Simulating the 3D Heisenberg magnet

For my final project in Computational Statistical Physics at ETH Zürich I simulated the classical Heisenberg model on a cubic lattice and measured its critical temperature and five critical exponents. The project received full marks. The numbers themselves are known to much higher precision already. The point of the project was to get them from scratch and to understand where every error bar comes from.

I have written this page for someone who knows some physics but has not met critical phenomena before, so it starts with the Ising model and works up to the measurements. The interactive simulations are small JavaScript versions I wrote for this page. The project code itself (Python with Numba, run on ETH’s Euler cluster) is not published here. The plots in the second half use the project’s actual data.

01Phase transitions

Iron stops being magnetic above 770 °C, its Curie temperature. Below that temperature the atomic magnetic moments line up and the metal has a net magnetisation \(m\). Above it they point in random directions and \(m = 0\).

What makes the transition interesting is how \(m\) goes to zero. It does not jump. As the temperature approaches the critical temperature \(T_c\) from below, it shrinks like a power law,

Close to \(T_c\) the magnet is full of fluctuating regions that are already ordered, with sizes up to a correlation length \(\xi\), and \(\xi\) diverges as

At \(T_c\) itself there are ordered patches of every size and no length scale stands out. The same thing happens at the liquid–gas critical point of a fluid, where density fluctuations become large enough to scatter light and the fluid turns milky (critical opalescence).

The exponents \(\beta\), \(\nu\) and a few others are the critical exponents. Measuring them for the 3D Heisenberg model is what this project was about.

02The Ising model

The simplest model of a magnet puts a spin \(s_i = \pm 1\) on each site of a lattice. Neighbouring spins lower the energy when they agree:

The first sum runs over nearest-neighbour pairs and \(h\) is an external field. At temperature \(T\) a configuration \(s\) appears with probability \(e^{-H(s)/T}/Z\), where \(Z\) normalises the distribution. I use units with \(k_B = J = 1\) throughout.

Wilhelm Lenz proposed the model in 1920 and gave it to his student Ernst Ising, who solved the one-dimensional chain in 1925. He found no transition and guessed, wrongly, that there would be none in higher dimensions either. Peierls showed in 1936 that the 2D model does order, and Onsager solved it exactly in 1944. On the square lattice

Even this model is far too big to sum over directly. The lattice below has \(128^2\) spins and so \(2^{16\,384}\) configurations. A Monte Carlo simulation instead produces a sequence of configurations in which each one turns up with its Boltzmann probability.

The Ising model gets used well outside magnetism, because “many two-state units that prefer to agree with their neighbours” describes a lot of systems:

• Fluids. Read each site as occupied or empty and you have a lattice gas (Lee and Yang, 1952). The liquid–gas critical point of real fluids has the same exponents as the 3D Ising model.
• Alloys. Order–disorder transitions in alloys such as β-brass (CuZn), where the two atom types arrange themselves on alternating sites.
• Neural networks. Hopfield networks (1982) and Boltzmann machines (1985) are Ising models whose couplings store patterns. This was the work behind the 2024 Nobel Prize in Physics.
• Images. Ising-type priors, in which neighbouring pixels tend to agree, were an early tool for denoising images (Geman and Geman, 1984).
• Optimisation. Problems such as max-cut can be written as finding the ground state of an Ising model, which is what quantum annealers are built to do.
• Biophysics. The helix–coil transition of polypeptides follows a one-dimensional Ising-type model (Zimm and Bragg, 1959).

03From arrows to vectors

In the classical Heisenberg model each spin is a unit vector in three dimensions instead of a sign:

Heisenberg’s 1928 quantum version explained ferromagnetism through the exchange interaction; the classical model is its limit of large spin. Ising (\(n=1\) component), XY (\(n=2\)) and Heisenberg (\(n=3\)) together are the O(n) models.

The difference that matters is the symmetry. Ising spins can only all be flipped, which is a discrete symmetry. Heisenberg spins can all be rotated together by any angle without changing the energy. A continuous symmetry allows spin waves, slow twists of the magnetisation across the sample that cost almost no energy, and in low dimensions these destroy order.

So the 2D Ising model orders but the 2D Heisenberg model never does. (The 2D XY model is a special case: it has no long-range order but still has the Berezinskii–Kosterlitz–Thouless transition, driven by vortices, which earned the 2016 Nobel Prize.) Three dimensions is the lowest dimension in which the Heisenberg model has an ordinary transition, at \(T_c \approx 1.443\) on the simple cubic lattice. Isotropic magnets such as EuO and EuS belong to the same universality class.

The widget below runs both models at the same ratio \(T/T_c\).

04Critical exponents

Near the critical point most quantities behave as power laws in the reduced temperature \(t = (T - T_c)/T_c\) or in the field \(h\). The five exponents I measured are:

They are related by scaling laws, for example \(\gamma = \nu(2-\eta)\) and \(\gamma = \beta(\delta - 1)\). The one I used as a check is hyperscaling, which in terms of the ratios I measured reads

The exponents are universal: they depend only on the dimension, on the number of spin components \(n\) and on the interactions being short-ranged, not on the lattice or the coupling strength. Kadanoff’s block-spin argument (1966) and Wilson’s renormalisation group (1971) explain why. The table shows how close the 3D classes are, which is what makes the measurement non-trivial.

The Heisenberg values are from Campostrini et al. (2002), the others from Pelissetto and Vicari (2002).

05Finite-size scaling

A simulation box of \(L \times L \times L\) spins with periodic boundaries cannot hold a correlation length much larger than \(L\), so nothing diverges. Finite-size scaling uses this. Near \(T_c\) the box size plays the role of \(\xi\), and at \(T_c\)

More generally an observable \(Q\) with exponent \(\rho\) satisfies

For dimensionless ratios such as the Binder cumulant \(U_4\) and \(\xi_L / L\), \(\rho = 0\), so their curves for different \(L\) cross at \(T_c\). These relations only hold up to corrections of order \(L^{-\omega}\), with \(\omega \approx 0.773\) for this model (Hasenbusch, 2001). In practice that meant the small lattices were not much use for the final numbers, and I needed lattices up to \(L = 256\).

06Monte Carlo and cluster algorithms

Metropolis

A Markov chain Monte Carlo simulation moves between configurations with transition probabilities \(W(a \to b)\) chosen so that the chain samples the Boltzmann distribution. The standard sufficient condition is detailed balance, \(P(a)\, W(a \to b) = P(b)\, W(b \to a)\), together with every configuration being reachable.

The Metropolis algorithm (1953) proposes a small change and accepts it with probability \(\min(1, e^{-\Delta E / T})\). For Heisenberg spins I rotated a spin by an angle \(\theta \in [-\theta_{\max}, \theta_{\max}]\) towards a random perpendicular direction \(\mathbf p\):

I tried three cone widths at \(T_c\) with \(L = 10\). A narrow cone (\(\theta_{\max} = 0.35\)) was accepted 92% of the time but barely moved the spins; a completely random new direction was accepted only 43% of the time. I used a wide cone (\(\theta_{\max} = 1.20\), 73% acceptance) as a compromise. One sweep is \(L^3\) attempted updates.

Critical slowing down

Near \(T_c\), Metropolis gets slow exactly where the measurements have to be made. The system consists of correlated regions as large as the box, and changing them one spin at a time takes many sweeps. The autocorrelation time grows as \(\tau \sim L^{z}\), and for my Metropolis runs \(z \approx 1.8\): about 300 sweeps at \(L = 16\) and about 3700 at \(L = 64\). At \(L = 256\) that would make every independent measurement far too expensive.

Wolff cluster updates

Swendsen and Wang (1987), building on Fortuin and Kasteleyn’s random-cluster representation, showed how to flip whole clusters of spins without ever rejecting a move. Wolff (1989) simplified this to growing one cluster at a time:

1. Pick a random seed spin.
2. Add each aligned neighbour to the cluster with probability \(P_{\text{add}} = 1 - e^{-2J/T}\).
3. Repeat from every newly added spin until nothing more joins.
4. Flip the cluster.

The reason this is exact is a short counting argument. Let \(n_{\text{same}}\) be the number of aligned bonds on the boundary of the cluster (all of which must have been rejected for the cluster to stop there) and \(n_{\text{opp}}\) the number of anti-aligned ones. Building the same cluster in the reverse move gives the factor \((1-P_{\text{add}})^{n_{\text{opp}}}\) instead of \((1-P_{\text{add}})^{n_{\text{same}}}\). The ratio, \(e^{-2J(n_{\text{same}} - n_{\text{opp}})/T}\), is exactly the Boltzmann factor of the energy change, so detailed balance holds with acceptance probability 1.

The widget grows one cluster bond by bond.

For Heisenberg spins Wolff’s idea needs one extra step. Choose a random unit vector \(\hat{\mathbf r}\) and use reflections through the plane perpendicular to it as the cluster move,

A neighbour joins the cluster with probability

The projections onto \(\hat{\mathbf r}\) behave like Ising spins with bond-dependent couplings, so the same argument applies.

Metropolis against Wolff

Both lattices below start from the same critical configuration, and each gets one sweep’s worth of work per frame.

In the real 3D runs the difference was large. One Wolff update flips only a fraction \(\langle \lvert C\rvert \rangle / L^3\) of the lattice, so to compare fairly I converted cluster updates to sweep-equivalents. In those units the Wolff autocorrelation time stayed between about 1 and 2 sweeps at every size I measured, while Metropolis grew as \(L^{1.8}\):

The reflection move only works because the zero-field energy is symmetric under \(\mathbf S \to -\mathbf S\). A field breaks that, so for the critical-isotherm runs in section 11 I went back to Metropolis. (There is a way around this: add a “ghost” spin coupled to every site, so that the field becomes an ordinary bond; Kent-Dobias and Sethna, 2018.)

07What gets measured

For each lattice size and temperature the simulation records:

• the magnetisation per spin \(\mathbf M = \frac{1}{N}\sum_i \mathbf S_i\), with \(N = L^3\), and its length \(m = \lVert \mathbf M\rVert\). In a finite box at zero field \(\langle \mathbf M\rangle = \mathbf 0\) because the direction wanders (you can see this in the cube above), so I work with \(m\);
• the susceptibility \(\chi = \frac{N}{T}\langle m^2\rangle\);
• the Binder cumulant \(U_4 = 1 - \dfrac{\langle m^4\rangle}{3\langle m^2\rangle^2}\), which for three-component spins goes from \(2/3\) in the ordered phase to \(4/9\) in the disordered phase;
• the correlation function \(G(r) = \frac{1}{N}\sum_i \langle \mathbf S_i \cdot \mathbf S_{i+r}\rangle\);
• the second-moment correlation length, from the structure factor \(S(\mathbf k) = \sum_{\mathbf x} e^{i\mathbf k \cdot \mathbf x}\langle \mathbf S_0 \cdot \mathbf S_{\mathbf x}\rangle\): \(\xi_L = \frac{1}{2\sin(\pi/L)}\sqrt{\frac{S(\mathbf 0)}{S(\mathbf k_{\min})} - 1}, \qquad \mathbf k_{\min} = \frac{2\pi}{L}\hat{\mathbf e}_\alpha .\)

Summing for \(S(\mathbf k)\) directly costs \(\mathcal O(L^6)\) per measurement, which is hopeless at \(L = 128\). With a fast Fourier transform, \(S(\mathbf k) = \lvert \mathrm{FFT}[\mathbf S(\mathbf x)]\rvert^2\) costs \(\mathcal O(L^3 \log L)\), and an inverse FFT of \(S(\mathbf k)\) gives \(G(r)\) for all displacements at once.

Every \((L, T)\) point was run with 8 independent seeds. At large \(L\) I started Wolff from a fully ordered state, because starting from random spins took a long time to build up the magnetisation near \(T_c\). After thermalisation I estimated autocorrelation times by blocking: the time series is grouped into bins of size \(b = 2^\ell\) and

levels off at the true \(\tau\) once \(b \gg \tau\).

08Finding T_c

Because \(U_4\) does not depend on \(L\) at \(T_c\), the curves for different lattice sizes cross there.

The five large sizes give ten pairwise crossings. Taking their mean, with a seed bootstrap for the error (section 12),

which agrees with \(T_c = 1/\beta_c = 1.4430(2)\) from Holm and Janke (1993). All the critical-point runs below use \(T = 1.443\).

09β/ν and γ/ν

At \(T_c\) the magnetisation and susceptibility are power laws in \(L\), so the slopes on a log–log plot give \(-\beta/\nu\) and \(\gamma/\nu\).

The small lattices clearly bend away from the line. I chose the lower cutoff from the fit quality: with \(L_{\min} \le 32\) the reduced \(\chi^2\) was between about 30 and 100, and with \(L_{\min} = 48\) it dropped to about 1.0 and 1.9. Fitting \(L = 48\) to \(256\) gives

both within \(1\sigma\) of the literature. The hyperscaling check gives \(2\beta/\nu + \gamma/\nu = 3.006 \pm 0.039\), which is \(0.15\sigma\) from 3.

10ν from a data collapse

The scaling form says that \(\xi_L / L\) depends on \(T\) and \(L\) only through \(x = (T - T_c)\, L^{1/\nu}\). If you plot the data against \(x\) with the right \(T_c\) and \(\nu\), all five lattice sizes should lie on one curve; with the wrong values they spread apart. The sliders below do exactly this with my data.

A collapse has no fitting function, so the usual count of degrees of freedom does not apply. I scored candidate \((T_c, \nu)\) pairs with the leave-one-out score above and minimised it on a grid. The result is

where 0.0064 of the error comes from a seed bootstrap and 0.0085 from how much \(\nu\) moved when I dropped the smallest lattice. That is \(0.6\sigma\) from the literature value \(0.7112(5)\). I tried the same collapse with \(U_4\) and got \(\nu = 0.719 \pm 0.013\), but its score was about 3.1 compared with 1.29 for \(\xi_L/L\), so I used \(\xi_L/L\) for the final number.

11η and δ

The correlation function at T_c

At \(T_c\) in three dimensions, \(G(r) \sim r^{-(1+\eta)}\). On a periodic lattice, though, two spins a distance \(r\) apart are also a distance \(L - r\) apart the other way round the box, and both paths contribute. That is why the measured curves bend upwards as \(r\) approaches \(L/2\):

I fitted all five lattice sizes together with one shared exponent and an amplitude and offset for each size:

This gives \(\eta = 0.036 \pm 0.022\), within \(0.1\sigma\) of the literature value \(0.0375(5)\). Most of the error comes from the choice of inner cutoff: \(\eta\) moved by about 0.04 as I varied \(r_{\min}\) from 4 to 8. As a cross-check, \(\chi \sim L^{2-\eta}\) on the same lattices gives \(\eta = 0.026 \pm 0.013\).

The critical isotherm

The last exponent describes the response to a field at \(T_c\), \(M_z \sim h^{1/\delta}\). A straight log–log fit on the largest lattice gives \(\delta_{\text{eff}} = 5.01\), well above the literature value. The problem is corrections to scaling in the field. The renormalisation group (derived in the appendix of my report) suggests

which I fitted jointly over \(L \ge 32\) and \(h \le 0.1\):

The corrected fit gives \(\delta = 4.79 \pm 0.09\) with \(\chi^2/\text{dof} = 0.70\), against \(4.783(3)\) in the literature.

12Error bars

Each uncertainty I quote has two parts added in quadrature,

The sample part estimates how much the answer would move with new random seeds. For the slope fits (β/ν, γ/ν, η, δ) this is the usual least-squares error, the change in the parameter that raises \(\chi^2\) by 1. \(T_c\) and \(\nu\) come from crossings and collapses rather than slopes, so for those I used a seed bootstrap: resample the eight seeds with replacement, rebuild every \((L, T)\) point from the same resampled seeds, redo the analysis, and repeat 2000 times. Using the same seeds at every temperature keeps the correlations between neighbouring temperatures intact.

The choice part covers decisions the data do not settle, mainly which lattice sizes go into a fit. I took the largest shift when the smallest lattices were dropped, or for \(\eta\) half the spread over \(r_{\min} \in \{4, 5, 6, 8\}\).

For most of the exponents the choice part was the larger of the two. Corrections to scaling, not statistical noise, were the main limitation.

13Results

All six results agree with the 3D O(3) values within one standard deviation, and the collapse rules out the neighbouring Ising and XY values of \(\nu\). Without the Wolff algorithm, lattices of \(L = 256\) would have been out of reach.

14Limitations

The literature values come from Campostrini et al. (2002), who used an improved φ⁴ Hamiltonian, tuned so that the leading \(L^{-\omega}\) corrections cancel, together with high-temperature series. Their error bars are far smaller than anything a direct simulation of the standard model can reach, so my results are best seen as a consistency check on the universality class, not a competing measurement.

The main limitations were:

• the standard model has sizeable \(L^{-\omega}\) corrections, which dominate several of my error bars;
• eight seeds per point was what the compute time allowed;
• the field runs used Metropolis, which is slow at large \(L\), so the isotherm stops at \(L = 128\);
• several observables come from the same Wolff runs, so their errors are correlated.

A more precise follow-up would use histogram reweighting around \(T_c\), a denser temperature grid, or an improved Hamiltonian.

15References

1. K. Binder, “Finite size scaling analysis of Ising model block distribution functions,” Z. Phys. B 43, 119 (1981).
2. M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi and E. Vicari, “Critical exponents and equation of state of the three-dimensional Heisenberg universality class,” Phys. Rev. B 65, 144520 (2002).
3. M. Hasenbusch, “Eliminating leading corrections to scaling in the three-dimensional O(N)-symmetric φ⁴ model: N = 3 and 4,” J. Phys. A 34, 8221 (2001).
4. C. Holm and W. Janke, “Critical exponents of the classical three-dimensional Heisenberg model: a single-cluster Monte Carlo study,” Phys. Rev. B 48, 936 (1993).
5. J. Kent-Dobias and J. P. Sethna, “Cluster representations and the Wolff algorithm in arbitrary external fields,” Phys. Rev. E 98, 063306 (2018).
6. J.-K. Kim, A. J. F. de Souza and D. P. Landau, “Numerical computation of finite size scaling functions: an alternative approach to finite size scaling,” Phys. Rev. E 54, 2291 (1996).
7. N. D. Mermin and H. Wagner, “Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic Heisenberg models,” Phys. Rev. Lett. 17, 1133 (1966).
8. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087 (1953).
9. A. Pelissetto and E. Vicari, “Critical phenomena and renormalization-group theory,” Phys. Rep. 368, 549 (2002).
10. R. H. Swendsen and J.-S. Wang, “Nonuniversal critical dynamics in Monte Carlo simulations,” Phys. Rev. Lett. 58, 86 (1987).
11. U. Wolff, “Collective Monte Carlo updating for spin systems,” Phys. Rev. Lett. 62, 361 (1989).