Phonons

A phonon is one quantum of a stable vibrational normal mode of a material. In a crystal, the modes are collective displacement patterns labeled by wave vector $\mathbf q$, branch $\nu$, and polarization. Their frequencies $\omega_{\nu}(\mathbf q)$ form a phonon dispersion, and their eigenvectors determine which atoms move, in which directions, and with what relative phases.

Harmonic force constants predict frequencies and atomic displacement patterns. Interactions can shift and broaden these modes, while a probe may suppress their intensity through destructive interference. Thus a calculated mode, a long-lived excitation and a visible spectral peak need not coincide.

The normal-mode quantization is developed in Phonons as Many-Body Excitations. Here we use it to calculate material dispersions, densities of states, thermodynamics and scattering intensities.

Required background. Crystals and Lattices supplies primitive cells, basis atoms, and equilibrium positions. Reciprocal Lattice supplies reciprocal vectors and the Fourier-phase language used in the dynamical matrix. Brillouin Zones supplies reduced-wave-vector zones and reciprocal-space paths used to read a phonon dispersion.

Helpful background. Phonons as Many-Body Excitations supplies normal-mode quantization and the bosonic creation and annihilation operators used here without rederivation.

From a Crystal Structure to a Spectrum

Let $\mathbf R_l$ label primitive cells and $\boldsymbol\tau_\kappa$ the equilibrium position of basis atom $\kappa$ within a cell. Its position is

$$\mathbf r_{l\kappa} = \mathbf R_l + \boldsymbol\tau_\kappa + \mathbf u_{l\kappa},$$

where $\mathbf u_{l\kappa}$ is a displacement from the chosen reference structure.

Within the adiabatic Born–Oppenheimer description, the electronic state adjusts to the ionic positions. Expand its resulting potential-energy surface:

$$V ={} V_0 + \frac{1}{2} \sum_{\substack{l\kappa\alpha\\l'\kappa'\beta}} u_{l\kappa\alpha} \Phi_{l\kappa\alpha,l'\kappa'\beta} u_{l'\kappa'\beta} + V^{(3)} + V^{(4)} + \cdots.$$

The linear term vanishes at a stationary reference structure. The harmonic interatomic force constants are

$$\Phi_{l\kappa\alpha,l'\kappa'\beta} = \frac{\partial^2V} {\partial u_{l\kappa\alpha} \partial u_{l'\kappa'\beta}} \biggr\rvert_{\mathbf u=0}.$$

They are a real-space Hessian, not merely pairwise spring constants. Electronic screening, chemical bonding, long-range electrostatics, and the chosen electronic-structure approximation are all encoded in them.

Dynamical matrix

This page uses a Fourier phase containing the full equilibrium atomic position. Translation symmetry lets one fix the first cell at $l=0$ and define

$$D_{\kappa\alpha,\kappa'\beta}(\mathbf q) ={} \frac{1} {\sqrt{M_\kappa M_{\kappa'}}} \sum_l \Phi_{0\kappa\alpha,l\kappa'\beta} \times \exp \left[ i\mathbf q\cdot \left( \mathbf R_l + \boldsymbol\tau_{\kappa'} - \boldsymbol\tau_\kappa \right) \right].$$

Some codes omit the basis positions from the Fourier phase and move the corresponding phases into the eigenvectors and structure factors. Frequencies are unchanged, but raw polarization vectors from the two conventions cannot be compared component by component without conversion.

The mass-weighted eigenproblem is

$$\sum_{\kappa'\beta} D_{\kappa\alpha,\kappa'\beta}(\mathbf q) e_{\kappa'\beta,\nu}(\mathbf q) = \omega_{\nu}^2(\mathbf q) e_{\kappa\alpha,\nu}(\mathbf q).$$

For a Hermitian dynamical matrix, the eigenvectors can be normalized as

$$\sum_{\kappa\alpha} e_{\kappa\alpha,\nu}^*(\mathbf q) e_{\kappa\alpha,\nu'}(\mathbf q) = \delta_{\nu\nu'}.$$

The eigenvalue is $\omega^2$, not $\omega$. A negative eigenvalue therefore gives exponential growth in the harmonic equations of motion. Plotting such a mode at a negative ordinate is a visualization convention for an imaginary frequency.

Translation invariance

A rigid translation cannot change the energy of a free bulk crystal. The force constants must obey the acoustic sum rule

$$\sum_{l\kappa'} \Phi_{0\kappa\alpha,l\kappa'\beta} = 0.$$

Consequently, an ordinary stable three-dimensional crystal has three zero-frequency translation modes at $\mathbf q=0$. Small nonzero values in a calculation can come from finite convergence, interpolation, or imperfect enforcement of invariance. Large violations are evidence of a defective force-constant model, not a harmless plotting issue.

Rotational invariance also matters. For an isotropic membrane with bending rigidity $\kappa_b$, areal mass density $\rho_{2D}$ and uniform tensile stress $\sigma$, the harmonic flexural dispersion obeys

$$\omega^2(q)=\frac{\kappa_b q^4+\sigma q^2}{\rho_{2D}}.$$

It is quadratic in $q$ when the membrane is tensionless, but linear at sufficiently small $q$ for $\sigma>0$. A freestanding sheet is not necessarily tensionless. Anharmonic fluctuations can further change the long-wavelength behavior; see Ochoa et al. (2012).

Quantized displacement modes

Fix the crystal’s center of mass and exclude its exact zero-frequency translations from the oscillator sums below. A rigid translation is a free coordinate, not an oscillator with a well-defined factor $1/\sqrt{\omega}$. Every remaining stable nonzero normal mode gives

$$H_{\mathrm{ph}} = \sum_{\mathbf q\nu} \hbar\omega_{\nu}(\mathbf q) \left( b_{\mathbf q\nu}^{\dagger} b_{\mathbf q\nu} + \frac{1}{2} \right).$$

The displacement operator is

$$\hat u_{l\kappa\alpha} ={} \frac{1} {\sqrt{N_cM_\kappa}} \sum_{\mathbf q\nu} \sqrt{ \frac{\hbar} {2\omega_\nu(\mathbf q)} } e_{\kappa\alpha,\nu}(\mathbf q) \times e^{i\mathbf q\cdot (\mathbf R_l+\boldsymbol\tau_\kappa)} \left( b_{\mathbf q\nu} + b_{-\mathbf q\nu}^{\dagger} \right),$$

using paired polarizations $\mathbf e_{\kappa\nu}(-\mathbf q)=\mathbf e_{\kappa\nu}(\mathbf q)^*$, with consistent choices within degenerate subspaces. This pairing makes the displacement Hermitian.

A one-phonon state is obtained by applying $b_{\mathbf q\nu}^{\dagger}$. It is not one atom oscillating, and it need not have nonzero mean displacement. A classical-looking lattice wave is generally a coherent state involving many occupation numbers.

Phonons carry energy $\hbar\omega_\nu(\mathbf q)$ and crystal momentum $\hbar\mathbf q$ modulo a reciprocal vector. Their number is not generally conserved: cubic and higher anharmonic interactions can create and annihilate phonons, and equilibrium phonons ordinarily have zero chemical potential.

Reading a Phonon Dispersion

A dispersion plot samples $\omega_\nu(\mathbf q)$ along a chosen path through the Brillouin zone. It is analogous to an electronic band plot, but the eigenvalues are squared frequencies and the eigenvectors are mass-weighted displacement patterns.

Branch counting

With $r$ basis atoms in a three-dimensional bulk crystal, there are

$$3r$$

branches at each $\mathbf q$. Near the zone center, a stable free crystal ordinarily has:

$$\begin{aligned} 3 &\quad \text{acoustic branches}, \\ 3r-3 &\quad \text{optical branches}. \end{aligned}$$

The names refer to long-wavelength displacement character, not directly to audibility or optical brightness. A monatomic primitive cell has no optical branches in the harmonic branch-counting sense. A two-dimensional monolayer embedded in three-dimensional space still has three displacement directions per atom, even though $\mathbf q$ is two-dimensional.

Acoustic branches and sound

For a regular acoustic branch at small wave vector,

$$\omega_\nu(\mathbf q) \approx v_\nu(\hat{\mathbf q}) |\mathbf q|.$$

The slope is a direction- and polarization-dependent sound velocity. In continuum elasticity, the Christoffel matrix

$$\Gamma_{ij}(\hat{\mathbf n}) = \frac{1}{\rho} C_{ikjl} n_k n_l$$

has eigenvalues $v_\nu^2(\hat{\mathbf n})$. Here $C_{ikjl}$ is the elastic tensor, $\rho$ is mass density, and $\hat{\mathbf n}$ is the propagation direction.

Longitudinal and transverse labels are exact in isotropic media and along suitable symmetry directions. At a generic wave vector in a low-symmetry crystal, a mode can have mixed polarization. The eigenvector, not the branch color chosen by plotting software, decides.

Optical branches

At a zone-center optical mode, basis atoms move against one another while the cell’s center-of-mass translation cancels. The frequency can remain finite as $\mathbf q\to0$.

An optical branch is:

• infrared active if its mode effective charge is nonzero;
• Raman active if the mode changes the polarizability tensor in an allowed symmetry channel;
• silent in a given first-order optical probe if neither selection rule is met.

These are symmetry and matrix-element statements. “Optical phonon” does not mean “visible in every optical spectrum.”

Polar crystals and directional limits

In a three-dimensional polar insulator, a longitudinal optical displacement creates a macroscopic electric field. The long-range Coulomb interaction adds a nonanalytic term as $\mathbf q\to0$. Define dimensionless Born charges by $Z^*_{\kappa,i\alpha}=(\Omega_c/e)\partial P_i/\partial u_{\kappa\alpha}$, where $\mathbf P$ is polarization, and let $\boldsymbol\epsilon_\infty$ be the relative electronic dielectric tensor. In SI units the bulk correction is

$$D_{\kappa\alpha,\kappa'\beta}^{\mathrm{NA}} (\mathbf q) ={} \frac{e^2} {\epsilon_0\Omega_c \sqrt{M_\kappa M_{\kappa'}}} \times \frac{ \left( \mathbf q\cdot\mathsf Z_\kappa^* \right)_\alpha \left( \mathbf q\cdot\mathsf Z_{\kappa'}^* \right)_\beta }{ \mathbf q\cdot \boldsymbol\epsilon_\infty \mathbf q }.$$

The limit depends on the direction of approach. It produces longitudinal–transverse optical splitting, so a frequency at $\Gamma$ may require a specified electrical boundary condition and direction. An isolated two-dimensional layer has a different Coulomb kernel; inserting a supercell volume into this bulk formula does not give its physical long-wavelength limit (Gonze and Lee, 1997).

Crossings, avoided crossings, and branch identity

Eigenvalues can cross when symmetry forbids mixing. Modes with compatible symmetry generally repel and exchange character through an avoided crossing. Sorting branches only by frequency can then create discontinuous colors or labels.

Robust mode tracking uses some combination of:

• eigenvector overlap between neighboring $\mathbf q$ points;
• irreducible representations at high-symmetry points and lines;
• projections onto atoms, directions, or molecular units;
• continuity of experimentally relevant structure factors.

At a degeneracy, individual eigenvectors may rotate arbitrarily inside the degenerate subspace. The subspace is meaningful even when a particular polarization vector is not.

Phonon Density of States

For $N_q$ wave vectors sampling a primitive Brillouin zone, define the density of states per primitive cell by

$$g(\omega) = \frac{1}{N_q} \sum_{\mathbf q\nu} \delta \left[ \omega-\omega_\nu(\mathbf q) \right].$$

In the thermodynamic limit, a three-dimensional crystal with $r$ basis atoms has

$$\int_0^\infty d\omega\, g(\omega) = 3r.$$

This normalization counts modes per cell per angular-frequency interval. For a finite periodic sample with $N_c$ cells, omitting the three rigid translations gives $3r-3/N_c$ internal modes per cell; this correction vanishes in the bulk limit. A density per physical volume is $g_V(\omega)=g(\omega)/\Omega_c$. If energy $E=\hbar\omega$ is used instead, the Jacobian must be included:

$$g_E(E) = \frac{1}{\hbar} g_\omega \left( \frac{E}{\hbar} \right).$$

Constant-frequency surfaces

In the continuum limit,

$$g(\omega) = \frac{\Omega_c}{(2\pi)^3} \sum_\nu \int_{\mathrm{BZ}} d^3q\, \delta \left[ \omega-\omega_\nu(\mathbf q) \right].$$

At a regular frequency,

$$g(\omega) = \frac{\Omega_c}{(2\pi)^3} \sum_\nu \int_{\mathcal S_{\nu,\omega}} \frac{dS} { \left| \nabla_{\mathbf q} \omega_\nu \right| }.$$

Large DOS appears where constant-frequency surfaces have large area or where group velocity is small. Critical points of $\omega_\nu(\mathbf q)$ produce van Hove features. A sharp DOS peak does not by itself imply a localized vibration or a long lifetime.

Projected and weighted densities

For mass-weighted eigenvectors, an atom-projected density can be defined as

$$g_\kappa(\omega) = \frac{1}{N_q} \sum_{\mathbf q\nu} \sum_\alpha \left| e_{\kappa\alpha,\nu}(\mathbf q) \right|^2 \delta \left[ \omega-\omega_\nu(\mathbf q) \right].$$

With the stated normalization,

$$\sum_\kappa g_\kappa(\omega) = g(\omega).$$

Other projections can resolve Cartesian direction, chemical species, layer, or molecular unit. Their meaning depends on the eigenvector convention and projector.

An experimentally extracted “phonon DOS” is often a generalized or probe-weighted density. In incoherent neutron measurements, species weights can scale roughly as scattering cross section divided by mass. Multiphonon contributions, Debye–Waller factors, detector response, and instrumental resolution must be removed or modeled before comparing with the unweighted theoretical $g(\omega)$.

Low-frequency laws

For three regular, linearly dispersing acoustic branches in a stable three-dimensional continuum, let $v_s(\hat{\mathbf n})>0$ be the acoustic phase speed in propagation direction $\hat{\mathbf n}$. Define the angular average

$$\mathcal A_v = \sum_{s=1}^{3} \frac{1}{4\pi} \int d\Omega_{\hat{\mathbf n}}\, \frac{1}{v_s(\hat{\mathbf n})^3}.$$

The radial wave-vector integral in the DOS supplies the inverse cube of the phase speed, giving the low-frequency form

$$g_V(\omega) = \frac{\omega^2}{2\pi^2} \mathcal A_v.$$

If each branch speed is direction independent, this reduces to

$$g_V(\omega) = \frac{\omega^2}{2\pi^2} \sum_{s=1}^{3} \frac{1}{v_s^3}.$$

In an isotropic solid this becomes

$$g_V(\omega) = \frac{\omega^2}{2\pi^2} \left( \frac{1}{v_L^3} + \frac{2}{v_T^3} \right).$$

The $\omega^2$ law is not universal across dimension or dispersion. A two-dimensional linear branch gives $g(\omega)\propto\omega$, while an ideal quadratic flexural branch gives a constant low-frequency contribution.

Vibrational Thermodynamics

For the nonzero harmonic modes specified above,

$$n_{\mathrm B}(\omega) = \frac{1} {e^{\beta\hbar\omega}-1}.$$

The phonon internal energy is

$$U_{\mathrm{ph}} = \sum_{\mathbf q\nu} \hbar\omega_{\nu}(\mathbf q) \left[ n_{\mathrm B} \left( \omega_\nu \right) + \frac{1}{2} \right].$$

The free energy is

$$F_{\mathrm{ph}} = \sum_{\mathbf q\nu} \biggl\{ \frac{\hbar\omega_\nu}{2} + k_{\mathrm B}T \ln \left[ 1-e^{-\beta\hbar\omega_\nu} \right] \biggr\}.$$

The mode heat capacity is

$$C_{\mathbf q\nu} = k_{\mathrm B} \frac{x_{\mathbf q\nu}^2e^{x_{\mathbf q\nu}}} {\left(e^{x_{\mathbf q\nu}}-1\right)^2}, \qquad x_{\mathbf q\nu} = \frac{\hbar\omega_{\mathbf q\nu}} {k_{\mathrm B}T}.$$

Summing over all modes gives $C_V$ in the fixed-frequency harmonic approximation.

Low and high temperature

For a three-dimensional crystal with the regular linear acoustic modes above, the leading low-temperature lattice heat capacity is

$$\frac{C_V}{V} \approx \frac{2\pi^2k_{\mathrm B}^4T^3} {15\hbar^3} \mathcal A_v.$$

For direction-independent branch speeds this becomes

$$\frac{C_V}{V} \approx \frac{2\pi^2k_{\mathrm B}^4T^3} {15\hbar^3} \sum_{s=1}^{3} \frac{1}{v_s^3}.$$

This is the Debye $T^3$ law. Optical modes and high-frequency acoustic states are exponentially suppressed at sufficiently low temperature.

At temperatures large compared with all relevant $\hbar\omega/k_{\mathrm B}$,

$$C_V \longrightarrow 3r k_{\mathrm B}$$

per primitive cell in a three-dimensional harmonic crystal. This is the Dulong–Petit limit. Electronic, magnetic, defect, and anharmonic contributions can prevent measured heat capacity from being purely phononic.

Zero-point energy affects total energies, isotope shifts, equilibrium structures, and phase competition, but it does not contribute to harmonic heat capacity when the frequencies are temperature independent.

Volume dependence and thermal expansion

The quasiharmonic approximation lets frequencies depend on volume while retaining harmonic occupation at each volume. Define the mode Grüneisen parameter

$$\gamma_{\mathbf q\nu} = - \frac{\partial \ln\omega_{\mathbf q\nu}} {\partial\ln V}.$$

Take this derivative at fixed reduced wave vector, so the same mode label is followed as the reciprocal lattice changes. For a volume-only quasiharmonic free energy in hydrostatic equilibrium, the vibrational contribution to volumetric expansion is

$$\alpha_V = \frac{1} {B_TV} \sum_{\mathbf q\nu} \gamma_{\mathbf q\nu} C_{\mathbf q\nu}.$$

Here $B_T=V(\partial^2F/\partial V^2)_T$ is the isothermal bulk modulus of the equilibrium free energy. Negative-Grüneisen modes can drive negative thermal expansion. General anisotropic expansion instead requires strain derivatives and elastic compliances; one scalar bulk modulus cannot determine every linear expansion coefficient.

Quasiharmonicity includes thermal expansion from volume-dependent harmonic spectra. It does not include every explicit phonon–phonon frequency shift or linewidth at fixed volume.

Thermal transport is not fixed by the DOS

In a relaxation-time description,

$$\kappa_{ij} = \frac{1}{V} \sum_{\mathbf q\nu} C_{\mathbf q\nu} v_{\mathbf q\nu,i} v_{\mathbf q\nu,j} \tau_{\mathbf q\nu}.$$

Here $\tau$ is an appropriate transport relaxation time, not automatically the lifetime inferred from a spectral linewidth. Momentum-conserving normal collisions can redistribute phonons without efficiently relaxing the collective heat current, so a full collision-operator calculation may be needed. The DOS supplies heat-capacity weighting; velocities, scattering and boundaries determine how that energy moves. See Anharmonicity for the collision physics.

How Scattering Sees Phonons

In a one-phonon event, the probe transfers physical momentum $\hbar\mathbf Q$ and energy $\hbar\Omega$. Thus $\mathbf Q$ is a wave-vector transfer and $\Omega$ an angular-frequency transfer. Crystal momentum conservation gives

$$\mathbf Q = \mathbf G+\mathbf q,$$

while energy conservation gives

$$\Omega = \pm\omega_\nu(\mathbf q).$$

Here transfers are defined into the sample, so positive $\Omega$ means energy lost by the probe. At fixed $\mathbf Q=\mathbf G+\mathbf q$, the plus sign creates a phonon of crystal momentum $\hbar\mathbf q$; the minus sign annihilates one of crystal momentum $-\hbar\mathbf q$. The real-force-constant harmonic model has $\omega_\nu(-\mathbf q)=\omega_\nu(\mathbf q)$ and equal equilibrium Bose occupations for these partners.

Neutron and nonresonant x-ray scattering

For coherent one-phonon neutron scattering, use the same full-position polarization convention as in the displacement operator. Expanding the scattering density $\sum_{l\kappa}b_\kappa e^{i\mathbf Q\cdot(\mathbf R_l+\boldsymbol\tau_\kappa+\mathbf u_{l\kappa})}$ to first order selects a creation amplitude with basis phase $e^{i(\mathbf Q-\mathbf q)\cdot\boldsymbol\tau_\kappa}$. Thus a representative mode structure factor is

$$F_\nu(\mathbf Q) = \sum_\kappa \frac{ b_\kappa e^{-W_\kappa(\mathbf Q)} }{ \sqrt{M_\kappa} } e^{i\mathbf G\cdot\boldsymbol\tau_\kappa} \times \mathbf Q\cdot \mathbf e_{\kappa\nu}(\mathbf q)^*.$$

$b_\kappa$ is the nuclear scattering length and $W_\kappa$ the Debye–Waller exponent. In the alternative cell-phase convention, $\mathbf e^{\mathrm{cell}}_{\kappa\nu}(\mathbf q)=e^{i\mathbf q\cdot\boldsymbol\tau_\kappa}\mathbf e_{\kappa\nu}(\mathbf q)$, the same amplitude contains $e^{i\mathbf Q\cdot\boldsymbol\tau_\kappa}\mathbf e^{\mathrm{cell}}_{\kappa\nu}(\mathbf q)^*$. Mixing the phase from one convention with eigenvectors from the other changes interference intensities. The Phonopy structure-factor documentation gives an explicit implementation with the full-position convention.

The paired annihilation channel has the same structure-factor magnitude in this harmonic equilibrium setting. Omitting overall kinematic prefactors, the one-phonon contribution has the form

$$S^{(1)}(\mathbf Q,\Omega) \propto \sum_\nu \frac{|F_\nu(\mathbf Q)|^2} {2\omega_\nu(\mathbf q)} \times\biggl\{ \left[ n_\nu+1 \right] \delta \left( \Omega-\omega_\nu \right) + n_\nu \delta \left( \Omega+\omega_\nu \right) \biggr\}.$$

This equation explains several experimental facts:

• the creation side remains at zero temperature through $n+1$;
• the annihilation side vanishes as $n\to0$;
• modes polarized perpendicular to $\mathbf Q$ can be dark through $\mathbf Q\cdot\mathbf e$;
• basis atoms can interfere constructively or destructively;
• heavier atoms have smaller displacement amplitudes for the same mass-weighted eigenvector;
• intensity changes between reciprocal zones even though $\omega_\nu(\mathbf q)$ is periodic.

Nonresonant inelastic x-ray scattering has closely related displacement and interference factors, with electronic form factors replacing neutron scattering lengths and different instrumental kinematics. Neutrons are often favorable for light atoms and isotope contrast; hard x-rays can access very small samples and extreme sample environments. Neither probe displays “the dispersion” without matrix-element and resolution effects.

Neutron Scattering and X-Ray Scattering explain the instrument kinematics, resolution and intensity normalization needed to compare these expressions with measurements.

Detailed balance requires

$$S(\mathbf Q,-\Omega) = e^{-\beta\hbar\Omega} S(-\mathbf Q,\Omega).$$

When symmetry makes the $\mathbf Q$ and $-\mathbf Q$ responses equivalent, the ratio of annihilation to creation intensity is $e^{-\beta\hbar\omega}$.

Raman, infrared, and Brillouin probes

Visible and infrared photons carry little crystal momentum compared with a Brillouin-zone dimension. First-order Raman and infrared spectra therefore probe modes near $\mathbf q=0$.

• Infrared absorption couples to a changing dipole and requires a nonzero mode effective charge.
• Raman scattering couples to a changing polarizability tensor and obeys polarization-dependent symmetry selection rules.
• Brillouin light scattering resolves long-wavelength acoustic modes and elastic properties.
• Ultrafast pump–probe measurements can launch coherent zone-center or finite-wave-vector lattice motion, depending on the excitation mechanism and geometry.

A Raman peak is not automatically an infrared peak. A zone-boundary phonon can enter optical spectra through disorder-assisted scattering, zone folding in a superstructure, or a multiphonon process with the required total momentum. Resonance can enhance a permitted process, but it does not by itself remove momentum conservation.

Frequencies, linewidths, and spectral weight

An ideal harmonic mode gives a delta-function line. Anharmonic decay, defects, boundaries and coupling to electronic or magnetic continua can broaden it. Hybridization between two closed harmonic modes merely produces two new sharp modes; a broad unresolved feature or damping into a continuum requires an additional explanation. Instrumental resolution also broadens the measured spectrum.

For a pole at $\omega_0-i\Gamma$ with $\Gamma\ll\omega_0$, the amplitude decays as $e^{-\Gamma t}$ and the Lorentzian angular-frequency full width is $2\Gamma$. Its energy population decays with time constant $1/(2\Gamma)$ in the simple damped-oscillator model. State the convention before converting a width into a lifetime. Spectral Functions explains the pole approximation and its limitations; Structure Factors develops detailed balance for general correlations.

Computing Phonons in Materials

Two common first-principles routes are:

1. finite displacement: displace atoms in a supercell, calculate forces, infer real-space force constants, and Fourier interpolate;
2. density-functional perturbation theory: calculate the linear response to periodic displacements directly at selected $\mathbf q$ points.

Both routes approximate the electronic ground-state energy surface and require convergence tests. Agreement between software packages is not meaningful unless structure, pseudopotentials, exchange-correlation approximation, magnetic order, spin–orbit treatment, supercell or $\mathbf q$ mesh, and nonanalytic corrections are aligned.

The important convergence questions follow the calculation itself: are the structure and electronic forces converged, do the force constants satisfy translation and rotation constraints, and are the supercell and full-zone mesh large enough? Polar materials also require a consistent long-range correction. Use a dense mesh for thermodynamics and a line path for visualization; these serve different purposes.

Finite differences need a window

Displacements that are too small amplify force noise. Displacements that are too large mix in anharmonicity. A reliable finite-displacement calculation repeats the extraction over a reasonable amplitude range and checks that the force constants and frequencies are stable.

A path is not the Brillouin zone

A standard high-symmetry path can miss:

• an instability at a generic wave vector;
• a narrow avoided crossing;
• anisotropic sound velocity away from displayed directions;
• a DOS feature generated by a broad off-path surface.

Stability and thermodynamics require a converged full-zone mesh, not only an attractive line plot.

Harmonic stability is conditional

If every harmonic $\omega^2$ is nonnegative, there is no exponentially growing harmonic mode. This harmonic test alone does not establish a local minimum: zero modes not enforced by exact symmetries require higher-order analysis. For a finite smooth potential at a stationary reference structure, a Hessian that is positive definite after removing exact symmetry directions is a sufficient test for a strict local minimum modulo those symmetries. The thermodynamic or continuum limit needs separate control. Even a local minimum does not prove global thermodynamic stability; another structure may have lower free energy.

Conversely, an imaginary zero-temperature harmonic mode can be stabilized at finite temperature by anharmonic fluctuations, disorder, or entropy. The correct response is to examine the eigenvector, convergence, reference symmetry, and temperature-dependent effective potential rather than deleting the point or replacing $\omega$ by $|\omega|$.

A soft mode identifies a low-cost distortion and its wave vector identifies the translation symmetry it would break. Electron screening can produce a Kohn anomaly, but the full electron–phonon matrix elements and susceptibility matter, not nesting alone. Electron–Phonon Coupling develops that connection.

A diatomic chain: branches and state counting

Consider a one-dimensional chain with cell length $a$, masses $m$ and $M$, and identical nearest-neighbor springs $K$. Each atom has one scalar displacement. In the full-position convention with basis separation $a/2$,

$$D(q)=\begin{pmatrix} 2K/m&-2K\cos(qa/2)/\sqrt{mM}\\ -2K\cos(qa/2)/\sqrt{mM}&2K/M \end{pmatrix}.$$

Solving this two-dimensional eigenproblem gives

$$\begin{aligned} \omega_\pm^2(q)&=A\pm\sqrt{A^2-B\sin^2(qa/2)},\\ A&=K\left(\frac1m+\frac1M\right),\qquad B=\frac{4K^2}{mM}. \end{aligned}$$

The minus branch is acoustic and the plus branch optical. For $M>m$, the acoustic branch ends at $\sqrt{2K/M}$ and the optical branch begins at $\sqrt{2K/m}$ at the zone boundary. Their gap closes when the masses become equal. In that limit the actual primitive translation is $a/2$; keeping a two-atom cell folds one monatomic branch into two.

The state count is especially transparent because each branch is monotonic on $0\le q\le\pi/a$. Both signs of $q$ contribute, so

$$g_\pm(\omega)=\frac{a}{\pi}\left|\frac{dq}{d\omega}\right|, \qquad \int g_\pm(\omega)\,d\omega=1.$$

Thus slow group velocity produces a large DOS without adding states. No modes occur between the two bands. A path through a higher-dimensional crystal cannot generally supply its DOS this way: the full-zone integral is essential.

The chain with $M=2m$, in units $\omega_0=\sqrt{K/m}$. The lower panel counts the same modes in bins of width $0.02\omega_0$; its total area is two modes per cell. The empty interval $\omega_0<\omega<\sqrt2\omega_0$ is a real gap. Finite bins round the integrable edge singularities and slightly spread the outermost occupied bins; they do not represent a lifetime or instrumental width.

Exercises

1. Acoustic motion in the diatomic chain

Find the normalized mass-weighted translation eigenvector at $q=0$. Derive the acoustic sound speed and the two branches’ total number of states per cell.

Solution

A rigid physical displacement has mass-weighted components

$$\mathbf e_-(0)=\frac{(\sqrt m,\sqrt M)^T}{\sqrt{m+M}}.$$

Multiplying by $D(0)$ gives zero. Expanding the smaller eigenvalue gives

$$\omega_-^2\simeq\frac{K a^2}{2(m+M)}q^2, \qquad v_s=a\sqrt{\frac{K}{2(m+M)}}.$$

The mass density is $(m+M)/a$ and two equal springs in series give the same long-wavelength elastic result. Integrating $(a/\pi)|dq/d\omega|$ along either monotonic branch yields $(a/\pi)(\pi/a)=1$. There are two internal displacement coordinates per cell and two modes in the bulk count.

2. Which sound branch dominates low-temperature heat capacity?

An isotropic solid has $v_L=2v_T$. Compare the longitudinal contribution to the combined contribution of the two transverse branches in the Debye $T^3$ coefficient.

Solution

Each branch contributes in proportion to $v^{-3}$. The longitudinal weight is

$$\frac{1}{v_L^3} = \frac{1}{8v_T^3}.$$

The two transverse branches contribute

$$\frac{2}{v_T^3}.$$

Their ratio is

$$\frac{ 2/v_T^3 }{ 1/(8v_T^3) } = 16.$$

The slower transverse branches dominate because low-frequency state counting scales as $v^{-3}$.

3. A dark optical mode

Two identical atoms in a primitive cell have full-position polarization vectors with equal magnitudes and opposite directions. Choose $\mathbf Q=\mathbf G+\mathbf q$ such that $e^{i\mathbf G\cdot\boldsymbol\tau_\kappa}$ is equal for the two atoms, and assume equal $b_\kappa e^{-W_\kappa}/\sqrt{M_\kappa}$. What happens to the one-phonon structure factor?

Solution

The two terms are equal in magnitude. Their opposite polarizations give opposite signs in $\mathbf Q\cdot\mathbf e_{\kappa\nu}^*$, while the stated basis phases and species weights are equal. Therefore they cancel:

$$F_\nu(\mathbf Q) = 0.$$

The mode is dark in that geometry even though it exists. Moving to a reciprocal zone with a different basis phase, changing polarization geometry, or using another probe can restore intensity.

4. A finite-displacement error window

For $V(u)=ku^2/2+\lambda u^4/4$, estimate $k$ using $k_{\mathrm{FD}}=-[F(h)-F(-h)]/(2h)$. Each computed force has independent zero-mean error with variance $\sigma_F^2$. Find the systematic bias, noise variance and the displacement that minimizes mean-square error when $\lambda\ne0$.

Solution

Since $F(u)=-ku-\lambda u^3$, the central difference gives $k_{\mathrm{FD}}=k+\lambda h^2$ before force noise. The bias grows with $h$, whereas subtracting two noisy forces and dividing by $2h$ gives variance $\sigma_F^2/(2h^2)$. Hence

$$\begin{aligned} \operatorname{MSE}(h)&=\lambda^2h^4+\frac{\sigma_F^2}{2h^2},\\ h_{\mathrm{opt}}&=\left(\frac{\sigma_F^2}{4\lambda^2}\right)^{1/6}. \end{aligned}$$

This toy model explains why arbitrarily small displacements are not best. Real force errors may be correlated or systematic, so a stable interval under changes of displacement and electronic convergence is more informative than this one optimum.

References

• N. W. Ashcroft and N. D. Mermin, Solid State Physics, Holt, Rinehart and Winston, 1976.
• A. Q. R. Baron, “Introduction to High-Resolution Inelastic X-Ray Scattering,” in Synchrotron Light Sources and Free-Electron Lasers, Springer, 2016, doi:10.1007/978-3-319-14394-1_52; updated author version, arXiv:1504.01098v7, 2020.
• S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, “Phonons and Related Crystal Properties from Density-Functional Perturbation Theory,” Reviews of Modern Physics 73, 515–562 (2001), doi:10.1103/RevModPhys.73.515.
• M. Born and K. Huang, Dynamical Theory of Crystal Lattices, Oxford University Press, 1954.
• W. Cochran, “Crystal Stability and the Theory of Ferroelectricity,” Advances in Physics 9, 387–423 (1960), doi:10.1080/00018736000101229.
• M. T. Dove, Introduction to Lattice Dynamics, Cambridge University Press, 1993, doi:10.1017/CBO9780511619885.
• B. Fultz, “Vibrational Thermodynamics of Materials,” Progress in Materials Science 55, 247–352 (2010), doi:10.1016/j.pmatsci.2009.05.002.
• F. Giustino, “Electron–Phonon Interactions from First Principles,” Reviews of Modern Physics 89, 015003 (2017), doi:10.1103/RevModPhys.89.015003.
• X. Gonze and C. Lee, “Dynamical Matrices, Born Effective Charges, Dielectric Permittivity Tensors, and Interatomic Force Constants from Density-Functional Perturbation Theory,” Physical Review B 55, 10355–10368 (1997), doi:10.1103/PhysRevB.55.10355.
• H. Ochoa, E. V. Castro, M. I. Katsnelson, and F. Guinea, “Scattering by Flexural Phonons in Suspended Graphene under Back Gate Induced Strain,” Physica E 44, 963–966 (2012), doi:10.1016/j.physe.2011.03.017; author version.
• Phonopy developers, “Dynamic structure factor”, official documentation, accessed 2026-10-14. An explicit one-phonon scattering implementation with full-position dynamical-matrix phases and paired creation/annihilation conventions.
• G. L. Squires, Introduction to the Theory of Thermal Neutron Scattering, 3rd ed., Cambridge University Press, 2012.
• A. Togo and I. Tanaka, “First Principles Phonon Calculations in Materials Science,” Scripta Materialia 108, 1–5 (2015), doi:10.1016/j.scriptamat.2015.07.021.