Ferromagnetic Coulomb phase in classical spin ice
FFerromagnetic Coulomb phase in classical spin ice
Stephen Powell
School of Physics and Astronomy, The University of Nottingham, Nottingham, NG7 2RD, United Kingdom
Spin ice is a frustrated magnetic system that at low temperatures exhibits a Coulomb phase, a classical spinliquid with topological order and deconfined excitations. This work establishes the presence of a Coulombphase with coexisting ferromagnetic order in a microscopic model of classical spin ice subject to uniaxial latticedistortion. General theoretical arguments are presented for the presence of such a phase, and its existence isconfirmed using Monte Carlo results. This example is used to illustrate generic properties of spin liquids withmagnetic order, including deconfinement of monopoles, signatures in the neutron-scattering structure factor, andcritical behavior at phase transitions. An analogous phase, a superfluid with spontaneously broken particle–holesymmetry, 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 freedomexhibit strong local correlations, despite the persistence oflarge fluctuations, of either quantum mechanical or thermalorigin. They occur at low temperature in certain frustratedsystems, where interactions are large compared to thermalfluctuations, but mutual competition between the interactionsprevents formation of a rigidly ordered configuration. Spinliquids have been of theoretical interest for several decades, but evidence for their existence in physical systems, andeven microscopic models, is considerably more recent.While spin liquids are often distinguished from conven-tional low-temperature phases, such as ferromagnets, by thefact that they lack magnetic order, their defining characteris-tics go beyond the mere absence of conventional order. A pre-cise definition of a quantum spin liquid (QSL) can be phrasedin terms of long-range entanglement, while the Coulombphase, the classical spin liquid (CSL) that is of primary in-terest here, can be defined through deconfinement of fraction-alized “monopole” excitations. Experimental evidence existsfor a Coulomb phase in the spin-ice compounds, which can betreated as classical at relevant temperatures. These definitions provide positive characterizations forQSL and CSL phases, and also make clear the possibility ofa magnetically ordered spin liquid, in which spin-liquid phe-nomena coexist with conventional symmetry-breaking order.Some examples of such phases have been reported in the the-oretical literature: Mean-field studies of quantum spin ice identified an ordered QSL, referred to as a “Coulomb ferro-magnet”, although quantum Monte Carlo simulations have notrevealed such a phase. Recent work has also demonstratedthe possibility of antiferromagnetic order coexisting with aCSL.The compatibility of magnetic order and spin-liquid phe-nomenology also allows for the existence of phase transitionsbetween ordered and disordered spin liquids. One might an-ticipate novel critical behavior at such transitions, since it isknown that transitions from spin liquids into conventional or-dered phases can transcend the usual Landau paradigm. This work demonstrates that a ferromagnetic Coulombphase can occur in a model of classical spin ice, and providesa detailed study of this phase and the associated transitions. Theoretical arguments, including mapping to a related quan-tum model, are used to show that such a phase exists and that itcan be reached through a continuous transition from the para-magnetic Coulomb phase. We present Monte Carlo (MC) re-sults that confirm both of these statements, and illustrate thegeneric properties of ordered spin liquids, including the struc-ture factor for elastic neutron scattering.We also consider the critical behavior at the ordering tran-sition and predict that, despite the Ising nature of the orderparameter and the presence of only short-range interactionsin the microscopic model, the transition should belong in themean-field universality class, as a result of coupling to the ef-fective gauge-field fluctuations of the spin liquid. While thenumerical results are consistent with this prediction, largersystem sizes would be required for a definitive confirmation ofthe universality class. This phase transition provides anotherinteresting example of the diversity of critical phenomena thatexists in the neighborhood of spin-liquid phases.
Outline
In Section II, the model of spin ice is introduced, and achoice of perturbations that lead to a ferromagnetic Coulombphase is motivated. The basic structure of the phase diagramis then illustrated using MC results, showing the appearanceof such a phase at intermediate temperatures for certain val-ues of the parameters. In Section III, the phase in question isstudied in detail, to confirm that it has nonzero magnetizationwhile simultaneously exhibiting the characteristic features ofa Coulomb phase. Sections IV and V consider in turn the crit-ical behavior at the higher- and lower-temperature transitionsout of this ferromagnetic Coulomb phase.We conclude in Section VI, by summarizing the featuresthat are expected to be generic to ordered spin liquids, bothquantum and classical, and discuss briefly the e ff ect of anonzero density of magnetic monopoles. In the Appendix, theclassical–quantum mapping developed in Ref. 14 is appliedto this system, and the resulting quantum model is related toa problem of hard-core quantum bosons studied by Rokhsarand Kotliar. a r X i v : . [ c ond - m a t . s t a t - m ec h ] M a r II. MODEL AND PHASE STRUCTUREA. Nearest-neighbor model of spin ice
The spin-ice materials Ho Ti O and Dy Ti O are welldescribed by a model of classical spins S i on the sites i of a py-rochlore lattice, a network of corner-sharing tetrahedra. Eachspin is subject to a strong easy-axis anisotropy constraining itto point parallel or antiparallel to the local (cid:104) (cid:105) axis joiningthe centres of adjacent tetrahedra, S i = ± ˆ n i . Including onlynearest-neighbor interactions, the Hamiltonian can be writtenas H nn = − (cid:88) (cid:104) i j (cid:105) J i j S i · S j , (1)where J i j is a ferromagnetic coupling between nearest-neighbor sites (cid:104) i j (cid:105) 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 eachtetrahedron, two point in and two point out. The latter condi-tion is referred to as the “ice rule” and selects a set of statesthat is degenerate in the nearest-neighbor model and whosenumber grows exponentially with the number of spins. Whilea more realistic microscopic model than H nn would also in-clude dipolar interactions, their e ff ect is primarily to renor-malize J , with only a small splitting of the ice-rule states. We will mostly concentrate on the limit where the ice ruleis enforced as a constraint, represented by Eq. (1) with tem-perature T (cid:28) J . Assuming ergodicity within the ice-rulemanifold, the system in this limit exhibits a Coulomb phase,in which the spins are disordered but highly correlated. Re-placing the spins S i by a continuous vector field B ( r ) and theice rule by ∇ · B = ff ective coarse-graineddescription for this phase. A quadratic action for the “mag-netic field” B correctly describes the long-wavelength neutronscattering at low temperature in spin ice, and predicts thatmonopoles in B , corresponding to single tetrahedra where theice rule is broken, are deconfined. Much of the physics is infact qualitatively unaltered by a small density of such defects(see Section VI), and their e ff ects on the critical propertiescan be understood by treating monopole fugacity as a relevantperturbation (in the renormalization-group sense). An important property of the ice-rule states for presentpurposes is that they obey a topological constraint on themagnetization:
Starting from any ice-rule state and flippinga spin S i breaks the ice rule on the two tetrahedra shared bysite i . The only updates that connect configurations within theice-rule manifold are those that involve flipping a set of spinsaligned head-to-tail along a closed loop. Any such update fora contractible loop preserves the magnetization density, M = N s (cid:88) i S i , (2)where N s = (cid:80) i E t = - J - p E t = - J - p E t = - J + p FIG. 1. Three configurations of a single tetrahedron, and their en-ergy E t in the nearest-neighbor Hamiltonian H nn , Eq. (1), with J ij given in Section II B. Pairs of spins situated in the same horizon-tal plane, indicated with dashed (red) lines, have reduced coupling J ij = J − p , while others have J ij = J . All three configurations obeythe ice rule, having two spins pointing in and two pointing out. Thefirst two are lowest-energy configurations for a single tetrahedron(since the antiferromagnetically aligned pairs are those with reducedcoupling), while the one on the right is one of the remaining four ice-rule configurations whose energy is higher by 4 p . Excitations abovethe ground state are described by strings of spins flipped relative toa fully polarized configuration, and increase the energy by 4 p pertetrahedron. When a pair of strings pass through the same tetrahe-dron, all four spins are inverted and the energy is again minimized;the strings therefore feel an attractive interaction. with the same magnetization therefore constitute “topolog-ical sectors”, disconnected by local updates. This topo-logical conservation law is broken by a nonzero density ofmonopoles, but remains approximately valid, and conceptu-ally useful, at low temperatures.Nonzero magnetic susceptibility χ requires that the systemfluctuates between di ff erent sectors; in the thermodynamiclimit, one can therefore distinguish “incompressible” phaseswith χ = χ > B. Uniaxial distortion
To split the energy of the six ice-rule states on a giventetrahedron requires breaking the cubic symmetry of the py-rochlore lattice. Following Ref. 21, we consider an explicituniaxial symmetry breaking, with J i j = J − p ( J > p > J i j = J for all others. (Such a perturbation could bee ff ected in experiment by application of uniaxial pressure. )As illustrated in Fig. 1, the result is to favor the two configura-tions where the total (vector) spin of the tetrahedron is alongthe [001] axis, whose energy is reduced by 4 p compared tothe other four. In contrast to the case of an applied field, an Ising symmetry remains; there are two degenerate lowest-energy states, with all spins on all tetrahedra maximally polar-ized, 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 23for more details. For p (cid:28) T (cid:28) J , the system remainsin the Coulomb phase, while below a critical temperature T c = p (ln 2) − it becomes a fully polarized ferromagnet. Atthe transition, the up–down symmetry is broken and the mag-netization along the [001] direction, M z , becomes nonzero.While an Ising order parameter can naturally be defined, thetransition in the ice-rules limit has quite di ff erent propertiesfrom the standard Ising universality class. Starting from ei-ther of the fully polarized states, the only closed loops are“strings” spanning the system in the [001] direction, whichcost energy proportional to the (linear) system size L . Thetransition occurs when the entropy of a single string, also ∝ L ,outweighs the energy, and so its free energy changes from pos-itive to negative; the string density then increases from zero tononzero.As a consequence, the system on the lower-temperatureside of the transition is fully polarized, with zero string den-sity, as in the related case of an applied field. A crucialdistinction in this case is that two strings feel an e ff ective at-traction when sharing a tetrahedron (see Fig. 1). At the criticalpoint, this exactly balances the entropy cost of the excludedvolume due to the hard-core interactions between strings. Infact, as Jaubert et al. have shown, the free energy at T = T c as a function of string density is precisely constant (in the ther-modynamic limit). Because each string consists of a fixednumber of flipped spins relative to the starting configuration,this implies that the free energy is independent of magneti-zation. As the temperature increases through the transition,the global minimum of F ( M z ), which can be interpreted as aLandau function, jumps from M z = ± M sat to M z =
0. (The re-sulting discontinuity in the magnetization is illustrated belowin Fig. 3.) Since all coe ffi cients in the Landau free energy van-ish at the transition, this has been referred to as “infinite-ordermulticriticality”. C. Additional interactions
Given the magnetization-independent free energy at thetransition, it is clear that any perturbation that produces a pos-itive fourth-order coe ffi cient in the Landau function shouldlead to an intermediate phase with 0 < | M z | < M sat . Whilethis argument does not provide a prescription for construct-ing appropriate perturbations, one expects on general groundsthat a su ffi ciently long-ranged four-spin interaction will havethis e ff ect. (As will also be demonstrated, a quartic coe ffi cientwith opposite sign should lead to a first-order transition.)As we detail in the following, MC results in fact demon-strate that it is su ffi cient to add a four-spin interaction actingbetween tetrahedra on opposite sides of a hexagonal loop, asillustrated in Fig. 2. The perturbation used throughout thiswork can be written explicitly as H = V (cid:88) { tt (cid:48) } [ Θ + ( S t , S t (cid:48) ) + Θ − ( S t , S t (cid:48) )] , (3)where Θ ± ( S , S (cid:48) ) = S = S (cid:48) = ± √ ˆ z FIG. 2. Illustration of the four-spin interaction H added to themodel to stabilize the ferromagnetic Coulomb phase. The arrowsrepresent spins on the sites of a pyrochlore lattice, a network ofcorner-sharing tetrahedra. (This configuration obeys the ice rule,with two spins pointing into and two pointing out of each tetrahe-dron.) The additional interaction couples pairs of tetrahedra at oppo-site sides of hexagons; one such pair and its hexagon are highlighted. and S t ≡ (cid:80) i ∈ t S i is the total (vector) spin on tetrahedron t . Thesum in Eq. (3) is over pairs of tetrahedra { tt (cid:48) } across a hexagon(see Fig. 2), and the summand is one if both tetrahedra haveall spins polarized in the same vertical direction, and zero oth-erwise. (Note that, while this expression apparently involveseight spins, it is equivalent to a four-spin interaction underprojection into the ice-rule states. This form of the interactionis partly motivated by the quantum mapping, described in theAppendix.)Regarding the choice of H , it is not the goal of thiswork to classify the various types of interactions according towhether they produce a ferromagnetic Coulomb phase, and weare not aware of a general argument that would allow for sucha classification. (The search for appropriate interactions isin any case better informed by experimental evidence aboutwhich interactions occur in particular materials.) Rather, thegoal here is to study a particular case where such a phase isknown to exist, and elucidate those properties of the phaseand its transitions that are expected to be universal, or at leastqualitatively generic.Plots of the magnetization as a function of temperature, for V positive, negative, and zero, are shown in Figs. 3 and 4.These results were produced using MC simulations based ona directed-loop algorithm. The lattice consists of L × L × L cubic unit cells, each containing 4 tetrahedra of each orienta-tion, and hence 16 spins. For V ≤
0, a step is observed in themagnetization, from essentially fully saturated, with small de-viations due to finite-size e ff ects, to zero within error bars. m a gn e ti za ti on X M z \ (cid:144) M s a t T (cid:144) p FIG. 3. Magnetization versus temperature, for fixed V / T = V / T = − .
01 (right), and L =
24 ( N s = L (cid:39) × spins).In both cases, there is a jump from saturation to zero magnetization,at a transition temperature indicated with a vertical line. For eachtemperature, the spontaneous magnetization is found by applying aweak field along the z direction and extrapolating to zero field. m a gn e ti za ti on X M z \ (cid:144) M s a t T (cid:144) p FIG. 4. Magnetization versus temperature, for fixed V / T = .
05 andsystem size L =
24. The vertical lines at T / p (cid:39) .
15 and 3 .
35 indi-cate positions of phase transitions, determined as described in Sec-tions V and IV respectively. Below the lower-temperature transition,the magnetization remains at its saturation value, apart from smallfinite-size corrections, while above the higher-temperature transition,it vanishes. The intermediate phase is a ferromagnet with nonzeroand continuously varying magnetization.
This step is accompanied by a single peak in the specific heat(not shown), whose height grows with system volume, indi-cating a single first-order transition.By contrast, when V > T is lowered. The high-temperature phase is paramagnetic, with M = , and is theusual Coulomb phase observed at T (cid:28) J in spin ice. Themagnetization first becomes nonzero at T c > before reaching itssaturation value at T c < . While the system is ferromagnetic forall 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 theenergy (proportional to the specific heat) in this case displaystwo peaks, both at most weakly diverging with L , consistentwith a pair of continuous transitions.Fig. 6 shows histograms of the energy and magnetizationfor L = V / T = .
05, and T / p = .
31, near the higher- e n e r gyv a r i a n ce T (cid:144) p L = FIG. 5. Variance of energy, ( (cid:104) E (cid:105) − (cid:104) E (cid:105) ) / N s , versus temperaturefor fixed V / T = .
05 and various system sizes. The vertical linesindicate the positions of phase transitions in the thermodynamic limit(determined by other means). The double-peak structure, with peakheights at most weakly diverging with L , is consistent with a pair ofcontinuous transitions. temperature peak of the energy variance. The unimodal struc-ture of the energy distribution confirms the continuous natureof the transition, while the two peaks of the magnetization in-dicate that this is a symmetry-breaking transition into a statewith nonzero but unsaturated magnetization. This should becontrasted with the case of V =
0, where the magnetizationhistogram is flat at the transition. Fig. 7 shows the case of V <
0, where the transition is of first order, with a bimodalstructure in the energy and coexisting peaks in the magnetiza-tion distribution, at both M z = M z = ± M sat . III. INTERMEDIATE PHASE
Having established the presence of a pair of phase transi-tions when V >
0, we now turn to the intermediate phasein the temperature range T c < < T < T c > . It will be arguedthat this phase shares the essential spin-liquid features of theCoulomb phase above T c > , but is distinguished by a nonzerospontaneous magnetization.The presence of a nonzero but unsaturated magnetizationin the intermediate phase is evident from Figs. 4 and 6. Con-tinuously changing magnetization implies that the magneticsusceptibility is nonzero, and hence that there are fluctua-tions between di ff erent topological (magnetization) sectors.This fact alone distinguishes the intermediate phase from thelow-temperature saturated ferromagnet, where the flux sti ff -ness vanishes in the thermodynamic limit, and there are notopological-sector fluctuations. Two phenomena that are characteristic of the Coulombphase 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 di ff ers) cor- - - - - - E (cid:144) N s - - - M z (cid:144) M sat FIG. 6. Energy and magnetization histograms for L = V / T = .
05, and T / p = .
31 (near the higher-temperature transition). The unimodalenergy distribution is indicative of a continuous transition, while the two peaks in the magnetization distribution show that this transition isassociated with magnetic ordering and breaking of spin-reversal symmetry. - - - E (cid:144) N s - - M z (cid:144) M sat FIG. 7. Energy and magnetization histograms for L = V / T = − .
01, and T / p = .
76. In this case, the transition is strongly first-order,as indicated by the bimodal energy distribution, and occurs directly between the saturated ferromagnet ( M z = ± M sat ) and the paramagnet( M z = L the competing states are metastable.) responds to a monopole in the continuous vector field B ( r ).Such defects are rare for T (cid:28) J , and, at least as a first approx-imation, we treat the density of thermally excited monopolesas vanishing.It is useful to consider, however, the introduction of a sin-gle pair of oppositely charged monopoles into an otherwisedefect-free background. The e ff ective interaction between thepair, induced by the fluctuations of the surrounding spins, al-lows one to distinguish spin-liquid phases from others suchas the saturated ferromagnet. In the Coulomb phase, themonopoles are subject to an e ff ective Coulomb interaction,with a finite limit for large separation. The saturated ferro-magnet is, by contrast, a confining phase, in which the freeenergy of a pair of monopoles grows without bound as theirseparation increases. 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 pairof monopoles of opposite charge at r ± . (More explicitly,the ensemble is constrained so that all tetrahedra obey theice rule, apart from those at r ± , where three spins point outand one points in, and vice versa.) This function, which isrelated to the e ff ective 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 confinedphase, it instead decays exponentially to zero.In a finite system, these asymptotic behaviors are observedonly for separations much less than the system size L . Finite-size e ff ects can be controlled by fixing the ratio | r + − r − | / L and observing the scaling with L . Fig. 8 shows the ratio of themonopole distribution function calculated at R max , the largestdisplacement possible for L cubic unit cells ( L even) withperiodic boundaries, and at R max, z , the maximum separationalong the z direction ( | R max | = √ | R max, z | ). The ratio ap-proaches unity with increasing system size for all T > T c < ,while it decays to zero below T c < , indicating confinement.No qualitative di ff erence is seen when crossing the higher-temperature phase boundary at T c > , demonstrating that theintermediate phase, in common with the standard Coulombphase above T c > , exhibits deconfinement of monopoles.The form of the e ff ective interaction U m ( r ) = − T ln G m ( r )is determined by the approach of G m ( r ) to its asymptotic valuefor large separation r . Fig. 9 shows U m for temperatureswithin the intermediate phase and above T c > . In both cases,the interaction is anisotropic, because the spatial symmetryis reduced by the applied pressure (and H ). Up to finite-size e ff ects, the interaction can be fit to the Coulomb form, (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231)(cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231)(cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231)(cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231)(cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231)(cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231)(cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231) (cid:231)
10 12 14 16 18 20 220.00.20.40.60.81.0 system size L G m (cid:72) R m a x (cid:76) (cid:144) G m (cid:72) R m a x , z (cid:76) T (cid:144) p M z (cid:144) M s a t FIG. 8. Ratio of monopole distribution function G m evaluated at themaximum displacement in the lattice, R max , and at the maximum dis-placement along the z direction, R max, z . The curves (from bottom totop) have T / p corresponding to the dashed lines (from left to right) inthe inset, which reproduces the magnetization curve of Fig. 4. The ra-tio approaches unity with increasing system size for all temperaturesexcept the lowest, which is within the low-temperature confiningphase. No qualitative change is seen across the higher-temperaturetransition, indicating that the two higher-temperature phases both ex-hibit deconfinement of monopoles. Explicitly, the curves have, frombottom to top, T / p = .
125 (black), 3 .
175 (red), 3 .
226 (orange),3 .
279 (yellow), 3 .
333 (green), 3 .
390 (light blue), and 3 .
448 (darkblue). In all cases, V / T = .
05 is fixed. ∝ / | r | , confirming the identification of the intermediate phaseas a Coulomb phase. The e ff ective interaction is stronger par-allel to the pressure axis at both temperatures, with largeranisotropy at the lower temperature. B. Neutron-scattering structure factor
The most direct experimental signature of the Coulombphase is the presence of “pinch points” in the neutron-scattering structure factor, which reflect the algebraic (dipolar)correlations between the spins.
These features are clear-est in the spin-flip component of polarized neutron-scatteringdata with incident polarization along [1¯10]. The correspond-ing structure factor is S SF ( Q ) = ˆ η µ ˆ η ν S µν ( Q ) , (5)where S µν ( Q ) is the Fourier transform of the two-spin corre-lation function (cid:104) S µ ( r ) S ν ( r (cid:48) ) (cid:105) andˆ η = Q × P | Q || P | (6)is a unit vector orthogonal to both the wavevector Q and theincident neutron polarization P .This structure factor is shown in Fig. 10, for Q in the ( hh (cid:96) )plane and P along [1¯10]. Pinch points are visible for all T > T c < , with no qualitative change at T c > , showing that the T (cid:144) p (cid:61) (cid:72)(cid:254)(cid:76) T (cid:144) p (cid:61) (cid:72) (cid:166) (cid:76) T (cid:144) p (cid:61) (cid:72)(cid:254)(cid:76) T (cid:144) p (cid:61) (cid:72) (cid:166) (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) (cid:160) r (cid:164) (cid:144) L U m (cid:144) T FIG. 9. E ff ective (dimensionless) interaction U m ( r ) / T = − ln G m ( r )between monopoles, for fixed system size L =
32. For each temper-ature T , the interaction is plotted for separations r parallel ( (cid:107) ) andperpendicular ( ⊥ ) to the axis of the applied pressure, and the zero ofinteraction is chosen as U m ( R max ) =
0. The lines show least-squarefits to the Coulomb form a − b / | r | in the region 0 . < | r | / L < . ff erent parameters a and b for each case. There are deviationsfrom the fit at large separation, because of finite-size e ff ects, and atsmall separation, because of lattice-scale e ff ects and because of thefinite range of the additional interactions. (The parameter b is givenby b (cid:107) = . b ⊥ = . T / p = . b (cid:107) = . b ⊥ = . T / p = . dipolar correlations of the Coulomb phase remain until thelower-temperature transition. In the intermediate phase, theycoexist with Bragg peaks at certain reciprocal lattice vectors,indicating ferromagnetic ordering. The di ff use scattering iscompletely suppressed in the low-temperature saturated ferro-magnet, 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, andthat it is connected to the high- and low-temperature phases bycontinuous transitions. In this section and the following, thecritical properties of these two transitions will be addressed inturn, using analytical arguments supported by numerical re-sults. A. Critical theory
Near the transition, at T = T c > , between the paramagneticand ferromagnetic Coulomb phases, the magnetization is farfrom saturation and so the discrete nature of the spins is pre-sumably not important. Replacing the discrete spins by a con-tinuous vector field B ( r ), the partition function can be ex-pressed as Z = (cid:90) D B δ ( ∇ · B ) exp − (cid:90) d r (cid:32) κ | B | − α B z (cid:33) . (7) - - - - T (cid:144) p = ¥ - - - - T (cid:144) p = - - - - T (cid:144) p = - - - - T (cid:144) p = - - - - T (cid:144) p = FIG. 10. Structure factor for (spin-flip) polarized neutron scattering S SF ( Q ), defined in Eq. (5), for scattering wavevector Q in the ( hh (cid:96) )plane and incident polarization P along [1¯10]. The first three plots are for temperatures above T c > , the fourth is in the intermediate phase, T c < < T < T c > , and the last is below T c < . Pinch points, characteristic of the dipolar correlations of the Coulomb phase, are visible at alltemperatures but the lowest. In the intermediate phase, there are also Bragg peaks at certain reciprocal lattice vectors, indicating spontaneousmagnetization. The system size is L =
16 and all plots have V / T = .
05. (Wavevectors are measured in units corresponding to the conventionalcubic unit cell.)
The coe ffi cient κ is the flux sti ff ness in directions transverse tothe applied pressure, while α > ff ect of theuniaxial pressure, enhancing fluctuations along the z direction.(Higher-order terms have been omitted.)Using a Hubbard–Stratonovich field Φ to decouple theanisotropy term, this can be replaced by Z ∝ (cid:90) D Φ (cid:90) D B δ ( ∇ · B ) × exp − (cid:90) d r (cid:32) α Φ + κ | B | + Φ B z (cid:33) . (8)The real scalar field Φ has Ising symmetry and provides anorder parameter for the transition, taking a nonzero value inthe ferromagnetic phase. Integrating out B induces dipolar in-teractions for Φ , giving an e ff ective description that is equiv-alent to that of Ising spins with dipolar couplings. A similarconnection between the dipolar correlations in the Coulombphase and e ff ective dipolar interactions at a critical point hasbeen noted in Ref. 28.The 3D Ising transition with dipolar interactions is atits upper critical dimension, and so shows mean-field criti-cal exponents with logarithmic corrections. In particular, the specific-heat, order-parameter, susceptibility, and correlation-length exponents take the values α = β = , γ =
1, and ν = , respectively.It should be noted that scaling remains isotropic at thistransition, in the sense that all spatial dimensions scale withthe same exponents. For example, the correlation lengths inthe directions parallel and perpendicular to the applied pres-sure diverge with the same exponent ν , though with di ff erent(nonuniversal) prefactors. This is in contrast to the anisotropicscaling at the lower-temperature transition (see Section V). B. Numerical results
To determine the critical temperature T c > and find valuesof the exponents, it is convenient to identify a quantity whosescaling dimension vanishes, for which curves with di ff erent L coincide at the transition. While the Binder cumulant of themagnetization provides such a quantity for this ordering tran-sition, it is di ffi cult to calculate accurately, as a result of thetopological constraints on the magnetization, which suppressfluctuations of the latter.We instead consider the quantity L (cid:104) M z (cid:105) , which, as a result L Y M z ] T (cid:144) p
16 18 2022 24 2628 30 32 L = FIG. 11. Plot of L (cid:104) M z (cid:105) versus temperature near the higher-temperature transition, for various system sizes L . (In each case V / T = .
05 is fixed.) This quantity has vanishing scaling dimen-sion for the mean-field universality class, consistent with the crossingpoint for large L , at ( T / p ) c > = . of the scaling form (cid:104) M z (cid:105) ≈ L − d + γ/ν Ψ (cid:32) L /ν T − T c > T c > (cid:33) , (9)where Ψ is a universal function, is expected to have zero scal-ing dimension for this transition. (This quantity is equal, upto powers of L , to the flux sti ff ness Υ , which is not expectedto have vanishing scaling dimension at a transition betweentwo spin liquids.) As shown in Fig. 11, L (cid:104) M z (cid:105) plotted as afunction of T / p indeed has a crossing point for large systemsizes. Using the crossings for successive L values, we estimate( T / p ) c = . V / T = . d − γ/ν = . ν , which, for the Ising class, takesthe value ν = . L (cid:104) M z (cid:105) evaluated at T = T c > , which is expectedto scale as dd T L (cid:104) M z (cid:105) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T = T c > ∼ L − d + γ/ν + + /ν = L /ν . (10)While not conclusive, the results are consistent with Ising-likecriticality for smaller system sizes, crossing over to the truemean-field universality class for L > è è è è è è l n £ d d TL Y M z ] §
22 24 26 28 30 32 L FIG. 12. Log–log plot of the temperature derivative (in arbitraryunits) of L (cid:104) M z (cid:105) , evaluated at T = T c > , versus system size L . Theslope gives the reciprocal of the correlation-length exponent ν . The(blue) solid line has ν = , as expected for the mean-field universal-ity class, while the (purple) dashed line has the Ising value ν = . V. LOWER-TEMPERATURE TRANSITION
The lower-temperature transition, at T = T c < , separatesthe ferromagnetic Coulomb phase from a conventional ferro-magnet. 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 sat-uration transition in an applied field. This is a Kasteleyntransition, which exhibits anisotropic scaling in the directionsparallel and perpendicular to the magnetization, with relativescaling exponent z =
2. The transition is consequently at itsupper critical dimension, and so shows mean-field exponentswith logarithmic corrections.
The Kasteleyn transition has the distinguishing characteris-tic that the magnetization is saturated in the low-temperaturephase (in the thermodynamic limit), and decreases continu-ously but nonanalytically across the transition.
The mag-netization is plotted in Fig. 13 for large systems near thelower-temperature transition, showing the development of akink as the system size grows and indicating that the depar-ture from saturation magnetization for T < T c < is a finite-sizee ff ect.Although there is no symmetry breaking at the transition,the quantity 1 − (cid:104) M z (cid:105) / M sat can be identified as an order pa-rameter, taking a nonzero value only on the high-temperatureside. The critical theory for the transition can be writtenusing a U(1)-symmetric complex field ψ , in terms of which theorder parameter is given by 1 − (cid:104) M z (cid:105) / M sat ∼ ψ ∗ ψ . It followsthat the scaling dimension of L (1 − (cid:104) M z (cid:105) / M sat ) vanishes, andso a crossing point is expected when this quantity is plottedfor di ff erent L . This crossing is shown in Fig. 14, enabling( T / p ) c < = . have shown that logarithmic corrections are visible for L (cid:38) X M z \ (cid:144) M s a t T (cid:144) p
30 40 50 60 L = FIG. 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 akink develops at T c < (compare also Fig. 4, for L = (cid:104) M z (cid:105) = M sat ; the low temperatures and large systemsizes ensure that ergodicity is broken and the order parameter takes anonzero value. L H - X M z \ (cid:144) M s a t L T (cid:144) p FIG. 14. Order parameter for the lower-temperature transition, 1 −(cid:104) M z (cid:105) / M sat , multiplied by L and plotted versus temperature for various L (using the same symbols as Fig. 13). This quantity has vanishingscaling dimension at the Kasteleyn transition; a crossing is observedat ( T / p ) c < = . the model make the MC simulations more computationallydemanding, and only mean-field behavior is observed at ac-cessible system sizes.The Kasteleyn transition occurs when the magnetizationfirst deviates from saturation, which occurs through the ap-pearance of strings of spins flipped relative to the fully po-larized state. One can therefore determine the exact transi-tion temperature by considering the free energy of a singlestring, and finding the point where this becomes negative.
When V =
0, such a string contributes free energy of ∆ f = p − T ln 2 per unit length, where the second term reflects theentropy associated with the possible paths. Including the four- spin interaction V modifies this to ∆ f = p − T ln (cid:16) e V / T + e V / T (cid:17) , (11)because the two possible routes for the string (following a (cid:104) (cid:105) chain or otherwise) have di ff erent energies. The Kaste-leyn transition occurs when ∆ f =
0, giving ( T / p ) c < = . V / T = .
05. The discrepancy with the numericalresults, which is small in absolute terms but several times theestimated statistical error, may result from logarithmic correc-tions to the scaling form for (cid:104) M z (cid:105) .Finally, it should be noted that, while both the higher- andlower-temperature transitions are at their upper critical di-mensions, and hence have rational exponents, they otherwisehave quite di ff erent properties. The Kasteleyn transition hasanistropic scaling in directions parallel and perpendicular tothe magnetization, while at the higher-temperature transitionthe system scales isotropically. A second distinction is that themagnetization is the critical field for the higher-temperaturetransition, 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 icewith uniaxial distortion, which has a ferromagnetic phase atlow temperature. Analytical arguments were used to show thatan additional four-spin interaction can lead to an intermedi-ate phase with coexisting ferromagnetic order and spin-liquidcharacteristics; the presence of this ferromagnetic Coulombphase (FCP) has been established using Monte Carlo simula-tions.Many features of the FCP are expected to occur more gener-ally in spin liquids, both classical and quantum, with magneticorder. A clear experimental signature of an ordered Coulombphase is coexistence of Bragg peaks, indicating magnetic or-der, with pinch points (see Fig. 10). On the theoretical side,a defining characteristic of spin-liquid phases is fractional-ized excitations, such as magnetic monopoles in spin ice andspinons in quantum antiferromagnets; these remain decon-fined across the transition into an ordered spin liquid (seeFig. 8). Finally, such phase transitions have conventional or-der parameters, but their critical properties are modified bycoupling to the soft modes of the spin liquid (see Section IV).The analysis here, including the Monte Carlo data, hastreated the ice rule (see Section II A) as a strict constraint onconfigurations of the model. With a nonzero but small den-sity of monopoles (i.e., defects in this constraint), as in thespin-ice compounds at low temperature, one expects most ofthe results to apply essentially unchanged: While the lower-temperature transition is immediately replaced by a crossover,this remains sharp for small monopole density. The FCP isno longer qualitatively distinct from a conventional ferromag-net, but there can be a clear regime where the system is e ff ec-tively described by a classical spin liquid with a small densityof monopoles. The higher-temperature transition remains,but is strictly in the Ising universality class at any nonzero0monopole density; as in the case of the cubic dimer model, however, the critical behavior is strongly a ff ected by the pres-ence of the unconventional critical point at zero monopoledensity. ACKNOWLEDGMENTS
I am grateful to John Chalker, Ludovic Jaubert, andMichael Levin for helpful discussions. The Monte Carlo sim-ulations used resources provided by the University of Notting-ham High-Performance Computing Service.
Appendix: Quantum mapping
In this Appendix, we consider a model of quantum bosonsin two spatial dimensions (2D), which shows closely analo-gous behavior to the model of spin ice discussed in the bodyof the paper. In fact, using the general mapping between clas-sical statistical mechanics in 3D and quantum mechanics in2D, which has previously been applied to phase transitionsfrom CSL phases, one expects the universal features ofthe phases and transitions to be equivalent in these two mod-els.The nearest-neighbor model for spin ice, H nn , can bemapped to a system of hard-core lattice bosons, withspin-reversal symmetry replaced by particle–hole symmetry.The Coulomb phase of the spin model is equivalent to asuperfluid, while the saturated ferromagnet with M z = ± M sat is equivalent to the vacuum and fully-occupied statesof the bosonic model, which spontaneously break particle–hole symmetry. The strings of flipped spins that proliferateat 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 , dis-playing infinite-order multicriticality at the transition betweenthese two phases, is given by H = − t (cid:88) (cid:104) i j (cid:105) ( b † i b j + b † j b i ) − V (cid:88) (cid:104) i j (cid:105) (cid:104) ( n i − )( n j − ) − (cid:105) ,(A.1)where b i = | (cid:105) i (cid:104) | i and n i = b † i b i = | (cid:105) i (cid:104) | i are annihilation andnumber operators for hard-core bosons. The first term rep-resents tunneling t between neighboring sites (cid:104) i j (cid:105) , while thesecond is an attractive interaction of strength V > H conserves particle number, the Hilbert spacecan be divided into sectors of fixed density ρ = (cid:104) n i (cid:105) ; let E gs ( ρ )be the ground-state energy in each. For t > V , the overallground state occurs in the sector with ρ = , and the system isa particle–hole-symmetric superfluid. For t < V , E gs is insteadminimized by either the vacuum, ρ =
0, or the fully-occupiedstate, ρ =
1. At t = V , the model has a Rokhsar–Kivelson(RK) point, at which H can be written as a projector, H = t (cid:88) (cid:104) i j (cid:105) ( | (cid:105) i | (cid:105) j − | (cid:105) i | (cid:105) j )( (cid:104) | i (cid:104) | j − (cid:104) | i (cid:104) | j ) , (A.2) ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë - - - V (cid:144) t g r ound - s t a t ee n e r gy E g s H N L ‘ t FIG. 15. Ground-state energy E gs , in the Hilbert-space sectorwith N particles, for hard-core bosons with particle–hole symme-try, Eq. (A.1). Results were obtained using exact diagonalizationon a 4 × E gs ( N ), for the labeled values of N , versus nearest-neighbor attraction V , both in units of the tunneling strength t . Themodel has an RK point at V = t , 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 bosondensity zero and one ( N = N =
16 respectively), which spon-taneously 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 su-perposition of all configurations, has E gs ( ρ ) =
0. As illus-trated in Fig. 15, which shows results of exact diagonaliza-tion (ED) on a small system, this leads to a transition withidentical properties to the ordering transition of spin ice underuniaxial pressure, with F ( M z ) replaced by E gs ( ρ ). (The exactequivalence is established by noting that the transfer matrixfor the classical problem can be written as a projector at thetransition. )Following similar logic to Section II C, one expects thatquartic interactions between bosons should open up an inter-mediate phase with density changing continuously between ρ = and ρ = ffi ces to add a (particle–hole-symmetrized) four-body repulsion H = V (cid:88) i jkl ∈ (cid:3) ( n i − )( n j − )( n k − )( n l − ) , (A.3)where the sum is over sites i jkl around a square plaquette.In this case, there are two continuous transitions, with den-sity changing from (cid:104) n i (cid:105) = to 0 < |(cid:104) n i (cid:105) − | < and then to |(cid:104) n i (cid:105)− | = , as V / t is increased. (The nature of the transitionsis clear even for the small system sizes accessible in ED, be-cause the order parameter commutes with the Hamiltonian.)The transition into the vacuum or fully-occupied (vacuum ofholes) state is described by the standard critical theory for the1 ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë ë - - - V (cid:144) t g r ound - s t a t ee n e r gy E g s H N L ‘ t FIG. 16. Ground-state energy E gs ( N ) versus nearest-neighbor at-traction V , as in Fig. 15, with an additional four-body repulsion V = . t , Eq. (A.3). As in the case with V =
0, the ground state isat half filling for V (cid:28) t and has density either zero or one for V (cid:29) t .In between, there is a regime where the minimum of E gs crosses overbetween the extremes, stepping through each intermediate density inturn. In the thermodynamic limit, this is expected to become a phasewith continuously varying density, separated from small- and large- V / t phases by continuous quantum phase transitions. vacuum transition of bosons, while the transition at lower V / t involves spontaneous breaking of particle–hole symmetrywithin the superfluid, and is described by the critical theoryof Section IV A. In cases where the total particle number isfixed, the latter transition would lead to phase separation intoregions with di ff ering densities.The additional interaction H in the classical spin model,defined in Eq. (3), may be viewed as the equivalent of H .To see this, recall that strings of flipped spins are equivalentto bosons, and that these occur at low density near the tran-sition to saturation (bosonic vacuum). Two strings passingthrough a tetrahedron t change its total spin from + √ ˆ z to − √ ˆ z , and so the interaction H (with V >
0) amounts toan energy penalty when four strings are in close proximity(passing through two tetrahedra on opposite sides of the samehexagon). L. Balents, Nature , 199 (2010). P. W. Anderson, Mater. Res. Bull. , 153 (1973). P. W. Anderson, Science , 1196 (1987) P. A. Lee, Science , 1306 (2008). T.-H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A.Rodriguez-Rivera, C. Broholm, and Y. S. Lee, Nature , 406(2012) R. Moessner and S. L. Sondhi, Phys. Rev. Lett. , 1881 (2001). L. Balents, M. P. A. Fisher, and S. M. Girvin, Phys. Rev. B ,224412 (2002). C. Castelnovo, R. Moessner, and S. L. Sondhi, Ann. Rev. Cond.Matt. Phys. , 35 (2012). C. L. Henley, Ann. Rev. Cond. Matt. Phys. , 179 (2010). S. T. Bramwell, M. J. Harris, B. C. den Hertog, M. J. P. Gingras, J.S. Gardner, D. F. McMorrow, A. R. Wildes, A. L. Cornelius, J. D.M. Champion, R. G. Melko, and T. Fennell, Phys. Rev. Lett. ,047205 (2001). L. Savary and L. Balents, Phys. Rev. Lett. , 037202 (2012);Phys. Rev. B , 205130 (2013). O. Benton, O. Sikora, and N. Shannon, Phys. Rev. B , 075154(2012). M. E. Brooks-Bartlett, S. T. Banks, L. D. C. Jaubert, A.Harman-Clarke, and P. C. W. Holdsworth, Phys. Rev. X , 011007(2014). S. Powell and J. T. Chalker, Phys. Rev. B , 024422 (2008). D. S. Rokhsar and B. G. Kotliar, Phys. Rev. B , 10328 (1991). S. V. Isakov, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. ,217201 (2005). T. Fennell, P. P. Deen, A. R. Wildes, K. Schmalzl, D. Prabhakaran,A. T. Boothroyd, R. J. Aldus, D. F. McMorrow, and S. T. Bramwell,Science , 415 (2009). S. Powell, Phys. Rev. Lett. , 065701 (2012). S. Powell, Phys. Rev. B , 064414 (2013). L. D. C. Jaubert, M. J. Harris, T. Fennell, R. G. Melko, S. T.Bramwell, and P. C. W. Holdsworth, Phys. Rev. X , 011014 (2013). L. D. C. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R.Moessner, Phys. Rev. Lett. , 087201 (2010). L. D. C. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R.Moessner, Phys. Rev. Lett. , 067207 (2008). L. D. C. Jaubert, Ph.D. thesis, Ecole Normale Sup´erieure de Lyon,2009, http://hal.archives-ouvertes.fr/docs/00/46/29/70/PDF/Thesis.pdf . This remark notwithstanding, we note that some other apparentlyfour-spin interactions reduce to a sum of two-spin interactions underprojection to the ice-rule states. For example, an interaction similarto H but giving ± V when the two tetrahedra are polarized in thesame ( + ) or opposite ( − ) directions along ˆ z can be written as S zt S zt (cid:48) and is therefore a sum of two-spin interactions. G. T. Barkema and M. E. J. Newman, Phys. Rev. E , 1155 (1998). A. W. Sandvik and R. Moessner, Phys. Rev. B , 144504 (2006). There are Bragg peaks at wavevector hh (cid:96) (in units corresponding tothe conventional cubic unit cell), for h , (cid:96) integers and h + (cid:96) even,except those where h = h and (cid:96) are even but h + (cid:96) is notdivisible by 4. T. S. Pickles, T. E. Saunders, and J. T. Chalker, EPL , 36002(2008). A. I. Larkin and D. E. Khmel’nitskii, Sov. Phys. JETP , 1123(1969). A. Aharony, Phys. Rev. B , 3363 (1973). E. Brezin and J. Zinn-Justin, Phys. Rev. B , 251 (1976). M. Hasenbusch, K. Pinn, and S. Vinti, Phys. Rev. B , 11471(1999). J. Cardy,
Scaling and Renormalization in Statistical Physics (Cambridge University Press, Cambridge, UK, 1996). P. W. Kasteleyn, J. Math. Phys. , 287 (1963). G. J. Sreejith and S. Powell, Phys. Rev. B , 014404 (2014). S. Powell and J. T. Chalker, Phys. Rev. Lett. , 155702 (2008);Phys. Rev. B , 134413 (2009). D. S. Rokhsar and S. A. Kivelson, Phys. Rev. Lett. , 2376 (1988). S. Sachdev,