# The Kasteleyn transition in three dimensions: 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. The Kasteleyn transition in three dimensions: spin ice in a [100] field. HIGHLY FRUSTRATED MAGNETISM 2008, Sep 2008, Braunschweig, Germany. 012024 (7pp), 10.1088/1742-6596/145/1/012024. ensl-00366133

HAL Id: ensl-00366133

https://ens-lyon.hal.science/ensl-00366133

Submitted on 5 Mar 2009

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# The Kasteleyn transition in three dimensions: spin ice in a [100] field

L. D. C. Jaubert$^{†}$, J.T. Chalker$^{‡}$, P. C. W. Holdsworth$^{†}$ and R. Moessner*

$^{†}$Université de Lyon, Laboratoire de Physique, École Normale

Supérieure de Lyon, 46 Allée d'Italie, 69364 Lyon cedex 07, France.

$^{‡}$Rudolf Peierls Centre for Theoretical Physics, 1 Keble Rd., Oxford, OX1 3NP, UK and

$^{*}$Max-Plack-Institut für Physik komplexer Systeme, 01189 Dresden, Germany.

(Dated: March 5, 2009)

We discuss the nearest neighbour spin ice model in the presence of a magnetic field placed along the cubic [100] direction. As recently shown in Phys. Rev. Lett. 100, 067207, 2008, the symmetry sustaining ordering transition observed at low temperature is a three dimensional Kasteleyn transition. We confirm this with numerical data using a non-local algorithm that conserves the topological constraints at low temperature and from analytic calculations from a Bethe lattice of corner sharing tetrahedra. We present a thermodynamic description of the Kasteleyn transition and discuss the relevance of our results to recent neutron scattering experiments on spin ice materials.

## I. INTRODUCTION

This transition was introduced by Kasteleyn [1] in the context of ordering of hard core 2D dimers lying on the bonds of a honeycomb lattice. Applying local chemical potentials, $\mu_i$, for dimers oriented in one of the $i = 1, 2, 3$ possible directions leads to a singular ordering transition for finite values of the fugacities $z_i = \exp(-\beta\mu_i)$. The transition can also occur for dimers on three dimensional equivalents of the honeycomb lattice such as the brick lattice [2], while the exponents can be determined in dimension $d$ through mapping onto a directed polymer problem [3]. It has been extensively used to describe trans-gauche structural transitions in polymerized lipid bilayers and the resulting theory provides a good qualitative description of the singular density change occurring in solvent bilayer systems as a function of temperature [4]. A Kasteleyn transition also occurs in the five vertex model on a square lattice [5], which maps onto the square ice model, satisfying the ice rules in two dimensions with direct relevance to spin ice materials [6]. It has been shown that the same Kasteleyn transition exists in the nearest neighbour model for spin ice materials on a pyrochlore lattice, with magnetic field directed close to the body centred cubic [111] direction [7]. Placing the field along this direction isolates layers orthogonal to the field, each of which forms a kagomé lattice of directed tetrahedra whose vertices lie alternately above and below the plane. The magnetic field fixes only one out of the four sublattices with the result that the extensive ground state degeneracy of zero field is only partially lifted. The ground state entropy remains finite and the kagomé spin problem maps exactly onto Kasteleyn’s original dimer problem on the Honeycomb lattice. Tilting the field off the [111] direction, giving a finite perpendicular field component $\vec{H}_{\perp}$, is equivalent to applying chemical potentials $\mu_i$ to the dimers and the spins order in the kagomé planes for a finite value of $\vec{H}_{\perp}/k_B T$ via Kasteleyn’s transition. Ground state configurations of this model require the imposition of local constraints, which lead to power law dipole like correlations [8, 9]. Recent neutron studies of $Ho_2Ti_{27}$ with this field configuration [10] have shown experimentally the existence of these correlations in the form of well defined pinch points in reciprocal space for the intensity of the elastic diffuse scattering. In two dimensions these points should include logarithmically divergent intensities. In the presence of a perpendicular field component the singularities are predicted [7] to drift with changing field strength, $\vec{H}_{\perp}$. It has recently been shown that the ordering transition seen when a field is placed along the [100] direction [11] is a three dimensional Kasteleyn transition [12]. Magnetic measurements on Holmium and Dysprosium Titanate [13, 14] indicate that experimental systems could show strong features of the transition, but as we will see below, accessing it while remaining in thermodynamic equilibrium is a considerable experimental challenge.

## II. THE KASTELEYN TRANSITION IN SPIN ICE IN A [100] FIELD

In this paper we review the behaviour of the nearest neighbour spin ice model, with magnetic field placed along the cubic [100] direction. Nearest neighbour moments (spins) of unit length are constrained to lie along the body centred crystal field directions of the tetrahedra: $\vec{S}_i = \pm \vec{d}_i$ and interact with a ferromagnetic exchange interaction $J$ and with an external magnetic field, $\vec{H}$ in reduced units. The Hamiltonian [6], is

$$ \mathcal{H} = - \sum_{\langle ij \rangle} J \vec{S}_i \cdot \vec{S}_j - \vec{H} \cdot \sum_i \vec{S}_i. \qquad (1) $$
The crystal fields frustrate the ferromagnetic exchange leading to a highly degenerate ground state satisfying the so called ice rules [15] with two spins pointing in and two spins pointing out of each tetrahedron. This is the topological constraint which takes the form $\nabla \cdot \vec{S}(\vec{r}) = 0$ [8] in a continuum description. In zero field there is no phase transition down to zero temperature but, in the presence of a magnetic field and in the limit where the topological constraints are rigorously applied the system orders at low temperature via a three dimensional [12] Kasteleyn transition.

The remarkable thing about the Kasteleyn transition is that it occurs when the system is rigorously constrained to the manifold of spin ice states such that the internal energy for each microstate in the reduced phase space has the same value, $U = -NJ/3$. This allows for a particularly simple thermodynamic description which we present below. Just as for a paramagnet the magnetic Helmholtz free energy for the constrained system has an entropic component only: $F(N, T, M) = -TS(N, T, M)$, where M is the magnetic order parameter. In a regular paramagnetic system S goes to zero, as M reaches its maximum value, $M_{max}$ but it does so with infinite slope and hence an infinite value for $H/k_BT$, where $H = |\vec{H}|$ is the field conjugate to M:

$$H = -\frac{1}{N} \frac{\partial F}{\partial M} \qquad (2)$$

That is, there is no phase transition. However the hidden divergence free constraint can drive the entropy to zero for finite $H/k_BT$ giving a singular longitudinal susceptibility, $\chi = \partial M/\partial H$-this is the Kasteleyn transition. The Gibbs potential is defined $G^* = F - NMH$ from which one can construct a Landau free energy, $G_L$ by expanding $G^*$ in small $x = M_{max} - M$ near the transition:

$$\frac{G_L}{N} = (H - H_K)x + \frac{\alpha_2}{2}Tx^2 + \frac{\alpha_3}{3}Tx^3 + \dots - M_{max}h, \qquad (3)$$

where $H_K$ is the critical field for T fixed and where $\alpha_i$ are constants. If $\alpha_2 > 0$ one finds the mean field critical behaviour for fixed H:

$$x = x_0(T - T_K)^\beta, \quad \beta = 1, \qquad (4)$$

where $T_K$ is the transition temperature for fixed H. One also finds that the specific heat, $C_H$, and susceptibility, $\chi$, jump discontinuously to zero at $T_K$. For $T < T_K$ $G^*(M)$ has no minimum for $0 < M < M_{max}$. $\partial G^*/\partial M$ is finite and negative at $M = M_{max}$, consistently with there being zero fluctuations and the entropy being strictly zero: the system is perfectly ordered and represents a complete vacuum for excitations. Hence the transition is violently asymmetric, with fluctuations above the transition but not below. Through this asymmetry the transition has often been described as having both 1st and 2nd order characteristics[16]. However from the above thermodynamics it is clearly second order, or continuous, within the Ehrenfest classification[17].

Zero entropy can only be maintained at finite temperature if the topological constraints are rigorously applied ($J \gg k_BT$) so that the lowest energy excitations off the perfectly ordered state are extended, with infinitely large energy in the thermodynamic limit. To see that this is the case consider a strong field at low temperatures which moves the system from a disordered spin ice state to the long range ordered state illustrated by the black spins in Fig 1(a). Relaxing back to a disordered state while maintaining the ice rules requires the flipping of a string of spins spanning the entire system (in the following discussion we impose periodic boundary conditions, in which case the string closes through the boundaries). To see this consider flipping a single spin out of the ordered state. This breaks the ice rules creating a pair of "topological defects" [18], tetrahedra with three spins in and one out and vise versa, which can only be annealed by interacting with each other. Creating the defects costs magnetic enthalpy $4J/3 + 2H/\sqrt{3}$, but moving them further apart costs no further exchange energy in the nearest neighbour model. The exchange energy can be regained if the ensuing string passes right through the system and the topological defects destroy each other, returning the system to the constrained manifold of spin ice states. Such a string is illustrated by the red spins in Figure 1(a). Ignoring the exchange energy, each step of the string costs Zeeman energy $2H/\sqrt{3}$ but also represents an entropy gain of $\ln(2)$ as the string can leave the tetrahedron in one of two directions. If we define L as the number of spins in the [100] direction, the total free enthalpy change for placing a string in the system is then

$$\begin{aligned} \delta G &= L \left( \frac{2H}{\sqrt{3}} - T \ln 2 \right) \\ &\propto L (T_K - T), \quad \text{where} \\ T_K &= \frac{2H/\sqrt{3}}{\ln 2}, \end{aligned} \qquad (5)$$

and strings appear above the temperature $T_K$ only. The transition could be first order if strings attracted each other but two strings passing through the same tetrahedron clearly have a reduced entropy over that for two separated strings. Hence, strings repel each other and the result is a continuous Kasteleyn transition at temperature $T_K$.
<figure><img src="image_2.png" /><figcaption>FIG. 1: (a) The pyrochlore lattice with spins satisfying the ice rules. $\vec{H}$ is given along the [100] direction. Black spins: long range ordered state, red spins: a string of flipped spins.</figcaption></figure>

**Figure 1** (b) (Color online) Magnetization per spin $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}(\bullet)$, 0.13($\triangle$) and 0.58($\blacksquare$).

A consequence of the topological constraints and the broken magnetic symmetry is that a mapping exists between the strings of returned spins and world lines for bosons living in a $d-1$ dimensional space perpendicular to the field [12, 19]. The pyrochlore lattice has six fold cubic symmetry, which is reduced by the magnetic field, giving a unique definition for the arrow of time of the world lines along an imaginary time axis. The bosonic nature of the particles shows up in the passage of two strings through a single tetrahedron. The entropic contribution of the tetrahedron is zero and it is impossible to index strings "1" and "2" either side of the encounter without a fixed convention. However, with such a convention an exact mapping exists between spin and string configurations. In the absence of a field, there are six equivalent definitions for the arrow of time and a further convention is required to define the strings. In reference [12, 19] it was shown that through the definition of strings as world lines, the Kasteleyn transition in dimension $d$ maps onto a quantum phase transition in dimension $d-1$. The applied field corresponds to the boson chemical potential, $\mu \sim (h-h_K)$ and the transition is from a Bose condensed to a vacuum state at zero temperature. The upper critical dimension for this transition is $d=3$[20], which is compatible with previous analytic studies of the Kasteleyn transition in three dimensional dimer models[2]. Hence, for the three dimensional Kasteleyn transition we expect near mean field behaviour, with logarithmic corrections to scaling, while in two dimensions fluctuations should be important. This is indeed the case: the analytic solution of the two dimensional problem[7] gives a value for the exponent, $\beta=1/2$.

We have tested these ideas using a non-local Monte Carlo algorithm which flips connected lines of spins [12, 21]. The lines can form closed loops, either within the sample or through the boundaries, or they can transport topological defects non-locally. In the case of closed loops only, the algorithm restricts the system to fully constrained spin ice states, allowing an ergodic exploration of this manifold without passing through energetically unfavorable unconstrained states. In Figure 1(b) we show $M = \frac{1}{N} \sum_i \vec{S}_i$ vs $T$ for fixed $H$. The points represent data generated by simulation using the non-local loop algorithm, while solid lines show analytic data from a calculation on a tree structure (Bethe lattice) of connected tetrahedron. The lack of closed loops here, allows for an essentially analytic solution of the problem [22]. For $H/J = 0.001$ the data show the characteristics of the Kasteleyn transition: asymmetric singular behaviour, with the magnetisation arriving at its saturation value $M_{max} = 1/\sqrt{3}$ precisely at $T = T_K$. This is confirmed in figure 2(a), where the magnetisation difference, $x = M_{max} - M$ and the specific heat (and hence the susceptibility) are both seen to be singular. As predicted through the mapping to the quantum phase transition, the simulation data for $C_H$ fit accurately to a logarithmic divergence above the transition, while $x$ approaches zero linearly with $T - T_K$ and has logarithmic corrections for larger $x$:

$$x = \Delta M \propto t(1 - a \ln(t)) \Leftrightarrow C_h \propto -\ln(t) \qquad (6)$$

Detailed calculations for the Bethe lattice confirm the mean field picture developed phenomenologically above [22].

Previous numerical work found indications of a $1^{st}$ order phase transition at finite field [11]. The observed first order nature is a consequence of the loss of ergodicity due to the single spin flip dynamics used at that time. Interestingly experimental results [13] reflected the same out-of-equilibrium features, as natural dynamics is also local. However,
<figure><img src="image_3.png" /><figcaption>FIG. 2: (a) (Color online) Specific heat $C_h$ and relative magnetisation $\Delta M$ as a function of $t$ obtained from simulations (● for $C_h$ and ■ for $\Delta M$) and analytically (solid line for $C_h$) for $H = 10^{-3} J$. The dashed lines are a fit to logarithmic singularities given in eq. (6).</figcaption></figure>

**Figure 2** (b) (Color online) Finite-size scaling of spin correlations for $H/H_c = 0.33(\blacksquare)$, $0.66(\triangle)$ and $0.93(\bullet)$, on log-log scale. Lines have slopes 3 and 1.

with the introduction of the non-local algorithm the numerical results very clearly confirm that the transition is in fact a Kasteleyn transition similar to that observed for dimers on the brick lattice[2].

For larger fields the transition is rounded, this is because the topological constraints are no longer rigorously imposed and the system can support a finite concentration of topological defects. Above a characteristic scale fixed by their separation the system will become regularly paramagnetic and this is just what is observed in the rounding of the transition.

Despite the fact that the constraints induce long range correlations the ferromagnetic susceptibility is not divergent away from the Kasteleyn transition. The coupling between spin and configurational space for dipolar interactions leads to both direction and sublattice dependence for the amplitude of the spin spin correlation function. This gives a finite susceptibility at zero wave vector, $\vec{Q} = 0$, but very strong variation at finite $\vec{Q}$, leading to the so called pinch points on the Brillouin zone boundaries [23]. Although logarithmically divergent in 2D the intensity of the pinch points in three dimensions remains finite in zero field. Applying the appropriate magnetic field in the 2D problem breaks the square symmetry of zero field, introducing a length scale related to the symmetry breaking. As a result, the pinch points move in towards the Brillouin zone centre. Their arrival at the zone centre corresponds to the Kasteleyn transition. In 3D the pinch points remain fixed in reciprocal space but their intensity evolves with applied field, with the point at the Brillouin zone centre developing a logarithmic divergence at the transition. This evolution can be seen directly in real space through spin-spin correlation functions $C(\vec{r}) = \langle \vec{S}_i(\vec{r}) \cdot \vec{S}_i(0) \rangle$. In Figure 2(b) we show finite size scaling plots of $C(z = L/2)$ for three dimensional systems of size $L$ for magnetic field along the $z$ axis, from which we can extract the power law behaviour of $C(\vec{r})$ in the thermodynamic limit. For small field we find a $1/r^3$ behaviour consistently with the dipolar correlations leading to the pinch points. As one approaches the transition we observe a crossover to $1/r$ dependence. This is the regime where spin fluctuations along $z$ are due to multiple encounters with a single string. Here, mapping either to the quantum problem or to a random walk problem lead to the $1/r$ behaviour observed in the figure[12, 19].

The loop algorithm allows particularly accurate studies of finite size effects close to the transition where the density of strings is small. In the ordered regime, the free energy cost for the excitation of a string (5), $\delta G$ is dominated by the energy $E_{string} = 2HL/\sqrt{3}$ hence the Boltzmann weight for string excitation decreases exponentially with system size: $w_{loop} \propto \exp(-\beta E_{string})$ [21]. As the magnetisation difference $x$ is proportional to the number of strings one can observe this exponential finite size scaling by making an accurate measure of $x$ in the transition region. In Fig. (3(a)) we show $x$ on a logarithmic scale, as a function of $T_K/T$, for different system sizes. The data, averaged over up to $10^7$ configurations fit the exponential scaling remarkably well down to $x \sim 10^{-7}$. This analysis shows convincingly that the mechanism for the transition is indeed the generation of extensive excitations developed in this section.
<figure><img src="image_4.png" /><figcaption>FIG. 3: (Color online) Relative magnetisation $\Delta M$ as a function of $H/H_K$ obtained numerically (●) for different sizes of the system: $L = 20, 40, 100, 200$. The straight lines are a fit to the Boltzmann factor $e^{(L \ln 2) H/H_K}$.</figcaption></figure>

**Figure 3** (b) *M* vs *B* obtained analytically (solid lines) and experimentally (dots) by neutron scattering on the compound Ho₂Ti₂O₇ at *T* = 1.2*K*. For fitting details, see main text.

### III. DISCUSSION

As the experimental systems fall out of equilibrium at low temperature, it is difficult to compare theory with experiment very close to the Kasteleyn transition. But at reasonably low temperatures, where the concentration of topological defects is low but sufficiently high to maintain ergodicity, it is possible to observe a rounded transition and to compare data with our theory. As an example we show, in Fig. 3(b) one of the [100] field scans performed at fixed low temperature in Holmium Titanate [13]. The data is for $T = 1.2$K. This temperature is about 2/3 of the effective nearest neighbour exchange constant $J_{eff}S^2/k_B \sim 1.8$K estimated for Ho₂Ti₂O₇ [24]. The magnetisation flattens off abruptly at the saturation value, $M_{exp} = 10\mu_B/\sqrt{3}$ in a way reminiscent of the Kasteleyn transition but is indeed rounded near saturation. Also shown in the figure is our data calculated analytically from the Bethe lattice approximation for the pyrochlore system with the same ratio of $J/k_BT$. The scale of the magnetic field (that is changing from $|\vec{H}|$ in reduced units to *B* in Tesla) was estimated after demagnetisation effects were taken into account on the experimental data. A final scale factor of 1.6 was required to make the comparison. The comparison is qualitatively good, suggesting that the data of Ref. [13] do provide evidence for a rounded Kasteleyn transition in Ho₂Ti₂O₇, which could be more closely approached in future experiments.

We have enjoyed many fruitful discussions with S.T. Bramwell and T. Fennell as well as with C. Castel-novo, M.J. Harris and R. Melko. PCWH acknowledges financial support from the European Science Foundation PESC/RNP/HFM, from the Royal Society and from the London Centre of Nanotechnology. LJ acknowledges sup-port from European Science Foundation PESC/RNP/HFM and from ANR grant 05-BLAN-0105. Both thank the Rudolph Peierls Centre for Theoretical Physics, University of Oxford, where this work was begun.

### References

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

[2] S. M. Bhattacharjee, J. F. Nagle, D. A. Huse, and M. E. Fisher, J. Stat. Phys. **32**, 361 (1983).

[3] S. M. Bhattacharjee and J. J. Rajasekaren, Phys. Rev. A **44**, 6202 (1991).

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

[5] G. I. Watson, J. Stat. Phys. **94**, 1045 (1999).

[6] M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, and K. W. Godfrey, Phys. Rev. Lett. **79**, 2554 (1997).

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

[8] S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. **93**, 167204 (2004).
[9] C. L. Henley, Phys. Rev. B **71**, 014424 (2005).

[10] T. Fennell, S. T. Bramwell, D. F. McMorrow, P. Manuel, and A. Wildes, Nature Physics **3**, 566 (2007).

[11] M. J. Harris, S. T. Bramwell, P. C. W. Holdsworth, and J. D. M. Champion, Phys. Rev. Lett. **81**, 4496 (1998).

[12] L. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R. Moessner, Phys. Rev. Lett **100**, 067207 (2008).

[13] T. Fennell, O. A. Pentrenko, B. Fak, J. S. Garnder, S. T. Bramwell, and B. Ouladdiat, Phys. Rev. B **72**, 224411 (2005).

[14] H. Fukazawa, R. G. Melko, R. Higashinaka, Y. Maeno, and M. P. J. Gingras, Phys. Rev. B **65**, 054410 (2002).

[15] J. D. Bernal and R. H. Fowler, J. Chem. Phys. **1**, 515 (1933).

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

[17] H. E. Stanley, *Introduction to Phase transitions and Critical Phenomena* (OUP, 1971).

[18] C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature **451**, 42 (2008).

[19] S. Powell and J. T. Chalker, cond-mat/0803.4204 (2008).

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

[21] S. V. Isakov, K. S. Raman, R. Moessner, and S. L. Sondhi, Phys. Rev. B **70**, 104418 (2004).

[22] L. Jaubert, P. C. W. Holdsworth, and R. Moessner, Unpublished (2008).

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

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