Chapter 2: How phase transitions occur in principlePart 2: symmetry breaking, ergodicity, and lattice gases

2.8Existence of phase transitions, continued

Part 1 ended with the one-dimensional argument: a single domain wall costs a fixed energy \(2J\) but gains entropy \(k_B\log N\), so long-range order is destroyed at any \(T>0\). What changes in two dimensions?

The Peierls argument in \(d = 2\)

In \(d=2\), a domain of flipped spins is bounded by a wall. If the wall crosses \(n\) bonds between flipped and unflipped spins, the energy cost relative to the ordered state is \(\Delta E \approx 2Jn\): it now grows with the length of the wall. What about the entropy? Starting from a point on the boundary, count the ways the wall can continue. With coordination number \(z\) (the "degree", in network-science terms), an upper bound on the number of wall configurations of length \(n\) is \(z^n\).

That is an overestimate. At each step the wall can go only in \(z-1\) directions, since it cannot retrace itself, and it also may not intersect itself, so the true number of allowed configurations is smaller still. Taking \((z-1)^n\), as in \(d=1\),

In the thermodynamic limit \(n\to\infty\):

• if \(T > T_c \equiv 2J/k_B\log(z-1)\), then \(\Delta F_n \to -\infty\). Walls proliferate, giving a disordered paramagnetic phase with \(M=0\). (Because \(\Delta S\) is overestimated, the true \(\Delta F_n\) is less negative than this.)
• if \(0\lt T\lt T_c\), \(\Delta F_n\) is minimized by \(n\to0\). The state with long-range order is stable: there is a net spontaneous magnetization \(M\ne0\) without an applied field.

So there is a phase transition at a temperature \(T_c>0\), and the lower critical dimension of the Ising model is \(d=1\). Note that \(T_c\) depends on the coordination number, and therefore on the type of lattice: it is not a universal quantity. For the square lattice (\(z=4\)) the estimate gives \(k_BT_c \approx 1.82J\), whereas the exact result \(k_BT_c = 2J/\sinh^{-1}(1) \approx 2.27J\) was obtained by Kramers and Wannier in 1941 without solving the model. The estimate is too low because the entropy was overestimated.

Since \(d\) is an integer in these results, it is reassuring that the conclusion also holds for non-integer \(d\): writing \(d = 1+\epsilon\), one finds \(T_c = O(\epsilon)\). The analysis used only the fact that domain walls exist, so the conclusions apply in general to Hamiltonians with a discrete symmetry.

A teaser: the impossibility of phase transitions

We can easily "show" that a phase transition is impossible. Time-reversal symmetry gives \(F_\Omega(H,J,T) = F_\Omega(-H,J,T)\), and the magnetization obeys

This "impossibility theorem" says the magnetization in zero field must vanish. What has gone wrong?

2.9Spontaneous symmetry breaking

The argument above is correct for a finite system. It fails in the thermodynamic limit because \(f(H)\) can develop a discontinuity in its first derivative \(\partial f/\partial H\). The symmetry \(f(H) = f(-H)\) does not imply \(M(0) = 0\) unless we also assume \(f\) is smooth at \(H=0\), so that its left and right derivatives agree. If \(f(H) = f(0) + O(H^{\sigma})\) with \(\sigma>1\), then

and \(f\) is smooth. But if \(f(H) = f(0) - M_s|H| + O(H^{\sigma})\), then

The finite-system free energy is a smooth parabola at \(H=0\); in the infinite system it has a cusp, which is exactly the slope discontinuity convexity allows (Part 1, §2.6). Note that the two limits

do not commute. The spontaneous magnetization depends on temperature: \(|M_s(0)| = 1\), \(|M_s| \lt 1\) for \(0 \lt T \lt T_c\), and \(M_s = 0\) for \(T \ge T_c\).

This set of phenomena is called spontaneous symmetry breaking (SSB). "Spontaneous" distinguishes it from magnetization induced by an external field, which itself breaks the symmetry. At \(H=0\) the Hamiltonian is invariant under \(\{S_i\}\to\{-S_i\}\), yet \(M = \lim_{N\to\infty}N^{-1}\sum_i\langle S_i\rangle \ne 0\): the state is not invariant under time reversal.

SSB seems paradoxical for another reason. With the Boltzmann distribution \(P_\Omega(\{S_i\}) = e^{-\beta\mathcal H_\Omega(\{S_i\})}/Z_\Omega\), if \(\mathcal H_\Omega\) is time-reversal invariant then so is \(P_\Omega\), an even function, while \(S_i\) is odd. Thus \(\langle S_i\rangle \equiv \operatorname{Tr}\,P_\Omega(\{S_i\})S_i = 0\). This is just a restatement of the impossibility theorem. In an infinite system \(Z\to\infty\), and we have to be more careful about the probability distribution.

Consider two configurations related by time reversal, A with magnetization \(M\) and B with \(-M\). The coupling term \(-H\sum_iS_i\) gives

Take the thermodynamic limit first. For \(H>0\), \(P_B/P_A\to0\): the system must be in state A, with \(+M\), and this still holds as \(H\to0^+\). For \(H\lt 0\), the system must be in B. (Take \(H\to0\) before \(N\to\infty\), and \(P_A/P_B = 1\).) Using the limit \(H\to0^+\) together with the thermodynamic limit is therefore equivalent to setting \(H=0\) but using a restricted ensemble in which the microstates contributing to \(-M_s\) are not included; \(H\to0^-\) excludes those contributing to \(+M_s\).

Continuous symmetry

Long-range order in the Ising model for \(d>1\) relies on its symmetry group being \(\mathbb Z_2\), \(\mathcal H_\Omega[\{S_i\}] = \mathcal H_\Omega[\{-S_i\}]\). Because this symmetry is discrete, domain walls have finite width: a wall separates \(\uparrow\) from \(\downarrow\), the only possible spin states, so its thickness is one lattice unit.

Things change when the spins obey a continuous symmetry, for example when the \(\mathbf S_i\) are vectors free to point in all directions, as in the Heisenberg ferromagnet:

This is in general a more realistic model of a ferromagnet, with no preferred direction for the spins. (In real crystals, spin–orbit coupling ties the spins to the lattice and can break the rotational symmetry, favoring crystallographic axes; in extreme cases the Ising model is the better description.) With \(\mathbf S_i = (S_i^x, S_i^y, S_i^z)\) and \(R(\theta)\) a rotation by \(\theta\) about any axis, the energy in zero field is invariant under rotating all spins simultaneously, \(\mathcal H_\Omega[\{R(\theta)\mathbf S_i\}] = \mathcal H_\Omega[\{\mathbf S_i\}]\). The model is rotationally invariant, with symmetry group \(O(3)\). The rotation acts on the spins, not on the lattice, which stays fixed in space.

Thermal fluctuations destroy long-range order more easily with a continuous symmetry, because there are "more directions to point in" (made precise in Chapter 11). To compensate for the larger entropy, ordering needs more energy, which a higher-dimensional lattice provides: spins have more neighbors. A detailed analysis shows that long-range order requires \(d>2\). These results reflect the fact that with a continuous symmetry, a domain wall can spread to the size of the system.

2.10Ergodicity breaking

What does SSB mean for the dynamics? Statistical mechanics is usually justified by the equivalence of time averages and ensemble averages. For an observable \(A\{\eta_i\}\), where the \(\eta_i(t)\) are the dynamical degrees of freedom (positions and momenta for a fluid),

where \(P_{\text{eqm}}\) is the equilibrium distribution on the phase space of the system. The equivalence rests on the ergodic hypothesis: as \(t\to\infty\), \(\{\eta_i(t)\}\) comes arbitrarily close to every configuration allowed by energy conservation. (For the Ising model, the Hamiltonian does not generate an equation of motion; we assume dynamics is supplied separately, as in our Monte Carlo demos.)

When \(H=0\) and \(T\lt T_c\), the regions of configuration space with \(M_\Omega = +M\) and \(M_\Omega = -M\) should be sampled equally. A system that starts with positive magnetization must eventually reverse, and the rate \(R\) at which this happens is expected to take the Arrhenius form

The barrier \(\Delta F\) is proportional to a surface area, so the lifetime \(\tau\) of a state with a given magnetization diverges as \(N\to\infty\). The initial condition decides whether the magnetization is positive or negative, and in the thermodynamic limit the system is effectively trapped in one region of phase space. This is ergodicity breaking. Each trapped region has a time-reversed partner, in one-to-one correspondence.

Ergodicity breaking and symmetry

As an illustration, take a system with degrees of freedom \(\{\vec C_i\}\), \(i=1,\dots,N\), which may be thought of as particle positions in a many-body system, and a Hamiltonian \(\mathcal H_\lambda\{\vec C_i\}\), where \(\lambda\) stands for all the coupling constants. Suppose it is translationally invariant: \(\mathcal H_\lambda\{\vec C_i\} = \mathcal H_\lambda\{\vec C_i+\vec a\}\) for any vector \(\vec a\). By varying \(\lambda\) we intend to drive a transition with SSB. The order parameter is

the Fourier transform of the density. For \(\lambda>\lambda_c\), \(\langle f_i(\vec k)\rangle = 0\); for \(\lambda\lt\lambda_c\), \(\langle f_i(\vec k)\rangle\ne0\), meaning the density is not simply constant in space. The particles have ordered: solidification, the spontaneous breaking of translational (and rotational) invariance. With \(\mathcal G\) the set of microstates included in the trace,

If all microstates are included, \(\mathcal G\) is invariant under \(\vec C_i\to\vec C_i+\vec a\). Shifting the variables and using translational invariance gives \(\langle f_i(\vec k)\rangle_{\mathcal G} = e^{-i\vec k\cdot\vec a}\langle f_i(\vec k)\rangle_{\mathcal G}\) for every \(\vec a\), so \(\langle f_i(\vec k)\rangle_{\mathcal G} = 0\) regardless of \(\lambda\). The order parameter can be nonzero only if \(\mathcal G\) does not include all the allowed microstates, and only if this residual set is not invariant under the symmetry group of \(\mathcal H_\lambda\). An example is a system constrained to have its center of mass at a given position. That is ergodicity breaking. It accompanies SSB, but not all ergodicity breaking comes from SSB: the liquid–gas coexistence line is a counterexample, with no symmetry broken.

The structure of phase space depends on \(\lambda\) and \(\beta\). In what manner is it broken? For \(\lambda\lt\lambda_c\), let \(\mathcal G_1\) be the largest set of microstates for which \(\langle f_i(\vec k)\rangle\ne0\), and let \(\mathcal G_2 = T(\mathcal G_1)\) for some translation \(T_a:\vec C_i\to\vec C_i+\vec a\). By definition \(\mathcal G_2\ne\mathcal G_1\), while \(\mathcal H(T\vec C_i) = \mathcal H(\vec C_i)\), so \(Z_2 = Z_1\) and

If ergodicity is broken so that an order parameter is nonzero in some set \(\mathcal G\) of microstates, the order parameter takes a (possibly different) nonzero value in every set related to \(\mathcal G\) by the spontaneously broken symmetry of the Hamiltonian. But the symmetry-related ergodic regions need not exhaust phase space. When they do not, the situation is known in disordered systems as replica symmetry breaking, believed to occur in spin glasses and, after vulcanization, in rubber, where cross-linked polymers can no longer pass through one another and the system is trapped in ergodic regions not related by translation.

Spin glasses and the overlap

An Ising spin glass is an Ising model whose exchange constants are random variables with some spread about a mean,

where each \(J_{ij}\) may be ferromagnetic (\(>0\)) or antiferromagnetic (\(\lt 0\)). Magnetic alloys, with magnetic moments randomly distributed in the material, are examples. The ground state is no longer simply all up or all down: around a loop with an odd number of negative bonds, not all bonds can be satisfied, which is frustration. If \(P(J) = P(-J)\), then \(\sum_im_i = 0\) even in the low-temperature spin-glass phase, so the net magnetization cannot distinguish the frozen spin-glass state from the paramagnet. It is not a good order parameter here.

Label the ergodic regions \(\sigma, \sigma', \dots\), and define the overlap between two of them,

with \(m_i^\sigma\) the average moment at site \(i\) in region \(\sigma\). The overlap is zero in the paramagnetic phase but can be nonzero in the spin-glass phase. In particular \(q^{\sigma\sigma}\ne0\), and if \(\sigma_1'\) is \(\sigma_1\) with all spins flipped, \(q^{\sigma_1\sigma_1} = q^{\sigma_1\sigma_1'}\), so the overlap does not favor a region over its time-reversed partner. If there is only one ergodic region (plus its partner), \(q\) has only one value: the replica-symmetric case. With replica symmetry breaking, \(q^{\sigma\sigma'}\) takes a range of values. Since \(q\) usually cannot be computed directly, consider its distribution

\(P(q)\) varies from sample to sample, as the \(J_{ij}\) do; what is important and potentially calculable is its average over the distribution of couplings, \([P(q)]\). It has three generic forms: a delta function at the origin when there is no order (paramagnet); \(\delta(q-\hat q)\) with a unique \(\hat q\ne0\) when there is only one way to order (ferromagnet); and a nontrivial function, not a single delta function, when ordering can occur in many ways (the putative spin-glass phase). The nature of the low-temperature phase of real spin glasses was still debated when the book was written. Calculating \([P(q)]\) as the couplings vary offers a way to map the phase diagram of a complex system, where the elementary order-parameter concept is not enough.

Quenched disorder and the replica trick

The SK model above is the infinite-range spin glass, with \(i,j\) running over all pairs and Gaussian couplings, \(P(J) = (2\pi(\delta J)^2)^{-1/2}\exp[-(J-J_0)^2/2(\delta J)^2]\); for it mean-field theory becomes exact. The couplings are quenched disorder: random, but fixed for a given system, which is described by \(\mathcal H(\{S_i\};X)\) with \(X\) the random variables. If \(X\) could change during the dynamics, the disorder would be annealed; quenched impurities are bound to lattice sites on the timescale of the experiment. The two free energies are

A large sample can be regarded as many macroscopic subsystems, each a realization of \(X\), so a measurement averages over realizations. Systems for which this holds are self-averaging, and the observed free energy is

Averaging a logarithm is hard. The replica trick uses \(\log Z = \lim_{n\to0}(Z^n-1)/n\) (by l'Hôpital's rule), so we must compute

the trace over \(n\) replicas of the system, each with the same realization of the disorder, with \(nN(\Omega)\) degrees of freedom \(\{S_i^\alpha\}\). The effective replica Hamiltonian \(\mathcal H_n\), defined by \(e^{-\beta\mathcal H_n(\{S^\alpha\})} = [e^{-\beta\sum_\alpha\mathcal H(\{S^\alpha\};X)}]\), no longer depends on \(X\). It is symmetric under permutation of the replicas, but it is not a sum of \(n\) single-replica Hamiltonians: the disorder average has coupled the replicas. This permutation symmetry can be spontaneously broken in the limit \(n\to0\).

The central result of replica theory is that the disorder-averaged overlap distribution can be computed from replicas. For the infinite-range spin glass,

with the average taken with respect to \(\mathcal H_n\). Since \(\mathcal H_n\) is symmetric under permuting replicas, one might expect \(q^{\alpha\beta} = \hat q\) for all \(\alpha\ne\beta\), giving \([P(q)] = \delta(q-\hat q)\): simple SSB with ordering at a single value. Below the zero-field spin-glass transition this is not correct. Replica permutation symmetry is spontaneously broken, \([P(q)]\) is not a single delta function, and ergodicity is broken with more than one family of ergodic regions not related by symmetry. Interpreting replica symmetry breaking as ergodicity breaking shows that the replica method, far from being a mere averaging trick, is a complex encoding of phase space. (Giving meaning to an \(n\times n\) matrix order parameter as \(n\to0\) is delicate; see Parisi's papers.) The method has proved valuable in disordered systems, neural-network models, graph theory, and optimization.

Two Nobel Prizes in Physics trace back to this circle of ideas: Giorgio Parisi in 2021, for the interplay of disorder and fluctuations in physical systems, and John Hopfield in 2024, whose associative-memory network is a spin-glass model of neurons (Mézard, Parisi, and Virasoro, Spin Glass Theory and Beyond, ch. XIII).

2.11Fluids

Is there a connection between fluids and magnets? Consider a continuum fluid of \(N\) particles in a region \(\Omega\) confined by hard walls, with coordinates \(\vec r_i\) and momenta \(\vec p_i\), and Hamiltonian

Here \(U_1\) is an external field acting on each particle (gravity, say; usually zero), and \(U_2, U_3, \dots\) are two-body, three-body, … potentials, translationally invariant so that they depend only on coordinate differences. In the grand canonical ensemble, where \(N\) varies, with chemical potential \(\mu\),

Planck's constant \(h\) sets the unit of phase space and provides continuity with quantum statistical mechanics. The momentum integrals are Gaussian, each giving \(1/\Lambda_T^d\) with the thermal wavelength \(\Lambda_T \equiv h/\sqrt{2\pi mk_BT}\), so

with fugacity \(e^{\beta\mu}\) and configurational sum \(Q_N\). The grand potential is \(F_\Omega(T,\mu,U_1,U_2,\dots) = -k_BT\log\Xi_\Omega = U-TS-\mu N\). It is nonsingular for finite \(\Omega\) provided the potentials have a hard core, so that at most \(N_{cp}\) particles fit in \(V(\Omega)\). Its bulk density

involves the entropy per unit volume \(\sigma\) and the number density \(\rho = \lim\langle N\rangle/V\). Since \(f_b = -p\), analytic properties follow as in the magnetic case, many from convexity: (a) \(p(T,\mu)\ge0\); (b) \(p\) is continuous; (c) \(\partial p/\partial T\) and \(\partial p/\partial\mu\) exist almost everywhere; (d) \(\partial p/\partial\mu = \rho\ge0\); (e) \(\partial p/\partial T = \sigma\ge0\) (true in quantum statistical mechanics, but false for classical systems); (f) \(\partial p/\partial T\) is monotonically non-decreasing, so the heat capacity \(C_\mu\ge0\); (g) \(\partial p/\partial\mu\) is monotonically non-decreasing, so the isothermal compressibility \(K_T\ge0\).

2.12Lattice gases

A simple model of a fluid is the lattice gas of Lee and Yang. The basic idea is to relate the local density of particles in a fluid to the local density of up spins in a magnet. The microscopic density is \(\rho(\vec r) \equiv \sum_{i=1}^N\delta(\vec r-\vec r_i)\), with \(\int_V\rho\,d^dr = N(V)\), and the potentials can be rewritten in terms of it:

Now take a \(d\)-dimensional lattice with coordination number \(z\), where each site holds at most one molecule. The occupation number \(n_i = 0\) or \(1\) plays the role of \(\rho(\vec r)\), with \(N = \sum_in_i\), and a suitable Hamiltonian in the grand canonical ensemble is

with the \(\tfrac12\) cancelling the double counting. There is no kinetic term; as we saw, it only contributes to the fugacity, and we really model \(Q_N\). Define \(n_i = \tfrac12(1+S_i)\) with \(S_i = \pm1\). For short-ranged forces and low density, drop three-body and higher terms and take \(U_2(i,j) = U_2\) for nearest neighbors, zero otherwise. With the sum over ordered pairs, \(\sum_{ij}U_2(i,j) = U_2Nz\), \(\sum_{ij}U_2(i,j)S_i = zU_2\sum_iS_i\), and \(\sum_{ij}U_2(i,j)S_iS_j = 2U_2\sum_{\langle ij\rangle}S_iS_j\), since each nearest-neighbor pair appears twice. Hence

and, setting \(U_1=0\),

This is the desired result: the thermodynamics of the lattice gas follows from that of the Ising model. An attractive potential (\(U_2\lt 0\)) is a ferromagnet, and the field vanishes at \(\mu = zU_2/2\): the liquid–gas coexistence line is the Ising first-order line \(H=0\), \(T\lt T_c\). Even the ideal gas law and its corrections can be derived from the lattice gas.

From the continuum to the lattice

The configurational sum \(Q_N\) can be approximated by dividing \(\Omega\) into \(N(\Omega) = V(\Omega)/a^d\) cells of linear size \(a\), small enough (\(a \lesssim\) hard-core radius) that more than one molecule per cell is negligible. Then \(\int_\Omega\prod_i d^dr_i \approx a^{dN}\sum_\alpha\), with the number of cells generally different from the number of particles. Replacing \(U_2(\vec r_i,\vec r_j)\) by \(U_2(\alpha,\beta)\) when \(\vec r_i\in\alpha\), \(\vec r_j\in\beta\), and keeping short-ranged forces (nonzero only for neighboring cells, with constant value \(U_2\)), \(U\{\vec r_i\}\to U_2\sum_{\langle\alpha\beta\rangle}n_\alpha n_\beta+\cdots\). Each configuration \(\{n_\alpha\}\) corresponds to \(N!\) configurations \(\{\vec r_i\}\), so \(Q_N = N!\,a^{dN}\sum'_{\{n_\alpha\}}e^{-\beta U\{n_\alpha\}}\), where the prime constrains \(\sum_\alpha n_\alpha = N\). That constraint makes \(Q_N\) hard, but the grand sum removes it: summing over all \(N\) of constrained sums gives the unrestricted sum over all \(2^{N(\Omega)}\) states of \(\{n_\alpha\}\), the advantage of the grand canonical ensemble. Thus

So, heuristically, the grand partition function of the lattice gas follows from that of the continuum fluid.

2.13Equivalence in statistical mechanics

The lattice-gas mapping illustrates an equivalence between two systems. The term is used in two ways:

1. Exact equivalence: an exact mapping between the partition functions of two models, such as Ising model \(\leftrightarrow\) lattice gas.
2. Approximate equivalence: partition functions that are not equal or exactly related, but behave identically near a critical point. Such models are in the same universality class. An example we will meet is Landau theory, whose effective Hamiltonian is roughly a power series in \(M\) to order \(M^4\), and the Ising model, whose free energy written in these terms involves \(\log\cosh M\).

The phase diagram, thermodynamics, and correlation functions depend on (a) the specific coupling constants, (b) the symmetry of the model (discrete, continuous, the actual symmetry group), (c) the lattice type, (d) the dimensionality, and (e) everything else you can think of. The critical behavior turns out to depend only on (a) the symmetry of the order parameter, (b) the dimensionality, and (c) the nature of the critical point, and not on coupling constants, lattice type, or the precise Hamiltonian. (A notable exception is Baxter's eight-vertex model, where an exponent varies continuously with the couplings; this is understood within the RG and occurs only for models with marginal operators.) Models related by exact equivalence obviously share critical behavior. That unrelated systems so often do too, in experiment and theory, is the phenomenon of universality.

2.14Miscellaneous remarks

History of the thermodynamic limit

1. 1936Hendrik Kramers realizes that the limit \(\lim_{N\to\infty}N^{-1}\log Z\) can give distinct functions not analytically joined: phase transitions require the thermodynamic limit.
2. 1937Most physicists are unaware of this; even Sommerfeld believes a partition function describes only one phase. At a congress in Amsterdam for the centenary of van der Waals's birth, Kramers, as chairman, puts the question of whether a partition function can give a sharp transition to a vote. The outcome is "inconclusive."
3. 1941Kramers and Wannier use the transfer-matrix method and duality to find \(T_c\) of the \(d=2\) Ising model, assuming a single transition in the thermodynamic limit.
4. 1944Lars Onsager solves the \(d=2\) Ising model exactly. This settles the relevance of the thermodynamic limit, and shows that critical behavior can disagree with mean-field theory: the exponents need not be those of the Weiss ferromagnet or the van der Waals gas.

Do quantum effects matter?

Condensed matter is held together by quantum mechanics, but its gross characteristics need not be determined by quantum considerations. The self-avoiding walk of Chapter 1 is an example: quantum mechanics explains the polymer's existence and bond energies, but its long-scale statistics are independent of microscopic details. For continuous transitions at \(T>0\) the answer is no: as \(\xi\to\infty\), thermodynamic behavior is determined by regions much larger than the length scales at which quantum effects matter. Intrinsically quantum systems such as ferromagnets and superconductors can therefore be modeled by classical ones such as the Ising model. (For first-order transitions a quantum effect is not ruled out, though Goldenfeld knew of no example.) At \(T=0\) quantum effects are important, but a quantum system at \(T=0\) in \(d\) dimensions is related to a classical system in \(d+1\) dimensions, most cleanly seen with path integrals, so classical statistical mechanics can still be used.

This completes Chapter 2. Chapter 3 puts these ideas to work with exact and approximate calculations.

Further reading

• N. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group (Addison-Wesley, 1992), Chapter 2.
• R. Peierls, "On Ising's model of ferromagnetism," Proc. Cambridge Philos. Soc. 32, 477 (1936).
• H. A. Kramers and G. H. Wannier, "Statistics of the two-dimensional ferromagnet," Phys. Rev. 60, 252 and 263 (1941).
• L. Onsager, "Crystal statistics. I," Phys. Rev. 65, 117 (1944).
• T. D. Lee and C. N. Yang, "Statistical theory of equations of state and phase transitions. II. Lattice gas and Ising model," Phys. Rev. 87, 410 (1952).
• D. Sherrington and S. Kirkpatrick, "Solvable model of a spin-glass," Phys. Rev. Lett. 35, 1792 (1975).
• M. Mézard, G. Parisi, and M. A. Virasoro, Spin Glass Theory and Beyond (World Scientific, 1987).
• M. Dresden, "Kramers's contributions to statistical mechanics," Physics Today 41(9), 26 (1988).