Notebooks · posts/metropolis-ising-against-onsager

Metropolis sampling finds Onsager's magnetization

Adam Zettel · September 25, 2026 · written with an agent

A 16 by 16 grid of spins, each flipped at random by the Metropolis rule, settles at $T = 1.8$ to an average magnetization of $0.9567 \pm 0.0013$ . The exact infinite-lattice answer, from Onsager and Yang, is $0.95686$ . The critical temperature $2/ ln (1 + 2) = 2.26919$ is where that formula reaches zero.

M (T) = (1 - sinh^{- 4} (2 J / k_{B} T))^{1/8}, T < T_{c} = \frac{2 J}{k _{B} ln ( 1 + 2 )}

The test's tolerance is three standard errors, measured from the run itself by batch means: 20 batches of 60 sweeps each. A separate Metropolis code (posts/refs/metropolis-ising-against-onsager.py, 64 chains of 20000 sweeps) supports that choice three ways. It finds the correlation time of $∣ m ∣$ to be 1.3 sweeps, so the batches are independent. It predicts a standard error of 0.00131 for a run of this length, and the run measures 0.00130. And it finds that the 16 by 16 lattice's average sits $1.1 \times 1 0^{- 5}$ below the infinite lattice's, a hundred times smaller than the statistical error, so the finite size does not need an allowance of its own at this temperature.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

The model

The Ising model is a grid of little magnets, each pointing up ( $s = + 1$ ) or down ( $s = - 1$ ). It shows how a material with only short-range forces ends up magnetized on its own, and why only below a certain temperature. Each spin interacts with its four neighbors, and the energy of an arrangement is

E = - J ⟨ ij ⟩ \sum s_{i} s_{j},

summed over neighboring pairs, so two aligned neighbors lower the energy by $J$ and two opposed ones raise it by $J$ . The program works in units where $J = 1$ and Boltzmann's constant $k_{B} = 1$ , so temperature is measured in units of $J / k_{B}$ .

At temperature $T$ , an arrangement with energy $E$ occurs with probability proportional to $e^{- E / T}$ . The quantity of interest is the magnetization $m = \frac{1}{N} \sum_{i} s_{i}$ , the average spin. Summing $e^{- E / T}$ over all $2^{256}$ arrangements of a 16 by 16 grid is impossible, so the program samples arrangements instead.

The Metropolis rule

Visit each spin in turn and propose flipping it. The flip changes the energy by $Δ E = 2 s_{i} \sum_{neighbors} s_{j}$ , which is one of $- 8, - 4, 0, 4, 8$ . Accept the flip if $Δ E \leq 0$ ; otherwise accept it with probability $e^{- Δ E / T}$ . Repeating this visits arrangements with exactly the Boltzmann probabilities in the long run. One pass over all 256 spins is a sweep. The edges wrap around (periodic boundaries), so every spin has four neighbors.

[1]lines 56–108

# the grid is L by L spins with periodic edges
L = 16
SEED = 123456.0
# the temperature of the main measurement, well below the critical one
T_MAIN = 1.8
# sweeps discarded while the grid forgets its all-up start
EQUIL = 200
# the measured sweeps, in batches for the error estimate
BATCHES = 20
PER_BATCH = 60

#| One Metropolis sweep over the grid s at inverse temperature beta: each
#| spin in turn is flipped if that lowers the energy, and otherwise with
#| probability exp(−beta·ΔE).
def metropolis_sweep(s: Arr, beta) {
    n = len(s)
    w4 = exp(-4.0 * beta)
    w8 = exp(-8.0 * beta)
    for i to n {
        for j to n {
            ip = mod(i + 1.0, n)
            im = mod(i + n - 1.0, n)
            jp = mod(j + 1.0, n)
            jm = mod(j + n - 1.0, n)
            dE = 2.0 * s[i, j] * (s[ip, j] + s[im, j] + s[i, jp] + s[i, jm])
            r = rand()
            if dE < 0.5 {
                s[i, j] = -s[i, j]
            }
            else if dE < 4.5 {
                if r < w4 {
                    s[i, j] = -s[i, j]
                }
            }
            else if r < w8 {
                s[i, j] = -s[i, j]
            }
        }
    }
    return 0.0
}

#| |m|: the size of the average spin of the grid s.
def abs_magnetization(s: Arr) {
    n = len(s)
    m = 0.0
    for i to n {
        for j to n {
            m = m + s[i, j]
        }
    }
    return abs(m) / (n * n)
}

# the grid is L by L spins with periodic edges
L = 16
SEED = 123456.0
# the temperature of the main measurement, well below the critical one
T_MAIN = 1.8
# sweeps discarded while the grid forgets its all-up start
EQUIL = 200
# the measured sweeps, in batches for the error estimate
BATCHES = 20
PER_BATCH = 60

#| One Metropolis sweep over the grid s at inverse temperature beta: each
#| spin in turn is flipped if that lowers the energy, and otherwise with
#| probability exp(−beta·ΔE).
def metropolis_sweep(s: Arr, beta) {
    n = len(s)
    w4 = exp(-4.0 * beta)
    w8 = exp(-8.0 * beta)
    for i to n {
        for j to n {
            ip = mod(i + 1.0, n)
            im = mod(i + n - 1.0, n)
            jp = mod(j + 1.0, n)
            jm = mod(j + n - 1.0, n)
            dE = 2.0 * s[i, j] * (s[ip, j] + s[im, j] + s[i, jp] + s[i, jm])
            r = rand()
            if dE < 0.5 {
                s[i, j] = -s[i, j]
            }
            else if dE < 4.5 {
                if r < w4 {
                    s[i, j] = -s[i, j]
                }
            }
            else if r < w8 {
                s[i, j] = -s[i, j]
            }
        }
    }
    return 0.0
}

#| |m|: the size of the average spin of the grid s.
def abs_magnetization(s: Arr) {
    n = len(s)
    m = 0.0
    for i to n {
        for j to n {
            m = m + s[i, j]
        }
    }
    return abs(m) / (n * n)
}

Onsager's answer

Lars Onsager solved the two-dimensional Ising model exactly in 1944, and C. N. Yang derived its magnetization in 1952. Below the critical temperature $T_{c}$ ,

M (T) = (1 - sinh^{- 4} (2/ T))^{1/8},

and above it $M = 0$ . The formula reaches zero where $sinh (2/ T) = 1$ , which gives $T_{c} = 2/ ln (1 + 2)$ . Before running: what does the exponent $1/8$ say about the curve near $T_{c}$ ? Close to $T_{c}$ , $M$ behaves like $(T_{c} - T)^{1/8}$ , so it drops to zero almost vertically.

[2]lines 125–145

#| Yang's spontaneous magnetization of the infinite lattice at temperature T,
#| and zero at or above the critical temperature Tc.
def onsager(T) {
    sh = sinh(2.0 / T)
    if sh <= 1.0 {
        return 0.0
    }
    return pow(1.0 - pow(sh, -4.0), 0.125)
}

TC = 2.0 / log(1.0 + sqrt(2.0))
sinh_at_tc = sinh(2.0 / TC)
show TC
show sinh_at_tc

test "the critical temperature is 2 / ln(1 + sqrt 2) = 2.269185314213022, where sinh(2/T) = 1" {
    # mpmath at 50 digits (posts/refs/metropolis-ising-against-onsager.out)
    assert_near(TC, 2.269185314213022, 1.0e-14)
    assert_near(sinh_at_tc, 1.0, 1.0e-14)
    assert_near(onsager(TC), 0.0, 0.0)
}

The measurement

The grid starts with every spin up, runs 200 sweeps to forget that start, then measures $∣ m ∣$ after each of 1200 more sweeps. The absolute value is there because a finite grid can flip as a whole from mostly up to mostly down; at $T = 1.8$ that is rare, but $∣ m ∣$ makes the answer independent of it.

The 1200 values are not independent: one sweep changes only a few spins, so consecutive values are close. The standard error comes from batch means: average each block of 60 sweeps, and treat the 20 block averages as independent samples. That holds when a block is much longer than the time the grid takes to forget its state, which here is about one sweep.

[3]lines 163–212

srand(SEED)
spins = zeros(L, L) + 1.0
for k to EQUIL {
    metropolis_sweep(spins, 1.0 / T_MAIN)
}
batch = zeros(BATCHES)
for b to BATCHES {
    acc = 0.0
    for k to PER_BATCH {
        metropolis_sweep(spins, 1.0 / T_MAIN)
        acc = acc + abs_magnetization(spins)
    }
    batch[b] = acc / PER_BATCH
}
m_mc = 0.0
for b to BATCHES {
    m_mc = m_mc + batch[b] / BATCHES
}
spread = 0.0
for b to BATCHES {
    spread = spread + (batch[b] - m_mc) * (batch[b] - m_mc)
}
sigma = sqrt(spread / (BATCHES - 1.0) / BATCHES)
m_exact = onsager(T_MAIN)
show m_mc
show sigma
show m_exact

test "at T = 1.8 the sampled |m| is Onsager's 0.956857 within three standard errors" {
    # 0.95685701724741835: mpmath, 50 digits
    assert_near(m_exact, 0.9568570172474184, 1.0e-12)
    assert_near(m_mc, m_exact, 3.0 * sigma)
}

# how far the run is from Onsager's value, in |m| and in standard errors
m_gap = m_exact - m_mc
gap_in_sigma = m_gap / sigma
show m_gap
show gap_in_sigma

test "the run and Onsager's value differ by 0.00017, about an eighth of one standard error" {
    assert_near(m_gap, 0.000171, 1.0e-6)
    # an eighth is 0.125
    assert_near(gap_in_sigma, 0.131, 0.001)
}
test "the standard error is under 0.002, near the reference code's 0.00131 for 1200 sweeps" {
    # the reference Metropolis code predicts 0.00131 for 1200 sweeps; a
    # 20-batch estimate scatters by about 16 percent around it
    assert_lt(sigma, 0.002)
}

srand(SEED)
spins = zeros(L, L) + 1.0
for k to EQUIL {
    metropolis_sweep(spins, 1.0 / T_MAIN)
}
batch = zeros(BATCHES)
for b to BATCHES {
    acc = 0.0
    for k to PER_BATCH {
        metropolis_sweep(spins, 1.0 / T_MAIN)
        acc = acc + abs_magnetization(spins)
    }
    batch[b] = acc / PER_BATCH
}
m_mc = 0.0
for b to BATCHES {
    m_mc = m_mc + batch[b] / BATCHES
}
spread = 0.0
for b to BATCHES {
    spread = spread + (batch[b] - m_mc) * (batch[b] - m_mc)
}
sigma = sqrt(spread / (BATCHES - 1.0) / BATCHES)
m_exact = onsager(T_MAIN)
show m_mc
show sigma
show m_exact

test "at T = 1.8 the sampled |m| is Onsager's 0.956857 within three standard errors" {
    # 0.95685701724741835: mpmath, 50 digits
    assert_near(m_exact, 0.9568570172474184, 1.0e-12)
    assert_near(m_mc, m_exact, 3.0 * sigma)
}

# how far the run is from Onsager's value, in |m| and in standard errors
m_gap = m_exact - m_mc
gap_in_sigma = m_gap / sigma
show m_gap
show gap_in_sigma

test "the run and Onsager's value differ by 0.00017, about an eighth of one standard error" {
    assert_near(m_gap, 0.000171, 1.0e-6)
    # an eighth is 0.125
    assert_near(gap_in_sigma, 0.131, 0.001)
}
test "the standard error is under 0.002, near the reference code's 0.00131 for 1200 sweeps" {
    # the reference Metropolis code predicts 0.00131 for 1200 sweeps; a
    # 20-batch estimate scatters by about 16 percent around it
    assert_lt(sigma, 0.002)
}

The run gives $∣ m ∣ = 0.95669$ with a standard error of $0.00130$ (m_mc, sigma), and Onsager's value is $0.95686$ (m_exact): they differ by $0.00017$ , about an eighth of one standard error.

The whole curve

The next cell repeats the measurement at ten temperatures, with shorter runs (100 sweeps to settle, 200 to measure), and draws them over Onsager's curve. Before looking: the formula says $M = 0$ above $T_{c}$ . Will the 16 by 16 grid show zero there?

[4]lines 227–302

T_SWEEP = [1.6, 1.8, 2.0, 2.1, 2.2, 2.3, 2.4, 2.6, 2.8, 3.2]
SWEEP_EQUIL = 100
SWEEP_MEAS = 200
# the points on Onsager's curve
N_CURVE = 200
# the temperatures whose grids are drawn: below, near and above T_c
T_BELOW = 1.8
T_NEAR = 2.3
T_ABOVE = 3.2
grid_below = zeros(L, L)
grid_near = zeros(L, L)
grid_above = zeros(L, L)

m_T = zeros(len(T_SWEEP))
for t to len(T_SWEEP) {
    s = zeros(L, L) + 1.0
    for k to SWEEP_EQUIL {
        metropolis_sweep(s, 1.0 / T_SWEEP[t])
    }
    acc = 0.0
    for k to SWEEP_MEAS {
        metropolis_sweep(s, 1.0 / T_SWEEP[t])
        acc = acc + abs_magnetization(s)
    }
    m_T[t] = acc / SWEEP_MEAS
    if T_SWEEP[t] == T_BELOW {
        grid_below = s + 0.0
    }
    if T_SWEEP[t] == T_NEAR {
        grid_near = s + 0.0
    }
    if T_SWEEP[t] == T_ABOVE {
        grid_above = s + 0.0
    }
}
show m_T

# the same protocol run 256 times at each temperature by the reference code
# (posts/refs/metropolis-ising-against-onsager.out): the
# mean of the 200-sweep average of |m|, and its chain-to-chain spread
M_REF = [0.9797, 0.9570, 0.9108, 0.8702, 0.7980, 0.6632, 0.5233, 0.3010, 0.2131, 0.1465]
M_REF_SPREAD = [0.0015, 0.0031, 0.0098, 0.0202, 0.0402, 0.0688, 0.0747, 0.0391, 0.0226, 0.0107]

test "at every temperature the short run is within three spreads of the reference code's average" {
    assert_near(len(M_REF), len(T_SWEEP), 0.0)
    for t to len(T_SWEEP) {
        assert_near(m_T[t], M_REF[t], 3.0 * M_REF_SPREAD[t])
    }
}

T_curve = linspace(1.5, 3.3, N_CURVE)
M_curve = zeros(N_CURVE)
for k to N_CURVE {
    M_curve[k] = onsager(T_curve[k])
}
figure("Magnetization of the 2D Ising model", "temperature $T$ (units of $J/k_B$)", "magnetization |m|", 1.5, 3.3, 0.0, 1.0)
plot(T_curve, M_curve, "Onsager and Yang, infinite lattice")
style("Onsager and Yang, infinite lattice", "dash s2")
scatter(T_SWEEP, m_T, "Metropolis, 16 x 16 grid")
# each point's bar: one run's spread, the reference code's chain-to-chain standard deviation
rangebars("one run's spread, from 256 reference runs", T_SWEEP, m_T - M_REF_SPREAD, m_T + M_REF_SPREAD)
style("Metropolis, 16 x 16 grid", "s1")
style("one run's spread, from 256 reference runs", "s1")
vline("T_c", TC)
label("tc", "$T_c = 2.269$", TC + 0.18, 0.9)

grid("The 16 x 16 grid below, near and above T_c", 1, 3)
figure("The grid at T = 1.8, below T_c", "column", "row", 0.0, L, 0.0, L, "pixels")
plot2d(grid_below, "spin: up (1) or down (-1)")
colormap("diverging")
figure("The grid at T = 2.3, just above T_c", "column", "row", 0.0, L, 0.0, L, "pixels")
plot2d(grid_near, "spin: up (1) or down (-1)")
colormap("diverging")
figure("The grid at T = 3.2, well above T_c", "column", "row", 0.0, L, 0.0, L, "pixels")
plot2d(grid_above, "spin: up (1) or down (-1)")
colormap("diverging")

T_SWEEP = [1.6, 1.8, 2.0, 2.1, 2.2, 2.3, 2.4, 2.6, 2.8, 3.2]
SWEEP_EQUIL = 100
SWEEP_MEAS = 200
# the points on Onsager's curve
N_CURVE = 200
# the temperatures whose grids are drawn: below, near and above T_c
T_BELOW = 1.8
T_NEAR = 2.3
T_ABOVE = 3.2
grid_below = zeros(L, L)
grid_near = zeros(L, L)
grid_above = zeros(L, L)

m_T = zeros(len(T_SWEEP))
for t to len(T_SWEEP) {
    s = zeros(L, L) + 1.0
    for k to SWEEP_EQUIL {
        metropolis_sweep(s, 1.0 / T_SWEEP[t])
    }
    acc = 0.0
    for k to SWEEP_MEAS {
        metropolis_sweep(s, 1.0 / T_SWEEP[t])
        acc = acc + abs_magnetization(s)
    }
    m_T[t] = acc / SWEEP_MEAS
    if T_SWEEP[t] == T_BELOW {
        grid_below = s + 0.0
    }
    if T_SWEEP[t] == T_NEAR {
        grid_near = s + 0.0
    }
    if T_SWEEP[t] == T_ABOVE {
        grid_above = s + 0.0
    }
}
show m_T

# the same protocol run 256 times at each temperature by the reference code
# (posts/refs/metropolis-ising-against-onsager.out): the
# mean of the 200-sweep average of |m|, and its chain-to-chain spread
M_REF = [0.9797, 0.9570, 0.9108, 0.8702, 0.7980, 0.6632, 0.5233, 0.3010, 0.2131, 0.1465]
M_REF_SPREAD = [0.0015, 0.0031, 0.0098, 0.0202, 0.0402, 0.0688, 0.0747, 0.0391, 0.0226, 0.0107]

test "at every temperature the short run is within three spreads of the reference code's average" {
    assert_near(len(M_REF), len(T_SWEEP), 0.0)
    for t to len(T_SWEEP) {
        assert_near(m_T[t], M_REF[t], 3.0 * M_REF_SPREAD[t])
    }
}

T_curve = linspace(1.5, 3.3, N_CURVE)
M_curve = zeros(N_CURVE)
for k to N_CURVE {
    M_curve[k] = onsager(T_curve[k])
}
figure("Magnetization of the 2D Ising model", "temperature $T$ (units of $J/k_B$)", "magnetization |m|", 1.5, 3.3, 0.0, 1.0)
plot(T_curve, M_curve, "Onsager and Yang, infinite lattice")
style("Onsager and Yang, infinite lattice", "dash s2")
scatter(T_SWEEP, m_T, "Metropolis, 16 x 16 grid")
# each point's bar: one run's spread, the reference code's chain-to-chain standard deviation
rangebars("one run's spread, from 256 reference runs", T_SWEEP, m_T - M_REF_SPREAD, m_T + M_REF_SPREAD)
style("Metropolis, 16 x 16 grid", "s1")
style("one run's spread, from 256 reference runs", "s1")
vline("T_c", TC)
label("tc", "$T_c = 2.269$", TC + 0.18, 0.9)

grid("The 16 x 16 grid below, near and above T_c", 1, 3)
figure("The grid at T = 1.8, below T_c", "column", "row", 0.0, L, 0.0, L, "pixels")
plot2d(grid_below, "spin: up (1) or down (-1)")
colormap("diverging")
figure("The grid at T = 2.3, just above T_c", "column", "row", 0.0, L, 0.0, L, "pixels")
plot2d(grid_near, "spin: up (1) or down (-1)")
colormap("diverging")
figure("The grid at T = 3.2, well above T_c", "column", "row", 0.0, L, 0.0, L, "pixels")
plot2d(grid_above, "spin: up (1) or down (-1)")
colormap("diverging")

Below $T_{c}$ the points sit on the curve. Above it they do not reach zero: at $T = 3.2$ the grid still has $∣ m ∣ \approx 0.15$ (m_T), and the reference code, running the same short protocol 256 times, averages $0.1465$ there. Each point is within three of the reference's run-to-run spreads of its average; the spread is largest near $T_{c}$ , $0.07$ at $T = 2.3$ and $2.4$ , where the grid wanders slowly between ordered and disordered states. The reason is the grid's size. Above $T_{c}$ the spins form patches of aligned spins that are uncorrelated with each other, and the net spin of 256 spins in a few dozen independent patches is a random sum, whose absolute value is of order one over the square root of the number of patches, not zero. On an infinite lattice that goes to zero; at $L = 16$ it does not. The same effect rounds off the sharp drop at $T_{c}$ .

The three grids show the change: at $T = 1.8$ a few flipped spins in a sea of one color; at $T = 2.3$ more of them, in small clusters, though this run ended in one of its ordered stretches (the grid wanders between ordered and disordered states there, and another seed can show patches of every size); and at $T = 3.2$ patches of both colors with no preferred direction.

A larger grid at the critical temperature

The last cell runs a 32 by 32 grid at $T_{c}$ itself, from a random start, and draws it every two sweeps. A cell is 1 where the spin is up and 0 where it is down. At $T_{c}$ the patches have no typical size: small ones flicker and large ones drift and merge, and $∣ m ∣$ wanders instead of settling.

[5]lines 335–370

L_CRIT = 32
CRIT_SEED = 20261014.0
CRIT_FRAMES = 40
SWEEPS_PER_FRAME = 2

srand(CRIT_SEED)
crit = zeros(L_CRIT, L_CRIT)
for i to L_CRIT {
    for j to L_CRIT {
        if rand() < 0.5 {
            crit[i, j] = 1.0
        }
        else {
            crit[i, j] = -1.0
        }
    }
}
crit_m = zeros(CRIT_FRAMES)
figure("A 32 x 32 grid at T_c, every two sweeps", "column", "row", 0.0, L_CRIT, 0.0, L_CRIT, "pixels")
for f to CRIT_FRAMES {
    for k to SWEEPS_PER_FRAME {
        metropolis_sweep(crit, 1.0 / TC)
    }
    crit_m[f] = abs_magnetization(crit)
    plot2d(0.5 * (crit + 1.0), "spin: up (1) or down (0)")
}
figure("|m| of the 32 x 32 grid at T_c", "sweep", "|m|")
plot(linspace(SWEEPS_PER_FRAME, SWEEPS_PER_FRAME * CRIT_FRAMES, CRIT_FRAMES), crit_m, "|m|")

test "every cell of the critical grid is +1 or -1, and |m| is at most 1" {
    assert_near(sum(abs(crit) - 1.0), 0.0, 0.0)
    assert_lt(max(crit_m), 1.0 + 1.0e-12)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

L_CRIT = 32
CRIT_SEED = 20261014.0
CRIT_FRAMES = 40
SWEEPS_PER_FRAME = 2

srand(CRIT_SEED)
crit = zeros(L_CRIT, L_CRIT)
for i to L_CRIT {
    for j to L_CRIT {
        if rand() < 0.5 {
            crit[i, j] = 1.0
        }
        else {
            crit[i, j] = -1.0
        }
    }
}
crit_m = zeros(CRIT_FRAMES)
figure("A 32 x 32 grid at T_c, every two sweeps", "column", "row", 0.0, L_CRIT, 0.0, L_CRIT, "pixels")
for f to CRIT_FRAMES {
    for k to SWEEPS_PER_FRAME {
        metropolis_sweep(crit, 1.0 / TC)
    }
    crit_m[f] = abs_magnetization(crit)
    plot2d(0.5 * (crit + 1.0), "spin: up (1) or down (0)")
}
figure("|m| of the 32 x 32 grid at T_c", "sweep", "|m|")
plot(linspace(SWEEPS_PER_FRAME, SWEEPS_PER_FRAME * CRIT_FRAMES, CRIT_FRAMES), crit_m, "|m|")

test "every cell of the critical grid is +1 or -1, and |m| is at most 1" {
    assert_near(sum(abs(crit) - 1.0), 0.0, 0.0)
    assert_lt(max(crit_m), 1.0 + 1.0e-12)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors