Ising Model Simulator

Take a grid of tiny arrows, each pointing up or down. Give neighbouring arrows a reason to agree, then shake the whole grid with heat. That is the entire model.

Wilhelm Lenz introduced it in 1920, and a survey of its first hundred years notes that he made no mention of ferromagnetism at the time. The work fell to his doctoral student in Hamburg, Ernst Ising, who graduated in 1924 and published in 1925. Ising supposed that only neighbours matter, and that two neighbours have the least energy when they point the same way. Then, to make the mathematics tractable, he placed his magnets along a line.

The line let him down. On a chain the magnetization falls to zero as the applied field is removed, whatever the temperature, so there is no permanent magnet left to explain. Ising concluded that the model would not account for ferromagnetism, even in higher dimensions. He left research to teach, later fled Nazi Germany for the United States, and by the same survey’s account never published again and learned only much later how famous his name had become.

The answer hiding in two dimensions

He was wrong about the higher dimensions, and it took twenty years to see how wrong. Rudolf Peierls gave an argument that on a square grid the model does behave like a magnet when it is cold enough. In 1941 Hans Kramers and Gregory Wannier found a symmetry relating the hot model to the cold one. That symmetry pins down where the change must happen, provided it happens at a single point.

Then Lars Onsager solved it. He announced the free energy at a meeting in 1942 and published the calculation in 1944: the partition function of the two-dimensional model, computed exactly, for zero applied field. In units where the coupling and Boltzmann’s constant are both 1, the critical temperature is 2 divided by the natural logarithm of 1 plus the square root of 2, about 2.269. His abstract states the consequences plainly. The energy is continuous through that temperature, and the specific heat becomes infinite there, growing as the logarithm of the distance from it.

A second exact result followed. Onsager announced a formula for the magnetization that survives below the critical temperature with no field at all, and C. N. Yang supplied the proof, published in 1952. Take the hyperbolic sine of 2 over the temperature, raise it to the fourth power, subtract its reciprocal from 1, and raise the result to the power 1/8. That small power makes the curve remarkably steep. At a temperature of 2 the formula gives 0.911, at 2.2 it still gives 0.785, and at 2.269 it gives zero.

Both formulas describe a grid with no edge and no end. Nobody has found the matching solution with a field switched on.

What a computer adds, and what it costs

A simulation does something a formula cannot: it produces states you can look at. The method came from a 1953 paper by Metropolis, the Rosenbluths and the Tellers, which described a modified Monte Carlo integration over configurations and ran it on the Los Alamos MANIAC. For spins the recipe is short. Pick a spin and work out how much energy flipping it would cost. Flip it always if the energy drops, and otherwise with a probability that shrinks exponentially with the cost divided by the temperature.

The rule is correct in a precise sense. It satisfies a condition called detailed balance, and that is what makes states turn up with the frequencies statistical mechanics assigns them, given enough time. The trouble is in those last three words.

Near the critical temperature the grid is full of patches of every size, and the magnetization of the whole grid wanders over a wide range. A single flip changes the total by two. One review of cluster methods places the mismatch between that wide distribution and those tiny steps at the core of critical slowing down. Successive states resemble each other for a long time, and the time grows quickly with the size of the grid. For single-spin dynamics in two dimensions the exponent that governs this growth has been measured as 2.1665, which means, roughly, that doubling the side of the grid multiplies the wait by more than four.

In 1987 Robert Swendsen and Jian-Sheng Wang presented a different kind of move, one that produced unusually small values of that exponent. Two years later Ulli Wolff published a version that builds a single cluster. Start from one spin. Add each neighbour that points the same way with probability 1 minus the exponential of minus 2 over the temperature, keep growing from every spin added, and then flip the whole cluster at once. With that particular probability no accept-or-reject step is needed at all. Wolff reported autocorrelation times of one or two steps per spin on grids up to 128 by 128 at the critical temperature, with no sign of slowing down.

The cluster is not a trick of bookkeeping. It moves the magnetization by amounts as large as the ones the physics is asking for.

Error bars that tell the truth

Measuring is harder than it looks, because neighbouring states in a run are not independent samples. Average ten thousand of them, compute the usual standard error, and the bar comes out too small. The correction involves the integrated autocorrelation time: the variance of the mean is larger by a factor of one plus twice that time. Near the critical point, with single flips, that factor can run to hundreds.

The practical cure is called binning, or blocking. Average the series in pairs, then pairs of pairs, and compute the naive error at each stage. While the blocks are shorter than the correlations the error keeps rising. Once they are longer it levels off, and that level is the honest error. Henrik Flyvbjerg and Henrik Gordon Petersen described the approach in 1989 as a way to obtain the true statistical error on an average of correlated data. A later paper on the method adds the warning that matters most: if the error never levels off, the run is shorter than its own memory, and only a much longer run will do.

Here is a check you can run in the tool above. Choose the 4 by 4 lattice and leave the temperature at 2.27. A grid that small has 65,536 states, few enough to add up one by one, so its averages are known exactly: the energy per spin is -1.56499 and the mean absolute magnetization per spin is 0.84360. Press Advance with 1,000 updates a few times. The measured values settle beside the exact ones, and the gap between them should be about the size of the error bar shown. Switch between Metropolis and Wolff and the answers agree, as they must.

Now look at what that exact magnetization says. On the infinite grid the value at 2.27, just above the critical temperature, is zero. On the 4 by 4 grid it is 0.84. Nothing is wrong. A grid that small simply cannot tell which side of the transition it is on.

What a finite grid cannot show

This is the limit of every simulation of the model, and it is a matter of principle, not of patience. With no field, turning every arrow over leaves the energy unchanged, so a finite grid visits mostly-up and mostly-down states equally often and its plain average magnetization is zero at every temperature. That is why the quantity measured is the absolute value, and why the absolute value never reaches zero on the hot side.

The specific heat tells the same story. Onsager’s own abstract notes that for strips of finite width the maximum of the specific heat grows with the logarithm of the width. Arthur Ferdinand and Michael Fisher later analysed the finite grid with joined edges and found a specific heat per spin whose leading term is proportional to the logarithm of the side. A peak that grows like a logarithm stays finite on any grid you can store. The infinity in the exact solution is approached and never reached.

The tool’s opening sweep shows both regimes. A 32 by 32 grid at a temperature of 1.5 gives an absolute magnetization of 0.9864, against 0.9865 from the exact formula. The agreement is that close because far below the critical temperature the correlations are short and the size of the grid hardly matters. At 2.3, just above, the same grid reports about 0.57 where the infinite grid has zero. Both numbers are right. They answer different questions.

So a careful simulation does not announce a critical temperature. It reports estimates with error bars for one size, compares them with the exact curves where those curves apply, and says where they do not. Recovering the behaviour of the infinite grid takes several sizes and a theory of how results drift with size, which is more than one picture can carry. The tool also stays in two dimensions. In three, the survey notes, it is not even known whether an exact solution exists.

Ising thought his model had failed. A century on, people still try their methods on it first, as Wolff did with his clusters, and the reason is that exact solution: it makes a simulation easy to check, and easy to get subtly wrong.