statistical-mechanics

# Statistical Mechanics ## What I do I provide comprehensive expertise in statistical mechanics, the bridge between microscopic physics and thermodynamics. I enable you to apply ensemble theory, calculate partition functions, analyze phase transitions, implement Monte Carlo simulations, and model non-equilibrium dynamics. My knowledge spans from foundational statistical foundations to modern computational methods essential for condensed matter physics, soft matter, biological physics, and materials science. ## When to use me Use statistical mechanics when you need to: derive thermodynamic properties from microscopic models, analyze critical phenomena and phase transitions, simulate many-body systems using Monte Carlo, calculate equation of state for gases and materials, understand transport phenomena and fluctuations, model magnetic systems and spin models, or connect molecular dynamics to macroscopic observables. ## Core Concepts - **Microstates and Macrostates**: Microscopic configurations and their macroscopic aggregations with entropy S = k_B ln Ω. - **Phase Space**: Complete specification of particle positions and momenta for classical systems. - **Ensembles**: Collections of microstates with different constraints (microcanonical, canonical, grand canonical). - **Partition Functions**: Z = Σ exp(-βE) encoding all thermodynamic information. - **Thermodynamic Relations**: Connecting partition functions to free energy, entropy, and equations of state. - **Phase Transitions**: Qualitative changes in system behavior with singularities in derivatives of free energy. - **Mean-Field Theory**: Approximate treatment replacing interactions with effective fields. - **Critical Phenomena**: Universal behavior near phase transitions characterized by critical exponents. - **Fluctuation-Dissipation**: Connecting response functions to equilibrium fluctuations. - **Non-Equilibrium Dynamics**: Time evolution toward equilibrium and current flows. ## Code Examples ### Ensembles and Partition Functions ```python import numpy as np def boltzmann_factor(E, T): """exp(-E/k_B T)""" kB = 1.38e-23 return np.exp(-E / (kB * T)) def canonical_partition_function(energies, T): """Z = Σ exp(-βE_i)""" beta = 1 / (1.38e-23 * T) return np.sum(np.exp(-beta * np.array(energies))) def canonical_expectation(E, Z, T): """⟨E⟩ = (1/Z) Σ E_i exp(-βE_i)""" beta = 1 / (1.38e-23 * T) weights = np.exp(-beta * np.array(E)) return np.sum(E * weights) / Z def heat_capacity(energies, T, Cv): """Cv = (⟨E²⟩ - ⟨E⟩²) / (k_B T²)""" beta = 1 / (1.38e-23 * T) Z = canonical_partition_function(energies, T) E_avg = canonical_expectation(energies, Z, T) E2_avg = canonical_expectation([E**2 for E in energies], Z, T) return (E2_avg - E_avg**2) / ((1.38e-23) * T**2) # Example: Two-level system E0, E1 = 0, 0.001 # Energy levels in J (corresponds to ~72 K) energies = [E0, E1] print("Two-level system:") for T in [10, 50, 100, 300, 1000]: Z = canonical_partition_function(energies, T) E_avg = canonical_expectation(energies, Z, T) print(f" T = {T:4d} K: Z = {Z:.4f}, ⟨E⟩ = {E_avg*1e21:.2f} × 10⁻²¹ J") # Partition function for harmonic oscillator def harmonic_oscillator_z(T, omega): """Z = exp(-βℏω/2) / (1 - exp(-βℏω))""" hbar = 1.055e-34 beta = 1 / (1.38e-23 * T) if omega * hbar * beta > 100: return np.exp(-omega * hbar * beta / 2) return np.exp(-omega * hbar * beta / 2) / (1 - np.exp(-omega * hbar * beta)) # Free particle partition function (classical limit) def free_partition_function_1D(V, T, m=1.67e-27): """Z = V / λ³ where λ = h/√(2πmkT)""" h = 6.626e-34 kB = 1.38e-23 lam = h / np.sqrt(2 * np.pi * m * kB * T) return V / lam**3 # 3D ideal gas def ideal_gas_pressure(N, T, V, m=1.67e-27): """P = NkT/V (ideal gas law)""" kB = 1.38e-23 return N * kB * T / V N_A = 6.022e23 V = 0.0224 # 1 mole at STP T = 273 P = ideal_gas_pressure(N_A, T, V) print(f"\nIdeal gas at STP: P = {P/1e5:.2f} bar") # Grand canonical ensemble def grand_canonical_partition_function(grand_mu, T, single_particle_states): """Ξ = Π_k (1 + exp(-β(ε_k - μ))) for fermions""" beta = 1 / (1.38e-23 * T) Xi = 1 for eps in single_particle_states: Xi *= (1 + np.exp(-beta * (eps - grand_mu))) return Xi # Fermi-Dirac and Bose-Einstein occupation def fermi_dirac_occupation(eps, mu, T): """f(ε) = 1/(exp((ε-μ)/kT) + 1)""" kB = 1.38e-23 return 1 / (np.exp((eps - mu) / (kB * T)) + 1) def bose_einstein_occupation(eps, mu, T): """n(ε) = 1/(exp((ε-μ)/kT) - 1)""" kB = 1.38e-23 if mu >= eps: return np.inf return 1 / (np.exp((eps - mu) / (kB * T)) - 1) print(f"\nFermi-Dirac occupation (ε-μ = 0.1 eV, T=300K):") delta_E = 0.1 * 1.6e-19 # 0.1 eV f = fermi_dirac_occupation(delta_E, 0, 300) print(f" f(ε_F + 0.1 eV) = {f:.4f}") ``` ### Ising Model and Phase Transitions ```python import numpy as np class IsingModel: def __init__(self, L, J=1, h=0): """2D Ising model on L×L lattice.""" self.L = L self.J = J self.h = h self.N = L * L self.spins = np.random.choice([-1, 1], size=(L, L)) def energy(self): """Calculate total energy.""" E = 0 for i in range(self.L): for j in range(self.L): # Periodic boundaries right = self.spins[(i+1) % self.L, j] down = self.spins[i, (j+1) % self.L] E -= self.J * self.spins[i, j] * (right + down) E -= self.h * self.spins[i, j] return E / 2 # Each bond counted twice def magnetization(self): """Total magnetization.""" return np.sum(self.spins) def metropolis_step(self, T): """Single Metropolis sweep.""" kB = 1.0 # Units where kB = 1 for _ in range(self.N): # Random spin flip i = np.random.randint(self.L) j = np.random.randint(self.L) # Calculate energy change neighbors = (self.spins[(i-1) % self.L, j] + self.spins[(i+1) % self.L, j] + self.spins[i, (j-1) % self.L] + self.spins[i, (j+1) % self.L]) dE = 2 * self.spins[i, j] * (self.J * neighbors + self.h) # Metropolis criterion if dE < 0 or np.random.random() < np.exp(-dE / (kB * T)): self.spins[i, j] *= -1 def equilibrate(self, T, n_sweeps=1000): """Equilibrate at given temperature.""" for _ in range(n_sweeps): self.metropolis_step(T) def simulate(self, T, n_sweeps=10000): """Run simulation and return observables.""" self.equilibrate(T, n_sweeps) mags = [] eners = [] for _ in range(n_sweeps): self.metropolis_step(T) mags.append(self.magnetization()) eners.append(self.energy()) return np.mean(mags)/self.N, np.std(mags)/self.N, np.mean(eners)/self.N, np.std(eners)/self.N # Critical temperature for 2D Ising (Onsager solution: kB Tc/J ≈ 2.269) L = 20 ising = IsingModel(L) print("2D Ising model simulation:") temperatures = [1.0, 2.0, 2.2, 2.269, 2.3, 2.5, 3.0] for T in temperatures: m, dm, E, dE = ising.simulate(T, n_sweeps=5000) print(f" T = {T:.2f}: m = {m:.4f} ± {dm:.4f}, E = {E/L**2:.4f} ± {dE/L**2:.4f}") # Binder cumulant def binder_cumulant(mags): """U = 1 - ⟨m⁴⟩/(3⟨m²⟩²)""" m2 = np.mean(np.array(mags)**2) m4 = np.mean(np.array(mags)**4) return 1 - m4 / (3 * m2**2) print(f"\nBinder cumulant at Tc ≈ 2.27:") for T in [2.2, 2.269, 2.3]: ising2 = IsingModel(20) ising2.equilibrate(T, 5000) mags = [ising2.magnetization() for _ in range(5000)] U = binder_cumulant(mags) print(f" U(T={T}) = {U:.4f}") ``` ### Mean-Field Theory ```python import numpy as np def mean_field_ising(T, J, h=0, max_iter=100, tol=1e-6): """Mean-field solution of Ising model.""" kB = 1.0 # Units where kB = 1 beta = 1 / T # Initial magnetization m = 0.1 for iteration in range(max_iter): m_new = np.tanh(beta * (J * 6 * m + h)) # z=6 for 3D if abs(m_new - m) < tol: break m = m_new return m def solve_self_consistent(T, J): """Solve MF equations numerically.""" from scipy.optimize import brentq def equation(m): kB = 1.0 return m - np.tanh(m * J * 6 / (kB * T)) # Below Tc, look for non-zero solution Tc_mf = J * 6 # Mean-field Tc if T > Tc_mf: return 0 # Find root try: m = brentq(equation, 1e-6, 1.0) return m except: return 0 # Critical exponents Tc_mf = 6.0 # For J=1 print("Mean-field Ising critical behavior:") print(f" Mean-field Tc = {Tc_mf:.3f} J/kB") # Below Tc, magnetization follows m ∝ (Tc - T)^β with β = 0.5 for T in [0.5, 1.0, 2.0, 4.0, 5.0, 5.5, 5.9, 5.99]: m = solve_self_consistent(T, 1.0) if T < Tc_mf: beta_mf = np.log(m) / np.log((Tc_mf - T) / Tc_mf) print(f" T = {T:.2f}: m = {m:.4f}, effective β = {beta_mf:.3f}") else: print(f" T = {T:.2f}: m = 0") # Landau theory def landau_free_energy(m, a, b, T): """F(m) = a(T-Tc)m² + bm⁴ + hm""" return a * (T - Tc_mf) * m**2 + b * m**4 # Susceptibility def susceptibility(T, Tc, chi0=1): """χ ∝ |T-Tc|^(-γ) with γ = 1 (MF)""" if T > Tc: return chi0 / (T - Tc) return chi0 / (T_c - T) print(f"\nMean-field susceptibility divergence:") for T in [2.0, 1.1, 1.01, 1.001]: chi = susceptibility(T, 1.0) print(f" T = {T:.3f}: χ = {chi:.1f}") ``` ### Monte Carlo Methods ```python import numpy as np def metropolis_sampling(pdf, proposal_std, n_samples, x0): """Metropolis-Hastings sampling from arbitrary PDF.""" samples = [x0] x = x0 for _ in range(n_samples): x_proposed = x + np.random.normal(0, proposal_std) alpha = pdf(x_proposed) / pdf(x) if np.random.random() < alpha: x = x_proposed samples.append(x) return np.array(samples) def wolff_cluster_update(spins, J, T): """Wolff single-cluster Monte Carlo update.""" L = len(spins) visited = np.zeros(L, dtype=bool) # Start new cluster site = np.random.randint(L) cluster = [site] visited[site] = True S_cluster = spins[site] # Grow cluster i = 0 while i < len(cluster): current = cluster[i] P_add = 1 - np.exp(-2 * J / T)