Bose gases: photons, phonons and Bose–Einstein condensation

Bosons crowd into the same states, and that tendency produces some of the most important phenomena in physics: the spectrum of thermal radiation, the heat capacity of solids at low temperature, superfluidity, and the condensation of a gas into a single quantum state. This lesson treats three Bose gases. Photons and phonons are bosons whose number is not conserved, which makes their chemical potential zero and their statistics especially simple. Massive bosons such as helium-4 or rubidium-87 atoms have a conserved number, and below a critical temperature a macroscopic fraction of them falls into the lowest orbital, a Bose–Einstein condensate, predicted by Einstein in 1925 and made in the laboratory seventy years later.

Photons: the Planck distribution

Electromagnetic radiation in a cavity at temperature \(T\) is a gas of photons. Photons are created and absorbed freely by the walls, so their number is not fixed; the equilibrium number is whatever minimizes the free energy, which requires \(\partial F/\partial N = \mu = 0\). With \(\mu = 0\), the Bose–Einstein occupancy of a mode of frequency \(\nu\) is the Planck distribution:

The same result follows from treating each mode of the field as a quantum harmonic oscillator, as Planck effectively did in 1900, since the mean number of quanta in an oscillator is \(1/(e^{h\nu/kT} - 1)\). Counting modes in a box, with two polarizations and \(\varepsilon = h\nu = cp\), gives a density of states proportional to \(\nu^2\), and multiplying by \(h\nu\,\bar n\) gives the energy density per unit frequency:

Integrating over frequency with \(x = h\nu/kT\) gives \(u = (8\pi k^4T^4/h^3c^3)\int_0^\infty x^3/(e^x - 1)\,dx\). The integral equals \(\pi^4/15\), so the total energy density is \(aT^4\), and the power radiated per unit area of a black surface, \(cu/4\), is \(\sigma T^4\) with

The Stefan–Boltzmann law, found experimentally by Josef Stefan in 1879 and derived thermodynamically by Boltzmann in 1884, now has its constant computed from \(h\), \(k\) and \(c\). The classical theory, which gives each mode \(kT\) by equipartition, assigns infinite energy to the infinitely many high-frequency modes, the ultraviolet catastrophe; the Bose–Einstein factor cuts the high modes off exponentially. The program evaluates both Planck integrals numerically and recovers \(\sigma\). It also counts the photons of the cosmic microwave background at 2.7255 K: about 411 per cubic centimeter, filling all of space, a relic of the hot early universe whose blackbody spectrum was measured with extraordinary precision by the COBE satellite around 1990. PHYS-261 plotted the Planck spectrum and used it for the Earth's energy balance; here the emphasis is the statistics behind it.

The thermodynamics of a photon gas

Because \(\mu = 0\), the grand potential of the photon gas equals its Helmholtz free energy, and the rest of its thermodynamics follows quickly. Summing \(kT\ln(1 - e^{-h\nu/kT})\) over modes gives \(F = -\frac13 aT^4V\), where \(u = aT^4\) is the energy density. Then the entropy is \(S = -\partial F/\partial T = \frac43 aT^3V\) and the pressure is \(P = -\partial F/\partial V = \frac13 u\). Radiation pressure is one third of the energy density, a result also obtained from the momentum of light, and it supports the most massive stars against gravity. If radiation expands without exchanging heat, its entropy stays constant, so \(T^3V\) is constant and \(T\) falls in proportion to the size of the container. That is exactly what has happened to the cosmic background radiation: released at about 3,000 K when the universe became transparent, it has cooled to 2.7 K as the universe expanded by a factor of about 1,100, keeping a perfect blackbody spectrum all the while.

Phonons: the Debye model

The vibrations of a crystal lattice are quantized into phonons, and like photons they are bosons with \(\mu = 0\). Einstein's 1907 model gave every vibration the same frequency, which correctly made the heat capacity vanish at low temperature but predicted an exponential fall-off, much too fast. Real crystals have vibrational modes of all frequencies, from long-wavelength sound waves at nearly zero frequency upward. In 1912 Peter Debye treated the lattice as an elastic solid: three sound polarizations with a density of states proportional to \(\nu^2\), as for photons, cut off at a maximum frequency \(\nu_D\) chosen so that the total number of modes is \(3N\), the number of degrees of freedom of \(N\) atoms. With the Debye temperature \(\theta_D = h\nu_D/k\), the heat capacity is

At high temperature the integral approaches \((\theta_D/T)^3/3\) and \(C_V \to 3Nk\), the Dulong–Petit value. At low temperature the upper limit can be taken to infinity, the integral becomes \(4\pi^4/15\), and

the Debye \(T^3\) law. The \(T^3\) arises for the same reason as the \(T^4\) of radiation: at low temperature only low-frequency phonons, which behave like massless particles with a \(\nu^2\) density of states, are excited. Debye temperatures range from about 105 K for soft lead to 343 K for copper and about 2,200 K for diamond, which is why diamond's heat capacity is far below Dulong–Petit at room temperature.

import math
k = 1.380649e-23; h = 6.62607015e-34; hbar = h / (2 * math.pi); c = 2.99792458e8; u = 1.66053907e-27

def simpson(f, a, b, n=20000):
    hh = (b - a) / n
    s = f(a) + f(b) + sum((4 if i % 2 else 2) * f(a + i * hh) for i in range(1, n))
    return s * hh / 3

# 1. Photon gas: the Planck integrals and the Stefan-Boltzmann constant
I3 = simpson(lambda x: x ** 3 / math.expm1(x) if x > 0 else 0.0, 0, 60)
I2 = simpson(lambda x: x ** 2 / math.expm1(x) if x > 0 else 0.0, 0, 60)
print(f"integral x^3/(e^x - 1) = {I3:.6f}, pi^4/15 = {math.pi ** 4 / 15:.6f}")
print(f"integral x^2/(e^x - 1) = {I2:.6f}, 2 zeta(3) = {2 * 1.2020569:.6f}")
sigma = 2 * math.pi * k ** 4 * I3 / (h ** 3 * c ** 2)
print(f"Stefan-Boltzmann constant from the integral: {sigma:.4e} W/(m^2 K^4)")
T = 2.7255
n_ph = 8 * math.pi * (k * T / (h * c)) ** 3 * I2
print(f"cosmic microwave background at {T} K: {n_ph / 1e6:.1f} photons per cm^3, "
      f"energy density {4 * sigma / c * T ** 4 * 1e15:.3f} x 10^-15 J/m^3")

# 2. Phonons: the Debye heat capacity C/(3Nk) = 3 (T/theta)^3 int_0^{theta/T} x^4 e^x/(e^x-1)^2 dx
def debye(t):                                  # t = T/theta
    if t == 0: return 0.0
    xm = 1 / t
    g = lambda x: x ** 4 * math.exp(x) / math.expm1(x) ** 2 if x > 0 else 0.0
    return 3 * t ** 3 * simpson(g, 0, xm, 4000)
theta = 343.0                                  # copper
print("\nDebye model for copper (theta_D = 343 K): C as a fraction of 3Nk")
for T in (5, 10, 20, 50, 100, 200, 300):
    lowT = 12 * math.pi ** 4 / 5 * (T / theta) ** 3 / 3
    extra = f"   T^3 law gives {lowT:.5f}" if T <= 50 else ""
    print(f"  T = {T:3} K: C/3Nk = {debye(T / theta):.5f}{extra}")

# 3. Bose-Einstein condensation of an ideal gas
zeta32 = 2.612375
def Tc(m, n): return 2 * math.pi * hbar ** 2 / (m * k) * (n / zeta32) ** (2 / 3)
m_He = 4.0026 * u
print(f"\nideal-gas Tc for liquid helium density 145 kg/m^3: {Tc(m_He, 145 / m_He):.2f} K (lambda point 2.17 K)")
m_Rb = 86.909 * u
for n in (1e18, 1e19, 1e20):
    print(f"  Rb-87 at n = {n:.0e} per m^3: Tc = {Tc(m_Rb, n) * 1e9:6.0f} nK")
print("condensate fraction 1 - (T/Tc)^1.5 (box) and 1 - (T/Tc)^3 (harmonic trap):")
for t in (0.25, 0.5, 0.75, 0.9):
    print(f"  T/Tc = {t}: box {1 - t ** 1.5:.3f}, trap {1 - t ** 3:.3f}")

integral x^3/(e^x - 1) = 6.493939, pi^4/15 = 6.493939
integral x^2/(e^x - 1) = 2.404114, 2 zeta(3) = 2.404114
Stefan-Boltzmann constant from the integral: 5.6704e-08 W/(m^2 K^4)
cosmic microwave background at 2.7255 K: 410.7 photons per cm^3, energy density 41.748 x 10^-15 J/m^3

Debye model for copper (theta_D = 343 K): C as a fraction of 3Nk
  T =   5 K: C/3Nk = 0.00024   T^3 law gives 0.00024
  T =  10 K: C/3Nk = 0.00193   T^3 law gives 0.00193
  T =  20 K: C/3Nk = 0.01545   T^3 law gives 0.01545
  T =  50 K: C/3Nk = 0.19984   T^3 law gives 0.24139
  T = 100 K: C/3Nk = 0.59199
  T = 200 K: C/3Nk = 0.86709
  T = 300 K: C/3Nk = 0.93757

ideal-gas Tc for liquid helium density 145 kg/m^3: 3.13 K (lambda point 2.17 K)
  Rb-87 at n = 1e+18 per m^3: Tc =     18 nK
  Rb-87 at n = 1e+19 per m^3: Tc =     86 nK
  Rb-87 at n = 1e+20 per m^3: Tc =    398 nK
condensate fraction 1 - (T/Tc)^1.5 (box) and 1 - (T/Tc)^3 (harmonic trap):
  T/Tc = 0.25: box 0.875, trap 0.984
  T/Tc = 0.5: box 0.646, trap 0.875
  T/Tc = 0.75: box 0.350, trap 0.578
  T/Tc = 0.9: box 0.146, trap 0.271

The program evaluates the Debye function for copper. The \(T^3\) law is accurate to all printed digits up to about 20 K, which is \(\theta_D/17\), and fails badly by 50 K. At room temperature copper's lattice heat capacity is 94 parts in 100 of the Dulong–Petit value. Plotted against \(T/\theta_D\), the heat capacities of many elements fall close to a single Debye curve, an early example of a universal law in statistical physics.

Bose–Einstein condensation

For massive bosons the number is conserved, so \(\mu\) is set by the density. Lesson 8 gave \(n\lambda^3 = f_{3/2}(z)\) with fugacity \(z = e^{\mu/kT}\) (for spinless bosons). As the gas is cooled at fixed density, \(\lambda\) grows, so \(f_{3/2}(z)\) must grow, and \(z\) rises toward 1, which is \(\mu \to 0\), the energy of the lowest orbital. But \(f_{3/2}(z)\) has a maximum value at \(z = 1\), namely \(\zeta(3/2) = 2.612\). When \(n\lambda^3\) reaches 2.612 the excited orbitals can hold no more particles. This happens at the critical temperature

Below \(T_c\), the excited orbitals hold only \(N(T/T_c)^{3/2}\) particles; the rest, a macroscopic number, crowd into the single lowest orbital. The condensate fraction is

What went wrong with the integral? Replacing the sum over orbitals by an integral with density of states \(g(\varepsilon) \propto \sqrt\varepsilon\) gives the ground state zero weight, which is harmless in every other case but not when one orbital holds a finite fraction of all the particles. Einstein's 1925 paper treated the ground state separately and found this transition, which has no counterpart in classical physics: it occurs in an ideal gas, with no forces at all, purely from the statistics. In a harmonic trap, as used in experiments, the density of states grows as \(\varepsilon^2\) instead, and the condensate fraction is \(1 - (T/T_c)^3\).

Real condensates

For decades the only candidate was liquid helium-4. Below 2.17 K it becomes a superfluid that flows without friction, discovered by Pyotr Kapitsa in Moscow and, independently, by John Allen and Don Misener, two Canadian physicists working at Cambridge, in 1937 and 1938. Fritz London proposed in 1938 that this was Bose–Einstein condensation in a liquid. The program supports him: an ideal gas of helium atoms at the density of the liquid would condense at 3.1 K, close to the observed 2.17 K. The difference is due to the strong interactions in a liquid, which also reduce the condensate fraction to about a tenth even at the lowest temperatures, although the whole liquid becomes superfluid.

A nearly ideal condensate required a dilute atomic gas cooled to well under a microkelvin, using laser cooling followed by evaporative cooling in a magnetic trap. In June 1995 Eric Cornell and Carl Wieman at JILA in Boulder produced a condensate of about 2,000 rubidium-87 atoms at about 170 nK, and a few months later Wolfgang Ketterle's group at MIT made a much larger one with sodium atoms. The three shared the 2001 Nobel Prize in Physics. The program shows why such low temperatures are needed: at a typical trap density of \(10^{19}\) atoms per cubic meter, \(T_c\) for rubidium is under 100 nK. Images of the expanding cloud showed a sharp spike of atoms with nearly zero velocity rising out of a broad thermal distribution as the temperature crossed \(T_c\), a direct picture of the lowest orbital filling up. Condensates are now routine tools for atom interferometry, precision measurement and the simulation of many-body quantum systems.