# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0
import numpy as np
import matplotlib.pyplot as plt
k_B = 1.380649e-23 # J/K
def fermi_dirac(epsilon, mu, T, k_B):
"""Fermi-Dirac distribution"""
beta = 1 / (k_B * T)
x = beta * (epsilon - mu)
# Avoid overflow
if x > 100:
return 0
elif x < -100:
return 1
return 1 / (np.exp(x) + 1)
def bose_einstein(epsilon, mu, T, k_B):
"""Bose-Einstein distribution"""
beta = 1 / (k_B * T)
x = beta * (epsilon - mu)
if x <= 0:
return np.inf # μ < ε required
if x < 0.01:
return 1 / x # Approximation
return 1 / (np.exp(x) - 1)
def maxwell_boltzmann(epsilon, mu, T, k_B):
"""Maxwell-Boltzmann distribution"""
beta = 1 / (k_B * T)
x = beta * (epsilon - mu)
if x > 100:
return 0
return np.exp(-x)
# Energy range (chemical potential set to 0)
mu = 0
epsilon_range = np.linspace(-5, 5, 500) # In units of k_B T
# Different temperatures (normalized in k_B T units)
k_B_T = 1 # Normalized unit
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Distribution at T = 1 (in k_B T units)
ax1 = axes[0, 0]
n_FD = [fermi_dirac(eps * k_B_T, mu, 1, k_B) for eps in epsilon_range]
n_BE = [bose_einstein(eps * k_B_T, mu, 1, k_B) if eps > 0 else 0 for eps in epsilon_range]
n_MB = [maxwell_boltzmann(eps * k_B_T, mu, 1, k_B) for eps in epsilon_range]
ax1.plot(epsilon_range, n_FD, 'b-', linewidth=2, label='Fermi-Dirac')
ax1.plot(epsilon_range, n_BE, 'r-', linewidth=2, label='Bose-Einstein')
ax1.plot(epsilon_range, n_MB, 'g--', linewidth=2, label='Maxwell-Boltzmann')
ax1.axvline(0, color='k', linestyle=':', linewidth=1, label='μ = 0')
ax1.set_xlabel('(ε - μ) / (k_B T)')
ax1.set_ylabel('Occupation number ⟨n⟩')
ax1.set_title('Quantum Statistical Distributions (T = 1)')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_ylim([0, 3])
# Fermi-Dirac at low temperature
ax2 = axes[0, 1]
temperatures_FD = [0.1, 0.5, 1.0, 2.0]
colors_FD = ['blue', 'green', 'orange', 'red']
for T_val, color in zip(temperatures_FD, colors_FD):
n_FD_T = [fermi_dirac(eps * k_B_T, mu, T_val, k_B) for eps in epsilon_range]
ax2.plot(epsilon_range, n_FD_T, color=color, linewidth=2,
label=f'k_B T = {T_val}')
ax2.axvline(0, color='k', linestyle=':', linewidth=1)
ax2.set_xlabel('(ε - μ) / (E_F)')
ax2.set_ylabel('⟨n⟩')
ax2.set_title('Temperature Dependence of Fermi-Dirac Distribution')
ax2.legend()
ax2.grid(True, alpha=0.3)
# Bose-Einstein distribution and Bose condensation
ax3 = axes[1, 0]
epsilon_BE = np.linspace(0.01, 5, 100)
temperatures_BE = [0.5, 1.0, 2.0, 5.0]
for T_val, color in zip(temperatures_BE, colors_FD):
n_BE_T = [bose_einstein(eps * k_B_T, mu * 0.9, T_val, k_B) for eps in epsilon_BE]
ax3.semilogy(epsilon_BE, n_BE_T, color=color, linewidth=2,
label=f'k_B T = {T_val}')
ax3.set_xlabel('(ε - μ) / (k_B T)')
ax3.set_ylabel('⟨n⟩ (log scale)')
ax3.set_title('Bose-Einstein Distribution')
ax3.legend()
ax3.grid(True, alpha=0.3, which='both')
# Classical limit verification
ax4 = axes[1, 1]
epsilon_classical = np.linspace(0, 10, 100)
T_high = 5 # High temperature
n_FD_high = [fermi_dirac(eps * k_B_T, mu, T_high, k_B) for eps in epsilon_classical]
n_BE_high = [bose_einstein(eps * k_B_T, mu * 0.5, T_high, k_B) for eps in epsilon_classical]
n_MB_high = [maxwell_boltzmann(eps * k_B_T, mu * 0.5, T_high, k_B) for eps in epsilon_classical]
ax4.semilogy(epsilon_classical, n_FD_high, 'b-', linewidth=2, label='Fermi-Dirac')
ax4.semilogy(epsilon_classical, n_BE_high, 'r-', linewidth=2, label='Bose-Einstein')
ax4.semilogy(epsilon_classical, n_MB_high, 'g--', linewidth=2, label='Maxwell-Boltzmann')
ax4.set_xlabel('(ε - μ) / (k_B T)')
ax4.set_ylabel('⟨n⟩ (log scale)')
ax4.set_title(f'Classical Limit (k_B T = {T_high})')
ax4.legend()
ax4.grid(True, alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('stat_mech_quantum_statistics.png', dpi=300, bbox_inches='tight')
plt.show()
print("=== Quantum Statistics ===\n")
print("Fermi-Dirac distribution:")
print(" - Fermions (electrons, protons, neutrons, etc.)")
print(" - Pauli exclusion principle: ⟨n⟩ ≤ 1")
print(" - Step function at T = 0 (Fermi surface)\n")
print("Bose-Einstein distribution:")
print(" - Bosons (photons, phonons, He-4, etc.)")
print(" - Multiple occupation possible: ⟨n⟩ ≥ 0")
print(" - Bose condensation at low temperature (μ → ε₀)\n")
print("Classical limit (high temperature):")
print(" - When e^{β(ε-μ)} >> 1")
print(" - FD ≈ BE ≈ MB")
print(" - Quantum statistics → classical statistics\n")
# Numerical example
epsilon_test = 2 * k_B_T
T_test = 300 # K
print(f"Numerical example (ε-μ = 2k_BT, T = {T_test} K):")
n_FD_test = fermi_dirac(epsilon_test, mu, T_test, k_B)
n_BE_test = bose_einstein(epsilon_test, mu * 0.5, T_test, k_B)
n_MB_test = maxwell_boltzmann(epsilon_test, mu * 0.5, T_test, k_B)
print(f" Fermi-Dirac: ⟨n⟩ = {n_FD_test:.6f}")
print(f" Bose-Einstein: ⟨n⟩ = {n_BE_test:.6f}")
print(f" Maxwell-Boltzmann: ⟨n⟩ = {n_MB_test:.6f}")