Spin ice

Water ice and Pauling's entropy

Look inside an ice cube and every oxygen sits at the centre of four hydrogen bonds. Each bond carries one proton, close to one end or the other. Bernal and Fowler argued in 1933 that each oxygen keeps exactly two protons close, so that the solid is built from intact water molecules [1]. Slater later used the same constraint for hydrogen-bonded crystals and called it the ice rule: two near, two far at every vertex [2; 3, Sec. 8.1].

Here's the surprise: the ice rule still leaves an enormous number of arrangements. Pauling counted them in 1935, and his argument still repays a slow reading [4]. Take N molecules and 2N bonds. Ignoring the rule, each bond has two choices, 2^{2N} in all. At each oxygen only 6 of the 16 arrangements of its four bonds obey the rule. Treating the oxygens as independent filters,

W \approx 2^{2N} \left(\frac{6}{16}\right)^{N} = \left(\frac{3}{2}\right)^{N}, \qquad S_0 = R \ln \tfrac{3}{2} \ \text{per mole of water}.

Giauque and Stout's calorimetry on ice found a residual entropy close to Pauling's number [5]. Nagle refined the count with a series expansion and found S_0 = 0.8145 \pm 0.0002 cal K^{-1} mol^{-1} for both hexagonal and cubic ice, slightly above Pauling's 0.806 [6].

From water to spin ice

Now swap the protons for magnetic moments. In the pyrochlore oxides Ho_2Ti_2O_7 and Dy_2Ti_2O_7 the rare-earth ions sit on the corners of a network of corner-sharing tetrahedra, and a strong crystal field forces each moment to point along the line joining the centres of its two tetrahedra, into one or out of the other. These local \langle 111 \rangle axes make each moment an Ising variable \sigma_i = \pm 1. Harris and collaborators found that Ho_2Ti_2O_7 has ferromagnetic interactions and still doesn't order. On these axes a ferromagnetic coupling is frustrated, and the best it can manage is two moments pointing in and two out of every tetrahedron, which is exactly the ice rule [7]. They called the state spin ice. Anderson had already noticed in 1956 that Ising spins on this lattice, the B sites of a spinel, map onto Pauling's problem [8].

Ramirez and collaborators then measured the specific heat of Dy_2Ti_2O_7 and integrated C/T. The entropy released on warming fell short of R \ln 2 by about Pauling's \tfrac{R}{2} \ln \tfrac{3}{2} = 1.68 J mol^{-1} K^{-1} per mole of Dy [9]. The factor \tfrac12 is just bookkeeping. There are two spins per tetrahedron, as there are two protons per oxygen. A magnet in a lab had reproduced the entropy of frozen water. Bramwell and Gingras reviewed the field in 2001 [10], and den Hertog and Gingras showed that the long-range dipolar interaction between the moments is what produces the ice rules in these materials [11].

The model in this lab

The lab simulates nearest-neighbour spin ice. On a vertex t, a tetrahedron in the pyrochlore or a site of the square lattice, define the charge

Q_t = \frac{n_{\text{in}} - n_{\text{out}}}{2} \in \{-2, -1, 0, 1, 2\}.

With \sigma_i = +1 meaning out of one end of bond i and into the other, the Hamiltonian is

H = J \sum_{\langle ij \rangle} \sigma_i \sigma_j - h \sum_i \sigma_i \, \hat e_i \cdot \hat h = 2J \sum_t Q_t^2 - 2J N_t - h \sum_i \sigma_i \, \hat e_i \cdot \hat h ,

where the pair sum runs over the six pairs on each tetrahedron and \hat e_i is the local axis. The second form is the one to look at: with J > 0 the energy simply counts violations of the ice rule. Every ice state has E = -J per spin, a vertex with Q = \pm 1 costs 2J and one with Q = \pm 2 costs 8J. In the materials J is an effective coupling that absorbs the nearest-neighbour part of the dipolar interaction [10, 11].

The square lattice version puts arrows on the bonds of a square grid. With all six ice vertices equally weighted it is Lieb's square ice, the six-vertex model, whose residual entropy he computed exactly in 1967 [12; 3, Ch. 8]:

W = \left(\frac{4}{3}\right)^{3/2} = 1.5396\ldots \ \text{per vertex}, \qquad s_0 = \tfrac{3}{4} \ln \tfrac{4}{3} = 0.2158 \ \text{per spin},

against Pauling's s_0 = \tfrac12 \ln \tfrac32 = 0.2027. People have built real arrays of nanomagnets in this geometry [13], but there the six ice vertices aren't degenerate, so the lab's square ice is the idealised model.

Monopoles and the Coulomb phase

Flip one spin in an ice state and two neighbouring vertices pick up Q = +1 and Q = -1. Flip a second spin that continues the path, and the charge moves on at no extra cost. You're left with two defects joined by a string of flipped spins, each free to wander. Castelnovo, Moessner and Sondhi showed that in spin ice these are magnetic monopoles, sources and sinks of the magnetisation field, and that the dipolar interaction makes them attract and repel by a magnetic Coulomb law [14]. The string plays the part of a Dirac string, and here you can actually see it, because it's made of real flipped spins; Morris and collaborators saw such strings in Dy_2Ti_2O_7 by neutron scattering [15]. In the nearest-neighbour model the monopoles carry no energetic interaction at all, only an entropic one from the ice around them [16, 17].

Why? Because of a conservation law. Coarse-grain the arrows into a field \mathbf{B}. The ice rule says no flux is created or destroyed at a vertex, \nabla \cdot \mathbf{B} = 0, and the number of ice states with a given coarse field falls like \exp(-\tfrac{K}{2} \int \mathbf{B}^2). A divergence-free field with Gaussian weight is a magnetostatic field: its correlations decay as r^{-3} with the angular form of a dipole, and in reciprocal space

\langle B_\alpha(\mathbf q) B_\beta(-\mathbf q) \rangle \propto \delta_{\alpha\beta} - \frac{q_\alpha q_\beta}{q^2}.

Henley called this state the Coulomb phase [16]. The transverse projector is singular at \mathbf q = 0 and at every reciprocal lattice vector that plays its role. Come in along one direction and the intensity is finite; come in along the perpendicular one and it vanishes. These bow ties are pinch points. Isakov and collaborators derived the dipolar correlations for spin ice [18], and Fennell and collaborators resolved pinch points in Ho_2Ti_2O_7 with polarised neutrons, the first direct proof of a Coulomb phase in a magnet [19]. Monopoles are charges in this field. At density n they screen it over a length that shrinks as n grows, and the pinch points blur over the inverse of that length.

How the simulation works

A single-spin Metropolis step proposes flipping spin i, which joins vertices a and b. Only those two charges change, so the cost is local:

\Delta E = 4J \left[1 + \sigma_i (Q_a - Q_b)\right] + 2h \sigma_i \, \hat e_i \cdot \hat h .

At low temperature this dynamics freezes solid. Any flip out of an ice state costs 4J, so the ice manifold, with its exponential number of states, is effectively cut into pieces you can't move between. Barkema and Newman introduced loop moves for ice models to fix this [20], and Melko, den Hertog and Gingras adapted them to spin ice [21]. The lab uses the short-loop version. Start at a random vertex, follow a randomly chosen outgoing arrow and keep walking until the path comes back to a vertex it has visited. Then reverse every arrow on the closed loop. Reversing a directed loop swaps one in and one out at each vertex on it, so every charge, and so the ice energy, is unchanged. With a field the loop is accepted with the Metropolis probability of its Zeeman energy. Each vertex has as many outgoing arrows after the reversal as before, so the reverse loop is proposed with the same probability, which is detailed balance. The tests check this by enumeration: the loops visit all 148 ice states of a 3 \times 3 torus and all 90 of a single pyrochlore cell equally often.

Each sweep is N single-flip attempts followed, when loops are on, by V/2 loop attempts on a lattice with V vertices. The square lattice has N = 2L^2 spins; the pyrochlore has N = 16 L^3 in L^3 cubic cells.

The entropy measurement cools through 40 temperatures from T = 12J to 0.2J on a geometric grid. Above T_0 = 12J the high-temperature series is exact enough to start from:

\frac{\ln Z}{N} = \ln 2 + 3 \ln \cosh \beta J + \tfrac12 \ln\left(1 - 4t^3 + 3t^4\right) + c\, t^4 + \ldots, \qquad t = \tanh \beta J ,

where the logarithm counts the triangles and squares inside each tetrahedron and c = \tfrac12 counts the square plaquettes of square ice (c = 0 on the pyrochlore, whose first loop between tetrahedra is a hexagon). Then

S(T) = S(T_0) - \int_T^{T_0} \frac{C(T')}{T'} \, dT' = S(T_0) - \int_{E(T)}^{E(T_0)} \frac{dE}{T'} .

The two integrals are equal because C\,dT = dE. The lab does the second with the trapezoid rule on the measured mean energies, which settle far faster than the energy fluctuations that give C. On an 8 \times 8 torus the result lands within 0.005 of the exact count from a transfer matrix, while Pauling's value is 0.02 away.

Things to try

• The ice rules. Loops redraw the arrows every frame, the overlap with 64 sweeps ago drops to zero and the energy stays at exactly -J per spin. Every frame is a different ground state.
• Heat it: monopoles appear. At T = J about 15% of vertices carry a charge. Compare the measured fraction with the single-tetrahedron estimate below it.
• Quench: free monopoles wander. Random arrows dropped to T = 0.45J. The charges hop at zero cost and annihilate when opposite charges meet, until the few that thermal pair creation sustains remain, wandering as free defects.
• Drag a monopole. Press any vertex and drag. A pair appears and the charge you hold flips every arrow along its path, drawn bold. The energy rises by 4J for the pair, however long the string. Release it and watch it wander, extending or erasing its own string.
• Frozen without loops. The same ground state with single flips only: nothing moves and the overlap stays at 1. Switch moves to loops and it melts at once. Then switch loops on in the drag experiment: the string wanders and deforms, but its ends stay put. Only the monopoles are physical.
• Pinch points in S(q). The averaged structure factor grows sharp bow ties at (2\pi, 0), (0, 2\pi) and (2\pi, 2\pi). Raise T to 1.5J and they fill in as monopoles screen the constraint.
• Pyrochlore spin ice. One panel per cubic cell seen down [001], and the (h\,h\,l) plane of S(\mathbf q) with bow ties at (0\,0\,2), (1\,1\,1) and (2\,2\,0), the plane in which Fennell and collaborators resolved pinch points [19]. Then press Measure the residual entropy. The result should sit near Pauling's 0.2027; Nagle's series gives 0.2049 for the infinite lattice [6], a difference this run cannot resolve. On square ice the same sweep lands near Lieb's 0.2158, plus a finite-size excess of about 0.54/L^2 that the transfer-matrix counts in the tests show.
• Kagome ice in a [111] field. The field pins the spin on each tetrahedron that lies along [111]. The other three spins live on kagome planes and still obey the ice rule, with entropy left over: kagome ice [22, 23, 24]. Its magnetisation is \tfrac13 per spin, two thirds of saturation. Sweep a [111] field shows the plateau and the jump to the fully polarised, all-monopole state at h = 6J, where the Zeeman gain of h/6 per spin pays for J per spin of broken ice rules.

Reading the lab's instruments

• The lattice. Arrows are the spins. Blue dots are vertices with three in and one out (Q = +1), red dots one in and three out (Q = -1); larger dots are double charges. Bold arrows differ from the snapshot taken when you first pressed the lattice: the Dirac strings. On the pyrochlore each panel is one cubic cell viewed down [001], and each crossed square is a tetrahedron, drawn in the panel that holds its centre.
• Energy / spin. H/N in units of J, Zeeman term included: -1 in any ice state at h = 0, 0 for independent random spins.
• Charged sites. The fraction of vertices with Q \neq 0.
• Single tetrahedron. Pauling's approximation applied at finite T: each tetrahedron independent, with 6 ice states at -2J, 8 singly charged at 0 and 2 doubly charged at 6J. Its charged fraction is (8 + 2e^{-6\beta J}) / (6e^{2\beta J} + 8 + 2e^{-6\beta J}), and it neglects all interactions between tetrahedra, including the entropic Coulomb one. Shown only at h = 0.
• Net charge. \sum_t Q_t, always 0 on a periodic lattice: every flip moves a unit of charge from one end of a bond to the other.
• Specific heat. C/N = N(\langle e^2 \rangle - \langle e \rangle^2) / T^2 from the energy fluctuations since the last change, after 50 sweeps of equilibration.
• Magnetization. \frac1N \sum_i \sigma_i \hat e_i \cdot \hat h along [11] for square ice and [111] for the pyrochlore.
• Overlap, 64 sweeps. \frac1N \sum_i \sigma_i(t) \sigma_i(t - 64): 1 when nothing moves, near 0 when the configuration has decorrelated.
• Structure factor. S(\mathbf q) = \frac1N \langle |\mathbf M_\perp(\mathbf q)|^2 \rangle, the part of the Fourier-transformed moment perpendicular to \mathbf q, which is what unpolarised neutrons see. Square ice shows q_x, q_y \in [-4\pi, 4\pi]; the pyrochlore shows the (h\,h\,l) plane used by Fennell and collaborators [19]. Only wavevectors that fit the periodic box are computed, so small lattices give coarse maps. The colour ramp takes equal steps in CIELAB, so equal colour differences mean equal changes in S.
• Live plots. Energy, overlap and charged fraction with the estimate dashed, sweep by sweep.
• Sweep plots. Dots are this run, lines the single-tetrahedron estimate, against \log_{10} T. The entropy panel adds Lieb's or Nagle's value dashed. The M(h) panel compares with the T = 0 ground state of this model.

Order by disorder, and what the model leaves out

This site is named after order by disorder, where fluctuations pick one ordered state out of a degenerate ground manifold. Villain and collaborators showed it for thermal fluctuations in a frustrated Ising model [25]. Spin ice is the counterexample that teaches the most. Its manifold is so large and so evenly connected that fluctuations pick nothing, and what's left is a liquid with an emergent gauge field. Compare the frustrated triangular antiferromagnet in the Ising lab, which also keeps entropy at T = 0 without ordering.

Nature is less tidy than the model. The dipolar interaction selects an ordered ground state at very low temperature [21], and Pomaranski and collaborators found that Dy_2Ti_2O_7, given long enough to equilibrate, does release its residual entropy [26]. Real spin ice also falls out of equilibrium as it cools, its spins slowing to a crawl [27], much as the single-flip dynamics here freezes without loops. Castelnovo, Moessner and Sondhi's later review collects these threads [17].

References

[1] J. D. Bernal and R. H. Fowler, "A Theory of Water and Ionic Solution, with Particular Reference to Hydrogen and Hydroxyl Ions," Journal of Chemical Physics 1, 515-548 (1933). DOI: 10.1063/1.1749327

[2] J. C. Slater, "Theory of the Transition in KH_2PO_4," Journal of Chemical Physics 9, 16-33 (1941). DOI: 10.1063/1.1750821

[3] R. J. Baxter, Exactly Solved Models in Statistical Mechanics (Academic Press, London, 1982). Ch. 8, "Ice-Type Models."

[4] L. Pauling, "The Structure and Entropy of Ice and of Other Crystals with Some Randomness of Atomic Arrangement," Journal of the American Chemical Society 57, 2680-2684 (1935). DOI: 10.1021/ja01315a102

[5] W. F. Giauque and J. W. Stout, "The Entropy of Water and the Third Law of Thermodynamics. The Heat Capacity of Ice from 15 to 273°K," Journal of the American Chemical Society 58, 1144-1150 (1936). DOI: 10.1021/ja01298a023

[6] J. F. Nagle, "Lattice Statistics of Hydrogen Bonded Crystals. I. The Residual Entropy of Ice," Journal of Mathematical Physics 7, 1484-1491 (1966). DOI: 10.1063/1.1705058

[7] 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

[8] P. W. Anderson, "Ordering and Antiferromagnetism in Ferrites," Physical Review 102, 1008-1013 (1956). DOI: 10.1103/PhysRev.102.1008

[9] 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

[10] 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

[11] B. C. den Hertog and M. J. P. Gingras, "Dipolar Interactions and Origin of Spin Ice in Ising Pyrochlore Magnets," Physical Review Letters 84, 3430-3433 (2000). DOI: 10.1103/PhysRevLett.84.3430

[12] E. H. Lieb, "Exact Solution of the Problem of the Entropy of Two-Dimensional Ice," Physical Review Letters 18, 692-694 (1967), DOI: 10.1103/PhysRevLett.18.692; "Residual Entropy of Square Ice," Physical Review 162, 162-172 (1967), DOI: 10.1103/PhysRev.162.162

[13] R. F. Wang, C. Nisoli, R. S. Freitas, J. Li et al., "Artificial 'spin ice' in a geometrically frustrated lattice of nanoscale ferromagnetic islands," Nature 439, 303-306 (2006). DOI: 10.1038/nature04447

[14] C. Castelnovo, R. Moessner and S. L. Sondhi, "Magnetic monopoles in spin ice," Nature 451, 42-45 (2008). DOI: 10.1038/nature06433. arXiv: 0710.5515

[15] D. J. P. Morris, D. A. Tennant, S. A. Grigera, B. Klemke et al., "Dirac Strings and Magnetic Monopoles in the Spin Ice Dy_2Ti_2O_7," Science 326, 411-414 (2009). DOI: 10.1126/science.1178868

[16] C. L. Henley, "The 'Coulomb Phase' in Frustrated Systems," Annual Review of Condensed Matter Physics 1, 179-210 (2010). DOI: 10.1146/annurev-conmatphys-070909-104138. arXiv: 0912.4531

[17] C. Castelnovo, R. Moessner and S. L. Sondhi, "Spin Ice, Fractionalization, and Topological Order," Annual Review of Condensed Matter Physics 3, 35-55 (2012). DOI: 10.1146/annurev-conmatphys-020911-125058. arXiv: 1112.3793

[18] S. V. Isakov, K. Gregor, R. Moessner and S. L. Sondhi, "Dipolar Spin Correlations in Classical Pyrochlore Magnets," Physical Review Letters 93, 167204 (2004). DOI: 10.1103/PhysRevLett.93.167204. arXiv: cond-mat/0407004

[19] T. Fennell, P. P. Deen, A. R. Wildes, K. Schmalzl et al., "Magnetic Coulomb Phase in the Spin Ice Ho_2Ti_2O_7," Science 326, 415-417 (2009). DOI: 10.1126/science.1177582. arXiv: 0907.0954

[20] G. T. Barkema and M. E. J. Newman, "Monte Carlo simulation of ice models," Physical Review E 57, 1155-1166 (1998). DOI: 10.1103/PhysRevE.57.1155

[21] R. G. Melko, B. C. den Hertog and M. J. P. Gingras, "Long-Range Order at Low Temperatures in Dipolar Spin Ice," Physical Review Letters 87, 067203 (2001). DOI: 10.1103/PhysRevLett.87.067203. arXiv: cond-mat/0009225

[22] K. Matsuhira, Z. Hiroi, T. Tayama, S. Takagi and T. Sakakibara, "A new macroscopically degenerate ground state in the spin ice compound Dy_2Ti_2O_7 under a magnetic field," Journal of Physics: Condensed Matter 14, L559-L565 (2002). DOI: 10.1088/0953-8984/14/29/101

[23] M. Udagawa, M. Ogata and Z. Hiroi, "Exact Result of Ground-State Entropy for Ising Pyrochlore Magnets under a Magnetic Field along [111] Axis," Journal of the Physical Society of Japan 71, 2365-2368 (2002). DOI: 10.1143/JPSJ.71.2365

[24] R. Moessner and S. L. Sondhi, "Theory of the [111] magnetization plateau in spin ice," Physical Review B 68, 064411 (2003). DOI: 10.1103/PhysRevB.68.064411

[25] J. Villain, R. Bidaux, J.-P. Carton and R. Conte, "Order as an effect of disorder," Journal de Physique 41, 1263-1272 (1980). DOI: 10.1051/jphys:0198000410110126300

[26] D. Pomaranski, L. R. Yaraskavitch, S. Meng, K. A. Ross et al., "Absence of Pauling's residual entropy in thermally equilibrated Dy_2Ti_2O_7," Nature Physics 9, 353-356 (2013). DOI: 10.1038/nphys2591

[27] J. Snyder, J. S. Slusky, R. J. Cava and P. Schiffer, "How 'spin ice' freezes," Nature 413, 48-51 (2001). DOI: 10.1038/35092516