Specific Heat Arises from the Uncertainty of Energy

When we put heat into a material, where does it go? It is stored in the vibrations of the atoms. But if we calculate this with classical mechanics, the constant-volume specific heat per atom turns out to be 3k_B at any temperature — the Dulong–Petit law. Experiment does not agree. As T → 0, C_V drops to zero. The key is the uncertainty of energy. In fact, from thermodynamics and statistical mechanics alone C_V = (ΔE)²/k_BT² — the specific heat is nothing other than how much the energy of the substance can fluctuate — follows (Section 3). All that remains is that quantum mechanics says “energy can be exchanged only in quanta of ħω,” so at low temperatures the system sticks to its lowest level, can no longer fluctuate, and the specific heat vanishes. That is all. Here we follow this single path in the order classical theory → fluctuations → Einstein → Debye.

0. What Is Wrong?

Let us first list what is known from experiment.

For many solids near room temperature, the molar heat capacity clusters around 25 J/(mol K) (Dulong–Petit, 1819).
Yet diamond has only about 6 J/(mol K) even at room temperature.
Every solid has C_V → 0 as T → 0, and in insulators the decrease follows C_V ∝ T³.

Classical theory can explain only the first of these three. Worse, it does not merely fail to explain the second and third facts: it flatly contradicts them. To find where it went wrong, let us first carry the classical calculation through to the end.

1. Classical Theory — Energy Equipartition and Dulong–Petit

1.1 One atom is a three-dimensional harmonic oscillator

The atoms in a crystal vibrate slightly around their equilibrium positions. At the equilibrium position the slope of the potential is 0, so if we take it as the origin and Taylor-expand, the first-order term vanishes and the second-order term plays the leading role. In other words, one atom is a harmonic oscillator in each of the three directions x, y, z.

H = \sum_{α = x, y, z} [\frac{p_{α}^{2}}{2 m} + \frac{1}{2} m ω^{2} x_{α}^{2}]

What matters is that this is a quadratic form in both p and x. Counting the degrees of freedom, there are three from the kinetic energy and three from the potential energy, six in all.

1.2 Computing the classical partition function

In classical statistical mechanics, the partition function is built from an integral over phase space. We write β = 1/k_BT.

Z = \frac{1}{h^{3}} \int e^{- β H} d^{3} p d^{3} x = {(\sqrt{\frac{2 π m}{β}} \cdot \frac{1}{h} \sqrt{\frac{2 π}{β m ω^{2}}})}^{3} = {(\frac{1}{β ħ ω})}^{3}

That was just six Gaussian integrals. The mean energy follows from ⟨E⟩ = −∂(ln Z)/∂β.

⟨ E ⟩ = - \frac{\partial}{\partial β} (- 3 \ln β + const.) = \frac{3}{β} = 3 k_{B} T

This is the equipartition theorem of energy — each degree of freedom that can be written as a quadratic form carries (1/2)k_BT. There are six of them, so 3k_BT. Differentiating with respect to temperature gives the specific heat.

(1)

C_{V} = \frac{\partial ⟨ E ⟩}{\partial T} = 3 k_{B} (per atom)

Per mole this becomes C_V = 3R = 3 × 8.3145 = 24.94 J/(mol K). This is the Dulong–Petit law, and it does reproduce well the room-temperature values of many metals — Pb, Cu, Al, Ag, Au, ….

⚠ But this result is suspicious in three ways ① The information about the material has disappeared. Both m and ω cancel along the way and do not survive. Whether the spring is stiff or soft, and whether the atom is heavy or light, the answer is the same 3k_B.
② The temperature has disappeared. It is 3k_B at 1 mK and at 3000 K alike.
③ It contradicts the third law of thermodynamics. Since S(T) = ∫₀^T C_VdT′/T′, if C_V remained finite as T → 0, the entropy would diverge logarithmically. For S to approach a finite constant value as T → 0, C_V → 0 is required.

2. Energy Takes Discrete Values in Steps of ħω

This section gives only the essentials — see the harmonic oscillator page for details The results E₀ = ħω/2 and E_n = (n+1/2)ħω that appear here are treated head-on in How the Harmonic Oscillator Is Solved, Why It Cannot Come to Rest, and Second Quantization. Where the two overlap, that page is far more detailed.

Section 2, Solving the equation — the place where the Schrödinger equation is solved by a series and, from the condition that the series must terminate at a finite order, quantization comes out
Section 4, What the zero-point energy really is — the minimization of (Δp)² below, done more carefully
Section 5, Absolute zero and the third law of thermodynamics — the story that “absolute zero = standing still” is not true. It also connects directly with Section 3 (fluctuations) of this page
Section 10, Climbing up and down the ladder — another route to the same result using creation and annihilation operators

Here we restate only the minimum needed for the discussion of the specific heat.

2.1 An atom cannot “sit still at a fixed position”

To pin down what was wrong with classical theory, let us go back to a single oscillator. Position and momentum are subject to the restriction

Δ x Δ p \geq \frac{ħ}{2}

Now consider the state of sitting quietly at the equilibrium position (x = 0, p = 0): then Δx = Δp = 0, which violates this restriction. An atom cannot stand still. So how quiet can it get? Writing the expectation value of the energy in terms of Δx and Δp, we have

⟨ E ⟩ = \frac{(Δ p)^{2}}{2 m} + \frac{1}{2} m ω^{2} (Δ x)^{2} \geq \frac{(Δ p)^{2}}{2 m} + \frac{m ω^{2} ħ^{2}}{8 (Δ p)^{2}}

Minimizing the right-hand side with respect to Δp (the AM–GM inequality gives the same result), the minimum is attained at (Δp)² = mωħ/2, where its value is

E_{0} = \frac{1}{2} ħ ω

This is the zero-point energy. From the uncertainty principle alone, it follows that the vibration cannot stop. If we solve the Schrödinger equation properly, the allowed energies are

(2)

E_{n} = (n + \frac{1}{2}) ħ ω (n = 0, 1, 2, \dots)

That is, they take equally spaced discrete values. That the levels are equally spaced does not itself follow from the uncertainty relation — all the uncertainty relation tells us is that the lowest level cannot go below ħω/2. The route to the same result using creation and annihilation operators is in Section 10 of the harmonic oscillator page.

This is the most important point of today Energy can go in and out only in quanta of ħω. Intermediate values are not allowed, so to give energy to this mode you must hand over at least one whole ħω. The energy that a heat bath can supply is of the order of k_BT, so

ħω ≪ k_BT … the quanta are fine-grained, so the energy is effectively continuous. Same as classical theory.
ħω ≫ k_BT … cannot afford even a single quantum. This mode falls asleep.

The specific heat is “how much energy the system can take up when the temperature is raised,” so a sleeping mode makes no contribution to the specific heat at all. It is also clear at the same time that letting ħ → 0 makes the steps disappear and brings us back to classical theory — the deviation from Dulong–Petit is a direct manifestation of the finiteness of ħ.

3. The Specific Heat Itself Is the Fluctuation of Energy

Let us take a small detour here. Or rather, this section is the real content of the title of this page. Neither Debye nor Einstein has appeared yet, and we do not even need quantum mechanics. From thermodynamics and statistical mechanics alone,

C_{V} = \frac{{(Δ E)}^{2}}{k_{B} T^{2}}

ΔE is the uncertainty (standard deviation) of the energy.

the relation above emerges. The specific heat is nothing other than how much the energy of the substance can “fluctuate” — once you see this, the rest of the story becomes much clearer. What cannot fluctuate cannot store heat.

3.1 First, let us define “fluctuation”

In probability theory, the variance σ² and the standard deviation σ of an observable x are defined by σ² = (1/N)Σ_i(x_i − ⟨x⟩)². It is exactly the same in quantum mechanics; for momentum, for example,

{(Δ p)}^{2} = \int ψ^{*} {(\hat{p} - ⟨ p ⟩)}^{2} ψ d x = ⟨ p^{2} ⟩ - {⟨ p ⟩}^{2}

Expanding gives ⟨p²⟩ − 2⟨p⟩⟨p⟩ + ⟨p⟩². In the last term we have used ∫ψ*ψdx = 1 (the normalization condition).

The same holds for position and for energy: (ΔE)² = ⟨E²⟩ − ⟨E⟩². The ΔE here is a different thing from the Heisenberg uncertainty relation of Section 2 — it is a statistical width: how much the energy of a system in thermal equilibrium scatters among the values it can take. What we are going to show is that this ΔE is the specific heat itself.

3.2 Thermodynamics: C_V = −T(∂²F/∂T²)_V

The constant-volume specific heat is defined by C_V ≡ (∂U/∂T)_V. Let us differentiate the Helmholtz free energy F ≡ U − TS with respect to temperature at constant volume.

{(\frac{\partial F}{\partial T})}_{V} = {(\frac{\partial U}{\partial T})}_{V} - T {(\frac{\partial S}{\partial T})}_{V} - S = C_{V} - T {(\frac{\partial S}{\partial T})}_{V} - S

But from a thermodynamic relation, (∂F/∂T)_V = −S. Substituting this into the left-hand side, S cancels on both sides and we get

(3)

- S = C_{V} - T {(\frac{\partial S}{\partial T})}_{V} - S ⟹ C_{V} = T {(\frac{\partial S}{\partial T})}_{V}

What it means for heat to be stored in matter Eq. (3) can be read as “at constant volume, raising the temperature increases the entropy, and that is why heat accumulates inside the substance.” Whether heat can be banked somewhere depends on how much the number of accessible states increases when the temperature is raised — that is one way to put it. Further on, we will see that this is the same thing as the “fluctuation of energy.”

Using S = −(∂F/∂T)_V once more, we can write it with only the second derivative of F.

(4)

C_{V} = T {(\frac{\partial S}{\partial T})}_{V} = T {(\frac{\partial}{\partial T})}_{V} {(\frac{- \partial F}{\partial T})}_{V} = - T {(\frac{\partial^{2} F}{\partial T^{2}})}_{V}

3.3 Statistical mechanics: differentiating F twice via the partition function

In statistical thermodynamics, we define the partition function Z = Σ_i e^−E_i/k_BT and from it define the internal energy and the free energy by

U \equiv \sum_{i} E_{i} e^{- E_{i} / k_{B} T} / Z, F \equiv - k_{B} T \ln Z

Now the first differentiation. Differentiating the T inside Z gives ∂Z/∂T = Σ(E_i/k_BT²)e^−E_i/k_BT and therefore

{(\frac{\partial F}{\partial T})}_{V} = - k_{B} \ln Z - k_{B} T \frac{1}{Z} \sum_{i} \frac{E_{i}}{k_{B} T^{2}} e^{- E_{i} / k_{B} T} = - k_{B} \ln Z - \frac{U}{T}

The second term is exactly of the form U/T.

Now differentiate once more. Here a pleasing cancellation occurs — differentiating the first term (−k_BlnZ) produces −U/T², and differentiating the second term (−U/T) produces +U/T², so these two cancel. What remains is only the part from differentiating U with respect to T.

{(\frac{\partial^{2} F}{\partial T^{2}})}_{V} = - \frac{1}{T} \frac{1}{k_{B} T^{2}} [\sum_{i} E_{i}^{2} e^{- E_{i} / k_{B} T} / Z - {(\sum_{i} E_{i} e^{- E_{i} / k_{B} T} / Z)}^{2}] = - \frac{1}{T} \frac{1}{k_{B} T^{2}} (⟨ E^{2} ⟩ - {⟨ E ⟩}^{2})

Can you see that what is inside the square brackets is exactly the variance of the energy defined in Section 3.1? It is ⟨E²⟩ weighted by the Boltzmann factor, minus ⟨E⟩².

3.4 Conclusion — the specific heat is the square of the uncertainty of the energy

We just substitute this into Eq. (4).

(5)

C_{V} = - T {(\frac{\partial^{2} F}{\partial T^{2}})}_{V} = \frac{1}{k_{B} T^{2}} (⟨ E^{2} ⟩ - {⟨ E ⟩}^{2}) = \frac{{(Δ E)}^{2}}{k_{B} T^{2}}

What cannot fluctuate cannot store heat Eq. (5) says C_V = 0 ⟺ ΔE = 0. A system whose energy has been fixed to a single value receives nothing even when the temperature is raised. Conversely, the wider the range of energies a system can take, the more heat it can store.
No quantum mechanics has entered here yet. Eq. (5) is a general law that holds for both classical and quantum systems. Quantum mechanics matters when the E_i become discrete as in Section 2, and then at low temperatures the system sticks to its lowest level, ΔE → 0, and the specific heat vanishes — from the next section on, we look at this concretely for a single mode.

Rereading the classical answer through Eq. (5) The classical theory of Section 1 gave C_V = 3Nk_B. Solving Eq. (5) in reverse gives ΔE = √(C_Vk_BT²) = √(3N) k_BT. In other words, the system is said to carry a fluctuation of k_BT per oscillator (mode), without fail, whatever the temperature. That this never disappears, however low the temperature, was the true content of the classical claim.

Fluctuations are invisible in macroscopic systems For 1 mol of a solid (C_V = 3R = 24.94 J/K) placed at 300 K, the energy fluctuation is ΔE = √(C_Vk_B)T = 5.57 nJ. Compared with the internal energy U ≈ 7.48 kJ, this is ΔE/U = 7.4×10⁻¹³, which is exactly 1/√(3N). In a macroscopic system the fluctuations become relatively negligible, so the energy ordinarily looks “definite.” And yet it can still be measured as a specific heat — we are looking at invisible fluctuations with a thermometer, and that is the fun part of Eq. (5).

4. Quantum Statistics of a Single Mode

4.1 The partition function

We build the Boltzmann distribution on the discrete levels of Eq. (2). Since the sum is a geometric series, it can be evaluated by elementary means.

Z = \sum_{n = 0}^{\infty} e^{- β (n + \frac{1}{2}) ħ ω} = \frac{e^{- β ħ ω / 2}}{1 - e^{- β ħ ω}}

Computing ⟨E⟩ = −∂(ln Z)/∂β gives

(6)

⟨ E ⟩ = \frac{1}{2} ħ ω + \frac{ħ ω}{e^{ħ ω / k_{B} T} - 1} = ħ ω (\frac{1}{2} + ⟨ n ⟩)

The second term, ⟨n⟩ = 1/(e^ħω/k_BT − 1), is the Bose–Einstein distribution itself, and gives the average number of phonons in that mode. Phonons are particles whose number is not conserved, so this is a Bose distribution with zero chemical potential.

The word “phonon” When the lattice vibrations of a crystal are decomposed into normal modes, each mode with a definite frequency ω becomes a harmonic oscillator independent of all the others (a single wave running through the whole crystal, labeled by its wavenumber q — counting them is the subject of Section 6). Its energy changes only in steps of ħω as in Eq. (2), so we regard one such quantum as a particle and call it a phonon, and rephrase “the mode is in the nth level” as “there are n phonons in that mode.” ⟨n⟩ is their average. A phonon is a name for counting the excitations of the lattice vibrations, not a real particle — when we said above that “the number is not conserved,” we meant that the level of the mode merely goes up and down, so phonons can be created and destroyed in any number.

4.2 The specific heat of one mode

Differentiate Eq. (6) with respect to T. The zero-point term contains no temperature, so it vanishes. With x ≡ ħω/k_BT, we obtain

(7)

C = k_{B} \frac{x^{2} e^{x}}{{(e^{x} - 1)}^{2}} (x = ħ ω / k_{B} T)

Let us look at the limits.

High temperature (x → 0): since e^x − 1 ≈ x, C → k_Bx²/x² = k_B. Equipartition is recovered.
Low temperature (x → ∞): C ≈ k_Bx²e^−x → 0. The exponential wins, so it drops off in no time.

4.3 The same answer also comes from the “fluctuation” of Section 3

Let us apply the conclusion of Section 3, Eq. (5), to this single mode. The energy is E_n = (n+1/2)ħω, so only n fluctuates.

{(Δ E)}^{2} = {(ħ ω)}^{2} (⟨ n^{2} ⟩ - {⟨ n ⟩}^{2}) = {(ħ ω)}^{2} ⟨ n ⟩ (⟨ n ⟩ + 1) = {(ħ ω)}^{2} \frac{e^{x}}{{(e^{x} - 1)}^{2}}

P_n ∝ e^−nx is a geometric distribution, so the variance comes out neatly as ⟨n⟩(⟨n⟩+1).

Substituting this into C_V = (ΔE)²/k_BT² of Eq. (5), and using ħω = xk_BT, we get

C = \frac{{(Δ E)}^{2}}{k_{B} T^{2}} = k_{B} \frac{x^{2} e^{x}}{{(e^{x} - 1)}^{2}} = Eq. (7) itself

We have arrived at the same equation by a completely different route. Eq. (7) is what you get by “differentiating Eq. (6) with respect to temperature,” whereas this one is what you get by “measuring the fluctuation of the energy” — that the two agree is the claim of Eq. (5). Incidentally, dividing Eq. (5) by k_B gives an extremely simple form.

\frac{C}{k_{B}} = {(\frac{Δ E}{k_{B} T})}^{2}

The specific heat of a single mode is the square of the energy fluctuation of that mode measured in units of k_BT.

x = ħω/k_BT	⟨n⟩	ΔE/k_BT	C/k_B	State
0.2	4.5167	0.9983	0.9967	Nearly classical: ΔE ≈ k_BT
0.5	1.5415	0.9897	0.9794	Still awake
1	0.5820	0.9595	0.9207	Starting to take effect
2	0.1565	0.8509	0.7241	Fluctuations are thinning out
4	0.0187	0.5514	0.3041	Dozing off
8	0.00034	0.1466	0.0215	Almost surely fixed at n = 0

“Falling asleep” meant “becoming unable to fluctuate” Look at the bottom row of the table. A mode with ħω ≫ k_BT has ⟨n⟩ ≈ 0, that is, it is almost certainly in n = 0. Being certain means that the energy does not fluctuate, so ΔE → 0. By Eq. (5) the specific heat is zero as well.
Conversely, for ħω ≪ k_BT the steps are so fine that the system can climb freely to any of them, and ΔE settles at exactly k_BT — just as classical theory says.
“The mode freezes” is just another way of saying “the mode can no longer fluctuate,” and this is how the levels in steps of ħω in Section 2 and the fluctuation formula in Section 3 are connected.

Zero-point energy does not affect the specific heat The ħω/2 in Eq. (6) contains no temperature, so it disappears when differentiated with respect to T. In other words, zero-point motion does not appear in the heat capacity. Even so, it shows up clearly in the temperature factor of X-ray diffraction and in the isotope effect (the lighter the atom, the larger the zero-point motion) — we will come back to this in Section 8.

5. The Einstein Model (1907) — First, Put In Only the “Discreteness”

Einstein considered that the N atoms in a crystal all vibrate at the same frequency ω_E, independently of one another. With three directions per atom there are 3N modes in all, each following Eq. (7). Writing θ_E ≡ ħω_E/k_B (the Einstein temperature), we have

(8)

C_{V} = 3 N k_{B} {(\frac{θ_{E}}{T})}^{2} \frac{e^{θ_{E} / T}}{{(e^{θ_{E} / T} - 1)}^{2}}

We have already seen the limits in Section 4.2.

T → ∞: C_V → 3Nk_B. Dulong–Petit is properly restored.
T → 0: C_V ≈ 3Nk_B(θ_E/T)²e^−θ_E/T → 0.

This was a major event, because it showed that the quantum hypothesis alone can explain the drop in the specific heat: it was the third example after blackbody radiation and the photoelectric effect, and the first to bring the quantum hypothesis into solids, matter itself. The small specific heat of diamond at room temperature is likewise explained by its large θ_E.

⚠ But it falls off too fast Experiment shows a power-law decrease, C_V ∝ T³, whereas Eq. (8) falls off as an exponential, e^−θ_E/T. The cause is clear: the assumption that “all the frequencies are the same” is too crude. A real crystal has sound waves, and by making the wavelength longer the frequency can be lowered without limit. Modes with ħω ≪ k_BT always exist, however low the temperature, so the specific heat cannot go to zero exponentially. The problem shifts to “how to count the modes.”

6. Counting the Modes — Periodic Boundary Conditions

6.1 The Born–von Kármán condition

A finite crystal has surfaces, but since the effect of the surfaces is small compared with the volume, it is usual to join the crystal up periodically instead. For the direction in which N₁ atoms line up along one edge, we require u(R + N₁a₁) = u(R). For a plane wave e^iq·R this becomes

e^{i N_{1} q \cdot a_{1}} = 1 ⟹ q \cdot a_{1} = \frac{2 π n_{1}}{N_{1}} (n_{1} \in ℤ)

The wavenumbers have become discrete. Taking all three directions together, the allowed q are lattice points evenly spaced in reciprocal space, and the volume occupied by one point is

Δ q^{3} = \frac{{(2 π)}^{3}}{V} (V = crystal volume)

The larger the crystal, the denser the points, but the number of points per unit volume does not change. This is the mechanism by which V later drops out neatly.

Fig. 1. The wavenumbers allowed by periodic boundary conditions, and the Debye cutoff. Drawn in two dimensions. Each lattice point is one wave (one mode), and the spacing is 2π/Na. The outer square is the first Brillouin zone. The Debye model replaces this zone with a sphere of the same volume (a circle in two dimensions) — its radius is q_D. Because the areas are equal, the total number of modes does not change, but it is not the same set: the corners of the zone (pale rings) are discarded, and the bulges outside the edges (orange rings) are picked up instead. That is the net content of the approximation.

Atoms per side N = 8

6.2 The density of states g(ω)

For long waves (q → 0) the lattice can be regarded as a continuum and becomes a sound wave. With sound velocity v, ω = vq. The number of lattice points inside a sphere of radius q is (V/(2π)³)(4π/3)q³, and we multiply this by the three polarizations (one longitudinal + two transverse). Differentiating with respect to ω gives the density of states.

(9)

g (ω) = \frac{3 V ω^{2}}{2 π^{2} v^{3}}

It increases in proportion to ω², and that is all there is to it — modes of low frequency are few, and the higher the frequency, the more modes there are. But the total number of modes is only 3N (three degrees of freedom per atom), so we have to cut it off somewhere. Debye chose the most straightforward way of cutting it off.

\int_{0}^{ω_{D}} g (ω) d ω = 3 N ⟹ ω_{D}^{3} = 6 π^{2} v^{3} \frac{N}{V}

V has dropped out. We define the Debye temperature as Θ_D ≡ ħω_D/k_B. The stiffer the crystal (the larger v) and the lighter the atoms, the higher Θ_D becomes — that is why diamond has 2230 K and lead 105 K.

7. The Debye Model (1912)

7.1 Internal energy and specific heat

All that remains is to weight Eq. (6) by g(ω) and add up.

E = \int_{0}^{ω_{D}} g (ω) ħ ω (\frac{1}{2} + \frac{1}{e^{ħ ω / k_{B} T} - 1}) d ω

Differentiating with respect to T and changing the variable to x = ħω/k_BT (ω² dω = (k_BT/ħ)³x²dx, with upper limit Θ_D/T), we obtain

(10)

C_{V} = 9 N k_{B} {(\frac{T}{Θ_{D}})}^{3} \int_{0}^{Θ_{D} / T} \frac{x^{4} e^{x}}{{(e^{x} - 1)}^{2}} d x

This integral cannot be written in elementary functions, but the limits at both ends can be obtained by hand. That is the most interesting part of Eq. (10).

7.2 High-temperature limit — Dulong–Petit returns

For T ≫ Θ_D the upper limit Θ_D/T is small, and x ≪ 1 throughout the interval of integration. The integrand becomes x⁴e^x/(e^x−1)² → x⁴/x² = x², so

\int_{0}^{Θ_{D} / T} x^{2} d x = \frac{1}{3} {(\frac{Θ_{D}}{T})}^{3} ⟹ C_{V} \to 3 N k_{B}

The Dulong–Petit law has been correctly reproduced as the high-temperature limit. What classical theory claimed to hold “always” actually came with a condition: “when k_BT is larger than the ħω_D of the stiffest mode by a sufficient margin.”

7.3 Low-temperature limit — the T³ law

For T ≪ Θ_D the upper limit is effectively ∞. This definite integral is a famous value that can be written with the Riemann zeta function. By integration by parts it reduces to ∫₀^∞x³/(e^x−1)dx = π⁴/15.

\int_{0}^{\infty} \frac{x^{4} e^{x}}{{(e^{x} - 1)}^{2}} d x = 4 \int_{0}^{\infty} \frac{x^{3}}{e^{x} - 1} d x = \frac{4 π^{4}}{15} \approx 25.976

(11)

C_{V} \approx \frac{12 π^{4}}{5} N k_{B} {(\frac{T}{Θ_{D}})}^{3} = 233.78 N k_{B} {(\frac{T}{Θ_{D}})}^{3}

This is the T³ law. Unlike Einstein’s exponential, it falls off as a power law because, thanks to g(ω) ∝ ω², “however low the temperature, some modes with ħω < k_BT always remain.” If you read this as saying that the number of modes that are awake is roughly proportional to (T/Θ_D)³, and each of them contributes k_B, the form of Eq. (11) becomes completely natural.

Consistency with the third law Using Eq. (11), S(T) = ∫₀^TC_VdT′/T′ = C_V(T)/3, which converges properly to 0 as T → 0. The worry of Section 1 is resolved. Note that the T³ law actually works only for T < Θ_D/50 or so, and at T = Θ_D/10 it is already off by about 3 %.

8. The Same Statistics Appears in the “Displacement” Too — ADP and Zero-Point Motion

Up to here we have been talking about energy. From the same Eq. (6) we can also obtain how much the atoms are shaking. Using exactly the same machinery as for the specific heat, we can derive the atomic displacement parameter (ADP) familiar from crystal structure analysis. Let us follow it step by step.

8.1 First, ⟨x²⟩ for a single harmonic oscillator

In a harmonic oscillator, the expectation values of the kinetic energy and the potential energy are exactly half each (the virial theorem, a general rule for V ∝ x²). For the nth level,

⟨ V ⟩ = \frac{1}{2} M ω^{2} ⟨ x^{2} ⟩ = \frac{E_{n}}{2} = \frac{1}{2} (n + \frac{1}{2}) ħ ω

All that remains is to solve this for ⟨x²⟩.

⟨ x^{2} ⟩_{n} = (n + \frac{1}{2}) \frac{ħ}{M ω}

In thermal equilibrium we replace n by the Bose-distribution average ⟨n⟩. Here, using the identity 2⟨n⟩ + 1 = coth(ħω/2k_BT) ((e^x+1)/(e^x−1) = coth(x/2)), everything comes together neatly.

⟨ u^{2} ⟩_{q ν} = \frac{ħ}{2 M ω_{q ν}} (2 ⟨ n ⟩ + 1) = \frac{ħ}{2 M ω_{q ν}} \coth \frac{ħ ω_{q ν}}{2 k_{B} T}

This coth will split into two later coth approaches 1 as x → ∞ and 2/x as x → 0. The former corresponds to low temperatures where only the zero point remains, the latter to high temperatures where we return to classical theory. If we let ħ → 0, the latter holds at every temperature, and the classical answer ⟨u²⟩ → k_BT/Mω² comes out — which is exactly the equipartition theorem of Section 1.

8.2 In a crystal, we sum over all phonon modes — the ADP tensor

In a crystal the frequency is not a single value but spreads out as the phonon dispersion ω_qν. The displacement parameter of atom κ becomes a second-rank tensor, and the expression actually evaluated in first-principles phonon calculations is the following.

U_{κ, α β} (T) = \frac{ħ}{2 N M_{κ}} \sum_{q ν} \frac{e_{κ α} e_{κ β}^{*}}{ω_{q ν}} \coth \frac{ħ ω_{q ν}}{2 k_{B} T}

e_κα(qν) is the “atom κ, direction α” component of the phonon eigenvector, and N is the number of q points taken in the Brillouin zone. The reason for the factor 1/N is physically clear: the absolute value of the plane wave e^iq·R is 1 in every unit cell, so the amplitude of one mode is shared among the N cells.

8.3 Taking the trace makes the eigenvectors disappear

The eigenvectors are normalized for each mode (Σ_α|e_κα|² = 1), so taking the trace of the tensor eliminates e.

⟨ u^{2} ⟩ = \sum_{α} U_{α α} = \frac{ħ}{2 N M} \sum_{q ν} \frac{1}{ω_{q ν}} \coth \frac{ħ ω_{q ν}}{2 k_{B} T}

No assumption of isotropy is needed here The trace is independent of the eigenvectors because it follows from the normalization |ê| = 1 alone. No matter how strong the anisotropy, this expression holds exactly. Only if we assume isotropy can we say U_xx = U_yy = U_zz, and write U_iso = ⟨u²⟩/3 (Section 8.6).
The sum runs over 3N modes (N q points × 3 branches), so 1/N × 3N = 3 — this 3 is the “three dimensions’ worth” that survives to the end.

8.4 Turning the sum into an integral — unlike the specific heat, a weight of 1/ω appears

Here we use the density of states built in Section 6, normalized per atom (∫g dω = 3, that is, three degrees of freedom per atom). For the Debye model,

g (ω) = \frac{9 ω^{2}}{ω_{D}^{3}} (0 \leq ω \leq ω_{D})

With this the sum becomes an integral.

⟨ u^{2} ⟩ = \frac{ħ}{2 M} \int_{0}^{ω_{D}} \frac{g (ω)}{ω} \coth \frac{ħ ω}{2 k_{B} T} d ω = \frac{9 ħ}{2 M ω_{D}^{3}} \int_{0}^{ω_{D}} ω \coth \frac{ħ ω}{2 k_{B} T} d ω

The ADP looks at the spectrum differently from the specific heat Since g(ω) ∝ ω² is multiplied by 1/ω, the integrand is the first power of ω. The specific heat, Eq. (10), has ω⁴e^…/(e^…−1)², where the high-frequency side is cut off exponentially. In other words, the ADP picks up soft (small-ω) modes strongly, while the specific heat discounts the contribution of stiff modes. Even for the same Debye model, the Θ_D fitted to the specific heat and the Θ_M fitted to diffraction differ by 10–20 %, because of this difference in weighting.

8.5 Splitting coth into “zero point” and “heat”

This is the crux, but all we do is rewrite one line.

\coth \frac{ħ ω}{2 k_{B} T} = 1 + \frac{2}{e^{ħ ω / k_{B} T} - 1}

The 1 on the left is the zero point (it contains no temperature), and the term on the right is the part excited by heat (= 2⟨n⟩).

The zero-point term contains no temperature, so it can be integrated as it is.

\frac{9 ħ}{2 M ω_{D}^{3}} \int_{0}^{ω_{D}} ω d ω = \frac{9 ħ}{2 M ω_{D}^{3}} \cdot \frac{ω_{D}^{2}}{2} = \frac{9 ħ}{4 M ω_{D}} = \frac{9 ħ^{2}}{4 M k_{B} Θ_{D}}

At the end we simply substituted ω_D = k_BΘ_D/ħ.

For the heat term, we substitute ξ = ħω/k_BT. Then dω = (k_BT/ħ)dξ, and the upper limit of the integral is ħω_D/k_BT = Θ_D/T.

\frac{9 ħ}{2 M ω_{D}^{3}} \cdot 2 {(\frac{k_{B} T}{ħ})}^{2} \int_{0}^{Θ_{D} / T} \frac{ξ}{e^{ξ} - 1} d ξ = \frac{9 ħ^{2}}{M k_{B} Θ_{D}} {(\frac{T}{Θ_{D}})}^{2} \int_{0}^{Θ_{D} / T} \frac{ξ}{e^{ξ} - 1} d ξ

Adding the two gives the three-dimensional mean-square displacement.

⟨ u^{2} ⟩ = \frac{9 ħ^{2}}{M k_{B} Θ_{D}} [\frac{1}{4} + {(\frac{T}{Θ_{D}})}^{2} \int_{0}^{Θ_{D} / T} \frac{ξ}{e^{ξ} - 1} d ξ]

The 1/4 in the square brackets is what remains after factoring the zero-point term out with the same prefactor as the heat term.

8.6 Per direction, it becomes the ADP

The U_iso of crystallography is 1/3 of the trace of the ADP tensor — the mean-square displacement per direction in the isotropic case.

U_{iso} = \frac{1}{3} (U_{x x} + U_{y y} + U_{z z}) = \frac{1}{3} ⟨ u^{2} ⟩

So we just divide by 3 — that is the only reason the 9 becomes a 3.

(12)

U_{iso} = \frac{1}{3} ⟨ u^{2} ⟩ = \frac{3 ħ^{2}}{M k_{B} Θ_{D}} [\frac{1}{4} + {(\frac{T}{Θ_{D}})}^{2} \int_{0}^{Θ_{D} / T} \frac{ξ}{e^{ξ} - 1} d ξ]

Correspondence with the way textbooks write it In crystallography books (Warren, Willis–Pryor), the Debye function Φ(x) = (1/x)∫₀^x ξ/(e^ξ−1)dξ is used, and the result is often written as B = (6h²/mk_BΘ)[Φ(x)/x + 1/4], x = Θ/T, but Φ(x)/x = (T/Θ)²∫₀^Θ/T…, and since h = 2πħ, 6h²/8π² = 3ħ² — exactly the same expression as Eq. (12).

8.7 The limits at both ends

	Value of the integral	U_iso	Meaning
T → 0	→ 0	3ħ²/(4Mk_BΘ_D)	No temperature dependence — zero-point motion
T → ∞	since ξ/(e^ξ−1) → 1, → Θ_D/T	3ħ²T/(Mk_BΘ_D²) = 3k_BT/(Mω_D²)	Proportional to T — classical equipartition

That U_iso ∝ T at high temperatures — that is, that the temperature factor B is proportional to the temperature — is well-known behavior in diffraction experiments. It comes straight out of Eq. (12). Conversely, at low temperatures B sticks to a constant value and gets no smaller — this is the reason the intensity of the diffraction lines does not fully recover however much you cool the sample.

Checking with numbers (Na) Putting in M = 22.990 and Θ_D = 158 K, the prefactor of Eq. (12) is 3ħ²/(Mk_BΘ_D) = 0.040064 Å². One quarter of that is the zero point, U_iso(0 K) = 0.010016 Å², and converted to the three-dimensional displacement, √(3U_iso) = 0.173 Å. The alternative route 3ħ/(4Mω_D) gives the same value. This is the number that tab ① of the simulator displays at 0 K.

8.8 Where does B = 8π²U_iso come from?

When there is thermal vibration, the scattering amplitude is multiplied by the Debye–Waller factor e^−W = exp[−½⟨(Q·u)²⟩] (e^−2W for the intensity). In the isotropic case ⟨(Q·u)²⟩ = Q²U_iso, so W = ½Q²U_iso. Inserting the magnitude of the scattering vector Q = 4π sinθ/λ, we obtain

W = 8 π^{2} U_{iso} {(\frac{\sin θ}{λ})}^{2} \equiv B {(\frac{\sin θ}{λ})}^{2} ⟹ B = 8 π^{2} U_{iso}

In other words, B is just U_iso dressed in the conventional units of diffraction; the content is the same. U_iso has the plain dimension of Å², so in recent structure analysis it has become usual to report U.

The 1/4 inside the square brackets is the zero-point motion It contains no temperature at all. It does not vanish even at T = 0, so the atoms keep shaking even at absolute zero. As we saw in Section 4.2, this term vanishes when differentiated with respect to T (it does not appear in Eq. (10) either), so it does not appear in the specific heat. Yet it remains in Eq. (12) — the same zero-point motion is invisible in the specific heat but visible in diffraction. It is what is drawn as the purple circle in tab ① of the simulator. In Na it amounts to √(3U_iso) = 0.173 Å even at 0 K, which is 36 % of the 0.480 Å at room temperature.

The Lindemann melting rule There is an empirical rule that a crystal melts when the amplitude reaches roughly 10–15 % of the nearest-neighbor distance. Computing the value at the melting point with Eq. (12), all five materials in the table below fall in the range 0.097–0.143 (roughly 10–14 %) — even though Θ_D ranges from 105 K to 2230 K, the atomic weight from 12 to 207, and the melting point from 371 K to 4100 K.

The ADP is very sensitive to “soft modes” At high temperatures coth(x) ≈ 1/x, so the sum in Section 8.3 takes the form U ∝ (k_BT/M) Σ 1/ω². With a weight of 1/ω², a single low-frequency mode is enough to make the ADP jump up. This is why the ADP comes out anomalously large in ferroelectrics with a soft mode, in filled skutterudites that exhibit rattling, and in negative thermal expansion materials.
However, a large ADP is not in itself evidence that the atom is responsible for the thermal expansion. Whether a mode depends on volume is determined by the Grüneisen parameter γ_qν = −∂lnω_qν/∂lnV, so a judgment can be made only after matching ω_qν, the eigenvector e_κα, and γ_qν mode by mode.

⚠ Cautions when dealing with real ADPs ① The isotropic approximation itself can be crude. The tensor U_αβ of Section 8.2 is the thermal vibration ellipsoid; for a hexagonal crystal, for example, symmetry gives U₁₁ = U₂₂, but U₃₃ can take a different value. Looking at U₃₃/U₁₁, the difference between the vibrations along the c axis and in the ab plane can be read off directly as an anisotropy. U_iso squashes that into a single number.
② Θ_D and Θ_M are different things. As seen in Section 8.4, the specific heat and the ADP weight the spectrum differently. When borrowing Θ_D from the literature, you need to check which of the two it was determined from.
③ The measured B includes more than vibration. Anharmonicity, static disorder (occupancy disorder, defects), residual errors of the absorption correction, and the like are added on top, so back-calculating Θ_D from B tends to give a low value.

9. Checking with Numbers

Results of numerically integrating Eqs. (10) and (12) for five materials. C_V is per atom.
Material	Θ_D / K	C_V/k_B at 300 K	% of 3k_B	Molar heat capacity J/(mol K)	B / Å² at 300 K	√(3U)/d at the melting point
Pb (lead)	105	2.9817	99.4	24.79	1.514	0.097
Na (sodium)	158	2.9588	98.6	24.60	6.052	0.143
Cu (copper)	343	2.8127	93.8	23.39	0.478	0.110
Al (aluminum)	428	2.7156	90.5	22.58	0.736	0.101
C (diamond)	2230	0.4971	16.6	4.13	0.120	0.113

Note that diamond alone is smaller by an order of magnitude. For diamond, room temperature is an extremely low temperature, T/Θ_D = 0.13, at which most of the stiff modes are asleep. For Pb, on the other hand, room temperature is a high temperature, T/Θ_D = 2.9, at which Dulong–Petit holds almost exactly. Whether room temperature counts as “high” is determined only by its ratio to Θ_D.

The measured molar heat capacities (C_P at 298 K) are about Pb 26.6, Na 28.2, Cu 24.4, Al 24.2, and diamond 6.1 J/(mol K). The orders of magnitude and the relative sizes (only diamond is smaller by an order of magnitude, Cu and Al are comparable, and so on) agree well. The measured values for metals are slightly larger than the calculated ones because C_P − C_V (the thermal expansion part), the contribution of the conduction electrons, and anharmonicity are added on top.

10. How Far Is It Correct?

⚠ What the Debye model leaves out ① The dispersion is approximated by sound waves. The true ω(q) levels off at the zone boundary, but Debye keeps ω = vq and chops it off abruptly at q_D. As a result it is off by a few % at intermediate temperatures (around T ≈ Θ_D/5), and Θ_D is not constant when it is measured as a function of temperature.
② There are no optical modes. If the unit cell contains two or more atoms, optical branches appear. Their frequencies are almost constant, so they are rather well represented by the Einstein model (in practice one treats the acoustic branches with Debye and the optical branches with Einstein, and adds them up).
③ The electronic specific heat is not included. In metals C = γT + AT³, and below a few K the electronic term γT wins. Plotting C/T against T² gives a straight line, from whose intercept you read γ and from whose slope you read Θ_D — a standard experimental method.
④ Anharmonicity is ignored. Since terms of third order and higher in the Taylor expansion are dropped, neither thermal expansion nor the finiteness of the thermal conductivity comes out. This is why the measured values go beyond 3R at high temperatures.

11. Summary

In classical theory the specific heat stays put at 3k_B. Equipartition merely hands (1/2)k_BT to each quadratic-form degree of freedom, and neither ħ, nor the material, nor the temperature remains in the formula (Dulong–Petit).
The specific heat is the fluctuation of the energy itself. Just by matching F = U − TS against F = −k_BTlnZ, we obtain C_V = −T(∂²F/∂T²)_V = (ΔE)²/k_BT². What cannot fluctuate cannot store heat. No quantum mechanics has entered here yet.
Quantum mechanics makes the energy take discrete values in steps of ħω. ΔxΔp ≥ ħ/2 demands the zero-point energy ħω/2, and solving the Schrödinger equation gives the levels E_n = (n+1/2)ħω.
Modes that cannot afford a quantum fall asleep. A mode with ħω ≫ k_BT stays at n = 0 and takes up no energy even when the temperature is raised — it does not contribute to the specific heat.
Periodic boundary conditions make the modes countable. q is discretized as 2πn/Na, g(ω) ∝ ω², and cutting off at a total of 3N fixes Θ_D.
Both limits are contained in Debye’s single formula. 3k_B for T → ∞ (the revival of Dulong–Petit), (12π⁴/5)k_B(T/Θ_D)³ for T → 0 (the T³ law).
The same statistics appears in the displacement too. The 1/4 in the ADP, Eq. (12), is the zero-point motion, which does not show up in the specific heat but does show up in diffraction. A crystal melts when the amplitude reaches about a tenth of the nearest-neighbor distance (Lindemann).

Corresponding simulator: Phonons and Specific Heat.
Section 3 follows the argument of “A detour: specific heat arises from the uncertainty of energy” in Lecture 3 of the lecture materials “Computational Materials Science.”
All equations are typeset in the browser’s standard MathML (no external libraries or network communication are used). Θ_D, atomic weights, lattice constants, metallic radii, and melting points are standard literature values. The integrals in Eqs. (10) and (12) are evaluated with 48-point Gauss–Legendre quadrature, and the agreement with the T³ law is to 8 digits at T = 0.01Θ_D.