Anharmonic phonons — QHA, SCHA & the SSCHA

The harmonic dynamical matrix is the curvature of the Born–Oppenheimer surface at the $T=0$ equilibrium geometry. It says nothing about thermal expansion, phonon lifetimes, or the fate of a structure whose harmonic well is actually a shallow double well — and it predicts a temperature- independent spectrum, which is wrong wherever anharmonicity matters (soft modes, hydrides, high-symmetry phases stabilized only by entropy). This note organizes the corrections: the cubic and quartic force constants, the quasi-harmonic approximation (QHA) as the cheapest non-harmonic fix, and the genuinely self-consistent SCHA / SSCHA renormalization with its finite-temperature soft-mode stabilization. The one subtlety the formalism must get right — physical dispersion is the free-energy Hessian, not the auxiliary SSCHA force constants — is stated carefully in §4.

Prerequisites

dynamical matrix & ASR · imaginary modes

1. The anharmonic expansion

Expand the Born–Oppenheimer energy in atomic displacements $u$ (collective index over atom/cell/Cartesian) about a reference geometry:

$$V(\mathbf{u})=V_0+\frac{1}{2}\sum _{ab}{\Phi }_{ab}{u}_a{u}_b+\frac{1}{3!}\sum _{abc}{\Phi }_{abc}^{(3)}{u}_a{u}_b{u}_c+\frac{1}{4!}\sum _{abcd}{\Phi }_{abcd}^{(4)}{u}_a{u}_b{u}_c{u}_d+\cdots ,$$

with ${\Phi }_{ab}={\partial }^{2}V/\partial u_a\partial u_b$ the harmonic force constants (the dynamical matrix after mass-weighting) and the higher tensors the cubic and quartic anharmonic force constants. Translational invariance imposes an acoustic sum rule on each tensor (e.g. ${\sum }_a{\Phi }_{abc}^{(3)}=0$), the higher-order analogue of the ASR. The cubic term drives three-phonon scattering (the leading contribution to lattice thermal conductivity and to phonon linewidths via the bubble self-energy); the quartic term gives the leading frequency shift and four-phonon processes. In standard many-body perturbation theory the phonon self-energy $\Pi (\mathbf{q},\omega )={\Pi }^{(\text{loop})}+{\Pi }^{(\text{bubble})}(\omega )$ collects the quartic tadpole (real, instantaneous shift) and the cubic bubble (complex, $\omega$-dependent — its imaginary part is the inverse lifetime), and the renormalized frequencies solve $\det [{\omega }^2-{\omega }_{\mathbf{q}\nu }^2-2{\omega }_{\mathbf{q}\nu }\mathrm{Re}\Pi (\mathbf{q},\omega )]=0$. Perturbation theory fails when the harmonic ${\omega }^2$ is small or negative — exactly the soft-mode regime — which is what motivates the non-perturbative schemes below.

2. Quasi-harmonic approximation (QHA) — anharmonicity through volume only

The cheapest correction keeps the harmonic form but lets the frequencies depend on volume (more generally on the strain/lattice parameters). One computes harmonic phonons ${\omega }_{\mathbf{q}\nu }(V)$ at several volumes and writes the vibrational free energy

$$F(V,T)=V_0(V)+\frac{1}{2}\sum _{\mathbf{q}\nu }\hslash {\omega }_{\mathbf{q}\nu }(V)+k_{B}T\sum _{\mathbf{q}\nu }\ln [1-e^{-\hslash {\omega }_{\mathbf{q}\nu }(V)/k_{B}T}],$$

then minimizes over $V$ at each $T$ to get the equilibrium volume $V(T)$ — i.e. thermal expansion, with the mode Grüneisen parameters ${\gamma }_{\mathbf{q}\nu }=-\partial \ln {\omega }_{\mathbf{q}\nu }/\partial \ln V$ controlling the sign and size. The QHA captures the implicit anharmonicity (volume dependence of frequencies) but not the explicit anharmonicity (intrinsic phonon–phonon coupling at fixed volume): frequencies remain those of independent harmonic oscillators, modes have infinite lifetime, and a mode that is harmonically unstable at every accessible volume stays unstable. QHA is the right tool for thermal expansion and the equation of state of weakly anharmonic solids; it is the wrong tool for soft-mode phase transitions, where the relevant physics is exactly the explicit renormalization QHA omits.

3. Self-consistent harmonic approximation (SCHA / SSCHA)

The self-consistent harmonic approximation (Hooton; the stochastic implementation is the SSCHA of Errea–Calandra–Mauri) is a variational, non-perturbative resummation. Choose a trial harmonic density matrix ${\tilde{\rho }}_{\mathbf{\Phi },\mathcal{R}}$ parametrized by auxiliary force constants $\mathbf{\Phi }$ and centroid positions $\mathcal{R}$ (the average nuclear positions). The Gibbs–Bogoliubov inequality bounds the true free energy from above,

$$F\le \mathcal{F}[\tilde{\rho }]=F_{\text{harm}}[\tilde{\rho }]+⟨V-V_{\text{harm}}[\tilde{\rho }]{⟩}_{\tilde{\rho }},$$

and one minimizes the right-hand side over $\mathbf{\Phi }$ and $\mathcal{R}$. The expectation values $⟨\cdot {⟩}_{\tilde{\rho }}$ are Gaussian averages over the trial ensemble, evaluated by importance sampling — drawing configurations from the trial harmonic distribution and averaging the true ab-initio (or MLIP-supplied) forces. The stationarity conditions are the self-consistency equations: the centroids sit where the ensemble-averaged force vanishes, $⟨\mathbf{f}{⟩}_{\tilde{\rho }}=0$, and the auxiliary force constants satisfy

$${\Phi }_{ab}=⟨\frac{{\partial }^{2}V}{\partial u_a\partial u_b}{⟩}_{\tilde{\rho }},$$

the thermal/quantum average of the true Hessian over the trial distribution — which, because the distribution itself depends on $\mathbf{\Phi }$, must be iterated to self-consistency. This dresses the harmonic frequencies with the full anharmonicity (all orders, not truncated at quartic) and is variationally stable even when the bare harmonic $\Phi$ has imaginary modes, because the averaging is over a finite-width distribution that samples the true (stable) well. It naturally includes quantum nuclear motion (zero-point sampling) and is the workhorse for hydrides, ferroelectrics, charge-density-wave parents, and other strongly anharmonic systems.

4. The dispersion is the free-energy Hessian — not $\mathbf{\Phi }$

This is the point a careful treatment must not blur. The auxiliary force constants $\mathbf{\Phi }$ are variational parameters of the trial density matrix; they are not the physical phonon dispersion. The frequencies you would measure (the poles of the one-phonon Green’s function, the peaks of the dynamical structure factor) come from the curvature of the free energy with respect to the centroid positions — the static ($\omega \to 0$) limit of the physical self-energy:

$$\boxed{D_{ab}^{(F)}=\frac{1}{\sqrt{M_a{M}_b}}\frac{{\partial }^{2}F(\mathcal{R})}{\partial {\mathcal{R}}_a\partial {\mathcal{R}}_b}}\ne \frac{1}{\sqrt{M_a{M}_b}}{\Phi }_{ab}.$$

The free-energy Hessian differs from $\mathbf{\Phi }$ by the third- and fourth-order terms that the self-consistency does not fold into $\mathbf{\Phi }$: schematically $D^{(F)}=\Phi +{\Phi }^{(3)}\Lambda {\Phi }^{(3)}$ (the static “bubble” built from the SSCHA cubic force constants), so $D^{(F)}$ can have a soft or imaginary eigenvalue even when $\mathbf{\Phi }$ is positive-definite. Quoting the auxiliary $\mathbf{\Phi }$ frequencies as “the SSCHA dispersion” overstabilizes the structure and misses the very instability one is hunting. The static free-energy Hessian is the correct stability criterion: a second-order (displacive) phase transition is signaled by an eigenvalue of $D^{(F)}$ going through zero as a function of $T$ (or pressure), not by $\mathbf{\Phi }$ doing so. For the full spectral line shape one needs the dynamic self-energy $\Pi (\omega )$; $D^{(F)}$ is its static limit.

5. Soft-mode stabilization and the imaginary-mode link

A genuine $T=0$ harmonic instability — a negative ${\omega }^2$ whose eigenvector, when followed, lowers the energy into a lower-symmetry phase — is the displacive picture of the imaginary-mode diagnostic: the harmonic well is a double well. Two distinct things can happen, and the SSCHA distinguishes them:

• Anharmonic (entropic) stabilization. The free-energy Hessian eigenvalue is negative at $T=0$ but turns positive above a transition temperature $T_c$: the high-symmetry phase is dynamically stabilized by anharmonicity/entropy. Many perovskite cubic phases, body-centered-cubic refractory metals, and high-$T_c$ hydrides are stable only in this sense — harmonic theory wrongly condemns them as unstable, QHA cannot rescue them, and the SSCHA free-energy Hessian recovers the correct $T_c$.
• Genuine instability. The free-energy Hessian eigenvalue stays negative up to melting: the high-symmetry structure is never the (meta)stable phase, and one must distort along the eigenvector.

The distinction is precisely the “artifact vs. instability” flowchart of the imaginary-mode note, extended to finite temperature: a harmonically imaginary mode that the SSCHA stabilizes above $T_c$ is physical (and temperature-dependent), not a convergence artifact — but one must still rule out the ASR/mesh/basis artifacts of that note before attributing a soft mode to anharmonic physics.

---

Builds toward: lattice thermal conductivity (three-/four-phonon) · displacive ferroelectric transitions · conventional high-$T_c$ hydride superconductors · [[methods/imaginary-modes|finite-$T$ extension of the imaginary-mode flowchart]]

Key references

• SCHA / SSCHA — N. R. Werthamer, Phys. Rev. B 1, 572 (1970); I. Errea, M. Calandra & F. Mauri, Phys. Rev. B 89, 064302 (2014); L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea & F. Mauri, J. Phys.: Condens. Matter 33, 363001 (2021).
• Free-energy Hessian / physical dispersion — R. Bianco, I. Errea, L. Paulatto, M. Calandra & F. Mauri, Phys. Rev. B 96, 014111 (2017).
• Anharmonic phonon self-energy — A. A. Maradudin & A. E. Fein, Phys. Rev. 128, 2589 (1962).
• Quasi-harmonic / Grüneisen — see Ashcroft & Mermin, Solid State Physics, Ch. 25; G. Grimvall, Thermophysical Properties of Materials (North-Holland, 1999).