Spin slush in an extended spin ice model

We present a new classical spin liquid on the pyrochlore lattice by extending spin ice with further neighbour interactions. We find that this disorder-free spin model exhibits a form of dynamical heterogeneity with extremely slow relaxation for some spins, while others fluctuate quickly down to zero temperature. We thus call this state spin slush, in analogy to the heterogeneous mixture of solid and liquid water. This behaviour is driven by the structure of the ground-state manifold which extends the celebrated two-in/two-out ice states to include branching structures built from three-in/one-out, three-out/one-in and all-in/all-out tetrahedra defects. Distinctive liquid-like patterns in the magnetic correlations serve as a signature of this intermediate range order. Possible applications to materials as well the effects of quantum tunnelling are discussed.

T he physics of glasses plays an important role in many types of physical systems; from its origins in the physics of liquids 1 further realizations have been found in disordered magnets 2 , superconductors 3 and metals 4 through to softcondensed matter 5 and even biophysics 6 . While ubiquitous, a complete understanding of glasses remains an important open problem in condensed matter physics. Connections between these vastly different contexts have proven fruitful in making progress. For example, studying conceptually and computationally simpler spin models, may inform the physics of super-cooled liquids and structural glasses 1 . However, there are complications-while spin glasses are driven by the combination of random quenched disorder and frustration 2 , glass-forming liquids are intrinsically disorder-free 1 . Finding a disorder-free spin model that realizes the diverse range of phenomena observed in glass formers, such as the dramatic slowing down of relaxation and emergence of spatially heterogeneous dynamics, is a serious challenge. Some examples of disorder-free spin models with strong freezing have been proposed [7][8][9][10][11][12] . Each of these proposals has some limitations, be it the lack of heterogeneous dynamics, the need for multi-spin interactions, the use of uncontrolled approximations or the introduction of non-local dynamics.
In this article, we introduce a cooperative paramagnet that we call spin slush (SS) which appears in an extended spin ice (ESI) model. This classical SS model is disorder-free and includes only first-, second-and third-neighbour Ising bilinear exchange interactions and thus lacks the shortcomings discussed above. We start from spin ice (SI), a well-studied magnetic analogue of common water ice 13 , magnetic moments pointing in or out of the corner-shared tetrahedra of the pyrochlore lattice embody the proton displacements of water ice 14 . Similar to water ice, SI displays an extensive ground-state degeneracy, and thus an associated extensive residual entropy, characterized by the two-in/two-out ice rule condition on each tetrahedron 13 . In SS, we find that the ground-state manifold of SS is larger than that of SI and contains a far richer set of states. In addition to the two-in/two-out tetrahedra of the SI ground-state manifold, there are spatially extended structures assembled from three-in/one-out, three-out/one-in and all-in/all-out tetrahedra. Built from SI defects, these structures are not simply loops or strings, but include branching tree-like objects. After characterizing the static thermodynamic and magnetic properties of SS, we turn to dynamics. Approaching zero temperature, we find freezing, as in SI 15,16 , with an exponentially increasing average relaxation time. However, unlike in SI where all of the spins freeze uniformly as the temperature is lowered, the spins in the SS exhibit highly heterogeneous dynamics reminiscent of glass formers 17 . While many of the spins strongly freeze with an extremely slow relaxation rate, a fraction of the spins, organized into spatially local clusters, remain completely dynamic, relaxing almost immediately. Since this model is disorder-free, the random distribution of these dynamical spins derives solely from the overall freezing behaviour. This dynamical heterogeneity in SS at low temperatures motivates the name SS, in analogy to slush where liquid water and solid ice coexist as a mixture. Finally, we speculate on the behaviour of quantum SS as well as possible experimental relevance in frustrated pyrochlore magnets.

Results
Model. We start with a review of some key results for the nearestneighbour SI model 18 to establish our notation and motivate the SS model. The SI model is a nearest-neighbour Ising antiferromagnet on the pyrochlore lattice, with Hamiltonian J P hiji s i s j , where s i ¼ ±1 are the Ising spins. This can be reformulated in terms of ice rule defects, or charges, defined on each tetrahedron. With each tetrahedron identified with a dual diamond lattice site I, one defines the charge Q I 1 2 ð À 1Þ I P i2I s i , where ( À 1) I is a sign reflecting the sublattice of I. In this language, the nearest-neighbour SI Hamiltonian simply penalizes non-zero charges, taking the form The ground states of this model are those with Q I ¼ 0 for all tetrahedra, i.e., the celebrated two-in/two-out states of the ice manifold. This manifold is macroscopically degenerate with a residual entropy given approximately by S SI $ ðNk B =2Þlogð3=2Þ $ 0:202Nk B (ref. 13). Because of this extensive ground-state degeneracy, addition of small perturbations will generically select an ordered state from this manifold at low temperatures 13 .
To explore the effects of such perturbations, we consider the addition of second-and third-neighbour Ising exchanges of the form For third-neighbour exchange, we include only those that are composed of two nearest-neighbour steps, as illustrated in Fig. 1 and Supplementary Fig. 1. For many mechanisms that generate such interactions in real materials, for example super-exchange or third-neighbour exchange can be present in real materials either intrinsically 19,20 , or via a partial cancellation of the leading terms 21 . One can show that for any SI state, X We thus see that two terms are not independent and when J 2 ¼ J 3a J 0 , they effectively cancel each other when in the SI manifold. Moving along the J 2 ¼ J 3a line, the model moves away from the nearest-neighbour SI regime, but without lifting the degeneracy of the SI manifold. While SI persists as the ground state at low temperature for sufficiently small J 0 /J, eventually it gives way when another set of states crosses the SI manifold. This line of degeneracy also exists for analogous models with N-components spins, although how the termination manifests depends on the precise value of N. We thus refer to the model along this line as ESI. It will prove useful to write this model in terms of the charges Q I as We see that J 0 40 generates an attraction between nearestneighbour charges of the same sign. This short-range attraction between charges will play a central role in understanding the ground and excited states of ESI . One can show (see Supplementary Note 1) that the SI manifold persists until J 0 ¼ J/4 for J 0 40 and to J 0 ¼ À J/2 for J 0 o0. Going along the J 2 ¼ J 3a line towards negative values (away from the SS) in our model, the end point has a similar manifold of states (distinct from the SS states) to that studied in ref. 22. Specifically, at J 0 ¼ À J/2, the ground states include all configurations with staggered charge Q I ¼ Q 0 ( À 1) I . These are the ice states (Q 0 ¼ 0), the single-charge states (Q 0 ¼ 1)) and the all-in, all-out states (Q 0 ¼ 2). The collapse of excited states when approaching J 0 ¼ J/4 is illustrated in Fig. 2. We show only the simplest examples that cross the ice manifold, but as we shall see, there are an infinite set of such states. We focus on the end point at J 0 ¼ J/4 which we refer to as the SS model. At this special point, one can write the model (Eq. 4) as In this form, one notes a strong similarity to the SI model of equations (2), except with the fundamental unit now being a pair of tetrahedra, indicated by B þ B þ , rather than a single tetrahedron.
Ground-state manifold. The ground-state manifold of SS is most easily characterized in terms of the variables Following equations (7), any state with P i ¼ ±½ for all sites has the minimal energy À 3NJ/4 and is in the ground-state manifold. Alternatively, we can write this in terms of the SI charges; associating each site i of the pyrochlore lattice with a nearest-neighbour bond hIJi of the dual diamond lattice, one has P i ¼ ( À 1) I (Q I À Q J ) À s i /2. From this expression for P i in terms of the SI charges, it is clear that any SI state with Q I ¼ 0 for all sites also belongs to the SS manifold. In addition to the familiar SI states, many more states satisfy P i ¼ ±½. A naive enumeration of states for an isolated pair of tetrahedra shows that beyond the 18 ice states, there are an additional 52 states, 70 in total, that belong to the SS manifold. We caution that a Pauling-like estimate severely underestimates the degeneracy of the SS manifold. Given that the number of tetrahedron pairs is equal to the number of sites, one would estimate a residual entropy of Nk B (log 2 þ log (70/2 7 ))B0.0896Nk B . This reflects that the constraints provided by P i ¼ ± ½ are much less independent than in SI where Pauling's estimate is accurate. These additional states include configurations with both single-charge (Q I ¼ ± 1) and doublecharge (Q I ¼ ±2) defects. The influence of the nearest-neighbour attraction of charges manifests here; pairs of like single charges can appear together, while double charges only appear with accompanying single charges of the same sign. One finds from equation (5) that the energy cost of having a charge can be compensated by the energy gain of having two neighbouring charges of the same sign. From these observations, we formulate rules for constructing states that satisfy P i ¼ ±½. We formulate these rules from the perspective of specifying non-ice tetrahedra (Q I a0) states first, then populating the remaining tetrahedra with any compatible ice states afterward. The first rule for placing the non-ice, charged tetrahedra, is the single-charge rule: this states that the minority spin of a single charge, Q I ¼ ± 1, must be connected to a tetrahedron carrying a single or double charge of the same sign. The second rule, the double-charge rule, states that a double charge Q I ¼ ±2 must have its four nearest-neighbour tetrahedra occupied by single charges of the same sign. Finally the neighbour rule requires that a single charge, Q I ¼ ±1, cannot have any single charges of opposite sign as nearest neighbours. Once single and double charges have been placed such that they satisfy the above three rules, one can assign the remaining tetrahedra any allowed ice rule, Q I ¼ 0, states. The first rule allows the singlecharge tetrahedra to form branching tree-like structures where the minority spin of a given charge also belongs to the next charge in the structure. Each branch must terminate in some way compatible with the rules, so the minority spin must end up on another single charge. The possibilities for terminating a branch include looping back to itself, ending on a different branch or on one of the single charges attached to a double charge. Note that these single-and double-charge structures must exist for both signs of the charge to satisfy the global neutrality requirement P I Q I ¼ 0. The third rule implies that charge structures of opposite sign must be separated by at least one ice rule obeying tetrahedron. An illustration of an SS state incorporating all of these features, restricted to a single [111] kagomé plane, is shown in Fig. 1.
Thermodynamic and magnetic properties. With the ground states of SS identified, we now explore its finite temperature properties via classical Monte Carlo simulations using singlespin-flip dynamics, augmented with parallel tempering when appropriate. Basic thermodynamic quantities are shown in Fig. 3. The specific heat exhibits a broad peak at T*B0.3J, reminiscent of the peak seen in SI. This peak signals the release of entropy as one begins to enter the SS ground-state manifold. This can be seen explicitly in the entropy in Fig. 3 where, below T*, the entropy approaches the constant value S SS B0.266Nk B . As expected from the rules derived in the previous section, this is significantly higher than S SI B0.202Nk B found in SI. At these low temperatures severe freezing is encountered, preventing the simulations from reaching equilibrium below TB0.15J. It is not obvious how to construct a non-local move that would sample the SS manifold efficiently. Including the SI loop move does aid equilibriation, but it is only effective in regions where no single-and double-charge defects are present, leaving the freezing problem for future work. The frozen states belong to the SS manifold and exhibit the single-and double-charge structures discussed in the previous section. We found no evidence of ordering in any of our simulations, be it conventional or via an order-by-disorder mechanism. Further, the specific heat and entropy are somewhat immune to this freezing problem, showing consistent behaviour between simulations. The magnetic properties are however more sensitive.
The simplest probe of the magnetic behaviour is the uniform susceptibility, w, shown in Fig. 4, for the moments m i s iẑi pointing in/out of the tetrahedra along the local [111] directionẑ i . At both low and high temperatures, one finds Curie-like behaviour, with 3wT constant, separated by a broad peak at TBO(J). The constant approached as T-0 depends on the details of how the system freezes. This varies between simulations, taking on a distribution of values clustered around 3wTB1, reflected in the large error bars in Fig. 4. A more detailed probe of the magnetic structure are the spin-spin correlations, as can be investigated via neutron scattering. Recall that in SI the appearance of sharp pinch-points 23 in the transverse momentmoment correlation function signals the development of long-range dipolar spin-spin correlations. In SS, one finds sharp features in I(k) distinct from such pinch points. As shown in Fig. 4, below T* I(k) develops into sharp rings centred on zone centres in a given plane of reciprocal space. In the full [hkl] space, these features lie approximately on spheres, reminiscent of an isotropic liquid. As discussed in Supplementary Note 2, this analogy is even more striking in the structure factor of the SI charges Q I where the intensity is approximately uniform across the sphere, as shown in Supplementary Figs 2 and 3. The wave-vector jkj $ 0:5ð2p=aÞ k Ã , where a is the size of the conventional cubic unit cell, indicates these correlations have a characteristic length of two cubic cells and thus represent intermediate scale correlations. These correlations are consistent with the typical size of the charged structures that appear in the ESI manifold. Indeed, as seen in Fig. 1, even the smallest of these structures can span several cubic unit cells. These simulations confirm that the SS model does not order and the SS manifold shows all the rich charge structures at intermediate length scales implied by the SS rules. Indeed, at low temperatures, a significant fraction of tetrahedra, B30-35%, carry single charges, while a smaller but finite fraction, a per cent or so, carry double charges. Similar to the susceptibility, the amount of single and double charges present at low temperatures  varies somewhat from run to run, a consequence of the severe freezing problem. To better understand this issue, we now look more closely at the low-temperature dynamics of SS.
Dynamics and spin slush. To reflect the physics of real systems with local dynamics, we employ only single-spin flip, Metropolis dynamics, although we expect any local dynamics to give qualitatively the same behaviour. We primarily consider the site-resolved auto-correlation functions, defining is the Ising spin at a given Monte Carlo sweep t at site i, averaging over many initial times t 0 . Generically, one would expect exponential relaxation A i ðtÞ $ e À t=ti with a characteristic relaxation time t i . Indeed this is found in SI, with the relaxation time being siteindependent, with t i Bt and increasing exponentially as temperature is lowered 16 .
In contrast to SI, the dynamics in SS vary strongly from site to site. As temperature is lowered, most of the sites freeze, with their relaxation times becoming very long, similar to what is found in SI 15,16 . This can be seen in the site-averaged auto-correlation function AðtÞ shown in Fig. 5. However, there are clear differences, namely in the initial decrease and plateau in AðtÞ at short times as well as in the larger site to site variance in A i (t) at low temperatures. We can understand this behaviour by looking at the T-0 limit; one finds that a fraction of sites remain highly dynamic down to very low temperatures. This is illustrated in Fig. 5, where the site-resolved auto-correlation functions are shown for T ¼ 10 À 4 J. The frozen spins have A i (t) ¼ 1 at all times, while the unfrozen spins have A i (t) relaxing in 10 1 -10 2 Monte Carlo sweeps to a constant value lim t-N A i (t) A i (N)o1. To be precise, for the long-time limit lim t-N A i (t), we mean tc1 but still much smaller than the slow timescale BO(e J/T ). A non-zero value of A i (N)o1 indicates that, while fluctuating, on average more time is spent in one of the states s i ¼ ± 1 than the other. For example, if s i is sampling uniformly from values s (1) ,y, s (m) as a function of time, then A i ð1Þ $ ð 1 m P m n¼1 s ðnÞ Þ 2 at long times. Figure 5 shows that the long-time values A i (N) cluster about the squares of rational numbers, as would be expected from the above discussion. In these annealed simulations, the frozen spins make up the bulk of the system, while the number of unfrozen, dynamic spins is on the order of a few per cent.
To better understand these dynamic spins, we examine their real space structure. We find that these spins are spatially correlated, forming clusters of varying size n c . A dynamical cluster is defined by a set of spins, where lim t-N A i (t)o1 and each spin is connected by a first, second or third neighbour bond to another spin in the cluster. The SS state at low temperature is thus a mixture where regions of frozen and unfrozen spins coexist. Dynamical clusters built from a small number of sites can be identified directly from the SS rules. Figure 6 shows an SS ground state containing several of these dynamical clusters. For example, one has a single site that can be flipped while preserving all of the SS rules, representing an n c ¼ 1 dynamical cluster. A larger n c ¼ 2 cluster is also shown, where two spins can be flipped, though not independently. For both these examples we note that a large number of surrounding frozen spins are needed to construct these dynamical clusters. A naive counting for the n c ¼ 1 case yields a fraction of unfrozen to frozen spins of B1/ 25B4%, comparable to the few per cent average of unfrozen spins observed in our annealed simulations. These examples represent only a small subset of the possible dynamical clusters that can be constructed in the SS manifold. As described in Supplementary Note 3, there are several ways to construct dynamical clusters of arbitrary size (see Supplementary Fig. 4) as well as illustrations of the time evolution of dynamical clusters in simulations of small systems. The presence of such dynamical clusters is not specific to the single-spin-flip dynamics used; for example, analogous dynamical clusters can be constructed for spin-exchange dynamics, as illustrated in Supplementary Fig. 5, and we expect the same holds true for any local dynamics.

Discussion
Outside of any pure theoretical interest, one may be concerned with the fine-tuning required to reach the SS phase. As in SI 13 , though the precise point in phase space may be difficult to reach in a material realization, the nearby regions in phase space may be controlled primarily by the SS physics. Understanding the SS manifold then allows one to understand the surrounding phases and their higher temperature properties as perturbations to the SS model. Here we discuss two types of such perturbations: deviations from the J 2 ¼ J 3a ESI line and quantum terms, such as transverse field or exchange.
While the effects of finite second-and third-neighbour exchange on similar models has been studied extensively [24][25][26][27][28] , the regime along the ESI line and near the SS point remains largely unexplored. We find four neighbouring phases; the simplest are a (½, ½, ½) ordered phase expected from the J 3aþ N limit that appears for J 3a 4J/4 and a ferromagnetic SI state expected from the J 3a -À N limit that appears for J 3a oJ/4. For J 2 oJ/4 one finds a set of layered states (These are related to, but not identical to the layered states discussed for . with subextensive degeneracy B2 L . For J 2 4J/4, one finds a complex incommensurate ordering with wave-vector along [h00] or equivalents. The SS manifold ties these phases together, all of which are drawn from the SS ground-state manifold, with P i ¼ ± ½ for all pairs of tetrahedra, and extend over large regions of parameter space. We leave the detailed investigation of these neighbouring phases and other perturbations (such as J 3b , dipolar interactions and so on) for future studies.
The effect of quantum non-Ising interactions on SS is potentially much richer than in SI. In the latter, the addition of transverse exchange or transverse field induces tunnelling within the SI manifold yielding a U(1) quantum spin liquid [29][30][31][32] . This quantum spin liquid is described by an emergent electrodynamics, complete with a gapless photon excitation 29 . However, the associated energy scale of the quantum spin liquid is very small, due to tunnelling only appearing at high order in perturbation theory, confining its effects to very low temperatures and close proximity to the SI point 30,33 . In the SS, quantum dynamics appear at first order in perturbation theory (see Supplementary Note 4), and thus we expect them to be more significant than in SI. The presence of these first-order matrix elements is a direct reflection of the single-spin-flip and spinexchange dynamics of the SS manifold. Even with such mixing, when the perturbed Hamiltonian is projected into the SS manifold it still breaks up into infinitely many disconnected blocks, representing sets of states reachable by such local moves. The simplest blocks correspond to a small number of dynamical clusters well-separated by frozen regions. For example, there can be many n c ¼ 1 clusters as in Fig. 6, each with two states, corresponding to the freely flippable spin j "i and j #i for each cluster. Application of a transverse field $ À G P i s x i mixes the two states and gives a ground state of ðj "i þ j #iÞ= ffiffi ffi 2 p with energy gain of À G per dynamical spin. Other blocks correspond to more complicated dynamical clusters with more spatially extended structures. For example, for the large linear clusters discussed in Supplementary Note 3, the energy gain per dynamical spin is smaller, approaching approximately À 2G/n c for clusters of size n c (see Supplementary Note 4 and Supplementary Fig. 6 for an explicit example of this). More exotically, one can even construct states where a single dynamical cluster of size n c BO(N) encompasses nearly all of the spins in the system. Examples of such larger and more complex dynamical clusters are illustrated in Supplementary Movies 1 and 2. Similar considerations apply for transverse exchange A key difference is that odd-sized dynamical clusters are guaranteed to have degenerate ground states due to Kramers' theorem. In the exchange case, the n c ¼ 1 clusters thus contain free spins and gain no energy.
We thus conclude that for quantum SS, the ground states favoured at first order in perturbation theory will depend on the ground-state energies of this zoo of clusters as well as their effective packing fractions. We leave the detailed resolution of these non-trivial questions to future work. As this model is free of the sign problem, some of these questions should be addressable through quantum Monte Carlo simulations for both a ferromagnetic transverse exchange (J ± 40) or an arbitrary transverse field. The physics of the above dynamical clusters and the heterogeneous freezing could potentially enlighten our understanding of the phenomena of persistent dynamics in highly frustrated magnets 34 . In a more concrete setting, one may speculate that the SS could be connected to the physics observed in the quantum spin liquid candidate Tb 2 Ti 2 O 7 . A tantalizing clue are the shortrange correlations 35 at wave-vector (½, ½, ½) seen in Tb 2 Ti 2 O 7 and the (½, ½, ½) phase obtained by perturbing SS .
In summary, we have identified SS, a cooperative paramagnet on the pyrochlore lattice found by extending SI with further neighbour exchanges. This classical Ising model serves as a simple example of freezing and dynamical heterogeneity in a clean, disorder-free system. The features present in the magnetic correlations and the unusual low-temperature dynamics could prove useful in understanding such physics in real materials. We note that during the review of this article ref. 36, which also studies the ESI model and the SS point, appeared.

Methods
Monte Carlo simulations. For all Monte Carlo simulations, we used the standard Metropolis updating scheme with single-spin flip moves. For thermodynamic quantities, we simulated systems of N ¼ 16L 3 spins in L 3 conventional cubic unit cells of the pyrochlore lattice under periodic boundary conditions with linear size up to L ¼ 10. Typically, we used O(10 6 ) sweeps to anneal the system to each temperature and thermalize, then an additional O(10 6 ) sweeps were used to compute observables. Error estimates were computed via the bootstrap method. For spin-spin and charge-charge correlation functions, we simulated larger systems of size up to L ¼ 24, but only O(10 5 ) sweeps were needed to obtain accurate results. In both cases, we also used parallel tempering moves after each sweep to aid equilibriation. Longer simulations on smaller system sizes, with O(10 7 ) to O(10 8 ) sweeps produce results consistent with the shorter simulations on the larger systems. For dynamical quantities, a comparable number of sweeps and system sizes were used, except without the use of parallel tempering. To access the very low-  Figure 6 | Dynamical clusters in spin slush. We illustrate some of the dynamical clusters that can appear in the spin slush ground-state manifold. In (a) we show an example of part of a state with two such clusters, one containing a single dynamical spin (n c ¼ 1) and the other containing two dynamical spins (n c ¼ 2), with the dynamical spins and the surrounding charges highlighted in gold. In (b) we show the accessible states of an n c ¼ 1 dynamical cluster where the two states yield an average spin of zero and thus A i (N) ¼ 0. In (c) we show the n c ¼ 2 case, where one finds three accessible states with an average spin of ±1/3 and thus A i (N) ¼ A j (N) ¼ (1/3) 2 . In (b,c) the flippable spins for each state are highlighted in gold.
temperature auto-correlation function, we first slowly annealed the system to T/J ¼ 10 À 4 , guaranteeing that an SS ground state was reached, then followed the same protocol as the higher temperature simulations. This was repeated many times; two of these simulations are described in the main text.
Data availability. Raw data for any of the results reported in the text are available from the authors upon request.