Mean-field theory for confinement transitions and magnetization plateaux in spin ice

TL;DR

This paper develops a mean-field theory to analyze phase transitions in classical spin ice under magnetic fields, revealing a phase diagram with Coulomb phases and magnetization plateaux, highlighting the generic lobe structure in Coulomb spin liquids.

Contribution

It introduces a novel mean-field approach to study confinement and topological effects in spin ice, connecting phase diagram features to Bose-Hubbard-like models.

Findings

Identification of a phase diagram with Coulomb and magnetization plateau phases

Discovery of a lobe structure in the phase diagram similar to Bose-Hubbard models

Demonstration of the generic nature of this structure in Coulomb spin liquids

Abstract

We study phase transitions in classical spin ice at nonzero magnetization, by introducing a mean-field theory designed to capture the interplay between confinement and topological constraints. The method is applied to a model of spin ice in an applied magnetic field along the [100] crystallographic direction and yields a phase diagram containing the Coulomb phase as well as a set of magnetization plateaux. We argue that the lobe structure of the phase diagram, strongly reminiscent of the Bose-Hubbard model, is generic to Coulomb spin liquids.

Equations

$${\mathcal{H}}_{\text{nn}}=-J\sum _{⟨ij⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j=\frac{J}{3}\sum _{⟨ij⟩}{\sigma }_i{\sigma }_j\text{.}$$

$${\mathcal{H}}_{\text{nn}}=\frac{2}{3}J\sum _t\left(Q_t^2-1\right)\text{,}$$

$$Q_t=\frac{1}{2}{ϵ}_t\sum _{i\in t}{\sigma }_i$$

$$U_{\text{m}}({\mathbf{r}}_{t,t^{\prime }})=-T\ln {\mathcal{Z}}_{\text{m}}({\mathbf{r}}_{t,t^{\prime }})$$

$${\mathcal{Z}}_{\text{m}}({\mathbf{r}}_{t,t^{\prime }})=\sum _{\{{\sigma }_i\}\in {\mathfrak{C}}_{t,t^{\prime }}}{\mathrm{e}}^{-V/T}\text{,}$$

$$\mathcal{Z}[{\zeta }_t]=\sum _{\{{\sigma }_i\}}{\mathrm{e}}^{-V/T}\prod _t{\zeta }_t^{|Q_t|}\text{,}$$

$$G_{\text{m}}({\mathbf{r}}_{t,t^{\prime }})\equiv \frac{{\mathcal{Z}}_{\text{m}}({\mathbf{r}}_{t,t^{\prime }})}{\mathcal{Z}}=\frac{1}{2}{\frac{\partial }{\partial {\zeta }_t}\frac{\partial }{\partial {\zeta }_{t^{\prime }}}\ln \mathcal{Z}[{\zeta }_t]|}_{{\zeta }_{\cdot }=0}\text{,}$$

$$\mathcal{H}={\mathcal{H}}_{\text{nn}}+{\mathcal{H}}_u+{\mathcal{H}}_{\mathbf{h}}\text{,}$$

$${\mathcal{H}}_{\mathbf{h}}=-\mathbf{h}\cdot \mathbf{M}=-\mathbf{h}\cdot \sum _i{\mathbf{S}}_i\text{,}$$

$${\mathcal{H}}_{\mathbf{h}}=-\frac{h}{\sqrt{3}}\sum _i{\eta }_i{\sigma }_i\text{,}$$

$$M_z=\sum _i{\mathbf{u}}_z\cdot {\mathbf{S}}_i=\frac{1}{\sqrt{3}}\sum _i{\eta }_i{\sigma }_i\text{.}$$

$${\mathcal{H}}_u=+\frac{u}{2}\sum _{⟨ij⟩\in \mathfrak{R}}(3{\mathbf{S}}_i\cdot {\mathbf{S}}_j-1)=-\frac{u}{2}\sum _{⟨ij⟩\in \mathfrak{R}}({\sigma }_i{\sigma }_j+1)\text{,}$$

$$\mathcal{Z}=\text{Tr}\prod _{z=L_z}^1{\mathcal{T}}_z=\text{Tr}\left({\mathcal{T}}_{L_z}{\mathcal{T}}_{L_z-1}\cdots {\mathcal{T}}_2{\mathcal{T}}_1\right)\text{,}$$

$${\mathcal{T}}_z=\prod _{c{c}^{\prime }\in {\mathcal{P}}_z}{\mathcal{T}}_{c{c}^{\prime }}\text{.}$$

$$\mathcal{Z}=\text{Tr}{\left[{\mathcal{T}}^{({\Lambda }_z)}\right]}^{L_z/{\Lambda }_z}\text{,}\text{where}{\mathcal{T}}^{({\Lambda }_z)}=\prod _{z={\Lambda }_z}^1{\mathcal{T}}_z\text{.}$$

$${\mathcal{Z}}_c^{\text{mf}}={\text{Tr}}_c\prod _{z=L_z}^1{\mathcal{T}}_{c,z}^{\text{mf}}\text{,}$$

$${\mathcal{T}}_{c,z}^{\text{mf}}={⟨v_{{\mathcal{C}}_{c,z},z}|}_{{\mathcal{C}}_{c,z}}{\mathcal{T}}_{c,{\mathcal{C}}_{c,z}}{|v_{{\mathcal{C}}_{c,z},z-1}⟩}_{{\mathcal{C}}_{c,z}}\text{,}$$

$$\prod _{z^{\prime }=z}^{z-{\Lambda }_z+1}{\mathcal{T}}_{c,z}^{\text{mf}}\text{.}$$

$${\mathcal{T}}_{c,z}^{\text{mf}}={⟨v_{{\mathcal{C}}_{c,z}}|}_{{\mathcal{C}}_{c,z}}{\mathcal{T}}_{c,{\mathcal{C}}_{c,z}}{|v_{{\mathcal{C}}_{c,z}}⟩}_{{\mathcal{C}}_{c,z}}\text{,}$$

$${\mathcal{T}}_{c{c}^{\prime }}=\begin{array}{rl}& {\mathrm{e}}^{\frac{2}{\sqrt{3}}h/T}{|\uparrow ⟩}_c{|\uparrow ⟩}_{c^{\prime }}{⟨\uparrow |}_c{⟨\uparrow |}_{c^{\prime }}+{\mathrm{e}}^{-\frac{2}{\sqrt{3}}h/T}{|\downarrow ⟩}_c{|\downarrow ⟩}_{c^{\prime }}{⟨\downarrow |}_c{⟨\downarrow |}_{c^{\prime }}\\ & +{\mathrm{e}}^{2u/T}{|\uparrow ⟩}_c{|\downarrow ⟩}_{c^{\prime }}{⟨\uparrow |}_c{⟨\downarrow |}_{c^{\prime }}+{\mathrm{e}}^{2u/T}{|\downarrow ⟩}_c{|\uparrow ⟩}_{c^{\prime }}{⟨\downarrow |}_c{⟨\uparrow |}_{c^{\prime }}\\ & +{|\uparrow ⟩}_c{|\downarrow ⟩}_{c^{\prime }}{⟨\downarrow |}_c{⟨\uparrow |}_{c^{\prime }}+{|\downarrow ⟩}_c{|\uparrow ⟩}_{c^{\prime }}{⟨\uparrow |}_c{⟨\downarrow |}_{c^{\prime }}\text{.}\end{array}$$

$${\mathcal{T}}_{c,z}^{\text{mf}}=\begin{array}{rl}& \left[\frac{1}{2}(1+m_{c^{\prime }}^{\text{mf}}){\mathrm{e}}^{\frac{2}{\sqrt{3}}h/T}+\frac{1}{2}(1-m_{c^{\prime }}^{\text{mf}}){\mathrm{e}}^{2u/T}\right]{|\uparrow ⟩}_c{⟨\uparrow |}_c\\ & +\left[\frac{1}{2}(1-m_{c^{\prime }}^{\text{mf}}){\mathrm{e}}^{-\frac{2}{\sqrt{3}}h/T}+\frac{1}{2}(1+m_{c^{\prime }}^{\text{mf}}){\mathrm{e}}^{2u/T}\right]{|\downarrow ⟩}_c{⟨\downarrow |}_c\\ & +{\zeta }_{c^{\prime }}^{\text{mf}}{|\uparrow ⟩}_c{⟨\downarrow |}_c+{\left({\zeta }_{c^{\prime }}^{\text{mf}}\right)}^{\ast }{|\downarrow ⟩}_c{⟨\uparrow |}_c\text{,}\end{array}$$

$$m_{c^{\prime }}^{\text{mf}}$$

$${\zeta }_{c^{\prime }}^{\text{mf}}$$

$$\Theta (|v_1⟩,|v_2⟩)={⟨v_1|}_1{⟨v_2|}_2{\mathcal{T}}_{12}{|v_1⟩}_1{|v_2⟩}_2$$

$$|v_c⟩=|\uparrow ⟩\cos \frac{{\theta }_c}{2}+|\downarrow ⟩{\mathrm{e}}^{\mathrm{i}{\varphi }_c}\sin \frac{{\theta }_c}{2}\text{,}$$

$$\begin{gathered} \Theta (|v_1⟩,|v_2⟩)=\frac{1}{2}(1+\cos {\theta }_1\cos {\theta }_2)\cosh \frac{2h}{\sqrt{3}T}+\frac{1}{2}(\cos {\theta }_1+\cos {\theta }_2)\sinh \frac{2h}{\sqrt{3}T} \\ +\frac{1}{2}(1-\cos {\theta }_1\cos {\theta }_2){\mathrm{e}}^{2u/T}+\frac{1}{2}\sin {\theta }_1\sin {\theta }_2\cos ({\varphi }_1-{\varphi }_2)\text{.} \end{gathered}$$

$$\frac{M_z}{M_{\text{sat}}}=\frac{1}{2}\left(m_1^{\text{mf}}+m_2^{\text{mf}}\right)=\frac{1}{2}(\cos {\theta }_1+\cos {\theta }_2)\text{.}$$

$$\frac{M_z}{M_{\text{sat}}}=\{\begin{array}{ll}\mathrm{sgn}h& \text{if T≤TK}\\ \frac{\sinh \frac{2h}{\sqrt{3}T}}{2-\cosh \frac{2h}{\sqrt{3}T}}& \text{if T>TK,}\end{array}$$

Full text

I Introduction

Classical spin liquids Balents2010 , such as the Coulomb phase Henley2010 in spin ice Bramwell2001 and related systems, are examples of phases whose behavior is not captured by the standard Landau picture of broken symmetries Landau1999 . Their two defining characteristics are fractionalization, the emergence of excitations not constructed from finite combinations of the elementary degrees of freedom, and topological order, the presence of structure that can only be discerned by observing the system globally Castelnovo2012 .

In the case of the classical spin ices Bramwell2001 , a family of magnetic pyrochlore oxides, the consequences of these properties have been of sustained interest, from both theoretical and experimental perspectives. The spins in these materials carry large magnetic moments, and their fractional excitations take the form of magnetic monopoles Castelnovo2008 , acting as sources for the physical magnetic field. Furthermore, transitions between phases where these excitations are confined and deconfined have particularly interesting properties that cannot be captured by the Landau–Ginzburg theory of phase transitions MonopoleScalingPRL ; MonopoleScalingPRB .

Spin ice also has the unusual feature that the magnetization is a topologically constrained quantity, and so topological order can be probed directly Castelnovo2012 . In particular, within the low-energy configuration space relevant at low temperatures, any local dynamics conserves the uniform magnetization Henley2010 . This fact is known to have interesting consequences for critical properties at certain confinement transitions, such as the Kasteleyn transition in an applied magnetic field Jaubert2008 ; SpinIceCQ .

In this work, we investigate the interplay between these two aspects of the Coulomb spin liquid by studying confinement phase transitions in spin ice across the full range of magnetization. As noted by Henley Henley2010 , one can distinguish three categories of phase transition from the Coulomb phase: The first, where the magnetization is zero across the transiton, was studied systematically in Ref. SpinIceHiggs . A second, where the ordered state has saturated polarization, was considered in Refs. Jaubert2008 ; SpinIceCQ . Here we consider the phase structure and transitions more generally, including the third case, where the magnetization is nonzero but unsaturated, and may or may not change across the transition.

To do so, we introduce a mean-field theory (MFT) designed to capture the distinction between confined and deconfined phases. Standard mean-field approaches, based on the Landau picture of ordered phases, are not well suited to describing confinement transitions. Instead, we take inspiration from the well-known MFT for the Bose–Hubbard model Fisher1989 , making use of a mapping from spin ice to a quantum model of bosons SpinIceCQ . We apply the method to a simplified model of classical spin ice, in order to illustrate the general approach and its physical interpretation. We also briefly discuss extensions to the method that could be used to treat more physically realistic perturbations.

In addition, we present arguments for the generality of our results, beyond mean-field theory and our particular choice of Hamiltonian. We argue that the phase diagram generically consists of a set of confined phases at low temperature in which the magnetization is, to a very good approximation, fixed. This plateau structure is strongly reminiscent of the lobes present in the phase diagram of the Bose–Hubbard model Fisher1989 ; Sachdev2011 , a connection that helps clarify the general phase structure.

While the focus here is on classical spin ice in a $[100]$ field, the distinction between transitions in different flux sectors is more general. For example, spin ice in a field along the $[111]$ direction also exhibits magnetization plateaux Matsuhira2002 ; Moessner2003 ; Isakov2004 , while it has been suggested that the zero-field ground state of quantum spin ice may have nonzero magnetization Shannon2012 .

In Section II, we introduce the model to be used and briefly outline the general properties of the Coulomb phase and confinement transitions. The mean-field method is then described and applied to the model in Section III. We present general arguments for the phase structure and the nature of the phase transitions in Section IV, before concluding in Section V.

II Model and basic properties

Our presentation will be based on the case of classical spin ice in a magnetic field applied along the $[100]$ crystallographic direction. We first introduce the nearest-neighbor model of spin ice (NNSI) that describes the physics of the Coulomb spin liquid, before discussing perturbations, such as an applied field, and the resulting ordering transitions.

II.1 Nearest-neighbor spin ice

The spin ices Bramwell2001 are a family of frustrated magnetic materials with moments on the sites $i$ of a pyrochlore lattice, a network of corner-sharing tetrahedra. Prominent examples are the “classical” spin ices, such as $\mathrm{Dy}_{2}\mathrm{Ti}_{2}\mathrm{O}_{7}$ and $\mathrm{Ho}_{2}\mathrm{Ti}_{2}\mathrm{O}_{7}$, which are well described in terms of classical spins $\boldsymbol{S}_{i}$. Each spin is subject to strong easy-axis anisotropy along the local $\langle 111\rangle$ direction $\hat{\boldsymbol{n}}_{i}$ which joins the centres of the two tetrahedra sharing site $i$, as shown in Fig. 1.

The anisotropy results from a large (of order $100\text{\,}\mathrm{K}$) crystal-field splitting between states with maximum projection along $\hat{\boldsymbol{n}}_{i}$ and all other states in the single-moment Hilbert space. This large gap ensures that tunneling processes between the low-energy states are strongly suppressed, rendering the moments effectively classical Bramwell2001 . For the low temperatures of interest, it is therefore sufficient to treat the spins as effectively constrained to $\boldsymbol{S}_{i}=\sigma_{i}\hat{\boldsymbol{n}}_{i}$, where $\sigma=\pm 1$ is a classical Ising variable.

Each tetrahedron $t$ in the pyrochlore lattice can be labeled according to its orientation $\epsilon_{t}=\pm 1$, and the lattice structure is such that all neighbors of $t$ have orientation $-\epsilon_{t}$. The fixed unit vector $\hat{\boldsymbol{n}}_{i}$ for each pyrochlore site $i$, directed along the line joining the centers of its two tetrahedra, can be chosen to point towards the tetrahedron with $\epsilon_{t}=-1$, and so $\epsilon_{t}\sigma_{i}$ indicates whether $\boldsymbol{S}_{i}$ points out of ($+1$) or into ($-1$) tetrahedron $t$. With this convention, $\hat{\boldsymbol{n}}_{i}\cdot\hat{\boldsymbol{n}}_{j}=-\frac{1}{3}$ and hence $\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{j}=-\frac{1}{3}\sigma_{i}\sigma_{j}$, for all nearest neighbors $i$ and $j$.

The minimal Hamiltonian that captures the spin-liquid physics of the classical spin ices Castelnovo2012 contains only interactions between nearest-neighbor pairs $\langle ij\rangle$,

[TABLE]

The coupling between moments $\boldsymbol{S}_{i}$ is ferromagnetic ($J>0$), incorporating the net effect of the dipolar interactions between nearest neighbors and (antiferromagnetic) exchange Bramwell2001 . The effective interactions between the Ising variables $\sigma_{i}$ are therefore antiferromagnetic, and hence frustrated. Further-neighbor interactions, both dipolar and exchange, are significant in spin ice, but can be treated as a relatively small perturbation to $\mathcal{H}_{\text{nn}}$ Castelnovo2012 ; we will return to these in Section II.3.

Because two sites are nearest neighbors if and only if they share a tetrahedron, the interactions can be rewritten as

[TABLE]

where $\sum_{t}$ includes all tetrahedra (of both orientations) and

[TABLE]

is referred to as the (fictional) “charge” on tetrahedron $t$ in terms of a sum over its sites $i\in t$. The energy is therefore minimized by configurations with $Q_{t}=0$, which are those where two spins point into and two out of each tetrahedron. Configurations that satisfy this condition are said to obey the ice rule Bramwell2001 ; Castelnovo2012 ; all configurations that do so at every tetrahedron are degenerate minimal-energy states of $\mathcal{H}_{\text{nn}}$.

Tetrahedron configurations with three spins in and one out, or vice versa, have energy $\Delta=+\frac{2}{3}J$ above the ice-rule configurations; they have charge $Q_{t}=\pm 1$ and are hence referred to as “monopoles”. (Note that $\sum_{t}Q_{t}=0$, and so all configurations are globally charge-neutral.) For $\mathrm{Dy}_{2}\mathrm{Ti}_{2}\mathrm{O}_{7}$, the effective nearest-neighbor coupling is $\frac{1}{3}J=$1.1\text{,}\mathrm{K} [Bramwell2001 ](#bib.bib5), corresponding to a monopole cost of $\Delta=$2.2\text{\,}\mathrm{K}. The corresponding Boltzmann weight $\zeta=\mathrm{e}^{-\Delta/T}$ is called the monopole fugacity. While the significance of the term “monopole” is here only that they carry effective “charge” $Q_{t}$, it has been argued that monopole excitations in the classical spin ice materials in fact carry true magnetic charge Castelnovo2008 .

For temperature $T$ significantly below $\Delta$, the system is effectively constrained by the ice rule, and the ensemble of ice-rule configurations therefore determines its behavior Bramwell2001 . The number of such configurations grows exponentially with the total number of spins $N$, and so the entropy density remains nonzero even well below $\Delta$, the temperature scale at which an unfrustrated system would order. In spite of the absence of order, the ensemble exhibits correlations that decay only algebraically with distance. The resulting correlated paramagnet is referred to as a Coulomb phase Henley2010 , and is an example of a classical spin liquid Balents2010 .

The long-wavelength properties of the Coulomb phase can be captured by coarse-graining the vector spins $\boldsymbol{S}$ to give a continuum field $\boldsymbol{\mathfrak{B}}(\boldsymbol{r})$ referred to as the effective “magnetic field” Henley2010 . The ice-rule constraint applied to $\boldsymbol{S}$ implies that $\boldsymbol{\mathfrak{B}}$ is divergenceless, while monopoles act as sources or sinks, according to an effective Gauss law. The resulting predictions for spin–spin correlation functions are in good agreement with neutron-scattering experiments in classical spin ice Henley2010 .

The ice-rule configurations also display interesting topological properties Castelnovo2012 . For our purposes, the most important is that any pair of ice-rule configurations is connected by flipping all spins along one or more loops. For a system with periodic boundary conditions (which we will assume), such a loop can only change the total magnetization $\boldsymbol{M}=\sum_{i}\boldsymbol{S}_{i}$ if it has nontrivial winding number. For other systems that host Coulomb phases Henley2010 , such as classical dimer models, an analogous quantity, called the “flux”, can be defined, but does not necessarily correspond to a quantity that is so directly accessible experimentally.

II.2 Confinement transitions

A second important property of the Coulomb phase is that monopoles, defects in the ice rule, are deconfined Castelnovo2008 . To make this statement precise, imagine imposing a pair of monopoles with opposite charges $\pm 1$ at tetrahedra $t$ and $t^{\prime}$ in a background that otherwise obeys the ice rule. Let $\mathfrak{C}_{t,t^{\prime}}$ denote the set of spin configurations $\{\sigma_{i}\}$ compatible with this charge configuration. Their number $\mathcal{Z}_{\text{m}}(\boldsymbol{r}_{t,t^{\prime}})=\lvert\mathfrak{C}_{t,t^{\prime}}\rvert$ is a function of the separation $\boldsymbol{r}_{t,t^{\prime}}$ that decreases with $\lvert\boldsymbol{r}_{t,t^{\prime}}\rvert$ but reaches a finite limit as $\lvert\boldsymbol{r}_{t,t^{\prime}}\rvert\rightarrow\infty$. The result is an effective entropic interaction

[TABLE]

between the charges that is bounded, and hence allows them to be separated to infinity at finite free-energy cost. In fact, since monopoles are charges in the continuum field $\boldsymbol{\mathfrak{B}}$, coarse-graining predicts that the interaction obeys the Coulomb law, $U_{\text{m}}(\boldsymbol{r})\sim\text{const}-\lvert\boldsymbol{r}\rvert^{-1}$, for large $\lvert\boldsymbol{r}\rvert$ Henley2010 .

More generally, consider the partition function

[TABLE]

where $V$ is a perturbation that splits the degeneracy of the ice-rule configurations. This will tend to suppress fluctuations, and can eventually, as the temperature is reduced compared to the scale of the perturbation, drive the system out of the Coulomb phase. A transition into a phase where $U_{\text{m}}(\boldsymbol{r}_{t,t^{\prime}})$ increases without limit as $\lvert\boldsymbol{r}_{t,t^{\prime}}\rvert\rightarrow\infty$, and so $\mathcal{Z}_{\text{m}}(\boldsymbol{r}_{t,t^{\prime}})\rightarrow 0$, is referred to as a confinement transition MonopoleScalingPRL ; MonopoleScalingPRB . (In terms of the effective continuum field $\boldsymbol{\mathfrak{B}}$, Gauss’ law implies that an imposed pair of charges must be joined by a fixed quantity of flux. Linear confinement, where $U_{\text{m}}(\boldsymbol{r})\sim\lvert\boldsymbol{r}\rvert$ for large $\rvert\boldsymbol{r}\rvert$, occurs when the flux forms a narrow tube connecting the monopoles, with finite tension.)

While confinement may occur simultaneously with appearance of conventional symmetry-breaking order Sreejith2014 , there is no local order parameter that provides a signature of confinement. The confinement criterion, that $\mathcal{Z}_{\text{m}}(\boldsymbol{r}_{t,t^{\prime}})$ approaches zero as $\lvert\boldsymbol{r}_{t,t^{\prime}}\rvert\rightarrow\infty$ in a confined phase and conversely has a finite limit in a deconfined phase, instead resembles the characterization of a conventional symmetry-broken phase in terms of long-range order. This connection can be pursued further by considering the partition function with an unconstrained charge distribution and a position-dependent fugacity $\zeta_{t}$,

[TABLE]

in terms of which FootnoteTotalChargeZero

[TABLE]

where $\mathcal{Z}=\left.\mathcal{Z}[\zeta_{t}]\right\rvert_{\zeta_{\cdot}=0}$ is the partition function restricted to ice-rule configurations. The relationship between the monopole distribution function $G_{\text{m}}$ and the fugacity $\zeta$ expressed by Eq. (7) is exactly analogous to the fluctuation–dissipation theorem relating, for example, a spin–spin correlation function and an applied field. This connection motivates a mean-field theory, analogous to standard Weiss theory Yeomans1992 , that describes the confinement transition in terms of an effective self-consistent fugacity.

It should be noted that the confinement distinction is only sharp in the limit $\zeta=0$, because any nonzero density of monopoles screens the effective interaction $U_{+1,-1}$ at large separation. There are nonetheless implications for the critical properties away from this limit MonopoleScalingPRL ; MonopoleScalingPRB . If the phase transition also involves a different type of order (such as spontaneous symmetry breaking), then it may survive for $\zeta>0$, but with a different universality class Sreejith2014 .

II.3 Perturbations to NNSI

Within the nearest-neighbor model $\mathcal{H}_{\text{nn}}$, all configurations that obey the ice rules are degenerate and the system exhibits a deconfining Coulomb phase. To drive a confinement transition, one must add perturbations that split the degeneracy of the ice-rule configurations. Here we consider a Hamiltonian of the form

[TABLE]

where $\mathcal{H}_{\boldsymbol{h}}$ is the coupling to an applied magnetic field and $\mathcal{H}_{u}$ contains additional interactions described below.

Coupling to a magnetic field $\boldsymbol{h}$ is described by a Zeeman term

[TABLE]

where a factor dependent on the size of the magnetic moments has been included in the definition of $\boldsymbol{h}$. The case that we consider here has $\boldsymbol{h}=h\boldsymbol{u}_{z}$ where $\boldsymbol{u}_{z}$ is a unit vector along $[100]$, and so $\boldsymbol{h}\cdot\boldsymbol{S}_{i}=h\boldsymbol{u}_{z}\cdot\sigma_{i}\hat{\boldsymbol{n}}_{i}=\pm\frac{1}{\sqrt{3}}h\sigma_{i}$. The coupling can therefore be written as

[TABLE]

where $\eta_{i}=\sqrt{3}\hat{\boldsymbol{n}}_{i}\cdot\boldsymbol{u}_{z}=\pm 1$. Using the choice described in Section II.1, where $\hat{\boldsymbol{n}}_{i}$ points towards the tetrahedron with $\epsilon_{t}=-1$, $\eta_{i}$ alternates in sign on successive $(100)$ planes of spins. The component of the magnetization along the field can likewise be expressed as

[TABLE]

The perturbation $\mathcal{H}_{\boldsymbol{h}}$ alone can drive a confinement transition, referred to as a Kasteleyn transition Jaubert2008 . For $\Delta\gg h\gg T$, all spins align with the field (to the extent allowed by the easy-axis anisotropy), an arrangement that satisfies the ice rules. Any excitation within the ice-rule states involves flipping a loop of spins, but starting from the fully polarized state, the only available loops span the system in the $[100]$ direction. Such excitations therefore have energy cost proportional to the linear system size $L_{z}$ and are suppressed in the thermodynamic limit.

The result is a strictly saturated magnetization $M_{z}=M_{\text{sat}}\equiv\frac{1}{\sqrt{3}}N$ even for nonzero $T/h$ (but in the limit $T/\Delta\rightarrow 0$). Above a temperature $T_{\text{K}}$, however, the energy cost is outweighed by the gain in entropy, also $\propto L_{z}$, associated with the multitude of available paths for a spanning loop, and a Kasteleyn transition occurs from the saturated paramagnet to the Coulomb phase. The transition is possible only because the ice rule restricts to loop-like excitations; away from the limit $T/\Delta=0$, the transition is replaced by a crossover Jaubert2008 ; MonopoleScalingPRB .

To study the interplay between the topological constraints on the magnetization, inherent in the ice-rule configurations, and conventional ordering transitions, we also consider an additional interaction $\mathcal{H}_{u}$ that induces spontaneous symmetry breaking. A natural choice would be the further-neighbor coupling present in the spin-ice compounds, due to a combination of further-neighbor exchange Henelius2016 and long-range dipolar interactions. The latter have been shown using Monte Carlo Melko2001 to lead to a phase transition into the ordered configuration shown in Fig. 2, which we will call the Melko–den Hertog–Gingras (MDG) state FootnoteNomenclature .

The actual perturbation $\mathcal{H}_{u}$ that we will use here is one that is simpler to treat within our mean-field theory but that is expected to lead to an ordered phase of the same type. It is given by

[TABLE]

where $\mathfrak{R}$ is the set of nearest-neighbor bonds highlighted in Fig. 2, which form a set of left-handed screw chains. The effect of $u>0$ is to reduce the strength of the interactions on these bonds and to favor the MDG configuration shown, or the one with all spins inverted. Fig. 3 shows the energy for each of the six configurations of a single tetrahedron that obey the ice rule, including both the applied field and the perturbation $\mathcal{H}_{u}$.

Unlike dipolar interactions, the perturbation $\mathcal{H}_{u}$ reduces the symmetry of the lattice by picking an axis and a chirality. The 12-fold degeneracy of the MDG state is therefore reduced by a factor of $3\times 2$, but the symmetry under inversion of all spins remains, and so an ordering transition into this state spontaneously breaks an Ising symmetry. This transition (at $h=0$) is in fact a type considered in Ref. SpinIceHiggs , where it was argued to be most likely of first order.

As we will show, this simplified model has a phase diagram that illustrates general features emerging from the interplay between deconfinement, ordering, and topological constraints. In Section V, we will briefly discuss the prospects for including more realistic perturbations in the mean-field theory.

II.4 Mapping to extended Bose–Hubbard model

Some aspects of the phase structure can be elucidated using a mapping between classical spin ice and an effective quantum model of hard-core bosons on a lattice Jaubert2008 ; SpinIceCQ . This proceeds by identifying the $[100]$ direction, along which the external field is applied, with the imaginary-time axis in a path-integral representation of the quantum problem. The thermodynamic limit of the original classical system therefore corresponds to the zero-temperature limit in the quantum problem.

This mapping has previously been studied in detail for the Hamiltonian $\mathcal{H}_{\text{nn}}+\mathcal{H}_{\boldsymbol{h}}$ SpinIceCQ . The fully polarized state (along the $+z$ direction, say) is identified with the vacuum; flipped spins relative to this configuration form system-spanning loops, which are viewed as the world lines of bosons. The result is an effective model of hard-core bosons on a lattice with density $\rho=\frac{1}{2}(1-M_{z}/M_{\text{sat}})$ and hence chemical potential $\mu\sim-h$. In general, one expects all other local terms compatible with symmetry, such as hopping between nearby sites and further-neighbor interactions, with parameters related in nontrivial ways to those of the spin model SpinIceCQ .

As $\mu$ is increased, the system passes from the vacuum through a quantum phase transition Sachdev2011 to a superfluid, which corresponds to the Coulomb phase of spin ice. (Off-diagonal long-range order in the superfluid can be associated with deconfinement in the Coulomb phase SpinIceCQ .) For sufficiently large positive $\mu$, the system is driven into a fully occupied Mott insulator, with $\rho=1$, corresponding to a saturated paramagnet with $M_{z}=-M_{\text{sat}}$ in an inverted applied field. Critical properties at the phase transitions in spin ice are, mutatis mutandis, equivalent to those in the quantum model SpinIceCQ .

The presence of interactions between bosons on different lattice sites also makes possible Mott insulators with fractional filling. These can be identified with confining phases of spin ice (because they lack superfluid order) in which the magnetization is fixed to a value less than $M_{\text{sat}}$. The MDG phase, for example, is an ordered phase with $M_{z}=0$, and corresponds to a Mott insulator at half filling.

In general, more exotic phases such as supersolids are also possible in the effective quantum model. A supersolid would correspond to a Coulomb phase with broken spatial symmetry; such phases have been proposed in spin ice BrooksBartlett2014 , but are not captured by the mean-field theory used here.

III Transfer-matrix mean-field theory

The phase structure of $\mathcal{H}$ in Eq. (8) is expected to contain: the Coulomb phase, characterized by deconfinement of monopoles; the saturated paramagnet, in which fluctuations are suppressed by the ice-rule constraints; and the MDG phase, with conventional order that breaks spin-inversion and spatial symmetries. A mean-field theory that describes all three must therefore capture the confinement–deconfinement distinction, unlike mean-field approaches based on finite clusters Reimers1991 , and also respect the spatial structure, unlike Bethe-lattice calculations Jaubert2008 .

Here we apply an approach that is analogous to the standard mean-field theory for the Bose–Hubbard model Fisher1989 , but is applied directly to the classical partition function, written in terms of a transfer matrix. (For other applications of the variational method to the transfer matrix, see Refs. Baxter ; Rujan .) We start by partitioning the lattice into chains that span the system, which are the minimal units that can obey the ice rule. An effective model is then defined for each chain, with the coupling between chains treated self-consistently.

In standard mean-field theory for a symmetry-breaking transition Yeomans1992 , coupling between lattice sites is replaced by an effective field related to the order parameter. Here, in the absence of a local order parameter, we capture the confinement–deconfinement transition by introducing an effective self-consistent fugacity $\zeta^{\text{mf}}$ for monopoles. Since, according to Eq. (7), the monopole distribution function can be interpreted as the correlation function corresponding to a monopole fugacity, a nonzero $\zeta^{\text{mf}}$ implies deconfinement.

This decoupling is closely analogous to that used for the Bose–Hubbard model, where superfluid order is captured by introducing a self-consistent source term for bosons. Following the mapping between this quantum problem and the classical statistical mechanics of spin ice, we apply the method at the level of the transfer matrix. It should be emphasized that this differs from the standard mean-field approach for classical statistical systems, where one chooses a trial distribution to minimize the free energy. (The latter method cannot be used, because a trial configuration with nonzero fugacity does not enforce the ice rule and therefore has infinite energy.)

While we treat deconfinement using an effective mean-field fugacity $\zeta^{\text{mf}}$, we set the true monopole fugacity $\zeta=\mathrm{e}^{-\Delta/T}$ to zero. Nonzero $\zeta$ could straightforwardly be included in the approximation and is analogous to standard Weiss mean-field theory in the presence of an applied magnetic field Yeomans1992 .

III.1 Decomposition of the pyrochlore lattice into chains

We choose to decompose the pyrochlore lattice into right-handed screw chains with axes parallel to the $[100]$ direction, as illustrated in Fig. 4.

For our purposes, this choice has two advantages: First, the chains span the system in the imaginary-time direction under the mapping to a bosonic model, and so the connection to the Bose–Hubbard mean-field theory is transparent. Second, the ordered states that we aim to describe have natural interpretations in terms of chain order: in the MDG state, shown in Fig. 4, the chains order in a checkerboard pattern, while in the fully polarized paramagnet, all chains are aligned. Other orderings, such as the partially polarized state observed by Lin and Kao Lin2013 , can be described in terms of ordered configurations of the same chains, while the general approach can likely be adapted for ordered states that are not captured by this decomposition of the lattice. (For example, McClarty et al. McClarty2015 consider ordered states described in terms of linear chains of spins.)

The division of the lattice into screw chains has the property that each chain $c$ has $4$ neighbors, which it meets in turn at tetrahedra in successive $(100)$ layers. The set of tetrahedra at a given layer, at vertical position $z$, amounts to a partition of the chains into pairs, which we denote $\mathcal{P}_{z}$. We also define $\mathcal{C}_{c,z}$ as the chain that meets chain $c$ in a tetrahedron at layer $z$. The pyrochlore lattice is symmetric under translations by $\Lambda_{z}=4$ layers in the $[100]$ direction SpinIceCQ , and so $\mathcal{P}_{c+\Lambda_{z}}=\mathcal{P}_{c}$ and $\mathcal{C}_{c,z+\Lambda_{z}}=\mathcal{C}_{c,z}$.

In the quantum mapping, the vertical direction is treated as imaginary time. Applying the mapping at a microscopic level would therefore lead to a time-evolution operator that couples different pairs of lattice sites at successive time steps. This “stroboscopic” Bukov2015 implementation of the hopping in the effective quantum Hamiltonian has been shown to be irrelevant for the critical properties SpinIceCQ , but must be handled with some care when implementing the mean-field theory.

III.2 General approach

The partition function $\mathcal{Z}$ can be expressed exactly in terms of a set of operators acting on degrees of freedom associated with the chains. To each chain $c$ is associated a pair of basis vectors, $\ket{\uparrow}_{c}$ and $\ket{\downarrow}_{c}$, representing the orientation, up or down (i.e., $\sigma_{i}\eta_{i}=\pm 1$), of the spin at a given vertical position $z$ along the chain. Denoting by $\mathfrak{H}_{c}$ the space spanned by these vectors, the full Hilbert space for a layer of the system is $\mathfrak{H}=\bigotimes_{c}\mathfrak{H}_{c}$. (Note that there is a single two-dimensional Hilbert space for each chain, and so the basis vectors are labeled only by their chain $c$, and not by their position $z$ along the chain.)

The partition function can then be written as

[TABLE]

where $L_{z}$ is the system size (number of layers) in the $z$ direction and

[TABLE]

is the transfer operator for the layer of tetrahedra at vertical position $z$. (We assume periodic boundary conditions in the $z$ direction.) The tetrahedron transfer operator $\mathcal{T}_{cc^{\prime}}$, acting on $\mathfrak{H}_{c}\otimes\mathfrak{H}_{c^{\prime}}$, accounts for the possible configurations of the chains $c$ and $c^{\prime}$ on entering and leaving this tetrahedron. Using the symmetry of the pyrochlore lattice under translations by $\Lambda_{z}=4$ layers in the $z$ direction, the product in Eq. (13) can be rewritten as

[TABLE]

In the thermodynamic limit, the partition function is therefore given by $\left[\rho\left(\mathcal{T}^{(\Lambda_{z})}\right)\right]^{L_{z}/\Lambda_{z}}$, where $\rho$ is the spectral radius (largest absolute eigenvalue). It should be noted that, even if $\mathcal{T}_{z}$ is hermitian, in general $[\mathcal{T}_{z},\mathcal{T}_{z^{\prime}}]\neq 0$ for $z\neq z^{\prime}\pmod{\Lambda_{z}}$, and so $\mathcal{T}^{(\Lambda_{z})}$ is not hermitian.

The approximation we propose is to replace the system with an effective model of decoupled chains, whose parameters are determined self-consistently. The partition function for chain $c$ is taken as

[TABLE]

where the effective transfer operator for chain $c$ is

[TABLE]

an operator on $\mathfrak{H}_{c}$. The normalized vector $\ket{v_{c,z}}_{c}\in\mathfrak{H}_{c}$, which is interpreted as an ansatz for the configuration of chain $c$ immediately above layer $z$, is determined self-consistently as the thermal state of the one-dimensional system defined by Eq. (16). Assuming a self-consistent solution exists with at least periodicity under translation in $z$ by $\Lambda_{z}$, it is given by the normalized (right) eigenvector with largest absolute eigenvalue of

[TABLE]

For the types of phases that we aim to describe, it will in fact be sufficient to consider $\ket{v_{c,z}}_{c}=\ket{v_{c}}_{c}$ uniform along each chain. This leads to the simplification that

[TABLE]

which is therefore hermitian. For such a solution to exist, $\ket{v_{c}}_{c}$ must be a simultaneous eigenvector of $\mathcal{T}_{c,z}^{\text{mf}}$ for all $z$.

III.3 Two-chain order in spin ice

For the case with an applied field $\boldsymbol{h}$ and a chain-favoring interaction $u$, the tetrahedron transfer operator can be written explicitly as

[TABLE]

Each term corresponds to one of the six ice-rules configurations of the tetrahedron, illustrated in Fig. 3, and is associated with a Boltzmann weight implementing the terms $\mathcal{H}_{\boldsymbol{h}}$ and $\mathcal{H}_{u}$ in the Hamiltonian. (Terms breaking the ice rule could also be included; we are here treating such configurations as having zero Boltzmann weight.) The first four terms correspond to cases where the two chains maintain their configurations after passing through the tetrahedron—i.e., the two spins on each chain are aligned—while in the last two terms the chains exchange orientations.

Taking the inner product in $\mathfrak{H}_{c^{\prime}}$ gives

[TABLE]

where the expectation values

[TABLE]

are mean fields describing the magnetization and density of monopoles on chain $c^{\prime}=\mathcal{C}_{c,z}$. The effective model for chain $c$, described by the transfer operator $\mathcal{T}^{\text{mf}}_{c,z}$, involves an applied field that is a function of $m^{\text{mf}}_{c^{\prime}}$ and an effective monopole fugacity $\zeta^{\text{mf}}_{c^{\prime}}$.

This is the crux of the approximation: The exchange of orientation between chains, described by the last two terms of Eq. (20), has been replaced by terms where a single chain flips orientation, thus creating or destroying a monopole. The criterion for deconfinement, that a distant pair of monopoles can be inserted with finite free-energy cost, is trivially satisfied whenever $\zeta^{\text{mf}}_{c^{\prime}}\neq 0$. (In the analogous mean-field theory for the Bose–Hubbard model Fisher1989 , the expectation value of the bosonic annihilation operator is used as a mean field.)

In general, the solution of the self-consistent approximation requires finding simultaneous eigenvectors of $\mathcal{T}^{\text{mf}}_{c,z}$, given in Eq. (21), for all the values of $z$. For the types of ordered phases expected in our model of spin ice, however, some further simplifications are possible. In particular, to describe the MDG and saturated phases, it is sufficient to divide the chains into two “chain sublattices”, which we label $1$ and $2$, with one chain of each type meeting in every tetrahedron, as illustrated in Fig. 4. The operator $\mathcal{T}^{\text{mf}}_{c,z}$ is therefore independent of $z$, and it is sufficient to find $\ket{v_{1}}$ and $\ket{v_{2}}$ such that $\ket{v_{c}}_{c}$ is the eigenvector of $\mathcal{T}^{\text{mf}}_{c}=\bra{v_{c^{\prime}}}_{c^{\prime}}\mathcal{T}_{cc^{\prime}}\ket{v_{c^{\prime}}}_{c^{\prime}}$ with largest eigenvalue FootnotePerronFrobenius , for both $c=1$, $c^{\prime}=2$ and vice versa. Because $\mathcal{T}^{\text{mf}}_{c}$ is hermitian, this is equivalent to maximizing

[TABLE]

with respect to normalized $\ket{v_{1}}$ and $\ket{v_{2}}$.

As an aside, we note that this amounts to maximizing the matrix element of the layer transfer operator $\mathcal{T}_{z}$ in the subspace of factorizable vectors $\ket{v}=\prod_{c}\ket{v_{c}}_{c}$. In these terms, it is clear how this is related to the mean-field theory for the Bose–Hubbard model, where one minimizes the matrix element of the Hamiltonian $\mathcal{H}$. (Recall the standard connection between classical and quantum statistical mechanics, where $\mathrm{e}^{-\tau\mathcal{H}}$ is the evolution operator in imaginary time $\tau$.) In both cases, the operator is hermitian and so the matrix element is bounded by its extremal eigenvalue, but only for the BH model is the Hamiltonian constant in time and the matrix element bounded by the true ground-state energy. In the present case, the product of extremal matrix elements $\prod_{z}\bra{v}\mathcal{T}_{z}\ket{v}$ provides an approximation to the exact partition function, but this approximation is not variational.

To maximize $\Theta(\ket{v_{1}},\ket{v_{2}})$, we express both vectors as

[TABLE]

where $0\leq\theta_{c}\leq\pi$ and $-\pi\leq\phi<\pi$, in terms of which

[TABLE]

The distinct maxima of this expression correspond to phases of the mean-field theory, shown in Fig. 5.

The phase boundaries can be found by expanding in small deviations from the ordered phases and identifying the points where they cease to be stable maxima. (We have verified numerically that the ordered states maximize the matrix element globally if and only if they are local maxima.)

At low temperature, there are three confined phases with different, fixed, values of the magnetization

[TABLE]

For $\mathrm{e}^{\frac{2}{\sqrt{3}}\lvert h\rvert/T}>\mathrm{e}^{2u/T}+1$, the maximum of Eq. (26) has $\theta_{c}=0$ ($h>0$) or $\theta_{c}=\pi$ ($h<0$) for all $c$, representing a fully polarized paramagnet with $M_{z}=\pm M_{\text{sat}}$. For $\mathrm{e}^{u/T}>2\cosh\frac{2h}{\sqrt{3}T}$, it has $\theta_{c}=0$ for $c=1$ and $\theta_{c}=\pi$ for $c=2$, or vice versa. The two chain sublattices therefore have opposite orientations, $m_{1}^{\text{mf}}=-m_{2}^{\text{mf}}$, and the total magnetization vanishes. Comparison with Fig. 2 confirms that this corresponds to the MDG phase.

At higher temperature, the maximum has $0<\theta_{c}<\pi$ for all $c$ and $\phi_{1}=\phi_{2}$, corresponding to a paramagnet with continuously varying magnetization. This last phase is deconfined, having $\zeta^{\text{mf}}_{c}=\mathrm{e}^{-\mathrm{i}\phi_{c}}\sin\theta_{c}\neq 0$, and is therefore identified as the Coulomb phase.

Fig. 5 also shows the magnetization $M_{z}$ evaluated at the maximum of $\Theta(\ket{v_{1}},\ket{v_{2}})$. The resulting phase diagram consists of a set of low-temperature confined phases with fixed magnetization, surrounded by the Coulomb phase in which the magnetization varies smoothly with the parameters. The MDG phase consists of a lobe, whose width, i.e., the range of fields for which it is stable, decreases with increasing temperature.

As a benchmark for the method, we note that our model reduces for $u=0$ to that considered by Jaubert et al. Jaubert2008 (in the limit of zero monopole fugacity). In this case, one can restrict to confined phases with all chains equivalent (i.e., $\ket{v_{1}}=\ket{v_{2}}$), and it is straightforward to show that the magnetization is given by

[TABLE]

where $T_{\text{K}}=\frac{2h}{\sqrt{3}\log 2}$. Our method therefore agrees, in this case, with the results of the Bethe-lattice calculation JaubertThesis .

IV Discussion

We have shown that applying mean-field theory to the model defined by $\mathcal{H}$ in Eq. (8) produces a phase diagram containing the Coulomb phase Henley2010 as well as a fully-polarized paramagnet Jaubert2008 and the MDG phase Melko2001 . In this section, we argue that the general features in the phase diagram of Fig. 5, and in particular the lobe structure, can be understood by considering the interplay of confinement and flux in the Coulomb phase. We first consider the limit of vanishing monopole density, approximately valid at temperatures well below $\Delta$, before discussing the qualitative effect of including monopoles.

IV.1 Structure of phase diagram

The crucial observation underlying the phase structure is that confinement implies suppression of fluctuations in the magnetization $\boldsymbol{M}$, and hence vanishing uniform magnetic susceptibility $\chi$. By contrast, phases where monopoles are deconfined, such as the Coulomb phase, are magnetizable (i.e., have nonzero susceptibility $\chi$). This relationship between confinement and vanishing susceptibility (or, more generally, flux variance) occurs in any local model with a confining phase FootnoteMonopoleWinding .

The lobe structure of the phase diagram follows immediately from this observation. Starting within a confining phase and changing the applied field will eventually tip the balance in favor of the Coulomb phase, whose free energy is quadratic in the field. The higher entropy in the Coulomb phase implies that the range of applied fields for which a given ordered state is favored narrows as the temperature increases. The confined phases therefore comprise lobes with fixed magnetization, each centered around favorable values of the applied field. (The fully polarized phases are exceptions, extending to arbitrarily large $\lvert h\rvert$ as a result of their saturated magnetization.) This qualitative structure is confirmed in the mean-field phase diagram, Fig. 5. With more realistic interactions, we expect more lobes with other commensurate values of magnetization to appear in the gaps between the lobes shown, as indeed seen in the MC results of Ref. Lin2013 .

The general phase structure can also be understood through the mapping to a quantum model of bosons, described in Section II.4 and Ref. SpinIceCQ . For an extended hard-core Bose–Hubbard model in the grand canonical ensemble, the phase diagram contains a set of Mott-insulator lobes at small hopping, consisting of different possible ordered states Fisher1989 ; Sachdev2011 ; Kovrizhin2005 ; Kurdestany2012 . Returning to spin ice, each of these corresponds to a confined phase with order in the spin degrees of freedom and fixed magnetization.

The susceptibility strictly vanishes in the ordered phases only in the thermodynamic limit at zero monopole fugacity $\zeta=\mathrm{e}^{-\Delta/T}$. With finite monopole cost, confinement of monopoles is no longer precisely defined and the susceptibility becomes nonzero, though remains exponentially small at low temperature. The plateaux are therefore rounded for nonzero $\zeta$. Those phases that have conventional order, such as the MDG phase, remain separated from the paramagnet by a symmetry-breaking transition.

IV.2 Phase transitions

Certain properties of the phase transitions can also be determined on general grounds, following similar arguments to those for the Bose–Hubbard model Fisher1989 ; Sachdev2011 , but with density replaced by magnetization. We first note that phase transitions can be classified by whether the magnetization changes across the transition. The magnetization necessarily changes in the case where the ordered phase has saturated magnetization, but it can also change across a transition into a magnetically ordered phase.

At the tip of the lobe, the magnetization is identical in the two phases (ordered within and Coulomb outside), and so this is the only point at which the magnetization stays the same across the transition. The argument is that, inside the ordered phase but sufficiently close to the tip, an arbitrarily small change in the applied field pushes the system into the Coulomb phase, regardless of its sign. If the two phases had different magnetizations at this point, their linear coefficients in the free energy would be different and this could not be the case. (This argument is essentially identical to the corresponding one for the Bose–Hubbard phase diagram Fisher1989 .)

For any continuous transitions with fixed magnetization, critical properties can be calculated using the approach of Ref. SpinIceHiggs , but with a modified projective symmetry group due to the nonzero magnetization. While MDG found a strongly first-order transition Melko2001 and the transition into the partially polarized phase was seen to be of first order in Ref. Lin2013 , it remains possible that transitions in the presence of different perturbations could be continuous. (The possibility of unconventional transitions between a superfluid and fractional Mott insulators has been discussed in the extended Bose–Hubbard model Balents2005 .)

When the magnetization changes across the transition, i.e., everywhere except at the tip, the transition is in the same universality class as the Kasteleyn transition Kasteleyn1963 ; Jaubert2008 . It should be noted, however, that only the transitions into the saturated phases, with $M_{z}=\pm M_{\text{sat}}$, are true Kasteleyn transitions, which are distinguished by the complete absence of fluctuations on the ordered side.

V Conclusions

We have presented a mean-field theory designed to study confinement transitions, based on the analogy between the confinement criterion and long-range order in conventional ordering transitions. The method has been applied to a simplified model of classical spin ice and shown to produce a phase diagram containing, besides the Coulomb spin liquid Henley2010 , a set of magnetization plateaux. These include the saturated paramagnet Jaubert2008 as well as the MDG phase expected to occur in dipolar spin ice at low temperature Melko2001 . We have argued that the general properties of the phase diagram follow from the interplay between confinement and magnetization that characterizes the ice-rule states and pointed out the analogy with the phase diagram of the extended Bose–Hubbard model Fisher1989 .

The simplified model that was studied here reproduces the MDG phase and is particularly amenable to the mean-field approach. It would be desirable to extend the approach to allow for more realistic further-neighbor interactions between spins, particularly the dipolar interactions that are known to be significant in the classical spin-ice materials. Further-neighbor interactions can also stabilize the partially magnetized phase observed using MC in Ref. Lin2013 . (This phase can be described in terms of the chains shown in Fig. 4, but distinguishes four chain sublattices, rather than two.)

A transfer operator written in the form of Eq. (20) has the limitation, however, that it can represent interactions only within a single tetrahedron. (One cannot even define the layer transfer operator for general interactions.) A potential route to include further-neighbor interactions is through an additional mean-field decoupling, replacing a pairwise interaction $V_{ij}\sigma_{i}\sigma_{j}$ with terms such as $V_{ij}\sigma_{i}\langle\sigma_{j}\rangle$, representing a mean field acting on site $i$. This extension, and others such as including nonzero monopole fugacity $\zeta$, are left to future work.

Acknowledgements.

This work was supported by EPSRC Grant No. EP/M019691/1.