# Magnetic-Moment Fragmentation and Monopole Crystallization

M. E. Brooks-Bartlett, S. T. Banks, L. D. C. Jaubert, A. Harman-Clarke, P. C. W. Holdsworth

► To cite this version:

M. E. Brooks-Bartlett, S. T. Banks, L. D. C. Jaubert, A. Harman-Clarke, P. C. W. Holdsworth. Magnetic-Moment Fragmentation and Monopole Crystallization. Physical Review X, 2014, 4, pp.011007. 10.1103/PhysRevX.4.011007. hal-01541937

HAL Id: hal-01541937

https://hal.science/hal-01541937

Submitted on 21 Jun 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Magnetic moment fragmentation and monopole crystallization

M. E. Brooks-Bartlett,¹ S. T. Banks,¹ L. D. C. Jaubert,² A. Harman-Clarke,¹,³ and P. C. W. Holdsworth³, *

¹Department of Chemistry, University College London,

20 Gordon Street, London WC1H 0AJ, United Kingdom.

²OIST, Okinawa Institute of Science and Technology, Onna-son, Okinawa, 904-0495, Japan

³Laboratoire de Physique, École Normale Supérieure de Lyon,

Université de Lyon, CNRS, 46 Allée d'Italie, 69364 Lyon Cedex 07, France.

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. To date, 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 behaviour 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 quasi-particles in graphene-like systems [4] to solitons in one-dimension [5, 6], have been prominent topics in theoretical 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 [7] and models of heavy fermion behaviour 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 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. 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[12]. 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 [11] 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 [13] 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 [14]. 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 [15]. 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

<footnote>* peter.holdsworth@ens-lyon.fr</footnote>
this one is rather exceptional, as it is far from empty. Rather, the magnetic moments constitute the curl of a lattice gauge field [16], the Coulomb phase, with manifest experimental consequences for spin ice materials. These include diffuse neutron scattering patterns showing the sharp pinch point features [17] characteristic of dipolar correlations, and the generation of large internal magnetic fields despite the status of vacuum [18, 19].

In an analogy with electrostatics, the monopole charge distribution obeys Gauss' law for the magnetic field, $\vec{\nabla} \cdot \vec{H} = \rho_m$, where $\rho_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, $\vec{U} = U_C - \mu N$, where $U_C$ is the Coulomb energy of the neutral gas of $N$ = monopoles and $\mu$ is the chemical potential such that $-2\mu$ is the energy cost of introducing a neutral pair of monopoles an infinite distance apart. For temperatures of order $-\mu$ the Coulomb gas picture of spin ice must be extended to include double charged monopoles (see Appendix A).

The divergence in $\vec{H}$ is related to the breaking of the ice rules through $\vec{M}$, the magnetic moment density which itself obeys Gauss' law, $\vec{\nabla} \cdot \vec{M} = -\rho_m$. This does not however, define the entire moment density which can have two contributions: the first, $\vec{M}_m$, is the gradient of a scalar potential and provides the magnetic charge; the second, $\vec{M}_d$, a dipolar field, can be represented as the curl of a vector potential and is divergence-free [20, 21]

$$ \vec{M} = \vec{M}_m + \vec{M}_d = \vec{\nabla}\psi(r) + \vec{\nabla} \wedge \vec{Q}. \quad (1) $$

In states obeying the ice rules, $\vec{M}_m = 0$, while $\vec{M}_d = \vec{M}$; this corresponds to an emergent Coulomb phase [10]. Breaking the ice rules and introducing magnetic charge leads to the conversion of $\vec{M}$ from the divergence-free $\vec{M}_d$ to the divergence-full field $\vec{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 [22].

The total charge inside a microscopic volume surrounding a diamond lattice site $i$ is $Q_i = - \int \vec{M} \cdot d\vec{S}$ with $d\vec{S}$ an outward pointing surface element. For the dumbbell model the integral reduces to a discrete sum $Q_i = - \sum_j M_{ij}$ where $M_{ij} = \vec{M} \cdot d\vec{S}_{ij} = \pm m/a$, $d\vec{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$ [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 $|M_{ij}| = m/a$:

$$ \begin{aligned} [M_{ij}]\left(\frac{a}{m}\right) &= (-1, -1, -1, 1) \\ &= \left(-\frac{1}{2}, -\frac{1}{2}, -\frac{1}{2}, -\frac{1}{2}\right) + \left(-\frac{1}{2}, -\frac{1}{2}, -\frac{1}{2}, \frac{3}{2}\right). \end{aligned} \quad (2) $$

The first set of fields satisfy Gauss' law for the charge at the origin; the second set satisfy a discrete divergence-free 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 $\vec{M}_m$ which would be subtracted from the $M_{ij}$ to find the discrete elements of $\vec{M}_d$. Singly charged monopoles leave a residual contribution to $\vec{M}_d$ at each vertex and it is only when the vertex is occupied by a double charged monopole for, which $[M_{ij}] = \pm(m/a)(1, 1, 1, 1)$, that the contribution to $\vec{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 [17, 23], 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 [24]. 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, $\vec{M}_d$, placing a dimer along the bond of strength $3/2$ provides a mapping between
<figure><img src="image_2.png" /><figcaption>FIG. 1. (Colour) Lattice structures: (a) Pyrochlore lattice, showing monopole crystal and magnetic needle (spin) orientations. The sub-lattice index $\Delta_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 [11].</figcaption></figure>

the emergent dipolar field and hard core dimers on the (bipartite) diamond lattice [25]. The extensive entropy of the dipolar field is thus associated with closed loops of dimer moves [26]. Introducing quantum loop dynamics give rise to a U(1) liquid phase close to the Rokhsar-Kivelson point [27, 28].

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, $-\mu N$. For the monopolar crystal $U_C = (N_0/2)\alpha u$, where $u = -\mu_0 Q^2 / 4\pi a$ is the Coulomb energy for a nearest neighbour pair of monopoles of charge $\pm Q$, $\mu_0$ is the permeability, $N_0$ the number of diamond lattice sites and $\alpha = 1.638$ is the Madelung constant. We define a reduced chemical potential $\mu^* = \mu/u$ and thus the ground state should be a monopole crystal for

$$ \mu^* < \mu_0^* = \frac{\alpha}{2} = 0.819. \qquad (3) $$

Monopole crystallization has also been studied recently within the dipolar spin ice model, in the canonical ensemble, that is, with fixed monopole number [29], 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 [13], 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 $\vec{M}_d$ is everywhere zero, corresponding to the “all-in-all-out” magnetic order observed in FeF$_3$ [30].

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 $\mu = -4.35$K, while direct simulations for dipolar spin ice give $\mu = -4.46$K [31], so that $\mu^* \approx 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 [32]. 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 $\mu^*$. $M_c$ is defined as

$$M_c = \left\langle \left| \frac{1}{N_0} \sum_{i=1}^{N_0} q_i \Delta_i \right| \right\rangle, \quad (4)$$

where $q_i = (Q_ia/2m) = \{-1, 0, +1\}$ is the topological charge on site $i$, $\Delta_i = \pm 1$ is a diamond sub-lattice index (see Fig. 1), $\langle \dots \rangle$ denotes a statistical average and $n = \langle N \rangle / N_0$. The data show clear evidence of monopole crystallization at a transition temperature, $T_c^*$ that varies with $\mu^*$. 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 $(\mu^*, T)$ plane, becoming first order via a tri-critical point as $\mu^*$ increases [33]. 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, $\mu_{tr}^* \approx 0.78$, $T_{tr}^* \approx 0.13$. This temperature is comparable to that obtained from numerical simulation of a cubic lattice Coulomb gas [34], although $\mu_{tr}^*$ is surprisingly close to $\mu_0^*$. For small values of $\mu^*$ a continuous transition takes the system from a high density fluid ($n > 4/7$) to the crystalline phase. However, as $\mu^*$ 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 $\mu^* = 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 [34]. 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 [33]. 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.

In Figure 3 we show the resulting phase diagram, as mapped out by the divergence of the specific heat, $C_\mu$, at the phase transition. The monopole crystal phase terminates for $\mu^* \approx 0.8$ in approximate agreement with our prediction of $\mu_0^* = 0.819$, the small difference being most likely due to finite size effects exacerbated by the long range interactions.

<figure><img src="image_3.png" /><figcaption>FIG. 2. (Color) (a) Order parameter $M_c$ for monopole crystallization as defined in equation (4) and (b) density of charges $n$, as a function of reduced temperature $T^*$ for chemical potential $\mu^* = 0.767(+), 0.778(\times), 0.784(\square), 0.794(\blacksquare), 0.801(\circ)$. The transition appears to pass from second order to first order via a tri-critical point for $\mu_{tr}^* \approx 0.78$. The very narrow first order region ends at $\mu^* \approx 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).</figcaption></figure>

In Figure 4 we show a simulated elastic neutron scattering map determined within the static approximation,
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

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ᵢ 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₃ [30]. 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 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 (2).

Magnetic charge crystallization also occurs for the dumbbell model on the kagome lattice (see Fig. 1) [36–40]. Breaking a Z₂ symmetry for needle orientations on each triangular vertex plunges the system from the so-called 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 [41, 42]. 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:

$$ 
\begin{align} 
[M_{ij}] \left(\frac{a}{m}\right) &= (-1, -1, 1) \tag{5} \\ 
&= \left(-\frac{1}{3}, -\frac{1}{3}, -\frac{1}{3}\right) + \left(-\frac{2}{3}, -\frac{2}{3}, \frac{4}{3}\right). \nonumber 
\end{align} 
$$

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 [36–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 non-local worm algorithm [44] extended to include long range interactions [45]. 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 [31, 46]. The one-site, two-time monopole and spin autocorrelation functions are defined as

$$ C_c(t) = \left\langle \frac{1}{nN_0} \sum_{i=1}^{N_0} q_i(0) q_i(t) \right\rangle \quad (6) $$

$$ C_s(t) = \left\langle \frac{1}{2N_0} \sum_{j=1}^{2N_0} \vec{S}_j(0) \cdot \vec{S}_j(t) \right\rangle . \quad (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 (5). 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 [26, 47] if the system is to
<figure><img src="image_5.png" /><figcaption>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\pi/a_0$ where $a_0 = 4a/\sqrt{3}$ is the lattice parameter of the cubic unit cell of the pyrochlore lattice. The kagome ice data is taken from [35].</figcaption></figure>

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\tilde{U} \sim -u\alpha + 2\mu$ 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, $\chi$, as a function of temperature for different values of $\mu^*$. 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 temperature

phase to the low temperature Coulomb phase [14, 48, 49]. It was demonstrated that taking the needles as scatterers of unit length gives rise to a crossover from $C = 3T\chi = 1$ to $C \approx 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 \sim 3/2$ on entering the monopole crystal phase and again to $C \approx 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.
<figure><img src="image_6.png" /></figure>

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

A. Magnetic field

An external magnetic field couples to both the monopoles and to $\vec{M}_d$ providing a gradient to the chemical potential, inducing transient monopole currents and ordering the dipolar field (the two processes being intimately related [11, 14, 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 sub-lattices, allowing experimental access to monopole crystallization [11]. This corresponds to the first order phase transition observed as spin ice materials leave the plateau region of the $(\vec{B}, T)$ phase diagram. The corresponding phase boundary terminates at a critical end point, $(\vec{B}_C, T_C)$ [51]. If one were to start from the charge ordered Coulomb phase in zero field, an external (111) field would couple uniquely to $\vec{M}_d$. The

<figure><img src="image_7.png" /></figure>

FIG. 6. (Color) (a) Magnetic susceptibility $\chi$ as a function of reduced temperature $T^*$ for chemical potential $\mu^* = 0.33(+), 0.65(x), 0.784(\triangle), 0.794(\Box) < \alpha/2$, where $\alpha$ is the Madelung constant. The inverse susceptibility is plotted in the inset on a log-log scale.
(b) Curie law crossover $3\chi T$ vs $T$ (same color labeling with two additional values of $\mu^* = 0.801(\blacksquare), 0.98(\circ) > \alpha/2$). The dashed lines are theoretical expectations for the spin liquid Curie law prefactor $C$ of the Coulomb phase ($C \approx 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^* \approx 0.007$) by the worm algorithm (See Appendix B for more details).

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 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 $\pm 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 [24, 56, 57]. For systems close to the spin ice-antiferromagnetic phase boundary [13] 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 [27, 28]. 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 $\langle 111 \rangle$ 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. $Tb_2Ti_2O_7$

Fifteen years of intense research have made $Tb_2Ti_2O_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 study case for magnetic moment fragmentation. $Tb_2Ti_2O_7$ has a negative Curie-Weiss temperature, $\Theta_{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 [13]. At ambient pressure it fails to develop magnetic order down to at least 50 mK [61]. Recently, diffuse neutron scattering from single crystal samples have 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 [17].

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 [64]. 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” [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 a yet unexplained long range order, accompanied by weak antiferromagnetism [67]. Finally, both strong magneto-elastic coupling linked to quantum spin liquid behaviour [68, 69] and splitting of the single ion ground state doublet [63, 70–72] have been reported.

It would be naive to suggest a quantitative connection between the rich behaviour of $Tb_2Ti_2O_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 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, $\vec{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_2Ti_2O_7$ or $Dy_2Ti_2O_7$ [75, 76], the double layer structure of monopoles observed in $Tb_2Ti_2O_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_2Ti_2O_7$ [61, 78, 79] only involves a fraction of the spins ($\gtrsim 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_2Ti_2O_7$ is beyond the scope of this paper, we do propose magnetic moment fragmentation as a promising route to understanding 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_2Ti_2O_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 ~400 mK. The magnetic anisotropy in $Yb_2Ti_2O_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 behaviour experimentally close to a quantum spin liquid [56, 80]. Within this context, a quantum spin liquid – classical spin gas transition has recently been proposed [86]. With such complex behaviour 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 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_2Sn_2O_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 (equation 5). Crystallites of the KII phase have recently been realized in permalloy nano-arrays 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 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 [11] 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 [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 which can also occur in experiments on rare earth 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 [24, 98, 99] 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.

## 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. Ben-ton, M. Faulkner, T. Fennell, M. J. P. Gingras, J. Gardner, V. Kaiser, A. C. Maggs, P. McClarty, R. Moess-ner 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 Theoret-ical 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 Professor-ship. M.E.B-B. thanks the STINT research fund from Uppsala Universitet for financial support.
Appendix A: The Coulomb Physics of Spin Ice

The dipolar spin ice Hamiltonian can be written (see SI of [11]), within the dumbbell approximation as

$$H - H_0 = \frac{1}{2} \sum_{i \neq j} \frac{\mu_0 Q_i Q_j}{4\pi r_{ij}} + \frac{1}{2} v_0 \sum_i Q_i^2, \quad (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_0Q^2$, with $Q = 2m/a$ the monopole charge. The ice rules and their consequent violation impose that $Q_i = 0, \pm Q, \pm 2Q$ only and the diagonal term provides the chemical potentials for both singly ($\mu$) and doubly charged ($\mu_2$) monopoles:

$$\frac{1}{2} v_0 \sum_i Q_i^2 = -\mu N - \mu_2 N_2, \quad (A2)$$

where $\mu = -v_0Q^2/2$, $\mu_2 = -2v_0Q^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 two-in-two out ground state, to a three in-one out monopole and finally to an all-in or all-out double monopole.

<figure><img src="image_8.png" /></figure>

Neglecting the double monopoles, the internal energy, $\tilde{U} = U_C - \mu 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 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 = \frac{\mu_0 \alpha}{4\pi a}. \quad (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 [11]

$$|\mu| = \frac{v_0 Q^2}{2} = - \left[ \frac{2J}{3} + \frac{8}{3} \left( 1 + \sqrt{\frac{2}{3}} \right) D \right] \quad (A4)$$

where $J$ is the nearest neighbor (antiferromagnetic) exchange constant and $D = \frac{\mu_0 m^2}{4\pi 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)| = \frac{8}{3}\sqrt{\frac{2}{3}}D$, so that equation (3) becomes

$$\frac{2J}{3} + \frac{8}{3} \left( 1 + \sqrt{\frac{2}{3}} \right) D = \frac{\alpha}{2} \frac{8}{3} \sqrt{\frac{2}{3}} D, \quad (A5)$$

and hence

$$\frac{J_{nn}}{D_{nn}} = -\frac{4}{5} \left[ 1 + \sqrt{\frac{2}{3}} \left( 1 - \frac{\alpha}{2} \right) \right] = -0.918, \quad (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 [32] find $J_{nn}/D_{nn} = -0.905$ with hysteresis down to $J_{nn}/D_{nn} \simeq -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, $\mu_2 = 4\mu$ 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, sub lattice dependent chemical potential. This point is addressed further in section IV B, where we discuss the relevance of our work to experiment.

Appendix B: Simulations

The numerical results in this paper were obtained by simulating the dumbbell model [11], 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 $\mu^*$ 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_2Ti_2O_7$ and $Ho_2Ti_2O_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); (iii) SSF with worm updates and parallel tempering [101] (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.
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 $\mu^*$) after which observations of the autocorrelation function were made every MCS/s for a total of $2.5 \times 10^4$ MCS/s ($\mu^* = 0.57$) and every 100 MCS/s for a total of $10^6$ MCS/s ($\mu^* = 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 $\mu^*$; 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.

## 1. 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 [13]. 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 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 [36] 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:

$$ \begin{aligned} [M_{ij}]_1 \left(\frac{a}{m}\right) &= (-1, -1, 1, -1) & (C1) \\ &= \left(-1, -\frac{1}{2}, 0, -\frac{1}{2}\right) + \left(0, -\frac{1}{2}, 1, -\frac{1}{2}\right), \end{aligned} $$

$$ \begin{aligned} [M_{ij}]_2 \left(\frac{a}{m}\right) &= (-1, -1, -1, 1) & (C2) \\ &= \left(0, -\frac{1}{2}, -1, -\frac{1}{2}\right) + \left(-1, -\frac{1}{2}, 0, \frac{3}{2}\right), \end{aligned} $$

where in each case the first term and second terms are the contributions to the divergence-full, $\vec{M}_{in}$, and divergence-free, $\vec{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 [26]. 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.
<figure><img src="image_9.png" /><figcaption>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$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 8. (Color) A sequence of Metropolis updates in the monopole crystal phase of a dumbbell model on a square lattice. The grey ellipses show the fictive dimer positions and the blue and red circles signify north and south poles respectively. The chevrons shows the needle orientations along the bonds. The sequence of moves simulates a hard core dimer flip on a square plaquet.</figcaption></figure>

## Appendix E: Finite-Size Scaling

The susceptibility of the monopole crystallization order parameter is defined as

$$ \chi_c = \frac{\langle M_c^2 \rangle - \langle M_c \rangle^2}{N_0 T} \qquad (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 \pm 0.01$ where it becomes 1st order before disappearing for $μ^* > 0.800 \pm 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 [33] 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 1st 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 1st 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 1st 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 1st order transition (see lower panels of Figure 9). In particular, the maximum in $C_μ$ develops scaling behaviour in this region – close to $α/ν = 3$ for $L \ge 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 1st 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.
<figure><img src="image_11.png" /><figcaption>FIG. 9. (Color) Top: Finite size scaling of the maxima of the specific heat $C_μ$ and susceptibility $χ_c$ of the order parameter $M_c$ as a function of linear system size $L$, for $μ^* = 0.490$ (+), 0.654 (x), 0.768 (□), 0.778 (■), 0.784 (○), 0.794 (●), 0.795 (△), 0.796 (▲), 0.801 (※). 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”[102], 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 $(\propto L^3)$ appearing in the 1st order regime for $L \ge 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.</figcaption></figure>

[1] P. W. Anderson, “Ordering and Antiferromagnetism in Ferrites,” Phys. Rev., **102**, 1008 (1956).

[2] J. Villain, “Spin glass with non-random interactions,” Journal of Physics C: Solid State Physics, **10**, 1717 (1977).

[3] H. L. Stormer, D. C. Tsui, and A. C. Gossard, “The fractional quantum Hall effect,” Rev. Mod. Phys., **71**, S298 (1999).

[4] C.-Y. Hou, C. Chamon, and C. Mudry, “Electron Fractionalization in Two-Dimensional Graphene-like Structures,” Phys. Rev. Lett., **98**, 186809 (2007).

[5] R. Jackiw and C. Rebbi, “Solitons with fermion number,” Phys. Rev. D, **13**, 3398 (1976).

[6] H. Steinberg, G. Barak, A. Yacoby, L. N. Pfeiffer, K. W. West, B. I. Halperin, and K. Le Hur, “Charge fractionalization in quantum wires,” Nat. Phys., **4**, 116 (2008), ISSN 1745-2473.

[7] J. Bernal and R. Fowler, “A theory of water and ionic

solution, with particular reference to hydrogen and hydroxyl ions,” J. Chem. Phys., **1**, 515 (1933).

[8] P. Fulde, A. N. Yaresko, A. A. Zvyagin, and Y. Grin, “On the origin of heavy quasiparticles in $LiV_2O_4$,” Europhys. Lett., **54**, 779 (2001).

[9] S. T. Banks and S. T. Bramwell, “Magnetic frustration in the context of pseudo-dipolar ionic disorder,” Europhys. Lett., **97**, 27005 (2012).

[10] C. L. Henley, “The Coulomb Phase in Frustrated Systems,” Annual Review of Condensed Matter Physics, **1**, 179 (2010).

[11] C. Castelnovo, R. Moessner, and S. L. Sondhi, “Magnetic monopoles in spin ice,” Nature, **451**, 42 (2008).

[12] The Coulomb phase can remain ergodic, either thanks to deconfined topological excitations, or to loop like quantum fluctuations.

[13] B. C. den Hertog and M. J. P. Gingras, “Dipolar interactions and origin of spin ice in Ising pyrochlore magnets,”
Phys. Rev. Lett., **84**, 3430 (2000).

[14] I. A. Ryzhkin, “Magnetic relaxation in rare-earth oxide pyrochlores,” Journal of Experimental and Theoretical Physics, **101**, 481 (2005).

[15] S. V. Isakov, R. Moessner, and S. L. Sondhi, “Why spin ice obeys the ice rules,” Phys. Rev. Lett., **95**, 217201 (2005).

[16] S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi, “Dipolar spin correlations in classical pyrochlore magnets,” Phys. Rev. Lett., **93**, 167204 (2004).

[17] T. Fennell, P. P. Deen, A. R. Wildes, K. Schmalzl, D. Prabhakaran, A. T. Boothroyd, R. J. Aldus, D. F. McMorrow, and S. T. Bramwell, “Magnetic Coulomb Phase in the Spin Ice $Ho_2Ti_2O_7$,” Science, **326**, 415 (2009).

[18] G. Sala, C. Castelnovo, R. Moessner, S. L. Sondhi, K. Kitagawa, M. Takigawa, R. Higashinaka, and Y. Maeno, “Magnetic Coulomb Fields of Monopoles in Spin Ice and Their Signatures in the Internal Field Distribution,” Phys. Rev. Lett., **108**, 217203 (2012).

[19] S. R. Dunsiger, A. A. Aczel, C. Arguello, H. Dabkowska, A. Dabkowski, M. H. Du, T. Goko, B. Javanparast, T. Lin, F. L. Ning, H. M. L. Noad, D. J. Singh, T. J. Williams, Y. J. Uemura, M. J. P. Gingras, and G. M. Luke, “Spin Ice: Magnetic Excitations without Monopole Signatures Using Muon Spin Rotation,” Phys. Rev. Lett., **107**, 207207 (2011).

[20] A. C. Maggs and V. Rossetto, “Local Simulation Algorithms for Coulomb Interactions,” Phys. Rev. Lett., **88**, 196402 (2002).

[21] S. T. Bramwell, “Generalized longitudinal susceptibility for magnetic monopoles in spin ice,” Phil. Trans. R. Soc. A, **370**, 5738 (2012).

[22] R. Moessner and S. L. Sondhi, “Irrational Charge from Topological Order,” Phys. Rev. Lett., **105**, 166401 (2010).

[23] A. Sen, R. Moessner, and S. L. Sondhi, “Coulomb Phase Diagnostics as a Function of Temperature, Interaction Range, and Disorder,” Phys. Rev. Lett., **110**, 107202 (2013).

[24] L. Savary and L. Balents, “Coulombic Quantum Liquids in Spin-1/2 Pyrochlores,” Phys. Rev. Lett., **108**, 37202 (2012).

[25] D. A. Huse, W. Krauth, R. Moessner, and S. L. Sondhi, “Coulomb and Liquid Dimer Models in Three Dimensions,” Phys. Rev. Lett., **91**, 167004 (2003).

[26] F. Alet, J. L. Jacobsen, G. Misguich, V. Pasquier, F. Mila, and M. Troyer, “Interacting Classical Dimers on the Square Lattice,” Phys. Rev. Lett., **94**, 235702 (2005).

[27] D. L. Bergman, G. A. Fiete, and L. Balents, “Ordering in a frustrated pyrochlore antiferromagnet proximate to a spin liquid,” Phys. Rev. B, **73**, 134402 (2006).

[28] O. Sikora, F. Pollmann, N. Shannon, K. Penc, and P. Fulde, “Quantum Liquid with Deconfined Fractional Excitations in Three Dimensions,” Phys. Rev. Lett., **103**, 247001 (2009).

[29] R. A. Borzi, D. Slobinsky, and S. A. Grigera, “Charge Ordering in a Pure Spin Model: Dipolar Spin Ice,” Phys. Rev. Lett., **111**, 147204 (2013).

[30] J. N. Reimers, J. E. Greedan, and M. Björgvinsson, “Critical properties of highly frustrated pyrochlore antiferromagnets,” Phys. Rev. B, **45**, 7295 (1992).

[31] L. D. C. Jaubert and P. C. W. Holdsworth, “Magnetic Monopole Dynamics in Spin Ice,” Journal of Physics-Condensed Matter (2011).

[32] R. G. Melko and M. J. P. Gingras, “Monte Carlo studies of the dipolar spin ice model,” Journal of Physics-Condensed Matter, **16**, R1277 (2004).

[33] V. Kobelev, B. Kolomeisky, and M. E. Fisher, “Lattice models of ionic systems,” J. Chem. Phys., **116**, 7589 (1992).

[34] R. Dickman and G. Stell, “Phase diagram of the lattice restricted primitive model,” AIP Conf. Proc., **492**, 225 (1999).

[35] A. Harman-Clarke, *Topological constraints and ordering in model frustrated magnets*, Ph.D. thesis, University College London and Ecole Normale Supérieure de Lyon (2011).

[36] G. Möller and R. Moessner, “Magnetic multipole analysis of kagome and artificial spin-ice dipolar arrays,” Phys. Rev. B, **80**, 140409 (2009).

[37] G.-W. Chern, P. Mellado, and O. Tchernyshyov, “Two-Stage Ordering of Spins in Dipolar Spin Ice on the Kagome Lattice,” Phys. Rev. Lett., **106**, 207202 (2011).

[38] G.-W. Chern and O. Tchernyshyov, “Magnetic charge and ordering in kagome spin ice,” Phil. Trans. R. Soc. A, **370**, 5718 (2012).

[39] A. J. Macdonald, P. C. W. Holdsworth, and R. G. Melko, “Classical topological order in kagome ice,” J. Phys: Cond. Mat., **23**, 164208 (2011).

[40] O. Cépas and B. Canals, “Heterogeneous freezing in a geometrically frustrated spin model without disorder: Spontaneous generation of two time scales,” Phys. Rev. B, **86**, 024434 (2012).

[41] K. Matsuhira, Z. Hiroi, T. Tayama, S. Takagi, and T. Sakakibara, “A new macroscopically degenerate ground state in the spin ice compound $Dy_2Ti_2O_7$ under a magnetic field,” J. Phys: Cond. Mat., **14**, L559 (2002).

[42] T. Fennell, S. T. Bramwell, D. F. McMorrow, P. Manuel, and A. R. Wildes, “Pinch points and Kasteleyn transitions in kagome ice,” Nat. Phys., **3**, 566 (2007).

[43] P. Fulde, K. Penc, and N. Shannon, “Fractional charges in pyrochlore lattices,” Ann. Phys. (Leipzig), **11**, 892 (2002).

[44] L. D. C. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R. Moessner, “Three-Dimensional Kasteleyn Transition: Spin Ice in a [100] Field,” Phys. Rev. Lett., **100**, 067207 (2008).

[45] R. G. Melko, B. C. den Hertog, and M. J. P. Gingras, “Long-Range Order at Low Temperatures in Dipolar Spin Ice,” Phys. Rev. Lett., **87**, 067203 (2001).

[46] L. D. C. Jaubert and P. C. W. Holdsworth, “Signature of magnetic monopole and Dirac string dynamics in spin ice,” Nat. Phys., **5**, 258 (2009).

[47] D. S. Rokhsar and S. A. Kivelson, “Superconductivity and the Quantum Hard-Core Dimer Gas,” Phys. Rev. Lett., **61**, 2376 (1988).

[48] S. V. Isakov, K. S. Raman, R. Moessner, and S. L. Sondhi, “Magnetization curve of spin ice in a [111] magnetic field,” Phys. Rev. B, **70**, 104418 (2004).

[49] L. D. C. Jaubert, M. J. Harris, T. Fennell, R. G. Melko, S. T. Bramwell, and P. C. W. Holdsworth, “Topological-Sector Fluctuations and Curie-Law Crossover in Spin Ice,” Phys. Rev. X, **3**, 011014 (2013).

[50] S. R. Giblin, S. T. Bramwell, P. C. W. Holdsworth, Monopole Dynamics in Spin Ice,” Journal of Physics-Condensed Matter (2011).
D. Prabhakaran, and I. Terry, "Creation and Mea-
surement of Long-Lived Magnetic Monopole Currents
in Spin Ice," Nat. Phys., **7**, 252 (2011).

[51] T. Sakakibara, T. Tayama, Z. Hiroi, K. Matsuhira, and S. Takagi, "Observation of a Liquid-Gas-Type Transition in the Pyrochlore Spin Ice Compound Dy₂Ti₂O₇ in a Magnetic Field," Phys. Rev. Lett., **90**, 207205 (2003).

[52] R. Moessner and S. L. Sondhi, "Theory of the [111] magnetization plateau in spin ice," Phys. Rev. B, **68**, 064411 (2003).

[53] A. P. Sazonov, A. Gukasov, I. Mirebeau, and P. Bonville, "Double-layered monopolar order in the Tb₂Ti₂O₇ spin liquid," Phys. Rev. B, **85**, 214420 (2012).

[54] H. D. Zhou, S. T. Bramwell, J. G. Cheng, C. R. Wiebe, G. Li, L. Balicas, J. A. Bloxsom, H. J. Silverstein, J. S. Zhou, J. B. Goodenough, and J. S. Gardner, "High pressure route to generate magnetic monopole dimers in spin ice," Nat Commun, **2**, 478 (2011).

[55] J. Gardner, (2013), private communication.

[56] K. A. Ross, L. Savary, B. D. Gaulin, and L. Balents, "Quantum Excitations in Quantum Spin Ice," Phys. Rev. X, **1**, 021002 (2011).

[57] L. Bovo, J. Bloxsom, D. Prabhakaran, G. Aeppli, and S. Bramwell, "Brownian motion and quantum dynamics of magnetic monopoles in spin ice," Nat. Commun., **4**, 1535 (2013).

[58] J. S. Gardner, S. R. Dunsiger, B. D. Gaulin, M. J. P. Gingras, J. E. Greedan, R. F. Kiefl, M. D. Lumsden, W. A. MacFarlane, N. P. Raju, J. E. Sonier, I. Swainson, and Z. Tun, "Cooperative paramagnetism in the geometrically frustrated pyrochlore antiferromagnet Tb₂Ti₂O₇," Phys. Rev. Lett., **82**, 1012 (1999).

[59] H. R. Molavian, M. J. P. Gingras, and B. Canals, "Dynamically induced frustration as a route to a quantum spin ice state in Tb₂Ti₂O₇ via virtual crystal field excitations and quantum many-body effects," Phys. Rev. Lett., **98**, 157204 (2007).

[60] M. J. P. Gingras, B. C. den Hertog, M. Faucher, J. S. Gardner, S. R. Dunsiger, L. J. Chang, B. D. Gaulin, N. P. Raju, and J. E. Greedan, "Thermodynamic and single-ion properties of Tb³⁺ within the collective paramagnetic-spin liquid state of the frustrated pyrochlore antiferromagnet Tb₂Ti₂O₇," Phys. Rev. B, **62**, 6496 (2000).

[61] J. S. Gardner, A. Keren, G. Ehlers, C. Stock, E. Segal, J. M. Roper, B. Fåk, M. B. Stone, P. R. Hammar, D. H. Reich, and B. D. Gaulin, "Dynamic frustrated magnetism in Tb₂Ti₂O₇ at 50 mK," Phys. Rev. B, **68**, 180401 (2003).

[62] T. Fennell, M. Kenzelmann, B. Roessli, M. K. Haas, and R. J. Cava, "Power-Law Spin Correlations in the Pyrochlore Antiferromagnet Tb₂Ti₂O₇," Phys. Rev. Lett., **109**, 017201 (2012).

[63] S. Petit, P. Bonville, J. Robert, C. Decorse, and I. Mirebeau, "Spin liquid correlations, anisotropic exchange, and symmetry breaking in Tb₂Ti₂O₇," Phys. Rev. B, **86**, 174403 (2012).

[64] I. Mirebeau, I. N. Goncharenko, P. Cadavez-Peres, S. T. Bramwell, M. J. P. Gingras, and J. S. Gardner, "Pressure-induced crystallization of a spin liquid," Nature, **420**, 54 (2002).

[65] S. Legl, C. Krey, S. R. Dunsiger, H. A. Dabkowska, J. A. Rodriguez, G. M. Luke, and C. Pfleiderer, "Vibrating-Coil Magnetometry of the Spin Liquid Properties of

Tb$_{2}$Ti$_{2}$O$_{7}$,” Phys. Rev. Lett., **109**, 047201 (2012).

[66] J. P. C. Ruff, B. D. Gaulin, K. C. Rule, and J. S. Gardner, "Superlattice correlations in Tb₂Ti₂O₇ under the application of [110] magnetic field," Phys. Rev. B, **82** (2010), doi:DOI 10.1103/PhysRevB.82.100401.

[67] T. Taniguchi, H. Kadowaki, H. Takatsu, B. Fåk, J. Ollivier, T. Yamazaki, T. J. Sato, H. Yoshizawa, Y. Shimura, T. Sakakibara, T. Hong, K. Goto, L. R. Yaraskavitch, and J. B. Kycia, "Long-range order and spin-liquid states of polycrystalline Tb$_{2+x}$Ti$_{2-x}$O$_{7+y}$," Phys. Rev. B, **87**, 060408 (2013).

[68] J. P. C. Ruff, Z. Islam, J. P. Clancy, K. A. Ross, H. Nojiri, Y. H. Matsuda, H. A. Dabkowska, A. D. Dabkowski, and B. D. Gaulin, "Magnetoplastics of a Spin Liquid: X-Ray Diffraction Studies of Tb$_{2}$Ti$_{2}$O$_{7}$ in Pulsed Magnetic Fields," Phys. Rev. Lett., **105**, 077203 (2010).

[69] S. Guitteny, J. Robert, P. Bonville, J. Ollivier, C. Decorse, P. Steffens, M. Boehm, H. Mutka, I. Mirebeau, and S. Petit, "Anisotropic Propagating Excitations and Quadrupolar Effects in Tb$_{2}$Ti$_{2}$O$_{7}$," Phys. Rev. Lett., **111**, 087201 (2013).

[70] P. Bonville, I. Mirebeau, A. Gukasov, S. Petit, and J. Robert, "Tetragonal distortion yielding a two-singlet spin liquid in pyrochlore Tb$_{2}$Ti$_{2}$O$_{7}$," Phys. Rev. B, **84**, 184409 (2011).

[71] B. D. Gaulin, J. S. Gardner, P. A. McClarty, and M. J. P. Gingras, "Lack of evidence for a singlet crystal-field ground state in the magnetic pyrochlore Tb$_{2}$Ti$_{2}$O$_{7}$," Phys. Rev. B, **84**, 140402 (2011).

[72] K. Fritsch, K. A. Ross, Y. Qiu, J. R. D. Copley, T. Guidi, R. I. Bewley, H. A. Dabkowska, and B. D. Gaulin, "Antiferromagnetic spin ice correlations at (1/2, 1/2, 1/2) in the ground state of the pyrochlore magnet Tb$_{2}$Ti$_{2}$O$_{7}$," Phys. Rev. B, **87**, 094410 (2013).

[73] M. Enjalran, M. J. P. Gingras, Y. J. Kao, A. Del Maestro, and H. R. Molavian, "The spin liquid state of the Tb$_{2}$Ti$_{2}$O$_{7}$ pyrochlore antiferromagnet: a puzzling state of affairs," Journal of Physics-Condensed Matter, **16**, S673 (2004).

[74] I. Mirebeau, I. N. Goncharenko, G. Dhalenne, and A. Revcolevschi, "Pressure and Field Induced Magnetic Order in the Spin Liquid Tb$_{2}$Ti$_{2}$O$_{7}$ as Studied by Single Crystal Neutron Diffraction," Phys. Rev. Lett., **93**, 187204 (2004).

[75] T. Fennell, O. A. Petrenko, G. Balakrishnan, S. T. Bramwell, J. D. M. Champion, B. Fak, M. J. Harris, and D. M. Paul, "Field-induced partial order in the spin ice dysprosium titanate," Applied Physics A-Materials Science & Processing, **74**, S889 (2002).

[76] Z. Hiroi, K. Matsuhira, and M. Ogata, "Ferromagnetic Ising spin chains emerging from the spin ice under magnetic field," Journal of the Physical Society of Japan, **72**, 3045 (2003).

[77] T. Fennell, M. Kenzelmann, B. Roessli, H. Mutka, J. Ollivier, M. Ruminy, U. Stuhr, O. Zaharko, L. Bovo, A. Cervellino, M. Haas, and R. Cava, "Magnetoplastics excitations in the pyrochlore spin liquid Tb$_{2}$Ti$_{2}$O$_{7}$," arXiv:1305.5405 (2013).

[78] Y. Yasui, M. Kanada, M. Ito, H. Harashina, M. Sato, H. Okumura, K. Kakurai, and H. Kadowaki, "Static correlation and dynamical properties of Tb$^{3+}$-moments in Tb$_{2}$Ti$_{2}$O$_{7}$ - Neutron scattering study," Journal of the Physical Society of Japan, **71**, 599 (2002).

[79] E. Lhotel, C. Paulsen, P. D. de Réotier, A. Yaouanc,
C. Marin, and S. Vanishri, “Low-temperature magnetization in geometrically frustrated Tb$_2$Ti$_2$O$_7$,” Phys. Rev. B, **86**, 020410 (2012).

[80] N. R. Hayre, K. A. Ross, R. Applegate, T. Lin, R. R. P. Singh, B. D. Gaulin, and M. J. P. Gingras, “Thermodynamic properties of Yb$_2$Ti$_2$O$_7$ pyrochlore as a function of temperature and magnetic field: Validation of a quantum spin ice exchange Hamiltonian,” Phys. Rev. B, **87**, 184423 (2013).

[81] S. T. Bramwell, M. N. Field, M. J. Harris, and I. P. Parkin, “Bulk magnetization of the heavy rare earth titanate pyrochlores - a series of model frustrated magnets,” J. Phys: Cond. Mat., **12**, 483 (2000).

[82] J. A. Hodges, P. Bonville, A. Forget, A. Yaouanc, P. Dalmas de Réotier, G. André, M. Rams, K. Królas, C. Ritter, P. C. M. Gubbens, C. T. Kaiser, P. J. C. King, and C. Baines, “First-Order Transition in the Spin Dynamics of Geometrically Frustrated Yb$_2$Ti$_2$O$_7$,” Phys. Rev. Lett., **88**, 077204 (2002).

[83] J. S. Gardner, G. Ehlers, N. Rosov, R. W. Erwin, and C. Petrovic, “Spin-spin correlations in Yb$_2$Ti$_2$O$_7$: A polarized neutron scattering study,” Phys. Rev. B, **70**, 180404 (2004).

[84] K. A. Ross, L. R. Yaraskavitch, M. Laver, J. S. Gardner, J. A. Quilliam, S. Meng, J. B. Kycia, D. K. Singh, T. Proffen, H. A. Dabkowska, and B. D. Gaulin, “Dimensional evolution of spin correlations in the magnetic pyrochlore Yb$_2$Ti$_2$O$_7$,” Phys. Rev. B, **84**, 174442 (2011).

[85] L.-J. Chang, S. Onoda, Y. Su, Y.-J. Kao, K.-D. Tsuei, Y. Yasui, K. Kakurai, and M. R. Lees, “Higgs transition from a magnetic Coulomb liquid to a ferromagnet in Yb$_2$Ti$_2$O$_7$,” Nat. Commun., **3**, 992 (2012).

[86] L. Savary and L. Balents, “Spin liquid regimes at nonzero temperature in quantum spin ice,” Phys. Rev. B, **87**, 205130 (2013).

[87] K. A. Ross, T. Proffen, H. A. Dabkowska, J. A. Quilliam, L. R. Yaraskavitch, J. B. Kycia, and B. D. Gaulin, “Lightly stuffed pyrochlore structure of single-crystalline Yb$_2$Ti$_2$O$_7$ grown by the optical floating zone technique,” Phys. Rev. B, **86**, 174424 (2012).

[88] R. M. D'Ortenzio, H. A. Dabkowska, S. R. Dunsiger, B. D. Gaulin, M. J. P. Gingras, T. Goko, J. B. Kycia, L. Liu, T. Medina, T. J. Munsie, D. Pomaranksi, K. A. Ross, Y. J. Uemura, T. J. Williams, and G. M. Luke, “Unconventional Magnetic Ground State in Yb$_2$Ti$_2$O$_7$,” ArXiv e-prints (2013), arXiv:1303.3850 [cond-mat.str-el].

[89] Y. Han, L. D. C. Jaubert, O. Benton, and N. Shannon, in preparation (2013).

[90] I. Mirebeau, A. Apetrei, J. Rodriguez-Carvajal, P. Bonville, A. Forget, D. Colson, V. Glazkov, J. P. Sanchez, O. Isnard, and E. Suard, “Ordered Spin Ice State and Magnetic Fluctuations in Tb$_2$Sn$_2$O$_7$,” Phys. Rev. Lett., **94**, 246402 (2005).

[91] M. L. Dahlberg, M. J. Matthews, P. Jiramongkolchai, R. J. Cava, and P. Schiffer, “Low-temperature dynamic freezing and the fragility of ordering in Tb$_2$Sn$_2$O$_7$,” Phys. Rev. B, **83**, 140410 (2011).

[92] P. M. McClarty, P. Stasiak, and M. J. P. Gingras, “Soft dipolar spin ice physics and the ordered phase of the frustrated Tb$_2$Sn$_2$O$_7$ pyrochlore magnet,” arXiv:1011.6346v1 (2010).

[93] S. Zhang, I. Gilbert, C. Nisoli, G.-W. Chern, M. J.

Erickson, L. O'Brien, C. Leighton, P. E. Lammert, V. H. Crespi, and P. Schiffer, “Crystallites of magnetic charges in artificial spin ice,” Nature, **500**, 553 (2013).

[94] O. A. Petrenko, C. Ritter, M. Yethiraj, and D. McK Paul, “Investigation of the Low-Temperature Spin-Liquid Behavior of the Frustrated Magnet Gadolinium Gallium Garnet,” Phys. Rev. Lett., **80**, 4570 (1998).

[95] S.-H. Lee, C. Broholm, T. H. Kim, W. Ratcliff, and S.-W. Cheong, “Local Spin Resonance and Spin-Peierls-like Phase Transition in a Geometrically Frustrated Antiferromagnet,” Phys. Rev. Lett., **84**, 3718 (2000).

[96] J. S. Gardner, M. J. P. Gingras, and J. E. Greedan, “Magnetic pyrochlore oxides,” Reviews of Modern Physics, **82**, 53 (2010).

[97] H. B. Cao, A. Gukasov, I. Mirebeau, and P. Bonville, “Anisotropic exchange in frustrated pyrochlore Yb$_2$Ti$_2$O$_7$,” J. Phys: Cond. Mat., **21**, 492202 (2009).

[98] M. Hermele, M. P. A. Fisher, and L. Balents, “Pyrochlore photons: The U(1) spin liquid in a S=1/2 three-dimensional frustrated magnet,” Phys. Rev. B, **69**, 64404 (2004).

[99] O. Benton, O. Sikora, and N. Shannon, “Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice,” Phys. Rev. B, **86**, 075154 (2012).

[100] R. Applegate, N. R. Hayre, R. R. P. Singh, T. Lin, A. G. R. Day, and M. J. P. Gingras, “Vindication of Yb$_2$Ti$_2$O$_7$ as a Model Exchange Quantum Spin Ice,” Phys. Rev. Lett., **109**, 097205 (2012).

[101] R. H. Swendsen and J.-S. Wang, “Replica Monte Carlo Simulation of Spin-Glasses,” Phys. Rev. Lett., **57**, 2607 (1986).

[102] I. Wolfram Research, Mathematica Edition: Version 9.0 (Wolfram Research, Inc., Champaign, Illinois, 2012).