Metropolis sampling finds Onsager's magnetization

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 10^{-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.

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

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 $\Delta E=2s_i{\sum }_{\text{neighbors}}{s}_j$, which is one of $-8,-4,0,4,8$. Accept the flip if $\Delta E\le 0$; otherwise accept it with probability $e^{-\Delta 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.

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$,

and above it $M=0$. The formula reaches zero where $\sinh (2/T)=1$, which gives $T_c=2/\ln (1+\sqrt{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.

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.

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?

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.