Part I: Ensembles | Chapter 4

Quantum Statistics

Identical particles, the spin-statistics theorem, Bose-Einstein and Fermi-Dirac distributions, density of states, and occupation number formalism

Historical Context

In 1924, Satyendra Nath Bose derived Planck's radiation law by treating photons as indistinguishable particles, sending his paper to Einstein who extended the method to massive particles. In 1926, Enrico Fermi and Paul Dirac independently developed the statistics for particles obeying the Pauli exclusion principle. The spin-statistics theorem, proved by Pauli in 1940, established that integer-spin particles (bosons) obey Bose-Einstein statistics while half-integer-spin particles (fermions) obey Fermi-Dirac statistics.

These quantum distributions reduce to Maxwell-Boltzmann statistics in the dilute (\(n\lambda_{dB}^3 \ll 1\)) limit but produce dramatically different behavior at low temperatures, giving rise to phenomena like superconductivity, superfluidity, and white dwarf stability.

1. Identical Particles and Symmetry

In quantum mechanics, identical particles are truly indistinguishable. The exchange operator \(\hat{P}_{12}\) satisfies \(\hat{P}_{12}^2 = 1\), so its eigenvalues are \(\pm 1\):

\[\hat{P}_{12}|\psi\rangle = \pm|\psi\rangle\]

Bosons (+1)

Integer spin (0, 1, 2, ...). Wave function symmetric under exchange. Any number can occupy the same state.

Examples: photons, phonons, \(^4\text{He}\), gluons, W/Z bosons, Higgs

Fermions (-1)

Half-integer spin (1/2, 3/2, ...). Wave function antisymmetric under exchange. At most one per quantum state (Pauli exclusion).

Examples: electrons, protons, neutrons, \(^3\text{He}\), quarks

Derivation 1: Counting States for Identical Particles

For \(N\) distinguishable particles distributed among \(G\) states, the number of arrangements is \(G^N\). For identical particles:

\[\text{Bosons: } \binom{G + N - 1}{N} = \frac{(G+N-1)!}{N!(G-1)!}\]

\[\text{Fermions: } \binom{G}{N} = \frac{G!}{N!(G-N)!}\]

In the dilute limit \(N \ll G\), both reduce to \(G^N/N!\), which is the Maxwell-Boltzmann result with the Gibbs correction.

2. Quantum Distribution Functions

Derivation 2: Bose-Einstein Distribution

For bosons in a single-particle state \(k\) with energy \(\epsilon_k\), the occupation number \(n_k\) can be 0, 1, 2, .... The grand partition function for this mode is:

\[\xi_k = \sum_{n_k=0}^{\infty} e^{-\beta(\epsilon_k - \mu)n_k} = \frac{1}{1 - e^{-\beta(\epsilon_k - \mu)}}\]

The mean occupation number:

\[\boxed{\langle n_k \rangle_{\text{BE}} = \frac{1}{e^{\beta(\epsilon_k - \mu)} - 1}}\]

For bosons, \(\mu \leq \epsilon_{\text{min}}\) (usually \(\mu \leq 0\)) to ensure non-negative occupation numbers.

Derivation 3: Fermi-Dirac Distribution

For fermions, \(n_k = 0\) or \(1\) only:

\[\xi_k = 1 + e^{-\beta(\epsilon_k - \mu)}\]

\[\boxed{\langle n_k \rangle_{\text{FD}} = \frac{1}{e^{\beta(\epsilon_k - \mu)} + 1}}\]

At \(T = 0\), this becomes a step function:\(\langle n_k \rangle = \theta(\mu - \epsilon_k)\). The chemical potential at\(T = 0\) defines the Fermi energy\(\epsilon_F = \mu(T=0)\).

3. Density of States

Derivation 4: Free Particle Density of States

For free particles in a box of volume \(V = L^3\) with periodic boundary conditions, the allowed momenta are \(\mathbf{p} = \hbar\mathbf{k}\) with\(\mathbf{k} = (2\pi/L)(n_x, n_y, n_z)\). The number of states with energy less than \(\epsilon\):

\[\Phi(\epsilon) = g_s \cdot \frac{V}{(2\pi)^3} \cdot \frac{4\pi}{3}k_{\max}^3 = g_s \frac{V}{6\pi^2}\left(\frac{2m\epsilon}{\hbar^2}\right)^{3/2}\]

where \(g_s\) is the spin degeneracy. The density of states is:

\[\boxed{g(\epsilon) = \frac{d\Phi}{d\epsilon} = g_s \frac{V}{4\pi^2}\left(\frac{2m}{\hbar^2}\right)^{3/2} \epsilon^{1/2}}\]

In \(d\) dimensions, \(g(\epsilon) \propto \epsilon^{d/2 - 1}\). For 2D: \(g(\epsilon) = \text{const}\); for 1D: \(g(\epsilon) \propto \epsilon^{-1/2}\).

4. Thermodynamic Quantities

All thermodynamic quantities for quantum gases can be expressed as integrals over the density of states weighted by the distribution function:

\[\langle N \rangle = \int_0^{\infty} g(\epsilon) f(\epsilon)\,d\epsilon\]

\[\langle E \rangle = \int_0^{\infty} \epsilon\, g(\epsilon) f(\epsilon)\,d\epsilon\]

\[PV = \frac{2}{3}\langle E \rangle \quad \text{(non-relativistic, 3D)}\]

Derivation 5: Equation of State \(PV = \frac{2}{3}E\)

Starting from the grand potential for a single mode:

\[PV = k_BT \ln\Xi = \mp k_BT \sum_k \ln(1 \mp ze^{-\beta\epsilon_k})\]

Converting to an integral and integrating by parts:

\[PV = \mp k_BT \int_0^{\infty} g(\epsilon)\ln(1 \mp ze^{-\beta\epsilon})\,d\epsilon\]

Using \(g(\epsilon) \propto \epsilon^{1/2}\) in 3D and integrating by parts with \(\Phi(\epsilon) \propto \epsilon^{3/2}\):

\[PV = \frac{2}{3}\int_0^{\infty}\epsilon\, g(\epsilon) f(\epsilon)\,d\epsilon = \frac{2}{3}\langle E \rangle\]

This is a purely kinematic result, valid for both bosons and fermions, depending only on the non-relativistic dispersion \(\epsilon = p^2/(2m)\).

5. Comparing the Three Statistics

Unified Formula

\[\langle n \rangle = \frac{1}{e^{\beta(\epsilon - \mu)} + a}, \qquad a = \begin{cases} +1 & \text{Fermi-Dirac} \\ 0 & \text{Maxwell-Boltzmann} \\ -1 & \text{Bose-Einstein}\end{cases}\]

In the limit \(\beta(\epsilon - \mu) \gg 1\) (dilute regime), all three converge to \(\langle n \rangle \approx e^{-\beta(\epsilon - \mu)}\). Quantum effects become important when \(n\lambda_{dB}^3 \gtrsim 1\), i.e., when the interparticle spacing becomes comparable to the de Broglie wavelength.

Degeneracy Parameter

The dimensionless quantity \(n\lambda_{dB}^3\) (where \(n = N/V\)) determines when quantum effects become important:

\[n\lambda_{dB}^3 = n\left(\frac{2\pi\hbar^2}{mk_BT}\right)^{3/2}\]

For air at room temperature: \(n\lambda_{dB}^3 \sim 10^{-6}\) (classical). For electrons in a metal: \(n\lambda_{dB}^3 \sim 10^3\) (deeply quantum). For liquid \(^4\text{He}\) at 2K: \(n\lambda_{dB}^3 \sim 7\) (strongly quantum).

6. Application: Photon Gas and Planck Distribution

Photons are massless bosons with \(\mu = 0\) (photon number is not conserved). The mean number of photons per mode with frequency \(\omega\) is:

\[\langle n(\omega) \rangle = \frac{1}{e^{\beta\hbar\omega} - 1}\]

The energy density per unit frequency is the Planck distribution:

\[u(\omega) = \frac{\hbar\omega^3}{\pi^2 c^3}\frac{1}{e^{\beta\hbar\omega} - 1}\]

Integrating gives the Stefan-Boltzmann law: \(U/V = (\pi^2 k_B^4 / 15\hbar^3 c^3)T^4\).

7. Computational Exploration

This simulation compares quantum distribution functions, demonstrates the density of states, shows Planck's radiation law, and visualizes the crossover from quantum to classical regimes.

Quantum Statistics: Distribution Functions, Planck Spectrum, and Density of States

Python

script.py146 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# Quantum Statistics: Distribution Functions and Applications
# ============================================================

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle('Quantum Statistics', fontsize=16, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#1a0a0a')
    ax.tick_params(colors='white')
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('#fca5a5')
    for spine in ax.spines.values():
        spine.set_color('#333')

# --- Panel 1: Comparison of distribution functions ---
ax1 = axes[0, 0]
eps_range = np.linspace(0, 5, 500)
mu = 1.0
kT = 0.5
beta = 1.0 / kT
x = beta * (eps_range - mu)
# Fermi-Dirac
n_FD = 1.0 / (np.exp(x) + 1)
# Bose-Einstein (only for eps > mu)
mask_BE = eps_range > mu + 0.01
n_BE = np.zeros_like(eps_range)
n_BE[mask_BE] = 1.0 / (np.exp(x[mask_BE]) - 1)
# Maxwell-Boltzmann
n_MB = np.exp(-x)

ax1.plot(eps_range, n_FD, color='#60a5fa', linewidth=2, label='Fermi-Dirac')
ax1.plot(eps_range[mask_BE], n_BE[mask_BE], color='#f87171', linewidth=2, label='Bose-Einstein')
ax1.plot(eps_range, n_MB, color='#34d399', linewidth=2, linestyle='--', label='Maxwell-Boltzmann')
ax1.axvline(x=mu, color='#fbbf24', linestyle=':', alpha=0.5, label='mu')
ax1.set_xlabel('Energy (epsilon)')
ax1.set_ylabel('<n(epsilon)>')
ax1.set_title('Distribution Functions (kT = 0.5, mu = 1)')
ax1.set_ylim(0, 4)
ax1.legend(fontsize=8, facecolor='#1a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: Planck distribution ---
ax2 = axes[0, 1]
x_planck = np.linspace(0.01, 15, 500)
for T_scale in [1.0, 1.5, 2.0, 3.0]:
    # u(x) ~ x^3 / (exp(x/T) - 1), x = hbar*omega/(k_B * T_ref)
    u = x_planck**3 / (np.exp(x_planck / T_scale) - 1)
    u /= np.max(u)  # normalize for display
    ax2.plot(x_planck, u, linewidth=1.5, label='T = ' + str(T_scale) + ' T_ref')

# Wien displacement
for T_scale in [1.0, 1.5, 2.0, 3.0]:
    x_peak = 2.82 * T_scale
    ax2.axvline(x=x_peak, color='white', linestyle=':', alpha=0.15)

ax2.set_xlabel('hbar * omega / (k_B * T_ref)')
ax2.set_ylabel('u(omega) [normalized]')
ax2.set_title('Planck Radiation Spectrum')
ax2.legend(fontsize=8, facecolor='#1a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Fermi function at various temperatures ---
ax3 = axes[1, 0]
eps_F = 5.0  # Fermi energy
for kT_val in [0.01, 0.2, 0.5, 1.0, 2.0]:
    eps = np.linspace(0, 10, 500)
    f = 1.0 / (np.exp((eps - eps_F) / kT_val) + 1)
    ax3.plot(eps, f, linewidth=1.5, label='kT = ' + str(kT_val))

ax3.axvline(x=eps_F, color='#fbbf24', linestyle='--', alpha=0.5, label='eps_F')
ax3.set_xlabel('Energy')
ax3.set_ylabel('f(epsilon)')
ax3.set_title('Fermi Function at Different Temperatures')
ax3.legend(fontsize=7, facecolor='#1a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: Density of states in different dimensions ---
ax4 = axes[1, 1]
eps = np.linspace(0.01, 5, 500)
# 3D: g ~ eps^{1/2}
g_3d = np.sqrt(eps)
g_3d /= np.max(g_3d)
# 2D: g ~ const
g_2d = np.ones_like(eps) * 0.5
# 1D: g ~ eps^{-1/2}
g_1d = 1.0 / np.sqrt(eps)
g_1d = np.clip(g_1d / 3.0, 0, 1)

ax4.plot(eps, g_3d, color='#f87171', linewidth=2, label='3D: g ~ eps^(1/2)')
ax4.plot(eps, g_2d, color='#60a5fa', linewidth=2, label='2D: g = const')
ax4.plot(eps, g_1d, color='#34d399', linewidth=2, label='1D: g ~ eps^(-1/2)')
ax4.set_xlabel('Energy (epsilon)')
ax4.set_ylabel('g(epsilon) [normalized]')
ax4.set_title('Density of States by Dimension')
ax4.legend(fontsize=8, facecolor='#1a0a0a', edgecolor='#333', labelcolor='white')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Plot saved to output.png")

# Numerical results
print("\n" + "=" * 60)
print("QUANTUM STATISTICS: NUMERICAL RESULTS")
print("=" * 60)

print("\n--- Occupation numbers at beta*(eps-mu) = 1 ---")
x_val = 1.0
n_be = 1.0 / (np.exp(x_val) - 1)
n_fd = 1.0 / (np.exp(x_val) + 1)
n_mb = np.exp(-x_val)
print("  Bose-Einstein: " + str(round(n_be, 6)))
print("  Fermi-Dirac:   " + str(round(n_fd, 6)))
print("  Maxwell-Boltz: " + str(round(n_mb, 6)))

print("\n--- Degeneracy parameter n * lambda_dB^3 ---")
# For electrons in copper
hbar = 1.055e-34
me = 9.109e-31
kB = 1.381e-23
n_cu = 8.5e28  # electrons/m^3
for T in [300, 1000, 10000]:
    lam = np.sqrt(2 * np.pi * hbar**2 / (me * kB * T))
    deg = n_cu * lam**3
    print("  Electrons in Cu, T=" + str(T) + "K: n*lam^3 = " + str(round(deg, 2)))

print("\n--- Planck distribution peak (Wien law) ---")
print("  x_peak = hbar*omega_peak/(kT) = 2.822")
for T in [3000, 5778, 10000]:
    omega_peak = 2.822 * kB * T / hbar
    lambda_peak = 2 * np.pi * 3e8 / omega_peak
    print("  T = " + str(T) + " K: lambda_peak = " + str(round(lambda_peak * 1e9, 1)) + " nm")

print("\n--- Stefan-Boltzmann law ---")
sigma = 2 * np.pi**5 * kB**4 / (15 * (2*np.pi*hbar)**3 * (3e8)**2)
print("  sigma = " + str(round(sigma, 4e-8)) + " W/(m^2 K^4)")
for T in [300, 5778, 40000]:
    power = sigma * T**4
    print("  T = " + str(T) + " K: radiated power = " + str(round(power, 1)) + " W/m^2")
print("\n" + "=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

8. Summary and Key Results

Core Formulas

BE: \(\langle n \rangle = 1/(e^{\beta(\epsilon-\mu)}-1)\)
FD: \(\langle n \rangle = 1/(e^{\beta(\epsilon-\mu)}+1)\)
DOS (3D): \(g(\epsilon) \propto \epsilon^{1/2}\)
Kinematic: \(PV = \frac{2}{3}E\)

Physical Insights

Spin-statistics theorem: bosons symmetric, fermions antisymmetric
Quantum effects when \(n\lambda_{dB}^3 \gtrsim 1\)
Photons: \(\mu = 0\), Planck spectrum
Classical limit recovered for dilute systems

← Grand Canonical Ensemble Ideal Quantum Gases →