Ferromagnetic Coulomb phase in classical spin ice

Spin ice is a frustrated magnetic system that at low temperatures exhibits a Coulomb phase, a classical spin liquid with topological order and deconfined excitations. This work establishes the presence of a Coulomb phase with coexisting ferromagnetic order in a microscopic model of classical spin ice subject to uniaxial lattice distortion. General theoretical arguments are presented for the presence of such a phase, and its existence is confirmed using Monte Carlo results. This example is used to illustrate generic properties of spin liquids with magnetic order, including deconfinement of monopoles, signatures in the neutron-scattering structure factor, and critical behavior at phase transitions. An analogous phase, a superfluid with spontaneously broken particle-hole symmetry, is demonstrated in a model of hard-core lattice bosons, related to spin ice through the quantum-classical correspondence.


I. INTRODUCTION
Spin liquids are phases where magnetic degrees of freedom exhibit strong local correlations, despite the persistence of large fluctuations, 1 of either quantum mechanical or thermal origin. They occur at low temperature in certain frustrated systems, where interactions are large compared to thermal fluctuations, but mutual competition between the interactions prevents formation of a rigidly ordered configuration. Spin liquids have been of theoretical interest for several decades, 2,3 but evidence for their existence in physical systems, 4,5 and even microscopic models, [6][7][8] is considerably more recent.
While spin liquids are often distinguished from conventional low-temperature phases, such as ferromagnets, by the fact that they lack magnetic order, their defining characteristics go beyond the mere absence of conventional order. A precise definition of a quantum spin liquid (QSL) can be phrased in terms of long-range entanglement, 1 while the Coulomb phase, 9 the classical spin liquid (CSL) that is of primary interest here, can be defined through deconfinement of fractionalized "monopole" excitations. 8 Experimental evidence exists for a Coulomb phase in the spin-ice compounds, which can be treated as classical at relevant temperatures. 10 These definitions provide positive characterizations for QSL and CSL phases, and also make clear the possibility of a magnetically ordered spin liquid, in which spin-liquid phenomena coexist with conventional symmetry-breaking order. Some examples of such phases have been reported in the theoretical literature: Mean-field studies of quantum spin ice 11 identified an ordered QSL, referred to as a "Coulomb ferromagnet", although quantum Monte Carlo simulations have not revealed such a phase. 12 Recent work 13 has also demonstrated the possibility of antiferromagnetic order coexisting with a CSL.
The compatibility of magnetic order and spin-liquid phenomenology also allows for the existence of phase transitions between ordered and disordered spin liquids. One might anticipate novel critical behavior at such transitions, since it is known that transitions from spin liquids into conventional ordered phases can transcend the usual Landau paradigm. 9 This work demonstrates that a ferromagnetic Coulomb phase can occur in a model of classical spin ice, and provides a detailed study of this phase and the associated transitions.
Theoretical arguments, including mapping to a related quantum model, are used to show that such a phase exists and that it can be reached through a continuous transition from the paramagnetic Coulomb phase. We present Monte Carlo (MC) results that confirm both of these statements, and illustrate the generic properties of ordered spin liquids, including the structure factor for elastic neutron scattering.
We also consider the critical behavior at the ordering transition and predict that, despite the Ising nature of the order parameter and the presence of only short-range interactions in the microscopic model, the transition should belong in the mean-field universality class, as a result of coupling to the effective gauge-field fluctuations of the spin liquid. While the numerical results are consistent with this prediction, larger system sizes would be required for a definitive confirmation of the universality class. This phase transition provides another interesting example of the diversity of critical phenomena that exists in the neighborhood of spin-liquid phases.

Outline
In Section II, the model of spin ice is introduced, and a choice of perturbations that lead to a ferromagnetic Coulomb phase is motivated. The basic structure of the phase diagram is then illustrated using MC results, showing the appearance of such a phase at intermediate temperatures for certain values of the parameters. In Section III, the phase in question is studied in detail, to confirm that it has nonzero magnetization while simultaneously exhibiting the characteristic features of a Coulomb phase. Sections IV and V consider in turn the critical behavior at the higher-and lower-temperature transitions out of this ferromagnetic Coulomb phase.
We conclude in Section VI, by summarizing the features that are expected to be generic to ordered spin liquids, both quantum and classical, and discuss briefly the effect of a nonzero density of magnetic monopoles. In the Appendix, the classical-quantum mapping developed in Ref. 14 is applied to this system, and the resulting quantum model is related to a problem of hard-core quantum bosons studied by Rokhsar and Kotliar. 15

A. Nearest-neighbor model of spin ice
The spin-ice materials 8,10 Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 are well described by a model of classical spins S i on the sites i of a pyrochlore lattice, a network of corner-sharing tetrahedra. Each spin is subject to a strong easy-axis anisotropy constraining it to point parallel or antiparallel to the local 111 axis joining the centres of adjacent tetrahedra, S i = ±n i . Including only nearest-neighbor interactions, the Hamiltonian can be written as where J i j is a ferromagnetic coupling between nearestneighbor sites i j of the lattice. In the unperturbed model, the interaction is uniform, J i j = J > 0, and favors those states where, of the four spins on each tetrahedron, two point in and two point out. The latter condition is referred to as the "ice rule" and selects a set of states that is degenerate in the nearest-neighbor model and whose number grows exponentially with the number of spins. While a more realistic microscopic model than H nn would also include dipolar interactions, 10 their effect is primarily to renormalize J, with only a small splitting of the ice-rule states. 16 We will mostly concentrate on the limit where the ice rule is enforced as a constraint, represented by Eq. (1) with temperature T J. Assuming ergodicity within the ice-rule manifold, the system in this limit exhibits a Coulomb phase, in which the spins are disordered but highly correlated. Replacing the spins S i by a continuous vector field B(r) and the ice rule by ∇ · B = 0 leads to an effective coarse-grained description for this phase. 9 A quadratic action for the "magnetic field" B correctly describes the long-wavelength neutron scattering at low temperature in spin ice, 17 and predicts that monopoles in B, corresponding to single tetrahedra where the ice rule is broken, are deconfined. 8 Much of the physics is in fact qualitatively unaltered by a small density of such defects (see Section VI), and their effects on the critical properties can be understood by treating monopole fugacity as a relevant perturbation (in the renormalization-group sense). 18,19 An important property of the ice-rule states for present purposes is that they obey a topological constraint on the magnetization: 8,9 Starting from any ice-rule state and flipping a spin S i breaks the ice rule on the two tetrahedra shared by site i. The only updates that connect configurations within the ice-rule manifold are those that involve flipping a set of spins aligned head-to-tail along a closed loop. Any such update for a contractible loop preserves the magnetization density, where N s = i 1 is the number of spins. Changing the magnetization while remaining within an ice-rule state in fact requires flipping spins along a loop that spans the system (assuming periodic boundary conditions). Sets of states Three configurations of a single tetrahedron, and their energy E t in the nearest-neighbor Hamiltonian H nn , Eq. (1), with J i j given in Section II B. Pairs of spins situated in the same horizontal plane, indicated with dashed (red) lines, have reduced coupling J i j = J − 3p, while others have J i j = J. All three configurations obey the ice rule, having two spins pointing in and two pointing out. The first two are lowest-energy configurations for a single tetrahedron (since the antiferromagnetically aligned pairs are those with reduced coupling), while the one on the right is one of the remaining four icerule configurations whose energy is higher by 4p. Excitations above the ground state are described by strings of spins flipped relative to a fully polarized configuration, and increase the energy by 4p per tetrahedron. When a pair of strings pass through the same tetrahedron, all four spins are inverted and the energy is again minimized; the strings therefore feel an attractive interaction.
with the same magnetization therefore constitute "topological sectors", 8 disconnected by local updates. This topological conservation law is broken by a nonzero density of monopoles, but remains approximately valid, and conceptually useful, at low temperatures. Nonzero magnetic susceptibility χ requires that the system fluctuates between different sectors; 20 in the thermodynamic limit, one can therefore distinguish "incompressible" phases with χ = 0 from those with χ > 0.

B. Uniaxial distortion
To split the energy of the six ice-rule states on a given tetrahedron requires breaking the cubic symmetry of the pyrochlore lattice. Following Ref. 21, we consider an explicit uniaxial symmetry breaking, with J i j = J − 3p (J > 3p > 0) for pairs of spins whose relative displacement lies in the (001) plane and J i j = J for all others. (Such a perturbation could be effected in experiment by application of uniaxial pressure. 21 ) As illustrated in Fig. 1, the result is to favor the two configurations where the total (vector) spin of the tetrahedron is along the [001] axis, whose energy is reduced by 4p compared to the other four. In contrast to the case of an applied field, 14,22 an Ising symmetry remains; there are two degenerate lowestenergy states, with all spins on all tetrahedra maximally polarized, consistent with the local easy axes, either along ("up") or against ("down") the [001] direction.
We first briefly review the phase structure of the model H nn with this distortion; readers are referred to Refs. 21 and 23 for more details. For p T J, the system remains in the Coulomb phase, while below a critical temperature T c = 4p(ln 2) −1 it becomes a fully polarized ferromagnet. At the transition, the up-down symmetry is broken and the magnetization along the [001] direction, M z , becomes nonzero. While an Ising order parameter can naturally be defined, the transition in the ice-rules limit has quite different properties from the standard Ising universality class. Starting from either of the fully polarized states, the only closed loops are "strings" spanning the system in the [001] direction, which cost energy proportional to the (linear) system size L. The transition occurs when the entropy of a single string, also ∝ L, outweighs the energy, and so its free energy changes from positive to negative; the string density then increases from zero to nonzero.
As a consequence, the system on the lower-temperature side of the transition is fully polarized, with zero string density, as in the related case of an applied field. 14,22 A crucial distinction in this case is that two strings feel an effective attraction when sharing a tetrahedron (see Fig. 1). At the critical point, this exactly balances the entropy cost of the excluded volume due to the hard-core interactions between strings. In fact, as Jaubert et al. 21 have shown, the free energy at T = T c as a function of string density is precisely constant (in the thermodynamic limit). Because each string consists of a fixed number of flipped spins relative to the starting configuration, this implies that the free energy is independent of magnetization. As the temperature increases through the transition, the global minimum of F(M z ), which can be interpreted as a Landau function, jumps from M z = ±M sat to M z = 0. (The resulting discontinuity in the magnetization is illustrated below in Fig. 3.) Since all coefficients in the Landau free energy vanish at the transition, this has been referred to as "infinite-order multicriticality". 21

C. Additional interactions
Given the magnetization-independent free energy at the transition, it is clear that any perturbation that produces a positive fourth-order coefficient in the Landau function should lead to an intermediate phase with 0 < |M z | < M sat . While this argument does not provide a prescription for constructing appropriate perturbations, one expects on general grounds that a sufficiently long-ranged four-spin interaction will have this effect. (As will also be demonstrated, a quartic coefficient with opposite sign should lead to a first-order transition.) As we detail in the following, MC results in fact demonstrate that it is sufficient to add a four-spin interaction acting between tetrahedra on opposite sides of a hexagonal loop, as illustrated in Fig. 2. The perturbation used throughout this work can be written explicitly as where and S t ≡ i∈t S i is the total (vector) spin on tetrahedron t. The sum in Eq. (3) is over pairs of tetrahedra {tt } across a hexagon (see Fig. 2), and the summand is one if both tetrahedra have all spins polarized in the same vertical direction, and zero otherwise. (Note that, while this expression apparently involves eight spins, it is equivalent to a four-spin interaction under projection into the ice-rule states. This form of the interaction is partly motivated by the quantum mapping, described in the Appendix.) Regarding the choice of H 4s , it is not the goal of this work to classify the various types of interactions according to whether they produce a ferromagnetic Coulomb phase, and we are not aware of a general argument that would allow for such a classification. 24 (The search for appropriate interactions is in any case better informed by experimental evidence about which interactions occur in particular materials.) Rather, the goal here is to study a particular case where such a phase is known to exist, and elucidate those properties of the phase and its transitions that are expected to be universal, or at least qualitatively generic.
Plots of the magnetization as a function of temperature, for V 4 positive, negative, and zero, are shown in Figs. 3 and 4. These results were produced using MC simulations based on a directed-loop algorithm. 25,26 The lattice consists of L × L × L cubic unit cells, each containing 4 tetrahedra of each orientation, and hence 16 spins. For V 4 ≤ 0, a step is observed in the magnetization, from essentially fully saturated, with small deviations due to finite-size effects, 21 to zero within error bars. In both cases, there is a jump from saturation to zero magnetization, at a transition temperature indicated with a vertical line. For each temperature, the spontaneous magnetization is found by applying a weak field along the z direction and extrapolating to zero field. This step is accompanied by a single peak in the specific heat (not shown), whose height grows with system volume, indicating a single first-order transition. By contrast, when V 4 > 0 (Fig. 4), there are clearly three distinct regimes as the temperature T is lowered. The hightemperature phase is paramagnetic, with M = 0, and is the usual Coulomb phase observed at T J in spin ice. The magnetization first becomes nonzero at T c> before reaching its saturation value at T c< . While the system is ferromagnetic for all T < T c> , it is a saturated ferromagnet, with M z = ±M sat , only below T < T c< . As shown in Fig. 5, the variance of the energy (proportional to the specific heat) in this case displays two peaks, both at most weakly diverging with L, consistent with a pair of continuous transitions. Fig. 6 shows histograms of the energy and magnetization for L = 16, V 4 /T = 0.05, and T/p = 3.31, near the higher- temperature peak of the energy variance. The unimodal structure of the energy distribution confirms the continuous nature of the transition, while the two peaks of the magnetization indicate that this is a symmetry-breaking transition into a state with nonzero but unsaturated magnetization. This should be contrasted with the case of V 4 = 0, where the magnetization histogram is flat at the transition. 21 Fig. 7 shows the case of V 4 < 0, where the transition is of first order, with a bimodal structure in the energy and coexisting peaks in the magnetization distribution, at both M z = 0 and M z = ±M sat .

III. INTERMEDIATE PHASE
Having established the presence of a pair of phase transitions when V 4 > 0, we now turn to the intermediate phase in the temperature range T c< < T < T c> . It will be argued that this phase shares the essential spin-liquid features of the Coulomb phase above T c> , but is distinguished by a nonzero spontaneous magnetization.
The presence of a nonzero but unsaturated magnetization in the intermediate phase is evident from Figs. 4 and 6. Continuously changing magnetization implies that the magnetic susceptibility is nonzero, and hence that there are fluctuations between different topological (magnetization) sectors. This fact alone distinguishes the intermediate phase from the low-temperature saturated ferromagnet, where the flux stiffness vanishes in the thermodynamic limit, and there are no topological-sector fluctuations. 8,20 Two phenomena that are characteristic of the Coulomb phase are deconfinement and algebraic spin-spin correlations; these are discussed in turn in the following subsections.

A. Monopole distribution function
A single tetrahedron at which the ice rule is broken (i.e., where the number of spins pointing in and out differs) cor- responds to a monopole in the continuous vector field B(r). Such defects are rare for T J, and, at least as a first approximation, we treat the density of thermally excited monopoles as vanishing.
It is useful to consider, however, the introduction of a single pair of oppositely charged monopoles into an otherwise defect-free background. The effective interaction between the pair, induced by the fluctuations of the surrounding spins, allows one to distinguish spin-liquid phases from others such as the saturated ferromagnet. In the Coulomb phase, the monopoles are subject to an effective Coulomb interaction, with a finite limit for large separation. The saturated ferromagnet is, by contrast, a confining phase, in which the free energy of a pair of monopoles grows without bound as their separation increases. 8,18,19 To determine directly whether monopoles are deconfined, one can define the monopole distribution function G m (r + , r − ) as the partition function calculated in the presence of a pair of monopoles of opposite charge at r ± . (More explicitly, the ensemble is constrained so that all tetrahedra obey the ice rule, apart from those at r ± , where three spins point out and one points in, and vice versa.) This function, which is related to the effective interaction between monopoles U m by G m = e −U m /T , has a nonzero limit for infinite separation |r + − r − | only when monopoles are deconfined. In a confined phase, it instead decays exponentially to zero. In a finite system, these asymptotic behaviors are observed only for separations much less than the system size L. Finitesize effects can be controlled by fixing the ratio |r + − r − |/L and observing the scaling with L. Fig. 8 shows the ratio of the monopole distribution function calculated at R max , the largest displacement possible for L 3 cubic unit cells (L even) with periodic boundaries, and at R max,z , the maximum separation along the z direction (|R max | = √ 3|R max,z |). The ratio approaches unity with increasing system size for all T > T c< , while it decays to zero below T c< , indicating confinement. No qualitative difference is seen when crossing the highertemperature phase boundary at T c> , demonstrating that the intermediate phase, in common with the standard Coulomb phase above T c> , exhibits deconfinement of monopoles.
The form of the effective interaction U m (r) = −T ln G m (r) is determined by the approach of G m (r) to its asymptotic value for large separation r. Fig. 9 shows U m for temperatures within the intermediate phase and above T c> . In both cases, the interaction is anisotropic, because the spatial symmetry is reduced by the applied pressure (and H 4s ). Up to finitesize effects, the interaction can be fit to the Coulomb form, The most direct experimental signature of the Coulomb phase is the presence of "pinch points" in the neutronscattering structure factor, which reflect the algebraic (dipolar) correlations between the spins. 9,10 These features are clearest in the spin-flip component of polarized neutron-scattering data with incident polarization along [110]. 17 The corresponding structure factor is where S µν (Q) is the Fourier transform of the two-spin correlation function S µ (r)S ν (r ) and is a unit vector orthogonal to both the wavevector Q and the incident neutron polarization P. This structure factor is shown in Fig. 10, for Q in the (hh ) plane and P along [110]. Pinch points are visible for all T > T c< , with no qualitative change at T c> , showing that the  27 The diffuse scattering is completely suppressed in the low-temperature saturated ferromagnet, and only the Bragg peaks remain.

IV. HIGHER-TEMPERATURE TRANSITION
The previous sections have established that the phase at T c< < T < T c> is a spin liquid with ferromagnetic order, and that it is connected to the high-and low-temperature phases by continuous transitions. In this section and the following, the critical properties of these two transitions will be addressed in turn, using analytical arguments supported by numerical results.

A. Critical theory
Near the transition, at T = T c> , between the paramagnetic and ferromagnetic Coulomb phases, the magnetization is far from saturation and so the discrete nature of the spins is presumably not important. Replacing the discrete spins by a continuous vector field 9 B(r), the partition function can be expressed as The coefficient κ is the flux stiffness in directions transverse to the applied pressure, while α > 0 represents the effect of the uniaxial pressure, enhancing fluctuations along the z direction.
(Higher-order terms have been omitted.) Using a Hubbard-Stratonovich field Φ to decouple the anisotropy term, this can be replaced by The real scalar field Φ has Ising symmetry and provides an order parameter for the transition, taking a nonzero value in the ferromagnetic phase. Integrating out B induces dipolar interactions for Φ, giving an effective description that is equivalent to that of Ising spins with dipolar couplings. A similar connection between the dipolar correlations in the Coulomb phase and effective dipolar interactions at a critical point has been noted in Ref. 28. The 3D Ising transition with dipolar interactions 29-31 is at its upper critical dimension, and so shows mean-field critical exponents with logarithmic corrections. In particular, the specific-heat, order-parameter, susceptibility, and correlationlength exponents take the values α = 0, β = 1 2 , γ = 1, and ν = 1 2 , respectively. It should be noted that scaling remains isotropic at this transition, in the sense that all spatial dimensions scale with the same exponents. For example, the correlation lengths in the directions parallel and perpendicular to the applied pressure diverge with the same exponent ν, though with different (nonuniversal) prefactors. This is in contrast to the anisotropic scaling at the lower-temperature transition (see Section V).

B. Numerical results
To determine the critical temperature T c> and find values of the exponents, it is convenient to identify a quantity whose scaling dimension vanishes, for which curves with different L coincide at the transition. While the Binder cumulant of the magnetization provides such a quantity for this ordering transition, it is difficult to calculate accurately, as a result of the topological constraints on the magnetization, which suppress fluctuations of the latter.
We instead consider the quantity L M 2 z , which, as a result of the scaling form where Ψ is a universal function, is expected to have zero scaling dimension for this transition. (This quantity is equal, up to powers of L, to the flux stiffness Υ, which is not expected to have vanishing scaling dimension at a transition between two spin liquids.) As shown in Fig. 11, L M 2 z plotted as a function of T/p indeed has a crossing point for large system sizes. Using the crossings for successive L values, we estimate (T/p) c = 3.3509(3) for V 4 /T = 0.05.
While the observed crossing is consistent with the meanfield exponents, it is also compatible with the Ising universality class, which has 32 d − γ/ν = 1.0366 (8). We can go some way to excluding this possibility by calculating the correlation-length exponent ν, which, for the Ising class, takes the value 32 ν = 0.6298(5). Fig. 12 shows the temperature derivative of L M 2 z evaluated at T = T c> , which is expected to scale as While not conclusive, the results are consistent with Ising-like criticality for smaller system sizes, crossing over to the true mean-field universality class for L > 25.
For the available system sizes, there is no evidence of the expected logarithmic corrections to scaling. We do not consider this to be strong evidence for their absence, however, since much larger systems are often required to observe logarithmic corrections. 33

V. LOWER-TEMPERATURE TRANSITION
The lower-temperature transition, at T = T c< , separates the ferromagnetic Coulomb phase from a conventional ferromagnet. Because the magnetization is nonzero on both sides, the spin-inversion symmetry of the Hamiltonian is immaterial, and the transition is in the same universality class as the saturation transition in an applied field. 14,22 This is a Kasteleyn transition, which exhibits anisotropic scaling in the directions parallel and perpendicular to the magnetization, with relative scaling exponent z = 2. The transition is consequently at its upper critical dimension, and so shows mean-field exponents with logarithmic corrections. 19,22 The Kasteleyn transition has the distinguishing characteristic that the magnetization is saturated in the low-temperature phase (in the thermodynamic limit), and decreases continuously but nonanalytically across the transition. 22,34 The magnetization is plotted in Fig. 13 for large systems near the lower-temperature transition, showing the development of a kink as the system size grows and indicating that the departure from saturation magnetization for T < T c< is a finite-size effect.
Although there is no symmetry breaking at the transition, the quantity 1 − M z /M sat can be identified as an order parameter, taking a nonzero value only on the high-temperature side. The critical theory for the transition 14,22 can be written using a U(1)-symmetric complex field ψ, in terms of which the order parameter is given by 1 − M z /M sat ∼ ψ * ψ. It follows that the scaling dimension of L(1 − M z /M sat ) vanishes, 19 and so a crossing point is expected when this quantity is plotted for different L. This crossing is shown in Fig. 14, enabling (T/p) c< = 3.15302(9) to be found and providing confirmation of the Kasteleyn universality class.
Previous MC simulations of spin ice in an applied field 19,22 have shown that logarithmic corrections are visible for L 100. In the present case, the further-neighbor interactions in  13. Magnetization near the lower-temperature transition at T c< , indicated with a vertical line, for large system sizes. As L increases, the magnetization approaches saturation below the transition, and a kink develops at T c< (compare also Fig. 4, for L = 24). In this case, the MC simulation is run starting from a state with saturated magnetization, M z = M sat ; the low temperatures and large system sizes ensure that ergodicity is broken and the order parameter takes a nonzero value. Order parameter for the lower-temperature transition, 1 − M z /M sat , multiplied by L and plotted versus temperature for various L (using the same symbols as Fig. 13). This quantity has vanishing scaling dimension at the Kasteleyn transition; a crossing is observed at (T/p) c< = 3.15302(9).
the model make the MC simulations more computationally demanding, and only mean-field behavior is observed at accessible system sizes.
The Kasteleyn transition occurs when the magnetization first deviates from saturation, which occurs through the appearance of strings of spins flipped relative to the fully polarized state. 22 One can therefore determine the exact transition temperature by considering the free energy of a single string, and finding the point where this becomes negative. 21,22 When V 4 = 0, such a string contributes free energy of ∆ f = 4p − T ln 2 per unit length, where the second term reflects the entropy associated with the possible paths. Including the four-spin interaction V 4 modifies this to ∆ f = 4p − T ln e 12V 4 /T + e 11V 4 /T , (11) because the two possible routes for the string (following a 011 chain or otherwise) have different energies. The Kasteleyn transition occurs when ∆ f = 0, giving (T/p) c< = 3.15343 for V 4 /T = 0.05. The discrepancy with the numerical results, which is small in absolute terms but several times the estimated statistical error, may result from logarithmic corrections to the scaling form for M z . Finally, it should be noted that, while both the higher-and lower-temperature transitions are at their upper critical dimensions, and hence have rational exponents, they otherwise have quite different properties. The Kasteleyn transition has anistropic scaling in directions parallel and perpendicular to the magnetization, while at the higher-temperature transition the system scales isotropically. A second distinction is that the magnetization is the critical field for the higher-temperature transition, while it is related to a bilinear of the critical field ψ for the Kasteleyn transition.

VI. DISCUSSION
This work has studied a nearest-neighbor model of spin ice with uniaxial distortion, which has a ferromagnetic phase at low temperature. Analytical arguments were used to show that an additional four-spin interaction can lead to an intermediate phase with coexisting ferromagnetic order and spin-liquid characteristics; the presence of this ferromagnetic Coulomb phase (FCP) has been established using Monte Carlo simulations.
Many features of the FCP are expected to occur more generally in spin liquids, both classical and quantum, with magnetic order. A clear experimental signature of an ordered Coulomb phase is coexistence of Bragg peaks, indicating magnetic order, with pinch points (see Fig. 10). On the theoretical side, a defining characteristic of spin-liquid phases is fractionalized excitations, such as magnetic monopoles in spin ice and spinons in quantum antiferromagnets; these remain deconfined across the transition into an ordered spin liquid (see Fig. 8). Finally, such phase transitions have conventional order parameters, but their critical properties are modified by coupling to the soft modes of the spin liquid (see Section IV).
The analysis here, including the Monte Carlo data, has treated the ice rule (see Section II A) as a strict constraint on configurations of the model. With a nonzero but small density of monopoles (i.e., defects in this constraint), as in the spin-ice compounds at low temperature, one expects most of the results to apply essentially unchanged: While the lowertemperature transition is immediately replaced by a crossover, this remains sharp for small monopole density. 22 The FCP is no longer qualitatively distinct from a conventional ferromagnet, but there can be a clear regime where the system is effectively described by a classical spin liquid with a small density of monopoles. 19 The higher-temperature transition remains, but is strictly in the Ising universality class at any nonzero monopole density; as in the case of the cubic dimer model, 35 however, the critical behavior is strongly affected by the presence of the unconventional critical point at zero monopole density.

ACKNOWLEDGMENTS
I am grateful to John Chalker, Ludovic Jaubert, and Michael Levin for helpful discussions. The Monte Carlo simulations used resources provided by the University of Nottingham High-Performance Computing Service.

Appendix: Quantum mapping
In this Appendix, we consider a model of quantum bosons in two spatial dimensions (2D), which shows closely analogous behavior to the model of spin ice discussed in the body of the paper. In fact, using the general mapping between classical statistical mechanics in 3D and quantum mechanics in 2D, which has previously been applied to phase transitions from CSL phases, 14,36 one expects the universal features of the phases and transitions to be equivalent in these two models.
The nearest-neighbor model for spin ice, H nn , can be mapped to a system of hard-core lattice bosons, with spin-reversal symmetry replaced by particle-hole symmetry. The Coulomb phase of the spin model is equivalent to a superfluid, 14,22 while the saturated ferromagnet with M z = ±M sat is equivalent to the vacuum and fully-occupied states of the bosonic model, which spontaneously break particlehole symmetry. The strings of flipped spins that proliferate at the transition (see Section II B) map to boson world lines (trajectories in space and imaginary time). An equivalent bosonic model to the Hamiltonian H nn , displaying infinite-order multicriticality at the transition between these two phases, is given by 15 i are annihilation and number operators for hard-core bosons. The first term represents tunneling t between neighboring sites i j , while the second is an attractive interaction of strength V > 0 between nearest-neighbor bosons, written in a manifestly particlehole-symmetric form.
Because H 0 conserves particle number, the Hilbert space can be divided into sectors of fixed density ρ = n i ; let E gs (ρ) be the ground-state energy in each. For t > V, the overall ground state occurs in the sector with ρ = 1 2 , and the system is a particle-hole-symmetric superfluid. For t < V, E gs is instead minimized by either the vacuum, ρ = 0, or the fully-occupied state, ρ = 1. At t = V, the model has a Rokhsar-Kivelson (RK) point, 37 at which H 0 can be written as a projector, at which E gs is independent of N. For V < t the ground state of the system occurs for half filling, N = 8. For V > t there are two generate ground states, with boson density zero and one (N = 0 and N = 16 respectively), which spontaneously break particle-hole symmetry. The insets show E gs versus N for fixed values of V/t, indicated by the arrows. and the ground state in each sector, an equal-amplitude superposition of all configurations, 15 has E gs (ρ) = 0. As illustrated in Fig. 15, which shows results of exact diagonalization (ED) on a small system, this leads to a transition with identical properties to the ordering transition of spin ice under uniaxial pressure, with F(M z ) replaced by E gs (ρ). (The exact equivalence is established by noting that the transfer matrix for the classical problem can be written as a projector at the transition. 21 ) Following similar logic to Section II C, one expects that quartic interactions between bosons should open up an intermediate phase with density changing continuously between ρ = 1 2 and ρ = 0 or 1. The precise form required is again unclear from these general arguments, but ED results, shown in Fig. 16, indicate that it suffices to add a (particle-holesymmetrized) four-body repulsion where the sum is over sites i jkl around a square plaquette. In this case, there are two continuous transitions, with density changing from n i = 1 2 to 0 < | n i − 1 2 | < 1 2 and then to | n i − 1 2 | = 1 2 , as V/t is increased. (The nature of the transitions is clear even for the small system sizes accessible in ED, because the order parameter commutes with the Hamiltonian.) The transition into the vacuum or fully-occupied (vacuum of holes) state is described by the standard critical theory for the . As in the case with V 4 = 0, the ground state is at half filling for V t and has density either zero or one for V t. In between, there is a regime where the minimum of E gs crosses over between the extremes, stepping through each intermediate density in turn. In the thermodynamic limit, this is expected to become a phase with continuously varying density, separated from small-and large-V/t phases by continuous quantum phase transitions. vacuum transition of bosons, 38 while the transition at lower V/t involves spontaneous breaking of particle-hole symmetry within the superfluid, and is described by the critical theory of Section IV A. In cases where the total particle number is fixed, the latter transition would lead to phase separation into regions with differing densities.
The additional interaction H 4s in the classical spin model, defined in Eq. (3), may be viewed as the equivalent of H 4b . To see this, recall that strings of flipped spins are equivalent to bosons, and that these occur at low density near the transition to saturation (bosonic vacuum). Two strings passing through a tetrahedron t change its total spin from + 4 √ 3ẑ to − 4 √ 3ẑ , and so the interaction H 4s (with V 4 > 0) amounts to an energy penalty when four strings are in close proximity (passing through two tetrahedra on opposite sides of the same hexagon).