The dynamical matrix & the acoustic sum rule

The harmonic lattice is, in the end, one matrix. Once the Born–Oppenheimer surface $E(\{\mathbf{R}\})$ is fixed, its curvature at the equilibrium geometry — the Hessian of the surface — is all there is to the vibrational spectrum, and packaging that curvature for a periodic crystal is exactly what the dynamical matrix $D_{\alpha \beta }^{\kappa {\kappa }^{\prime }}(\mathbf{q})$ does. This note builds $D(\mathbf{q})$ from the interatomic force constants, states the eigenproblem it solves, shows that translational invariance forces the three acoustic branches to zero at $\Gamma$ through the acoustic sum rule (ASR) — and why that exact statement is only ever numerically approximate, so it has to be re-imposed by hand — and then points to the one place the picture is genuinely incomplete: the non-analytic LO–TO term in polar crystals.

Prerequisites

Born–Oppenheimer PES · Hellmann–Feynman forces · Harmonic approximation · crystal symmetry

From the PES Hessian to the force constants

Expand the BO energy to second order in displacements $u_{\alpha }(l\kappa )$ of atom $\kappa$ (in the unit cell at lattice vector ${\mathbf{R}}_l$) along Cartesian direction $\alpha$ about equilibrium. The linear term vanishes at a minimum — the Hellmann–Feynman forces are zero there — and the quadratic coefficient is the interatomic force constant, the second derivative of the surface,

$${\Phi }_{\alpha \beta }^{l\kappa ,l^{\prime }{\kappa }^{\prime }}={\frac{{\partial }^{2}E}{\partial u_{\alpha }(l\kappa )\partial u_{\beta }(l^{\prime }{\kappa }^{\prime })}|}_0=-{\frac{\partial F_{\alpha }(l\kappa )}{\partial u_{\beta }(l^{\prime }{\kappa }^{\prime })}|}_0.$$

So $\Phi$ is the rate at which the force on one atom changes when another is displaced — which is what a finite-displacement calculation literally measures, and what DFPT computes analytically by linear response. Lattice periodicity means $\Phi$ depends on the cell indices only through the difference ${\mathbf{R}}_{l^{\prime }}-{\mathbf{R}}_l$, so it is enough to fix the reference atom in the home cell ($l=0$) and tabulate ${\Phi }_{\alpha \beta }^{0\kappa ,l^{\prime }{\kappa }^{\prime }}$.

The classical equation of motion for the harmonic crystal is then $M_{\kappa }{\ddot{u}}_{\alpha }(l\kappa )=-{\sum }_{l^{\prime }{\kappa }^{\prime }\beta }{\Phi }_{\alpha \beta }^{l\kappa ,l^{\prime }{\kappa }^{\prime }}{u}_{\beta }(l^{\prime }{\kappa }^{\prime })$, a coupled system over the whole crystal. Bloch’s theorem decouples it.

The dynamical matrix and its eigenproblem

Seek normal-mode solutions of Bloch form, $u_{\alpha }(l\kappa )=(M_{\kappa }{)}^{-1/2}{e}_{\alpha \kappa }(\mathbf{q})e^{i(\mathbf{q}\cdot {\mathbf{R}}_l-\omega t)}$ — the mass weighting is the standard device that turns the generalized eigenproblem $\Phi \mathbf{e}={\omega }^{2}M\mathbf{e}$ into an ordinary, Hermitian one. Substituting collapses the crystal-wide coupling into a $3n_{\text{at}}\times 3n_{\text{at}}$ problem at each wavevector $\mathbf{q}$, governed by the dynamical matrix — the mass-weighted Fourier transform of the force constants:

$$D_{\alpha \beta }^{\kappa {\kappa }^{\prime }}(\mathbf{q})=\frac{1}{\sqrt{M_{\kappa }{M}_{{\kappa }^{\prime }}}}\sum _{l^{\prime }}{\Phi }_{\alpha \beta }^{0\kappa ,l^{\prime }{\kappa }^{\prime }}{e}^{i\mathbf{q}\cdot \left({\mathbf{R}}_{l^{\prime }}-{\mathbf{R}}_0\right)},\det \left[D(\mathbf{q})-{\omega }_{\mathbf{q}\nu }^{2}1\right]=0,$$

where $\kappa$ indexes the $n_{\text{at}}$ basis atoms, $\alpha ,\beta$ the Cartesian directions, and $l^{\prime }$ runs over lattice vectors. The eigenvalues are the squared phonon frequencies ${\omega }_{\mathbf{q}\nu }^2$ ($\nu =1,\dots ,3n_{\text{at}}$ the branch index) and the eigenvectors $e_{\alpha \kappa }(\mathbf{q}\nu )$ are the (mass-weighted) polarization patterns. Diagonalizing $D(\mathbf{q})$ on a path through the Brillouin zone is the phonon dispersion; on a grid, with the density of states, it is everything the harmonic approximation has to say.

Three structural facts are worth stating because everything downstream leans on them. (i) Because $\Phi$ is a real symmetric Hessian and obeys ${\Phi }_{\alpha \beta }^{0\kappa ,l^{\prime }{\kappa }^{\prime }}={\Phi }_{\beta \alpha }^{l^{\prime }{\kappa }^{\prime },0\kappa }$, the matrix is Hermitian, $D(\mathbf{q})=D^{†}(\mathbf{q})$, so the ${\omega }^2$ are real — but not necessarily positive. A negative eigenvalue is plotted as an imaginary frequency and flags a direction of negative curvature of the PES; whether that is a genuine instability or a numerical artifact is the subject of imaginary modes. (ii) $D(\mathbf{q})$ is periodic in reciprocal space up to a gauge, and its transformation under the space group is what makes phonons at symmetry-related $\mathbf{q}$ degenerate and what blocks the matrix by irreducible representation; Maradudin and Vosko worked out the full symmetry theory of $D(\mathbf{q})$. (iii) The whole construction is the periodic-crystal specialization of “phonons are the Hessian of the BO surface” — $D$ is that Hessian, mass-weighted and Fourier-transformed, nothing more.

The acoustic sum rule as translational invariance

The single most important exact constraint on $\Phi$ is not a symmetry of the crystal but a symmetry of space: the total energy cannot change if the entire crystal is rigidly translated. Displace every atom by the same constant vector $\mathbf{c}$, i.e. $u_{\beta }(l^{\prime }{\kappa }^{\prime })=c_{\beta }$ for all $l^{\prime }{\kappa }^{\prime }$. The force on every atom must remain zero, which forces

$$\sum _{l^{\prime }{\kappa }^{\prime }}{\Phi }_{\alpha \beta }^{0\kappa ,l^{\prime }{\kappa }^{\prime }}=0⟹{\omega }_{\Gamma ,\nu }^{\text{acoustic}}=0.$$

This is the acoustic sum rule (ASR). Read the implication directly off $D$: at $\mathbf{q}=\Gamma$ the phase factor is unity, so $D_{\alpha \beta }^{\kappa {\kappa }^{\prime }}(\Gamma )=(M_{\kappa }{M}_{{\kappa }^{\prime }}{)}^{-1/2}{\sum }_{l^{\prime }}{\Phi }_{\alpha \beta }^{0\kappa ,l^{\prime }{\kappa }^{\prime }}$, and the ASR says the mass-weighted rigid-translation vector $e_{\alpha \kappa }\propto \sqrt{M_{\kappa }}{\hat{\mathbf{c}}}_{\alpha }$ is annihilated by $D(\Gamma )$. Three orthogonal choices of $\hat{\mathbf{c}}$ give three zero eigenvalues — the three acoustic branches going to $\omega =0$ as $\mathbf{q}\to 0$. Equivalently: a uniform shift is a Goldstone mode of broken continuous translational symmetry, and the acoustic phonons are its finite-$\mathbf{q}$ dispersion, $\omega \sim c_s|\mathbf{q}|$ with $c_s$ the sound velocities.

Translational invariance is in fact only the first of the Born–Huang invariance conditions on $\Phi$ that follow from conserving total linear and angular momentum; the rotational-invariance and vanishing-stress (Huang) conditions become binding — not automatic — in low-dimensional systems, where they govern the flexural branch. I deliberately stop here and defer the full low-dimensional treatment to imaginary modes; for bulk crystals the translational ASR is the condition that matters in practice.

Why the ASR has to be re-imposed numerically

The ASR is exact for the true force constants and violated by essentially every set you actually compute. The reason is mundane and unavoidable: ${\sum }_{l^{\prime }{\kappa }^{\prime }}\Phi$ is a sum of independently evaluated second derivatives, each carrying its own numerical error, and there is no mechanism that makes the errors cancel. The usual culprits:

• Finite-displacement noise. In the frozen-phonon approach $\Phi$ comes from finite-differencing forces in displaced supercells. Residual SCF error, grid discretization, and any Pulay force left uncorrected all leak into the sum. (An uncorrected Pulay term is especially insidious because it biases forces systematically, not randomly.)
• Supercell truncation. A finite supercell truncates the range of $\Phi$; force-constant tails that fall outside the cell are simply dropped from the sum, which by itself breaks the ASR.
• Brillouin-zone / basis incompleteness. A sparse $\mathbf{q}$-grid in DFPT, or an incomplete basis, shifts individual $\Phi$ entries and unbalances the sum.

The signature is unmistakable and lives precisely at $\Gamma$: the three acoustic frequencies do not come back to zero. The residual is often a harmless few cm${}^{-1}$, but with a poor grid or basis it can be a deceptively large few THz of spurious frequency — small imaginary values that masquerade as a dynamical instability. The standard fix is to enforce the sum rule on the computed force constants, reapportioning the residual so that ${\sum }_{l^{\prime }{\kappa }^{\prime }}{\Phi }_{\alpha \beta }^{0\kappa ,l^{\prime }{\kappa }^{\prime }}=0$ holds exactly. The common recipe subtracts the self-term: for each $\kappa ,\alpha ,\beta$,

$${\Phi }_{\alpha \beta }^{0\kappa ,0\kappa }⟵{\Phi }_{\alpha \beta }^{0\kappa ,0\kappa }-\sum _{l^{\prime }{\kappa }^{\prime }}{\Phi }_{\alpha \beta }^{0\kappa ,l^{\prime }{\kappa }^{\prime }},$$

pushing the entire drift into the on-site block, which is where it physically belongs (the on-site force constant is the one quantity not independently constrained by a force measurement). In phonopy this is FC_SYMMETRY together with the explicit acoustic-sum-rule correction; Quantum ESPRESSO’s q2r.x/matdyn.x apply the analogous correction to the real-space force constants before interpolation. One caveat carried over from imaginary modes: FC_SYMMETRY imposes the translational sum rule (and permutation symmetry) and nothing more — it does not enforce the rotational-invariance or vanishing-stress conditions, so a re-imposed ASR is necessary but not, in low dimensions, sufficient.

Polar crystals: the non-analytic LO–TO term

There is one regime where the construction above is not merely numerically imperfect but physically incomplete: ionic (infrared-active) crystals, where the atoms carry nonzero Born effective charges $Z_{\kappa ,\alpha \beta }^{\ast }$. A long-wavelength longitudinal-optical (LO) phonon sets up a macroscopic polarization, hence a macroscopic electric field, and the resulting long-range dipole–dipole interaction makes the force constants not absolutely summable. The Fourier transform defining $D(\mathbf{q})$ then fails to be analytic at $\mathbf{q}=\Gamma$: the limit $\mathbf{q}\to 0$ depends on the direction $\hat{\mathbf{q}}$ of approach. The clean fix, due to Pick, Cohen and Martin and cast in DFPT form by Gonze and Lee, is to split $D$ into an analytic short-range part plus a non-analytic correction (NAC),

$$D_{\alpha \beta }^{\kappa {\kappa }^{\prime }}(\mathbf{q}\to 0)=D_{\alpha \beta }^{\kappa {\kappa }^{\prime },\text{an}}(\mathbf{q}\to 0)+\frac{1}{\sqrt{M_{\kappa }{M}_{{\kappa }^{\prime }}}}\frac{4\pi e^2}{\Omega }\frac{(\hat{\mathbf{q}}\cdot {\mathbf{Z}}_{\kappa }^{\ast }{)}_{\alpha }(\hat{\mathbf{q}}\cdot {\mathbf{Z}}_{{\kappa }^{\prime }}^{\ast }{)}_{\beta }}{\hat{\mathbf{q}}\cdot {\mathbf{\epsilon }}^{\infty }\cdot \hat{\mathbf{q}}},$$

with $\Omega$ the cell volume, ${\mathbf{Z}}_{\kappa }^{\ast }$ the Born-charge tensor of atom $\kappa$, and ${\mathbf{\epsilon }}^{\infty }$ the electronic (high-frequency, clamped-ion) dielectric tensor. The NAC is the dipole field of the displaced ions screened by the electrons. Its effect is to split the LO from the TO modes at $\Gamma$: only the longitudinal mode couples to the macroscopic field, is stiffened by it, and is pushed up to ${\omega }_{\text{LO}}>{\omega }_{\text{TO}}$ — the cubic-crystal limit being the LST relation, ${\omega }_{\text{LO}}^2/{\omega }_{\text{TO}}^2={\epsilon }^0/{\epsilon }^{\infty }$. Omit the correction and the optical branches arrive at $\Gamma$ as a single spurious LO–TO–degenerate point, the dispersion has the wrong shape near the zone center, and any phonon-derived dielectric or infrared quantity is wrong.

The ingredients ${\mathbf{Z}}^{\ast }$ and ${\mathbf{\epsilon }}^{\infty }$ are exactly what DFPT returns as the response to a homogeneous electric field; equivalently ${\mathbf{Z}}^{\ast }$ is a derivative of the macroscopic polarization, which is what the Berry-phase modern theory of polarization makes well-defined in a periodic solid in the first place. The Born charge is the bridge: it is simultaneously the force per unit field and the polarization per unit displacement, and it is the object that ties the dynamical matrix to the dielectric response of the crystal. (For 2D polar materials the $1/\Omega$ Coulomb structure changes qualitatively and the splitting must be handled with a 2D cutoff — another thread I leave to imaginary modes.)

Builds toward: finite-displacement phonons · imaginary modes · anharmonic phonons & SSCHA · electron–phonon coupling

Key references

• Foundational lattice dynamics — M. Born & K. Huang, Dynamical Theory of Crystal Lattices (Oxford University Press, 1954). The harmonic expansion, the dynamical matrix, and the invariance (sum-rule) conditions on $\Phi$.
• Symmetry of the dynamical matrix — A. A. Maradudin & S. H. Vosko, “Symmetry properties of the normal vibrations of a crystal,” Rev. Mod. Phys. 40, 1 (1968).
• Macroscopic field & LO–TO splitting — R. M. Pick, M. H. Cohen & R. M. Martin, “Microscopic theory of force constants in the adiabatic approximation,” Phys. Rev. B 1, 910 (1970).
• ASR and the non-analytic term in DFPT — X. Gonze & C. Lee, “Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory,” Phys. Rev. B 55, 10355 (1997).
• DFPT review — S. Baroni, S. de Gironcoli, A. Dal Corso & P. Giannozzi, “Phonons and related crystal properties from density-functional perturbation theory,” Rev. Mod. Phys. 73, 515 (2001).