# A Three Dimensional Kasteleyn Transition: Spin Ice in a [100] Field

Ludovic D. C. Jaubert, J. T. Chalker, Peter C.W. Holdsworth, R. Moessner

► To cite this version:

Ludovic D. C. Jaubert, J. T. Chalker, Peter C.W. Holdsworth, R. Moessner. A Three Dimensional Kasteleyn Transition: Spin Ice in a [100] Field. Physical Review Letters, 2008, 100, pp.067207. 10.1103/PhysRevLett.100.067207. ensl-00176660v4

HAL Id: ensl-00176660

https://ens-lyon.hal.science/ensl-00176660v4

Submitted on 3 Mar 2008

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A Three Dimensional Kasteleyn Transition: Spin Ice in a [100] Field

L. D. C. Jaubert,¹ J. T. Chalker,² P. C. W. Holdsworth,¹ and R. Moessner²,³

¹Laboratoire de Physique, École Normale Supérieure de Lyon,
Université de Lyon, CNRS, 46 Allée d'Italie, 69364 Lyon Cedex 07, France.

²Theoretical Physics, Oxford University, 1 Keble Road, Oxford, OX1 3NP, United Kingdom.

³Max-Plack-Institut für Physik komplexer Systeme, 01189 Dresden, Germany

(Dated: March 3, 2008)

We examine the statistical mechanics of spin-ice materials with a [100] magnetic field. We show that the approach to saturated magnetisation is, in the low-temperature limit, an example of a 3d Kasteleyn transition, which is topological in the sense that magnetisation is changed only by excitations that span the entire system. We study the transition analytically and using a Monte Carlo cluster algorithm, and compare our results with recent data from experiments on Dy₂Ti₂O₇.

PACS numbers: 75.10.Hk, 75.40.Cx, 75.40.Mg

Topological phases occur in connection with new liquid states and concurrent exotic transitions out of them. A central idea is the emergence of degrees of freedom which cannot be defined or manipulated locally. Their presence has profound ramifications for the static and dynamic properties of such phases. They can exhibit algebraic correlations in the absence of criticality [1] and unusually slow relaxation towards thermal equilibrium [2, 3].

A magnetic system whose topological properties have a particularly simple and transparent origin is known as spin ice [4, 5] (Fig. 1): in its macroscopically degenerate ground state, a magnetic version of the ice rules apply. These stipulate that for each tetrahedron of the corner-sharing pyrochlore lattice, two of the four spins point in, and two point out of the tetrahedron along the local body centred diagonal axes. As a result, the lattice divergence of the spin field vanishes everywhere: $\nabla \cdot \vec{S} = 0$. As a consequence of this local constraint, the magnetisation of each (100) plane is the same [6] and is therefore a topological quantity which can be changed only by making a simultaneous change throughout the system: the smallest such magnetic excitation involves a set of spins on a string spanning the system. An attractive feature of spin-ice materials is that one can couple directly to this topological quantity using a uniform magnetic field, which lifts the macroscopic degeneracy of the zero-field ground states.

In this paper, we focus on the behaviour near saturation, when the ordered $q = 0$ state [5, 7] is reached by applying a magnetic field along the [100] direction as illustrated in Fig. 1. In the saturated state at large fields, the ice rules prohibit local fluctuations in the limit of low temperature $T$. The resulting phase transition as the field is lowered is a three dimensional (3d) example of a Kasteleyn transition [8], first proposed in the context of 2d dimer models [9]. Its hallmark is an asymmetric character, appearing to be first-order on one side, and continuous on the other. The Kasteleyn transition has in the past been associated with soft matter: the trans-gauche disordering transition for polymers embedded in

<figure><img src="image_2.png" /><figcaption>FIG. 1: (Color online) Left: Pyrochlore lattice showing $q = 0$ spin structure (black narrow spins) and a string defect (red thick spins). Right: Bethe lattice with a central tetrahedron (red) and first layer of nearest-neighbour tetrahedra (blue).</figcaption></figure>

lipid bilayer system, observed via density measurements in dipalmitoyl lecithin (DPL) [10]. To our knowledge this is the first example of a three-dimensional Kasteleyn transition in a magnetic system.

We analyse the transition using a combination of analytical and numerical tools. By mapping the strings onto bosonic worldlines, we obtain an analytical expression for thermodynamic quantities close to the critical point, and for correlation functions at all values of the magnetic field. These predictions are confirmed to high precision by numerical simulation using a Monte Carlo algorithm incorporating a non-local cluster update. In addition, we find that a calculation on a Bethe lattice provides an excellent description of the thermodynamics of the magnet away from the critical point. A crucial characteristic of this transition is that it is symmetry-sustaining, as recognised in Ref. [11]; the first-order liquid-gas transition reported there, however, is an artefact of the local dynamics employed at that time. We compare our results with magnetisation measurements on the spin ice compound Dy₂Ti₂O₇ [12]: they represent well the observed behaviour at large field and temperature, but for lower values the experimental system falls out of equilibrium.

We consider a spin ice model with nearest-neighbour
magnetic interactions only. The Hamiltonian reads

$$H = -J \sum_{\langle ij \rangle} \vec{S}_i \cdot \vec{S}_j - \vec{h} \cdot \sum_i \vec{S}_i, \quad (1)$$

where $\vec{S}_i$ is a vector of unit length, encoding an Ising spin constrained to lie along the local body centred cubic axes of the tetrahedron; $J > 0$ is the ferromagnetic exchange constant and $\vec{h}$ is an external field of strength $h$ in reduced units, aligned along the [100] direction. This model is a simplification in that the real Hamiltonian is known to incorporate dipolar interactions of at least equal importance [13, 14]. However, it has been demonstrated that a nearest-neighbour description is adequate as long as the ice rules are not violated, since corrections fall off with separation $r$ as $r^{-5}$ [15]. Our results should therefore be robust for $T \ll J$, when behaviour depends only on the ratio $h/T$, whereas interactions between further neighbours may become important at higher $T$.

Useful insights into the Kasteleyn transition come from discussing in detail the string excitations introduced above: start from the $q = 0$ state, flip spins in successive [100] layers to generate another configuration, and link flipped spins in adjacent layers with segments of the strings. This mapping between spin and string configurations is one-to-one provided we adopt a fixed convention for the paths followed by two strings through the same tetrahedron. Strings of finite length do not occur in ice-rule states as their ends cost exchange energy $O(J)$. In contrast, strings that span the system cost only Zeeman energy, of $2h/\sqrt{3}$ per segment. An isolated string also has an entropy, of $\ln 2$ per segment, and a finite density of strings can account for the macroscopic ground state entropy of spin ice. The total free energy for a single spanning string of $L$ segments is

$$G = L \left( 2h/\sqrt{3} - T \ln 2 \right). \quad (2)$$

From this one obtains the exact value of the Kasteleyn transition temperature, $T_K = 2h/\sqrt{3}\ln 2$, and for $h, T \ll J$, properties are a function only of $T/T_K$. For $T < T_K$ the equilibrium phase is the $q = 0$ state without any thermal fluctuations, while for $T > T_K$ strings are present. Depending on their interactions, different scenarios are possible for the variation of their density, and hence sample magnetisation, with $T/T_K$. If strings attract, their density has a jump at the transition. At finite string density, if they form a commensurate crystal, there is a magnetisation plateau; and if they form a liquid, their magnetisation varies smoothly. In fact, by construction, strings have hard core repulsion, and the simulations we report below show only a liquid phase.

By thinking of the strings as world lines for hard core bosons moving in two dimensions at zero temperature, with the [100] direction as imaginary time, we can obtain thermodynamic behaviour near the critical point.

Under this mapping, the reduced temperature difference $t \equiv (T/T_K - 1)$ is proportional to the chemical potential for bosons, the deviation of the magnetisation $M$ per site from its saturation value $M_{sat}$ is proportional to boson density, and the spin-ice free energy is proportional to the boson ground-state energy. The Kasteleyn transition is between a vacuum state and a Bose condensate as the chemical potential changes sign. Here, three (2+1) is the upper critical dimension, since in the quantum description the dynamical critical exponent takes the value $z=2$, a reflection of the fact that the [100] field renders the classical model highly anisotropic. Transcription of established results for the boson problem [16] gives for $t > 0$

$$M = M_{\text{sat}}[1 - at - bt \ln(1/t)], \quad (3)$$

where $a$ and $b$ are constants. It follows from this that the heat capacity $C_h$ and the differential susceptibility $\chi$ diverge logarithmically in $t$ as the critical point is approached from the high-temperature side. Similar conclusions for a dimer model were reached using arguments based on the interaction between two strings [8, 17].

Correlation functions are of great interest throughout the partially magnetised phase, and at $h = 0$ are known to have an algebraic dipolar form that is a consequence of the ice rules [18]. The bosonic description of the system that we have introduced provides both a new way of viewing this result and also a means of finding how it changes with $h$. For a complete treatment, it is necessary to take into account the distinction between the four sites in the basis of the pyrochlore lattice when constructing the mapping between spins and strings. Moreover, a precise representation of the transfer matrix for strings in terms of a boson Hamiltonian would include multiparticle interactions, which are not expected to influence long-distance behaviour or critical properties, since both of these should be universal. We defer these details to a future publication and present here a simplified calculation, in which bosons move on a square lattice and the imaginary time direction is approximated as continuous. We start from a boson Hamiltonian written in terms of spin $S = 1/2$ operators as

$$\mathcal{H} = -\mathcal{J} \sum_{\langle ij \rangle} (S_i^x S_j^x + S_i^y S_j^y) - \mathcal{B} \sum_i S_i^z, \quad (4)$$

where $\mathcal{J}$ is an effective coupling constant and $\mathcal{B} \propto h$. In this incarnation, the Kasteleyn transition is a quantum phase transition in which the two-dimensional XY ferromagnet is fully polarized by a field $\mathcal{B}$ coupled to the $z$-component of magnetisation, at a critical field strength $\mathcal{B} = \mathcal{B}_c$. Moreover, $M$ for spin ice is proportional to $\langle S_i^z \rangle$ for the quantum magnet. We use a large $S$ approximation to evaluate the connected correlation function, defined for spin ice as

$$C(\mathbf{r}, z) = \langle \vec{S}_i \cdot \vec{S}_j \rangle - \langle \vec{S}_i \rangle \cdot \langle \vec{S}_j \rangle, \quad (5)$$
<figure><img src="image_3.png" /><figcaption>FIG. 2: (Color online) M vs. T/T$_K$ obtained from simulations for the pyrochlore lattice (dots) and analytically on the Bethe lattice (solid lines) for h/J = 10$^{-3}$(•), 0.13(▲) and 0.58(■).</figcaption></figure>

where the separation between lattice sites i and j is r within the (100) plane and z along the [100] direction. Since we have treated the lattice structure in a simplified way, our results apply only to intra-sublattice correlations. Writing B/B$_c$ = cos θ, we find for large r and z

$$
C(\mathbf{r}, z) = \frac{S \sin \theta [3[z \sin \theta]^2 - (r^2 + [z \sin \theta]^2)]}{8\pi (r^2 + [z \sin \theta]^2)^{5/2}} . \quad (6)
$$

For B = 0 the form of this expression reduces to that of Ref. [18] and has the cubic symmetry of the zero field problem. For finite field this symmetry is broken and distances in the z direction acquire a scale factor 1/√1 − ( B/B$_{c}$ ) $^{2}$. Close to the Kasteleyn transition, there is a second regime of behaviour, which also follows from Eq. (4) but is not displayed in Eq. (6). It arises because correlations at distances shorter than the string separation are controlled by the behaviour of an isolated string, and hence reduce to the correlations along a 2d random walk, giving for example, C(0, z) ∼ z − $^{1}$.

We have investigated the transition via Monte Carlo (MC) simulations, using a cluster algorithm [19, 20, 21]. This allows us to avoid the apparent first-order discontinuities [11] that arise from loss of ergodicity with single spin flip dynamics for T ≪ J. Each MC move consists of the reversal of all spins on a string, which in the limit T/J → 0 and with periodic boundary conditions must close on itself. The algorithm is efficient because in this limit all attempted flips are accepted below saturation. A typical simulation involves 4 × 10 6 spins, studied for 10 6 MC steps at each temperature. Monte Carlo data for magnetisation (Fig. 2), differential susceptibility and heat capacity (Fig. 3) show conclusively a 3d Kasteleyn transition in the limit h ≪ J. For T < T$_K$, the system is completely frozen, with M = M$_{sat}$ and C$_h$ = χ = 0, while for T > T$_K$, M varies continuously with T/T$_K$. The behaviour of C$_h$ and χ close to T$_K$, related by C$_h$ = h $^{2}$χ/T , is consistent with the logarithmic singularity expected

<figure><img src="image_4.png" /><figcaption>FIG. 3: C$_h$ and χ vs t ≡ (T/T$_K$ − 1) obtained from simulations (•) and analytically (solid lines) for h = 10 − $^{3}$J . The dashed line is a fit to a logarithmic singularity derived from Eq. (3).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4: (Color online) Finite-size scaling of spin correlations for h/h c = 0.33(■), 0.66(▲) and 0.93(•), on log-log scale. Lines have slopes 3 and 1.</figcaption></figure>

from Eq. (3). Away from the limit h ≪ J, there is no sharp transition: instead M varies smoothly with T .

Data for the correlation function, illustrated in Fig. 4, show the behaviour expected from Eq. (6). At the maximum separation along the [100] direction in a system of linear size L (r=0 and z=L/2) correlations decrease with increasing L as L − $^{3}$, and increase in amplitude with h/h$_{c}$ . Close to the transition, a slower decrease is observed at short distances, consistent with the L − 1 law expected for isolated strings.

We supplement our simulations with analytical calculations for a related model: a Bethe lattice (BL) of tetrahedra, whose central element is illustrated in Fig. 1. As on the pyrochlore lattice, each tetrahedron is connected to four others. Crucially, however, there are no closed loops of tetrahedra [22]: the statistical mechanics can then be solved exactly and a Kasteleyn transition appears at T$_{K}$ . Analytical results for the BL are very close to the simulation data for the pyrochlore lattice – in the
<figure><img src="image_6.png" /><figcaption>FIG. 5: Magnetisation *M* vs the internal magnetic strength *μ*₀*H*int obtained analytically for the BL (solid lines) and experimentally (points) for the compound Dy₂Ti₂O₇ at *T*exp = 1.8K; data from Ref. [12].</figcaption></figure>

absence of adjustable parameters (Fig. 2 and 3) – provid-
ing a *posteriori* justification for the BL approach. Small
differences are apparent only for *h* ≪ *J*: these arise be-
cause the transition on the BL is necessarily mean field,
and so does not have the logarithmic contributions to
critical behaviour of Eq. (3).

What about the experimental situation? Studies of magnetization as a function of [100] field strength at fixed low temperature have been reported for the spin ice materials Dy₂Ti₂O₇ [12] and Ho₂Ti₂O₇ [23]. Here we present a comparison of our theoretical results with data for Dy₂Ti₂O₇. The measured magnetisation flattens off abruptly at the saturation value, M<sub>exp</sub> = 10μ<sub>B</sub>/√3, in a way reminiscent of the Kasteleyn transition. For temperatures well below the scale set by the exchange constant, experiments indicate that spin ice materials are out of equilibrium, as expected if the actual dynamics is local. We therefore make a comparison at a temperature comparable with the exchange energy. At such temperatures the behaviour of the pyrochlore lattice and BL are indistinguishable. We show in Fig. 5 data for Dy₂Ti₂O₇ at T<sub>exp</sub> = 1.8K (Ref. [12]) together with our BL results at a fitting temperature T = 0.87 T<sub>exp</sub>. This temperature is about 1.6 times the effective nearest-neighbour exchange constant J<sub>eff</sub>S<sup>2</sup>/k<sub>B</sub> ~ 1.11K estimated for Dy₂Ti₂O₇ [12]. The agreement between experiment and our theory is clearly good, indicating that the former is indeed close to the transition discussed here.

In conclusion, we propose that spin ice materials pro-
vide a new opportunity for the study of the Kasteleyn
transition. The 3d transition presented here differes from
that in 2d, predicted to occur [24] for spin ice materials in
a [111] field, both in its critical properties and in the be-
havior of correlation functions. In 2d, it takes the form of
a commensurate-incommensurate phase transition which

can be evidenced by a divergent peak in the correlation
function (absent in 3d), the wavevector of which is field-
dependent [24], as reported in recent measurements [25].
We believe that the data of Refs. [12] and [23] provide ev-
idence for rounded 3d Kasteleyn transitions in Dy2Ti2O7
and Ho2Ti2O7 respectively, which could be more closely
approached in future experiments.

It is a pleasure to thank S. T. Bramwell, C. Castelnovo,
T. Fennell and M. J. Harris for useful discussions. We ac-
knowledge financial support from the European Science
Foundation PESC/RNP/HFM (PCWH and LJ), from
the Royal Society and the London Centre for Nanotech-
nology (PCWH), from ANR grant 05-BLAN-0105 (LJ),
and from EPSRC Grant No. GR/83712/01 (JTC). LJ
and PCWH thank the Rudolph Peierls Centre for The-
oretical Physics, University of Oxford, where this work
was completed.

[1] R. W. Youngblood and J. D. Axe, Phys. Rev. B **23**, 232 (1981).

[2] C. L. Henley, J. Stat Phys. **89**, 483 (1997).

[3] B. Chakraborty, D. Das, and J. Kondev, Eur. Phys. J. E **9**, 227 (2002).

[4] P. W. Anderson, Phys. Rev. **102**, 1008 (1956).

[5] M. J. Harris et al., Phys. Rev. Lett. **79**, 2554 (1997).

[6] R. Moessner and J. T. Chalker, Phys. Rev. B **58**, 12049 (1998).

[7] A. P. Ramirez et al., Nature **399**, 333 (1999).

[8] S. M. Bhattacharjee et al., J. Stat. Phys. **32**, 361 (1983).

[9] P. W. Kasteleyn, J. Math. Phys. **4**, 287 (1963).

[10] J. F. Nagle, Proc. Nat. Acad. Sci. USA **70**, 3443 (1973).

[11] M. J. Harris et al., Phys. Rev. Lett. **81**, 4496 (1998).

[12] H. Fukazawa et al., Phys. Rev. B **65**, 054410 (2002).

[13] R. Siddharthan et al., Phys. Rev. Lett. **83**, 1854 (1999).

[14] B. C. den Hertog and M. J. P. Gingras, Phys. Rev. Lett. **84**, 3430 (2000).

[15] S. V. Isakov, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. **95**, 217201 (2005).

[16] D. S. Fisher and P. C. Hohenberg, Phys. Rev. B **37**, 4936 (1988).

[17] J. F. Nagle, C. S. O. Yokoi, and S. M. Bhattacharjee, *Phase Transitions and Critical Phenomena, Ed. C. Domb and J. L. Lebowitz* (Academic Press, 1989), chap. 2.

[18] S. V. Isakov et al., Phys. Rev. Lett. **93**, 167204 (2004).

[19] G. T. Barkema and M. E. J. Newman, Phys. Rev. E **57**, 1155 (1998).

[20] R. G. Melko, B. C. den Hertog, and M. J. P. Gingras, Phys. Rev. Lett. **87**, 067203 (2001).

[21] S. V. Isakov et al., Phys. Rev. B **70**, 104418 (2004).

[22] P. Chandra and B. Doucot, J. Phys. A **27**, 1541 (1994).

[23] T. Fennell et al., Phys. Rev. B **72**, 224411 (2005).

[24] R. Moessner and S. L. Sondhi, Phys. Rev. B **68**, 064411 (2003).

[25] T. Fennell et al., Nature Physics **3**, 566 (2007).