of 23
PHYSICAL REVIEW B
94
, 195440 (2016)
Nonlinear damping and dephasing in nanomechanical systems
Juan Atalaya,
1
,
*
Thomas W. Kenny,
2
M. L. Roukes,
3
and M. I. Dykman
1
1
Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824, USA
2
Department of Mechanical Engineering, Stanford University, Stanford, California 94305, USA
3
Kavli Nanoscience Institute and Departments of Physics, Applied Physics, and Bioengineering, Caltech, Pasadena, California 91125, USA
(Received 28 September 2016; published 28 November 2016)
We present a microscopic theory of nonlinear damping and dephasing of low-frequency eigenmodes in
nanomechanical and micromechanical systems. The mechanism of the both effects is scattering of thermally
excited vibrational modes off the considered eigenmode. The scattering is accompanied by energy transfer of
2

ω
0
for nonlinear damping and is quasielastic for dephasing. We develop a formalism that allows studying both
spatially uniform systems and systems with a strong nonuniformity, which is smooth on the typical wavelength
of thermal modes but is pronounced on their mean free path. The formalism accounts for the decay of thermal
modes, which plays a major role in the nonlinear damping and dephasing. We identify the nonlinear analogs of the
Landau-Rumer, thermoelastic, and Akhiezer mechanisms and find the dependence of the relaxation parameters
on the temperature and the geometry of a system.
DOI:
10.1103/PhysRevB.94.195440
I. INTRODUCTION
Nanomechanical and micromechanical vibrational systems
have been attracting much attention in recent years. Vibrational
modes in such systems are not only interesting on their own, but
also provide a platform for studying physics far from thermal
equilibrium and quantum physics at the macroscale [
1
11
].
The typical spatial scale of the low-lying modes is the size
of the system, and their frequencies are
10
5
–10
8
Hz. They
have enabled record high force resolution [
12
15
] and mass
sensitivity [
16
18
] and have found applications in modern
electronic devices such as accelerometers and gyroscopes.
One of the most important characteristics of nanoscale
vibrational modes is their decay rate, which can be much
smaller than the vibration frequency [
1
3
,
12
18
]. The relax-
ation mechanisms include radiation into the bulk modes of
the medium surrounding the nanosystem [
19
23
], nonlinear
coupling between the modes localized inside the system
[
10
,
24
28
], and coupling with electronic excitations [
29
32
]
and two-level defects [
33
35
].
Conventionally, decay of a vibrational mode is associated
with an exponential falloff of the vibration amplitude in
time. This behavior is often observed in the experiment. In
the phenomenological description of the mode dynamics, it
comes from the friction force
2

̇
q
, where
q
is the vibration
coordinate and

is the friction coefficient. Such friction is
called linear. Recent experiments have shown that, already
for moderate vibration amplitudes, in various nanomechan-
ical and micromechanical systems the amplitude decay is
nonexponential in time and/or the maximal amplitude of
forced vibrations nonlinearly depends on the force amplitude
[
28
,
36
39
]; respectively, the friction force is a nonlinear
function of the mode coordinate and velocity. Nonlinear fric-
tion (nonlinear damping) has recently attracted considerable
attention also in the context of quantum information processing
with microwave cavity modes, where it was engineered to be
strong [
40
,
41
].
*
Current address: Department of Electrical Engineering, UC River-
side, Riverside, CA 92521.
Nonlinear damping is well understood in quantum optics,
where this is a major mechanism that limits the laser
intensity [
42
]. Significantly less is known about mechanical
systems. Here, a simple microscopic mechanism is a decay
process in which two quanta of the considered mode scatter
into a quantum of a mode of a continuous spectrum [
43
].
The corresponding process is sketched in Fig.
1(a)
.This
mechanism was implemented in [
40
,
41
]. For graphene-based
nanoresonators, it was discussed for nonlinear leakage of
the resonator modes into bulk modes [
23
], and it should be
relevant for the decay into bulk modes (clamping losses) in
other systems. However, as we show, in many cases other
mechanisms become more important.
In this paper, we develop a theory of nonlinear damping
due to the anharmonic coupling between vibrational modes of
a nanosystem/microsystem, which serves as a resonator. We
consider the damping of a low-lying eigenmode. The mode
eigenfrequency
ω
0
is assumed to be small compared to
k
B
T/

and the Debye frequency of a crystal. At the same time, the
mode decay rate is much smaller than
ω
0
, i.e., the mode has a
high quality factor.
An important feature of nanosystems and microsystems is
nonuniformity. Familiar examples are ripples of the graphene
sheets [
44
], bending and twisting of nanotubes, or layer
thickness variations in multilayer nanobeams. The scale of the
nonuniformity
l
sm
should be compared with the wavelength
λ
T
and the mean free path
l
T
of thermal phonons, i.e., of phonons
with frequencies
k
B
T/

. If both
λ
T
and
l
T
in the absence
of the nonuniformity are small,
λ
T

l
T

l
sm
, phonon
transport can be described in terms of coordinate-dependent
transport coefficients, for example, a coordinate-dependent
thermal conductivity [
45
]. A qualitatively different situation
occurs if
l
sm
<l
T
, as in rippled graphene membranes for low
temperature [
44
,
46
]. One has to be careful in defining
l
T
in
this case, as the normal vibrational modes are no longer plain
waves, and
l
T
should be defined as the mean free path of
these modes due to their nonlinear coupling with each other.
In this case, one can no longer use locally defined transport
coefficients even where
l
sm

λ
T
.
The aim of this paper is twofold: first, to develop a theory of
damping of low-frequency vibrational modes in nonuniform
2469-9950/2016/94(19)/195440(23)
195440-1
©2016 American Physical Society
ATALAYA, KENNY, ROUKES, AND DYKMAN
PHYSICAL REVIEW B
94
, 195440 (2016)
(a)(b)
0
0
0
0
FIG. 1. Elementary phonon-phonon interaction processes leading
to nonlinear damping. (a) Three-phonon scattering: two quanta of the
mode of interest, which has frequency
ω
0
, scatter into a phonon of
a continuous spectrum with frequency
ω
κ
2
ω
0
. (b) Four-phonon
scattering: a phonon of a continuous spectrum is scattered off the
considered mode into another phonon, with the frequency change
ω
κ

ω
κ
2
ω
0
. The modes
κ,κ

are not assumed to be plane waves.
systems and, second, to study the nonlinear damping and
dephasing. The theory is based on the eikonal formulation
of the kinetics of thermal phonons. We use the term “phonon”
loosely to describe normal modes of a nanosystem that belong
to the quasicontinuous spectrum, even though they are not
plane waves. The phonons are enumerated by a quantum
number
κ
which, for an ideal crystal, includes the wave vector
and the branch number, and have frequencies
ω
κ
.
The nonlinear damping with decay into a single phonon [see
Fig.
1(a)
] comes primarily from the cubic anharmonicity of the
vibrations [
43
]. As mentioned earlier, this mechanism can be
important for the nonlinear coupling to substrate phonons.
However, it is less likely to work for decay into the modes
localized within the system. This is because the localized
low-frequency modes have a discrete spectrum, and the decay
of the considered mode into one of these modes
κ
requires
resonance
ω
κ
2
ω
0
. It also requires that the decay rate of
the mode
κ
be much larger than the coupling and the decay
rate of the considered mode. Such fine tuning does not happen
generically, but if it happened, the resulting nonlinear damping
would be strong.
A moreimportant contributiontothenonlinear dampingcan
come from the processes described by the quartic anharmonic-
ity. In such processes, a phonon with frequency
ω
κ

2
ω
0
is
scattered off the considered low-frequency mode into another
high-frequency phonon [see Fig.
1(b)
]. It is implied that the
spectrum of the involved phonons
κ,κ

is quasicontinuous.
This condition is met if the phonons are weakly confined
to the resonator or
l
T
is small compared to at least one of
the resonator dimensions. The frequency change 2
ω
0
can be
smaller than the decay rate of the involved high-frequency
modes, which complicates the analysis.
The conventional mechanisms of linear damping for bulk
modes [
47
,
48
] are based on the processes similar to those
in Fig.
1(b)
, except that the scattering of the high-frequency
phonons is accompanied by absorption/emission of one quan-
tum of the considered low-frequency mode. Linear damping
comes from the cubic nonlinearity. The results of this paper
describe such damping in nonuniform systems.
The rates of linear and nonlinear damping can become
comparable even for comparatively small vibration amplitudes
as a consequence of symmetry restrictions or the restrictions
imposed on linear damping by the energy and momentum con-
servation. The best-known example is linear dissipation due to
quartic anharmonicity proposed by Landau and Khalatnikov
[
49
,
50
] to describe phonon scattering in liquid helium. Another
example is suppression of the cubic-anharmonicity-induced
linear damping for flexural modes in flat atomically thin
membranes. It is rooted in the reflection symmetry, which
allows nonlinear coupling that involves only an even number
of flexural modes. Ultimately, the linear damping is due to
four-wave processes, with participation of four flexural modes
or two flexural modes and two acoustic phonons [
25
]. The
nonlinear damping shown in Fig.
1(b)
occurs in the same order
in the anharmonicity and therefore has no additional smallness
associated with the higher-order anharmonicity. This might be
behind the numerically observed strong nonlinear damping in
graphene nanomembranes [
51
].
Along with damping, scattering of phonons off the con-
sidered mode leads to its dephasing, i.e., loss of coherence,
as was first discussed for high-frequency vibrations localized
near defects in solids [
52
,
53
]. In contrast to damping, the
scattering that leads to dephasing is quasielastic, the energy of
the mode is not changed. To the lowest order, the process comes
from the quartic nonlinearity in the mode coupling, although
the parameters are renormalized by the cubic nonlinearity,
which is also the case for nonlinear friction. As we show, the
formalism developed to describe nonlinear friction extends to
the analysis of the dephasing rate. The dephasing we consider
is a
T
2
process, in the language of nuclear magnetic resonance.
We will not consider 1
/f
-type frequency noise, which usually
does not come from nonlinear mode coupling.
The problem of linear damping of low-frequency eigen-
modes in nanoresonators and microresonators overlaps with
the problem of sound absorption in dielectrics, which has
been intensely studied since 1930s [
54
,
55
] (see [
45
,
56
60
]
and references therein for more recent work). However, to the
best of our knowledge, the nonlinear damping has not been
discussed nor has there been developed a theory for systems
with the nonuniformity on the scale small compared to the
mean free path of thermal phonons.
Outline
Our analysis of damping and dephasing of an underdamped
low-frequency mode is based on the expressions that relate
the relaxation rates to the correlators of the bath to which
the mode is coupled (see Sec.
II
). These correlators are the
main object of interest. We consider the case where the bath
is formed by phonons and the coupling to the bath is due to
the vibration nonlinearity. Then, the correlators of interest are
those of the pairs of the phonon variables and are related to
the phonon density matrix weighted with the interaction. They
are easy to calculate and to find the nonlinear damping where
the frequency
ω
0
of the considered mode is large compared to
the decay rates of thermal phonons, the Landau-Rumer limit
(see Sec.
III
). The problem becomes more complicated where
the phonon decay cannot be disregarded, and this is the central
part of the paper.
The analysis of the phonon correlators is done for nonuni-
form systems with an arbitrarily strong nonuniformity, which
is however smooth on the typical wavelength
λ
T
of thermal
modes. It is based on the eikonal approximation, as mentioned
195440-2
NONLINEAR DAMPING AND DEPHASING IN . . .
PHYSICAL REVIEW B
94
, 195440 (2016)
above. This approximation is developed in Sec.
IV
.The
thermal modes (phonons) are characterized by coordinate-
dependent wave vectors, and there holds approximate local
momentum conservation in phonon scattering. The analysis of
the dynamics can be conveniently done by writing the sought
correlators in the Wigner representation. We formulate it in
Sec.
V
in terms of the eigenmodes of the nonuniform system
(not plane waves). The central result of this section, Eq. (
33
),
relates the nonlinear decay rate of the considered system to the
Wigner transform of the thermal phonon correlator

α
(
R
,
k
,t
),
which is a function of the position
R
, the wave vector
k
, and
the phonon mode index
α
.
In Sec.
VI
, we discuss the quantum transport equation for
function

α
. The derivation of this equation is somewhat
cumbersome and is described in Appendixes A–D. Our
analysis suggests that, for smooth nonuniformity, this equation
is local in space even where the characteristic nonuniformity
length
l
sm
is small compared to the phonon mean free path
l
T
, but still
l
sm

(
l
T
λ
T
)
1
/
2
. The locality means that the
scattering couples function

α
(
R
,
k
,t
) to functions

α

(
R
,
k

,t
)
with different wave vectors and mode indices, but the same
position. In Sec.
VII
and Appendix
E
, we briefly consider the
effect on the evolution of function

α
of a weak short-range
disorder.
The analysis of the quantum transport equation in uniform
systems is done in Sec.
VIII
. Function

α
is expanded
in the eigenfunctions of the collision operator, which is
advantageous, since this operator is independent of the
coordinates in this case. The results are used in Sec.
IX
to
study nonlinear thermoelastic relaxation of a low-frequency
mode. We show that this relaxation is of particular interest for
Lam
́
e modes in thin plates and also for flexural modes in thin
nanobeams, where thermal diffusion transverse to the beam
is fast compared to the mode vibrations; this limit is opposite
to the one where the Zener linear thermoelastic damping is
dominating. Along with the microscopic theory we develop a
phenomenological theory of nonlinear thermoelastic damping
and show that they lead to the same results.
In Sec.
X
, we briefly discuss the nonlinear Akhiezer
relaxation. As in the case of linear Akhiezer relaxation, it
is important where the mode frequencies are large compared
to the reciprocal time of thermal diffusion. We show that the
expressions for the Akhiezer decay rate are grossly simplified
in the
τ
approximation for the dynamics of thermal phonons;
however, the value of
τ
is generally different from the one used
to characterize thermal diffusion.
In Sec.
XI
, we discuss the dephasing of low-frequency
vibrational modes due to their coupling to thermal phonons.
For high-frequency modes, where the decay of thermal
phonons can be disregarded, the theory is similar to the theory
of dephasing for localized vibrations in solids [
52
,
53
]. For
smaller
ω
0
, the decay of thermal phonons has to be taken into
account. As we show, this can be done by simply extending the
results on nonlinear friction. Finally, in Sec.
XII
we summarize
the main conclusions of this work.
Notations.
Throughout the paper, we use a hat to indicate
tensors; for example, ˆ
is the strain tensor with components
ij
and
ˆ
θ
is an auxiliary tensor with components
θ
ij
.Weuse
a wide hat, like
̂
, for operators that describe the coupling
between different phonon correlators.
II. GENERAL EXPRESSION FOR THE NONLINEAR
DAMPING RATE
We consider a low-frequency mode of a mesoscopic
resonator interacting with a thermal bath. The bath excitations
are vibrational modes with a quasicontinuous frequency
spectrum, which we call phonons. The total Hamiltonian
H
=
H
0
+
H
b
+
H
i
is the sum of the Hamiltonians of the
mode
H
0
, the bath
H
b
, and their interaction
H
i
. We write
H
in
terms of the creation and annihilation operators of the mode,
a
and
a
, and the phonons,
b
κ
and
b
κ
, where
κ
enumerates the
phonons.
The displacement field of the mode of interest can be written
as
q
(
t
)
u
(
r
). Function
u
(
r
) describes the spatial structure of
the mode, whereas
q
(
t
) characterizes the displacement at the
antinode. We choose the normalization of
u
(
r
) in such a way
that the Lagrangian of the mode in the harmonic approximation
is (
M
0
/
2)(
̇
q
2
ω
2
0
q
2
), where
M
0
and
ω
0
are the effective mass
and the eigenfrequency of the mode. We then introduce the
operator of the mode coordinate as
q
=
q
0
(
a
+
a
)
,q
0
=
(

/
2
M
0
ω
0
)
1
/
2
.
(1)
Length
q
0
is the scaled amplitude of the zero-point vibrations
at the mode antinode. The characteristic wavelength of the
mode is of the order of the size of the system, and therefore
q
0
usually (but not necessarily) decreases with the system size.
The mode Hamiltonian
H
0
is
H
0
=

ω
0
a
a
+
1
2

Va
2
a
2
.
(2)
The term
V
comes from the internal nonlinearity of the
mode and from its nonlinear coupling to phonons (cf. [
43
]). It
plays an important role in the dynamics of mesoscopic modes
[
2
], even though typically
|
V
|
ω
0
. The energy decay rate is
not affected by this term, except where a mode is very strongly
driven [
28
,
61
]; such driving is not considered here.
In the Hamiltonian of the bath
H
b
, along with the harmonic
part it is necessary [
47
,
48
] to take into account the nonlinear
mode coupling responsible for the phonon decay, in particular,
the decay of the high-frequency phonons
κ,κ

involved in the
processes shown in Fig.
1(b)
. Respectively, we write
H
b
as
H
b
=

κ
ω
κ
b
κ
b
κ
+
H
ph-ph
,
H
ph-ph
=
1
2
κκ

κ

v
κκ

κ

b
κ
b
κ

b
κ

+
1
3
κκ

κ

v

κκ

κ

b
κ
b
κ

b
κ

+
H.c
.
(3)
The phonons
κ
are generally not plane waves. This is a
consequence of the finite size of the resonator and also of the
presence of built-in static deformation and/or stress and maybe
defects. However, in a deformed finite-size system one can still
find normal vibrational mode, and they are enumerated by
κ
.
In the coupling of the considered mode to the bath, we will
take into account the terms linear and quadratic in the mode
displacement
q
a
+
a
. The expression for
H
i
in terms of
the operators
a,a
can be written as
H
i
=
ah
1
+
a
h
1
+
a
2
h
2
+
a
2
h
2
+
a
ah
3
.
(4)
Operators
h
1
,
2
,
3
depend on the variables of the bath.
195440-3
ATALAYA, KENNY, ROUKES, AND DYKMAN
PHYSICAL REVIEW B
94
, 195440 (2016)
For weak coupling
H
i
, one can describe the dynamics of
the considered mode by switching to the rotating frame with
the canonical transformation exp(
0
a
at
). In slow time
compared to the vibration period 2
π/ω
0
one then obtains a
Markovian quantum kinetic equation [
41
,
43
,
62
]
̇
ρ
=−
i
1
2
V
[
a
2
a
2
]
+

D
[
a
]
ρ
+

(nl)
D
[
a
2
]
ρ
+

(
φ
)
D
(
φ
)
ρ,
(5)
where
ρ
is the mode density matrix and the superoperators
D
[
a
] and
D
[
a
2
] describe linear and nonlinear damping of the
mode
D
[
a
m
]
ρ
=−
(
̄
n
(
m
)
+
1)(
a
m
a
m
ρ
2
a
m
ρa
m
+
ρa
m
a
m
)
̄
n
(
m
)
(
a
m
a
m
ρ
2
a
m
ρa
m
+
ρa
m
a
m
)
.
(6)
Here,
̄
n
(
m
)
̄
n
(
0
) is the Planck number,
̄
n
(
ω
)
=
[exp(

ω/k
B
T
)
1]
1
. To the lowest order in
H
i
, the linear
and nonlinear damping rates

and

(nl)
are

=

2
[
̄
n
(
ω
0
)
+
1]
1
Re
0
dt e
0
t
h
1
(
t
)
h
1
(0)
,

(nl)
=

2
[
̄
n
(2
ω
0
)
+
1]
1
Re
0
dt e
2
0
t
h
2
(
t
)
h
2
(0)
.
(7)
These rates are determined by the processes in which the
considered mode makes a transition to the nearest energy
level, in the case of linear damping,
|
n
→|
n
1
,orto
the next-nearest level, in the case of nonlinear damping,
|
n
→|
n
2
, with the energy going into the thermal bath
(
|
n
is the mode Fock state).
The superoperator
D
(
φ
)
describes dephasing,
D
(
φ
)
ρ
=
2
a
aρa
a
(
a
a
)
2
ρ
ρ
(
a
a
)
2
.
(8)
The dephasing rate

(
φ
)
is

(
φ
)
=

2
Re
0
dt
h
3
(
t
)
h
3
(0)
.
(9)
As indicated above, the dephasing we consider is a
T
2
-type
process; it is caused by quasielastic scattering of the bath
excitations off the mode, with no transitions between the mode
energy levels.
In Eqs. (
7
) and (
9
), the averaging
...
is done over the
thermal state of the bath. We assume that
h
1
,
2
,
3
=
0. The
quantum kinetic equation applies for the relaxation rates small
compared to
ω
0
and to the reciprocal correlation time of the
bath excitations, i.e., to the reciprocal time over which the
correlators in Eqs. (
7
) and (
9
) decay. Also, it is required that
the derivatives of
,
(nl)
over
ω
0
be small. These conditions
are met if the density of states of the bath weighted with the ap-
propriate interaction
h
1
,
2
,
3
is smooth near
ω
0
,2
ω
0
, and
ω
=
0.
Along with decay and dephasing, the coupling (
4
) leads
to a renormalization of the mode eigenfrequency
ω
0
and the
parameter
V
that describes the nonequidistance of the mode
energy levels. We assume this renormalization to have been
done; we note that the renormalization of
ω
0
is temperature
dependent even to the second order in
H
i
. In addition, the
linear in the oscillator coordinates coupling
h
1
generally
renormalizes the operator of the nonlinear coupling
h
2
.The
corresponding contribution to nonlinear friction can be thought
of as a result of a nonlinear response of the bath to the
perturbation
ah
1
+
a
h
1
. It does not lead to new effects.
A. Relation to the phenomenological nonlinear friction
The nonlinear friction term of the form (
5
) corresponds
to the nonlinear friction in the van der Pol equation, which
is broadly used for describing various types of self-excited
modes, in particular, in laser physics [
42
]. This can be seen
from the equation for
a
=
Tr(
) that follows from Eq. (
5
):
̇
a
=−
(

+

(
φ
)
4
̄
n
(2)

(nl)
)
a
(
iV
+
2

(nl)
)
a
a
2
.
(10)
The standard van der Pol equation in the rotating wave
approximation is obtained from this equation if one ignores
fluctuations, i.e., formally replaces
a
a
2
a
a
2
and
disregards the thermal term
̄
n
(2)
.Theterm
a
a
2
describes
losses which nonlinearly depend on the mode amplitude;
because of this term, the decay of
a
(
t
)
is nonexponential. We
note that, in contrast to the nonlinear friction usually discussed
in lasers, in our case the dissipation and the fluctuations in
the quantum kinetic equation (
5
) come from the coupling
to a bath in thermal equilibrium and are related by the
fluctuation-dissipation relation.
Parameters
V
and

(nl)
have dimension of frequency and
describe, respectively, the changes of the oscillator frequency
and the decay rate, which occur for the oscillator displacement
on the order of the zero-point vibration amplitude
q
0
[Eq. (
1
)].
In the laboratory frame, in the phenomenological classical
description of the oscillator dynamics, with account taken of
the fluctuations that come along with the linear and nonlinear
damping, Eq. (
10
) corresponds to the equation of motion of
the form
̈
q
+
ω
2
0
q
+
γq
3
=−
2

̇
q
4

(nl)
(
q/q
0
)
2
̇
q
+
ξ
lin
(
t
)
+
ξ
(
t
)
q/q
0
.
(11)
As seen from this equation, the friction force nonlinearly
depends on the displacement
q
(
t
). We note that,
phenomenologically, the friction force could have been written
as
(4
/
3
ω
2
0
q
2
0
)

(nl)
̇
q
3
. This would lead to the same equation
of motion in the rotating frame. Parameter
γ
=
(2
ω
0
/
3
q
2
0
)
V
in
Eq. (
11
) is the coefficient of nonlinearity of the restoring force.
Both
V
and

(nl)
are proportional to the corresponding classical
nonlinearity parameters in Eq. (
11
) multiplied by
q
2
0

.
The term
ξ
lin
(
t
)inEq.(
11
) represents the additive thermal
Gaussian noise with the intensity
T
related to the linear
friction coefficient by the fluctuation-dissipation theorem. The
term
ξ
(
t
) is also noise, but it has two distinct parts. One is
Gaussian noise at frequencies close to 2
ω
0
. It is related to
the nonlinear friction, has intensity

(nl)
T
, and the power
spectrum, which is flat around 2
ω
0
in the region much broader
than
,
(nl)
. The other part is noise with a bandwidth smaller
than
ω
0
yet white on the slow time scale
1
/,
1
/
(nl)
; its
intensity is given by

(
φ
)
(cf. Ref. [
63
] and papers cited
therein). We note that, in contrast to the phenomenological
model, Eq. (
5
) does not assume that the dynamics of the
oscillator is Markovian in the laboratory frame; it is Markovian
only on the time scale that largely exceeds
ω
1
0
and the
correlation time of the thermal reservoir.
195440-4
NONLINEAR DAMPING AND DEPHASING IN . . .
PHYSICAL REVIEW B
94
, 195440 (2016)
III. NONLINEAR DAMPING IN THE
LANDAU-RUMER LIMIT
The goal of this paper is to calculate the correlators (
7
)
and (
9
) for the case where operators
h
1
,
2
,
3
are determined by
the coupling to phonons. In what follows, we concentrate on
the analysis of nonlinear damping. We then briefly mention
the relation to linear damping, and also discuss dephasing.
To calculate the rate of nonlinear damping, we keep in
h
2
[Eq. (
4
)] the relevant terms of the first and second order in the
phonon coordinates,
h
2
=
κ
v
κ
b
κ
+
κκ

v
κ

κ
b
κ

b
κ
.
(12)
Here, the first term describes the decay of two quanta of the
considered mode into a phonon of the continuous spectrum
[Fig.
1(a)
]. The second term describes scattering of phonon
κ
and two quanta of the considered mode into a phonon
κ

[Fig.
1(b)
].
In the single-phonon decay, the considered mode pri-
marily emits substrate phonons, as mentioned earlier. Pa-
rameters
v
κ
in Eq. (
12
) determine nonlinear coupling
to such phonons at the boundary of the resonator. In
the case of the Fermi resonance between the con-
sidered mode and another weakly damped intracavity
mode, 2
ω
0
ω
κ
[
64
],
v
κ
determines the strength of the mode
coupling. The analysis of nonlinear dynamics in this case is
beyond the scope of this paper.
The two-mode coupling parameters
v
κ

κ
in Eq. (
12
)areof
particular interest for coupling to high-frequency phonons,
which usually have a high density of states. Since the
considered mode is a standing wave and the coupling depends
on the displacements of the modes, operator
v
κκ

b
κ
b
κ

can
be made Hermitian,
v
κκ

=
v
κ

κ
.InEq.(
12
), we disregarded
processes where two quanta of the considered mode decay into
two phonons
ω
κ
+
ω
κ

2
ω
0
. For small
ω
0
, such phonons
have low density of states, and decay into them is generally
less probable than the single-phonon decay or the decay via
scattering of high-frequency phonons.
The expression for the nonlinear damping rate (
7
) takes
a simple form in the Landau-Rumer limit where one can
disregard the decay of the phonons of the continuous spectrum,
i.e., their decay rate is small compared to 2
ω
0
. The Landau-
Rumer nonlinear damping rate is

(nl)
LR
=
π

2
κ
v
2
κ
δ
(
ω
κ
2
ω
0
)
+
π

2
κκ

v
2
κ

κ
(
̄
n
κ
̄
n
κ

)
×
δ
(
ω
κ

ω
κ
2
ω
0
)
,
(13)
where
̄
n
κ
̄
n
(
ω
κ
). The first term in Eq. (
13
) describes the
single-phonon decay and is independent of temperature.
The second term describes the decay due to scattering of
high-frequency phonons; the scattering rate is determined by
the thermal population of these phonons and thus depends on
temperature.
If the high-frequency phonons
κ,κ

that contribute to the
second term in Eq. (
13
) are of acoustic type and
ω
κ
κ


ω
0
,
the coupling parameters are
v
κ

κ

ω
κ
. This relation applies
also if one thinks of the coupling
v
κκ

as coming from
a parametric modulation of the energy of high-frequency
phonons by the strain from the considered mode [
55
]. The
squared displacement of the considered mode gives a factor

/M
0
ω
0
in
v
κκ

.For

ω
0

k
B
T
,inEq.(
13
)
̄
n
κ
̄
n
κ

(2

ω
0
/k
B
T
)
̄
n
κ
(
̄
n
κ
+
1). Therefore, for scattering by acoustic
phonons, the second term in the nonlinear damping rate
(
13
)is
C
v
T/ρ
0
c
2
s
, where
C
v
is the specific heat of the
nanomechanical resonator,
c
s
is the speed of sound, and
ρ
0
is the resonator mass density. Typically, this factor is small,
C
v
T/ρ
0
c
2
s

1.
In obtaining Eq. (
13
) we disregarded the rate of inelastic
scattering of phonons
κ
compared to
ω
0
, but the static disorder
was not assumed weak. This disorder makes high-frequency
vibrational modes different from plane waves. Equation (
13
)
applies both for short- and long-wavelength disorder, one just
has to use the normal modes
κ
calculated with the disorder
taken into account. However, the inelastic decay rate of the
modes
κ
should exceed the spacing of their frequencies.
Equation (
13
) also applies if modes
κ,κ

are weakly decaying
bulk modes coupled to the considered mode at the resonator
boundary.
IV. EIKONAL APPROXIMATION
We now consider how the nonlinear coupling between the
modes of the quasicontinuous spectrum (phonons) affects the
nonlinear damping rate

(nl)
. Of primary interest is the case of
smooth spatial nonuniformity of the resonator, with the spatial
scale
l
sm
large compared to the wavelength
λ
T
of thermal
modes, i.e., the modes with energy
k
B
T
. We will concentrate
on the nontrivial case where the mean free path
l
T
of the
thermal modes due to their nonlinear coupling is large, so that
for acoustic-type modes
l
T

l
sm

λ
T
T
=
2
π

c
s
/k
B
T
(14)
(
c
s
is the speed of sound). We will further assume that the
thermal phonons are not affected by a magnetic field and the
crystal is nonpolar. To simplify notations, we will assume that
there is just one atom per unit cell. The analysis immediately
extends to several atoms per cell.
For smooth nonuniformity, thermal phonons can be de-
scribed in the eikonal approximation (see Appendix
A
). The
displacement operator of an
n
th atom (at equilibrium position
r
n
) can be written as
u
n
=
κ
[

/
2
M
n
ω
κ
]
1
/
2
u
κ
(
r
n
)
b
κ
+
H
.
c
.,
(15)
where
M
n
is the atomic mass. Instead of being plane waves,
the long-wavelength normal modes
u
κ
(
r
)havetheform
u
κ
(
r
)
=
N
1
/
2
e
κ
(
r
)exp[
iS
κ
(
r
)]
,
S
κ
(
r
)
=
k
κ
(
r
)
.
(16)
One can think of Eq. (
16
) as the long-wavelength continuation
of the displacement field of the mode from the discrete values
u
κ
(
r
n
) defined on the lattice sites.
Vector
k
κ
plays the role of the coordinate-dependent wave
vector of the mode
κ
, and
e
κ
describes the polarization of the
mode;
N
is the number of the unit cells. Generally, not only
|
e
κ
|
2
depends on the coordinate
r
, but also the relation between
different components of
e
κ
varies with
r
. To the leading order,
in the eikonal approximation the equation of motion for
u
κ
(
r
)
is an algebraic equation that gives
k
κ
(
r
) for a given mode
frequency and branch, cf. [
65
] (see Appendix
A
).
195440-5
ATALAYA, KENNY, ROUKES, AND DYKMAN
PHYSICAL REVIEW B
94
, 195440 (2016)
The eikonal approximation implies smoothness of the
vectors
k
κ
(
r
) and
e
κ
(
r
). The spatial derivatives like
|
k
κ
|
and
(
k
κ
)
e
κ
are small compared to
k
2
κ
and
|
e
κ
|
k
2
κ
, respectively.
An important case is where one or two of the nanores-
onator dimensions are comparable to or smaller than the
thermal wavelength
λ
T
. This is of interest for nanowires or
thin membranes, including carbon nanotubes and graphene
membranes. Here, modes propagate within the resonator in
one or two directions spanned by vector
r

, whereas in
the transverse direction the modes are “quantized,” forming
different transverse branches. In Eq. (
16
), the eikonal now is a
function of the propagation direction,
S
κ
S
κ
(
r

), with
r

S
κ
(
r

)
=
k
κ
(
r

)
.
The polarization vector
e
κ
(
r
) smoothly depends on the coor-
dinate
r

.
The mode number
κ
has discrete and quasicontinuous com-
ponents. For an ideal bulk system these are the phonon branch
α
κ
and the wave vector
k
κ
. For thin nanoresonators, the discrete
components of
κ
also enumerate the transverse branches. It is
important that
κ
has a quasicontinuous component and, as
mentioned above, the frequency spectrum is quasicontinuous,
at least in the range where

ω
κ
k
B
T
.
One can find the quasicontinuous values of
κ
by formally
extending the system to make it possible to impose periodic
boundary conditions. One then requires that the difference of
S
κ
(
r
)[or
S
κ
(
r

)] on the boundaries be a multiple of 2
π
.The
quasicontinuous part of
κ
can be associated with the value of
vector
k
κ
at some arbitrary point
r
0
:
[
κ
]
quasicontinuous
⇐⇒
k
κ
(
r
0
)
.
(17)
Normal vibrational modes introduced in Eq. (
15
) satisfy
the standard orthonormality and completeness conditions; in
particular,
n
u
κ
(
r
n
)
·
u
κ

(
r
n
)
=
δ
κκ

.
(18)
The operator of the nonlinear phonon coupling
H
ph-ph
(
3
)
is obtained from the cubic in atomic displacements term in
the potential energy of the system. This term is a sum of local
contributions that depend on the difference of displacements
of a few neighboring atoms. From the expansion (
16
)wesee
that, if we disregard umklapp processes, the matrix elements
of the nonlinear mode-mode coupling in Eq. (
3
) can be
calculated in the stationary phase approximation. We change
from summation over lattice sites to integration
n
v
1
c
(
r
)
d
r
,
where
v
c
(
r
) is the volume of the unit cell that smoothly depends
on
r
. Then,
v
κ
1
κ
2
κ
3
d
r
e
i
[
S
κ
3
(
r
)
S
κ
1
(
r
)
S
κ
2
(
r
)]
ˆ
U
κ
1
κ
2
κ
3
(
r
)
e
i
[
S
κ
3
(
r
)
S
κ
1
(
r
)
S
κ
2
(
r
)]
ˆ
U
κ
1
κ
2
κ
3
(
r
)
.
(19)
Here,
r
r
(
κ
1
2
3
) is given by equation
k
κ
1
(
r
)
+
k
κ
2
(
r
)
k
κ
3
(
r
)
=
0
(20)
and
ˆ
U
κ
1
κ
2
κ
3
(
r
) is a tensor that describes the nonlinear phonon
coupling. In systems with translational symmetry,
ˆ
U
is
independent of
r
as are also vectors
k
κ
[
48
], and then
Eq. (
20
) becomes just the condition of the quasimomentum
conservation in phonon scattering. The matrix elements
v

κ
1
κ
2
κ
3
in Eq. (
3
) are given by an expression similar to (
19
). In Eq. (
19
),
we skipped the smooth mode polarization factors and
v
c
(
r
);
matrix elements
v
κ
1
κ
2
κ
3
are determined by the values of these
factors at
r
.
The stationary phase approximation applies provided the
nonuniformity is sufficiently strong. The range
δr
q
of the
values of
r
that contribute to the integral over
r
in the first line
of Eq. (
19
) can be estimated by expanding
S
κ
n
(
r
)(
n
=
1
,
2
,
3)
in Eq. (
19
) to second order in
r
r
. We assume that, for the
values of
κ
that correspond to thermal phonons, there holds
the condition
λ
T

δr
q

l
sm
,δr
q
=
2
r
S
κ
1
/
2
(
λ
T
l
sm
)
1
/
2
.
(21)
The above estimate of
δr
q
is obtained by noting that, for suf-
ficiently strong nonuniformity, the position-dependent wave
vector of thermal phonons
k
κ
(
r
) changes by
λ
1
T
on the
length
l
sm
. Parameter
δr
q
gives the momentum uncertainty in
the momentum conservation condition (
20
),
δk
(
δr
q
)
1
.
Similar arguments apply to the matrix elements
v
κκ

of the coupling of the considered low-frequency mode to
high-frequency modes
κ,κ

[see Eq. (
12
)]. For the typical
wavelength of the low-frequency mode larger than
l
sm
,from
Eq. (
21
) we can estimate the range of the difference in the
quasimomenta of the high-frequency phonons as
|
k
κ
(
r
)
k
κ

(
r
)
|

(
δr
q
)
1
.
The transition
n
v
1
c
(
r
)
d
r
and the stationary-phase
approximation used in Eq. (
19
) apply also if
k
κ
1
(
r
)
+
k
κ
2
(
r
)
k
κ
3
(
r
)
=
K
(
r
)
,
where
K
(
r
) is a smoothly dependent on
r
reciprocal lattice
vector. This expression describes umklapp processes in a
nonuniform medium. However, we concentrate on the tem-
perature range where typical
|
k
κ
(
r
)
||
K
(
r
)
|
.
V. TWO-PHONON CORRELATOR
It is convenient to express the nonlinear damping rate

(nl)
[Eqs. (
7
) and (
12
)] in terms of the two-phonon correlation
function
φ
κκ

(
t
). For

ω
0

k
B
T
,

(nl)
=−
2
ω
0
k
B
T
Im
κ,κ

v
κ

κ
0
dt e
i
(2
ω
0
+
)
t
φ
κκ

(
t
) (22)
(
ε
→+
0), where
φ
κκ

(
t
)
=
i

κ
0

0
v
κ

0
κ
0
b
κ
(
t
)
b
κ

(
t
)
b
κ

0
(0)
b
κ
0
(0)
.
(23)
Functions
φ
κκ

that contribute to the nonlinear damping
correspond to mode pairs
κ,κ

with close eigenfrequencies
|
ω
κ
ω
κ

|∼
2
ω
0

ω
κ
. For the values of
r
where
v
κ

κ
are large for given
κ,κ

, the wave vectors of the modes
κ,κ

are close,
|
k
κ
(
r
)
k
κ

(
r
)
||
k
κ
(
r
)
|
,
|
k
κ

(
r
)
|
.Themajor
contribution to

(nl)
comes from the modes
κ,κ

that belong
to the same vibrational branch
α
κ
=
α
κ

. Functions
φ
κκ

for
195440-6
NONLINEAR DAMPING AND DEPHASING IN . . .
PHYSICAL REVIEW B
94
, 195440 (2016)
different branches are fast oscillating on the scale
ω
0
, except
for the region where the branches cross or touch, but this region
is small, and therefore the overall contribution of terms with
α
κ

=
α
κ

is small.
The correlator
φ
κκ

is the weighted off-diagonal phonon
density matrix
b
κ
(
t
)
b
κ

(
t
)
. The weighting factor is indepen-
dent of time, as seen from Eq. (
23
). The reduction of the prob-
lem to a calculation of the weighted thermal-phonon density
matrix significantly simplifies the analysis. Our formulation
allows us to study spatially nonuniform systems, and, in the
first place, the nonlinear friction and the dephasing. Our results
for linear friction in the limit of a spatially uniform system and
weak disorder coincide with the results based on the linearized
Boltzmann equation.
A natural approach to the analysis of the phonon dynamics
for close
k
κ
(
r
) and
k
κ

(
r
) is to use the Wigner transformation
W
. This is a linear integral transformation. The Wigner
transform of the correlator
φ
κκ

for a nonuniform system is
tensor
ˆ

α
:
ˆ

α
(
R
,
k
,t
)
=
W
R
,
k
[
φ
κκ

]
κ,κ

ˆ
θ
(
R
,
k
;
κ,κ

;
α
)
φ
κκ

(
t
)
.
(24)
The transformation should be constructed in such a way
that correlator [

α
(
R
,
k
)]
ij
describes the dynamics of a wave
packet centered at
R
with the typical wave vector
k
and that
ˆ

α
smoothly depends on
R
on the scale
λ
T
. The kernel of
such transformation for a nonuniform system is tensor
ˆ
θ
with
components
θ
ij
(
R
,
k
;
κ,κ

;
α
)
=
d
ρ
v
c
(
R
)
u
κ

i
(
R
+
1
2
ρ
)
e
i
k
ρ
×
u
κ
j
(
R
1
2
ρ
)
δ
α
κ
α
δ
α
κ

α
.
(25)
Alternatively, in the eikonal approximation (
16
) for the lattice
displacements, tensor
ˆ
θ
could be defined as a sum over lattice
sites
n,m
of
u
κ

i
(
r
n
)
u
κj
(
r
m
)exp[
i
k
(
r
n
r
m
)] calculated for
R
=
(
r
n
+
r
m
)
/
2 and for
κ,κ

belonging to a branch
α
.The
typical value of
|
ρ
|
in Eq. (
25
)is
λ
T
; it is small compared to
the nonuniformity length
l
sm
.Wehaveset
v
c
(
R
+
1
2
ρ
)
v
c
(
R
1
2
ρ
)
v
2
c
(
R
).
In nanowires and thin membranes, vectors
R
,
k
,
ρ
in Eq. (
25
)
are one or two dimensional, respectively. In this case, Eq. (
25
)
implies summation over the lattice sites in the transverse
direction for the corresponding
α
; function
θ
ij
is independent
of the transverse coordinates. In a spatially uniform system
where
k
κ
is a good quantum number, from Eq. (
25
)
ˆ
θ
δ
[
1
2
(
k
κ
+
k
κ

)
k
]exp[
i
(
k
κ

k
κ
)
R
]; it has the form of the
kernel of the standard Wigner transformation in the continuous
limit.
The orthonormality of functions
u
κ
(
r
)[Eq.(
18
)] leads to
the relation
d
R
d
k
(2
π
)
d
Tr[
ˆ
θ
(
R
,
k
;
κ
1

1
;
α
1
)
ˆ
θ
(
R
,
k
;
κ,κ

;
α
)]
=
δ
κκ
1
δ
κ

κ

1
δ
α
κ
α
δ
α
κ

α
δ
αα
1
,
(26)
where
d
is the dimension of vectors
R
and
k
and the trace is
taken over the tensor indices.
Scalar Wigner function in the eikonal approximation
For smooth nonuniformity, one can simplify tensor
ˆ
θ
by
integrating over
ρ
in Eq. (
25
) with the account taken of
the eikonal form of function
u
κ
[Eq. (
16
)]. Expanding the
exponents
S
κ
(
R
1
2
ρ
) and
S
κ

(
R
+
1
2
ρ
) to first order in
ρ
,we
obtain
ˆ
θ
(
R
,
k
;
κ,κ

;
α
)
̃
δ
{
1
2
[
k
κ
(
R
)
+
k
κ

(
R
)]
k
}
,
(27)
where
̃
δ
(
k
)isasharp
δ
-like peak that describes the approximate
momentum conservation; its typical width is
(
δr
q
)
1

λ
1
T
and is determined by the second-order term in the expansion
of
S
κ
,S
κ

in
ρ
.
If
ˆ
θ
were defined as a lattice sum, which is of interest
for high temperatures, one would have quasimomentum
conservation in the form
k
κ
(
R
)
+
k
κ

(
R
)
k
+
K
(
R
). The
analysis of this case is a straightforward extension of the
present analysis.
Since
|
k
κ
(
R
)
k
κ

(
R
)
||
k
κ
(
R
)
|
and
α
κ
=
α
κ

,weintro-
duce a quasicontinuous component ̃
κ
, which is the “center” of
κ,κ

for given
R
,
k
and is defined by equation
k
̃
κ
(
R
)
=
k
,
̃
κ
̃
κ
(
R
,
k
)
.
(28)
Keeping in mind that the polarization vectors
e
κ
(
r
)inEq.(
16
)
are smooth functions of the quasicontinuous part of
κ
, we can
change from tensor
ˆ
θ
to a scalar
θ
:
ˆ
θ
(
R
,
k
;
κ,κ

;
α
)
ˆ
M
(
R
,
k
;
α
)
θ
α
(
R
,
k
;
κ,κ

)
,
θ
α
(
R
,
k
;
κ,κ

)
=
(2
π
)
d
N
v
c
(
R
)
̃
δ
[
k
κ
(
R
)
+
k
κ

(
R
)
2
k
]
×|
e
̃
κ
|
2
exp[
iS
κ

(
R
)
iS
κ
(
R
)]
δ
α
κ
δ
α
κ

.
(29)
Here,
e
̃
κ
e
̃
κ
(
R
), whereas
ˆ
M
(
R
,
k
;
α
) is a unit-trace tensor
M
ij
(
R
,
k
;
α
)
=
e
̃
κi
e
̃
κ
j
/
|
e
̃
κ
|
2
.
(30)
Using Eq. (
24
), we can now introduce a scalar correlation
function

α
(
R
,
k
,t
):
ˆ

α
(
R
,
k
,t
)]
ˆ
M
(
R
,
k
;
α
)

α
(
R
,
k
,t
)
,

α
(
R
,
k
,t
)
=
κκ

θ
α
(
R
,
k
;
κ,κ

)
φ
κκ

(
t
)
.
(31)
Similarly we introduce a scalar parameter of the coupling to
the low-frequency mode
V
α
(
R
,
k
) (the vertex in the Wigner
representation):
V
α
(
R
,
k
)
=
κ,κ

θ
α
(
R
,
k
;
κ,κ

)
v
κ

κ
.
(32)
Using the orthonormality of tensors
ˆ
θ
[Eq. (
26
)], after simple
algebra one obtains the expression for the nonlinear damping
195440-7
ATALAYA, KENNY, ROUKES, AND DYKMAN
PHYSICAL REVIEW B
94
, 195440 (2016)
rate that reads as

(nl)
=−
2
ω
0
k
B
T
Im
d
R
d
k
(2
π
)
d
α
V
α
(
R
,
k
)

α,
2
ω
0
(
R
,
k
)
,

α,ω
(
R
,
k
)
=
0
dt 
α
(
R
,
k
,t
)exp(
iωt
)
.
(33)
The correlator

α
(
R
,
k
,t
) is the scalar Wigner transform
of
φ
κκ

in a spatially nonuniform system, i.e., the Wigner
transform of the off-diagonal phonon density matrix weighted
with a time-independent operator. It is much more convenient
for the analysis than the tensor. Equation (
33
) expresses the
nonlinear damping parameter

(nl)
in terms of this scalar. We
note that
V
α
in Eq. (
33
) is real for a nonlinear coupling to a
standing mode. Interestingly, as shown in Appendix
C
, starting
with the dynamical equations for the tensor
ˆ

α
, one can obtain
a closed-form equation for the scalar function

α
(
R
,
k
,t
). It is
discussed in the next section.
VI. TRANSPORT EQUATION FOR THE
WIGNER FUNCTION
Function

α
(
R
,
k
,t
) smoothly depends on the coordinate
R
and the wave vector
k
. As a function of
R
,itvariesonthe
distance determined by the scale of the spatial nonuniformity
and the phonon mean free path. This distance is much longer
than
|
k
|
1
λ
T
. The separation of the spatial scales makes it
possible to derive a transport equation for function

α
.The
derivation is somewhat involved since the standard momentum
conservation in phonon scattering does not apply. It is given in
the Appendix for the Fourier transform of

α
(
R
,
k
). In the
time domain the transport equation reads as
t

α
+{

α
,
α
}=
St[

α
]
,

α

α
(
R
,
k
,t
)
,
α

α
(
R
,
k
)
.
(34)
Here,

α
(
R
,
k
)
=
κ
ω
κ
θ
α
(
R
,
k
;
κ,κ
)
(35)
is the effective position- and momentum-dependent frequency
of the branch
α
, and
{
A,B
}=
k
A∂
R
B
k
B∂
R
A
is the Poisson bracket.
The effective collision term St

α
is
St[

α
(
R
,
k
,t
)]
=
α
0
d
R
0
d
k
0
[
λ
α
0
α
(
R
,
R
0
,
k
,
k
0
)
]
ω
→+
i
0
×

α
0
(
R
0
,
k
0
,t
)
.
(36)
The kernel
λ
for phonon-phonon scattering is given by Eq. (
C6
)
with operator
̂
calculated for
ω
→+
i
0. The leading-order
term in [
λ
]
ω
→+
i
0
is real. It contains
δ
functions of energy
conservation in the processes where modes
κ,κ

decay into
other modes or are scattered by other modes. The imaginary
terms in
λ
, which are quadratic in the small parameters of the
phonon nonlinearity, lead to a small renormalization of

α
,
which we assume to have been done.
The initial condition for function

α
follows from Eq. (
31
)
in which one replaces
φ
κκ

(
t
) with
φ
κκ

(0):
φ
κκ

(0)
=−
i

1
v
κ

κ
̄
n
κ
(
̄
n
κ

+
1)
,

α
(
R
,
k
,
0)
=−
i

1
V
α
(
R
,
k
)
̄
n
̃
κ
(
̄
n
̃
κ
+
1)
,
(37)
where ̃
κ
̃
κ
(
R
,
k
) is given by Eq. (
28
).
A feature of central importance in Eq. (
36
)isthe
local-
ity
. As shown in Appendix
C1
, the kernel
λ
α
0
α
(
R
,
R
0
,
k
,
k
0
)
becomes small for
|
R
R
0
|
exceeding the quantum scale
δr
q
=
(
λ
T
l
sm
)
1
/
2
. On the other hand, for a given
R
R
0
such that
|
R
R
0
|

δr
q
, function
λ
α
0
α
varies with
R
on the
nonuniformity scale
l
sm
. It is this variation that determines
the spatial variation of

α
(
R
,
k
,t
) (recall that
l
sm
is small
compared to the phonon mean free path
l
T
). Therefore,
function

α
remains nearly constant on the quantum spatial
scale
δr
q

l
sm
, and one can rewrite Eq. (
36
)as
St[

α
(
R
,
k
,t
)]
=
α
0
d
k
0
̃
λ
α
0
α
(
R
,
k
,
k
0
)

α
0
(
R
,
k
0
,t
)
,
̃
λ
α
0
α
(
R
,
k
,
k
0
)
=
d
R
0
[
λ
α
0
α
(
R
,
R
0
,
k
,
k
0
)]
ω
→+
i
0
,
(38)
which means that the evolution of

α
(
R
,
k
,t
) described by
Eq. (
34
) depends on

α
0
(
R
,
k
0
,t
), i.e., on the Wigner function
calculated for the same
R
but different
k
0
0
.
The physics behind the locality is as follows. In phonon
scattering, the total phonon momentum is conserved to within

/δr
q
[cf. Eq. (
21
)]. Respectively, the region in which the
scattering occurs is
δr
q
. On the other hand, the scattering rate
changes on the nonuniformity scale
l
sm
that largely exceeds
δr
q
. On the scale
l
sm
, phonon scattering events look “local.”
In the Markovian approximation used to derive Eq. (
34
)
(the ladder diagram approximation), this equation has the
same form as the kinetic equation for the Wigner transform
of the phonon density matrix. The difference is in the initial
condition, which has the form (
37
) for function

α
.The
decoupling used to derive Eq. (
34
) disregards many-phonon
correlations. This equation is an analog of the linearized Boltz-
mann equation for the phonon distribution function. However,
here it is derived for a phonon correlation function starting
from a microscopic model of a nonuniform system with the
nonuniformity, which is smooth on phonon wavelength, but
not on the phonon mean free path.
VII. SHORT-RANGE DISORDER
The above analysis directly extends to the case where,
along with smooth nonuniformity, the crystal has a short-range
disorder, examples being mass disorder or impurities (cf.
Refs. [
45
,
57
,
60
] and references therein). This disorder leads to
scattering of thermal phonons, but plays practically no role in
the dynamics of the considered long-wavelength mode itself,
except that it may slightly change its waveform. We will
describe the disorder by the extra term in the Hamiltonian
of thermal phonons
H
d
=
κκ

v
(
d
)
κκ

b
κ
b
κ

;
(39)
195440-8
NONLINEAR DAMPING AND DEPHASING IN . . .
PHYSICAL REVIEW B
94
, 195440 (2016)
we disregard the terms with
b
κ
b
κ

and
b
κ
b
κ

, which lead to
frequency renormalization. We will further assume that the
disorder is weak. The simplest model is a weak zero-mean
Gaussian field, with correlations limited to a single lattice
site (a
δ
-correlated field, in the continuous limit). In this
case, parameters
v
(
d
)
κκ

for thermal phonons are sums over the
lattice sites
n
of the factors
u
κ

i
(
r
n
)
u
κj
(
r
n
) weighted with
the perturbation on site
n
. Then, the quadratic in the disorder
terms, which determine phonon scattering in the second order
of the perturbation theory, have the structure
v
(
d
)
κκ

v
(
d
)
κ
1
κ

1
n
v
(
d
)
)
κκ
1
κ

κ

1
(
r
n
)exp[
iS
κ

(
r
n
)
iS
κ
(
r
n
)
+
iS
κ

1
(
r
n
)
iS
κ
1
(
r
n
)]
,
(40)
where tensors ˆ
v
(
d
)
smoothly depend on
n
and on the mode
indices.
Similar to Eq. (
19
), the sum (
40
) can be written as an integral
over
r
n
because we are interested in the correlation functions
φ
κκ

with close
κ
and
κ

. This integral can be calculated in
the stationary phase approximation. The major contribution to
Eq. (
40
) comes from the vicinity of a point
r
∗∗
in which the
total local quasimomentum is conserved:
k
κ
(
r
∗∗
)
k
κ

(
r
∗∗
)
+
k
κ
1
(
r
∗∗
)
k
κ

1
(
r
∗∗
)
=
0
.
(41)
We note that, because of the smooth nonuniformity, the
quasimomentum of the system of phonons is not conserved
even without short-range disorder and umklapp processes.
Phonons do not propagate along straight lines. However,
Eq. (
20
) indicates that, locally, the quasimomentum of the
system of phonons is conserved. Short-range disorder breaks
down this conservation; in each scattering event, the local
momentum of a phonon is transferred to a short-range scatterer.
It is shown in Appendix
E
that, in the presence of short-
range disorder, one still obtains a transport equation for the
Wigner function of the form of Eqs. (
34
) and (
38
), but now
the kernel of the collision term
λ
α
0
α
has an extra term. This
term is quadratic in parameters
v
(
d
)
κκ

[see Eq. (
E1
)]. Effectively,
the approximation corresponds to a generalized Matthiessen
rule where the phonon-phonon and the short-range disorder
scattering processes do not affect each other. However, both
types of processes are affected by the smooth nonuniformity.
Short-range scattering is particularly important for lower tem-
peratures, where the rate of the umklapp processes becomes
exponentially small.
VIII. NONLINEAR FRICTION IN SPATIALLY
UNIFORM SYSTEMS
Solutions of the transport equation can be obtained for
specific types of spatial nonuniformity, and in most cases this
requires numerical calculations. In what follows, we illustrate
the power of our approach by applying it to spatially uniform
systems and calculating nonlinear damping in such systems.
In uniform systems
k
κ
(
r
) and

α
(
r
,
k
) are independent of
the coordinate
r
. The spatial dependence remains only in
the displacement of the considered mode and in its Wigner
function

α
. We therefore switch to using
r
instead of
R
for
the coordinate of the Wigner function. This will make other
notations more conventional; in particular, the displacement
of the mode will be a function of
r
, not
R
.
In a uniform system function,
θ
α
[Eq. (
29
)] becomes
θ
α
(
r
,
k
;
κ,κ

)
=
(2
π
)
d
V
δ
[
1
2
(
k
κ
+
k
κ

)
k
]
×
exp[
i
(
k
κ

k
κ
)
r
]
δ
α
κ
δ
α
κ

,
whereas in Eq. (
36
) the collision operator takes the form
St[

α
(
r
,
k
,t
)]
=
V
(2
π
)
d
α
0
d
k
0
̃
k
0
α
0
k
α

α
0
(
r
,
k
0
,t
)
.
(42)
Here,
V
=
N
v
c
is the volume of the
d
-dimensional system;
operator
̃
is given by Eq. (
C10
) and is expressed in terms
of the rates of phonon scattering (
C11
), (
C12
), and (
E2
). In
addition, in the transport equation (
34
) the Poisson bracket
becomes
{

α
,
α
}=
v
k
α
r

α
, where
v
k
α
=
k
ω
k
α
is the
group velocity of phonons of the branch
α
with the wave
vector
k
. We note that the phonon-phonon collision term in
Eq. (
42
) can be written in the same form as the collision term
in the linearized phonon kinetic equation for the correction to
the phonon occupation number [
48
].
We are looking for a nonstationary solution of Eq. (
34
) with
the initial condition (
37
) and for the Fourier transform of this
solution convoluted with function
V
α
(
r
,
k
). The expression
for
V
α
(
r
,
k
) and the initial condition for

α
are simplified
in a spatially uniform system. Indeed, in the considered
problem of decay of smoothly varying in space standing-wave
type modes, of interest are the parameters of coupling to
phonons
v
κ

κ
=
v
k

α

,
k
α
with
α

=
α
and
|
k
k

||
k
|
.They
are given by a Fourier transform
ζ
α
(
r
) of the squared strain
tensor components of the considered mode convoluted with
the polarization vectors of thermal phonons. Function
ζ
α
(
r
)
is calculated for the considered mode being in the ground
vibrational state. It smoothly depends on
r
on the scale
|
k
|
1
.
Thus,
v
k

α,
k
α
V
1
̃
v
α
[(
k
+
k

)
/
2]
d
r
exp[
i
(
k
k

)
r
]
ζ
α
(
r
),
and from Eq. (
32
) the coupling vertex in the Wigner represen-
tation is
V
α
(
r
,
k
)
=
̃
v
α
(
k
)
ζ
α
(
r
)
.
(43)
The coefficient ̃
v
α
(
k
) is quadratic in
k
for acoustic modes and
is inversely proportional to
ω
k
α
. One can think of
V
α
(
r
,
k
)
/

as a scaled local (at a point
r
) change of the frequency of
the mode (
k
) due to the squared strain from the considered
low-frequency mode in its ground vibrational state.
In a uniform system the density
ρ
ρ
0
is constant.
Therefore, if the displacement field of the considered mode
is
q
(
t
)
u
(
r
), we have the effective mass of the mode
M
0
=
ρ
0
d
r
[
u
(
r
)]
2
[note that
u
(
r
) is dimensionless].
A. Eigenvalue problem for the collision term
The collision operator St[

α
(
r
,
k
,t
)] is independent of
r
,
it couples the values of

α
(
r
,
k
,t
) with different
k
and
α
,but
with the same
r
. One can therefore introduce right and left
eigenvectors of operator St, which we denote by
ψ
ν
(
k
) and
195440-9