Statistical Mechanics & Thermodynamics — Boltzmann, Ising Model, Phase Transitions and the Arrow of Time

Why Physics Needs Statistics

Classical mechanics tells you exactly where a single mass will be given its initial position and velocity. But one mole of gas contains 6 × 10²³ molecules. You cannot track them individually. What you can do is ask: given the laws of mechanics, what is the most probable distribution of speeds? What fraction of molecules have enough energy to break a chemical bond? When does a magnetic material spontaneously order itself?

These are questions of statistical mechanics — and the answers turn out to be deeply tied to one concept: entropy, the measure of how many microscopic configurations are compatible with the macroscopic state you observe. The entropy formulation of thermodynamics (Boltzmann's S = k_B ln Ω) was so revolutionary that it was carved on his tombstone.

The six simulations in this guide move through the field from molecular speeds to phase transitions to heat engines. Each one makes a different statistical phenomenon visible and measurable.

Part 1: Molecular Speeds & the Maxwell-Boltzmann Distribution

Why molecules don't all move at the same speed

When you heat a gas, you add kinetic energy — but that energy does not distribute uniformly. Molecules collide continuously, exchanging momentum in every collision. The result is a spread of speeds that follows the Maxwell-Boltzmann distribution, derived by Maxwell in 1860 and given a statistical mechanics foundation by Boltzmann in 1872.

The Maxwell-Boltzmann simulation runs a 2D elastic hard-sphere gas in real time. Molecules start with equal speeds; after a few hundred collisions, the histogram of speeds converges to the theoretical MB curve. Adjust temperature or molecular mass and watch the distribution shift — higher temperature broadens and shifts the peak right; heavier molecules narrow and shift it left.

Probability density of speed v at temperature T for molecules of mass m:

           ┌  m  ┐³/²            ┌   mv²  ┐
f(v) = 4π │ ─── │    v²  exp  − │ ──────  │
           └ 2πk_BT ┘            └  2k_BT ┘

Key statistics:
  Most probable speed:  v_p  = √(2k_BT / m)
  Mean speed:           v̄   = √(8k_BT / πm)
  RMS speed:            v_rms = √(3k_BT / m)  ← enters kinetic pressure P=nk_BT

Boltzmann factor:
  Probability of state with energy E: P(E) ∝ exp(−E / k_BT)
  Partition function: Z = Σ exp(−Eᵢ / k_BT)

The Boltzmann factor exp(−E/k_BT) appears everywhere in physics and chemistry: it governs the fraction of molecules with enough energy to surmount a reaction barrier (Arrhenius equation), the population of atomic energy levels (spectroscopy), and the probability of magnetic spin alignment (next section). Understanding it in the context of gas speeds is the cleanest entry point.

Part 2: Molecular Structure & Lennard-Jones Dynamics

From gas to liquid to solid — the same force law

Molecules do not just collide like billiard balls. At short range, electron clouds repel (Pauli exclusion); at medium range, induced dipole forces attract (London dispersion); beyond a few nanometres the interaction vanishes. The Lennard-Jones 12-6 potential captures this two-regime interaction with one elegant equation.

The Lennard-Jones simulation runs 2D molecular dynamics with Velocity Verlet integration, an Andersen thermostat for temperature control, and periodic boundary conditions. Lower the temperature slider from the gas preset and watch phase transitions happen: disordered gas → dense liquid → crystalline solid, all emerging from the same pairwise potential. The simulation draws bonds between atoms closer than 1.3σ, so the lattice structure of the solid phase is immediately visible.

LJ 12-6 potential between atoms i and j at separation r:

  U(r) = 4ε [ (σ/r)¹² − (σ/r)⁶ ]

  ε  = depth of potential well (energy scale)
  σ  = finite distance where U = 0  (size scale)
  r_min = 2^(1/6) σ ≈ 1.122 σ       (equilibrium separation)

Forces: F = −dU/dr (Lennard-Jones force, applied via Newton's 3rd law)

Velocity Verlet integrator (symplectic, conserves energy):
  x(t+dt) = x(t) + v(t)dt + ½a(t)dt²
  a(t+dt) = F(t+dt) / m
  v(t+dt) = v(t) + ½[a(t) + a(t+dt)]dt

Phase transitions (LJ system):
  Gas:    k_BT/ε ≳ 1.3,  ρσ² ≪ 1
  Liquid: k_BT/ε ≈ 0.7–1.3
  Solid:  k_BT/ε ≲ 0.7,  ρσ² ≈ 0.85

Part 3: Ferromagnetism & the Ising Model

When a phase transition comes from nothing but neighbour interactions

The Ising model is arguably the most studied model in all of statistical mechanics. Each site on a lattice carries a spin (+1 or −1), and the energy of the system is lowered when neighbouring spins align. At high temperature, thermal fluctuations dominate and the spins point randomly in all directions — zero net magnetisation. Below the Curie temperature T_C, the exchange interaction wins: domains form and the material spontaneously magnetises.

The Ising Model simulation uses the Metropolis algorithm to sample the Boltzmann distribution of the spin lattice. Run it at high temperature: a chaotic salt-and-pepper pattern. Gradually lower the temperature past T_C and watch spontaneous symmetry breaking — one spin orientation wins and domains grow. At exactly T_C, the critical point, domain structure appears at every length scale (fractal geometry — this is a second-order phase transition).

Hamiltonian:
  H = −J Σ⟨i,j⟩ sᵢsⱼ − B Σᵢ sᵢ
  J > 0: ferromagnetic coupling   B: external field
  ⟨i,j⟩: sum over nearest-neighbour pairs

Metropolis MCMC step:
  1. Pick a spin sᵢ at random
  2. Compute ΔE = 2J·sᵢ·Σneighbour sⱼ + 2B·sᵢ
  3. If ΔE < 0: flip (energy decreases) — always accept
     Else: flip with probability exp(−ΔE / k_BT)

Order parameter and critical point (2D square lattice):
  Curie temperature: k_BT_C / J ≈ 2.269  (Onsager exact solution)
  Magnetisation near T_C: m ∝ (T_C − T)^β,  β = 1/8 (2D Ising)
  Correlation length:      ξ ∝ |T − T_C|^−ν, ν = 1  (diverges at T_C!)

Part 4: Blackbody Radiation & the Quantum Bridge

Where classical thermodynamics ran into trouble — and Planck fixed it

Every hot object emits electromagnetic radiation. Classical thermodynamics (the Rayleigh-Jeans law) predicted that an ideal blackbody should radiate infinite power — the "ultraviolet catastrophe." Max Planck resolved the crisis in 1900 by quantising energy: E = nhν. The resulting Planck distribution matches observation exactly.

The Blackbody Radiation simulation plots the spectral radiance from 100 K to 30 000 K. Star presets cover the full range: cool red dwarfs through Sun-like G-stars to blue-white O-stars. Watch Wien's displacement law in action — the peak wavelength shifts as 1/T — and observe Stefan-Boltzmann total power growing as T⁴.

Planck spectral radiance (power per area per wavelength per steradian):

         2hc²      1
B(λ,T) = ────  · ──────────────
          λ⁵    exp(hc/λk_BT)−1

Wien's displacement law: λ_max · T = 2.898 × 10⁻³ m·K
  → Sun (T≈5778 K): λ_max ≈ 502 nm (green/yellow)
  → Human body (T≈310 K): λ_max ≈ 9.3 μm (mid-infrared)

Stefan-Boltzmann total radiated power (per unit area):
  P = σ T⁴      σ = 5.67 × 10⁻⁸ W/(m²·K⁴)

Classical Rayleigh-Jeans (fails at short λ):
  B_RJ(λ,T) = 2ck_BT / λ⁴     (diverges as λ → 0)

Part 5: The Carnot Cycle & Thermodynamic Limits

The best any heat engine can ever do

The second law of thermodynamics imposes a hard limit on how efficiently any heat engine can convert thermal energy to work. Carnot's theorem (1824) shows that the maximum efficiency depends only on the temperatures of the hot and cold reservoirs — not on the working fluid, not on the engineering details.

The Carnot Cycle simulation animates all four reversible steps: isothermal expansion (absorbing Q_H from the hot reservoir), adiabatic expansion (cooling the gas), isothermal compression (rejecting Q_C to the cold reservoir), and adiabatic compression back to the start. The P-V diagram is drawn in real time; the area enclosed equals the net work W done per cycle.

Four Carnot steps on an ideal gas:
  1. Isothermal expansion  (T_H): ΔU=0; Q_H = nRT_H ln(V₂/V₁) = W₁
  2. Adiabatic expansion:  Q=0;  T_H → T_C with TV^(γ-1) = const
  3. Isothermal compression (T_C): Q_C = −nRT_C ln(V₄/V₃) (heat rejected)
  4. Adiabatic compression:        T_C → T_H

Net work per cycle:  W = Q_H − Q_C

Carnot efficiency (theoretical maximum):
  η_C = W / Q_H = 1 − T_C / T_H       (temperatures in Kelvin)

Example — coal power plant:
  T_H ≈ 800 K,  T_C ≈ 300 K  →  η_C = 1 − 300/800 ≈ 62.5%
  Real plant efficiency: ~35–45% (irreversibility, friction, heat loss)

Second law entropy statement:
  ΔS_universe = ΔS_system + ΔS_surroundings ≥ 0
  Equality holds only for reversible (Carnot) processes

Part 6: Brownian Motion & the Einstein Connection

Visible proof of molecules

In 1827 Robert Brown observed pollen grains jiggling erratically under a microscope. He could not explain it. In 1905, the same year that defined modern physics, Albert Einstein published a quantitative theory: the pollen's random walk is driven by thermal collisions with water molecules too small to see individually. Jean Perrin confirmed the theory experimentally in 1908, providing direct evidence for the existence of atoms.

The Brownian Motion simulation implements the Langevin equation: a large tracer particle subject to both viscous drag and random thermal kicks. The mean-squared displacement ⟨Δr²⟩ grows linearly with time (diffusive regime), with slope 4D in 2D. The simulation plots the MSD in real time and overlays Einstein's theoretical prediction. Change temperature or tracer radius and watch D respond as k_BT/(6πηr).

Stokes-Einstein diffusion coefficient:
  D = k_BT / (6πηr)
  k_B = 1.38 × 10⁻²³ J/K  (Boltzmann constant)
  η  = fluid dynamic viscosity  (Pa·s)
  r  = tracer particle radius  (m)

Mean-squared displacement (2D):
  ⟨Δr²(t)⟩ = 4Dt      (slope gives D experimentally)

Langevin equation (overdamped):
  m·dv/dt = −6πηr·v + ξ(t)
  ⟨ξ(t)⟩ = 0,  ⟨ξ(t)ξ(t′)⟩ = 2·6πηr·k_BT·δ(t−t′)

Perrin's measurement (1908):
  Measured D of gamboge resin spheres in water
  → derived k_B → first experimental measurement of Avogadro's number

Statistical Mechanics Across the Simulation Collection

The same Boltzmann factor and entropy arguments that appear in this guide connect to many other categories:

• Maxwell-Boltzmann + Reaction Kinetics — the Arrhenius equation k = A·exp(−E_a/k_BT) is a direct application of the Boltzmann factor to chemical rate constants.
• Crystal Growth + Crystal Diffusion — solid-state diffusion is thermally activated via the same Boltzmann factor; grain boundary energy and nucleation theory are statistical mechanics problems.
• Big Bang Nucleosynthesis — the abundance of hydrogen and helium from the first minutes of the universe is determined by thermal equilibrium rates and Boltzmann factors at T ~ 10¹⁰ K.
• Debye Shielding — the Debye length over which a charge is screened in a plasma follows from the Boltzmann distribution of electrons in the electrostatic potential: n_e ∝ exp(eφ/k_BT).

Algorithms & Methods in This Collection