Order by disorder

The model

Place a spin s_i = \pm 1 on every site of a square lattice. The energy of a configuration \{s\} is

H = -J \sum_{\langle ij \rangle} s_i s_j \; - \; h \sum_i s_i ,

where \langle ij \rangle runs over nearest-neighbour pairs, each pair counted once, and h is an external field. With J > 0 the interaction is ferromagnetic: a pair of aligned neighbours contributes -J, an anti-aligned pair +J, so the system lowers its energy by lining up. The simulation sets J = 1 and h = 0, which leaves the Hamiltonian symmetric under flipping every spin at once. That symmetry is what gets spontaneously broken below the transition.

The lattice here is 96 \times 96 with periodic boundaries, so N = 9216 spins and 2N bonds, each spin with four neighbours. The two ground states (all up, all down) have energy -2J per spin.

In equilibrium at temperature T, each configuration appears with the Boltzmann weight e^{-\beta H}/Z, where \beta = 1/k_B T and

Z = \sum_{\{s\}} e^{-\beta H(\{s\})}

is a sum over all 2^N configurations. Everything thermodynamic follows from Z: the free energy is F = -k_B T \ln Z, the mean energy is \langle E \rangle = -\partial \ln Z / \partial \beta and the specific heat is the energy variance, C = (\langle E^2 \rangle - \langle E \rangle^2)/k_B T^2 [1, Ch. 1].

The two readouts under the sim are

m = \frac{1}{N} \sum_i s_i , \qquad e = \frac{H}{N} = -\frac{J}{N} \sum_{\langle ij \rangle} s_i s_j .

The magnetization m is the order parameter: zero on average in the disordered phase, nonzero in the ordered phase. The energy e is a local observable that measures the density of unsatisfied bonds; it runs from -2J (perfect order) toward 0 (independent random spins at T \to \infty). Rarely for an interacting model, its equilibrium value is known exactly.

A short history

The model was proposed by Wilhelm Lenz in 1920 [2], who handed it to his student Ernst Ising as a thesis problem. Ising solved the one-dimensional chain exactly and found no ordering at any T > 0. In the English translation of his 1925 paper, the conclusion for the chain reads: "thus we find no appearances of hysteresis and no spontaneous magnetization" [3]. He then argued, using an idealized limit in which whole layers flip together, that a three-dimensional version would not order either [3]. That extrapolation was wrong, and for about a decade the model was widely regarded as a failed theory of ferromagnetism. Stephen Brush's 1967 review is the standard account of this period and what followed [4].

Why can't the chain order? Domain walls are cheap. Breaking a perfectly ordered chain into two domains costs a fixed energy 2J but gains entropy k_B \ln N from where the wall can sit, so at any T > 0 walls proliferate and destroy long-range order [1, Ch. 2]. In two dimensions a domain of linear size \ell is bounded by a wall whose energy grows like \ell. Rudolf Peierls used this in 1936 to argue that the 2D model has a nonzero spontaneous magnetization at low enough temperature [5]: the number of closed contours of length \ell grows only exponentially in \ell, so at small T the energy cost wins. Peierls' argument had gaps; Griffiths repaired them in 1964 and made the proof rigorous [6].

Kramers and Wannier (1941) found a duality that maps the high-temperature expansion of the square-lattice model onto its low-temperature expansion [7]. If the transition is unique, it must sit at the self-dual point, \sinh(2J/k_B T_c) = 1. (Their quoted J/kT_c = 0.8814 uses a different normalization; it equals \ln(1+\sqrt 2) = 2J/k_B T_c here.) Baxter gives a clean modern version [1, Sec. 6.2].

Then in 1944 Lars Onsager computed the zero-field free energy exactly, by diagonalizing the transfer matrix [8]. Per spin,

-\beta f = \ln 2 + \frac{1}{8\pi^2} \int_0^{2\pi}\!\! \int_0^{2\pi} \ln\!\left[ \cosh^2 2\beta J - \sinh 2\beta J \,(\cos\theta_1 + \cos\theta_2) \right] d\theta_1\, d\theta_2 .

His abstract states the headline result plainly: "The energy is a continuous function of T; but the specific heat becomes infinite as -\log|T-T_c|." [8] This was the first exact demonstration that the partition function of a short-range model can produce a genuine singularity, and it overturned mean-field theory's prediction of a finite jump in C.

The spontaneous magnetization took longer. As Baxter records, "five years elapsed between Onsager's derivation of the free energy f and his announcement at a conference in Florence of the result for M_0" [1, Sec. 7.1]. The first published derivation was C. N. Yang's in 1952 [9]:

M_0(T) = \left[ 1 - \sinh^{-4}\!\left( \frac{2J}{k_B T} \right) \right]^{1/8}, \qquad T < T_c ,

and M_0 = 0 above T_c.

The phase transition

Setting \sinh(2J/k_B T_c) = 1 gives

\frac{k_B T_c}{J} = \frac{2}{\ln(1+\sqrt 2)} \approx 2.269185 .

Mean-field theory, which replaces each spin's neighbours by their average, predicts k_B T_c = 4J for the square lattice [1, Ch. 3]. The real transition is more than 40% lower because fluctuations, which mean-field theory ignores, are strong in two dimensions.

Near T_c, with reduced temperature t = (T - T_c)/T_c, thermodynamic quantities follow power laws whose exponents for the 2D Ising model are known exactly [1, Sec. 7.12; 10, Ch. 4 and Ch. 7]:

Quantity	Behaviour	2D Ising exponent
Specific heat	C \sim \lvert t \rvert^{-\alpha}	\alpha = 0 (logarithmic divergence)
Magnetization (T < T_c)	M_0 \sim (-t)^{\beta}	\beta = 1/8
Susceptibility	\chi \sim \lvert t \rvert^{-\gamma}	\gamma = 7/4
Critical isotherm	M \sim h^{1/\delta} at T_c	\delta = 15
Correlation length	\xi \sim \lvert t \rvert^{-\nu}	\nu = 1
Correlations at T_c	\langle s_0 s_r \rangle \sim r^{-(d-2+\eta)}	\eta = 1/4

You can check that they satisfy the scaling relations \alpha + 2\beta + \gamma = 2, \gamma = \beta(\delta - 1), \gamma = \nu(2-\eta) and d\nu = 2 - \alpha with d = 2 [10, Ch. 4]. The exponent \beta = 1/8 is unusually small: Yang's magnetization stays close to 1 until very near T_c and then drops almost vertically.

These numbers are universal. They do not depend on lattice details, on anisotropic couplings or on adding further short-range interactions, only on dimension, the symmetry of the order parameter and the range of the forces [1, Secs. 1.3 and 1.6]. Baxter notes that Onsager's own solution supplied the first evidence: with different couplings J and J' along the two axes the specific heat still diverges as \ln|T-T_c|, independent of J'/J [1, Sec. 6.1].

The physical origin is the correlation length \xi, the size of the typical correlated region, which diverges at T_c. At the critical point there is no characteristic length at all between the lattice spacing and the system size. You can see it in the sim: at T_c the picture is a nested hierarchy of clusters within clusters, with no typical domain size. Geometric clusters of like spins at T_c are fractals; Stella and Vanderzande showed their fractal dimension is exactly 187/96 \approx 1.948 [11].

Finite-size effects on a 96 \times 96 lattice. A finite system cannot have a true singularity, since Z is a finite sum of analytic functions. Once \xi reaches the system size L = 96, it stops growing, and every singularity is rounded over a temperature window of relative width roughly 1/L, here about 1%. Ferdinand and Fisher worked this out exactly for the specific heat of a periodic m \times n lattice: the peak is rounded on a scale \Delta T/T_c \sim 1/n and displaced from T_c by an amount of the same order; for a square torus it lies slightly above T_c [12].

By the up/down symmetry, \langle m \rangle = 0 exactly in any finite system at T > 0: given enough time, the lattice flips between the + and - states. What you see depends on time scales:

• Well below T_c, a flip requires growing a domain that spans the periodic box, which costs two system-spanning walls, with a free energy of order 2\sigma L (the interfacial tension \sigma is known exactly [8; 1, Sec. 7.10]). On 96 \times 96 at, say, T = 2.0 this almost never happens in practice, so m sits near +M_0 or -M_0.
• Close to T_c, the barrier is small and m wanders and changes sign often.

This is why simulations measure \langle |m| \rangle rather than \langle m \rangle. Even \langle |m| \rangle does not drop to zero at T_c on a finite lattice; finite-size scaling predicts \langle |m| \rangle \sim L^{-\beta/\nu} = L^{-1/8} at T_c [10, Ch. 4]. Doubling L reduces it by only 2^{-1/8} \approx 0.92, so on 96 \times 96 you will see a sizable |m| at T_c, decaying smoothly above it.

How the simulation works

We want samples from \pi(\{s\}) = e^{-\beta H}/Z without computing Z, whose 2^{9216} terms are out of reach. Markov chain Monte Carlo builds a random walk through configurations whose long-run distribution is \pi. If P(a \to b) is the probability of moving from configuration a to b, a sufficient condition is detailed balance,

\pi(a)\, P(a \to b) = \pi(b)\, P(b \to a),

together with ergodicity (every configuration reachable from every other). Then the chain converges to \pi from any starting state [13, Ch. 2]. Only ratios \pi(b)/\pi(a) = e^{-\beta \Delta E} appear, so Z cancels.

The Metropolis algorithm [14], generalized by Hastings [15], satisfies detailed balance with this rule:

1. Pick a spin i and propose flipping it.
2. Compute the energy change \Delta E.
3. Accept the flip with probability A = \min\!\left(1, e^{-\Delta E / k_B T}\right). Otherwise leave the spin as is.

Moves that lower the energy are always accepted; uphill moves are accepted with the Boltzmann factor. Those uphill steps are the thermal fluctuations. Without them the lattice would just slide down to the ground state.

For a single flip, only the four bonds touching spin i change sign, so

\Delta E = 2 J s_i \sum_{j \in \text{nn}(i)} s_j .

The neighbour sum takes only the values -4, -2, 0, 2, 4, so \Delta E \in \{-8, -4, 0, 4, 8\}\,J. Only two Boltzmann factors are ever needed, e^{-4J/k_BT} and e^{-8J/k_BT}, and they can be precomputed whenever the slider moves [13, Ch. 3]. At T_c these are about 0.17 and 0.03.

A sweep is N attempted flips, one per spin on average. It is the natural unit of simulation time, and every time quoted below is in sweeps. Metropolis is one of several local rules with the same equilibrium; the heat-bath (Glauber) rule is another [13, Ch. 3; 16].

Equilibration and autocorrelation. After a change of temperature the chain needs time to forget its starting point. Even in equilibrium, successive configurations are strongly correlated. The relevant scale is the autocorrelation time \tau, defined from the decay of \langle m(t)m(t+\Delta t)\rangle; a run of length t_{\text{run}} contains only about t_{\text{run}}/2\tau statistically independent samples [13, Ch. 3].

Critical slowing down. Near T_c, the system must reorganize correlated regions of size \xi using only single-spin moves, and \tau diverges as \tau \sim \xi^z. On a finite lattice \xi is capped by L, so \tau \sim L^z at T_c. Nightingale and Blöte put it this way: "The onset of criticality is marked by a divergence of both the correlation length ξ and the correlation time τ. While the former divergence yields singularities in static quantities, the latter manifests itself notably as critical slowing down." [17] Using single-spin-flip heat-bath dynamics and an eigenvalue method that sidesteps much of the statistical noise, they obtained z = 2.1665(12) for the 2D Ising model [17]. Local Metropolis dynamics is expected to share this exponent; Newman and Barkema discuss z \approx 2.17 for the Metropolis algorithm in 2D [13, Ch. 4]. With L = 96, L^{z} \approx 2 \times 10^4, so up to a prefactor of order one, independent samples at T_c arrive only every several thousand to tens of thousands of sweeps.

Cluster algorithms. The fix is to flip correlated regions all at once. Swendsen and Wang (1987) used the Fortuin and Kasteleyn mapping to the random-cluster model [18]: bonds between aligned neighbours are activated with probability p = 1 - e^{-2\beta J}, and each resulting cluster is flipped with probability 1/2. They showed that this dynamics falls outside the dynamic universality class of local updates and has a much smaller dynamic exponent [19]. Wolff (1989) simplified this to growing and flipping a single cluster from a random seed, which preferentially picks the large clusters that matter most near T_c [20]. He reported: "On lattices up to 128^2 no sign of critical slowing down is visible with autocorrelation times of 1-2 steps per spin for estimators of long-range quantities." [20] Both satisfy detailed balance for the same Boltzmann distribution, so the equilibrium physics is identical, only reached far faster [13, Ch. 4]. The price is that their dynamics no longer resembles physical relaxation.

Things to try in the sim

Reference values computed from Onsager's and Yang's formulas (energy and |m| per spin, J = k_B = 1) [8, 9, 21]:

T	Exact e	Exact M_0
1.5	-1.951	0.987
2.0	-1.746	0.911
2.2	-1.547	0.785
2.269 (T_c)	-\sqrt 2 \approx -1.414	0
2.5	-1.106	0
3.0	-0.817	0

The energy column uses the closed form that follows from differentiating Onsager's free energy [21, Ch. 13]:

e(T) = -J \coth 2\beta J \left[ 1 + \frac{2}{\pi}\left( 2\tanh^2 2\beta J - 1 \right) K_1(\kappa) \right], \qquad \kappa = \frac{2 \sinh 2\beta J}{\cosh^2 2\beta J},

where K_1 is the complete elliptic integral of the first kind. At T_c, \kappa = 1 and K_1 diverges, but its prefactor vanishes, leaving e(T_c) = -\sqrt 2\, J.

1. Watch domains grow. The preset starts hot, a salt-and-pepper picture with m \approx 0, and drops the temperature to T = 1.5. Domains of both signs appear; walls move to reduce their curvature, small domains vanish and the typical domain size grows as t^{1/2} [22]. Watch e fall quickly while m lags: e only measures wall length, while m needs one sign to win. Sometimes neither wins: quenches to very low T can end in two straight stripes spanning the periodic box, stable at T = 0 and extremely long-lived at small T. Spirin, Krapivsky and Redner found such stripe states roughly a quarter to a third of the time [23, 24]. A stripe is not equilibrium; its energy per spin sits above -2J by the two walls, about 2 \cdot 2J \cdot 96 / 9216 \approx 0.04J.

2. Sit at the critical point. T = 2.269, from a hot start. Wait. Clusters appear on every scale up to the box size, and m swings widely and changes sign. Over a few thousand sweeps, note how slowly the large-scale structure evolves compared with the small-scale flicker: critical slowing down in action.

3. Flip whole clusters at once. The same temperature with the Wolff algorithm. Each move grows one cluster with bond probability 1 - e^{-2\beta J} and flips it whole [20], so the large structure that barely moved under Metropolis now changes every frame. Wolff works only at h = 0.

4. Check the exact answer at T = 2. The preset starts ordered at T = 2.0; at that temperature, or at T = 3.0, let the system equilibrate for a few hundred sweeps and read off e. The instantaneous readout jitters by about 0.01 to 0.02, but its average over many sweeps should match the table closely. Away from T_c finite-size corrections are exponentially small in L/\xi; near T_c the agreement degrades.

5. Follow Yang. Below T_c, compare |m| with the M_0 column. At T = 1.5 and 2.0 the agreement should be good. At T = 2.2 you will notice |m| fluctuating more and lying slightly off the exact value, because \xi is already a sizable fraction of the box. Above T_c, |m| does not reach zero but shrinks steadily as T rises.

6. Paint a domain. Ordered at T = 1.5. Drag on the lattice to paint a blob of the other sign. Its wall costs energy in proportion to its length, so the blob shrinks and vanishes, the curvature-driven motion that also drives coarsening [22].

7. Flip the field. All spins up at T = 1.5 in a field h = -0.3. The up state is metastable: it survives until a droplet of down spins grows past the critical size and sweeps the lattice. Lifetimes of this kind depend on the field and the lattice size [33].

8. Antiferromagnet. J < 0 on the square lattice. The lattice orders as a checkerboard: m stays near zero while the staggered magnetization approaches one.

9. Frustration on triangles. J < 0 on the triangular lattice at T = 0.3. No arrangement satisfies all three bonds of a triangle, so the lattice never orders and keeps a finite entropy at T = 0, as Wannier showed [29].

10. Feel the equilibration time. From a fully ordered start, heat to T = 3.0: e settles within tens of sweeps. Do the same to T = 2.3: the approach takes far longer, and m keeps drifting after e looks settled. This is why cluster algorithms were invented.

Reading the lab's instruments

Fluctuations are measurements. The panel keeps running averages over the sweeps since you last changed a setting, after discarding the first 50 while the lattice adjusts. From them it computes two response functions using the fluctuation formulas of equilibrium statistical mechanics [1, Ch. 1; 13, Ch. 3]:

C = \frac{N}{T^2}\left(\langle e^2 \rangle - \langle e \rangle^2\right), \qquad \chi = \frac{N}{T}\left(\langle m^2 \rangle - \langle |m| \rangle^2\right).

The susceptibility uses |m| because on a finite lattice \langle m \rangle averages to zero once the system has flipped between its two ordered states. Both quantities peak near T_c and the peaks sharpen as the lattice grows: compare size 32 with size 128.

The Binder cumulant is

U = 1 - \frac{\langle m^4 \rangle}{3\langle m^2 \rangle^2},

which tends to 2/3 deep in the ordered phase and to 0 in the disordered phase, where m is Gaussian. Its value at T_c depends only weakly on the lattice size, so curves of U(T) for different sizes cross close to the critical point. Binder introduced it for exactly that purpose [32].

The field breaks the symmetry. With h \ne 0 the two ordered states are no longer equal. Below T_c, reversing h is a first-order transition: an all-up lattice in a weak downward field is metastable, and only decays once a droplet of down spins grows past a critical size. How long that takes depends strongly on the field, the temperature and the lattice size [33]. The hysteresis sweep makes this visible as a loop in m(h).

Three algorithms, one equilibrium. Metropolis and heat bath both flip single spins and satisfy detailed balance, so they sample the same distribution with different local dynamics [13, Ch. 3; 16]. Wolff flips whole clusters [20] and only applies at h = 0; with a field on, the lab uses Metropolis. Run each at T_c and watch how differently the large-scale structure moves.

Antiferromagnets and frustration. With J < 0 the square lattice orders as a checkerboard: m stays near zero while the staggered magnetization, which counts spins on alternating sublattices with opposite signs, approaches \pm 1. On the triangular lattice no checkerboard exists, the staggered magnetization is undefined (shown as n/a) and at best two of the three bonds of every triangle are satisfied, so the energy cannot fall below -|J| per spin [29].

Why it matters beyond magnets

Lattice gas and the liquid-gas critical point. Let n_i = (1 + s_i)/2 be the occupancy of site i by a particle. The Ising model in a field becomes a lattice gas with nearest-neighbour attraction, the field playing the role of the chemical potential. Lee and Yang proved that the two problems are mathematically equivalent [25]. Up and down domains become liquid and vapour, and T_c becomes a liquid-gas critical point. By universality, the critical exponents of real fluids such as carbon dioxide and xenon match those of the three-dimensional Ising model [1, Sec. 1.3].

Binary alloys. In an AB alloy such as brass, s_i = \pm 1 labels which atom sits on site i. Bragg and Williams (1934) treated the resulting order-disorder transition in a mean-field approximation that is equivalent to mean-field Ising theory [26].

Neural networks. Hopfield (1982) built an associative memory from binary units with symmetric couplings and asynchronous single-unit updates [27]. The network has an Ising-like energy function, stored patterns are its local minima and recall is relaxation downhill from a partial or corrupted input. Its update rule is the zero-temperature limit of the sim's single-spin dynamics, with couplings set by the stored patterns.

Frustration. Flip the sign to J < 0 (antiferromagnet). On the square lattice nothing deep changes: the lattice is bipartite, so flipping every other spin maps the antiferromagnet back onto the ferromagnet. On the triangular lattice that fails. Around each triangle, at most two of the three bonds can be satisfied. Bramwell and Gingras put the definition this way: "A frustrated system is one whose symmetry precludes the possibility that every pairwise interaction (“bond”) in the system can be satisfied at the same time." [28] Wannier (1950) solved the triangular Ising antiferromagnet exactly and found a ground state with macroscopic degeneracy, nonzero entropy at absolute zero and no ordering transition at any temperature [29]. His original numerical value for the zero-point entropy was corrected in a 1973 erratum to S_0 \approx 0.323\, N k_B [29].

Ising-like moments pointing along local axes on the corner-sharing tetrahedra of the pyrochlore lattice give spin ice. Harris and coworkers found in Ho_2Ti_2O_7 that ferromagnetic coupling plus local Ising anisotropy produces a frustrated, disordered ground state, the first experimental evidence for geometrical frustration in a ferromagnet [30]. Ramirez and coworkers then measured a zero-point entropy in Dy_2Ti_2O_7 consistent with Pauling's estimate for the proton disorder of water ice [31]. The local rule "two in, two out" on each tetrahedron maps exactly onto the ice rules of water ice, which is where the name comes from [28].

References

[1] R. J. Baxter, Exactly Solved Models in Statistical Mechanics (Academic Press, London, 1982). Ch. 1 (basic statistical mechanics, universality, lattice gas), Ch. 2 (one-dimensional Ising model), Ch. 3 (mean field model), Ch. 6 (duality), Ch. 7 (square lattice Ising model). Free PDF from the author's institution: https://physics.anu.edu.au/research/ftp/_files/Exactly.pdf

[2] W. Lenz, "Beitrag zum Verständnis der magnetischen Erscheinungen in festen Körpern," Physikalische Zeitschrift 21, 613 (1920). No DOI.

[3] E. Ising, "Beitrag zur Theorie des Ferromagnetismus," Zeitschrift für Physik 31, 253-258 (1925). DOI: 10.1007/BF02980577. English translation "Contribution to the Theory of Ferromagnetism," digital version by J. Ising, T. Cummings and U. Harsch (2001/02), Bibliotheca Augustana: https://tha.de/~harsch/anglica/Chronology/20thC/Ising/isi_fm00.html

[4] S. G. Brush, "History of the Lenz-Ising Model," Reviews of Modern Physics 39, 883-893 (1967). DOI: 10.1103/RevModPhys.39.883

[5] R. Peierls, "On Ising's model of ferromagnetism," Mathematical Proceedings of the Cambridge Philosophical Society 32, 477-481 (1936). DOI: 10.1017/S0305004100019174

[6] R. B. Griffiths, "Peierls Proof of Spontaneous Magnetization in a Two-Dimensional Ising Ferromagnet," Physical Review 136, A437-A439 (1964). DOI: 10.1103/PhysRev.136.A437

[7] H. A. Kramers and G. H. Wannier, "Statistics of the Two-Dimensional Ferromagnet. Part I," Physical Review 60, 252-262 (1941), DOI: 10.1103/PhysRev.60.252; "Part II," Physical Review 60, 263-276 (1941), DOI: 10.1103/PhysRev.60.263

[8] L. Onsager, "Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition," Physical Review 65, 117-149 (1944). DOI: 10.1103/PhysRev.65.117

[9] C. N. Yang, "The Spontaneous Magnetization of a Two-Dimensional Ising Model," Physical Review 85, 808-816 (1952). DOI: 10.1103/PhysRev.85.808

[10] M. Kardar, Statistical Physics of Fields (Cambridge University Press, 2007). Ch. 4 "The scaling hypothesis"; Ch. 7 "Series expansions" (includes the exact free energy of the square lattice Ising model).

[11] A. L. Stella and C. Vanderzande, "Scaling and fractal dimension of Ising clusters at the d=2 critical point," Physical Review Letters 62, 1067-1070 (1989). DOI: 10.1103/PhysRevLett.62.1067

[12] A. E. Ferdinand and M. E. Fisher, "Bounded and Inhomogeneous Ising Models. I. Specific-Heat Anomaly of a Finite Lattice," Physical Review 185, 832-846 (1969). DOI: 10.1103/PhysRev.185.832

[13] M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical Physics (Oxford University Press, 1999). Ch. 2 (principles of equilibrium thermal Monte Carlo); Ch. 3 "The Ising model and the Metropolis algorithm"; Ch. 4 "Other algorithms for the Ising model" (critical slowing down, Wolff and Swendsen-Wang algorithms).

[14] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, "Equation of State Calculations by Fast Computing Machines," Journal of Chemical Physics 21, 1087-1092 (1953). DOI: 10.1063/1.1699114

[15] W. K. Hastings, "Monte Carlo sampling methods using Markov chains and their applications," Biometrika 57, 97-109 (1970). DOI: 10.1093/biomet/57.1.97

[16] R. J. Glauber, "Time-Dependent Statistics of the Ising Model," Journal of Mathematical Physics 4, 294-307 (1963). DOI: 10.1063/1.1703954

[17] M. P. Nightingale and H. W. J. Blöte, "Dynamic Exponent of the Two-Dimensional Ising Model and Monte Carlo Computation of the Sub-Dominant Eigenvalue of the Stochastic Matrix," Physical Review Letters 76, 4548-4551 (1996). DOI: 10.1103/PhysRevLett.76.4548. arXiv: cond-mat/9601059

[18] C. M. Fortuin and P. W. Kasteleyn, "On the random-cluster model: I. Introduction and relation to other models," Physica 57, 536-564 (1972). DOI: 10.1016/0031-8914(72)90045-6

[19] R. H. Swendsen and J.-S. Wang, "Nonuniversal critical dynamics in Monte Carlo simulations," Physical Review Letters 58, 86-88 (1987). DOI: 10.1103/PhysRevLett.58.86

[20] U. Wolff, "Collective Monte Carlo Updating for Spin Systems," Physical Review Letters 62, 361-364 (1989). DOI: 10.1103/PhysRevLett.62.361

[21] R. K. Pathria and P. D. Beale, Statistical Mechanics, 3rd ed. (Academic Press / Elsevier, 2011). Ch. 13 "Phase Transitions: Exact (or Almost Exact) Results for Various Models."

[22] A. J. Bray, "Theory of phase-ordering kinetics," Advances in Physics 43, 357 (1994). DOI: 10.1080/00018739400101505. arXiv: cond-mat/9501089

[23] V. Spirin, P. L. Krapivsky and S. Redner, "Fate of zero-temperature Ising ferromagnets," Physical Review E 63, 036118 (2001). DOI: 10.1103/PhysRevE.63.036118. arXiv: cond-mat/0008019

[24] V. Spirin, P. L. Krapivsky and S. Redner, "Freezing in Ising ferromagnets," Physical Review E 65, 016119 (2002). DOI: 10.1103/PhysRevE.65.016119. arXiv: cond-mat/0105037

[25] T. D. Lee and C. N. Yang, "Statistical Theory of Equations of State and Phase Transitions. II. Lattice Gas and Ising Model," Physical Review 87, 410-419 (1952). DOI: 10.1103/PhysRev.87.410

[26] W. L. Bragg and E. J. Williams, "The effect of thermal agitation on atomic arrangement in alloys," Proceedings of the Royal Society of London A 145, 699-730 (1934). DOI: 10.1098/rspa.1934.0132

[27] J. J. Hopfield, "Neural networks and physical systems with emergent collective computational abilities," Proceedings of the National Academy of Sciences USA 79, 2554-2558 (1982). DOI: 10.1073/pnas.79.8.2554

[28] S. T. Bramwell and M. J. P. Gingras, "Spin Ice State in Frustrated Magnetic Pyrochlore Materials," Science 294, 1495-1501 (2001). DOI: 10.1126/science.1064761. arXiv: cond-mat/0201427

[29] G. H. Wannier, "Antiferromagnetism. The Triangular Ising Net," Physical Review 79, 357-364 (1950), DOI: 10.1103/PhysRev.79.357; erratum, Physical Review B 7, 5017 (1973), DOI: 10.1103/PhysRevB.7.5017

[30] M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske and K. W. Godfrey, "Geometrical Frustration in the Ferromagnetic Pyrochlore Ho_2Ti_2O_7," Physical Review Letters 79, 2554-2557 (1997). DOI: 10.1103/PhysRevLett.79.2554

[31] A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan and B. S. Shastry, "Zero-point entropy in 'spin ice'," Nature 399, 333-335 (1999). DOI: 10.1038/20619

[32] K. Binder, "Finite size scaling analysis of Ising model block distribution functions," Zeitschrift für Physik B 43, 119-140 (1981). DOI: 10.1007/BF01293604

[33] P. A. Rikvold, H. Tomita, S. Miyashita and S. W. Sides, "Metastable lifetimes in a kinetic Ising model: Dependence on field and system size," Physical Review E 49, 5080-5090 (1994). DOI: 10.1103/PhysRevE.49.5080