Quantum Statistics in a Nutshell

From Classical to Quantum: When Does It Matter?¶

In the canonical ensemble we derived the partition function for $N$ indistinguishable non-interacting particles:

Z = \frac{Z_1^N}{N!}, \qquad Z_1 = \frac{V}{\lambda^3}

(1)

where $\lambda = h/\sqrt{2\pi m k_B T}$ is the thermal de Broglie wavelength.

This classical (Maxwell-Boltzmann) approximation is valid when the mean occupation of every single-particle state is much less than one:

\bar{n}_i \ll 1 \quad \Longleftrightarrow \quad n\lambda^3 \ll 1

(2)

where $n = N/V$ is the number density. The dimensionless quantity $n\lambda^3$ compares the thermal wavelength to the mean interparticle spacing:

Regime	Condition	Statistics	Examples
Classical	$n\lambda^3 \ll 1$	Maxwell-Boltzmann	Gases at room $T$ , solutions, biomolecules
Quantum	$n\lambda^3 \gtrsim 1$	Bose-Einstein or Fermi-Dirac	Electrons in metals, $^4$ He at 2K, neutron stars

Quantum statistics matters when particles are light (small $m \Rightarrow$ large $\lambda$ ), cold (small $T \Rightarrow$ large $\lambda$ ), or dense (large $n$ ).

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import h, k, N_A, m_e, m_p

# Thermal de Broglie wavelength
def lambda_dB(m, T):
    return h / np.sqrt(2 * np.pi * m * k * T)

T = np.linspace(1, 1000, 500)

systems = {
    'Electron':       (m_e, 8.5e28, '#d62728'),     # free electron in copper
    r'$^4$He atom':   (4*m_p, 2.2e28, '#1f77b4'),   # liquid helium density
    r'$H_2$ gas (1 atm)': (2*m_p, 2.7e25, '#2ca02c'),  # H2 at STP
    r'$N_2$ gas (1 atm)': (28*m_p, 2.7e25, '#ff7f0e'),  # N2 at STP
}

fig, ax = plt.subplots(figsize=(9, 5))
for label, (m, n, color) in systems.items():
    lam = lambda_dB(m, T)
    ax.semilogy(T, n * lam**3, lw=2.5, color=color, label=label)

ax.axhline(1, color='gray', ls='--', lw=1.5)
ax.fill_between(T, 1, 1e6, alpha=0.08, color='blue')
ax.fill_between(T, 1e-30, 1, alpha=0.08, color='green')
ax.text(50, 30, 'Quantum regime', fontsize=13, color='blue', fontweight='bold')
ax.text(500, 1e-8, 'Classical regime', fontsize=13, color='green', fontweight='bold')
ax.set_xlabel('Temperature (K)', fontsize=12)
ax.set_ylabel(r'$n\lambda^3$', fontsize=14)
ax.set_title('Classical vs Quantum Regime', fontsize=14, fontweight='bold')
ax.set_ylim(1e-15, 1e5)
ax.set_xlim(1, 1000)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3, which='both')
plt.tight_layout()
plt.show()

Occupation Number Representation¶

For non-interacting particles, a microstate is specified by the occupation numbers $\{n_1, n_2, \ldots\}$ of single-particle energy levels $\{\epsilon_1, \epsilon_2, \ldots\}$ :

E = \sum_i n_i \epsilon_i, \qquad N = \sum_i n_i

(3)

The key question is: what values can $n_i$ take?

Particle type	Spin	Allowed $n_i$	Example
Bosons	Integer (0, 1, 2, ...)	$0, 1, 2, 3, \ldots$	Photons, phonons, $^4$ He
Fermions	Half-integer (1/2, 3/2, ...)	0 or 1 only	Electrons, protons, $^3$ He

This single constraint — whether particles can share quantum states — leads to profoundly different statistical behavior.

The Three Distribution Laws¶

Working in the grand canonical ensemble (where the $N = \sum_i n_i$ constraint is relaxed), the partition function factorizes over single-particle levels:

\Xi = \prod_i \Xi_i, \qquad \Xi_i = \sum_{n=0}^{n_{\max}} e^{-\beta(\epsilon_i - \mu)n}

(4)

Why the grand canonical ensemble?

In the canonical ensemble, the constraint $N = \sum_i n_i$ couples all occupation numbers and prevents the partition function from factorizing. In the grand canonical ensemble, $N$ is allowed to fluctuate, and the sum over all $\{n_i\}$ becomes unrestricted:

\Xi = \sum_{n_1, n_2, \ldots} e^{-\beta \sum_i (\epsilon_i - \mu) n_i} = \prod_i \left[\sum_{n=0}^{n_{\max}} e^{-\beta(\epsilon_i - \mu)n}\right] = \prod_i \Xi_i

(5)

Each factor $\Xi_i$ is a geometric series:

Bosons ( $n_{\max} = \infty$ ): $\Xi_i = (1 - e^{-\beta(\epsilon_i - \mu)})^{-1}$
Fermions ( $n_{\max} = 1$ ): $\Xi_i = 1 + e^{-\beta(\epsilon_i - \mu)}$

The mean occupation follows from $\langle n_i \rangle = \frac{1}{\beta}\frac{\partial \ln \Xi_i}{\partial \mu}$ , and fluctuations from $\sigma^2_{n_i} = \frac{1}{\beta^2}\frac{\partial^2 \ln \Xi_i}{\partial \mu^2}$ .

The resulting mean occupation number takes a universal form:

Quantum distribution functions

\boxed{\langle n_i \rangle = \frac{1}{e^{\beta(\epsilon_i - \mu)} + c}}

(6)

$c$	Distribution	Partition function per level
-1	Bose-Einstein	$\Xi_i = (1 - e^{-\beta(\epsilon_i - \mu)})^{-1}$
+1	Fermi-Dirac	$\Xi_i = 1 + e^{-\beta(\epsilon_i - \mu)}$
0	Maxwell-Boltzmann	$\Xi_i = e^{e^{-\beta(\epsilon_i - \mu)}}$

Bose-Einstein: Occupation is unbounded. Bosons tend to bunch into the same state. At $T \to 0$ , macroscopic occupation of the ground state leads to Bose-Einstein condensation. Note: for the ground state ( $\epsilon_0 = 0$ ), $\langle n_0 \rangle = (e^{-\beta\mu} - 1)^{-1}$ requires $\mu \leq 0$ to keep occupation non-negative. The chemical potential of a boson gas is always less than or equal to the ground state energy.
Fermi-Dirac: Occupation is 0 or 1 (Pauli exclusion). Fermions avoid each other. At $T = 0$ , all states below the Fermi energy $\epsilon_F = \mu(T=0)$ are filled, forming a sharp step function.
Maxwell-Boltzmann: The classical limit recovered when $e^{\beta(\epsilon_i - \mu)} \gg 1$ (i.e. $\mu \to -\infty$ , low density), where both quantum distributions converge to $\langle n_i \rangle \approx e^{-\beta(\epsilon_i - \mu)}$ .

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-2, 8, 500)  # x = beta*(epsilon - mu)

# Avoid division by zero
be = 1.0 / (np.exp(x) - 1)  # Bose-Einstein
fd = 1.0 / (np.exp(x) + 1)  # Fermi-Dirac
mb = np.exp(-x)              # Maxwell-Boltzmann

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Left panel: all three distributions
ax = axes[0]
ax.plot(x, be, lw=2.5, color='#1f77b4', label='Bose-Einstein ($c=-1$)')
ax.plot(x, fd, lw=2.5, color='#d62728', label='Fermi-Dirac ($c=+1$)')
ax.plot(x, mb, lw=2.5, color='#2ca02c', ls='--', label='Maxwell-Boltzmann ($c=0$)')
ax.axvline(0, color='gray', ls=':', lw=1)
ax.set_xlabel(r'$\beta(\epsilon_i - \mu)$', fontsize=13)
ax.set_ylabel(r'$\langle n_i \rangle$', fontsize=13)
ax.set_title('Quantum vs Classical Distributions', fontsize=13, fontweight='bold')
ax.set_xlim(-2, 6)
ax.set_ylim(0, 4)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Right panel: Fermi-Dirac at different temperatures
ax = axes[1]
eps = np.linspace(0, 15, 500)  # energy in units of some scale
mu_val = 7.0  # Fermi energy

for kT, color in zip([0.2, 0.5, 1.0, 2.0, 5.0], 
                      ['#1a1a2e', '#16213e', '#0f3460', '#e94560', '#ff9a00']):
    n_fd = 1.0 / (np.exp((eps - mu_val) / kT) + 1)
    ax.plot(eps, n_fd, lw=2, color=color, label=rf'$k_BT = {kT}$')

ax.axvline(mu_val, color='gray', ls='--', lw=1.5)
ax.text(mu_val + 0.2, 0.85, r'$\epsilon_F = \mu$', fontsize=11, color='gray')
ax.set_xlabel(r'$\epsilon$', fontsize=13)
ax.set_ylabel(r'$\langle n \rangle$', fontsize=13)
ax.set_title('Fermi-Dirac: Sharpness at Low $T$', fontsize=13, fontweight='bold')
ax.legend(fontsize=10)
ax.set_ylim(-0.05, 1.1)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Density of States¶

To compute thermodynamic quantities, we convert sums over discrete levels to integrals using the density of states $g(\epsilon)$ :

\sum_i (\cdots) \rightarrow \int_0^\infty g(\epsilon)\, (\cdots)\, d\epsilon

(7)

For free particles in 3D with $\epsilon = \hbar^2 k^2/2m$ :

g(\epsilon) = \frac{V}{4\pi^2} \left(\frac{2m}{\hbar^2}\right)^{3/2} \epsilon^{1/2} \propto \sqrt{\epsilon}

(8)

Derivation of $g(\epsilon)$ from k-space counting

For a particle in a box with periodic boundary conditions, the allowed wavevectors are $k_\alpha = 2\pi n_\alpha / L_\alpha$ , so each state occupies a volume $(2\pi)^3/V$ in k-space. The number of states with $|k| < k$ is:

N(k) = \frac{4\pi k^3/3}{(2\pi)^3/V} = \frac{Vk^3}{6\pi^2}

(9)

Differentiating and converting to energy via $\epsilon = \hbar^2 k^2/2m$ :

g(k) = \frac{dN}{dk} = \frac{Vk^2}{2\pi^2} \qquad \Rightarrow \qquad g(\epsilon) = g(k)\frac{dk}{d\epsilon} = \frac{V}{4\pi^2}\left(\frac{2m}{\hbar^2}\right)^{3/2} \sqrt{\epsilon}

(10)

For photons ( $\epsilon = \hbar\omega = \hbar ck$ ) with 2 polarizations: $g(\omega) = V\omega^2/(\pi^2 c^3)$ .

The total energy and particle number become:

\langle E \rangle = \int_0^\infty \epsilon\, g(\epsilon)\, \langle n(\epsilon) \rangle\, d\epsilon, \qquad \langle N \rangle = \int_0^\infty g(\epsilon)\, \langle n(\epsilon) \rangle\, d\epsilon

(11)

Application 1: Blackbody Radiation (Photon Gas)¶

Photons are massless bosons with $\mu = 0$ (photon number is not conserved). Applying Bose-Einstein statistics with the photon density of states $g(\omega) = V\omega^2/(\pi^2 c^3)$ gives:

Key results:

Wien’s displacement law: Peak frequency shifts $\propto T$
Stefan-Boltzmann law: Total energy density $u \propto T^4$
Classical limit ( $\hbar \to 0$ ): Rayleigh-Jeans formula $u \propto \omega^2 T$ , which diverges at high $\omega$ (ultraviolet catastrophe)

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import c, h, k

def planck(wav, T):
    """Planck spectral radiance as a function of wavelength."""
    return 2.0 * h * c**2 / (wav**5 * (np.exp(h * c / (wav * k * T)) - 1.0))

wavs = np.arange(100e-9, 3000e-9, 5e-9)

fig, ax = plt.subplots(figsize=(9, 5))
colors = ['#1a1a2e', '#0f3460', '#e94560', '#ff9a00', '#ffcc00']
for T_val, color in zip([3000, 4000, 5000, 6000, 7000], colors):
    ax.plot(wavs * 1e9, planck(wavs, T_val) * 1e-12, lw=2.5, color=color,
            label=f'{T_val} K')

ax.set_xlabel(r'Wavelength (nm)', fontsize=12)
ax.set_ylabel(r'Spectral radiance (TW/m$^3$/sr)', fontsize=12)
ax.set_title('Planck Blackbody Spectrum', fontsize=14, fontweight='bold')
ax.set_xlim(0, 2500)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

# Shade visible range
ax.axvspan(380, 750, alpha=0.1, color='orange', label='Visible')
ax.text(550, ax.get_ylim()[1]*0.9, 'Visible', ha='center', fontsize=10, 
        color='orange', fontweight='bold')

plt.tight_layout()
plt.show()

Application 2: Heat Capacity of Solids (Phonon Gas)¶

Atoms in a crystal vibrate as coupled harmonic oscillators. In normal-mode coordinates, the solid becomes a gas of phonons — quantized lattice vibrations that obey Bose-Einstein statistics with $\mu = 0$ .

Why $\mu = 0$ for phonons (and photons)?

The partition function for a quantum harmonic oscillator with frequency $\omega_i$ is:

Z_i = \sum_{n=0}^{\infty} e^{-\beta \hbar \omega_i (n + 1/2)} = e^{-\beta \hbar \omega_i / 2} \cdot \frac{1}{1 - e^{-\beta \hbar \omega_i}}

(13)

The second factor is exactly the bosonic grand canonical partition function $\Xi_i$ with $\mu = 0$ . This is because the number of phonons (or photons) is not conserved — they are created and destroyed freely as the crystal vibrates. Since there is no constraint on total phonon number, no chemical potential is needed to enforce it.

The mean occupation of mode $i$ is:

\langle n_i \rangle = \frac{1}{2} + \frac{1}{e^{\beta \hbar \omega_i} - 1}

(14)

where the $\frac{1}{2}$ is the quantum zero-point contribution.

The Debye model approximates the phonon density of states as $g(\omega) = 9N\omega^2/\omega_D^3$ up to a cutoff $\omega_D$ (set by requiring $\int_0^{\omega_D} g(\omega)\,d\omega = 3N$ total modes). The resulting heat capacity is:

C_V(T) = 9Nk_B \left(\frac{T}{\Theta_D}\right)^3 \int_0^{\Theta_D/T} \frac{x^4 e^x}{(e^x - 1)^2}\, dx

(15)

where $\Theta_D = \hbar \omega_D / k_B$ is the Debye temperature.

Limit	$C_V$	Physics
$T \gg \Theta_D$	$3Nk_B$ (Dulong-Petit)	All modes classically excited
$T \ll \Theta_D$	$\propto T^3$	Only low-frequency acoustic modes are excited

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.constants import k, N_A

def debye_cv(T_arr, T_D):
    """Debye heat capacity per mole."""
    def integrand(x):
        return x**4 * np.exp(x) / (np.exp(x) - 1)**2
    
    cv = np.zeros_like(T_arr)
    for i, T in enumerate(T_arr):
        if T < 0.01:
            cv[i] = 0
        else:
            cv[i] = 9 * N_A * k * (T / T_D)**3 * quad(integrand, 0, T_D / T)[0]
    return cv

T = np.linspace(1, 600, 500)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

# Left: C_v vs T for different Debye temperatures
materials = {'Diamond (2230 K)': 2230, 'Silicon (645 K)': 645, 
             'Copper (343 K)': 343, 'Lead (105 K)': 105}
colors = ['#1f77b4', '#2ca02c', '#ff7f0e', '#d62728']
for (name, T_D), color in zip(materials.items(), colors):
    cv = debye_cv(T, T_D)
    ax1.plot(T, cv, lw=2.5, color=color, label=name)

ax1.axhline(3 * N_A * k, color='gray', ls='--', lw=1.5)
ax1.text(450, 3 * N_A * k + 0.5, 'Dulong-Petit: $3Nk_B$', fontsize=10, color='gray')
ax1.set_xlabel('Temperature (K)', fontsize=12)
ax1.set_ylabel(r'$C_V$ (J/mol/K)', fontsize=12)
ax1.set_title('Debye Heat Capacity', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# Right: Universal curve C_v vs T/T_D
T_reduced = np.linspace(0.01, 3, 500)
cv_universal = debye_cv(T_reduced * 343, 343) / (3 * N_A * k)  # normalized
ax2.plot(T_reduced, cv_universal, lw=3, color='steelblue')
ax2.set_xlabel(r'$T / \Theta_D$', fontsize=13)
ax2.set_ylabel(r'$C_V / 3Nk_B$', fontsize=13)
ax2.set_title('Universal Debye Curve', fontsize=13, fontweight='bold')
ax2.axhline(1, color='gray', ls='--', lw=1.5)

# Show T^3 behavior
T_low = np.linspace(0.01, 0.3, 100)
cv_low = debye_cv(T_low * 343, 343) / (3 * N_A * k)
ax2.plot(T_low, cv_low, lw=3, color='steelblue')
ax2.annotate(r'$\propto T^3$', xy=(0.15, 0.02), fontsize=14,
             xytext=(0.5, 0.15), color='red', fontweight='bold',
             arrowprops=dict(arrowstyle='->', color='red', lw=1.5))
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

/tmp/ipykernel_3048/2195249607.py:9: RuntimeWarning: overflow encountered in exp
  return x**4 * np.exp(x) / (np.exp(x) - 1)**2
/tmp/ipykernel_3048/2195249607.py:9: RuntimeWarning: invalid value encountered in scalar divide
  return x**4 * np.exp(x) / (np.exp(x) - 1)**2
/tmp/ipykernel_3048/2195249607.py:9: RuntimeWarning: overflow encountered in scalar power
  return x**4 * np.exp(x) / (np.exp(x) - 1)**2
/tmp/ipykernel_3048/2195249607.py:9: RuntimeWarning: overflow encountered in scalar multiply
  return x**4 * np.exp(x) / (np.exp(x) - 1)**2
/tmp/ipykernel_3048/2195249607.py:16: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  cv[i] = 9 * N_A * k * (T / T_D)**3 * quad(integrand, 0, T_D / T)[0]

Application 3: Free Electron Gas (Fermi Gas)¶

Electrons in a metal are fermions. At $T = 0$ , they fill all available states up to the Fermi energy:

\epsilon_F = \frac{\hbar^2}{2m_e} \left(\frac{3\pi^2 N}{V}\right)^{2/3}

(16)

For typical metals, $\epsilon_F \sim 5{-}10$ eV, corresponding to a Fermi temperature $T_F = \epsilon_F / k_B \sim 50{,}000{-}100{,}000$ K. Since room temperature $T \ll T_F$ , the electron gas is deeply quantum even at 300 K.

Key consequences:

Only electrons within $\sim k_BT$ of $\epsilon_F$ contribute to heat capacity: $C_V^{\text{el}} \propto T$ (linear, not the classical $\frac{3}{2}Nk_B$ )
Even at $T = 0$ , electrons exert degeneracy pressure (Pauli exclusion prevents collapse)
This explains why metals are good conductors but have surprisingly small electronic heat capacity

import numpy as np
import matplotlib.pyplot as plt

# Fermi gas: occupation * density of states at different T
eps = np.linspace(0, 12, 500)
mu = 7.0  # eV, typical Fermi energy

fig, ax = plt.subplots(figsize=(9, 5))

# Density of states g(e) ~ sqrt(e)
dos = np.sqrt(eps)

for kT, alpha, color in [(0.05, 0.9, '#1a1a2e'), (0.2, 0.7, '#0f3460'), 
                          (0.5, 0.5, '#e94560'), (1.5, 0.4, '#ff9a00')]:
    fd = 1.0 / (np.exp((eps - mu) / kT) + 1)
    occupied_dos = dos * fd
    ax.fill_between(eps, 0, occupied_dos, alpha=alpha * 0.4, color=color)
    ax.plot(eps, occupied_dos, lw=2, color=color, label=rf'$k_BT = {kT}$ eV')

# Show full DOS for reference
ax.plot(eps, dos, lw=1.5, color='gray', ls='--', label=r'$g(\epsilon) \propto \sqrt{\epsilon}$')
ax.axvline(mu, color='black', ls=':', lw=1.5)
ax.text(mu + 0.15, ax.get_ylim()[1]*0.55, r'$\epsilon_F$', fontsize=14, fontweight='bold')

ax.set_xlabel(r'Energy $\epsilon$ (eV)', fontsize=12)
ax.set_ylabel(r'$g(\epsilon) \langle n(\epsilon) \rangle$', fontsize=12)
ax.set_title('Occupied States in a Free Electron Gas', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.set_xlim(0, 12)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Summary: When to Use What¶

System	Statistics	Why
Gas molecules at room $T$	Maxwell-Boltzmann	$n\lambda^3 \sim 10^{-6}$
Molecules in solution	Maxwell-Boltzmann	Dilute, heavy, warm
Photons (radiation)	Bose-Einstein ( $\mu = 0$ )	Massless bosons
Phonons (lattice vibrations)	Bose-Einstein ( $\mu = 0$ )	Quantized collective modes
Electrons in metals	Fermi-Dirac	$T_F \sim 10^5$ K $\gg T_{\text{room}}$
Liquid $^4$ He below 2.17 K	Bose-Einstein	Superfluid transition
Liquid $^3$ He	Fermi-Dirac	Fermionic helium isotope

Problems¶

Show that both Bose-Einstein and Fermi-Dirac distributions reduce to the Maxwell-Boltzmann distribution in the limit $e^{\beta(\epsilon_i - \mu)} \gg 1$ . What physical condition does this correspond to?
For a 2D free electron gas (relevant to graphene), the density of states is constant: $g(\epsilon) = Am/(\pi\hbar^2)$ where $A$ is the area. Derive the Fermi energy and show that $C_V \propto T$ at low temperatures.
Use the Planck distribution to show that the total energy density of a photon gas scales as $u \propto T^4$ (Stefan-Boltzmann law). Hint: substitute $x = \beta\hbar\omega$ and use $\int_0^\infty x^3/(e^x - 1)\, dx = \pi^4/15$ .
Compute the occupation number fluctuations $\sigma^2_{n_i}$ for bosons and fermions. Show that:
- Bosons: $\sigma^2_{n_i} = \langle n_i \rangle(1 + \langle n_i \rangle)$ — super-Poissonian (bunching)
- Fermions: $\sigma^2_{n_i} = \langle n_i \rangle(1 - \langle n_i \rangle)$ — sub-Poissonian (antibunching)
- Classical: $\sigma^2_{n_i} = \langle n_i \rangle$ — Poissonian
Estimate $n\lambda^3$ for: (a) nitrogen gas at STP, (b) electrons in copper ( $n_e \approx 8.5 \times 10^{28}$ m $^{-3}$ ), (c) water molecules in liquid water. Which systems require quantum statistics?