Magnetic-Moment Fragmentation and Monopole Crystallization

The Coulomb phase, with its dipolar correlations and pinch-point-scattering patterns, is central to discussions of geometrically frustrated systems, from water ice to binary and mixed-valence alloys, as well as numerous examples of frustrated magnets. The emergent Coulomb phase of lattice-based systems has been associated with divergence-free fields and the absence of long-range order. Here, we go beyond this paradigm, demonstrating that a Coulomb phase can emerge naturally as a persistent fluctuating background in an otherwise ordered system. To explain this behavior, we introduce the concept of the fragmentation of the field of magnetic moments into two parts, one giving rise to a magnetic monopole crystal, the other a magnetic fluid with all the characteristics of an emergent Coulomb phase. Our theory is backed up by numerical simulations, and we discuss its importance with regard to the interpretation of a number of experimental results.


I. INTRODUCTION
The interplay between frustration, topological order, fractionalization, and emergent physics has been the focus of a rapidly increasing body of work in recent years. In themselves, these concepts are not new. Frustration underpins theories of glassiness and has been much discussed since the seminal studies of Anderson and Villain [1,2]. Likewise, concepts of topological order and fractionalization, from the fractional quantum Hall effect [3] and quasiparticles in graphenelike systems [4] to solitons in one dimension [5,6], have been prominent topics in theoretical and experimental condensed-matter physics for over 30 years.
The subtleties of the interdependence between these concepts have been elucidated only recently, in the context of the so-called emergent Coulomb phase of highly frustrated magnetic models. However, the phenomenology of the Coulomb phase, with its characteristic dipolar correlations and emergent gauge structure, is not limited to the realms of frustrated magnetism. Indeed, similar behavior is observed in water ice [7] and models of heavy-fermion behavior in spinels [8], as well as in models of binary and mixed-valence alloys [9].
Henley [10] succinctly sets out three requirements for the emergence of a Coulomb phase on a lattice: (i) Each microscopic variable can be mapped onto a signed flux directed along a bond in a bipartite lattice, (ii) the sum of the incoming fluxes at each lattice vertex is zero, and (iii) the system has no long-range order. Elementary excitations out of the divergence-free manifold are seen to fractionalize, giving rise to effective magnetic monopoles that can interact via a Coulomb potential [11] and (in three dimensions) may be brought to infinite separation with finite energy cost.
In this paper, we demonstrate how the lattice-based magnetic Coulomb phase emerges in a considerably wider class of models than those covered by the conditions stipulated above. More specifically, we show that such a phase exists for a model ''magnetolyte'' irrespective of the magnetic monopole density or of monopole ordering [12]. We introduce the concept of magnetic-moment fragmentation, whereby the magnetic-moment field undergoes a novel form of fractionalization into two parts: a divergence-full part representing magnetic monopoles and a divergencefree part corresponding to the emergent Coulomb phase with independent and ergodic spin fluctuations [13]. Our results apply even for a monopole crystal that is shown to exist in juxtaposition with mobile spin degrees of freedom: a previously unseen coexistence between a spin liquid and long-range order induced by magnetic monopole crystallization. The implications of this field fragmentation are wide reaching and relevant for the interpretation of a number of experiments, as discussed below.

II. THE MODEL
Magnetic monopoles [11] emerge as quasiparticle excitations from the ground-state configurations of the dumbbell model of spin ice. Here, the point dipoles of the dipolar spin-ice model [14] are extended to infinitesimally thin magnetic needles lying along the axes linking the centers of adjoining tetrahedra of the pyrochlore lattice (see Fig. 1) that constitute a diamond lattice of nodes for magnetic charge [15]. The needles touch at the diamondlattice sites so that, by construction, the long-range parts of the dipolar interactions are perfectly screened for the ensemble of ice-rules states in which two needles point in and two out of each tetrahedron [16]. These ground states form a vacuum from which monopoles are excited by reversing the orientation of a needle, breaking the ice rules on a pair of neighboring sites. As vacuums go, this one is rather exceptional, as it is far from empty. Rather, the magnetic moments constitute the curl of a lattice gauge field [17], the Coulomb phase, with manifest experimental consequences for spin-ice materials. These effects include diffuse neutron-scattering patterns showing the sharp pinch-point features [18] characteristic of dipolar correlations, and the generation of large internal magnetic fields despite the status of vacuum [19,20].
In an analogy with electrostatics, the monopole-charge distribution obeys Gauss's law for the magnetic field r ÁH ¼ m , where m is the monopole density. The monopole number is not conserved, and the energetics of the dumbbell model at low temperature correspond to a Coulomb gas in the grand canonical ensemble, in which the phase space of monopole configurations is constrained by the underlying ice rules. One can define a Landau energyŨ ¼ U C À N, where U C is the Coulomb energy of the neutral gas of N monopoles and is the chemical potential such that À2 is the energy cost of introducing a neutral pair of monopoles an infinite distance apart. For temperatures of order À, the Coulomb-gas picture of spin ice must be extended to include doubly charged monopoles (see Appendix A).
The divergence inH is related to the breaking of the ice rules throughM, the magnetic-moment density that itself obeys Gauss's lawr ÁM ¼ À m . This constraint does not, however, define the entire moment density, which can have two contributions through a Helmholtz decomposition. The firstM m is the gradient of a scalar potential and provides the magnetic charge; the secondM d is a dipolar field, can be represented as the curl of a vector potential, and is divergence free [21,22]: In states obeying the ice rules,M m ¼ 0, whileM d ¼M; this decomposition corresponds to an emergent Coulomb phase [10]. Breaking the ice rules and introducing magnetic charge lead to the conversion ofM from the divergence-freeM d to the divergence-full fieldM m . Central to the discussions in this paper is the concept that this conversion is, in general, partial-the divergence-free part is not  4) is þ1 (respectively, À1) on the diamond sites where the blue (respectively, red) monopoles sit. The minority spins (red arrows) are equivalent to dimers positioned on the diamond lattice. (b) Two-in-one-out spin configurations in the confined KII kagome spin-ice phase. (c) Spin and needle configurations for a three-in-one-out vertex carrying an isolated north pole showing magneticmoment fractionalization into divergence-full and divergence-free elements. We emphasize that the divergence-free field emerging here, despite being a Coulomb phase, is different from the standard one in the absence of monopoles described in Ref. [11].
completely suppressed by the presence of the monopoles. Thus, in almost all circumstances, one has the coexistence of two complementary fields, reminiscent of the Hamiltonian splitting in topologically ordered phases [23]. The total charge inside a microscopic volume surrounding a diamond-lattice site i is Q i ¼ À RM Á dS, with dS an outward-pointing surface element. For the dumbbell model, the integral reduces to a discrete sum Q i ¼ À P j M ij , where M ij ¼M Á dS ij ¼ ±m=a, dS ij is an infinitesimal surface element cutting the needle, m is the strength of the magnetic moment, and a is the distance between diamond-lattice sites. The positive (negative) sign is for the needle pointing away from (into) the site i, so that M ij ¼ ÀM ji . An ice-rules configuration with two needles into and two out of the site indeed gives Q i ¼ 0, while three in and one out (three out and one in) gives Q i ¼ þðÀÞ2m=a [11].
The fields for an isolated three-in-one-out vertex can be split into divergence-full and divergence-free parts subject to the constraint that the amplitude of each field element is jM ij j ¼ m=a: The first set of fields satisfies Gauss's law for the charge at the origin; the second set satisfies a discrete divergencefree condition and constitutes a residual dipolar field dressing the monopole [see Fig. 1(c)]. Decomposition further away from the charge could be made by solving the lattice Poisson equation to find the field sets belonging toM m that would be subtracted from the M ij to find the discrete elements ofM d . Singly charged monopoles leave a residual contribution toM d at each vertex, and it is only when the vertex is occupied by a doubly charged monopole for which ½M ij ¼ ±ðm=aÞð1; 1; 1; 1Þ that the contribution tõ M d is totally suppressed. Hence, a fluid of singly charged monopoles should be accompanied by a correlated random dipolar field whose detailed structure is updated by the monopole dynamics and only destroyed on the temperature scale at which doubly charged monopoles proliferate. Indeed, the pinch points in diffuse neutron scattering from spin-ice materials are maintained up to surprisingly high temperatures [18,24], indicating the presence of such a dipolar field. The emergence of the dipolar field is further illustrated in Appendix C, where we show how the ½M ij are divided around a pair of isolated nearest-neighbor charges of the same sign. The random fluctuations in the underlying gauge field can be ironed out by breaking the two-sublattice translational symmetry of the diamond lattice, creating a monopole crystal with north and south poles localized on different sublattices.
For an ideally ordered array, the divergence-free fields on alternate sites are perfectly satisfied by the sets ½M ij d ¼ þðÀÞðm=aÞð1=2; 1=2; 1=2; À3=2Þ. Thus, one sees the emergence of a new Coulomb phase with extensive entropy superimposed on monopole order, in which each vertex has three contributions to the dipolar field of strength 1=2 [in units of ðm=aÞ] and one of strength 3=2, which is shared between a pair of neighboring sites on opposite sublattices. This fragmented state could be termed an ''antiferromagnetic Coulomb magnet'' by analogy with the ''ferromagnetic Coulomb magnet,'' predicted in the gauge mean-field theory of quantum magnets on a pyrochlore lattice [25]. The ordered component corresponds to a broken symmetry of the Ising spins described in the local-axis reference frame into the three-up-one-down or three-down-one-up sector. For the divergence-free partM d , placing a dimer along the bond of strength 3=2 provides a mapping between the emergent dipolar field and hard-core dimers on the (bipartite) diamond lattice [26]. The extensive entropy of the dipolar field is thus associated with closed loops of dimer moves [27]. Introducing quantum-loop dynamics gives rise to a U(1) liquid phase close to the Rokhsar-Kivelson point [28,29].
A monopole-crystal ground state can be induced in the dumbbell model by modifying the chemical potential so that the total Coulomb energy U C outweighs the energy cost for creating the particles ÀN. For the monopolar crystal is the Coulomb energy for a nearest-neighbor pair of monopoles of charge ±Q, 0 is the permeability, N 0 is the number of diamond-lattice sites, and ¼ 1:638 is the Madelung constant. We define a reduced chemical potential Ã ¼ =u, and thus the ground state should be a monopole crystal for Monopole crystallization has also recently been studied within the dipolar spin-ice model, in the canonical ensemble, that is, with fixed monopole number [30], leading to a region of phase separation between the crystalline and the fluid phases. For classical spin ice, crossing this phase boundary would correspond to leaving the spin-ice phase [14], at which point double charges become favored. The ordering is then to a structure in complete analogy with zinc blende: a crystal of doubly charged monopoles for whichM d is everywhere zero, corresponding to the ''all-in-all-out'' (AIAO) magnetic order observed in FeF 3 [31].
In the modeling of spin ice, the chemical potential can be extracted from the dumbbell approximation to the dipolar spin-ice model (see the Supplementary Information of Ref. [11]). To put spin-ice materials in the context of the present work, we note that the dumbbell approximation for dysprosium titanate yields ¼ À4:35 K, while direct simulations for dipolar spin ice give ¼ À4: 46 K [32], so that Ã % 1:42, well away from the monopole-crystal phase boundary. In Appendix A, we return to this modeling and show that the magnetostatics of monopoles leads to a prediction for the phase boundary where the ratio of nearest-neighbor exchange to dipolar energy scale reaches J nn =D nn ¼ À0:918, in excellent agreement with direct analysis of the dipolar model [33]. Further, we show that, for a pair of doubly charged monopoles, both the chemical potential and the Coulomb energy are scaled by a factor of 4 compared with the values for the singly charged monopoles discussed in the main text. Hence, the zerotemperature spin-ice-monopole-crystal phase boundary occurs at the same point, whether one includes or excludes the doubly charged monopoles that complete the magnetic charge description of spin ice. In this paper, we explicitly exclude double charges, taking us away from traditional spin-ice modeling for large monopole concentrations, a point we return to in Sec. IV when discussing the experimental relevance of our results.

III. RESULTS
We have tested these ideas directly through Monte Carlo simulations of the dumbbell model (details are given in Appendix B). In Fig. 2, we show the evolution of an order parameter for monopole crystallization M c and monopole number density n, as a function of reduced temperature where q i ¼ ðQ i a=2mÞ ¼ fÀ1; 0; þ1g is the topological charge on site i, Á i ¼ ±1 is a diamond-sublattice index (see Fig. 1), h. . .i denotes a statistical average, and n ¼ hNi=N 0 . The data show clear evidence of monopole crystallization at a transition temperature T Ã C that varies with Ã . At this temperature, a lattice fluid gives way to a phase with reduced translational symmetry, in which M c approaches unity. Debye-Hückel theory for an unconstrained Coulomb gas on a bipartite lattice predicts a line of second-order transitions in the ð Ã ; TÞ plane, becoming first order via a tricritical point as Ã increases [34]. Our data are consistent with this scenario, despite the additional constraints of the dumbbell model. From the finite-size scaling analysis shown in Appendix E, we estimate a tricritical point Ã tr % 0:78, T Ã tr % 0:13. This temperature is comparable to that obtained from the numerical simulation of a cubic-lattice Coulomb gas [35], although Ã tr is surprisingly close to Ã 0 . For small values of Ã , a continuous transition takes the system from a high-density fluid (n > 4=7) to the crystalline phase. However, as Ã increases toward the phase boundary at approximately 0.8, the fluid density is able to reach lower values near T c [shown by the data for Ã ¼ 0:794 in Fig. 2(b)], indicating that while the crystalline ground state is energetically favored, the finite monopole density of the fluid phase is stabilized by entropy. It is this minimum in the density as the transition is approached that drives the transition to first order, as in the Blume-Capel model for spin-one systems [35]. More work is required to extract the effect of the constraints in detail and to establish the tricritical parameters with precision. There could, in principle, also be a liquid-gas transition at higher temperature, between lowand high-density fluids, but at the level of Debye-Hückel theory, this transition is suppressed by the monopole ordering [34]. There is no strong evidence of this liquidgas transition in our simulations, although the crossover from a high-density fluid (region II of the phase diagram in The transition appears to pass from second order to first order via a tricritical point for Ã tr % 0:78. The very narrow first-order region ends at Ã % 0:80, above which spin-ice physics is recovered (e.g., the open black circles). The alternating positive or negative charges on the diamond lattice can also be seen as stacked monolayers of monopoles of the same charge in all three cubic directions. Note that the configurational constraints of the singly charged monopole fluid result in a high-temperature limit for the density of n ¼ 4=7, rather than the n ¼ 0:5 expected for an unconstrained bipartite lattice gas. A generalization of the worm algorithm has been developed to ensure equilibration at low temperature (Appendix B). monopole-crystal phase boundary is quite sharp and could be considered as a vestige of such a transition.
In Fig. 3, we show the resulting phase diagram, as mapped out by the divergence of the specific heat C at the phase transition. The monopole-crystal phase terminates for Ã % 0:8, in approximate agreement with our prediction of Ã 0 ¼ 0:819, the small difference being most likely due to finite-size effects exacerbated by the longrange interactions. In Fig. 4, we show a simulated elastic neutron-scattering map determined within the static approximation, setting the magnetic form factor equal to unity and by averaging over 2000 randomly selected configurations of the ideal monopole crystal. For scattering purposes, the magnetic needles are once again taken as unit vectors (spins)S i on the sites of the pyrochlore lattice. The structure factor SðQÞ is dominated by Bragg peaks at the (220) positions, characteristic of the all-in-all-out structure observed in FeF 3 [31]. The intensity of these peaks is precisely one quarter of that expected for complete all-in-all-out ordering and is consistent with scattering from monopoles constructed from fragmented spins with effective length 1=2. SðQÞ also reveals diffuse scattering with the clearly defined pinch points of a Coulomb phase. We have compared the intensity of the diffuse scattering with that found when the length of the minority spin at each vertex is extended to three while the length of the majority spins remains fixed at unity. The resulting structure factor has no Bragg peaks and has 4 times the intensity of the true diffuse scattering for a given wave vector. The ensemble of Bragg peaks plus emergent Coulomb phase therefore appears in excellent agreement with the predicted moment fragmentation of Eq. (2).
Magnetic charge crystallization also occurs for the dumbbell model on the kagome lattice (see Fig. 1) [36][37][38][39][40]. Breaking a Z 2 symmetry for needle orientations on each triangular vertex plunges the system from the unconstrained ''KI'' phase into the constrained ''KII'' Coulomb phase, realized in spin-ice materials by applying a field along the (111) direction of the pyrochlore lattice [41,42]. The mo- FIG. 3. The position of the singularity in the specific heat C in the ð Ã ; T Ã Þ plane traces out the phase boundary for the monopole-crystal phase. The most intense peak signals the passage of the transition from second to first order via a tricritical point. Indexes I, II, and III mark the monopole crystal, highdensity fluid, and low-density fluids, respectively. Regions II and III are separated by a continuous crossover, and the gas of monopoles (III) corresponds to the low-temperature phase for spin ice.
FIG. 4. Simulated unpolarized neutron-scattering structure factors SðQÞ for the pyrochlore monopole crystal (left) and for in-plane scattering from kagome ice (right). The pyrochlore SðQÞ has been calculated by averaging over 2000 distinct monopole-crystal ground states of a lattice with L ¼ 8. In order to reveal the diffuse scattering, the Bragg peaks in the pyrochlore data are plotted as contours in gray scale superimposed on the contribution to SðQÞ from the dipolar field. The wave vector Q is in units of 2=a 0 , where a 0 ¼ 4a= ffiffiffi 3 p is the lattice parameter of the cubic unit cell of the pyrochlore lattice. The kagome-ice data are taken from Ref. [102].
ments of a vertex with configurations constrained to two needles in and one out can, as above, be decomposed into divergence-full and divergence-free parts: Consequently, simulated neutron-scattering plots for the ensemble of constrained states have both Bragg peaks and pinch points characteristic of a two-dimensional Coulomb phase (Fig. 4). The Bragg peaks have an intensity of 1=9 of those for a fully ordered all-in-all-out phase. While deconfined monopole excitations away from these states carry a magnetic charge Q ¼ 2m=a, as in spin ice, the charge ordering corresponds to a crystal of objects with charge Q=2 ¼ m=a, providing a simple example of frustration-driven charge fractionalization [36][37][38]43].
Returning to the three-dimensional system, the persistent background fluctuations are further evidenced by studying local dynamics. We have collected two sets of data, one using local single spin-flip Metropolis dynamics, the other using a nonlocal worm algorithm [44] extended to include long-range interactions [45]. While the nonlocal algorithm is extremely powerful for extracting equilibrium properties in the highly constrained monopole-crystal phase, the local dynamics is of great interest, as it provides insight into how real systems might evolve with time [32,46]. The one-site, two-time monopole-and spinautocorrelation functions are defined as In region I of our phase diagram, the charge-autocorrelation function C c ðtÞ remains close to unity over the whole time window, reflecting the broken charge symmetry and the localization of monopoles, as shown in Fig. 5. The spinautocorrelation function, however, shows a decay over a modest simulation time, from unity to an asymptote of C s ðt ¼ 1Þ ¼ 1=4, reflecting the random projection of the spin onto the four orientations of a three-in-one-out or three-out-one-in vertex. This spin ergodicity, superimposed on a background of magnetic Bragg peaks, illustrates the collective nature of the monopole excitations: singular nodes in a fluctuating magnetic fluid, rather than static microscopic objects. We have evidence for the validity of these conclusions deep into region I, with simulations down to approximately T C =2 for each Ã . C s ðtÞ would be accessible experimentally, as it provides the diagonal contribution to the ac magnetic susceptibility ð!Þ, so that this finite time scale should show up as a characteristic frequency. For the perfect monopole crystal, the dynamics are restricted to dimer-loop moves [27,47] if the system is to remain on the ground-state manifold. Single spin-flip dynamics lead to excitations away from these ground states, initially by monopole-pair annihilation. In region I, single spin-flip dynamics are dominated by needle flips that destroy and recreate nearest-neighbor pairs. Short-lived excitations of this kind collectively displace the fictive dimer positions approximating to the loop dynamics of dimers (Appendix D). The energy scale for an isolated excitation of this kind dŨ $ Àu þ 2 goes to zero at the ground-state phase boundary, allowing for extensive local dynamics in this region. Hence, for M c close to unity, the system is ergodic while at the same time retaining Coulomb-phase correlations.
We now turn to an important consequence of the field fragmentation for the interpretation of experimental data. The fluctuating background of the Coulomb phase appears The time required for the autocorrelation function to reach its longtime asymptote is significantly greater at Ã ¼ 0:41 than at Ã ¼ 0:57. This behavior is indicative of a slowing of the dynamics as the monopole density increases. As Ã increases, n decreases [see Fig. 2(b)] and the t ¼ 1 asymptote of C s ðtÞ is slightly reduced relative to the value for the ideal monopole crystal. We note that both panels correspond to points far into region I. to obscure the phase transition from bulk magnetic measurements. In Fig. 6(a), we show the magnetic susceptibility as a function of temperature for different values of Ã . At this level of analysis, the susceptibility is virtually featureless through the transition, showing no evidence of the characteristic cusp that one might expect at an antiferromagnetic transition in an Ising system. As one moves toward the tricritical point, a very weak feature does appear, driven by the huge monopole-density change as the system passes through the transition; however, the unusual characteristic (for a magnetic phase transition) of being virtually transparent to the bulk susceptibility remains essentially intact. A more detailed analysis of the susceptibility does, however, yield interesting information. It was recently shown that spin-ice models show a crossover in the Curie constant C as the system moves from the uncorrelated high-temperature phase to the low-temperature Coulomb phase [15,48,49]. It was demonstrated that taking the needles as scatterers of unit length gives rise to a crossover from C ¼ 3T ¼ 1 to C % 2. A similar Curie-law crossover is observed in the current work, as shown in Fig. 6(b), where C is seen to evolve from 4=3 at high temperature, as expected for a paramagnetic 14-vertex model, to a value C $ 3=2 on entering the monopole-crystal phase and again to C % 2 on entering the constrained monopole vacuum. The change from a second-to a firstorder transition can clearly be seen from this evolution. In the first-order region, C evolves above 3=2 as the monopole density drops in the fluid phase before falling discontinuously at the transition onto the 3=2 plateau.

IV. RELATION TO EXPERIMENT
One of our goals has been to construct a minimal model, based on spin ice, in which monopole order can be shown to coexist with (Coulomb-phase) spin-liquid physics. Spin ice is a central pillar in the ever-expanding field of frustrated magnetism. More generally, models of, and experiments on, systems based on pyrochlore and kagome lattices have resulted in a wealth of often puzzling results stemming from inherent geometrical frustration. We propose our model as a step toward answering some of these questions. The rest of this paper is dedicated to the consideration of experimental systems in relation to our model and its effects, including rare-earth oxides, spin-ice candidates, artificial spin ice, and the use of magnetic fields to induce a staggered chemical potential for magnetic charge.

A. Magnetic field
An external magnetic field couples to both the monopoles and toM d , providing a gradient to the chemical potential, inducing transient monopole currents and ordering the dipolar field (the two processes being intimately related [11,15,46,50]). Applying a field in the [111] direction imposes kinetic constraints to monopole movement, restricting them to planes perpendicular to the field direction. In this scenario, the field constitutes a staggered chemical potential, breaking the Z 2 symmetry between the two sublattices, allowing experimental access to monopole crystallization [11]. This evolution corresponds to the first-order phase transition observed as spin-ice materials leave the plateau region of the ðB; TÞ phase diagram. The corresponding phase boundary terminates at a critical end point ðB C ; T C Þ [51]. If one were to start from the chargeordered Coulomb phase in zero field, an external (111) field would couple uniquely toM d . The system would then order via a three-dimensional Kasteleyn transition in complete analogy with that in two dimensions, driven by tilting the field off the (111) axis [52]. An external magnetic field can also stabilize a single-charge monopole in Tb 2 Ti 2 O 7 [53].

B. Spin-ice candidates
Experimental relevance in zero field is a more open question. The chemical potential can be reduced within the confines of classical spin ice by reducing the lattice parameter, as is the case for the material Dy 2 Ge 2 O 7 [54]. Here, the monopole number certainly increases but their proliferation is accompanied by the generation of double monopoles with charge ±2Q. Close to the spin-ice, all-in-all-out phase boundary, the energy and entropy balance of the monopole crystal could possibly stabilize it ahead of either the monopole fluid or the fully ordered doublemonopole crystal and it would certainly be interesting to study this problem both numerically and experimentally through high-pressure experiments or further substitution of smaller ions. Replacing germanium with silicon is, for example, a challenging possibility [55]. A further route to stabilization could be quantum fluctuations [25,56,57]. For systems close to the spin-ice-antiferromagnetic phase boundary [14], one might hope that zero-point fluctuations of the fragmented dipolar field could stabilize the monopole crystal over the classical all-in-all-out spin structure or the spin-ice manifold [28,29]. This situation would require both transverse spin fluctuations and dipole interactions between the moments, allowing the chemical potential to vary while permitting perturbative quantum spin fluctuations about the local h111i axes of spin ice. Another possibility is the generation of a staggered monopole chemical potential through a distortion of the lattice structure and the breaking of the crystal electric-field symmetry. Lifting the doublet degeneracy corresponding to the Ising-like spin-ice degrees of freedom in an ordered manner could thus lead to a perturbation that couples to monopoles but not to the dipolar field, favoring monopole crystallization.

C. Tb 2 Ti 2 O 7
Fifteen years of intense research have made Tb 2 Ti 2 O 7 one of the most intriguing rare-earth frustrated magnets, sitting somewhere between a spin liquid [58] and quantum spin ice [59]-or maybe spanning both. It is this dual nature that makes it an interesting case study for magnetic-moment fragmentation. Tb 2 Ti 2 O 7 has a negative Curie-Weiss temperature Â CW ¼ À14 K [60]. As such, it can be considered as an antiferromagnet and numerical simulations of the corresponding dipolar spin-ice model give a phase transition to the all-in-all-out state at 1.2 K [14]. At ambient pressure, it fails to develop magnetic order down to at least 50 mK [61]. Recently, diffuse neutron scattering from single-crystal samples has exposed pinch-point-scattering patterns, indicating the presence of Coulomb-phase correlations [62,63], albeit somewhat deformed compared with those observed for classical spin-ice materials [18].
This picture is, however, very sensitive to perturbations. Under high pressure, the liquidlike phase transforms into a partially ordered structure where each vertex has one spin along a spin-ice axis and three collinear spins [64]. Vertices on the two-diamond sublattices have spin configurations with mirror symmetry, and the ordered state coexists with a fluctuating magnetic background. There is also evidence that a small magnetic field imposed along the [111] direction induces weak all-in-all-out order, leading to the material being described as an ''incipient AIAO antiferromagnet'' [65]. A strong field in the [110] direction stabilizes a double-layered structure of singly charged monopoles [53,66]. Recent experiments on polycrystalline samples have also shown that a very small amount of Tb 3þ stuffing produces an as yet unexplained long-range order, accompanied by weak antiferromagnetism [67]. Finally, both strong magnetoelastic coupling linked to quantum spin-liquid behavior [68,69] and splitting of the singleion ground-state doublet [63,[70][71][72] have been reported.
It would be naive to suggest a quantitative connection between the rich behavior of Tb 2 Ti 2 O 7 and the simple classical model presented here, but at a qualitative level, the similarities are striking. In the monopole-crystal phase, the model is an antiferromagnet with Coulomb-phase correlations. Placing it in the fluid phase but close to the monopole-crystal phase boundary, one could have residual pinch-point correlations from the underlying gauge field emanating from the correlated monopole fluid. Higher pressure could then send the system over the line into the crystalline state through a modification of the chemical potential [73,74]. The experimentally observed spin configuration in the high-pressure phase is quite different from the three-in-one-out spin-ice vertices of the monopole crystal, but it does share the same two-sublattice structure, to which one would have to add transverse spin relaxation off the spin-ice axes. The background of magnetic fluctuations is consistent with the fluctuating dipolar degrees of freedomM d . The incipient, or partial, all-in-all-out behavior in the presence of a field along [111] is exactly what one would expect of our model when sat close to the monopolecrystal phase boundary. As for the influence of the [110] field, since it does not give rise to monopoles at very low temperature in single crystals of Ho 2 Ti 2 O 7 or Dy 2 Ti 2 O 7 [75,76], the double-layer structure of monopoles observed in Tb 2 Ti 2 O 7 is a strong indication that such singly charged monopoles can be stabilized via internal couplings, even if the microscopic mechanism remains unknown. Spin-lattice coupling could be particularly relevant in this context [53,77].
Turning to dynamics, the freezing observed in Tb 2 Ti 2 O 7 [61,78,79] only involves a fraction of the spins (* 10%) and has been shown to be different from spin-glass physics [79]. Such partial spin freezing seems consistent with magnetic-moment fragmentation where only a fraction of the degrees of freedom order; the precise value of this fraction could then be mediated by quantum fluctuations. Hence, while microscopic modeling of Tb 2 Ti 2 O 7 is beyond the scope of this paper, we do propose magneticmoment fragmentation as a promising route to understanding the apparent coexistence in this material of antiferromagnetism with the fluctuating Coulomb-phase physics of a frustrated ferromagnet [62].

D. Other materials
A second quantum spin-ice candidate is Yb 2 Ti 2 O 7 [56,80]. This material, with a Curie-Weiss temperature estimated at around 600 mK [81], shows an unusual phase transition at 200 mK [82], with apparently no accompanying magnetic order [82,83] and magnetic dimensional reduction [84] in the high-temperature phase. There are, however, reports of a (partial) ferromagnetic ordering [84,85] at approximately 250 mK. The magnetic anisotropy in Yb 2 Ti 2 O 7 was initially considered to be XY-like [82], but more recent analysis has suggested that it could in fact be considered as a spin ice with the low-temperature behavior experimentally close to a quantum spin liquid [56,80]. Within this context, a quantum spin-liquidclassical spin-gas transition has recently been proposed [86]. With such complex behavior and sample dependence [87,88], a quantitative understanding requires a detailed microscopic approach [56,80,89]. That being said, the joint features of a low-temperature magnetically fluctuating phase and the presence of a phase transition that is partially transparent to magnetic probes are not unlike the fragmentation-driven transition in this paper, and the concepts developed here could be of use in understanding this complex material.
The understanding of Tb 2 Sn 2 O 7 also remains incomplete. This material orders in a ferromagnetic structure with spins canted off the local spin-ice axes [90,91], as predicted by combining dipolar interactions and spin relaxation [92]. However, the ordering is accompanied by an as yet unexplained fluctuating magnetic background, while the correlations above the transition appear to be antiferromagnetic. Although the details will almost certainly be different, magnetic-moment fragmentation does seem to be at play here and the concept could be of use in understanding this ordered, yet fluctuating, system.

E. Artificial spin ice
There are immediate experimental consequences for our results for charge ordering in two dimensions. We have shown here that the KII phase on a kagome lattice, which was previously believed to be magnetically disordered, actually has partial all-in-all-out order [Eq. (5)]. Crystallites of the KII phase have recently been realized in permalloy nanoarrays with a honeycomb structure [93]. A simulated neutron-scattering analysis of the dipole orientations of the sample should therefore yield Bragg peaks of reduced intensity, similar to those observed in Fig. 4. Artificial spin-ice systems could therefore provide direct experimental realizations of magnetic-moment fragmentation. The ð2; 2; 0Þ Bragg peaks characteristic of two-dimensional charge ordering should, in principle, also occur in spin-ice materials with a field along the ð1; 1; 1Þ direction, although they may be masked by the fieldinduced magnetic order.

V. CONCLUSIONS
In conclusion, we have shown how, through the presence of singly charged monopoles, a gauge field emerges from the dumbbell model of spin ice [11] that only partially maps onto the physical degrees of freedom, the magnetic needles. As a consequence, the intrinsic moments fragment into two parts. The first part satisfies the discrete Poisson equation on a diamond lattice, giving the magnetic monopoles, but does not exhaust the magnetic resources associated with each vertex. What remains forms an emergent dipolar field that evolves through monopole dynamics. By varying the chemical potential for monopole-pair creation, one can observe a monopole-crystallization transition, below which the gauge field provides a fluctuating and ergodic magnetic background with Coulomb-phase correlations. An analogous description exists for magnetic charge ordering in the KII phase of magnetic needles on a kagome lattice.
Order, or partial freezing in the presence of a fluctuating magnetic background, is a recurring phenomenon in frustrated magnetism (see, for example, Refs. [65,82,90,94,95]). Here, the magnetic-moment fragmentation leads naturally to persistent spin fluctuations within a purely classical model based on spin-ice physics. These background fluctuations mask the magnetic phase transition from view in susceptibility measurements, a phenomenon that can also occur in experiments on rareearth pyrochlores [84,96,97]. It will be interesting to see if this concept of partial emergence can provide a more generic mechanism for persistent spin fluctuations in other situations. Finally, recent studies of quantum spin ices [25,98,99] have revealed a complete model for QED with magnetic monopoles, conjugate electric poles, and photon excitations. Including magnetic-moment decomposition within this model for QED opens the possibility for new levels of fractionalization, such as fractional charge and spin-charge separation.

ACKNOWLEDGMENTS
We thank S. T. Bramwell for discussions concerning the simulated neutron-scattering plots and issues related to Ref. [22]. It is also a pleasure to thank O. Benton, M.

APPENDIX A: THE COULOMB PHYSICS OF SPIN ICE
The dipolar spin-ice Hamiltonian can be written (see the Supplementary Information of Ref. [11]) within the dumbbell approximation as where Q i is the total magnetic charge on diamond-lattice site i and v 0 is an on-site term whose value is calculated from estimating spin-flip energies in the dipolar model. H 0 is the ground-state energy for a Pauling state within the dumbbell approximation H 0 ¼ ÀðN 0 =2Þv 0 Q 2 , with Q ¼ 2m=a the monopole charge. The ice rules and their consequent violation impose that Q i ¼ 0; ±Q; ±2Q only, and the diagonal term provides the chemical potentials for both singly () and doubly charged ( 2 ) monopoles: where ¼ Àv 0 Q 2 =2, 2 ¼ À2v 0 Q 2 and where the numbers of single and double monopoles are N and N 2 , respectively. The sketch below illustrates the energy scale, taking a single vertex or a diamond-lattice site from a twoin-two out ground state to a three-in-one-out monopole and finally to an all-in or all-out double monopole.
Neglecting the double monopoles, the internal energỹ U ¼ U C À N and we have the lattice Coulomb gas studied in the main text.
This analysis shows that both the chemical potential and the Coulomb energy scale by a factor of 4 when moving from single to double monopoles. Hence, the zerotemperature phase transition from monopole vacuum to monopole crystal [Eq. (3)] holds for both single-and double-monopole crystals, giving in both cases Below this threshold, if both species are present, the excess Coulomb energy of the double monopoles wins, ensuring the predicted all-in-all-out ground state outside the classical spin-ice phase. One can use the monopole crystallization as a criterion to estimate the position of the spin-ice phase boundary: Following Ref. [11], where J is the nearest-neighbor (antiferromagnetic) exchange constant and D ¼ 0 m 2 4k B a 3 is the strength of the dipole interaction between the moments of dipolar spin ice. The Coulomb interaction between nearest-neighbor monopoles can be written as juðaÞj ¼ 8 and hence where J nn ¼ J=3 and D nn ¼ 5D=3. This Coulomb-gas estimate is in excellent agreement with numerical estimates for dipolar spin ice. Melko et al. [33] find J nn =D nn ¼ À0:905 with hysteresis down to J nn =D nn ' À1, the origin of the difference being the small bandwidth for the Pauling states that is neglected in moving to the magnetic charge description.
In the present paper, the double charges are suppressed, leading to the monopole crystal with finite zero-point entropy and magnetic-moment fragmentation. In model spin ice, any perturbation that displaces the equality 2 ¼ 4 in favor of single monopoles will generate the monopole-crystal phase in a band between the spin ice and the double-monopole-crystal phases. It is possible that the phase could be stabilized near this phase boundary by quantum fluctuations [100], thermal fluctuations, or a staggered, sublattice-dependent chemical potential. This point is addressed further in Sec. IV B, where we discuss the relevance of our work to experiment.

APPENDIX B: SIMULATIONS
The numerical results in this paper are obtained by simulating the dumbbell model [11], equivalent to a gas of singly charged magnetic monopoles. The energy scale of the Coulomb interactions is entirely determined by uðaÞ, the energy scale for nearest-neighbor monopoles. The chemical potential Ã is then a free parameter.
Three variants of Monte Carlo simulations are employed: (i) a single spin-flip Metropolis update (SSF), to reproduce the local dynamics relevant for classical spin-ice materials such as Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 (used for the results in Figs. 3 and 6); (ii) joint dynamics of SSF and a worm algorithm specifically designed for the current model (used for the results in Fig. 5; see below for more details on the worm algorithm); and (iii) SSF with worm updates and parallel tempering [101] (used for the results in Fig. 2).
The results presented in Fig. 3 are obtained using a system comprised of 8L 3 ¼ 1000 diamond-lattice sites where L ¼ 5 is the number of cubic unit cells in each spatial dimension. The system is equilibrated over t eq ¼ 10 4 Monte Carlo steps per diamond-lattice site (MCS/s) with data collected over a further 10 5 MCS=s.
For the results in Fig. 6, we use a system with L ¼ 7 (N 0 ¼ 2744) with t eq ¼ 10 4 MCS=s (for both values of Ã ) after which observations of the autocorrelation function are made every MCS/s for a total of 2:5 Â 10 4 MCS=s ( Ã ¼ 0:57) and every 100 MCS=s for a total of 10 6 MCS=s ( Ã ¼ 0:41). Each point on these plots is the result of averaging over 100 consecutive observations. The density of monopoles n in Eq. (5) is chosen as that at t ¼ 0. Figure 5 shows data obtained for a system with L ¼ 4 (N 0 ¼ 512). After annealing from high temperature to the temperature T of interest over a period of 10 5 MCS=s, the system is equilibrated at temperature T for a further t eq ¼ 10 5 MCS=s prior to the data-collection period lasting 10 6 MCS=s, during which observations are made every 10 MCS/s. Fifty worm updates are performed every 10 MCS=s to facilitate thermalization. We run six independent simulations for each value of the parameter Ã ; the error bars are the standard deviations of these six samples at each temperature. Similarly, for the data in Fig. 2, the parameters are L ¼ 8, t eq ¼ 10 4 MCS=s, with an observation period of 10 5 MCS=s and averaging over four independent simulations. Again, 50 worm updates are performed every 10 MCS=s; 100 different temperatures between 0.2 and 0.6 K are used for parallel tempering.

Worm algorithm
In the absence of interactions between particles, the free energy of a system in the grand canonical ensemble only depends on the density of charges. With the addition of Coulomb interactions, the free energy also depends on the position of the charges. Hence, an update that (i) does not modify the number of charges or their positions and (ii) respects detailed balance, i.e., has the same probability flux to be formed and erased, will necessarily be rejection free. In the absence of double charges, if we randomly choose an initial tetrahedron and spin (say, pointing ''in''), it will always be possible to move forward and start a worm by flipping an ''out'' spin on the chosen tetrahedron. The number of out spins can be 1 (three in and one out), 2 (two in and two out), or 3 (three out and one in). Given that these choices remain the same irrespective of whether the worm is being created or destroyed, detailed balance is obeyed. When the worm closes on itself, it can be flipped at no energy cost while respecting detailed balance; the update can be accepted with probability 1. The strength of this algorithm is that it is rejection free in the Coulomb phase (two in and two out), in the dimer covering of the diamond lattice (alternating three in and one out and three out and one in), and for all densities of monopoles in between. This worm algorithm could also be developed for the dipolar spin-ice model [14]. In this case, an additional global Metropolis argument would be needed to take into account the degeneracy lifting between states due to corrections of quadrupolar order when the needles of the dumbbell model are replaced by point dipoles on the nodes of the pyrochlore lattice.

APPENDIX C: FIELD DISTRIBUTIONS
As an example of magnetic-moment fragmentation in the monopole fluid phase, we show, in Fig. 7, two isolated neighboring north poles on a square lattice. It is useful to consider this case (even though we do not consider in detail the dumbbell model on a square lattice [36] in this paper), as the fields can be easily visualized. Starting at nine o'clock and turning clockwise, the fields for sites 1 (on the left) and 2 (on the right) can be decomposed as follows: (C2) where in each case the first and second terms are the contributions to the divergence-fullM m and divergence-freeM d fields, respectively.

APPENDIX D: THE DYNAMICS OF DIMER FLIPS
Single spin-flip Metropolis dynamics below the crystallization transition are dominated by needle flips that create and destroy north-south monopole charges. In Fig. 8, we show a typical sequence of moves for the monopole-crystal phase on a square lattice. In the first instance, the vertical needles flip independently, thus destroying the two neutral pairs of monopoles and reestablishing the ice rules on the four vertices of the square plaquette. These moves are followed by flipping of the horizontal needles that reestablishes locally the monopole crystal. A net consequence of such a sequence is to flip the fictive dimers from a horizontal to a vertical arrangement, as in a dimer-loop move [27]. Similar sequences occur for the pyrochlore and kagome lattices considered in this article, for which the shortest loop is a hexagon comprised of six needles or, equivalently, three dimers.

APPENDIX E: FINITE-SIZE SCALING
The susceptibility of the monopole-crystallization order parameter is defined as In order to characterize more quantitatively the nature of the phase transitions, we perform finite-size scaling on the maxima of the specific heat C and susceptibility c , plotted in Fig. 9. As explained in the paper, two regimes clearly appear. The transition is continuous up to a tricritical point at Ã tr ¼ 0:78±0:01, where it becomes first order before disappearing for Ã > 0:800±0:05. Our simulations suggest the continuous transition line to be of the 3D Ising universality class-consistent with a Debye screening of the long-range interactions-however, distinguishing between Ising and mean-field exponents is a difficult task [34] that would require further numerical and/or theoretical effort.
The discontinuity of the order parameter in Fig. 2 of the main text strongly supports the first-order nature of the phase transition, but its quantitative signature in finite-size scaling is rather challenging. That is, the diverging correlation length on either side of the tricritical point could lead to overestimates of Ã tr and T Ã tr . The first-order regime in our system occupies only a small region of parameter space, making it difficult to separate first-order from tricritical behavior. Increasing the system size beyond the correlation length rapidly becomes very time consumingthe long-range nature of the Coulomb interactions makes the CPU time scale as L 6 . This computational cost is compounded by the relative inefficiency of the parallel tempering algorithm for first-order transitions, especially in large systems. Nonetheless, it has been possible to show a sharp increase of the scaling exponents close to the lowtemperature phase boundary, their values approaching those of a first-order transition (see the lower panels of Fig. 9). In particular, the maximum in C develops scaling behavior in this region-close to = ¼ 3 for L ! 4 and Ã ¼ 0:796. An interesting consequence of our theory is the appearance of critical correlations even in the spin-ice regime where there is no phase transition. In Fig. 9, the green data points ( Ã ¼ 0:801) are constant for L > 4, as expected for a spin-ice crossover into the two-in-two-out Coulomb phase. However, for very small systems, both C and c seem to scale the same way as in the first-order region. This behavior suggests that spin-ice materials and models Divergence-full and divergence-free field distributions for two isolated nearest-neighbor north poles (particle 1 on the left and 2 on the right) for the dumbbell model on a square lattice. Each chevron corresponds to a magnetic-moment field strength of m=2a and a dotted line to zero-field strength. The blue circles represent north poles of charge 2m=a. close enough to the low-temperature phase boundary can exhibit correlations inherited from the monopole crystallization, which should be visible using a local probe such as neutron scattering.