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

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.


I. INTRODUCTION
Classical spin liquids [1], such as the Coulomb phase [2] in spin ice [3] and related systems, are examples of phases whose behavior is not captured by the standard Landau picture of broken symmetries [4]. 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 [5].
In the case of the classical spin ices [3], 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 [6], 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 [7,8].
Spin ice also has the unusual feature that the magnetization is a topologically constrained quantity, and so topological order can be probed directly [5]. In particular, within the low-energy configuration space relevant at low temperatures, any local dynamics conserves the uniform magnetization [2]. 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 [9,10].
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 [2], 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. [11]. A second, where the ordered state has saturated polarization, was considered in Refs. [9,10]. 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 [12], making use of a mapping from spin ice to a quantum model of bosons [10]. 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 [12,13], 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 [14][15][16], while it has been suggested that the zero-field ground state of quantum spin ice may have nonzero magnetization [17].
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.

A. Nearest-neighbor spin ice
The spin ices [3] 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 Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 , which are well described in terms of classical spins S i . Each spin is subject to strong easy-axis anisotropy along the local 111 directionn 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 K) crystal-field splitting between states with maximum projection alongn 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 [3]. For the low temperatures of interest, it is therefore sufficient to treat the spins as effectively constrained to S i = σ ini , where σ = ±1 is a classical Ising variable.
Each tetrahedron t in the pyrochlore lattice can be labeled according to its orientation t = ±1, and the lattice structure is such that all neighbors of t have orientation − t . The fixed unit vectorn i for each FIG. 1. A pair of tetrahedra forming part of the pyrochlore lattice. In the classical spin-ice compounds, magnetic moments S i (shown as green arrows) lie at the sites i of the lattice and are constrained to point along the local 111 axes (dashed lines) that join the centers of adjacent tetrahedra. pyrochlore site i, directed along the line joining the centers of its two tetrahedra, can be chosen to point towards the tetrahedron with t = −1, and so t σ i indicates whether S i points out of (+1) or into (−1) tetrahedron t. With this convention,n i ·n j = − 1 3 and hence S i · S j = − 1 3 σ i σ j , for all nearest neighbors i and j.
The minimal Hamiltonian that captures the spin-liquid physics of the classical spin ices [5] contains only interactions between nearest-neighbor pairs i j , The coupling between moments S i is ferromagnetic (J > 0), incorporating the net effect of the dipolar interactions between nearest neighbors and (antiferromagnetic) exchange [3]. The effective interactions between the Ising variables σ 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 H nn [5]; we will return to these in Section II C.
Because two sites are nearest neighbors if and only if they share a tetrahedron, the interactions can be rewritten as where t includes all tetrahedra (of both orientations) and is referred to as the (fictional) "charge" on tetrahedron t in terms of a sum over its sites i ∈ 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 [3,5]; all configurations that do so at every tetrahedron are degenerate minimal-energy states of H nn .
Tetrahedron configurations with three spins in and one out, or vice versa, have energy ∆ = + 2 3 J above the ice-rule configurations; they have charge Q t = ±1 and are hence referred to as "monopoles". (Note that t Q t = 0, and so all configurations are globally charge-neutral.) For Dy 2 Ti 2 O 7 , the effective nearestneighbor coupling is 1 3 J = 1.1 K [3], corresponding to a monopole cost of ∆ = 2.2 K. The corresponding Boltzmann weight ζ = e −∆/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 [6].
For temperature T significantly below ∆, the system is effectively constrained by the ice rule, and the ensemble of ice-rule configurations therefore determines its behavior [3]. The number of such configurations grows exponentially with the total number of spins N, and so the entropy density remains nonzero even well below ∆, 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 [2], and is an example of a classical spin liquid [1].
The long-wavelength properties of the Coulomb phase can be captured by coarse-graining the vector spins S to give a continuum field B(r) referred to as the effective "magnetic field" [2]. The ice-rule constraint applied to S implies that 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 [2].
The ice-rule configurations also display interesting topological properties [5]. 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 M = i S i if it has nontrivial winding number. For other systems that host Coulomb phases [2], 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.

B. Confinement transitions
A second important property of the Coulomb phase is that monopoles, defects in the ice rule, are deconfined [6]. To make this statement precise, imagine imposing a pair of monopoles with opposite charges ±1 at tetrahedra t and t in a background that otherwise obeys the ice rule. Let C t,t denote the set of spin configurations {σ i } compatible with this charge configuration. Their number Z m (r t,t ) = |C t,t | is a function of the separation r t,t that decreases with |r t,t | but reaches a finite limit as |r t,t | → ∞. The result is an effective entropic interaction 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 B, coarse-graining predicts that the inter- More generally, consider the partition function 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 m (r t,t ) increases without limit as |r t,t | → ∞, and so Z m (r t,t ) → 0, is referred to as a confinement transition [7,8]. (In terms of the effective continuum field B, Gauss' law implies that an imposed pair of charges must be joined by a fixed quantity of flux. Linear confinement, where U m (r) ∼ |r| for large |r|, 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 [18], there is no local order parameter that provides a signature of confinement. The confinement criterion, that Z m (r t,t ) approaches zero as |r t,t | → ∞ 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 ζ t , in terms of which [19] where Z = Z[ζ t ]| ζ · =0 is the partition function restricted to ice-rule configurations. The relationship between the monopole distribution function G m and the fugacity ζ expressed by Eq. (7) is exactly analogous to the fluctuation-dissipation theorem relating, for example, a spin-spin correlation function and an applied field. confinement transition in terms of an effective self-consistent fugacity.
It should be noted that the confinement distinction is only sharp in the limit ζ = 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 [7,8]. If the phase transition also involves a different type of order (such as spontaneous symmetry breaking), then it may survive for ζ > 0, but with a different universality class [18].

C. Perturbations to NNSI
Within the nearest-neighbor model H 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 where H h is the coupling to an applied magnetic field and H u contains additional interactions described below.
Coupling to a magnetic field h is described by a Zeeman term where a factor dependent on the size of the magnetic moments has been included in the definition of h. The case that we consider here has h = hu z where u z is a unit vector along [100], and so h · S i = hu z · σ ini = ± 1 √ 3 hσ i . The coupling can therefore be written as where Using the choice described in Section II A, wheren i points towards the tetrahedron with t = −1, η i alternates in sign on successive (100) planes of spins. The component of the magnetization along the field can likewise be expressed as The perturbation H h alone can drive a confinement transition, referred to as a Kasteleyn transition [9].
For ∆ h 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 sat ≡ 1 √ 3 N even for nonzero T/h (but in the limit T/∆ → 0). Above a temperature T K , however, the energy cost is outweighed by the gain in entropy, also ∝ 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/∆ = 0, the transition is replaced by a crossover [8,9].
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 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 [21] and long-range dipolar interactions. The latter have been shown using Monte Carlo [22] 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 [23].
The actual perturbation 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 where 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 H u .
Unlike dipolar interactions, the perturbation 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 × 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. [11], 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. between nearby sites and further-neighbor interactions, with parameters related in nontrivial ways to those of the spin model [10].
As µ is increased, the system passes from the vacuum through a quantum phase transition [13] 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 [10].) For sufficiently large positive µ, the system is driven into a fully occupied Mott insulator, with ρ = 1, corresponding to a saturated paramagnet with M z = −M 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 [10].
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 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 [26], but are not captured by the mean-field theory used here.
The phase structure of 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 [27], and also respect the spatial structure, unlike Bethe-lattice calculations [9].
Here we apply an approach that is analogous to the standard mean-field theory for the Bose-Hubbard model [12], 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. [28,29].) 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 selfconsistently.
In standard mean-field theory for a symmetry-breaking transition [20], 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 ζ 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 ζ 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 ζ mf , we set the true monopole fugacity ζ = e −∆/T to zero. Nonzero ζ could straightforwardly be included in the approximation and is analogous to standard Weiss mean-field theory in the presence of an applied magnetic field [20]. 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 P z . We also define C c,z as the chain that meets chain c in a tetrahedron at layer z. The pyrochlore lattice is symmetric under translations by Λ z = 4 layers in the [100] direction [10], and so P c+Λ z = P c and C c,z+Λ z = 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" [30] implementation of the hopping in the effective quantum Hamiltonian has been shown to be irrelevant for the critical properties [10], but must be handled with some care when implementing the mean-field theory. The partition function can then be written as where L z is the system size (number of layers) in the z direction and is the transfer operator for the layer of tetrahedra at vertical position z.
In the thermodynamic limit, the partition function is therefore given by ρ T (Λ z ) L z /Λ z , where ρ is the spectral radius (largest absolute eigenvalue). It should be noted that, even if T z is hermitian, in general [T z , T z ] 0 for z z (mod Λ z ), and so T (Λ 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 where the effective transfer operator for chain c is an operator on H c . The normalized vector v c,z c ∈ 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 Λ z , it is given by the normalized (right) eigenvector with largest absolute eigenvalue of For the types of phases that we aim to describe, it will in fact be sufficient to consider v c,z c = |v c c uniform along each chain. This leads to the simplification that which is therefore hermitian. For such a solution to exist, |v c c must be a simultaneous eigenvector of T mf c,z for all z.

C. Two-chain order in spin ice
For the case with an applied field h and a chain-favoring interaction u, the tetrahedron transfer operator can be written explicitly as 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 H h and 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 H c gives where the expectation values are mean fields describing the magnetization and density of monopoles on chain c = C c,z . The effective model for chain c, described by the transfer operator T mf c,z , involves an applied field that is a function of m mf c and an effective monopole fugacity ζ mf c . 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 ζ mf c 0. (In the analogous mean-field theory for the Bose-Hubbard model [12], 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 T 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 T mf c,z is therefore independent of z, and it is sufficient to find |v 1 and |v 2 such that |v c c is the eigenvector of T mf c = v c | c T cc |v c c with largest eigenvalue [31], for both c = 1, c = 2 and vice versa. Because T mf c is hermitian, this is equivalent to maximizing with respect to normalized |v 1 and |v 2 .
As an aside, we note that this amounts to maximizing the matrix element of the layer transfer operator where e −τH is the evolution operator in imaginary time τ.) 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 z v| T z |v provides an approximation to the exact partition function, but this approximation is not variational.
To maximize Θ(|v 1 , |v 2 ), we express both vectors as where 0 ≤ θ c ≤ π and −π ≤ φ < π, in terms of which The distinct maxima of this expression correspond to phases of the mean-field theory, shown in Fig , it has θ c = 0 for c = 1 and θ c = π for c = 2, or vice versa. The two chain sublattices therefore have opposite orientations, m mf 1 = −m mf 2 , and the total magnetization vanishes. Comparison with Fig. 2 confirms that this corresponds to the MDG phase.
At higher temperature, the maximum has 0 < θ c < π for all c and φ 1 = φ 2 , corresponding to a paramagnet with continuously varying magnetization. This last phase is deconfined, having ζ mf c = e −iφ c sin θ c 0, and is therefore identified as the Coulomb phase. As a benchmark for the method, we note that our model reduces for u = 0 to that considered by Jaubert et al. [9] (in the limit of zero monopole fugacity). In this case, one can restrict to confined phases with all chains equivalent (i.e., |v 1 = |v 2 ), and it is straightforward to show that the magnetization is given by where Our method therefore agrees, in this case, with the results of the Bethe-lattice calculation [32].

IV. DISCUSSION
We have shown that applying mean-field theory to the model defined by H in Eq. (8) produces a phase diagram containing the Coulomb phase [2] as well as a fully-polarized paramagnet [9] and the MDG phase [22]. 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 ∆, before discussing the qualitative effect of including monopoles.

A. Structure of phase diagram
The crucial observation underlying the phase structure is that confinement implies suppression of fluctuations in the magnetization M, and hence vanishing uniform magnetic susceptibility χ. By contrast, phases where monopoles are deconfined, such as the Coulomb phase, are magnetizable (i.e., have nonzero susceptibility χ). This relationship between confinement and vanishing susceptibility (or, more generally, flux variance) occurs in any local model with a confining phase [33].
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 |h| 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. [24].
The general phase structure can also be understood through the mapping to a quantum model of bosons, described in Section II D and Ref. [10]. 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 [12,13,35,36]. 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 ζ = e −∆/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 ζ. Those phases that have conventional order, such as the MDG phase, remain separated from the paramagnet by a symmetry-breaking transition.

B. 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 [12,13], 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 [12].) For any continuous transitions with fixed magnetization, critical properties can be calculated using the approach of Ref. [11], but with a modified projective symmetry group due to the nonzero magnetization.
While MDG found a strongly first-order transition [22] and the transition into the partially polarized phase was seen to be of first order in Ref. [24], 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 [37].) 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 [9,38]. It should be noted, however, that only the transitions into the saturated phases, with M z = ±M 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 [2], a set of magnetization plateaux. These include the saturated paramagnet [9] as well as the MDG phase expected to occur in dipolar spin ice at low temperature [22]. 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 [12].
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 furtherneighbor 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. [24]. (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 i j σ i σ j with terms such as V i j σ i σ j , representing a mean field acting on site i. This extension, and others such as including nonzero monopole fugacity ζ, are left to future work.

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