of 26
Nonlinear damping and dephasing in nanomechanical systems
Juan Atalaya
Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824, USA
Thomas W. Kenny
Department of Mechanical Engineering, Standford University, Stanford, California 94305, USA
M. L. Roukes
Kavli Nanoscience Institute and Departments of Physics, Applied Physics,
and Bioengineering, Caltech, Pasadena, California 91125, USA
M. I. Dykman
Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824, USA
(Dated: September 29, 2016)
We present a microscopic theory of nonlinear damping and dephasing of low-frequency eigenmodes
in nano- and micro-mechanical 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 quasieleastic 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 not 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.
PACS numbers: 05.40.-a
I. INTRODUCTION
Nano- and micro-mechanical vibrational systems have
been attracting much attention in recent years. Vibra-
tional modes in such systems are not only interesting
on their own, but also provide a platform for study-
ing 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 nano-scale
vibrational modes is their decay rate, which can be much
smaller than the vibration frequency [1–3, 12–18]. The
relaxation mechanisms include radiation into the bulk
modes of the medium surrounding the nanosystem [19–
23], nonlinear coupling between the modes localized in-
side 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 asso-
ciated with an exponential fall-off of the vibration am-
plitude 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
,
Current address: Department of Electrical Engineering, UC
Riverside, Riverside, CA 92521
where
q
is the vibration coordinate and Γ is the fric-
tion coefficient. Such friction is called linear. Recent ex-
periments have shown that, already for moderate vibra-
tion amplitudes, in various nano- and micromechanical
systems the amplitude decay is nonexponential in time
and/or the maximal amplitude of forced vibrations non-
linearly depends on the force amplitude [28, 36–39]; re-
spectively, the friction force is a nonlinear function of the
mode coordinate and velocity. Nonlinear friction (non-
linear damping) has recently attracted considerable at-
tention also in the context of quantum information pro-
cessing with microwave cavity modes, where it was engi-
neered to be strong [40, 41].
Nonlinear damping is well understood in quantum op-
tics, where this is a major mechanism that limits the laser
intensity [42]. Significantly less is known about mechan-
ical 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). It was implemented in [40, 41]. For graphene-
based nanoresonators, it was discussed for nonlinear leak-
age of the resonator modes into bulk modes [23], and it
should be relevant for the decay into bulk modes (clamp-
ing 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 damp-
ing due to the anharmonic coupling between vibrational
modes of a nano/micro system, which serves as a res-
onator. We consider the damping of a low-lying eigen-
arXiv:1609.08714v1 [cond-mat.mes-hall] 28 Sep 2016
2
(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 con-
tinuous 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.
mode. 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 nano- and microsystems is
nonuniformity.
Familiar examples are ripples of the
graphene sheets [44], bending and twisting of nanotubes,
or layer thickness variations in multi-layer nanobeams.
The scale of the nonuniformity
l
sm
should be compared
with the wavelength
λ
T
and the mean free path
l
T
of ther-
mal phonons, i.e., of phonons with frequencies
k
B
T/
~
.
If both
λ
T
and
l
T
in the absence of the nonuniformity
are both 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 oc-
curs if
l
sm
< l
T
, as in rippled graphene membranes for
low temperature [44, 46]. One has to be careful in defin-
ing
l
T
in this case, as the normal vibrational modes are no
longer plain waves, and
l
T
should be defined as the scat-
tering rate 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 systems and second, to study the non-
linear 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 quasi-
continuous spectrum, even though they are not plane
waves. The phonons are enumerated by a quantum num-
ber
κ
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 ear-
lier, 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 sys-
tem. This is because the localized low-frequency modes
have a discrete spectrum, and the decay of the consid-
ered 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 de-
cay rate of the considered mode. Such fine tuning does
not happen generically, but if it happened, the resulting
nonlinear damping would be strong.
A more important contribution to the nonlinear damp-
ing can come from the processes described by the quartic
anharmonicity. In such processes, a phonon with fre-
quency
ω
κ

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 quasi-continuous. 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 dimen-
sions. 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 simi-
lar to those in Fig. 1(b), except that the scattering of
the high-frequency phonons is accompanied by absorp-
tion/emission of one quantum of the considered low-
frequency mode. Linear damping comes from the cu-
bic 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 am-
plitudes as a consequence of symmetry restrictions or the
restrictions imposed by the energy and momentum con-
servation. The best known example where dissipation
is due to quartic anharmonicity was proposed by Lan-
dau and Khalatnikov [49, 50] to describe phonon scat-
tering 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 nonlin-
ear coupling that involves only an even number of flex-
ural 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 ob-
served strong nonlinear damping in graphene nanomem-
branes [51].
Along with damping, scattering of phonons off the con-
sidered mode leads to its dephasing, i.e., loss of coher-
ence, 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 quasi-
elastic, the energy of the mode is not changed. To the
lowest order, the process comes from the quartic non-
linearity 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 formal-
3
ism 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 mag-
netic 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 nano- and micro-resonators 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 the-
ory for systems with the nonuniformity on the scale small
compared to the mean free path of thermal phonons.
A. Outline
Our analysis of damping and dephasing of an under-
damped 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 den-
sity 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 compli-
cated 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
nonuniform systems with an arbitrarily strong nonunifor-
mity, which is however smooth on the typical wavelength
λ
T
of thermal modes. It is based on the eikonal ap-
proximation, as mentioned 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 correla-
tor Φ
α
(
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 some-
what cumbersome and is described in Appendices A to
D. Our analysis suggests that, for smooth nonuniformity,
this equation is local in space even where the character-
istic 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 uni-
form systems is done in Sec. VIII. Function Φ
α
is ex-
panded in the eigenfunctions of the collision operator,
which is advantageous, since this operator is indepen-
dent 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 phe-
nomenological theory of nonlinear thermoelastic damp-
ing and show that they lead to the same results.
In Sec. X we briefly discuss the nonlinear Akhiezer re-
laxation. As in the case of linear Akhiezer relaxation, it
is important where the mode frequencies are large com-
pared 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 dy-
namics of thermal phonons; however, the value of
τ
is
generally different from that used to describe 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 sim-
ilar 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 nonlin-
ear friction. Finally, in Sec. XII we summarize the main
conclusions of this work.
Notations:
Throughout the paper we use a hat to in-
dicate tensors; for example, ˆ

is the strain tensor with
components

ij
and
ˆ
θ
is an auxiliary tensor with compo-
nents
θ
ij
. We use a wide hat, like
̂
Λ, for operators that
describe the coupling between different phonon correla-
tors .
II. GENERAL EXPRESSION FOR THE
NONLINEAR DAMPING RATE
We consider a low-frequency mode of a mesoscopic res-
onator interacting with a thermal bath. The bath ex-
citations are vibrational modes with a quasi-continuous
frequency spectrum, which we call phonons. The total
Hamiltonian
H
=
H
0
+
H
b
+
H
i
is the sum of the Hamilto-
nians of the mode
H
0
, the bath
H
b
, and their interaction
4
H
i
. We write
H
in terms of the creation and annihila-
tion 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 dis-
placement 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 vi-
brations at the mode antinode. The characteristic wave-
length 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
~
V a
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 ac-
count 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.
For weak coupling
H
i
, one can describe the dynam-
ics 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 super-
operators
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
dte
0
t
h
1
(
t
)
h
1
(0)
,
Γ
(nl)
=
~
2
[ ̄
n
(2
ω
0
) + 1]
1
Re
0
dte
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
, or to
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 super-operator
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 quasi-elastic 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 corre-
lation time of the bath excitations, i.e., to the reciprocal
time over which the correlators in Eqs. (7) and (9) de-
cay. Also, it is required that the derivatives of Γ
,
Γ
(nl)
over
ω
0
be small. These conditions are met if the den-
sity of states of the bath weighted with the appropriate
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
5
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 fol-
lows 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)
. The
term
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 change of the oscillator
frequency and the decay rate, which occur for the oscilla-
tor displacement on the order of the zero-point vibration
amplitude
q
0
, Eq. (1). In the lab frame, in the phe-
nomenological classical description of the oscillator dy-
namics, 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
(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, phe-
nomenologically, the friction force could have been writ-
ten 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 non-
linearity 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
) in Eq. (11) represents the additive
thermal Gaussian noise with the intensity
Γ
T
re-
lated 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 frequen-
cies 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 phenomenolog-
ical model, Eq. (5) does not assume that the dynamics of
the oscillator is Markovian in the lab frame; it is Marko-
vian only on the time scale that largely exceeds
ω
1
0
and
the correlation time of the thermal reservoir.
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 deter-
mined 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 scatter-
ing 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.
Parameters
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 considered mode
and another intracavity mode, 2
ω
0
ω
κ
[64],
v
κ
deter-
mines 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)
are of 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, op-
erator
v
κκ
b
κ
b
κ
can be made Hermitian,
v
κκ
=
v
κ
κ
.
In Eq. (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
6
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 con-
tinuous 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 tempera-
ture. The second term describes the decay due to scat-
tering 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
, in Eq. (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 nanome-
chanical 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 in-
elastic 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 nor-
mal modes
κ
calculated with the disorder taken into ac-
count. However, the inelastic decay rate of the modes
κ
should exceed the spacing of their frequencies. Equa-
tion (13) also applies if modes
κ,κ
are weakly decaying
bulk modes coupled to the considered mode at the res-
onator boundary.
IV. THE EIKONAL APPROXIMATION
We now consider how the nonlinear coupling between
the modes of the quasi-continuous spectrum (phonons)
affects the nonlinear damping rate Γ
(nl)
. Of primary in-
terest 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 en-
ergy
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
described in the eikonal approximation, see Appendix A.
The displacement operator of an
n
th atom (at equilib-
rium poistion
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, which smoothly depends
on
n
. Instead of being plane waves, the normal modes
u
κ
(
r
) have the form
u
κ
(
r
) =
N
1
/
2
e
κ
(
r
) exp[
iS
κ
(
r
)]
,
S
κ
(
r
) =
k
κ
(
r
)
.
(16)
One can think of Eq. (16) as the analytic 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 polar-
ization 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 approx-
imation 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).
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 “quan-
tized”, 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
coordinate
r
.
The mode number
κ
has discrete and quasi-continuous
components. For an ideal bulk system these are the
7
phonon branch
α
κ
and the wave vector
k
κ
. For thin
nanoresonators, the discrete components of
κ
also enu-
merate the transverse branches. It is important that
κ
has a quasi-continuous component and, as mentioned
above, the frequency spectrum is quasi-continuous, at
least in the range where
~
ω
κ
k
B
T
.
One can find the quasi-continuous values of
κ
by for-
mally 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 quasi-continuous part of
κ
can be
associated with the value of vector
k
κ
at some arbitrary
point
r
0
,
[
κ
]
quasi
continuous
⇐⇒
k
κ
(
r
0
)
.
(17)
Normal vibrational modes introduced in Eq.(15) satisfy
the standard orthonormality and completeness condi-
tions; 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 (15) we see that, if we disregard umklapp pro-
cesses, 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 symme-
try,
ˆ
U
is independent of
r
as are also vectors
k
κ
[48], and
then Eq. (20) becomes just the condition of the quasi-
momentum conservation in phonon scattering. The ma-
trix 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
sufficiently 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 momen-
tum uncertainty in the momentum conservation condi-
tion (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 quasi-momenta 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 lat-
tice vector. This expression describes umklapp processes
in a nonuniform medium.
V. THE TWO-PHONON CORRELATOR
It is convenient to express the nonlinear damping rate
Γ
(nl)
, Eqs. (7) and (12), in terms of the two-phonon cor-
relation function
φ
κκ
(
t
). For
~
ω
0

k
B
T
Γ
(nl)
=
2
ω
0
k
B
T
Im
κ,κ
v
κ
κ
0
dte
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 eigenfrequen-
cies,
|
ω
κ
ω
κ
| ∼
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
)
|
.
The major contribution to Γ
(nl)
comes from the modes
κ,κ
that belong to the same vibrational branch,
α
κ
=
α
κ
. Functions
φ
κκ
for different branches are fast os-
cillating 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
α
κ
6
=
α
κ
is small.
8
The correlator
φ
κκ
is the weighted off-diagonal phonon
density matrix
b
κ
(
t
)
b
κ
(
t
)
. The weighting factor is in-
dependent of time, as seen from Eq. (23). The reduction
of the problem to a calculation of the weighted thermal-
phonon density matrix significantly simplifies the analy-
sis. Our formulation allows us to study spatially nonuni-
form systems, and, in the first place, the nonlinear fric-
tion and the dephasing. Our results for linear friction
in the limit of a spatially uniform system and weak dis-
order coincide with the results based on the linearized
Boltzmann equation.
A natural approach to the analysis of the phonon dy-
namics for close
k
κ
(
r
) and
k
κ
(
r
) is to use the Wigner
transformation
W
. This is a linear integral transforma-
tion. 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
κ,κ
belong-
ing to a branch
α
. The typical value of
|
ρ
|
in Eq. (25) is
λ
T
; it is small compared to the nonuniformity length
l
sm
. We have set
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
α
; func-
tion
θ
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.
A. 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
) is a sharp
δ
-like peak that describes the ap-
proximate 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 quasi-momentum
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
α
κ
=
α
κ
, we
introduce a quasi-continuous component ̃
κ
, which is the
“center” of
κ,κ
for given
R
,
k
and is defined by equa-
tion
k
̃
κ
(
R
) =
k
,
̃
κ
̃
κ
(
R
,
k
)
.
(28)
Keeping in mind that the polarization vectors
e
κ
(
r
) in
Eq. (16) are smooth functions of the quasi-continuous
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 correla-
tion 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 sim-
ple algebra one obtains the expression for the nonlinear
9
damping rate that reads
Γ
(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 trans-
form of
φ
κκ
in a spatially nonuniform system, i.e., the
Wigner transform of the off-diagonal phonon density ma-
trix weighted with a time-independent operator. It is
much more convenient for the analysis than the tensor.
Equation (33) expresses the nonlinear damping param-
eter Γ
(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, start-
ing 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 coor-
dinate
R
and the wave vector
k
. As a function of
R
,
it varies on the distance determined by the scale of the
spatial nonuniformity and by the phonon scattering rate.
This distance is much longer than
|
k
|
1
λ
T
. The sep-
aration of the spatial scales makes it possible to derive
a transport equation for function Φ
α
. The derivation is
somewhat involved, since the standard momentum con-
servation 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
t
Φ
α
+
{
α
,
Φ
α
}
= St[Φ
α
]
,
Φ
α
Φ
α
(
R
,
k
,t
)
,
α
α
(
R
,
k
)
.
(34)
Here,
α
(
R
,
k
) =
κ
ω
κ
θ
α
(
R
,
k
,
;
κ,κ
)
(35)
is the effective position- and momentum-dependent fre-
quency 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 non-
linearity, 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) is the
locality
.
As shown in Appendix C 1, the kernel
λ
α
0
α
(
R
,
R
0
,
k
,
k
0
) becomes small for
|
R
R
0
|
exceed-
ing 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 vari-
ation 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 equa-
tion has the same form as the kinetic equation for the
Wigner transform of the phonon density matrix. The dif-
ference 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 Boltzmann equation for the
phonon distribution function. However, here it is derived
for a phonon correlation function starting from a micro-
scopic model of a nonuniform system with the nonunifor-
mity, which is smooth on phonon wavelength, but not on
the phonon mean free path.
10
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 impu-
rities, 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)
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 contin-
uous 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 de-
termine phonon scattering in the second order of the per-
turbation 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 correla-
tion 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 quasi-momentum
is conserved,
k
κ
(
r
∗∗
)
k
κ
(
r
∗∗
) +
k
κ
1
(
r
∗∗
)
k
κ
1
(
r
∗∗
) = 0
.
(41)
We note that, because of the smooth nonuniformity,
the quasi-momentum of the system of phonons is not
conserved even without short-range disorder and umk-
lapp processes. Phonons do not propagate along straight
lines.
However, Eq. (20) indicates that, locally, the
quasi-momentum 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 equa-
tion 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 scat-
tering is particularly important for lower temperatures,
where the rate of the umklapp processes becomes expo-
nentially 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 fol-
lows we illustrate the power of our approach by applying
it to spatially uniform systems and calculating nonlin-
ear 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 sys-
tem; 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 trans-
form 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
11
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 repre-
sentation 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. The 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 Ψ
ν
(
k
), respectively, and which are
independent of
r
,
St[
ψ
ν
(
k
)] =
ε
ν
ψ
ν
(
k
)
,
V
(2
π
)
d
α
d
k
Ψ
ν
(
k
)
ψ
ν
(
k
) =
δ
ν,ν
.
(44)
Here,
ε
ν
are eigenvalues of the collision operator St. On
physical grounds, since this operator describes phonon
relaxation, Re
ε
ν
0. The left eigenvectors Ψ
ν
have
the same eigenvalues as the right eigenvectors and are
orthogonal to the right eigenvectors (or can be made or-
thogonal, for degenerate
ε
ν
), the condition indicated in
Eq. (44). The eigenvector with zero eigenvalue,
ε
0
= 0,
is
ψ
0
(
k
) =
c
0
ω
k
α
̄
n
k
α
( ̄
n
k
α
+ 1) (for concreteness, in
what follows we set
c
0
=
~
); the left eigenvector is
Ψ
0
(
k
) = const
×
ω
k
α
, where the constant is determined
by the normalization condition (44).
Except for special cases, vectors
ψ
ν
(
k
) provide a
complete set, and therefore the solution of the transport
equation can be sought in the form
Φ
α
(
r
,
k
,t
) =
ν
T
ν
(
r
,t
)
ψ
ν
(
k
)
.
(45)
It follows from Eqs. (34) and (37) that functions
T
ν
sat-
isfy equation
t
T
ν
(
r
,t
) +
ν
v
νν
r
T
ν
(
r
,t
) =
ε
ν
T
ν
(
r
,t
)
,
v
νν
=
V
(2
π
)
d
α
d
k
Ψ
ν
(
k
)
v
k
α
ψ
ν
(
k
)
.
(46)
with the initial condition
T
ν
(
r
,
0) =
i
V
(2
π
)
d
~
α
d
k
Ψ
ν
(
k
)
×
V
α
(
r
,
k
) ̄
n
k
α
( ̄
n
k
α
+ 1)
.
(47)
Both linear and nonlinear decay of low-frequency
standing waves is determined by functions
T
ν
, which
correspond to even in
k
eigenfunctions
ψ
ν
(
k
), as
Φ
α
(
r
,
k
,
0) in this case is even in
k
. The separation of
eigenfunctions in even and odd in
k
is a consequence of
the symmetry of operator
̃
Λ
k
0
α
0
k
α
with respect to the si-
multaneous inversion
k
→−
k
,
k
0
→−
k
0
.
B. Thermal diffusion
We are interested in the evolution of Φ
α
(
r
,
k
,t
) on
the time scale
ω
1
0
, which largely exceeds the recip-
rocal frequencies of thermal phonon
~
/k
B
T
. Another
time scale is the characteristic relaxation time of thermal
phonons
τ
r
=
l
T
/c
s
. Yet another relevant time scale is
the relaxation time of the spatial nonuniformity, which
characterizes the decay of functions
T
ν
(
r
,t
). Of partic-
ular interest is the decay of function
T
0
(
r
,t
), since this
may be the slowest decay which determines the long-time
behavior of Φ
α
(
r
,
k
,t
).
If the decay of
T
0
(
r
,t
) is slow on the time scale
(Re
ε
ν>
0
)
1
, one can describe it in the adiabatic approx-
imation familiar from the analysis of spatial diffusion in
phonon or electron systems. Physically, it corresponds to
the picture in which there is locally (for given
r
) formed
a thermal distribution of phonons with a coordinate-
dependent temperature. Formally, in Eq. (46) one as-
sumes that functions
T
ν>
0
adiabatically follow function
T
0
, so that
T
ν
(
r
,t
)
≈−
ε
1
ν
v
ν
0
r
T
0
(
r
,t
) for
ν >
0. Then
the equation for
T
0
reads
t
T
0
(
r
,t
) =
ij
r
i
D
ij
r
j
T
0
(
r
,t
)
,
D
ij
=
ν>
0
(
v
0
ν
)
i
(
v
ν
0
))
j
ν
.
(48)
Function
T
0
(
r
,t
) can be thought of as a scaled
coordinate-dependent correction to the temperature of
high-frequency phonons. This interpretation follows from
Φ
α
(
r
,
k
,t
) being the weighted Wigner phonon density
matrix, which comes to local thermal equilibrium for a
given
T
0
(
r
,t
). In our perturbation theory
T
0
(
r
,t
) is gen-
erally complex and small in the absolute value; the oc-
cupation numbers of thermal phonons are close to their
equilibrium values ̄
n
k
α
.