Thermal properties: phonon heat capacity and thermal conductivity

Lesson 5 showed that the vibrations of a crystal are a collection of harmonic oscillators, the phonon modes. Thermal physics follows almost immediately. The energy stored in those modes at temperature \(T\) determines the heat capacity; their anharmonic interactions produce thermal expansion; and their motion and scattering determine how well an insulating crystal conducts heat. PHYS-261 met the Einstein and Debye models in outline, and PHYS-361 derives the statistical mechanics of bosons in general. This lesson applies those tools to real crystals: it derives the Debye \(T^3\) law, compares the models with copper, and explains why diamond conducts heat five times better than copper while glass barely conducts at all.

The classical result and its failure

Classically, each vibrational mode is a harmonic oscillator with average energy \(k_BT\), half kinetic and half potential, by the equipartition theorem. A crystal of \(N\) atoms has \(3N\) modes, so its vibrational energy is \(3Nk_BT\) and its heat capacity is \(3Nk_B\), or \(3R = 24.9\) J/(mol K) per mole of atoms. This is the law of Dulong and Petit, found empirically in 1819, and it holds well for most elements at room temperature. But it fails badly at low temperatures. Measurements from the late nineteenth century onward showed that heat capacities fall toward zero as \(T \to 0\), as the third law of thermodynamics requires, and that diamond, with its stiff bonds and light atoms, falls short of \(3R\) even at room temperature.

Einstein's model

In 1907, Albert Einstein applied Planck's quantum hypothesis to the problem. He assumed that all \(3N\) modes oscillate at the same frequency \(\omega_E\). The average energy of a quantum oscillator is \(\hbar\omega_E/(e^{\hbar\omega_E/k_BT} - 1)\), plus the constant zero-point energy, so the heat capacity is

At high temperatures, \(x \to 0\) and \(C_E \to 3Nk_B\), recovering Dulong and Petit. At low temperatures, \(C_E \approx 3Nk_Bx^2e^{-x}\) falls exponentially, because the thermal energy \(k_BT\) is too small to excite even one quantum of the oscillator. This was one of the first successes of quantum theory outside radiation, and with \(\theta_E \approx 1300\) K it explained diamond's anomaly. But the exponential decay is too fast: real heat capacities fall as \(T^3\), much more slowly. The flaw is the single frequency. Real crystals have acoustic modes of arbitrarily low frequency, and those modes stay excited at temperatures where the Einstein oscillators are frozen out.

Debye's model

Peter Debye fixed this in 1912. He replaced the true phonon spectrum with that of an elastic continuum, \(\omega = v_sk\), for all three acoustic polarizations, with a single average sound speed \(v_s\). In a crystal of volume \(V\), the number of modes with wave vector smaller than \(k\) is \(3 \cdot V\,\tfrac43\pi k^3/(2\pi)^3\), so the density of states in frequency is

This grows as \(\omega^2\), as for photons. To keep the total number of modes equal to \(3N\), Debye cut the spectrum off at a maximum frequency \(\omega_D\), chosen so that \(\int_0^{\omega_D}g\,d\omega = 3N\). That gives \(\omega_D^3 = 6\pi^2v_s^3N/V\), and defines the Debye temperature \(\theta_D = \hbar\omega_D/k_B\). Effectively, the Brillouin zone is replaced by a sphere of the same volume. The energy is then

and differentiating with respect to \(T\) gives the heat capacity

At high temperatures, the integrand is approximately \(x^2\), the integral is \(\tfrac13(\theta_D/T)^3\), and \(C_D \to 3Nk_B\). At low temperatures, the upper limit can be taken to infinity, the integral equals \(4\pi^4/15\), and

This is the Debye \(T^3\) law. Its physical content is simple. At temperature \(T\), only modes with \(\hbar\omega \lesssim k_BT\) are excited; these have wave vectors up to \(k_T = k_BT/\hbar v_s\), occupying a sphere of volume proportional to \(T^3\) in \(k\)-space; each contributes about \(k_B\) to the heat capacity. Because the law depends only on the low-frequency acoustic modes, which the continuum describes exactly, it is exact at low temperatures for any crystal, not just an approximation.

import math
# Debye and Einstein heat capacities, and a kinetic estimate of thermal conductivity.
R = 8.314462618
def debye(T, th, n=400):
    if T == 0: return 0.0
    xD = th / T
    h = xD / n
    f = lambda x: x ** 4 * math.exp(x) / (math.exp(x) - 1) ** 2 if x > 0 else 0.0
    s = f(0) + f(xD) + sum((4 if i % 2 else 2) * f(i * h) for i in range(1, n))
    return 9 * R * (T / th) ** 3 * s * h / 3
def einstein(T, thE):
    x = thE / T
    return 3 * R * x * x * math.exp(x) / (math.exp(x) - 1) ** 2
thD = 343.0               # copper
thE = 0.806 * thD         # a common matching choice
print(" T (K)   Debye C (J/mol K)   Einstein C   Debye T^3 law")
for T in (5, 10, 20, 50, 100, 200, 300, 600):
    low = 12 * math.pi ** 4 / 5 * R * (T / thD) ** 3
    lows = f"{low:12.4f}" if T <= 50 else "        (n/a)"
    print(f"{T:5d}   {debye(T, thD):12.4f}   {einstein(T, thE):12.4f}   {lows}")
print(f"Dulong-Petit limit 3R = {3 * R:.3f} J/mol K")
print()
# Kinetic theory: kappa = (1/3) C v l, with C per unit volume
for name, C_mol, Vm, v, kappa in (("diamond", 6.1, 3.42e-6, 13400, 2200), ("silicon", 19.8, 12.06e-6, 6400, 148),
                                   ("quartz glass", 44.0, 27.3e-6, 4100, 1.4)):
    C_vol = C_mol / Vm
    l = 3 * kappa / (C_vol * v)
    print(f"{name:13s}: measured kappa {kappa:7.1f} W/m K implies phonon mean free path {l * 1e9:8.1f} nm")

 T (K)   Debye C (J/mol K)   Einstein C   Debye T^3 law
    5         0.0060         0.0000         0.0060
   10         0.0482         0.0000         0.0482
   20         0.3853         0.0047         0.3853
   50         4.9847         3.0510         6.0211
  100        14.7661        13.6803           (n/a)
  200        21.6283        21.3242           (n/a)
  300        23.3862        23.2507           (n/a)
  600        24.5405        24.5067           (n/a)
Dulong-Petit limit 3R = 24.943 J/mol K

diamond      : measured kappa  2200.0 W/m K implies phonon mean free path    276.1 nm
silicon      : measured kappa   148.0 W/m K implies phonon mean free path     42.3 nm
quartz glass : measured kappa     1.4 W/m K implies phonon mean free path      0.6 nm

The program integrates the Debye formula numerically with Simpson's rule. For copper, with \(\theta_D = 343\) K, the heat capacity at 300 K is 23.4 J/(mol K), about 4 parts in 100 below the measured 24.4; most of the difference is the extra heat capacity at constant pressure, which includes the work of thermal expansion, whereas the model gives the value at constant volume. It approaches 24.9 at higher temperatures. Below about 20 K the Debye curve and the \(T^3\) law agree exactly, while the Einstein model, matched to the same mean frequency, falls far below: at 20 K, 0.0047 against 0.385 J/(mol K). Above \(\theta_D/50\) or so, the \(T^3\) law overshoots, because it ignores the cutoff. Measured Debye temperatures span a wide range: 2230 K for diamond, 645 K for silicon, 428 K for aluminum, 343 K for copper, 105 K for lead and 92 K for argon. High \(\theta_D\) signals stiff bonds and light atoms.

In metals there is a second contribution, from the conduction electrons, which lesson 7 shows is linear in \(T\). Because it falls more slowly than \(T^3\), it dominates below a few kelvin. Low-temperature measurements are usually plotted as \(C/T\) against \(T^2\), which gives a straight line whose intercept measures the electronic term and whose slope gives \(\theta_D\).

How good is the Debye model?

The Debye model replaces the true phonon density of states, with its van Hove peaks and separate branches, by a smooth parabola cut off at \(\omega_D\). It works far better than such a crude replacement has any right to, for two reasons. At low temperatures, only the long-wavelength acoustic modes matter, and for those the continuum description is exact. At high temperatures, every mode contributes \(k_B\) regardless of its frequency, so only the total number of modes matters, and Debye's cutoff gets that right by construction. The model is least accurate in between, around \(\theta_D/10\) to \(\theta_D/3\), where the detailed shape of the spectrum matters.

A sensitive way to see the discrepancy is to fit the Debye formula to measured heat capacities at each temperature separately and plot the resulting effective \(\theta_D(T)\). If the model were perfect, the curve would be flat. For real crystals it typically dips by 10 to 20 parts in 100 at intermediate temperatures, revealing the extra low-frequency modes of the transverse acoustic branches, whose dispersion flattens near the zone boundary. Crystals with optical branches, such as silicon or sodium chloride, are often described more accurately by a Debye term for the acoustic branches plus an Einstein term for each optical branch, since optical branches are nearly flat and have a narrow range of frequencies. That combination is how chemists and materials scientists usually tabulate heat capacities in practice.

The Debye temperature also gives a quick estimate of other properties. Since \(k_B\theta_D = \hbar\omega_D\) and \(\omega_D \approx v_s(6\pi^2N/V)^{1/3}\), the Debye temperature can be computed from measured sound speeds; for most elements the elastic and calorimetric values agree to within about 5 parts in 100. And the Lindemann criterion, proposed in 1910, says a crystal melts when the root-mean-square vibration amplitude reaches roughly a tenth of the spacing between atoms; combined with the Debye model, it predicts that the melting point scales as \(M\theta_D^2a^2\), a rule that holds within a factor of about two across many elements.

Thermal expansion

In a perfectly harmonic crystal, the average position of each atom would not change with temperature, and the crystal would not expand. Expansion is an anharmonic effect. The real interatomic potential of lesson 4 is steeper on the compressed side than on the stretched side, so as atoms vibrate more vigorously, their average separation increases. The Grüneisen relation connects the volume expansion coefficient \(\beta\) to the heat capacity: \(\beta = \gamma C_V/(BV)\), where \(B\) is the bulk modulus and \(\gamma\), the Grüneisen parameter, measures how phonon frequencies change with volume, \(\gamma = -d\ln\omega/d\ln V\), typically between 1 and 3. So thermal expansion follows the heat capacity, rising as \(T^3\) at low temperature and leveling off above \(\theta_D\). Silicon is unusual: some of its transverse acoustic modes soften under compression, giving negative \(\gamma\), and silicon contracts slightly on warming between about 20 and 120 K.

Thermal conductivity of insulators

In an electrical insulator, heat is carried by phonons. Kinetic theory, the same argument used for gases, gives the thermal conductivity

where \(C\) is the heat capacity per unit volume, \(v\) the average phonon speed (roughly the sound speed) and \(\ell\) the phonon mean free path, the average distance a phonon travels before scattering. The program inverts this formula for three materials. Diamond, with \(\kappa \approx 2200\) W/(m K), the highest of any bulk material at room temperature, has a phonon mean free path of nearly 300 nm. Silicon's is about 40 nm. Fused quartz, a glass, has a mean free path of under 1 nm, about the spacing between atoms: in an amorphous solid, phonons with wavelengths near the atomic spacing scatter almost immediately from the disorder.

What limits \(\ell\) in a perfect crystal? In a harmonic crystal, phonons would never scatter, and \(\kappa\) would be infinite. Anharmonic forces allow phonons to collide, combining two into one or splitting one into two. But Rudolf Peierls pointed out in 1929 that collisions conserving crystal momentum, called normal processes, cannot by themselves produce thermal resistance, because they redistribute momentum among phonons without destroying the net heat current. Only umklapp processes, in which the total crystal momentum changes by a reciprocal lattice vector, reverse the flow. Umklapp requires phonons with wave vectors of order half the zone, which have energies near \(k_B\theta_D/2\); their number falls as \(e^{-\theta_D/2T}\) at low temperature, so \(\ell\) and \(\kappa\) grow rapidly on cooling.

At still lower temperatures, \(\ell\) grows until phonons scatter mainly from the surfaces of the sample or from defects and isotopes, and \(\kappa\) follows the heat capacity down as \(T^3\). So the thermal conductivity of a pure crystal rises on cooling, peaks typically between 10 and 100 K, and then falls. Isotopically purified diamond and silicon conduct heat noticeably better than natural ones, because the random mixture of isotopes, such as about 1 part in 100 of carbon-13 in natural carbon, scatters phonons. The very high conductivity of diamond, cubic boron arsenide and graphite in its planes makes them candidates for spreading heat in high-power electronics.