Statistical Mechanics Deep Dives#
Twenty-one worked examples, one per conceptual breakthrough in
Breakthroughs in Statistical Mechanics, each pairing a schematic diagram of
the underlying lattice, chain, or ensemble with a short, runnable
physicskit.statphys script reproducing the phenomenon’s textbook
observable signature. Every reference here has a verified counterpart in
the Examples example gallery; this page exists to
put the physics of each breakthrough next to its diagram and its code, in
one place.
1870s – Boltzmann’s Statistical Entropy and the H-Theorem#
Boltzmann’s \(S = k_B \ln W\) and his H-theorem are two sides of one claim: macroscopic irreversibility is nothing but the overwhelmingly probable behavior of reversible microscopic collisions. The clearest way to see the H-theorem is to watch a real N-body gas relax: start every particle with the same speed (about as far from equilibrium as a fixed total energy allows) and watch the speed distribution spread into the Maxwell-Boltzmann bell curve while \(H = \int f \ln f\, dv\) falls.
Periodic box, N Lennard-Jones particles, all |v| equal at t=0:
+-------------------------------+
| o-> o-> o-> o-> | <- every particle: same speed,
| o-> o-> o-> | random direction (delta-fn
| o-> o-> o-> o-> | speed distribution)
| o-> o-> o-> |
+-------------------------------+
| collisions randomize speeds
v
+-------------------------------+
| o^ o<- o> ov o\ o/ | <- speeds now spread out;
| o> o\ o^ o<- o/ | H(t) has decreased toward
| o<- o/ ov o> o\ | its Maxwell-Boltzmann minimum
+-------------------------------+
physicskit.statphys.chapters.molecular_dynamics.LennardJonesGas
integrates exactly this system with Velocity-Verlet dynamics, and
h_function()
estimates \(H(t)\) from the instantaneous speed histogram at every
recorded step.
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.molecular_dynamics import LennardJonesGas
from physicskit.statphys.visualizers.particle_render import plot_velocity_histogram
# Start every particle with the same speed -- far from Maxwell-Boltzmann.
gas = LennardJonesGas(
n_particles=225, box_size=25.0, initial_velocity_distribution="delta", seed=0
)
history = gas.run(n_steps=4000, steps_per_record=40)
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].plot(history["time"], history["H"])
axes[0].set_xlabel("time")
axes[0].set_ylabel(r"Boltzmann $H(t) = \int f \ln f\, dv$")
axes[0].set_title("H-theorem: monotonic relaxation")
plot_velocity_histogram(gas, ax=axes[1], show_theory=True)
axes[1].set_title("Final speed distribution vs. Maxwell-Boltzmann")
plt.tight_layout()
plt.show()
1902 – Gibbs’s Ensemble Formalism#
The canonical ensemble ties macroscopic response functions directly to microscopic fluctuations. Sweeping a lattice model’s temperature and measuring the variance of its energy and magnetization at each point reproduces the specific heat and susceptibility peaks that every later lattice-model chapter in this package relies on.
................
: Heat bath : T fixed, energy exchanged freely
: (reservoir) :
''''''''|''''''''
| <- fluctuating energy E
+------v------+
| L x L | P(configuration) proportional to
| spin | exp(-beta * E)
| lattice |
+-------------+
Var(E) -> C_v , Var(M) -> chi (fluctuation-dissipation)
import numpy as np
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.ising_lattice import Ising2D
from physicskit.statphys.visualizers.lattice_render import plot_thermodynamics
model = Ising2D(L=32, seed=0)
temperatures = np.linspace(model.T_C + 1.2, max(model.T_C - 1.2, 0.5), 20)
result = model.run_temperature_sweep(
temperatures, n_equil=200, n_measure=300, algorithm="wolff"
)
plot_thermodynamics(result, T_c=model.T_C)
plt.suptitle(r"Canonical-ensemble $C_v$ and $\chi$ from energy/magnetization variance")
plt.show()
1905 – Einstein’s Explanation of Brownian Motion#
An ensemble of independent random walkers, each taking many small independent kicks, has a mean squared displacement growing linearly in time – the diffusion law Einstein derived to explain Brownian motion and that Perrin used to measure Avogadro’s number.
Single walker trajectory (2D): Ensemble MSD(t):
MSD
* <- start | . '
\ | .'
\___ | .'
\ | .' slope = 2 d D
____/ | .'
/ | .'
* <- t=1000 +------------------ t
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.random_walk import RandomWalk
from physicskit.statphys.visualizers.random_walk_render import plot_msd, plot_trajectories_2d
walk = RandomWalk(n_walkers=500, n_steps=2000, dim=2, kind="gaussian", seed=0)
D = walk.diffusion_coefficient()
print(f"Estimated diffusion coefficient D = {D:.4f}")
fig, axes = plt.subplots(1, 2, figsize=(10, 4.5))
plot_trajectories_2d(walk, ax=axes[0], n_show=25)
plot_msd(walk, ax=axes[1], show_theory=True)
plt.tight_layout()
plt.show()
1907 – The Ehrenfest Urn Model#
The Ehrenfest urn is exactly reversible and does recur, yet a system started far from equilibrium (all \(N\) balls in one box) relaxes to near-equal occupancy almost immediately, while the recurrence time is astronomically longer than the simulation could ever run.
t = 0 (maximally non-equilibrium) t large (near equilibrium)
+-----------+ +-----------+ +-----------+ +-----------+
|o o o o o o| | | |o o o o o | |o o o o o|
|o o o o o o| | | --> | o o o | |o o o o |
|o o o o o o| | | |o o o o | | o o o o|
+-----------+ +-----------+ +-----------+ +-----------+
n_left = N n_right = 0 n_left ~ N/2 n_right ~ N/2
one ball, chosen uniformly at random, crosses over at every time step
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.ehrenfest_urn import EhrenfestUrn
from physicskit.statphys.visualizers.urn_render import plot_ehrenfest_history
urn = EhrenfestUrn(n_balls=200, seed=0) # all 200 balls start in the left box
history = urn.run(n_steps=6000)
plot_ehrenfest_history(history, n_balls=200)
plt.suptitle("Reversible dynamics, irreversible-looking relaxation to n_left ~ N/2")
plt.show()
1920-1944 – The Ising Model and Onsager’s Exact Solution#
Onsager’s exact solution locates the 2D Ising model’s continuous phase transition at \(T_C = 2J/(k_B \ln(1+\sqrt2)) \approx 2.269\, J/k_B\). Sweeping temperature across \(T_C\) and snapshotting the lattice below, at, and above it makes the transition – from an ordered, mostly-aligned phase to a disordered one – directly visible.
T << T_C (ordered) T ~ T_C (critical) T >> T_C (disordered)
^ ^ ^ ^ ^ ^ ^ ^ ^ v ^ ^ v v ^ ^ ^ v v ^ v ^ v v
^ ^ ^ ^ ^ ^ ^ ^ v ^ ^ v v ^ v ^ v ^ v v ^ v ^ ^
^ ^ ^ ^ ^ ^ ^ ^ ^ ^ v v v ^ ^ ^ v v ^ ^ v ^ v v
^ ^ ^ ^ ^ ^ ^ ^ v v ^ ^ ^ v v ^ ^ v ^ v v ^ ^ v
bulk magnetized scale-invariant zero net
(spontaneous domains, all sizes magnetization
symmetry breaking)
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.ising_lattice import Ising2D
from physicskit.statphys.visualizers.lattice_render import plot_spin_grid
print(f"Onsager T_C = {Ising2D().T_C:.4f} J/kB")
fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))
for ax, dT in zip(axes, [-1.0, 0.0, 1.0]):
model = Ising2D(L=64, seed=1)
T = model.T_C + dT
model.sweep(beta=1.0 / T, algorithm="wolff", n_sweeps=200)
plot_spin_grid(model.spins, ax=ax, title=f"T = {T:.2f}")
plt.tight_layout()
plt.show()
1937 – Landau Mean-Field Theory#
Expanding the free energy in powers of an order parameter \(m\), with no microscopic input beyond symmetry, already reproduces spontaneous symmetry breaking: one minimum above \(T_C\), a split double well below it.
F(m) F(m)
| \ / | \ /
| \ / T > T_C | \ /\ / T < T_C
| \ / (m=0) | \/ \/ (m = +-m_eq)
| \/ | /\ /\
+------------- m +------------- m
0 -m_eq 0 +m_eq
import numpy as np
import matplotlib.pyplot as plt
from physicskit.statphys.utils.landau_theory import (
landau_free_energy, landau_equilibrium_magnetization,
)
Tc = 2.0
m = np.linspace(-1.5, 1.5, 300)
fig, ax = plt.subplots(figsize=(6, 4.5))
for T in [Tc + 0.8, Tc, Tc - 0.8]:
ax.plot(m, landau_free_energy(m, T=T, Tc=Tc), label=f"T={T:.1f}")
ax.set_xlabel("order parameter m"); ax.set_ylabel("F(m)"); ax.legend()
plt.tight_layout()
plt.show()
print(f"m(T_C - 0.8) = {landau_equilibrium_magnetization(Tc - 0.8, Tc):.4f} "
f"(mean-field beta=1/2 power law)")
1952 – The Yang-Lee Circle Theorem#
Continuing a finite Ising system’s partition function into the complex fugacity plane, every zero of a ferromagnetic system lands exactly on the unit circle – and a phase transition, in the infinite-size limit, is where these zeros pinch the positive real axis.
Im(z)
| x x
| x x every zero of Z(z) lies
----x------+------x---- Re(z) exactly on |z|=1
| x x (ferromagnetic J > 0)
| x x
import numpy as np
import matplotlib.pyplot as plt
from physicskit.statphys.utils.partition_function import (
ising_partition_polynomial, yang_lee_zeros,
)
g = ising_partition_polynomial(N=18, beta=0.4, J=1.0, periodic=True)
zeros = yang_lee_zeros(g)
print(f"max deviation from |z|=1: {np.max(np.abs(np.abs(zeros) - 1)):.2e}")
theta = np.linspace(0, 2 * np.pi, 200)
plt.figure(figsize=(5, 5))
plt.plot(np.cos(theta), np.sin(theta), "--", color="gray", linewidth=0.8)
plt.scatter(zeros.real, zeros.imag, color="tab:red")
plt.gca().set_aspect("equal")
plt.xlabel("Re(z)"); plt.ylabel("Im(z)")
plt.tight_layout()
plt.show()
1924-1926 – Bose-Einstein and Fermi-Dirac Quantum Statistics#
Bosons pile without limit into a shared ground state below a critical temperature; fermions, forbidden by the Pauli exclusion principle from ever doubling up, fill states one at a time up to a sharp Fermi level as \(T \to 0\).
Bosons, T < T_c: Fermions, T -> 0:
energy energy
| | (empty states)
| o | ----------------
| o | ---------------- <- Fermi level mu
| o | ================ (filled, 1 per
| o o o | ================ state, Pauli-
|__ooooooooooooo__ ground state |__================_ excluded above)
macroscopic condensate hard step at T=0
import numpy as np
import matplotlib.pyplot as plt
from physicskit.statphys.utils.thermodynamics import (
bose_einstein_occupation, fermi_dirac_occupation, bec_condensate_fraction,
)
energy = np.linspace(0.01, 5.0, 400)
fig, axes = plt.subplots(1, 2, figsize=(10, 4.5))
for T, style in zip([0.5, 1.0, 2.0], ["-", "--", ":"]):
axes[0].plot(energy, fermi_dirac_occupation(energy, mu=0.0, temperature=T),
style, color="steelblue", label=f"Fermi-Dirac, T={T}")
axes[0].plot(energy, bose_einstein_occupation(energy, mu=-0.01, temperature=T),
style, color="crimson", label=f"Bose-Einstein, T={T}")
axes[0].set_ylim(0, 5)
axes[0].set_xlabel(r"energy $\varepsilon$"); axes[0].legend(fontsize=8)
T = np.linspace(0.0, 2.0, 200)
axes[1].plot(T, bec_condensate_fraction(T, critical_temperature=1.0))
axes[1].axvline(1.0, color="k", linestyle="--", linewidth=1, label=r"$T_c$")
axes[1].set_xlabel(r"$T/T_c$"); axes[1].set_ylabel(r"condensate fraction $N_0/N$")
axes[1].legend()
plt.tight_layout()
plt.show()
1952 – The Potts Model#
The q-state Potts model rewards neighboring sites only for matching exactly. On the square lattice, watching the order parameter as temperature sweeps through \(T_C\) shows the transition sharpen from a smooth, continuous drop at small \(q\) into an abrupt jump once \(q > 4\).
q = 3 Potts configuration (letters = discrete states, not spin signs):
A B B C A A H = -J * sum_<i,j> delta(s_i, s_j)
A A C C B B (energy lowered only by *exact* neighbor matches,
C A A B B C not by a continuous alignment as in Ising/XY)
B B A A C A
C C A B B A
import numpy as np
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.ising_lattice import PottsModel2D
fig, ax = plt.subplots(figsize=(6, 4.5))
for q in [2, 4, 8]:
model = PottsModel2D(L=24, q=q, seed=0)
temperatures = np.linspace(model.T_C + 0.8, max(model.T_C - 0.8, 0.1), 15)
result = model.run_temperature_sweep(temperatures, n_equil=150, n_measure=150)
ax.plot(result["T"] / model.T_C, result["m"], marker="o", label=f"q={q}")
ax.set_xlabel(r"$T / T_C(q)$"); ax.set_ylabel("Potts order parameter $m$")
ax.set_title("Transition sharpens from continuous (q <= 4) to first-order (q > 4)")
ax.legend()
plt.tight_layout()
plt.show()
1953 – The Metropolis Algorithm#
Every lattice chapter in this package is, underneath, the same accept/reject rule: propose a local change, accept downhill moves always and uphill moves with probability \(e^{-\beta \Delta E}\).
current configuration s
|
v
propose s' (flip one randomly chosen spin)
|
v
dE = E(s') - E(s)
| |
dE <= 0 dE > 0
| |
v v
always accept accept with prob. exp(-beta * dE)
| |
+----------+----------+
v
repeat for every site (one sweep)
import numpy as np
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.ising_lattice import Ising2D
model = Ising2D(L=32, seed=0)
beta = 1.0 / model.T_C
energies = [model.energy() / model.n_sites]
for _ in range(200):
model.sweep(beta=beta, algorithm="metropolis", n_sweeps=1)
energies.append(model.energy() / model.n_sites)
plt.figure(figsize=(6, 4))
plt.plot(energies)
plt.xlabel("Metropolis sweep"); plt.ylabel("energy per site")
plt.title(r"Metropolis relaxation toward equilibrium at $T_C$")
plt.tight_layout()
plt.show()
1957 – Molecular Dynamics Simulation#
Alder and Wainwright integrated Newton’s equations directly for hundreds of interacting particles rather than sampling configurations stochastically – the direct-simulation complement to Metropolis Monte Carlo.
+-------------------------------+ <- periodic boundary: a particle
| o o o o | leaving the right edge re-enters
| o o o o o | on the left (and similarly top/
| o o o o o | bottom)
| o o o o o |
| o o o o o | Velocity Verlet: x(t+dt), v(t+dt)
+-------------------------------+ from Lennard-Jones pairwise forces
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.molecular_dynamics import LennardJonesGas
from physicskit.statphys.visualizers.particle_render import plot_particle_snapshot
gas = LennardJonesGas(n_particles=144, box_size=20.0, seed=0)
E0 = gas.total_energy()
gas.step(n_steps=2000)
E1 = gas.total_energy()
print(f"Energy drift over 2000 Verlet steps: {abs(E1 - E0) / abs(E0):.2e} (relative)")
fig, ax = plt.subplots(figsize=(5, 5))
plot_particle_snapshot(gas, ax=ax, color_by_speed=True)
plt.tight_layout()
plt.show()
1957 – Percolation Theory#
Below the percolation threshold \(p_c\), only finite clusters exist; above it, an infinite spanning cluster appears with probability one – a purely geometric phase transition with a sharp threshold and a power-law cluster-size distribution right at \(p_c\).
p < p_c (subcritical) p = p_c (critical) p > p_c (spanning)
. X X . X . X . X X . X X . X X X X X X X X X X
X . X . . X . X X . X X . X . X X X X X X X X X
. X . X X . X . . X X . X X X . X X X X X X X X
X . X . X . X . X . X X . X X X X X X X X X X X
(small clusters only) (spans, fractal) (fully occupied)
import numpy as np
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.percolation import Percolation2D
from physicskit.statphys.visualizers.lattice_render import (
plot_percolation_clusters, plot_cluster_size_distribution,
)
perc = Percolation2D(L=80, seed=0)
p_values = np.linspace(0.3, 0.85, 25)
p, P_span = perc.spanning_probability(p_values, n_trials=60)
fig, axes = plt.subplots(1, 3, figsize=(13, 4))
axes[0].plot(p, P_span, marker="o")
axes[0].axvline(perc.p_c, color="k", linestyle="--", label=r"$p_c$")
axes[0].set_xlabel("p"); axes[0].set_ylabel(r"$P_{\rm span}(p)$"); axes[0].legend()
perc.generate(perc.p_c)
plot_percolation_clusters(perc, ax=axes[1])
plot_cluster_size_distribution(perc.cluster_size_distribution(), ax=axes[2])
plt.tight_layout()
plt.show()
1966-1971 – Kadanoff / Wilson Renormalization Group#
Coarse-graining a critical (\(T = T_C\)) configuration by repeatedly replacing \(2 \times 2\) blocks with a single majority-rule spin produces a self-similar sequence of lattices; away from \(T_C\), the same procedure flows visibly toward a fully ordered or fully disordered fixed point.
Fine lattice L 2x2 majority rule Coarse L/2
+ + - - + -
+ + - - each 2x2 block -> sign(sum of block) ->
- - + + - +
- - + +
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.renormalization import BlockSpinRG
from physicskit.statphys.visualizers.rg_render import plot_rg_flow
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
for ax, (label, T) in zip(axes, [("T << T_C", 1.0), ("T = T_C", None), ("T >> T_C", 6.0)]):
rg = BlockSpinRG(L=64, T=T, seed=0)
grids = rg.iterate()
order_params = [BlockSpinRG.order_parameter(g) for g in grids]
ax.plot(order_params, marker="o")
ax.set_title(label); ax.set_xlabel("coarse-graining step")
ax.set_ylabel("|<s>|"); ax.set_ylim(-0.05, 1.05)
plt.tight_layout()
plt.show()
1973 – Vortex Unbinding in the 2D XY Model#
Below the Kosterlitz-Thouless temperature, vortices remain bound in tight pairs of opposite topological charge; above it, thermal fluctuations tear these pairs into a free “vortex plasma” – a phase transition organized by topology rather than a local order parameter.
Bound vortex-antivortex pair (T < T_KT):
^ -> charge +1 (vortex): spins wind
\ \ counterclockwise once around the core
(+1) ->
/ \ charge -1 (antivortex): spins wind
v <- (-1) <- clockwise once around the core
/ \
<- v tightly bound below T_KT; unbind into a
free plasma above T_KT
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.ising_lattice import XYModel2D
from physicskit.statphys.visualizers.vortex_render import plot_vortices
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
for ax, dT in zip(axes, [-0.3, 0.5]):
model = XYModel2D(L=48, seed=0)
T = model.T_KT + dT
model.sweep(beta=1.0 / T, n_sweeps=300)
nv, na = model.vortex_count()
plot_vortices(model.theta, ax=ax)
ax.set_title(f"T = {T:.2f} ({'below' if dT < 0 else 'above'} $T_{{KT}}$): "
f"{nv} vortices, {na} antivortices")
plt.tight_layout()
plt.show()
1975 – The Edwards-Anderson Spin Glass#
A generic plaquette with an odd number of antiferromagnetic bonds cannot satisfy all four of its bonds at once – it is frustrated. The Edwards-Anderson order parameter, the overlap between two replicas sharing the same disorder, stays nonzero as long as both replicas freeze into the same disorder-selected (if disordered) pattern.
o --+J-- o --+J-- o Plaquette bonds (clockwise): +J, -J, +J, -J
| | | product of four bonds < 0 ==> FRUSTRATED:
-J +J -J no single spin assignment on its four sites
| | | satisfies all four bonds simultaneously
o --+J-- o --+J-- o
import numpy as np
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.spin_glass import EdwardsAndersonSpinGlass2D
from physicskit.statphys.visualizers.lattice_render import plot_spin_grid
model = EdwardsAndersonSpinGlass2D(L=32, seed=0)
print(f"Frustrated plaquette fraction: {model.frustration_density():.3f} (expect ~0.5)")
temperatures = np.linspace(0.3, 3.0, 10)
q2_values = [
model.edwards_anderson_order_parameter(beta=1.0 / T, n_equil=150, n_measure=150)
for T in temperatures
]
fig, axes = plt.subplots(1, 2, figsize=(10, 4.5))
model.sweep(beta=2.0, n_sweeps=500)
plot_spin_grid(model.spins, ax=axes[0], title="Frozen configuration, T=0.5 J/kB")
axes[1].plot(temperatures, q2_values, marker="o")
axes[1].set_xlabel("Temperature"); axes[1].set_ylabel(r"$\langle q^2 \rangle$")
plt.tight_layout()
plt.show()
1975-1980 – The Sherrington-Kirkpatrick Model and Parisi’s Replica Symmetry Breaking#
Every spin coupled to every other spin by an independent Gaussian bond. Parisi’s replica-symmetry-breaking solution predicts the pooled replica-overlap distribution P(q) is a single narrow peak above the transition, but broad and structured below it.
P(q) P(q)
| ___ | __
| / \ T > T_SG | _/ \_ T < T_SG
| / \ (paramagnetic) | / \ (many pure states,
+--+---+---+-- q +---------+-- q broad, non-Gaussian)
-1 0 1 -1 0 1
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.spin_glass import SherringtonKirkpatrick
from physicskit.statphys.visualizers.spin_glass_render import plot_overlap_distribution
fig, ax = plt.subplots(figsize=(6, 4))
for beta, label in [(0.3, "T=3.3 J/kB (paramagnetic)"), (2.0, "T=0.5 J/kB (spin-glass)")]:
model = SherringtonKirkpatrick(N=64, seed=0)
samples = model.overlap_distribution(beta=beta, n_disorder=25, n_equil=300, n_measure=60)
plot_overlap_distribution(samples, ax=ax, label=label)
plt.tight_layout()
plt.show()
1981 – The Binder Cumulant and Finite-Size Scaling#
The fourth-order cumulant \(U_4 = 1 - \langle M^4\rangle / (3\langle M^2\rangle^2)\) is dimensionless, so at \(T_C\) its curves for different lattice sizes \(L\) all cross at (very nearly) the same temperature – letting \(T_C\) be pinned down without ever simulating an infinite system.
U4
| L=8 \ / L=8
| \ X / X marks the common crossing:
| \ | / the finite-size-independent
| \ | / estimate of T_C
| L=16 / | \ L=16
| / | \
+------------------------------------- T
T_C
import numpy as np
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.ising_lattice import Ising2D
from physicskit.statphys.utils.thermodynamics import binder_cumulant
from physicskit.statphys.utils.finite_size_scaling import binder_cumulant_crossing
temperatures = np.linspace(1.8, 2.8, 12)
U4_by_L = {}
for L in [8, 16, 24]:
model = Ising2D(L=L, seed=0)
U4 = []
for T in temperatures:
beta = 1.0 / T
model.sweep(beta=beta, algorithm="wolff", n_sweeps=100) # equilibrate
mags = []
for _ in range(200):
model.sweep(beta=beta, algorithm="wolff", n_sweeps=1)
mags.append(model.magnetization())
U4.append(binder_cumulant(np.array(mags)))
U4_by_L[L] = np.array(U4)
T_c_estimate = binder_cumulant_crossing(temperatures, U4_by_L)
print(f"Binder-cumulant crossing estimate of T_C: {T_c_estimate:.3f} "
f"(Onsager: {Ising2D().T_C:.3f})")
plt.figure(figsize=(6, 4.5))
for L, U4 in U4_by_L.items():
plt.plot(temperatures, U4, marker="o", label=f"L={L}")
plt.axvline(T_c_estimate, color="k", linestyle="--", label=r"$T_C$ estimate")
plt.xlabel("T"); plt.ylabel(r"$U_4$"); plt.legend()
plt.tight_layout()
plt.show()
1986 – The Kardar-Parisi-Zhang Equation#
A local growth rule with no adjustable surface-tension or coupling constant at all – the Restricted Solid-On-Solid model – is enough to place a growing interface in the KPZ universality class: early-time width growth \(w \sim t^{1/3}\) crossing over to a size-limited saturation \(w_{\text{sat}} \sim L^{1/2}\).
height log w
| . . | _________ w_sat ~ L^0.5
| . . .. . rough, | _/
| .. . . uncorrelated | _/ w ~ t^(1/3)
| . .. noise + local | _/
+------------------- x +---------------------- log t
growing interface growth then saturation
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.kpz_growth import KPZInterface
from physicskit.statphys.utils.finite_size_scaling import power_law_exponent
from physicskit.statphys.visualizers.kpz_render import plot_width_growth
interface = KPZInterface(L=400, seed=0)
times, widths = interface.run_growth_curve(t_max=1500, n_points=30)
beta_fit, _ = power_law_exponent(times, widths)
print(f"fitted growth exponent beta = {beta_fit:.3f} (KPZ theory: 1/3)")
plot_width_growth(times, widths, beta=1.0 / 3.0)
plt.tight_layout()
plt.show()
1987 – Self-Organized Criticality and the BTW Sandpile#
Grains added one at a time tune the pile itself to a critical state, with no parameter fine-tuning: avalanche sizes follow a power law with no characteristic scale, the geometric signature of self-organized criticality.
1 2 1 3 2
2 4*2 1 3 site (1,1) reaches height 4 = threshold -> topples,
1 3 1 2 1 sending one grain to each of its four neighbors,
3 2 4 1 2 which may push a neighbor over threshold too --
2 1 2 3 1 an avalanche, size = total topplings triggered
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.sandpile import BTWSandpile
from physicskit.statphys.visualizers.sandpile_render import (
plot_sandpile_heights, plot_avalanche_size_distribution,
)
pile = BTWSandpile(L=64, seed=0)
sizes = pile.run(n_grains=20000, warmup=20000) # warmup reaches the SOC state
fig, axes = plt.subplots(1, 2, figsize=(10, 4.5))
plot_sandpile_heights(pile, ax=axes[0])
plot_avalanche_size_distribution(sizes, ax=axes[1])
plt.tight_layout()
plt.show()
1989 – The Wolff Cluster Algorithm#
A Wolff update grows one cluster of same-sign spins from a random seed – each bond added independently with probability \(1 - e^{-2\beta J}\) – and flips the whole cluster in a single global Monte Carlo move, sidestepping the critical slowing down that limits single-spin Metropolis dynamics right at \(T_C\).
+ + + + + + [+ + +]+ cluster grown from a random seed
+ + + + + --> +[+ + + +] (bracketed): each same-sign neighbor
+ - + + + +[+]- + + bond joins with prob. 1-exp(-2*beta*J)
+ + + - + + + + - + then the entire cluster is flipped
as ONE Monte Carlo move
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.ising_lattice import Ising2D
model = Ising2D(L=64, seed=0)
beta = 1.0 / model.T_C # right at T_C: hardest regime for single-spin dynamics
cluster_sizes = [model.sweep(beta=beta, algorithm="wolff", n_sweeps=1) for _ in range(300)]
plt.figure(figsize=(6, 4))
plt.plot(cluster_sizes)
plt.axhline(model.n_sites, color="k", linestyle="--", linewidth=1, label="full lattice")
plt.xlabel("Wolff update"); plt.ylabel("cluster size")
plt.title(r"Wolff cluster sizes at $T_C$: a global, not local, Monte Carlo move")
plt.legend()
plt.tight_layout()
plt.show()
1997 – The Jarzynski Equality#
Averaging \(e^{-\beta W}\), not \(W\), over repeated realizations of a nonequilibrium protocol recovers the exact equilibrium free energy difference – however fast or dissipative each individual realization is.
density
| <W> <W> sits above the true Delta F
| v (second law), but the exponential
| _ _|_ average of the SAME work samples
| / \_/ \___ Delta F lands exactly back on Delta F
+------|----------- W (Jensen's inequality, saturated)
Delta F
import matplotlib.pyplot as plt
from physicskit.statphys.chapters.nonequilibrium_work import JarzynskiHarmonicTrap
from physicskit.statphys.visualizers.jarzynski_render import plot_work_distribution
model = JarzynskiHarmonicTrap(seed=0)
work = model.run_protocol(lambda_0=0.0, lambda_1=2.0, tau=0.3, n_steps=400, n_trajectories=4000)
dF_estimate = model.jarzynski_free_energy_estimate(work)
print(f"<W> = {work.mean():.3f} Jarzynski dF_estimate = {dF_estimate:.3f} (true dF = 0)")
plot_work_distribution(work, dF_true=0.0, dF_estimate=dF_estimate)
plt.tight_layout()
plt.show()
See Also#
Breakthroughs in Statistical Mechanics for the full chronology and the historical context behind each entry above.
Examples for the complete, gallery-rendered example suite these snippets are drawn from.