The Ising Model

In an ideal paramagnet, each dipole responds only to an external field. In real materials, atomic dipoles are influenced by their neighbors: the energy depends on whether neighboring dipoles are parallel or antiparallel. When they prefer to align parallel, even with no external field, the material is a ferromagnet (iron is the most familiar example); when they prefer antiparallel, it is an antiferromagnet (Cr, NiO, FeO). The spontaneous, long-range order shows up as a net nonzero magnetization, which a paramagnet can’t have.

Raising the temperature causes random fluctuations that reduce the magnetization. For every ferromagnet there is a critical temperature, the Curie temperature (about 1043 K for iron), above which the net magnetization vanishes and the material becomes a paramagnet. Even below it, a piece of iron needn’t be magnetized: it breaks up into domains, each with billions of aligned dipoles but pointing in different directions. Heating it, applying a field, and cooling it again aligns the domains into a permanent magnet.

The Ising model describes a single domain. It keeps only the tendency of neighbors to align and assumes a preferred axis, so each dipole points either parallel or antiparallel to it. It is not an accurate model of a real ferromagnet at low temperature, where the quantum mechanics is subtler and the excitations are long-wavelength magnons, but it turns out to be much more accurate near the Curie temperature.

Let \(s_i = +1\) when dipole \(i\) points up and \(s_i = -1\) when it points down. A pair of neighboring dipoles contributes \(-\varepsilon\) when parallel and \(+\varepsilon\) when antiparallel, which is \(-\varepsilon s_i s_j\) in either case. The total energy, summed over neighboring pairs, and the partition function, summed over all sets of alignments, are

With \(N\) dipoles, that sum has \(2^N\) terms, so brute force is hopeless for any real sample. For a 4 × 4 lattice, though, there are only \(2^{16} = 65{,}536\) states, few enough for the computer to add up exactly. Click dipoles to flip them; the dots mark parallel and antiparallel bonds (Problem 8.15).

Problem 8.17. With just two dipoles, the states \(\uparrow\uparrow\) and \(\downarrow\downarrow\) have energy \(-\varepsilon\), and \(\uparrow\downarrow\) and \(\downarrow\uparrow\) have \(+\varepsilon\), so

Both dipoles pointing up is more likely than one up and one down when \(\tfrac12 P_{\text{parallel}} > P_{\text{antiparallel}}\), which works out to \(e^{2\beta\varepsilon} > 2\), that is \(T < 2\varepsilon/(k_{\mathrm{B}}\ln 2) \approx 2.885\,\varepsilon/k_{\mathrm{B}}\).

A real ferromagnet is three-dimensional and not exactly solvable, but the one-dimensional chain is. With \(U = -\varepsilon(s_1s_2 + s_2s_3 + \cdots + s_{N-1}s_N)\), the last sum gives \(\sum_{s_N = \pm1} e^{\beta\varepsilon s_{N-1}s_N} = 2\cosh\beta\varepsilon\) whatever \(s_{N-1}\) is, and so on down the chain:

These are exactly the results for a two-state paramagnet, with \(\mu B\) replaced by \(\varepsilon\), except that here the dipoles line up with each other rather than with a field. As \(T \to 0\), \(\bar{U} \to -N\varepsilon\) (perfectly aligned); as \(T \to \infty\), \(\bar{U} \to 0\). The order sets in gradually, with no abrupt transition at any nonzero temperature: each dipole has only two neighbors, too few to sustain long-range order. The chain below runs the Metropolis algorithm described further down; the domains grow as it cools, but never take over completely.

A very crude approximation can “solve” the Ising model in any dimension and shows why dimensionality matters. Focus on one dipole with \(n\) nearest neighbors (2 in one dimension, 4 on a square lattice, 6, 8, or 12 for simple, body-centered, or face-centered cubic). If \(\bar{s}\) is the average alignment of its neighbors, \(E_\uparrow = -\varepsilon n\bar{s}\) and \(E_\downarrow = +\varepsilon n\bar{s}\), so

Up to here this is exact. The approximation is to assume there are no fluctuations, so every neighborhood is typical and \(\bar{s}_i = \bar{s}\), the same idea used for the van der Waals equation:

Solve it graphically. When \(\beta\varepsilon n < 1\), the slope of the tanh at the origin is less than 1, and the only solution is \(\bar{s} = 0\), which is stable. When \(\beta\varepsilon n > 1\), \(\bar{s} = 0\) becomes unstable and two stable solutions appear, equally likely to be positive or negative. The system must choose one: the symmetry is spontaneously broken. The critical temperature is

proportional to \(\varepsilon\) and to \(n\): more neighbors take more thermal energy to disorder. In one dimension it predicts \(T_c = 2\varepsilon/k_{\mathrm{B}}\), which we know is wrong, but the approximation improves with dimension (it becomes exact above four). Problem 8.22 adds a field, giving each dipole \(\mp\mu_B B\) and \(\bar{s} = \tanh[\beta(\varepsilon n\bar{s} + \mu_B B)]\).

Problem 8.24. Using \(\tanh x \approx x - x^3/3\), the magnetization just below \(T_c\) is \(\bar{s} \approx \sqrt{3/T_c}\,(T_c - T)^{1/2}\), a critical exponent \(\beta = 1/2\). The susceptibility \(\chi = (\partial M/\partial B)_T\) diverges as \(\mu_B/k_{\mathrm{B}}(T - T_c)\) above \(T_c\) and half that below, so \(\gamma = 1\) on both sides. The exact exponents differ: \(\beta = 1/8\) and \(\gamma = 7/4\) in two dimensions (from Onsager’s solution), and \(\beta \approx 1/3\), \(\gamma \approx 1.24\) in three. Here is the mean-field magnetization for a square lattice (\(n = 4\), \(T_c = 4\varepsilon/k_{\mathrm{B}}\)) against Onsager’s exact result (\(T_c \approx 2.27\,\varepsilon/k_{\mathrm{B}}\)).

A 10 × 10 lattice already has \(2^{100} \approx 10^{30}\) states, impossible to enumerate. Sampling a million states at random wouldn’t work well either, since most would be irrelevant. The better idea is importance sampling, using the Boltzmann factors themselves as a guide. The Metropolis algorithm (1953): pick a dipole at random and compute the energy change \(\Delta U\) its flip would cause. If \(\Delta U \le 0\), flip it; if \(\Delta U > 0\), flip it with probability \(e^{-\Delta U/k_{\mathrm{B}}T}\). Repeat.

For two states 1 and 2 that differ by one flip, with \(U_1 \le U_2\), the transition probabilities are \(P(2\to1) = 1/N\) and \(P(1\to2) = (1/N)\,e^{-\beta(U_2-U_1)}\), so their ratio is exactly the ratio of Boltzmann factors. This detailed balance carries over to any sequence of steps, so states occur with the frequencies Boltzmann statistics demands. On a square lattice \(\Delta U = 2\varepsilon s_i(\text{sum of the four neighbors})\) takes only a few values, so the needed exponentials are computed in advance. The edges wrap around (periodic boundary conditions).

The lecture’s Python code scans the temperature, equilibrates, and then averages the energy, magnetization, heat capacity, and susceptibility. This does the same in your browser on a smaller lattice, using fluctuations: \(C = (\langle U^2\rangle - \langle U\rangle^2)/k_{\mathrm{B}}T^2\) and \(\chi \propto (\langle M^2\rangle - \langle |M|\rangle^2)/k_{\mathrm{B}}T\).

Problem 8.29 quantifies clustering with the correlation function: for dipoles a distance \(r\) apart,

where \(\xi\) is the correlation length. At \(T = T_c\), \(\xi \to \infty\): there is no characteristic scale, and near \(T_c\) it follows a power law, \(\xi \sim |T - T_c|^{-\nu}\). This is computed from the live simulation above, averaged over all pairs along rows and columns.

Problem 8.32. Divide the lattice into 3 × 3 blocks and replace each by a single dipole that follows the majority rule. The new lattice is a third as wide, and it looks like an Ising lattice at a different temperature. Away from \(T_c\), repeating the transformation drives the system toward one of the trivial fixed points, zero or infinite temperature. At \(T_c\), the system is invariant under the transformation, because of its scale-free nature: the critical point is itself a fixed point. Averaging away the small-scale details leaves the critical behavior unchanged, which is why very different systems, such as magnets and fluids near their critical points, share the same critical exponents. This is one version of the renormalization group, for which Kenneth Wilson received the 1982 Nobel Prize in Physics.