Lattice vibrations: phonon dispersion in one and three dimensions

The atoms of a crystal are not frozen in place. They vibrate about their equilibrium positions, and because each atom is coupled to its neighbors by the forces of lesson 4, the vibrations travel through the crystal as waves. At long wavelengths these waves are ordinary sound. At short wavelengths, comparable to the spacing between atoms, the discreteness of the lattice changes them in important ways: their frequency stops growing with wave vector, a maximum frequency appears, and in crystals with more than one atom per cell, new optical branches emerge. Quantized, these vibrations are phonons, the particles of sound, and they carry heat, scatter electrons, limit electrical conductivity in metals and, in conventional superconductors, bind electrons into pairs. PHYS-230 treated coupled oscillators and normal modes for a few masses; this lesson extends the method to an infinite periodic lattice.

The monatomic chain

Consider a line of identical atoms of mass \(m\), equally spaced by \(a\), with each atom connected to its two neighbors by springs of constant \(C\). This is the harmonic approximation: we expand the energy of lesson 4 to second order in the displacements, so that the restoring forces are linear. Let \(u_n\) be the displacement of atom \(n\) along the chain. Newton's second law gives

Because every atom is equivalent, we try a traveling wave \(u_n = A\,e^{i(kna - \omega t)}\). Substituting and canceling the common factor gives \(-m\omega^2 = C(e^{ika} + e^{-ika} - 2) = -2C(1 - \cos ka) = -4C\sin^2(ka/2)\), so

This is the dispersion relation of the chain. For small \(k\), \(\sin(ka/2) \approx ka/2\) and \(\omega \approx a\sqrt{C/m}\,k\): the waves travel without dispersion at the speed of sound \(v_s = a\sqrt{C/m}\), as for an elastic continuum. As \(k\) increases, \(\omega\) grows more slowly than linearly, and it reaches a maximum \(\omega_{\max} = 2\sqrt{C/m}\) at \(k = \pi/a\). There the group velocity \(d\omega/dk\) vanishes: the wave is a standing wave in which neighboring atoms move in opposite directions.

A striking feature is that wave vectors differing by \(2\pi/a\) describe exactly the same motion. Replacing \(k\) by \(k + 2\pi/a\) multiplies \(u_n\) by \(e^{i2\pi n} = 1\). A wave with wavelength shorter than \(2a\) cannot be distinguished, at the positions of the atoms, from one with longer wavelength; only the displacements of the atoms matter, not what happens between them. So all distinct modes are captured by \(k\) in the range \(-\pi/a < k \le \pi/a\), which is exactly the first Brillouin zone of the one-dimensional lattice from lesson 3. For a chain of \(N\) atoms with periodic boundary conditions, \(u_{N+n} = u_n\), the allowed wave vectors are \(k = 2\pi j/(Na)\), and the zone contains exactly \(N\) of them, one mode per atom, as it must.

The diatomic chain

Now let the chain alternate between two masses, \(m_1\) and \(m_2\), with the same springs. The repeat length is now \(a\), containing two atoms. Writing \(u_n\) and \(v_n\) for the displacements of the two atoms in cell \(n\), with plane-wave amplitudes \(A\) and \(B\), the equations of motion become a \(2 \times 2\) eigenvalue problem

Setting the determinant to zero gives two solutions for each \(k\):

The lower branch, \(\omega_-\), starts at zero frequency and is linear at small \(k\): it is the acoustic branch, in which both atoms of a cell move together, like a sound wave. The upper branch, \(\omega_+\), starts at a finite frequency \(\sqrt{2C(1/m_1 + 1/m_2)}\) at \(k = 0\), where the two atoms of each cell move against each other with their center of mass fixed. In an ionic crystal such as NaCl, this motion separates positive and negative charges, creating an oscillating dipole that couples strongly to infrared light: hence the name optical branch. At the zone boundary, the branches reach \(\sqrt{2C/m_2}\) and \(\sqrt{2C/m_1}\) (taking \(m_2 > m_1\)), and between them lies a gap in which no wave can propagate.

import math, random
# Phonon dispersion of a one-dimensional chain, monatomic and diatomic.
C, m, a = 1.0, 1.0, 1.0
wmax = 2 * math.sqrt(C / m)
print("monatomic chain, omega(k) = 2 sqrt(C/m) |sin(ka/2)|")
for frac in (0.05, 0.25, 0.5, 0.75, 1.0):
    k = frac * math.pi / a
    w = wmax * abs(math.sin(k * a / 2))
    vg = wmax * a / 2 * math.cos(k * a / 2)
    print(f"  k = {frac:4.2f} pi/a: omega = {w:.4f}, phase velocity = {w / k:.4f}, group velocity = {vg:.4f}")
print("  long-wave sound speed a sqrt(C/m) =", f"{a * math.sqrt(C / m):.4f}")
print()
# Diatomic chain: masses M1, M2, spring C, cell length a (two atoms)
M1, M2 = 1.0, 3.0
def branches(k):
    s = C * (1 / M1 + 1 / M2)
    r = math.sqrt(s * s - 4 * C * C * math.sin(k * a / 2) ** 2 / (M1 * M2))
    return math.sqrt(s - r), math.sqrt(s + r)
ac, op = branches(math.pi / a)
print(f"diatomic chain M2/M1 = {M2 / M1:.0f}: at zone edge acoustic = {ac:.4f}, optical = {op:.4f}")
print(f"  expected sqrt(2C/M2) = {math.sqrt(2 * C / M2):.4f}, sqrt(2C/M1) = {math.sqrt(2 * C / M1):.4f};"
      f" optical at k = 0: sqrt(2C(1/M1+1/M2)) = {math.sqrt(2 * C * (1 / M1 + 1 / M2)):.4f}")
print(f"  frequency gap between branches: {op - ac:.4f}")
print()
# Density of states of the monatomic chain by histogram of N allowed k values
Nk = 200000
bins = [0] * 10
for j in range(Nk):
    k = -math.pi + 2 * math.pi * (j + 0.5) / Nk
    w = wmax * abs(math.sin(k / 2))
    bins[min(int(w / wmax * 10), 9)] += 1
print("monatomic DOS histogram (fraction of modes per tenth of the band):")
print("  " + " ".join(f"{b / Nk:.3f}" for b in bins))
print("  the last bin is largest: the DOS diverges as 1/sqrt(omega_max - omega) at the zone edge (a van Hove singularity)")

monatomic chain, omega(k) = 2 sqrt(C/m) |sin(ka/2)|
  k = 0.05 pi/a: omega = 0.1569, phase velocity = 0.9990, group velocity = 0.9969
  k = 0.25 pi/a: omega = 0.7654, phase velocity = 0.9745, group velocity = 0.9239
  k = 0.50 pi/a: omega = 1.4142, phase velocity = 0.9003, group velocity = 0.7071
  k = 0.75 pi/a: omega = 1.8478, phase velocity = 0.7842, group velocity = 0.3827
  k = 1.00 pi/a: omega = 2.0000, phase velocity = 0.6366, group velocity = 0.0000
  long-wave sound speed a sqrt(C/m) = 1.0000

diatomic chain M2/M1 = 3: at zone edge acoustic = 0.8165, optical = 1.4142
  expected sqrt(2C/M2) = 0.8165, sqrt(2C/M1) = 1.4142; optical at k = 0: sqrt(2C(1/M1+1/M2)) = 1.6330
  frequency gap between branches: 0.5977

monatomic DOS histogram (fraction of modes per tenth of the band):
  0.064 0.064 0.066 0.068 0.071 0.076 0.084 0.097 0.123 0.287
  the last bin is largest: the DOS diverges as 1/sqrt(omega_max - omega) at the zone edge (a van Hove singularity)

The program confirms the formulas: for the monatomic chain, the phase and group velocities agree at small \(k\) and the group velocity falls to zero at the zone boundary. For the diatomic chain with mass ratio 3, the branches at the zone edge sit at 0.8165 and 1.4142 in units of \(\sqrt{C/m_1}\), with a gap of 0.60 between them. The density of states, the number of modes per unit frequency, is computed by counting modes on a fine grid of \(k\). It rises steeply toward \(\omega_{\max}\), where the flat dispersion packs many modes into a narrow frequency range; in one dimension it diverges as \(1/\sqrt{\omega_{\max} - \omega}\). Such features, called van Hove singularities after Léon van Hove, who analyzed them in 1953, appear wherever the group velocity vanishes, and they show up in heat capacities and optical spectra.

Three dimensions

In a three-dimensional crystal with \(p\) atoms in the primitive cell, each wave vector \(\mathbf k\) in the first Brillouin zone has \(3p\) modes, because each atom can move in three directions. Three of them are acoustic branches, whose frequency vanishes linearly as \(\mathbf k \to 0\): one longitudinal (atoms move along \(\mathbf k\)) and two transverse (atoms move perpendicular to it). The remaining \(3p - 3\) are optical branches. Copper, with one atom per primitive cell, has only the three acoustic branches. Silicon and diamond, with two, have three acoustic and three optical. A crystal of \(N\) primitive cells has \(3pN\) vibrational modes in total, matching the \(3pN\) degrees of freedom of its atoms.

The general method replaces the chain's spring constant with a matrix of force constants, \(\Phi_{\alpha\beta}(\mathbf R)\), giving the force in direction \(\alpha\) on one atom when another, a lattice vector \(\mathbf R\) away, moves in direction \(\beta\). Fourier transforming gives the dynamical matrix \(D(\mathbf k)\), a \(3p \times 3p\) Hermitian matrix whose eigenvalues are \(\omega^2\) for each branch. Computing force constants from first principles, using the electronic structure methods that grew out of the band theory of lessons 8 and 9, is now routine, and the results agree closely with experiment.

Measuring dispersion: neutron scattering

How do we know phonon dispersion curves are real? The decisive technique is inelastic neutron scattering, pioneered by Bertram Brockhouse in the 1950s at Chalk River. Thermal neutrons have wavelengths comparable to atomic spacings and energies comparable to phonon energies, so a neutron scattering from a crystal can create or absorb a single phonon and change both its momentum and its energy noticeably. Conservation of energy and crystal momentum gives

where \(\mathbf q\) is the phonon wave vector. By measuring the energy and direction of neutrons scattered from a single crystal, Brockhouse built the triple-axis spectrometer and traced complete dispersion curves, point by point. His measurements of aluminum in 1955 were among the first. Today, neutron spectrometers at the Spallation Neutron Source in Tennessee and the NIST Center for Neutron Research and inelastic x-ray scattering at synchrotrons map phonons in hours.

The appearance of \(\mathbf G\) in the momentum condition reflects a subtle point: \(\hbar\mathbf q\) is called crystal momentum, but it is not the true mechanical momentum of the atoms. Because \(\mathbf q\) and \(\mathbf q + \mathbf G\) describe the same phonon, crystal momentum is conserved only modulo \(\hbar\mathbf G\). Processes in which a reciprocal lattice vector is absorbed, called umklapp processes after the German for 'flipping over', are essential to the thermal resistance of insulators in lesson 6.

Phonons: quantizing the vibrations

Each normal mode is an independent harmonic oscillator with frequency \(\omega_s(\mathbf k)\), where \(s\) labels the branch. Quantum mechanics, from PHYS-372, gives its energy levels as \((n + \tfrac12)\hbar\omega_s(\mathbf k)\), with \(n = 0, 1, 2, \dots\). We describe a mode in its \(n\)th excited state as containing \(n\) phonons, each carrying energy \(\hbar\omega\) and crystal momentum \(\hbar\mathbf k\). Phonons are bosons: any number can occupy the same mode, and in thermal equilibrium the average number is given by the Bose–Einstein distribution,

which PHYS-361 derives. Unlike photons in vacuum, phonons have a maximum frequency and a finite number of modes, a difference that dominates the low- and high-temperature behavior of the heat capacity in lesson 6. Typical phonon energies range up to a few tens of meV, or a few THz in frequency: 30 meV in copper, 165 meV for the optical phonons of diamond. The zero-point energy, \(\tfrac12\hbar\omega\) per mode, never disappears and is responsible for the expansion of light crystals such as argon and helium noted in lesson 4.