Ising Model Simulation

What the simulation shows

Each site of an L × L grid holds a spin s = +1 or −1. Neighboring spins that agree lower the energy, so at low temperature the lattice orders into one color; heat randomizes it. The lab runs the Metropolis algorithm and plots the magnetization per spin m = Σs_i/L² and the energy per spin as it goes.

• Below T_c (try T = 2): the lattice stays magnetized, with small islands of flipped spins flickering in and out. The measured |m| sits close to Onsager's exact curve.
• Near T_c (try 2.27 with a random start): clusters of every size appear, from single spins to the whole lattice. The magnetization swings wildly and the simulation slows down, because large clusters take a long time to rearrange.
• Above T_c (try 3.5): a salt-and-pepper pattern with short-lived small clusters, and m fluctuating around 0.
• A quench (random start at T = 1.5): domains of both colors form at once and coarsen slowly. Sometimes two stripes wrap around the periodic lattice and survive for thousands of sweeps.

Each finished run adds a point to the phase diagram, so moving the temperature slider step by step traces the magnetization curve. With a finite lattice the points round off the sharp drop at T_c; larger L makes them follow it more closely.

The model and Onsager’s exact solution

With coupling J = 1 and no external field, the energy of a configuration is a sum over nearest-neighbor pairs, and at temperature T (in units of J/k_B) configurations occur with Boltzmann probability:

Ernst Ising solved the one-dimensional chain in 1925 and found no phase transition. In 1944 Lars Onsager solved the two-dimensional square lattice exactly. The critical temperature is

and below it the spontaneous magnetization of the infinite lattice is

a formula Onsager announced in 1949 and C. N. Yang first derived in print in 1952. Near T_c it vanishes like (T_c − T)^1/8, a steep drop no mean-field theory predicts. The lab also compares the energy per spin with Onsager's exact expression, which involves a complete elliptic integral; at T_c it equals −√2 ≈ −1.414. Both exact results are for an infinite lattice; a finite one deviates most near T_c, where the correlation length exceeds L.

The Metropolis algorithm

1. Pick a site i uniformly at random.
2. Compute the energy change of flipping it: ΔE = 2s_i × (sum of its four neighbors), one of −8, −4, 0, 4, 8.
3. If ΔE ≤ 0, flip. Otherwise flip with probability e^−ΔE/T, using the two factors e^−4/T and e^−8/T computed once.
4. Repeat; L² attempts make one sweep. Update the magnetization and energy incrementally with each accepted flip.

The acceptance rule satisfies detailed balance with respect to the Boltzmann distribution, so after enough sweeps the lattice visits configurations with the right probabilities (Metropolis–Hastings generalizes the idea to any distribution). “Enough” is the catch: successive sweeps are correlated, and near T_c the correlation time grows sharply with L (critical slowing down). Cluster algorithms such as Wolff's flip whole clusters at once and avoid most of it. The lab discards the first half of the run as equilibration and averages the second half.

It averages |m|, not m: a finite lattice has no true symmetry breaking, and given enough time it flips between the up and down states, so ⟨m⟩ would be 0 at every temperature. The acceptance rate is a useful diagnostic too; it falls toward 0 at low temperature, where almost every proposed flip costs energy.

Example result

On a 48 × 48 lattice at T = 2 from a cold start (400 sweeps, seed 42), the magnetization averaged over the second half of the run is 0.909 against Onsager’s 0.911, and the energy per spin is −1.741 against the exact −1.746. 8.0% of the proposed flips were accepted; the critical temperature is T_c = 2.269.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator, so a given seed produces the same result in both languages. The excerpts show the core of each.

// Metropolis acceptance factors: flipping spin i changes the energy by
// ΔE = 2 s_i (sum of its four neighbours) ∈ {−8, −4, 0, 4, 8} (units of J).
export function boltzmann(T) {
  return { w4: Math.exp(-4 / T), w8: Math.exp(-8 / T) };
}

// One sweep = N attempted single-spin flips at random sites.
export function sweep(g, w, rng) {
  const { s, right, left, up, down, N } = g;
  for (let k = 0; k < N; k++) {
    const i = rng.int(N);
    const dE = 2 * s[i] * (s[right[i]] + s[left[i]] + s[up[i]] + s[down[i]]);
    // accept downhill moves always, uphill ones with probability exp(−ΔE/T)
    if (dE <= 0 || rng.uniform() < (dE === 4 ? w.w4 : w.w8)) {
      s[i] = -s[i];
      g.M += 2 * s[i];
      g.E += dE;
      g.accepted++;
    }
  }
  g.attempted += N;
}

def sweep(g, w4, w8, rng):
    """N attempted single-spin flips at random sites (Metropolis)."""
    s, right, left, up, down, N = g["s"], g["right"], g["left"], g["up"], g["down"], g["N"]
    for _ in range(N):
        i = rng.int(N)
        dE = 2 * s[i] * (s[right[i]] + s[left[i]] + s[up[i]] + s[down[i]])
        # accept downhill moves always, uphill ones with probability exp(-dE/T)
        if dE <= 0 or rng.uniform() < (w4 if dE == 4 else w8):
            s[i] = -s[i]
            g["M"] += 2 * s[i]
            g["E"] += dE
            g["accepted"] += 1
    g["attempted"] += N

Run the full model from a terminal:

node simulation.js --size 48 --temperature 2 --seed 42
python simulation.py --size 48 --temperature 2 --seed 42

Download the model

Assumptions & limitations

• Nearest-neighbor coupling, no field. Only the four neighbors interact, with equal strength; real magnets have longer-range and anisotropic interactions and quantum spins.
• Finite lattice with periodic boundaries. The exact results are for an infinite lattice; expect visible finite-size effects near T_c and for small L.
• Dynamics is a sampling device. Metropolis sweeps are not physical time; they are a way to draw configurations from the Boltzmann distribution. Coarsening after a quench is only qualitatively like real kinetics.
• Fixed equilibration rule. The first half of each run is discarded. From a random start below T_c that is often not enough, and the lab says so.
• Statistical error not shown. Correlated samples make the naive standard error too small; a binning or autocorrelation analysis is needed for error bars (see Monte Carlo error).

Related simulations

• OS-005Metropolis–HastingsA Markov chain that samples any density known up to a constant; tune the step size and watch it mix.
• OS-014Markov chainA process that remembers only its current state: stationary distribution, mixing and absorption, live.
• OS-022Game of LifeTwo rules on a grid of cells: gliders, guns, oscillators and soups that settle or keep going.
• OS-037DiffusionThe heat equation solved two ways, by finite differences and by random walkers, against the exact Gaussian.

Background: How Accurate Is a Monte Carlo Estimate? · Pseudo-Random Numbers and Seeds in Simulation

References

1. Ising, E. (1925). Beitrag zur Theorie des Ferromagnetismus. Zeitschrift für Physik, 31, 253–258. doi.org/10.1007/BF02980577
2. Onsager, L. (1944). Crystal statistics. I. A two-dimensional model with an order-disorder transition. Physical Review, 65, 117–149. doi.org/10.1103/PhysRev.65.117
3. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21(6), 1087–1092. doi.org/10.1063/1.1699114
4. Newman, M. E. J. & Barkema, G. T. (1999). Monte Carlo Methods in Statistical Physics. Oxford University Press.