Anharmonic Lattice Dynamics

Anharmonic lattice dynamics describes what changes when restoring forces are not linear in atomic displacement. Cubic and higher terms couple the harmonic normal modes: they can shift frequencies, permit decay, redistribute energy and alter the preferred structure. Volume-dependent harmonic frequencies also produce thermal expansion. These effects answer different questions: a frequency shift is not a decay rate, and a finite spectral width is not by itself the rate that relaxes a heat current.

We work with equilibrium crystals on an adiabatic ionic potential-energy surface. The normal modes of Phonons provide the starting coordinates; interactions among them determine how the harmonic predictions change.

Required background. Phonons supplies force constants, mass-weighted polarizations, Fourier phases and the meaning of a harmonic instability.

Helpful background. Diagrammatic Methods supplies perturbative self-energy organization, and Spectral Functions supplies the distinction between a pole and a measured peak. Boltzmann Transport supplies collision-operator reasoning; Heat Capacity and Thermodynamics supplies the measurement and ensemble distinctions used in thermal expansion.

Cubic and quartic forces

Let $i=(l,\kappa,\alpha)$ collect a cell, atom and Cartesian component, and let $u_i$ measure displacement from a stationary reference structure. At a fixed cell volume and strain, expand the Born–Oppenheimer energy as

$$\begin{aligned} U&=U_0+U_2+U_3+U_4+\cdots,\\ U_n&=\frac{1}{n!}\sum_{i_1\ldots i_n} \Phi^{(n)}_{i_1\ldots i_n} u_{i_1}\cdots u_{i_n}. \end{aligned}$$

The derivatives $\Phi^{(n)}=\partial^n U/\partial u_{i_1}\cdots\partial u_{i_n}$ are evaluated at that structure. They have units of energy divided by length to the power $n$ and are symmetric under permutation of their complete indices. The factorial belongs to the Taylor expansion; an interaction vertex defined with that factorial absorbed uses a different convention.

The harmonic Hessian is only the $n=2$ member. Cubic forces describe asymmetric local distortions and three-mode couplings; quartic forces change the curvature sampled at finite vibration amplitude and also generate four-phonon processes. A centrosymmetric crystal does not generically have zero cubic force constants: inversion transforms atom labels and displacement components together. A particular mode combination can be forbidden by its symmetry, while other cubic combinations remain allowed.

For an isolated bulk crystal, a rigid translation leaves $U$ invariant. Differentiating that identity gives higher-order acoustic sum rules, for example

$$\sum_{l\kappa} \Phi^{(3)}_{l\kappa\alpha,j,k}=0$$

for each fixed Cartesian direction $\alpha$ and remaining indices. Rotational invariance imposes further constraints. Truncating interaction range or fitting forces without these constraints can generate artificial low-frequency scattering. A positive quartic coefficient in one coordinate does not prove that a multicomponent truncated potential is bounded in every direction.

Normal-mode vertices with a fixed normalization

Use $N_c$ primitive cells with periodic boundaries and positive harmonic frequencies $\omega_\lambda$, where $\lambda=(\mathbf q,\nu)$ and $-\lambda=(-\mathbf q,\nu)$ in a paired polarization gauge. Remove the exact rigid-translation coordinates. In the full-position Fourier convention of Phonons, write

$$\begin{aligned} u_i&=\sum_\lambda d_{i\lambda}X_\lambda,\\ X_\lambda&=b_\lambda+b^\dagger_{-\lambda},\\ d_{i\lambda} &=\sqrt{\frac{\hbar}{2N_cM_\kappa\omega_\lambda}} e_{\kappa\alpha,\nu}(\mathbf q) e^{i\mathbf q\cdot(\mathbf R_l+\boldsymbol\tau_\kappa)}. \end{aligned}$$

Thus $d_{i,-\lambda}=d_{i\lambda}^*$ and $X_{-\lambda}=X_\lambda^\dagger$. Substitution gives an unambiguous finite-crystal interaction:

$$\begin{aligned} U_n&=\frac{1}{n!}\sum_{\lambda_1\ldots\lambda_n} \mathcal V^{(n)}_{\lambda_1\ldots\lambda_n} X_{\lambda_1}\cdots X_{\lambda_n},\\ \mathcal V^{(n)}_{\lambda_1\ldots\lambda_n} &=\sum_{i_1\ldots i_n} \Phi^{(n)}_{i_1\ldots i_n} \prod_{a=1}^n d_{i_a\lambda_a}. \end{aligned}$$

The vertex $\mathcal V^{(n)}$ has energy units and already contains the factors of $\hbar$, mass, frequency and $N_c$ shown above. Do not insert another normalization by cell count when using it in a transition matrix element. Translating all real-space indices together supplies a sum over $N_c$ cells. For a translationally invariant crystal, this enforces

$$\sum_{a=1}^n\mathbf q_a=\mathbf G, \qquad \mathcal V^{(n)}\propto N_c^{1-n/2},$$

with $\mathbf G$ reciprocal and the proportionality understood for allowed fixed wavevectors and a fixed bulk interaction. Cubic vertices therefore scale as $N_c^{-1/2}$ and quartic vertices as $N_c^{-1}$. Momentum sums in rates compensate these factors to give an intensive lifetime. Togo, Chaput and Tanaka (2015) give a material implementation of the force-constant construction; their placement of combinatorial factors must be translated before comparing a printed vertex with this one.

Decay channels and thermal occupation

A cubic operator contains, among other terms, one phonon annihilating into two and two phonons combining into one. For decay $1\to2+3$, energy and crystal momentum require

$$\omega_1=\omega_2+\omega_3, \qquad \mathbf q_1=\mathbf q_2+\mathbf q_3+\mathbf G.$$

For combination $1+2\to3$, change the corresponding signs. These are necessary kinematic conditions. A vanishing polarization contraction, symmetry selection rule or small force constant can suppress a channel even when its phase space is large.

The transition rate between normalized many-phonon states fixes the normalization:

$$W_{I\to F} =\frac{2\pi}{\hbar} \bigl|\langle F|U_3|I\rangle\bigr|^2 \delta(E_F-E_I).$$

Here the delta function has inverse-energy units. Replacing it by a delta function of angular frequency contributes another $1/\hbar$. Count each distinct final state once; summing ordered daughter labels instead requires their symmetry factor. This form avoids guessing an overall factor from a differently normalized three-phonon vertex.

For distinct modes, the forward decay occupation factor is $n_1(1+n_2)(1+n_3)$; its inverse has $(1+n_1)n_2n_3$. Holding the daughter bath at equilibrium, the coefficient that relaxes a small excess of $n_1$ is proportional to

$$(1+n_2)(1+n_3)-n_2n_3 =1+n_2+n_3.$$

Spontaneous decay can therefore survive at zero temperature if energy, momentum and the vertex permit it. In a combination process the corresponding difference is $n_2-n_3$, with $\omega_3=\omega_1+\omega_2$. Detailed balance cancels the full collision term in the equilibrium Bose distribution. It does not set every individual transition rate to zero.

At weak anharmonicity, the lowest cubic contribution to a dynamical shift or damping is usually second order in $\Phi^{(3)}$. The first-order quartic contraction gives a real frequency renormalization; quartic scattering can contribute a width at higher perturbative order. Static force or stress corrections also shift the equilibrium positions and must be treated consistently. Near a soft mode, calling a term “higher order” does not ensure it is small: the displacement fluctuations and small denominators can invalidate expansion about the original harmonic state.

Frequency shifts, linewidths and lifetimes

For an isolated nondegenerate branch, define the retarded coordinate propagator without a factor $1/\hbar$:

$$D^R_\lambda(t) =-i\theta(t)\langle[X_\lambda(t),X_{-\lambda}(0)]\rangle.$$

Fourier transformation with $e^{izt}$ is analytic for $\operatorname{Im}z>0$. In this convention $D^R$ has time units and may be parametrized by

$$D^R_\lambda(z) =\frac{2\omega_\lambda} {z^2-\omega_\lambda^2-2\omega_\lambda\Pi^R_\lambda(z)}.$$

Write $\Pi^R=\Delta-i\Gamma$ on the positive real-frequency axis; both $\Delta$ and $\Gamma$ have angular-frequency units. For a passive stable narrow resonance, the positive-frequency pole is $z_*\simeq\Omega-i\gamma$ with $\gamma>0$. To leading order in damping, solve

$$\Omega^2-\omega_\lambda^2 -2\omega_\lambda\Delta(\Omega)=0.$$

Linearizing the real denominator around $\Omega$ then gives

$$\begin{aligned} Z_\lambda &=\frac{\omega_\lambda} {\Omega-\omega_\lambda\partial_\omega\Delta(\Omega)},\\ \gamma&\simeq Z_\lambda\Gamma(\Omega),\\ D^R_\lambda(\omega) &\simeq\frac{Z_\lambda}{\omega-\Omega+i\gamma}. \end{aligned}$$

This local approximation assumes a simple pole, positive residue there and a self-energy varying slowly across its width. Only when the shift and derivative correction are small is $Z_\lambda\simeq1$ and $\gamma\simeq\Gamma(\omega_\lambda)$. In a degenerate or strongly mixed subspace, retain the matrix self-energy and locate its coupled poles instead of assigning a scalar correction to an arbitrary polarization basis.

The pole contribution to $-\operatorname{Im}D^R$ is Lorentzian. Its angular-frequency half width at half maximum is $\gamma$, full width is $2\gamma$, and coherent amplitude decays as $e^{-\gamma t}$. A fit using energy transfer has full width $2\hbar\gamma$.

For a single weakly damped Markov oscillator with no independent dephasing, its excess occupation decays twice as fast as its amplitude:

$$\tau_{\rm amp}=\gamma^{-1}, \qquad \tau_{\rm pop}=(2\gamma)^{-1}.$$

Additional dephasing, coupled population relaxation or a non-Lorentzian continuum breaks that conversion. Instrumental resolution, unresolved branches, static inhomogeneity and electronic damping can also broaden a measured phonon line. Use the spectral/probe forward model before identifying its width with intrinsic phonon–phonon scattering. A shifted real frequency without an imaginary self-energy supplies no lifetime at all.

A frequency shift without decay

For a single stable oscillator with $U=M\omega_0^2x^2/2+\lambda_4x^4/4$ and small positive $\lambda_4$, first-order stationary perturbation theory changes its discrete energy levels. It produces no irreversible decay: an isolated finite oscillator has discrete transitions. At finite temperature, different transition frequencies can create several unresolved lines. Treating their envelope as one Lorentzian would manufacture a lifetime. Exercise 1 calculates the shifts using the harmonic oscillator matrix elements.

Quasiharmonic expansion and explicit anharmonicity

In the quasiharmonic approximation (QHA), each fixed-volume calculation remains harmonic, but its frequencies depend on volume. For a fixed number of cells and a single stable structural branch, take

$$\begin{aligned} F(V,T)&=E_{\rm BO}(V)+\sum_\lambda f_\lambda(V,T), \\ f_\lambda &=\frac{\hbar\omega_\lambda(V)}{2} +k_{\rm B}T\ln\!\left(1-e^{-\hbar\omega_\lambda(V)/(k_{\rm B}T)}\right). \end{aligned}$$

The volume $V$ here is the total crystal volume and the mode sum contains the whole crystal. Per-cell free energies and volumes give the same relations if normalized together. At pressure $P$, minimize $F(V,T)+PV$. This construction, documented in Phonopy’s QHA formulation, includes zero-point pressure as well as thermal pressure. It requires real, stable phonon frequencies throughout the sampled branch; inserting imaginary harmonic frequencies into its Bose free energy is not a stabilization procedure.

Define the mode Grüneisen parameter at fixed reduced wavevector and branch by

$$\gamma_\lambda^{\rm G} =-\frac{\partial\ln\omega_\lambda}{\partial\ln V}.$$

The superscript distinguishes this volume derivative from the damping rate $\gamma$. Since $\partial f_\lambda/\partial\omega_\lambda=U_\lambda/\omega_\lambda$, the vibrational pressure is

$$\begin{aligned} P_{\rm vib}&=\frac1V\sum_\lambda \gamma_\lambda^{\rm G}U_\lambda,\\ U_\lambda&=\hbar\omega_\lambda(n_\lambda+\tfrac12). \end{aligned}$$

At equilibrium, let $B_T=-V(\partial P/\partial V)_T>0$ be the isothermal bulk modulus of this same free energy. Differentiating the fixed-pressure condition gives the volumetric expansion coefficient

$$\begin{aligned} \alpha_V &=\frac1V\left(\frac{\partial V}{\partial T}\right)_P =\frac1{B_T}\left(\frac{\partial P}{\partial T}\right)_V,\\ \alpha_V &=\frac{1}{VB_T}\sum_\lambda \gamma_\lambda^{\rm G}C_\lambda. \end{aligned}$$

Here $C_\lambda=(\partial U_\lambda/\partial T)_V$ is the harmonic heat capacity of that mode. The second line assumes the static energy and frequencies have no explicit temperature dependence at fixed volume; electronic, magnetic or explicit anharmonic free energies need their own derivatives. A negative weighted sum of mode Grüneisen parameters gives negative expansion without requiring negative heat capacity. Linear expansion is $\alpha_V/3$ only in the isotropic case. Anisotropic expansion requires strain derivatives and the elastic-compliance tensor.

The thermodynamic identity

$$C_P-C_V=TVB_T\alpha_V^2$$

uses heat capacities of the same sample and applies to a smooth stable phase under hydrostatic conditions. It is nonnegative even when $\alpha_V<0$. It is not permission to add an arbitrary constant to a measured heat capacity.

An observed frequency at fixed pressure changes both because the volume changes and because the interacting dynamics at a fixed volume changes:

$$\left(\frac{d\Omega}{dT}\right)_P =\left(\frac{\partial\Omega}{\partial T}\right)_V +\alpha_VV\left(\frac{\partial\Omega}{\partial V}\right)_T.$$

In QHA alone the first term vanishes and the second is $-\alpha_V\gamma_\lambda^{\rm G}\omega_\lambda$. Explicit cubic/quartic self-energies generate the fixed-volume term and intrinsic damping. Fitting a temperature-dependent frequency does not separate those mechanisms without independent volume or pressure information.

Heat-current relaxation is a collision problem

In a quasiparticle regime with slowly varying local temperature, the population contribution to heat current is

$$j^Q_i=\frac1V\sum_\lambda \hbar\omega_\lambda v_{\lambda i}\,\delta n_\lambda, \qquad \mathbf v_\lambda=\nabla_{\mathbf q}\omega_\lambda.$$

If the linearized collision operator can be replaced by independent transport relaxation times, then

$$\begin{aligned} \delta n_\lambda &=-\tau_\lambda^{\rm tr}\, \mathbf v_\lambda\cdot\nabla T\, \frac{\partial n_\lambda^{\rm B}}{\partial T},\\ \kappa_{ij}^{\rm RTA} &=\frac1V\sum_\lambda C_\lambda v_{\lambda i}v_{\lambda j}\tau_\lambda^{\rm tr}. \end{aligned}$$

This relaxation-time approximation (RTA) gives three-dimensional bulk conductivity in $\mathrm{W\,m^{-1}\,K^{-1}}$. A sheet calculation needs a declared thickness or an explicitly sheet-normalized coefficient. Frequencies, heat capacities, velocities and collision rates must belong to compatible approximations; inserting a renormalized spectrum into an inconsistent scattering kernel need not preserve detailed balance.

The central issue is what the collision operator conserves. In a normal three-phonon event, the chosen first-zone momenta obey conservation with $\mathbf G=0$. An Umklapp event has a nonzero reciprocal-vector transfer. Both can relax the occupation of one particular mode and contribute to its linewidth. Normal events nevertheless conserve the crystal pseudomomentum of the phonon distribution, while Umklapp can transfer it to the lattice.

For normal collisions alone, a small displaced Bose distribution

$$n_\lambda^{\mathbf u} =\left[ e^{\beta\hbar(\omega_\lambda-\mathbf u\cdot\mathbf q)}-1 \right]^{-1}$$

has vanishing collision integral wherever this distribution is well defined. Energy and momentum conservation cancel the forward and reverse factors. It can carry heat current even though individual phonons have finite scattering rates. Normal scattering redistributes heat and can relax components of current orthogonal to conserved momentum; it cannot relax the protected component that overlaps that conserved quantity. Therefore a linewidth inserted as $1/(2\gamma)$ into a diagonal RTA can predict resistance where the full operator has a nondecaying current component.

Umklapp, isotopic disorder, defects and boundaries remove or alter these constraints in different ways. Ward and collaborators (2009) illustrate why solving the linearized collision equation matters in a material with strong normal scattering. The geometry and competing scattering lengths also decide whether diffusive Fourier conduction is an appropriate boundary-value model. A ballistic device or hydrodynamic flow needs more than a bulk scalar mean free path.

When branches overlap strongly or mode coherences contribute substantially, even a population-only Boltzmann equation may be insufficient. Simoncelli, Marzari and Mauri (2019) develop a thermal-transport formulation including coherence between vibrational modes. This is a specific extension with its own assumptions, not a license to interpret every broad vibrational feature as a phonon with a fitted lifetime.

Soft modes and finite-temperature stability

A negative eigenvalue of the bare harmonic Hessian means that the chosen reference is unstable within that quadratic expansion. A zero eigenvalue is only inconclusive: $+x^4$ and $-x^4$ have the same vanishing Hessian at the origin but opposite local stability. Finite-temperature stabilization requires a statement about the appropriate thermodynamic free energy, not a change in how imaginary frequencies are plotted.

Self-consistent harmonic and related variational methods introduce a positive trial covariance and optimize an approximate free energy. Three objects must be distinguished:

• the frequencies of the auxiliary Gaussian distribution;
• the curvature of the optimized free energy with respect to a structural order parameter or centroid;
• the complex poles and continuum structure of a dynamical response.

Bianco and collaborators (2017) show how free-energy curvature and a dynamical extension are constructed in the self-consistent harmonic approximation. Positive auxiliary frequencies do not alone establish positive thermodynamic curvature or sharp dynamical quasiparticles.

A general variational calculation makes the distinction concrete. Let $Q$ be a centroid or order parameter and $s$ the adjustable covariance parameters, with $F_{\rm opt}(Q)=\min_s F_{\rm var}(Q,s)$. At a smooth stationary solution with invertible positive $F_{ss}$, eliminating the response of $s$ gives

$$\frac{d^2F_{\rm opt}}{dQ^2} =F_{QQ}-F_{Qs}F_{ss}^{-1}F_{sQ}.$$

Thus the curvature at fixed variational parameters can overestimate the relaxed curvature. In dimensionless local coordinates, $F_{QQ}=F_{ss}=1$ and $F_{Qs}=2$ give a relaxed curvature $-3$, despite the two positive diagonal curvatures. This is a local derivative example; $F_{QQ}$ is not being identified with the auxiliary oscillator frequency squared. A positive Gaussian covariance by itself says still less about this Schur complement.

Positive curvature, after removing exact symmetry zero modes, establishes local stability within the chosen approximation. Comparing competing free energies is a separate question, particularly at a first-order transition. Static curvature and a dynamical peak also describe different responses and need not vary together.

Connecting calculations to measurements

Electronic-force errors, displacement size, interaction-range truncation and momentum sampling affect different parts of the calculation. Numerical delta-function broadening is particularly easy to confuse with physics: a smooth curve at one broadening may conceal a sparse momentum mesh. Frequencies and rates should stabilize as the mesh is refined and the broadening reduced together. Translation and rotation constraints provide additional tests. A soft reference that invalidates perturbation theory requires a different reference or method; numerical smoothing cannot stabilize it.

Compare temperature and pressure dependence with the same structural state and probe selection rules. Inelastic neutron and x-ray measurements weight polarizations and atomic species; optical measurements select different channels. The Phonons scattering discussion and Spectral Functions explain the forward models. Thermodynamic and transport comparisons then provide independent constraints: fitting expansion does not validate a decay rate, and fitting a linewidth does not validate heat-current relaxation.

Exercises

1. A quartic shift is not a lifetime

For $H=p^2/(2M)+M\omega_0^2x^2/2+\lambda_4x^4/4$, use $x=x_0(a+a^\dagger)$ with $x_0=\sqrt{\hbar/(2M\omega_0)}$. Derive the first-order level and adjacent-transition shifts. State the smallness condition for the levels being used and whether this isolated system acquires a decay width.

Solution

The oscillator algebra gives

$$\langle n|x^4|n\rangle =3x_0^4(2n^2+2n+1).$$

Therefore

$$\begin{aligned} \delta E_n &=\frac{3\lambda_4x_0^4}{4}(2n^2+2n+1),\\ \delta\omega_{n\to n+1} &=\frac{3\lambda_4x_0^4}{\hbar}(n+1). \end{aligned}$$

The fractional transition shift scales as $\lambda_4x_0^4(n+1)/(\hbar\omega_0)$. Small mixing of the harmonic eigenstates is a stronger requirement at large $n$. For example,

$$\langle n+2|(a+a^\dagger)^4|n\rangle =(4n+6)\sqrt{(n+1)(n+2)}.$$

Comparing this off-diagonal perturbation with the gap $2\hbar\omega_0$ gives the sufficient condition $\lambda_4x_0^4(n+1)^2/(\hbar\omega_0)\ll1$ over the occupied levels. A small relative shift alone does not guarantee an accurate unperturbed eigenvector. The exact isolated oscillator still has real, discrete transition energies and no irreversible decay width. A bath, continuum limit or specified coarse-graining is additional physics.

2. Convert a linewidth with its convention intact

A resolved phonon has a pole energy near $20\,\mathrm{meV}$ and intrinsic energy FWHM $1.00\,\mathrm{meV}$. Assume a Lorentzian isolated Markov mode without additional dephasing. Find the amplitude and population lifetimes using $\hbar=0.658212\,\mathrm{meV\,ps}$. Would these data alone determine its contribution to thermal conductivity?

Solution

The energy half width is $\hbar\gamma=0.500\,\mathrm{meV}$, hence

$$\begin{aligned} \gamma&=0.759634\,\mathrm{ps^{-1}},\\ \tau_{\rm amp}&=1.316424\,\mathrm{ps},\\ \tau_{\rm pop}&=0.658212\,\mathrm{ps}. \end{aligned}$$

The half-width to pole-energy ratio is $0.025$, consistent with a narrow local resonance if no other spectral complication is present. A frequency in cycles per second would require the extra factor $2\pi$ when translated to this angular-frequency convention. The heat capacity and group velocity are still needed, and the collision operator may not reduce to this population lifetime. In particular, normal collisions can broaden the line while preserving a drifting heat-current component.

3. Negative expansion with positive heat capacity

At one temperature, two groups of modes contribute $C_1=2k_{\rm B}$ and $C_2=k_{\rm B}$ per primitive cell, with Grüneisen parameters $-2$ and $+1$, respectively. Let the cell volume be $60\,\mathrm{\mathring A^3}$, $B_T=100\,\mathrm{GPa}$ and $T=300\,\mathrm K$. Within the stated isotropic QHA, calculate $\alpha_V$ and $(C_P-C_V)/k_{\rm B}$ per cell. Which group is weighted more strongly than in an unweighted average of the Grüneisen parameters?

Solution

The heat-capacity-weighted sum is $(-2)(2k_{\rm B})+(1)(k_{\rm B})=-3k_{\rm B}$. With $B_TV_c=6.00\times10^{-18}\,\mathrm J$,

$$\alpha_V =-6.903245\times10^{-6}\,\mathrm{K^{-1}}.$$

Using the same per-cell normalization in the thermodynamic identity gives

$$\frac{C_P-C_V}{k_{\rm B}} =\frac{TB_TV_c\alpha_V^2}{k_{\rm B}} \simeq0.006213.$$

The negative-Grüneisen group carries twice the heat capacity of the other group. The weighted average is $-1$, while averaging the two listed Grüneisen parameters without their capacities gives $-1/2$. Both $C_V=3k_{\rm B}$ and $C_P$ are positive. The result describes the chosen temperature and input coefficients; it does not assert that either group has a temperature-independent capacity over an entire experimental range.

4. A finite collision rate with a stationary drift

Consider a one-dimensional first Brillouin zone and a reciprocal vector $G>0$. Suppose available branches permit energy-conserving combination events $1+2\to3$ at the following reduced wavevectors: event N has $(q_1,q_2,q_3)=(0.20,0.25,0.45)G$, and event U has $(0.35,0.40,-0.25)G$. No common linear dispersion is assumed. For $n_i=[e^{\beta\hbar(\omega_i-uq_i)}-1]^{-1}$, compare the forward and reverse occupation factors. Set $\beta\hbar uG=0.10$ for a numerical comparison, with $u$ small enough that all occupations are well defined.

Solution

Since $n_i/(1+n_i)=e^{-\beta\hbar(\omega_i-uq_i)}$ and $\omega_1+\omega_2=\omega_3$,

$$\frac{n_1n_2(1+n_3)}{(1+n_1)(1+n_2)n_3} =e^{\beta\hbar u(q_1+q_2-q_3)}.$$

For N, the exponent is zero, so forward and reverse terms cancel for any allowed small drift. Each term can be nonzero: cancellation of the collision integral is not absence of scattering. For U, the wavevector difference is $G$ and the ratio is $e^{0.10}\simeq1.10517$. A generic drift is therefore not stationary under that Umklapp channel. This calculation establishes the conservation distinction, not the absolute rate; the force-constant matrix element and phase space are still required. It explains why a mode linewidth and the eigenvalues relevant to heat-current relaxation need not coincide.

References

• R. Bianco, I. Errea, L. Paulatto, M. Calandra and F. Mauri, “Second-order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: Theory and stochastic implementation,” Physical Review B 96, 014111 (2017), doi:10.1103/PhysRevB.96.014111, author manuscript.
• Phonopy developers, “Quasi harmonic approximation,” official documentation, accessed 14 October 2026 — volume-dependent harmonic free energy and constant-pressure thermodynamics.
• M. Simoncelli, N. Marzari and F. Mauri, “Unified theory of thermal transport in crystals and glasses,” Nature Physics 15, 809–813 (2019), doi:10.1038/s41567-019-0520-x.
• A. Togo, L. Chaput and I. Tanaka, “Distributions of phonon lifetimes in Brillouin zones,” Physical Review B 91, 094306 (2015), doi:10.1103/PhysRevB.91.094306, author manuscript.
• A. Ward, D. A. Broido, D. A. Stewart and G. Deinzer, “Ab initio theory of the lattice thermal conductivity in diamond,” Physical Review B 80, 125203 (2009), doi:10.1103/PhysRevB.80.125203.