The Ising model
The 2D Ising model is a prototypical example of a system with a higher-order phase transition. Introduced by Wilhelm Lenz in 1920 as an idealization of ferromagnetic materials (and studied by Ernst Ising) it involves a square array s of spins, each either up or down (+1 or -1), corresponding to two orientations for magnetic moments of atoms. The magnetic energy of the system is taken to be
e[s_] := -1/2 Apply[Plus, s ListConvolve[{{0, 1, 0}, {1, 0, 1}, {0, 1, 0}}, s, 2], {0, 1}]
so that each pair of adjacent spins contributes -1 when they are parallel and +1 when they are not. The overall magnetization of the system is given by m[s_] := Apply[Plus, s, {0, 1}].
In physical ferromagnetic materials what is observed is that at high temperature, corresponding to high internal energy, there is no overall magnetization. But when the temperature goes below a critical value, spins tend to line up, and an overall magnetization spontaneously develops. In the context of the 2D Ising model this phenomenon is associated with the fact that those configurations of a large array of spins that have high total energy are overwhelmingly likely to have near zero overall magnetization, while those that have low total energy are overwhelmingly likely to have nonzero overall magnetization. For an n × n array s of spins there are a total of 2n2 possible configurations. The pictures below show the results of picking all configurations with a given energy e[s] (cyclic boundary conditions are assumed) and then working out their distribution of magnetization values m[s]. Even for small n the pictures demonstrate that for large e[s] the magnetization m[s] is likely to be close to zero, but for smaller e[s] two branches approaching +1 and -1 appear. In the limit n ∞ the distribution of magnetization values becomes sharp, and a definite discontinuous phase transition is observed.
Following the work of Lars Onsager around 1944, it turns out that an exact solution in terms of traditional mathematical functions can be found in this case. (This seems to be true only in 2D, and not in 3D or higher.) Almost all spin configurations with e[s] > -√2 (where here and below all quantities are divided by the total number of spins, so that -2 ≤ e[s] ≤ 2 and -1 ≤ m[s] ≤ + 1) yield m[s] 0. But for smaller e[s] one can show that
Abs[m[s]] (1 - Sinh[2 β]-4)1/8
where β can be deduced from
e[s] -(Coth[2 β](1 + 2 EllipticK[4 Sech[2 β]2 Tanh[2 β]2] (-1 + 2 Tanh[2 β]2)/π))
This implies that just below the critical point e0 = -√2 (which corresponds to β = Log[1 + √2]/2) Abs[m] ~ (e0 - e)1/8, where here 1/8 is a so-called critical exponent. (Another analytical result is that for e~e0 correlations between pairs of spins can be expressed in terms of Painlevé functions.)
Despite its directness, the approach above of considering sets of configurations with specific energies e[s] is not how the Ising model has usually been studied. Instead, what has normally been done is to take the array of spins to be in thermal equilibrium with a heat bath, so that, following standard statistical mechanics, each possible spin configuration occurs with probability Exp[-β e[s]], where β is inverse temperature. It nevertheless turns out that in the limit n ∞ this so-called canonical ensemble approach yields the same results for most quantities as the microcanonical approach that I have used; β simply appears as a parameter, as in the formulas above.
About actual spin systems evolving in time the Ising model itself does not make any statement. But whenever the evolution is ergodic, so that all states of a given energy are visited with equal frequency, the average behavior obtained will at least eventually correspond to the average over all states discussed above.
In Monte Carlo studies of the Ising model one normally tries to sample states with appropriate probabilities by randomly flipping spins according to a procedure that can be thought of as emulating interaction with a heat bath. But in most actual physical spin systems it seems unlikely that there will be so much continual interaction with the environment. And from my discussion of intrinsic randomness generation it should come as no surprise that even a completely deterministic rule for the evolution of spins can make the system visit possible states in an effectively random way.
Among the simplest possible types of rules all those that conserve the energy e[s] turn out to have behavior that is too simple and regular. And indeed, of the 4096 symmetric 5-neighbor rules, only identity and complement conserve e[s]. Of the 232 general 5-neighbor rules 34 conserve e[s]—but all have only very simple behavior. (Compositions of several such rules can nevertheless yield complex behavior. Note that as indicated on page 1022, 34 of the 256 elementary 1D rules conserve the analog of e[s].) Of the 262,144 9-neighbor outer totalistic rules the only ones that conserve e[s] are identity and complement. But among all 2512 9-neighbor rules, there are undoubtedly examples that show effectively random behavior. One marginally more complicated case effectively involving 13 neighbors is
IsingEvolve[list_, t_Integer] := First[Nest[IsingStep, {list, Mask[list]}, t]]
IsingStep[{a_, mask_}] := {MapThread[If[#2 2 && #3 1, 1 - #1, #1]&, {a, ListConvolve[{{0, 1, 0}, {1, 0, 1}, {0, 1, 0}}, a, 2], mask}, 2], 1 - mask}
where
Mask[list_] := Array[Mod[#1 + #2, 2]&, Dimensions[list]]
is set up so that alternating checkerboards of cells are updated on successive steps.
One can see a phase transition in this system by looking at the dependence of behavior on conserved total energy e[s]. If there are no correlations between spins, and a fraction p of them are +1, then m[s] p and e[s] -2 (1 - 2p)2. And since the evolution conserves e[s] changing the initial value of p allows one to sample different total energies. But since the evolution does not conserve m[s] the average of this after many steps can be expected to be typical of all possible states of given e[s].
The pictures at the top of the next page show the values of m[s] (densities of +1 cells) after 0, 10, 100 and 1000 steps for a 500*500 system as a function of the initial values of m[s] and e[s]. Also shown is the result expected for an infinite system at infinite time. (The slow approach to this limit can be viewed as being a consequence of smallness of finite size scaling exponents in Ising-like systems.)
The phase transition in the Ising model is associated with a lack of smoothness in the dependence of the final m value on e or the initial value p of m in limiting cases of the pictures above. The transition occurs at e = -√2, corresponding to p = (1 ± 2-1/4)/2. The pictures show typical configurations generated after 1000 steps from various initial densities, as well as slices through their evolution.
And what one sees at least roughly is that right around the phase transition there are patches of black and white of all sizes, forming an approximately nested random pattern. (See also pages 989 and 1149.)