M E Brooks-Bartlett

S T Banks

L D C Jaubert

A Harman-Clarke

P C W Holdsworth
email: peter.holdsworth@ens-lyon.fr

Magnetic-Moment Fragmentation and Monopole Crystallization

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 [START_REF] Anderson | Ordering and Antiferromagnetism in Ferrites[END_REF][START_REF] Villain | Spin glass with non-random interactions[END_REF]. Likewise, concepts of topological order and fractionalization, from the fractional quantum Hall effect [START_REF] Stormer | The fractional quantum Hall effect[END_REF] and quasi-particles in graphene-like systems [START_REF] Hou | Electron Fractionalization in Two-Dimensional Graphene-like Structures[END_REF] to solitons in one-dimension [START_REF] Jackiw | Solitons with fermion number[END_REF][START_REF] Steinberg | Charge fractionalization in quantum wires[END_REF], have been prominent topics in theroretical and experimental condensed matter physics for over thirty years.

The subtleties of the inter-dependence between these concepts has 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 realms of frustrated magnetism. Indeed similar behaviour is observed in water ice [START_REF] Bernal | A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions[END_REF] and models of heavy fermion behaviour in spinels [START_REF] Fulde | On the origin of heavy quasiparticles in LiV2O4[END_REF], as well as in models of binary and mixed valence alloys [START_REF] Banks | Magnetic frustration in the context of pseudo-dipolar ionic disorder[END_REF].

Henley [START_REF] Henley | The Coulomb Phase in Frustrated Systems[END_REF] 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 interact via a Coulomb potential [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF] 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. 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 divergence-free part corresponding to the emergent Coulomb phase with independent and ergodic spin fluctuations [START_REF]The Coulomb phase can remain ergodic, either thanks to deconfined topological excitations[END_REF]. Our results apply even for a monopole crystal which 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 moment 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 [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF] emerge as quasi-particle excitations from the ground state configurations of the dumbbell model of spin ice. Here, the point dipoles of the dipolar spin ice model [START_REF] Hertog | Dipolar interactions and origin of spin ice in Ising pyrochlore magnets[END_REF] are extended to infinitesimally thin magnetic needles lying along the axes linking the centres of adjoining tetrahedra of the pyrochlore lattice (see Fig. 1) which constitute a diamond lattice of nodes for magnetic charge [START_REF] Ryzhkin | Magnetic relaxation in rare-earth oxide pyrochlores[END_REF]. The needles touch at the diamond lattice sites so that, by construction, the long range part 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 [START_REF] Isakov | Why spin ice obeys the ice rules[END_REF]. 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 neighbouring 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 [START_REF] Isakov | Dipolar spin correlations in classical pyrochlore magnets[END_REF], the Coulomb phase, with manifest experimental consequences for spin ice materials. These include diffuse neutron scattering patterns showing the sharp pinch point features [START_REF] Fennell | Magnetic Coulomb Phase in the Spin Ice Ho2Ti2O7[END_REF] characteristic of dipolar correlations, and the generation of large internal magnetic fields despite the status of vacuum [START_REF] Sala | Magnetic Coulomb Fields of Monopoles in Spin Ice and Their Signatures in the Internal Field Distribution[END_REF][START_REF] Dunsiger | Spin Ice: Magnetic Excitations without Monopole Signatures Using Muon Spin Rotation[END_REF].

In an analogy with electrostatics, the monopole charge distribution obeys Gauss' law for the magnetic field, ∇. 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 double charged monopoles (see Appendix A).

The divergence in H is related to the breaking of the ice rules through M , the magnetic moment density which itself obeys Gauss' law, ∇. M = -ρ m . This does not however, define the entire moment density which can have two contributions: the first, M m , is the gradient of a scalar potential and provides the magnetic charge; the second, M d , a dipolar field, can be represented as the curl of a vector potential and is divergence-free [START_REF] Maggs | Local Simulation Algorithms for Coulomb Interactions[END_REF][START_REF] Bramwell | Generalized longitudinal susceptibility for magnetic monopoles in spin ice[END_REF]

M = M m + M d = ∇ψ(r) + ∇ ∧ Q. ( 1
)
In states obeying the ice rules, M m = 0, while M d = M ; this corresponds to an emergent Coulomb phase [START_REF] Henley | The Coulomb Phase in Frustrated Systems[END_REF].

Breaking the ice rules and introducing magnetic charge leads to the conversion of M from the divergence-free M d to the divergence-full field M m . Central to the discussions in this paper is the concept that this conversion is, in general, partial -the divergence-free part is not completely suppressed by the presence of the monopoles. Thus, in almost all circumstances, one has the coexistence of two complimentary fields, reminiscent of the Hamiltonian splitting in topologically ordered phases [START_REF] Moessner | Irrational Charge from Topological Order[END_REF]. The total charge inside a microscopic volume surrounding a diamond lattice site i is Q i = -M .d S with d S an outward pointing surface element. For the dumbbell model the integral reduces to a discrete sum

Q i = -j M ij where M ij = M .d S ij = ±m/a, d S 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-one-out (three-out-one-in) gives Q i = +(-)2m/a [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF].

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

|M ij | = m/a: [M ij ] a m = (-1, -1, -1, 1) (2) = - 1 2 , - 1 2 , - 1 2 , - 1 2 + - 1 2 , - 1 2 , - 1 2 , 3 2 .
The first set of fields satisfy Gauss' law for the charge at the origin; the second set satisfy a discrete divergencefree condition and constitute a residual dipolar field dressing the monopole (see Fig. 1c). Decomposition further away from the charge could be made by solving the lattice Poisson equation to find the field sets belonging to M m which would be subtracted from the M ij to find the discrete elements of M d . Singly charged monopoles leave a residual contribution to M d at each vertex and it is only when the vertex is occupied by a double charged monopole for, which [M ij ] = ±(m/a)(1, 1, 1, 1), that the contribution to 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 [START_REF] Fennell | Magnetic Coulomb Phase in the Spin Ice Ho2Ti2O7[END_REF][START_REF] Sen | Coulomb Phase Diagnostics as a Function of Temperature, Interaction Range, and Disorder[END_REF], indicating the presence of such a dipolar field. The emergence of the dipolar field is further illustrated in the Appendix C where we show how the [M ij ] are divided around a pair of isolated nearest neighbour charges of the same sign. The random fluctuations in the underlying gauge field can be 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 sub lattices. 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 neighbouring sites on opposite sub-lattices. This fragmented state could be termed an "Antiferromagnetic Coulomb Magnet" (AFCM) by analogy with the "Ferromagnetic Coulomb Magnet" (FCM), predicted in the gauge mean field theory of quantum magnets on a pyrochlore lattice [START_REF] Savary | Coulombic Quantum Liquids in Spin-1/2 Pyrochlores[END_REF]. 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 part, M d , placing a dimer along the bond of strength 3/2 provides a mapping between The sub-lattice index ∆i defined in equation ( 4) is +1 (resp. -1) on the diamond sites where the blue (resp. red) monopoles sit. The minority spins (red) 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 magnetic moment fractionalization into divergence-full and divergence-free elements. We shall emphasize that the divergence-free field emerging here, despite being a Coulomb phase, is different from the standard one in absence of monopoles described in [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF].

the emergent dipolar field and hard core dimers on the (bipartite) diamond lattice [START_REF] Huse | Coulomb and Liquid Dimer Models in Three Dimensions[END_REF]. The extensive entropy of the dipolar field is thus associated with closed loops of dimer moves [START_REF] Alet | Interacting Classical Dimers on the Square Lattice[END_REF]. Introducing quantum loop dynamics give rise to a U(1) liquid phase close to the Rokhsar-Kivelson point [START_REF] Bergman | Ordering in a frustrated pyrochlore antiferromagnet proximate to a spin liquid[END_REF][START_REF] Sikora | Quantum Liquid with Deconfined Fractional Excitations in Three Dimensions[END_REF].

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 U C = (N 0 /2)αu, where u = -µ 0 Q 2 /4πa is the Coulomb energy for a nearest neighbour pair of monopoles of charge ±Q, µ 0 is the permeability, N 0 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

µ * < µ * 0 = α 2 = 0.819. ( 3
)
Monopole crystallization has also been studied recently within the dipolar spin ice model, in the canonical ensemble, that is, with fixed monopole number [START_REF] Borzi | Charge Ordering in a Pure Spin Model: Dipolar Spin Ice[END_REF], 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 [START_REF] Hertog | Dipolar interactions and origin of spin ice in Ising pyrochlore magnets[END_REF], at which point double charges become favoured. The ordering is then to a structure in complete analogy with zinc blende: a crystal of doubly charged monopoles for which M d is everywhere zero, corresponding to the "allin-all-out" magnetic order observed in FeF 3 [START_REF] Reimers | Critical properties of highly frustrated pyrochlore antiferromagnets[END_REF].

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. [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF]). To put spin ice materials in the context of the present work we note that the dumbbell approximation for dysprosium titanate yields µ = -4.35K, while direct simulations for dipolar spin ice give µ = -4.46K [START_REF] Jaubert | Magnetic Monopole Dynamics in Spin Ice[END_REF], so that µ * ≈ 1.42, well away from the monopole crystal phase boundary. In the appendix 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 [START_REF] Melko | Monte Carlo studies of the dipolar spin ice model[END_REF]. 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 four compared with values for the monopoles discussed in the main text. Hence, the zero temperature 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 from traditional spin ice modeling for large monopole concentrations; a point we return to in section 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 the Appendix B). In Figure 2 we show the evolution of an order parameter for monopole crystallization and the monopole number density, M c and n respectively, as a function of reduced temperature T * = T /|u| for different values of µ * . M c is defined as

M c = 1 N 0 N0 i=1 q i ∆ i , (4)
where 1), ... denotes a statistical average and n = N /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 tri-critical point as µ * increases [START_REF] Kobelev | Lattice models of ionic systems[END_REF]. Our data is consistent with this despite the additional constraints of the dumbbell model. From the finite size scaling analysis shown in the Appendix E we estimate a tri-critical point, µ * tr ≈ 0.78, T * tr ≈ 0.13. This temperature is comparable to that obtained from numerical simulation of a cubic lattice Coulomb gas [START_REF] Dickman | Phase diagram of the lattice restricted primitive model[END_REF], 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 towards 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 favoured, 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 first order, as in the Blume-Capel model for spin-one systems [START_REF] Dickman | Phase diagram of the lattice restricted primitive model[END_REF]. More work is required to extract the effect of the constraints in detail and to establish the tri-critical parameters with precision. There could, in principle also be a liquid-gas transition at higher temperature, between a low and high density fluid, but at the level of Debye-Hückel theory this transition is suppressed by the monopole ordering [START_REF] Kobelev | Lattice models of ionic systems[END_REF]. There is no strong evidence of this in our simulations, although the crossover from high density fluid (region II of the phase diagram in Fig. 3) to a low density fluid (region III) just outside the monopole crystal phase boundary is quite sharp and could be considered as a vestige of such a transition.

q i = (Q i a/2m) = {-1, 0, +1} is the topological charge on site i, ∆ i = ±1 is a diamond sub-lattice in- dex (see Fig.
In Figure 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 long range interactions. The transition appears to pass from second order to first order via a tri-critical 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 black open circles). The alternating positive/negative charges on the diamond lattice can also be seen as stacked monolayers of same charge monopoles 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).

In Figure 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 were 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 [START_REF] Reimers | Critical properties of highly frustrated pyrochlore antiferromagnets[END_REF]. The intensity of these peaks is precisely one quarter of that expected for complete allin-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 whilst the length of the majority spins remain fixed at unity. The resulting structure factor has no Bragg peaks and has four 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 equation [START_REF] Villain | Spin glass with non-random interactions[END_REF].

Magnetic charge crystallization also occurs for the dumbbell model on the kagome lattice (see Fig. 1) [START_REF] Möller | Magnetic multipole analysis of kagome and artificial spin-ice dipolar arrays[END_REF][START_REF] Chern | Two-Stage Ordering of Spins in Dipolar Spin Ice on the Kagome Lattice[END_REF][START_REF] Chern | Magnetic charge and ordering in kagome spin ice[END_REF][START_REF] Macdonald | Classical topological order in kagome ice[END_REF][START_REF] Cépas | Heterogeneous freezing in a geometrically frustrated spin model without disorder: Spontaneous generation of two time scales[END_REF]. Breaking a Z 2 symmetry for needle orientations on each triangular vertex plunges the system from the socalled KI phase into a constrained Coulomb phase, the KII phase, realized in spin ice materials by applying a field along the (111) direction of the pyrochlore lattice [START_REF] Matsuhira | A new macroscopically degenerate ground state in the spin ice compound Dy2Ti2O7 under a magnetic field[END_REF][START_REF] Fennell | Pinch points and Kasteleyn transitions in kagome ice[END_REF]. The moments 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:

[M ij ] a m = (-1, -1, 1) (5) = - 1 3 , - 1 3 , - 1 3 + - 2 3 , - 2 3 , 4 3 .
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 (Figure 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 [START_REF] Möller | Magnetic multipole analysis of kagome and artificial spin-ice dipolar arrays[END_REF][START_REF] Chern | Two-Stage Ordering of Spins in Dipolar Spin Ice on the Kagome Lattice[END_REF][START_REF] Chern | Magnetic charge and ordering in kagome spin ice[END_REF][START_REF] Fulde | Fractional charges in pyrochlore lattices[END_REF].

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 non-local worm algorithm [START_REF] Jaubert | Three-Dimensional Kasteleyn Transition: Spin Ice in a [100] Field[END_REF] extended to include long range interactions [START_REF] Melko | Long-Range Order at Low Temperatures in Dipolar Spin Ice[END_REF]. While the non-local 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 [START_REF] Jaubert | Magnetic Monopole Dynamics in Spin Ice[END_REF][START_REF] Jaubert | Signature of magnetic monopole and Dirac string dynamics in spin ice[END_REF]. The one-site, two-time monopole and spin autocorrelation functions are defined as

C c (t) = 1 nN 0 N0 i=1 q i (0) q i (t) (6)
C s (t) = 1 2N 0 2N0 j=1 S j (0). S j (t) . (7)
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 Figure [START_REF] Jackiw | Solitons with fermion number[END_REF]. The spin autocorrelation function, however, shows a decay over a modest simulation time, from unity to an asymptote of C s (t = ∞) = 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) is accessible experimentally as it provides the diagonal contribution to the a.c 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 [START_REF] Alet | Interacting Classical Dimers on the Square Lattice[END_REF][START_REF] Rokhsar | Superconductivity and the Quantum Hard-Core Dimer Gas[END_REF] 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 neighbour 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 to obscure the phase transition from bulk magnetic measurements. In Figure 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 towards the tri-critical 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 tem-perature phase to the low temperature Coulomb phase [START_REF] Ryzhkin | Magnetic relaxation in rare-earth oxide pyrochlores[END_REF][START_REF] Isakov | Magnetization curve of spin ice in a [111] magnetic field[END_REF][START_REF] Jaubert | Topological-Sector Fluctuations and Curie-Law Crossover in Spin Ice[END_REF]. 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 Figure 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 first order 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 towards 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. ) and the t = ∞ asymptote of Cs(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.

A. Magnetic field

An external magnetic field couples to both the monopoles and to M d providing a gradient to the chemical potential, inducing transient monopole currents and ordering the dipolar field (the two processes being intimately related [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF][START_REF] Ryzhkin | Magnetic relaxation in rare-earth oxide pyrochlores[END_REF][START_REF] Jaubert | Signature of magnetic monopole and Dirac string dynamics in spin ice[END_REF][START_REF] Giblin | Creation and Measurement of Long-Lived Magnetic Monopole Currents in Spin Ice[END_REF]). 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 sub-lattices, allowing experimental access to monopole crystallization [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF]. This 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 ) [START_REF] Sakakibara | Observation of a Liquid-Gas-Type Transition in the Pyrochlore Spin Ice Compound Dy 2 Ti2O7 in a Magnetic Field[END_REF]. If one were to start from the charge ordered Coulomb phase in zero field, an external (111) field would couple uniquely to M 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 [START_REF] Moessner | Theory of the [111] magnetization plateau in spin ice[END_REF]. An external magnetic field can also stabilize single charge monopole in Tb 2 Ti 2 O 7 [START_REF] Sazonov | Double-layered monopolar order in the Tb2Ti2O7 spin liquid[END_REF].

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 [START_REF] Zhou | High pressure route to generate magnetic monopole dimers in spin ice[END_REF].

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, entropy balance of the monopole crystal could possibly stabilize it ahead of either the monopole fluid, or the fully ordered double monopole 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 [START_REF] Savary | Coulombic Quantum Liquids in Spin-1/2 Pyrochlores[END_REF][START_REF] Ross | Quantum Excitations in Quantum Spin Ice[END_REF][START_REF] Bovo | Brownian motion and quantum dynamics of magnetic monopoles in spin ice[END_REF]. For systems close to the spin ice -antiferromagnetic phase boundary [START_REF] Hertog | Dipolar interactions and origin of spin ice in Ising pyrochlore magnets[END_REF] one might hope that zero point fluctuations of the fragmented dipolar field could stabilize the monopole crystal over the classical all-in-allout spin structure or the spin ice manifold [START_REF] Bergman | Ordering in a frustrated pyrochlore antiferromagnet proximate to a spin liquid[END_REF][START_REF] Sikora | Quantum Liquid with Deconfined Fractional Excitations in Three Dimensions[END_REF]. 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 111 ) 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, favouring monopole crystallization.

C. Tb2Ti2O7

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 [START_REF] Gardner | Cooperative paramagnetism in the geometrically frustrated pyrochlore antiferromagnet Tb2Ti2O7[END_REF] and quantum spin ice [START_REF] Molavian | Dynamically induced frustration as a route to a quantum spin ice state in Tb2Ti2O7 via virtual crystal field excitations and quantum many-body effects[END_REF] -or maybe spanning both. It is this dual nature that makes it an interesting study case for magnetic moment fragmentation. Tb 2 Ti 2 O 7 has a negative Curie-Weiss temperature, Θ CW = -14 K [START_REF] Gingras | Thermodynamic and single-ion properties of Tb 3+ within the collective paramagnetic-spin liquid state of the frustrated pyrochlore antiferromagnet Tb2Ti2O7[END_REF]. 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 [START_REF] Hertog | Dipolar interactions and origin of spin ice in Ising pyrochlore magnets[END_REF]. At ambient pressure it fails to develop magnetic order down to at least 50 mK [START_REF] Gardner | Dynamic frustrated magnetism in Tb2Ti2O7 at 50 mK[END_REF]. Recently, diffuse neutron scattering from single crystal samples have exposed pinch point scattering patterns indicating the presence of Coulomb phase correlations [START_REF] Fennell | Power-Law Spin Correlations in the Pyrochlore Antiferromagnet Tb2Ti2O7[END_REF][START_REF] Petit | Spin liquid correlations, anisotropic exchange, and symmetry breaking in Tb2Ti2O7[END_REF], albeit somewhat deformed compared with those observed for classical spin ice materials [START_REF] Fennell | Magnetic Coulomb Phase in the Spin Ice Ho2Ti2O7[END_REF].

This picture is however very sensitive to perturbations. Under high pressure the liquid-like phase transforms into a partially ordered structure where each vertex has one spin along a spin ice axis and three collinear spins [START_REF] Mirebeau | Pressure-induced crystallization of a spin liquid[END_REF]. Vertices on the two diamond sub-lattices 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" [START_REF] Legl | Vibrating-Coil Magnetometry of the Spin Liquid Properties of Tb2Ti2O7[END_REF]. A strong field in the [110] direction stabilizes a double-layered structure of singly charged monopoles [START_REF] Sazonov | Double-layered monopolar order in the Tb2Ti2O7 spin liquid[END_REF][START_REF] Ruff | Superlattice correlations in Tb2Ti2O7 under the application of [110] magnetic field[END_REF]. Recent experiments on polycrystalline samples have also shown that a very small amount of Tb 3+ stuffing produces a yet unexplained long range order, accompanied by weak antiferromagnetism [START_REF] Taniguchi | Long-range order and spin-liquid states of polycrystalline Tb2+xTi2-xO7+y[END_REF]. Finally, both strong magneto-elastic coupling linked to quantum spin liquid behaviour [START_REF] Ruff | Magnetoelastics of a Spin Liquid: X-Ray Diffraction Studies of Tb2Ti2O7 in Pulsed Magnetic Fields[END_REF][START_REF] Guitteny | Anisotropic Propagating Excitations and Quadrupolar Effects in Tb2Ti2O7[END_REF] and splitting of the single ion ground state doublet [START_REF] Petit | Spin liquid correlations, anisotropic exchange, and symmetry breaking in Tb2Ti2O7[END_REF][START_REF] Bonville | Tetragonal distortion yielding a two-singlet spin liquid in pyrochlore Tb2Ti2O7[END_REF][START_REF] Gaulin | Lack of evidence for a singlet crystal-field ground state in the magnetic pyrochlore Tb2Ti2O7[END_REF][START_REF] Fritsch | Antiferromagnetic spin ice correlations at ( 1 2 , 1 2 , 1 2 ) in the ground state of the pyrochlore magnet Tb2Ti2O7[END_REF] have been reported.

It would be naive to suggest a quantitative connection between the rich behaviour 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 [START_REF] Enjalran | The spin liquid state of the Tb2Ti2O7 pyrochlore antiferromagnet: a puzzling state of affairs[END_REF][START_REF] Mirebeau | Pressure and Field Induced Magnetic Order in the Spin Liquid Tb2Ti2O7 as Studied by Single Crystal Neutron Diffraction[END_REF]. 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 sub-lattice 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 freedom, M d . The incipient, or partial, all-in-all-out behaviour in the presence of a field along [111] is exactly what one would expect of our model when sat close to the monopole crystal 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 [START_REF] Fennell | Field-induced partial order in the spin ice dysprosium titanate[END_REF][START_REF] Hiroi | Ferromagnetic Ising spin chains emerging from the spin ice under magnetic field[END_REF], 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. Spinlattice coupling could be particularly relevant in this context [START_REF] Sazonov | Double-layered monopolar order in the Tb2Ti2O7 spin liquid[END_REF][START_REF] Fennell | Magnetoelastic excitations in the pyrochlore spin liquid Tb2Ti2O7[END_REF].

Turning to dynamics, the freezing observed in Tb 2 Ti 2 O 7 [START_REF] Gardner | Dynamic frustrated magnetism in Tb2Ti2O7 at 50 mK[END_REF][START_REF] Yasui | Static correlation and dynamical properties of Tb 3+ -moments in Tb2Ti2O7 -Neutron scattering study[END_REF][START_REF] Lhotel | Low-temperature magnetization in geometrically frustrated Tb2Ti2O7[END_REF] only involves a fraction of the spins ( 10%) and has been shown to be different from spin glass physics [START_REF] Lhotel | Low-temperature magnetization in geometrically frustrated Tb2Ti2O7[END_REF]. 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 magnetic moment fragmentation as a promising route to understanding apparent co-existence in this material of antiferromagnetism with the fluctuating Coulomb phase physics of a frustrated ferromagnet [START_REF] Fennell | Power-Law Spin Correlations in the Pyrochlore Antiferromagnet Tb2Ti2O7[END_REF].

D. Other materials

A second quantum spin ice candidate is Yb 2 Ti 2 O 7 , [START_REF] Ross | Quantum Excitations in Quantum Spin Ice[END_REF][START_REF] Hayre | Thermodynamic properties of Yb2Ti2O7 pyrochlore as a function of temperature and magnetic field: Validation of a quantum spin ice exchange Hamiltonian[END_REF]. This material, with a Curie-Weiss temperature estimated at around 600 mK [START_REF] Bramwell | Bulk magnetization of the heavy rare earth titanate pyrochlores -a series of model frustrated magnets[END_REF], shows an unusual phase transition at 200 mK, [START_REF] Hodges | First-Order Transition in the Spin Dynamics of Geometrically Frustrated Yb2Ti2O7[END_REF] with apparently no accompanying magnetic order [START_REF] Hodges | First-Order Transition in the Spin Dynamics of Geometrically Frustrated Yb2Ti2O7[END_REF][START_REF] Gardner | Spin-spin correlations in Yb2Ti2O7: A polarized neutron scattering study[END_REF] and magnetic dimensional reduction [START_REF] Ross | Dimensional evolution of spin correlations in the magnetic pyrochlore Yb2Ti2O7[END_REF] in the high temperature phase. There are, however, reports of a (partial) ferromagnetic ordering [START_REF] Ross | Dimensional evolution of spin correlations in the magnetic pyrochlore Yb2Ti2O7[END_REF][START_REF] Chang | Higgs transition from a magnetic Coulomb liquid to a ferromagnet in Yb2Ti2O7[END_REF] at ∼ 400 mK. The magnetic anisotropy in Yb 2 Ti 2 O 7 was initially considered to be XY like [START_REF] Hodges | First-Order Transition in the Spin Dynamics of Geometrically Frustrated Yb2Ti2O7[END_REF] but more recent analysis has suggested that it could in fact be considered as a spin ice with the low temperature behaviour experimentally close to a quantum spin liquid [START_REF] Ross | Quantum Excitations in Quantum Spin Ice[END_REF][START_REF] Hayre | Thermodynamic properties of Yb2Ti2O7 pyrochlore as a function of temperature and magnetic field: Validation of a quantum spin ice exchange Hamiltonian[END_REF]. Within this context, a quantum spin liquidclassical spin gas transition has recently been proposed [START_REF] Savary | Spin liquid regimes at nonzero temperature in quantum spin ice[END_REF]. With such complex behaviour and sample dependence [START_REF] Ross | Lightly stuffed pyrochlore structure of singlecrystalline Yb2Ti2O7 grown by the optical floating zone technique[END_REF][START_REF] D'ortenzio | Unconventional Magnetic Ground State in Yb2Ti2O7[END_REF], a quantitative understanding requires a detailed microscopic approach [START_REF] Ross | Quantum Excitations in Quantum Spin Ice[END_REF][START_REF] Hayre | Thermodynamic properties of Yb2Ti2O7 pyrochlore as a function of temperature and magnetic field: Validation of a quantum spin ice exchange Hamiltonian[END_REF][START_REF] Han | [END_REF]. That being said, the joint features of a low temperature magnetically fluctuating phase and the presence of a phase transition which is partially transparent to magnetic probes is 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 [START_REF] Mirebeau | Ordered Spin Ice State and Magnetic Fluctuations in Tb2Sn2O7[END_REF][START_REF] Dahlberg | Low-temperature dynamic freezing and the fragility of ordering in Tb2Sn2O7[END_REF], as predicted by combining dipolar interactions and spin relaxation [START_REF] Mcclarty | Soft dipolar spin ice physics and the ordered phase of the frustrated Tb2Sn2O7 pyrochlore magnet[END_REF]. 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 (equation 5). Crystallites of the KII phase have recently been realized in permalloy nano-arrays with a honeycomb structure [START_REF] Zhang | Crystallites of magnetic charges in artificial spin ice[END_REF]. 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 Figure 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 field along the (1, 1, 1) direction, although they may be masked by the field induced 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 [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF] which only partially maps onto the physical degrees of freedom, the magnetic needles. As a consequence, the intrinsic moments fragment into two parts. The first 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 which 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 [START_REF] Legl | Vibrating-Coil Magnetometry of the Spin Liquid Properties of Tb2Ti2O7[END_REF][START_REF] Hodges | First-Order Transition in the Spin Dynamics of Geometrically Frustrated Yb2Ti2O7[END_REF][START_REF] Mirebeau | Ordered Spin Ice State and Magnetic Fluctuations in Tb2Sn2O7[END_REF][START_REF] Petrenko | Investigation of the Low-Temperature Spin-Liquid Behavior of the Frustrated Magnet Gadolinium Gallium Garnet[END_REF][START_REF] Lee | Local Spin Resonance and Spin-Peierlslike Phase Transition in a Geometrically Frustrated Antiferromagnet[END_REF]). 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 which can also occur in experiments on rare earth pyrochlores [START_REF] Ross | Dimensional evolution of spin correlations in the magnetic pyrochlore Yb2Ti2O7[END_REF][START_REF] Gardner | Magnetic pyrochlore oxides[END_REF][START_REF] Cao | Anisotropic exchange in frustrated pyrochlore Yb2Ti2O7[END_REF]. 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 [START_REF] Savary | Coulombic Quantum Liquids in Spin-1/2 Pyrochlores[END_REF][START_REF] Hermele | Pyrochlore photons: The U(1) spin liquid in a S=1/2 three-dimensional frustrated magnet[END_REF][START_REF] Benton | Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice[END_REF] have revealed a complete model for quantum electrodynamics (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.

Appendix A: The Coulomb Physics of Spin Ice

The dipolar spin ice Hamiltonian can be written (see SI of [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF]), within the dumbbell approximation as

H -H 0 = 1 2 i =j µ 0 Q i Q j 4πr ij + 1 2 v 0 i Q 2 i , (A1)
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:

1 2 v 0 i Q 2 i = -µN -µ 2 N 2 , ( A2
)
where µ = -v 0 Q 2 /2, µ 2 = -2v 0 Q 2 and where the number of single and double monopoles are N and N 2 respectively. The sketch below illustrates the energy scale taking a single vertex, or 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. 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 zero temperature phase transition from monopole vacuum to monopole crystal, equation (3), holds for both single and double monopole crystals, giving in both cases

v 0 = µ 0 α 4πa . ( A3
)
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 [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF]

|µ| = v 0 Q 2 2 = - 2J 3 + 8 3 1 + 2 3 D (
J nn D nn = - 4 5 1 + 2 3 1 - α 2 = -0.918, (A6)
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 [START_REF] Melko | Monte Carlo studies of the dipolar spin ice model[END_REF] find J nn /D nn = -0.905 with hysteresis down to J nn /D nn ≃ -1, the origin of the difference being the small band width for the Pauling states which 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 which 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 [START_REF] Applegate | Vindication of Yb2Ti2O7 as a Model Exchange Quantum Spin Ice[END_REF], thermal fluctuations, or a staggered, sub lattice dependent chemical potential. This point is addressed further in section IV B, where we discuss the relevance of our work to experiment.

For the results in Figure 6, we used 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 were 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 equation ( 5) was 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 was 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 were made every 10 MCS/s. Fifty worm updates were performed every 10 MCS/s to facilitate thermalization. We ran 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 Figure 2, the parameters are L = 8, t eq = 10 4 MCS/s, with an observation period of 10 5 MCS/s and averaging over 4 independent simulations. Again, 50 worm updates were performed every 10 MCS/s; 100 different temperatures between 0.2 and 0.6 K were 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 which:

• does not modify the number of charges or their positions and

• 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-one-out), 2 (two-in-two-out) or 3 (three-out-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 whilst 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-two-out), the dimer covering of the diamond lattice (alternating three-in-one-out and three-out-one-in) and for all densities of monopoles in between.

This worm algorithm could also be developed for the dipolar spin ice model [START_REF] Hertog | Dipolar interactions and origin of spin ice in Ising pyrochlore magnets[END_REF]. 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. The susceptibility of the monopole crystallization order parameter is defined as

χ c = M 2 c -M c 2 N 0 T (E1)
In order to characterize more quantitatively the nature of the phase transitions, we performed finite size scaling on the maxima of the specific heat C µ and susceptibility χ c , plotted in Figure 9. As explained in the paper, two regimes clearly appear. The transition is continuous up to a tri-critical point at µ * tr = 0.78 ± 0.01 where it becomes 1 st 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 [START_REF] Kobelev | Lattice models of ionic systems[END_REF] which would require further numerical and/or theoretical effort.

The discontinuity of the order parameter in Figure 2 of the main text strongly supports the 1 st order nature of the phase transition, but its quantitative signature in finite size scaling is rather challenging. That is, the diverging correlation length either side of the tri-critical point could lead to over estimates of µ * tr and T * tr . The 1 st order regime in our system occupies only a small region of parameter space making it difficult to separate first order from tri-critical behaviour. Increasing the system size beyond the correlation length rapidly becomes very time consuming -the 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 1 st order transitions, especially in large systems. Nonetheless, it has been possible to show a sharp increase of the scaling exponents close to the low temperature phase boundary, their values approaching those of a 1 st order transition (see lower panels of Figure 9). In particular, the maximum in C µ develops scaling behaviour in this regionclose 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 Figure 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 1 st order region. This suggests that spin ice materials and models close enough to the low temperature phase boundary can exhibit correlations inherited from the monopole crystallization, which should be visible using local probe such as neutron scattering. as a function of linear system size L, for µ * = 0.490 (+), 0.654 (×), 0.768 ( ), 0.778 ( ), 0.784 (•), 0.794 (•), 0.795 (△), 0.796 ( ), 0801 ( ). The error bars are the standard deviation σ over 4 independent simulation outcomes. Each solid line is the best fit obtained using "Wolfram Mathematica v9.0" [START_REF] Research | Mathematica Edition: Version 9.0[END_REF], including all data points for a given µ * and weighting each data point by 1/σ 2 . The dashed line is a guide-to-the-eye for the cubic power law (∝ L 3 ) appearing in the 1 st order regime for L 4. Bottom: Scaling Exponents ratio α/ν and γ/ν as a function of µ * . The error bars represent a confidence level of 90%, based on the statistical uncertainty of the data plotted in the top panels.

FIG. 1 .

 1 FIG. 1. (Colour) Lattice structures: (a) Pyrochlore lattice, showing monopole crystal and magnetic needle (spin) orientations.The sub-lattice index ∆i defined in equation (4) is +1 (resp. -1) on the diamond sites where the blue (resp. red) monopoles sit. The minority spins (red) 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 magnetic moment fractionalization into divergence-full and divergence-free elements. We shall emphasize that the divergence-free field emerging here, despite being a Coulomb phase, is different from the standard one in absence of monopoles described in[START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF].

FIG. 2 .

 2 FIG. 2. (Color) (a) Order parameterMc for monopole crystallization as defined in equation (4) and (b) density of charges n, as a function of reduced temperature T * for chemical potential µ * = 0.767(+), 0.778(×), 0.784( ), 0.794( ), 0.801(•). The transition appears to pass from second order to first order via a tri-critical 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 black open circles). The alternating positive/negative charges on the diamond lattice can also be seen as stacked monolayers of same charge monopoles 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).

FIG. 3 .

 3 FIG. 3. (Color) 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 tri-critical point. Index I, II and III mark the monopole crystal, high density 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 .

 4 FIG. 4. (Color) Simulated unpolarised 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 grayscale superimposed on the contribution to S(Q) from the dipolar field. The wavevector Q is in units of 2π/a0 where a0 = 4a/ √ 3 is the lattice parameter of the cubic unit cell of the pyrochlore lattice. The kagome ice data is taken from[START_REF] Harman-Clarke | Topological constraints and ordering in model frustrated magnets[END_REF].

FIG. 5 .

 5 FIG. 5. (Color) Charge and spin autocorrelation functions for µ * = 0.57 (top) and µ * = 0.41 (bottom) for T * = 0.20.The dotted line in each plot indicates the asymptote at Cs(t = ∞) = 1/4. The time required for the autocorrelation function to reach its long-time asymptote is significantly greater at µ * = 0.41 than at µ * = 0.57. This is indicative of a slowing of the dynamics as the monopole density increases. As µ * increases, n decreases (see Fig.2(b)) and the t = ∞ asymptote of Cs(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.

FIG. 6 .

 6 FIG. 6. (Color) (a) Magnetic susceptibility χ as a function of reduced temperature T * for chemical potential µ * = 0.33(+), 0.65(×), 0.784(△), 0.794( ) < α/2, where α is the Madelung constant. The inverse susceptibility is plotted in the inset on a log-log scale. (b) Curie law crossover 3χT vs T (same color labeling with two additional values of µ * = 0.801( ), 0.98(•) > α/2). The dashed lines are theoretical expectations for the spin liquid Curie law prefactor C of the Coulomb phase (C ≈ 2) and the singly charged monopole fluid (C = 4/3), while the one for C = 1.52 is a guide to the eye. Equilibration has been ensured down to T = 20 mK (T * ≈ 0.007) by the worm algorithm (See Appendix B for more details).

  Neglecting the double monopoles, the internal energy, Ũ = U C -µN and we have the lattice Coulomb gas studied in the main text.

  A4) where J is the nearest neighbor (antiferromagnetic) exchange constant and D = µ 0 m 2 4πk 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 |u(a)| =

FIG. 7 .

 7 FIG. 7. (Color) 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.

FIG. 9 .

 9 FIG. 9. (Color)Top: Finite size scaling of the maxima of the specific heat Cµ and susceptibility χc of the order parameter Mc as a function of linear system size L, for µ * = 0.490 (+), 0.654 (×), 0.768 ( ), 0.778 ( ), 0.784 (•), 0.794 (•), 0.795 (△), 0.796 ( ), 0801 ( ). The error bars are the standard deviation σ over 4 independent simulation outcomes. Each solid line is the best fit obtained using "Wolfram Mathematica v9.0"[START_REF] Research | Mathematica Edition: Version 9.0[END_REF], including all data points for a given µ * and weighting each data point by 1/σ 2 . The dashed line is a guide-to-the-eye for the cubic power law (∝ L 3 ) appearing in the 1 st order regime for L 4. Bottom: Scaling Exponents ratio α/ν and γ/ν as a function of µ * . The error bars represent a confidence level of 90%, based on the statistical uncertainty of the data plotted in the top panels.

ACKNOWLEDGMENTS

We thank S. T. Bramwell for discussions concerning the simulated neutron scattering plots and issues related to reference [21]. It is also a pleasure to thank O. Benton, M. Faulkner, T. Fennell, M. J. P. Gingras, J. Gardner, V. Kaiser, A. C. Maggs, P. McClarty, R. Moessner and K. Penc for useful discussions. P.C.W.H. thanks the Institut Universitaire de France for financial support and is also grateful to the Yukawa Institute for Theoretical Physics, Kyoto for hospitality during the workshop NQS2011. S.T.B. thanks the Ecole Normale Supérieure de Lyon for support in the form of a visiting Professorship. M.E.B-B. thanks the STINT research fund from Uppsala Universitet for financial support.

Appendix B: Simulations

The numerical results in this paper were obtained by simulating the dumbbell model [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF], equivalent to a gas of singly charged magnetic monopoles. The energy scale of the Coulomb interactions was entirely determined by u(a), the energy scale for nearest neighbour monopoles. The chemical potential µ * is then a free parameter.

Three variants of Monte Carlo simulations were employed: (i) 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 and6); (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); (iii) SSF with worm updates and parallel tempering [START_REF] Swendsen | Replica Monte Carlo Simulation of Spin-Glasses[END_REF] (used for the results in Fig. 2).

The results presented in Figure 3 were 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 was 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.

Appendix C: Field distributions

As an example of magnetic moment fractionalization in the monopole fluid phase, we show, in Figure 7, two isolated neighbouring 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 [START_REF] Möller | Magnetic multipole analysis of kagome and artificial spin-ice dipolar arrays[END_REF] in this paper) as the fields can be easily visualized. Starting at 9 o'clock and turning clockwise the fields for sites 1 (on the left) and 2 (on the right) can be decomposed as follows:

where in each case the first term and second terms are the contributions to the divergence-full, M m , and divergencefree, M 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 Figure S8 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 re-establishing the ice rules on the four vertices of the square plaquet. These moves are followed by flipping of the horizontal needles which re-establishes 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 [START_REF] Alet | Interacting Classical Dimers on the Square Lattice[END_REF]. 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.