# Tunable non-equilibrium dynamics: field quenches in spin ice

###### Abstract

We present non-equilibrium physics in spin ice as a novel setting which combines kinematic constraints, emergent topological defects, and magnetic long range Coulomb interactions. In spin ice, magnetic frustration leads to highly degenerate yet locally constrained ground states. Together, they form a highly unusual magnetic state – a “Coulomb phase” – whose excitations are pointlike defects – magnetic monopoles – in the absence of which effectively no dynamics is possible. Hence, when they are sparse at low temperature, dynamics becomes very sluggish. When quenching the system from a monopole-rich to a monopole-poor state, a wealth of dynamical phenomena occur the exposition of which is the subject of this article. Most notably, we find reaction diffusion behaviour, slow dynamics due to kinematic constraints, as well as a regime corresponding to the deposition of interacting dimers on a honeycomb lattice. We also identify new potential avenues for detecting the magnetic monopoles in a regime of slow-moving monopoles. The interest in this model system is further enhanced by its large degree of tunability, and the ease of probing it in experiment: with varying magnetic fields at different temperatures, geometric properties – including even the effective dimensionality of the system – can be varied. By monitoring magnetisation, spin correlations or zero-field Nuclear Magnetic Resonance, the dynamical properties of the system can be extracted in considerable detail. This establishes spin ice as a laboratory of choice for the study of tunable, slow dynamics.

## I Introduction

Nature and origin of unusual – in particular, slow – dynamics in disorder-free systems Stinchcombe2001 () are amongst the most fascinating aspects of disciplines as diverse as the physics of structural glasses and polymers Angell1995 (), chemical reactions and biological matter Tokuyama2004 (). Kinetically constrained models, following the original idea by Fredrickson and Andersen Fredrickson1984 (), represent one paradigm in which unusual dynamics is generated by short-distance ingredients alone without disorder Ritort2003 (). Another is provided by reaction-diffusion systems, in which spatial and temporal fluctuations feed off each other to provide a wide variety of dynamical phenomena Toussaint1983 () especially due to the slow decay of long wavelength fluctuations.

Spin ice systems Bramwell2001 () allow to combine both aspects – thanks to the nature of their emergent topological excitations, which take the form of magnetic monopoles Castelnovo2008 (); Castelnovo2012 () with long range Coulomb interactions. The ground state correlations in these localised spin systems lead to kinematic constraints in the reaction-diffusion behaviour of these mobile excitations Castelnovo2010 (); Levis2012 ().

Understanding the dynamics of spin ice systems, and in particular proposing new ways to probe their out-of-equilibrium properties, is of direct experimental relevance. For instance, modelling the emergent excitations near equilibrium Ryzhkin2005 (); Jaubert2009 () allowed to gain insight on the observed spin freezing at low temperatures Matsuhira2000 (); Snyder2004 (), and to explain in part the time scales measured in zero-field NMR Henley2013 (). Despite the fast paced progress, several open questions remain unanswered in particular concerning the behaviour of spin ice materials far from equilibrium, as evidenced for instance by recent low-temperature magnetic relaxation experiments Matsuhira2011 (); Yaraskavitch2012 (); Revell2012 ().

In this article, we study the strongly out-of-equilibrium behaviour following a quench from a monopole-rich to a monopole-poor regime by means of varying an applied magnetic field. We uncover a wide range of dynamical regimes and we provide a theoretical understanding of their origin. We show that the initial evolution maps onto a deposition of dimers on a honeycomb lattice, in presence of long range Coulomb interactions between the ‘vacancies’, i.e., the uncovered sites. The long time behaviour can instead be understood as a dynamical arrest due to the appearance of field-induced energy barriers to monopole motion. This regime can be seen as similar to conventional spin ice, but with a monopole hopping time exponentially sensitive to temperature: we have a Coulomb liquid in which time can pass arbitrarily slowly.

We discuss how to probe these phenomena in experiment, showing how by monitoring magnetisation, spin correlations or zero-field Nuclear Magnetic Resonance (NMR), the dynamical properties of the system can be extracted in detail.

Overall, this richness and versatility establishes spin ice as a laboratory of choice for the study of slow dynamics arising from an interplay of frustration (local constraints), topological defects (monopoles) and magnetic long-range Coulomb interactions.

## Ii Phase diagram and quench protocols

Spin ice materials consist of magnetic rare earth ions arranged on a corner-sharing tetrahedral (pyrochlore) lattice. Their magnetic moments are well-modeled by classical Ising spins constrained to point in the local direction (either into or out of their tetrahedra) Bramwell2001 (), as illustrated in the top left panel of Fig. 2 (see also the Suppl. Info.).

The frustration of the magnetic interactions manifests itself in the fact that there are not just a handful of ground states – an unfrustrated Ising ferromagnet, for instance, has only two ground states. Rather, spin ice has an exponentially large number of ground states, namely all those configurations which obey the ice rule that each tetrahedron has two spins pointing into it, and two out (2in-2out).

At any finite temperature, there will be excitations in the form of tetrahedra which violate the ice rule, having three spins pointing in and one out or vice versa. Crucially, these defects can be thought of as magnetic monopoles which are deconfined and free to carry their magnetic charge across the system Castelnovo2008 (). The energy cost associated with the creation of a monopole is of the order of a few degrees Kelvin for spin ice materials DyTiO and HoTiO, leading to an exponential suppression of their density for .

We refer the reader interested in the detailed background to Ref. Castelnovo2012 (); Henley2010 (). Here, we emphasize that spin ice is the only experimentally accessible magnetic system that realises a Coulomb phase described by a low-energy emergent gauge field. The new phenomena we describe here can ultimately be traced back to this fundamentally new feature.

### ii.1 The equilibrium phase diagram of spin ice

Four distinct regimes are encountered in spin ice in presence of a field along a [111] crystallographic direction, as displayed in Fig. 1 Bramwell2001 (). In order to understand this phase diagram, it is convenient to divide the spin lattice (pyrochlore, or corner-sharing tetrahedral lattice) into alternating kagome and triangular layers, as illustrated in Fig. 2. Further details on spin ice systems as well as on their phase diagram in a field can be found in the Suppl. Info.

In the limit of strong fields, in the saturated (SatIce) regime, all the spins point along the field direction, while respecting the local easy axes (see Fig. 2, top middle panel). The ice rules are violated everywhere and each tetrahedron hosts a monopole; the monopoles form an “ionic crystal”.

As the field strength is reduced, violations of the ice rules are no longer sufficiently offset by a gain in Zeeman energy and a regime where most tetrahedra obey the ice rule is recovered for . This necessarily requires some of the spins to point against the applied field. At intermediate field strengths, they will dominantly be the spins in the kagome planes, as their Zeeman energy is smaller by a factor of compared to the spins in the triangular planes. This leads to an extensively degenerate regime known as kagome ice (KagIce), illustrated in the bottom middle panel in Fig. 2.

At low field strengths, the kagome ice regime becomes entropically unstable to the ‘conventional’ spin ice regime, namely the ensemble of all configurations satisfying the ice rules irrespective of the polarisation of the triangular spins (top right panel in Fig. 2). All of these regimes crossover at sufficiently large temperatures into a conventional paramagnetic regime.

These different regimes appear in the field-temperature phase diagram of spin ice as illustrated in Fig. 1. Note that only the KagIce and SatIce phases are separated at low temperatures by an actual (first order) phase transition Hiroi2003 (); Sakakibara2003 () (see Suppl. Info. for further details).

### ii.2 Field quench protocols and non-equilibrium phenomena

In this article we discuss field quenches from SatIce to KagIce, as exemplified by the blue vertical arrows in Fig. 1. We focus on field quenches at (constant) temperature , as the system appears to relax straightforwardly to equilibrium at higher temperatures (blue dash-shaded region in Fig. 1) on time scales of the order of a typical microscopic spin-flip time scale.

We are interested here mostly in the low-temperature dynamics that follows quenches across the first order phase transition in Fig. 1. As we shall see in the following, these quenches are characterised by interesting slow relaxation phenomena and dynamical arrest. By comparison, quenches at higher temperatures appear to be fast and featureless. The presence of a critical point in the phase diagram at these higher temperatures raises a number of separately interesting questions and ideas for further work, such as the possibility of quenching to / across a critical point, and testing whether the critical correlations bear any signatures in the resulting dynamical behaviour. Whilst interesting in their own right, these issues are beyond the scope of the present work.

The dynamics of spin ice at K appears to be based on incoherent spin reversals, well modelled by a Monte Carlo (MC) single-spin flip dynamicsJaubert2009 (). From AC susceptibility measurements Snyder2004 (), a Monte Carlo time step corresponds to approximately 1 ms in DyTiO, which we use to convert MC time into physical time in this work.

We used Monte Carlo simulations with single spin flip dynamics, with DyTiO spin ice parameters as in Ref. Melko2004 (). We employed Ewald summation techniques deLeeuw1980 () to treat the long range dipolar interactions, and the Waiting Time Method WTM_refs () to access long simulation times at low temperatures (see also Ref. Castelnovo2010 ()).

Observables that we monitor after a field quench are monopole density and magnetisation, as illustrated in Fig. 3.

In addition to gaining insight on non-equilibrium dynamical phenomena in spin ice – in particular on the role of the emergent magnetic monopoles – our aim is to investigate whether a protocol can be devised that allows to prepare spin ice samples in a low-temperature yet monopole-rich state by means of a rapid reduction of an applied magnetic field. Attaining such state would be a significant step towards a more direct experimental detection of the magnetic monopoles, for instance using zero-field NMR techniques Sala2012 (), as we discuss towards the end of the article.

The use of field quenches instead of thermal quenches Castelnovo2010 () offers two practical advantages: (i) a fast variation in the externally applied magnetic field is experimentally less challenging than performing a controlled thermal quench; and (ii) in quenches from SatIce to KagIce, the monopole density is to an excellent approximation proportional to the magnetisation of the sample, a thermodynamic quantity, and therefore directly and simply accessible in experiments.

In order to make contact with a broader range of experimental settings, we briefly review towards the end of the paper alternative experimental protocols, such as quenches to zero field, as well as slow “quenches” (i.e., finite-rate field ramps). For completeness, these are discussed in greater detail in the Supporting Information.

We point out that in the following every reference to the magnetic field is intended as the externally applied field . Furthermore, we note that Monte Carlo simulations are devoid of demagnetisation corrections, which ought to be taken into account when comparing to possible experimental results.

## Iii Saturation to kagome ice

As illustrated in Fig. 2, KagIce comprises the hardcore dimer coverings of the honeycomb lattice. The equilibrium density of monopoles in this regime is exponentially small at low temperatures, and they correspond to monomers in the dimer language. SatIce instead corresponds to an ionic crystal of monopoles, i.e., a configuration without any dimers. Microscopically, in SatIce pairs of neighbouring positive and negative monopoles on the same kagome plane annihilate by inverting the intervening spin, which is equivalent to the deposition of a dimer on the corresponding bond in the honeycomb lattice.

Field quenches from SatIce to KagIce can thus be understood as a stochastic deposition process of dimers, with non-vanishing desorption probability, in presence of long range (3D) Coulomb interactions between the monomers (i.e., the monopoles). Eventually the dimers become fully packed, up to an exponentially small concentration of thermally excited monomers at equilibrium.

The physics of dimer deposition processes with only short range interactions was recently considered in the experimental context of molecular random tilings Blunt2008 (), which were shown to exhibit a variety of different (stochastic) dynamical regimes akin to kinetically constrained models Stinchcombe2001 (); Garrahan2009 (). In the following we show how the presence of long range interactions, as well as additional kinematic constraints in the motion of monopoles due to the underlying spin configuration, give rise to remarkably rich behaviour.

### iii.1 Short-times: dimer deposition

The initial spin flip events following a quench from SatIce to KagIce lower the energy of spin interactions by removing monopoles, but they increase the Zeeman energy. As the latter is proportional to the quench field, at weak fields, all flips are initially downhill in energy and thus take place on a time scale of the order of 1 MC step, independently of temperature. A mean-field equation capturing the initial decay of the monopole density is

(1) | |||||

(2) |

where is the monopole density per tetrahedron and the factor 3 appears because each monopole can equally annihilate with three of its neighbours. This is in good agreement with the simulation results in Fig. 4.

For quench fields above a threshold of Tesla, the initial decay becomes activated, as the Zeeman energy for reversing a kagome spin K exceeds the spin-spin interaction energy gain of about K (DyTiO parameters) footnote1 (). For fields beyond this threshold, the initial decay is described by a dimer deposition process with field-dependent energy barriers.

The near-neighbour monopole pair annihilation continues until no suitable pairs are left. In dimer deposition language, this corresponds to the random packing limit where it is no longer possible to deposit futher dimers without first rearranging those already present on the lattice.

For a wide range of fields and temperatures, the initial decay takes approximately MC steps (see Fig. 3) and reaches a monopole density in the range 0.05 - 0.08, depending on field and temperature.

We compare this value to that obtained from numerical simulations of random depositions of dimers on the honeycomb lattice, where we measure the density of vacancies at the end of each deposition process when it is no longer possible to add a (hard-core) dimer to the system without rearranging the ones already deposited. The density of leftover vacancies in the non-interacting simulations equals . It is interesting to notice that this value is at least % larger than that following field quenches from SatIce to KagIce. We find that the difference is due to the combination of a non-vanishing desorption probability and long range Coulomb attraction between oppositely charged monopoles, which are present in the spin ice quenches and not in the random deposition process.

### iii.2 Long times: from noncontractible pairs to dynamically-obstructed monopole diffusion

At the end of the initial decay, no neighbouring monopoles are left in any of the kagome planes. The majority of monopoles are ‘isolated’ and only a minority forms (noncontractible) pairs across adjacent planes (3).

Two oppositely-charged monopoles on neighbouring tetrahedra form a noncontractible pair whenever inverting the direction of the shared spin leads to an even higher energy configuration with the two tetrahedra in the 4in and 4out state, respectively Castelnovo2010 (). For instance, every monopole pair across a triangular spin in SatIce is noncontractible (see Fig. 2). Indeed, in KagIce, where the triangular spins are polarised in the field direction, noncontractible pairs can only form between positive and negative monopoles belonging to adjacent kagome planes.

For small values of the applied field ( Tesla), we see from Fig. 3 that the decay in the density of the isolated monopoles is fast in comparison to the decay of noncontractible pairs, and eventually only the latter are left to determine the long time behaviour of the system. This regime can be understood in analogy to thermal quenches Castelnovo2010 () (see Suppl. Info.).

Whereas isolated monopoles are essentially free to move, noncontractible pairs have an activation barrier to move/decay, and thus their lifetime increases exponentially with inverse temperature. This explains why the isolated monopoles decay faster than noncontractible pairs; it also explains the mechanism responsible for the long time decay in the overall monopole density of the system. Whereas lattice-scale kinematic constraints control the formation (or possibly the survival from the initial SatIce state) of noncontractible pairs, it is their decay that ultimately controls the long time behaviour of the monopole density.

The life time of a noncontractible pair is determined by the energy barriers of the allowed decay processes (namely, processes where the two charges separate and annihilate one another elsewhere on the lattice; or they may annihilate with other charges in the system). For quenches to zero field, the energy barriers are solely due to spin-spin interactions and the corresponding decay occurs via a process, identified in thermal quenches Castelnovo2010 (), whereby the noncontractible pair recombines by hopping the long way around a hexagon containing the obstructing spin.

The presence of a small applied field does not alter the picture qualitatively, but it introduces two important differences. Firstly, the mechanism and barrier to the decay of noncontractible pairs is altered. Once the majority of the triangular spins are polarised, it is straightforward to see that the hexagonal processes invoked in Ref. Castelnovo2010, are no longer available: decay must proceed by separation of the monopoles until they meet oppositely charged ones from other noncontractible pairs. This process incurs a Coulomb energy barrier for separating the positive and negative charges in a noncontractible pair, which grows with the distance to the next monopole in the plane. Moreover, the presence of an applied magnetic field also affects the value of the barrier to separation, as the Zeeman energy for spin flips needs to be taken into consideration along with the Coulomb attraction of the oppositely charged monopoles.

Secondly, isolated monopoles hopping within kagome planes exhibit activated behaviour of their own due to the Zeeman energy (see below for details). This explains the opening of an intermediate time window between the end of the initial decay, discussed in the previous section, and the final decay controlled by noncontractible pairs. This time window is absent in thermal quenches Castelnovo2010 () as well as in field quenches down to zero field (see Suppl. Info.).

This intermediate-time regime is governed by diffusion-annihilation processes of isolated monopoles across larger and larger distances as their density is reduced. Depending on the value of the field and temperature, our simulations show different behaviours following the initial decay, including an apparent power-law at low fields (see left panel in Fig. 3). This regime is controlled by “finite time” processes rather than by any asymptotic behaviour. Albeit challenging to model in the absence of an asymptotic regime, it is a novel and interesting example of a reaction-diffusion process in presence of long range Coulomb interactions and kinematic constraints that is experimentally accessible in spin ice materials.

As the field strength increases, the decay of the isolated monopoles gets comparatively slower, until for Tesla they remain the majority with respect to noncontractible pairs at all time. In this regime, a new mechanism controls the long time behaviour: the system remains always in a Coulomb liquid phase, rather than condensing its monopoles into noncontractible pairs. It is the isolated monopoles themselves that, by becoming exceedingly slow, determine the long time decay of the system.

In order to better understand and confirm this scenario, let us look in detail at how such slowing down of the monopole time scales takes place. Consider the motion of a monopole within a kagome layer (the illustrative Fig. S3 provided in the Suppl. Info. may be of help here). Ordinary monopole motion involves alternatively flipping pairs of spins with negative and positive projection onto the applied field. This yields a periodic energy landscape due to the Zeeman energy change for the kagome spins, where is measured in Tesla and the energy is measured in degrees Kelvin. As the field increases from zero, this barrier progressively slows down the monopoles.

For sufficiently large applied fields, another process becomes energetically more convenient: A new monopole-antimonopole pair is first excited neighbouring the existing monopole, which then annihilates with the opposite member of the pair (pair-assisted hopping – see figure in the Supplementary Material). The energy barrier for this process is , where K is the cost of flipping a spin in spin ice and K is the Coulomb energy difference between two monopoles at nearest-neighbour vs second-neighbour distance (DyTiO parameters). The new process becomes energetically more favourable than the ordinary Zeeman barrier if Tesla. Notice that the sign of the Zeeman contribution in the new process is reversed: if the fields were increased significantly beyond the threshold, monopole diffusion would once again become a fast process. However, this is forbidden by an intervening first order phase transition Castelnovo2008 () which already affects field quenches to Tesla.

The largest quench field we could study was Tesla. At this value, the two processes above give rise to barriers of K and K, respectively. These values ought to be corrected for and broadened by quadrupolar terms that are missing in the monopole picture used in the estimates Castelnovo2008 (). Indeed, the long time decay of the curves in Fig. 3, collapses upon rescaling the time axis by a Boltzmann factor , as illustrated in Fig. 5.

The corresponding value of the activation barrier, K, is in reasonable agreement with the estimates for the barrier to isolated monopole motion.

We stress that for Tesla, the system enters a long-lived Coulomb liquid phase with an enhanced (out-of-equilibrium) density of exceedingly slow monopoles. This is a regime that may be experimentally promising to probe and detect single monopoles, as we discuss farther below.

## Iv Experimental probes

We have argued that different quench regimes lead to the appearance of a variety of non-equilibrium dynamical phenomena. The next question to ask is how best to access such quantities experimentally. On the one hand, a solid state system like spin ice does not permit direct imaging of the magnetic state in the bulk, unlike two-dimensional systems such as artificial spin ices artificial_spin_ice_refs (). On the other hand, we have access to a number of well-developed probes; this being a spin system, natural options are magnetisation and NMR.

### iv.1 Magnetisation

One of the advantages of using field quenches from SatIce to KagIce is the fact that we have direct access to the monopole density per tetrahedron by measuring the magnetisation in the direction (in units of Dy): the equation is extremely accurate for Tesla, and it remains very useful also at lower field strengths, with deviations due to reversed triangular spins of less than % down to Tesla (see Supporting Information), for K.

### iv.2 Nuclear magnetic resonance

A promising approach to measure the monopole density in spin ice samples follows from the zero-field NMR measurements pioneered by the Takigawa group TakigawaNMR (); Sala2012 (). This technique uses DyTiO samples with NMR active O nuclei at the centres of the tetrahedra. The large internal fields induced by the rare earth moments at these locations ( Tesla) result in a peak in the NMR signal from the O ions in zero external field.

When a monopole is present in a tetrahedron, the violation of the ice rules results in a sizeable reduction of the field at the centre of a tetrahedron ( Tesla). This in turn ought to give rise to an additional shifted peak in the NMR signal. The relative intensity of the two peaks can be used for instance to extract the magnetic monopole density in the sample Sala2012 ().

One of the challenging aspects of the zero-field NMR approach to detect monopoles is the trade-off between having a large number of monopoles to enhance the signal and having sufficiently long persistence times (the time a monopole spends in a given tetrahedron) to reduce the noise in the signal. High-temperatures cannot be used to increase the monopole population because they lead to fast magnetic fluctuations. On the other hand, the system cannot be allowed to thermalise at low temperature as this would result in an exponentially small – thus undetectable – monopole density.

The field quenches discussed in this paper can be used to prepare a spin ice sample in a strongly-out-of-equilibrium state that offers great potential to NMR measurements. Indeed, the long-time regime at low temperatures ( K) has a sizeable density of long-lived, static monopoles (see Fig. 3).

In Fig. 6 we show the density of monopoles that are present in the system at a given time and that have remained in their current tetrahedron for longer than (persistence time).

Clearly, and , . Fig. 6 confirms that the long time decay of the monopole density is indeed due to long persistence times of the residual (essentially static) monopoles rather than mobile monopoles which are somehow unable to annihilate each other.

### iv.3 Alternative experimental protocols

In practice, depending on details of the set-up, other protocols may be easier to implement than sudden quenches to non zero field values. For instance, it may be more practical to remove a field altogether – quenching to SpinIce rather than to KagIce – or to lower the field gradually to allow time to dissipate the field energies. We present here only a brief summary of the results. Details can be found in the relevant sections of the Supporting Information.

We find that zero-field quenches lack the intermediate-time regime between the end of the dimer deposition process and the onset of the long time decay controlled by noncontractible pairs (this intermediate-time regime is evident for instance in Fig. 3). Once again, one concludes that zero field quenches can be used to prepare spin ice samples in a metastable state rich in monopoles forming noncontractible pairs, albeit the process appears to be about one order of magnitude less efficient than SatIce to KagIce quenches with respect to the resulting density.

Replacing an instantaneous field quench (on the scale of the microscopic time MC step ms) with a constant-rate field sweep clearly has a profound effect on the short time behaviour – it corresponds to an annealed dimer deposition process, which would itself be an interesting topic for further study. However, shortly after the field reaches its final value, MC simulations show that the evolution of the monopole density agrees quantitatively with the behaviour following a field quench to the same field value.

Therefore, any experimental limitations to achieve fast variations of applied magnetic fields affect only the behaviour of the system while the field is varied, whereas the main features discussed above occurring on a longer time scale remain experimentally accessible.

## Acknowledgments

This work was supported in part by EPSRC grants EP/G049394/1 and EP/K028960/1 (C.C.). This work was in part supported by the Helmholtz Virtual Institute “New States of Matter and Their Excitations”. C. Castelnovo acknowledges hospitality and travel support from the MPIPKS in Dresden. We are grateful to M. Takigawa and K. Kitagawa, R. Kremer, S. Zherlytsin and to T. Fennell for discussions on experimental probes of field quenches in spin ice.

## V Conclusions

Spin ice offers a realization of several paradigmatic concepts in non-equilibrium dynamics: dimer absorption, Coulombic reaction-diffusion physics and kinetically constrained slow dynamics. There is an unusually high degree of tuneability: the timescale of the elementary dynamical move through a Zeeman energy barrier; the dimensionality of the final state (d = 2 KagIce vs d = 3 SpinIce); and the relative importance of dimer desorption compared to Coulomb interactions. The additional availability of a range of experimental probes thus allows broad and detailed experimental studies.

With the advent of incipient thermal ensembles of artificial spin ice, where the constituent degrees of freedom can even be imaged individually in real space, there even appears an entirely new setting for the study of these phenomena on the horizon thermal_ASI ().

## References

- (1) Stinchcombe R (2001) Stochastic non-equilibrium systems. Adv Phys 50: 431-496.
- (2) Angell C A (1995) Formation of Glasses from Liquids and Biopolymers. Science 267: 1924-1935.
- (3) 8th Tohwa University International Symposium, Fukuoka, Japan, November 1998. Slow Dynamics in Complex Systems, eds Tokuyama M and Oppenheim I (Springer), AIP Conference Proceedings, Vol. 469 (2004).
- (4) Fredrickson G H, Andersen H C (1984) Kinetic Ising-model of the glass-transition. Phys. Rev. Lett. 53: 1244-1247.
- (5) Ritort F, Sollich P (2003) Glassy dynamics of kinetically constrained models. Adv Phys 52: 219-342.
- (6) Toussaint D, and Wilczek F (1983) Particle-antiparticle annihilation in diffusive motion. J Chem Phys 78: 2642-2647.
- (7) Bramwell S T and Gingras M J P (2001) Spin ice state in frustrated magnetic pyrochlore materials. Science 294: 1495-1501.
- (8) Castelnovo C, Moessner R, and Sondhi S L (2008) Magnetic monopoles in spin ice. Nature 451: 42-45.
- (9) Castelnovo C, Moessner R, and Sondhi S L (2012) Spin Ice, Fractionalization, and Topological Order. Annu Rev Condens Matter Phys 3: 35-55.
- (10) Castelnovo C, Moessner R, and Sondhi S L (2010) Thermal quenches in spin ice. Phys. Rev. Lett. 104: 107201-4.
- (11) Levis D and Cugliandolo L F (2012) Out of equilibrium dynamics in the bidimensional spin-ice model. EPL 97: 30002-6; Levis D and Cugliandolo L F (2013) Defects dynamics following thermal quenches in square spin-ice. Phys. Rev. B 87: 214302-14.
- (12) Ryzhkin I A (2005) Magnetic Relaxation in Rare-Earth Oxide Pyrochlores. JETP 101: 481-486.
- (13) Jaubert L D C and Holdsworth P C W (2009) Signature of magnetic monopole and Dirac string dynamics in spin ice. Nat Phys 5: 258-261.
- (14) Matsuhira K, Hinatsu Y, Tenya K, and Sakakibara T (2000) Low temperature magnetic properties of frustrated pyrochlore ferromagnets HoSnO and HoTiO. J Phys: Condens Matter 12: L649-656.
- (15) Snyder J et al. (2004) Low-temperature spin freezing in the DyTiO spin ice. Phys. Rev. B 69: 064414-6.
- (16) Henley C L (2013) NMR relaxation in spin ice due to diffusing emergent monopoles. arXiv:1210.8137.
- (17) Matsuhira K, et al. (2011) Spin Dynamics at Very Low Temperature in Spin Ice DyTiO. J Phys Soc Japan 80: 123711-4
- (18) Yaraskavitch L R, et al. (2012) Spin dynamics in the frozen state of the dipolar spin ice material DyTiO. Phys. Rev. B 85: 020410(R)-5.
- (19) Revell H M, et al. (2012) Evidence of impurity and boundary effects on magnetic monopole dynamics in spin ice. Nat Phys 9: 34-47.
- (20) Henley C L (2010) The Coulomb phase in frustrated systems. Ann Rev Condens Matter Phys 1: 179-210.
- (21) Hiroi Z, Matsuhira K, Takagi S, Tayama T, and Sakakibara T (2003) Specific heat of kagome ice in the pyrochlore oxide DyTiO. J Phys Soc Japan 72: 411-418.
- (22) Sakakibara T, Tayama T, Hiroi Z, Matsuhira K, and Takagi S (2003) Observation of a liquid-gas-type transition in the pyrochlore spin ice compound DyTiO in a magnetic field. Phys. Rev. Lett. 90: 207205-4.
- (23) Melko R G and Gingras M J P (2004) Monte Carlo studies of the dipolar spin ice model. J Phys: Condens Matter 16: R1277-R1319.
- (24) de Leeuw S W, Perram J W, and Smith E R (1980) Simulation of electrostatic systems in periodic boundary conditions. I: Lattice sums and simulation results. Proc R Soc London Ser A 373: 27-56.
- (25) Dall J and Sibani P (2001) Faster Monte Carlo simulations at low temperatures. The waiting time method. Computer Phys Comm 141: 260-267.
- (26) Sala G, et al. (2012) Magnetic Coulomb Fields of Monopoles in Spin Ice and Their Signatures in the Internal Field Distribution. Phys. Rev. Lett. 108: 217203-4.
- (27) Blunt M O, et al. (2008) Random tiling and topological defects in a two-dimensional molecular network. Science 322: 1077-1081.
- (28) Garrahan J P, Stannard A, Blunt M O, and Beton P H (2009) Molecular random tilings as glasses. Proc Natl Acad Sci USA 106: 15209-15213.
- (29) For reviews, see Heyderman L J and Stamps R L (2013) Artificial ferroic systems: novel functionality from structure, interactions and dynamics. J Phys: Condens Matter 25, 363201; Nisoli C, Moessner R, and Schiffer P (2013) Artificial spin ice: Designing and imaging magnetic frustration. Rev Mod Phys 85, 1473-1490.
- (30) Kitagawa K and Takigawa M, poster at the ESF-HFM meeting “Topics In the Frustration of Pyrochlore Magnets” in Abingdon, UK (2009).
- (31) Porro J M, Bedoya-Pinto A, Berger A, and Vavassori P (2013) Exploring thermally induced states in square artificial spin-ice arrays. New J Phys 15: 055012-12; Zhang S, et al. (2013) Crystallites of magnetic charges in artificial spin ice. Nature 500: 553-557; Farhan A, et al. (2013) Exploring hyper-cubic energy landscapes in thermally active finite artificial spin-ice systems. Nat Phys 9: 375-382; Farhan A, et al. (2013) Direct Observation of Thermal Relaxation in Artificial Spin Ice. Phys. Rev. Lett. 111: 057204-5.
- (32) The energy gain of K for reversing a kagome spin in SatIce was obtained by extrapolating MC values to infinite system size. It is in reasonable agreement with the value K that obtains from the effective description in terms of magnetic monopoles, where is the Madelung constant of the Zincblende structure (an ionic diamond lattice), K is the energy cost of a monopole and K is the energy of two neighbouring monopoles (DyTiO parameters).

Supplementary Information

## Appendix A Dipolar spin ice model

We present here a brief summary of the model and properties of spin ice. For a more thorough review, see for instance Ref. Bramwell2001, and Ref. Castelnovo2012, .

Spin ice models and materials are systems where the magnetic degrees of freedom are localised on a lattice of corner sharing tetrahedra (pyrochlore lattice), illustrated in Fig. S1.

Given a choice of direction (i.e., one of the major diagonals of the cube), the pyrochlore lattice can be seen as a layered structure of alternating triangular and kagome planes of spins. This is particularly important in the present manuscript where we will study the behaviour of the system in presence of a magnetic field applied in the direction, which results in a crucial difference between the corresponding triangular and kagome spins.

The largest energy scale in spin ice systems is a single ion anisotropy that forces the spins to lie along the axis connecting the centre of a tetrahedron to the corresponding vertex, as illustrated in Fig. S2. This energy scale is usually so large compared to interactions as well as the relevant temperature and field energy scales that one can model the spins as strictly easy-axis (Ising).

The spins interact via exchange and dipolar coupling, which in these systems happen to be of approximately the same magnitude at nearest-neighbour distance. The model used in this paper is known as the dipolar spin ice model Shastry1999 (). If we label the sites of the three-dimensional pyrochlore lattice, and the local easy-axis , each spin can be represented as a classical vector , where is the size of the magnetic moment and .

The thermodynamic properties of spin ice are described with good accuracy by a Hamiltonian energy that encompasses a uniform nearest-neighbour exchange term of strength K and long-range dipolar interactions,

(3) |

where Å is the pyrocholore lattice constant, and K.

Due to a peculiar interplay between interactions, lattice geometry and local easy-axis anisotropy, it is not possible to minimise simultaneously each term in the Hamiltonian and the system is frustrated. To a first approximation, the energy is minimised by configurations where each tetrahedron has two spins pointing in and two pointing out (2in-2out, Fig. S2) Isakov2005 (). These configuration are in one-to-one correspondence with proton disorder in (cubic) water ice, hence the name spin ice and the fact that the 2in-2out rules are referred to as ice rules.

There are many ways to fulfill the ice rule condition on the pyrochlore lattice, which result in an extensive degeneracy. The ensamble of these configurations is neither ordered nor trivially disordered; it is perhaps best understood as an instance of classical topological order, whereby the system develops an additional (gauge) symmetry at low temperatures Castelnovo2012 ().

Excitations over these low energy states take the form of defective tetrahedra with one spin pointing in and three out, or vice versa (3out-1in, Fig. S2). In recent years it was shown that these defects behave as Coulomb-interacting point like magnetic charges (i.e., a Coulomb liquid) free to move in three dimensions Castelnovo2008 ().

At closer inspection, the 2in-2out configurations are not exactly degenerate (in presence of dipolar interactions). Their splitting is however much smaller than the strength of the interactions , and a conventional ordering transition is not observed in the system down to a correspondingly small temperature (only numerical evidence of this transition is available to date, though see Ref. Pomaranski2013, for some preliminary experimental results). Therefore, these systems exhibit a relatively broad temperature range between the cost of a monopole excitation, , down to the transition temperature , in which their behaviour is captured at a coarse grained level by a Coulomb liquid of magnetic charges. This description has indeed gone a long way in allowing us to understand the behaviour observed experimentally in spin ice materials.

## Appendix B Phase diagram in a magnetic field

The effect of an external magnetic field applied on a spin ice system along one of the axes is not the same on all spin sublattices. As it is evident from Fig. S1 and Fig. S2, the field is aligned with the local easy axis of one spin per tetrahedron and canted with respect to the other three spins, which thus incur a lesser Zeeman coupling. In the alternating kagome-triangular layer interpretation of the pyrochlore lattice, the parallel spins live on the triangular planes, whereas the canted spins live on the kagome planes (see Fig. S1).

When the field energy is stronger than spin-spin interactions, the lowest energy state where each spin has a positive component in the direction of the field has each tetrahedron either in a 1in-3out or in a 3in-1out configuration (an example of such configuration is shown in the right panel of Fig. S2).

This configuration corresponds to a perfect ionic crystal of magnetic charges living at the centres of the tetrahedra (which form a bipartite diamond lattice).

As the field is reduced, the interactions favour the restoration of the 2in-2out ice rules. However, since one out of four spins is more strongly Zeeman coupled to the field, the ice rules will be restored predominantly via re-arrangement of the canted spins. Namely, the spins in the triangular planes remain essentially fully polarised whereas one every three spins in the kagome layer gets reversed to point opposite to the field direction (left panel in Fig. S2).

These new configurations that occur at intermediate field values are not as disordered as generic 2in-2out spin ice, yet they are only partially polarised and they remain extensively degenerate. They are referred to as kagome ice states. These states are, up to thermal fluctuations, devoid of monopole excitations.

It is interesting to re-interpret the role of the applied field in terms of the Coulomb liquid picture of the excitations of the spin ice ground states. From vanishing to intermediate fields, the system is able to respond without violating the ice rules, and a continuous crossover is observed from spin ice to kagome ice. When the field is further increased, it begins to compete with the ice rules, as it favours an ionic crystal of magnetic charged. Namely, the applied field in the spin language translates into a staggered chemical potential in the monopole picture.

The phase diagram of a liquid of Coulomb interacting charges as a function of temperature and staggered chemical potential is the well-known liquid-gas phase diagram. This is indeed a rather distinctive phase diagram, characteristic of itinerant systems and long range interactions, which is quite at odds (unprecedented!) with the behaviour of a localised spin system. The experimental observation of such phase diagram in spin ice materials was indeed one of the first smoking guns that the theoretical proposal for emergent magnetic monopole excitations in these systems was indeed correct Castelnovo2008 ().

Details about this phase diagram can be found in Ref. Castelnovo2008, (Fig. 4) and it is schematically represented in Fig. 1 in the main text. At low temperatures, the high-field polarised state (saturated ice) is separated from the intermediate field kagome ice state by a first order transition. In the monopole language, this transition connects a low-density (kagome) to a high-density (saturated) monopole phase. In the spin language, the order parameter describing the transition is the magnetisation, which jumps discontinuously from a low to a high moment. As discussed in the main text, up to thermal fluctuations in the triangular spins (which remain essentially polarised throughout the transition), the monopole density and the magnetisation of the system are directly related ().

At some finite temperature, the line of first order points terminates at a characteristic critical end point (Coulombic criticality), and at higher temperatures the transition is replaced by a continuous crossover from kagome to saturated ice.

## Appendix C Field quenches in nearest-neighbour spin ice

In view of the advent of “exchange spin ice” materials, where the Coulomb interactions between monopoles is relatively weak Kimura2013 () on account of the super exchange between neighbouring spins dominating over their dipolar interactions, it is instructive to study the long time behaviour of the monopole density in field quenches in the simpler context of nearest-neighbour spin ice. There, the lack of dipolar (or further range exchange) interactions implies that the monopoles are trivially deconfined and no metastable bound states can form. For field quenches down to zero field, this removes all possible activation barriers and the monopole density is controlled by the asymptotics of reaction diffusion processes (see Ref. Castelnovo2010 ()).

The situation is different in field quenches to finite fields, where there exist Zeeman energy barriers to monopole motion which thus affect the long time decay of the mononopole density. In the discussion hereafter we shall assume for convenience that the triangular spins in the system are fully polarized, which is indeed the relevant scenario for the long time decay of the monopole density.

Straightforward hopping of a monopole is object to an alternating energy landscape due to the Zeeman energy to flip a kagome spin, (see main text and Fig. S3). As this energy scale grows with increasing field, an alternative process becomes energetically preferable: the creation of a new monopole pair neighbouring the existing monopole, with subsequent annihilation of the latter with its opposite member of the pair. The energy scale for this process in nearest-neighbour spin ice is , which decreases with increasing field. The two processes are illustrated schematically in Fig. S3.

For the typical value of K denHertog2000 (), the two energy scales cross over at Tesla. For we expect that the long time decay of the monopole density is controlled by the energy barrier ; for we expect that the relevant barrier is instead . Fig. S4 confirms that the corresponding Boltzmann factors lead to a very good collapse of the long time behaviour of the monopole density in nearest-neighbour spin ice.

Note that the very good collapse of the data from nearest-neighbour simulations suggests that the entropic long range Coulomb interactions that are present even in nearest-neighbour spin ice do not play a significant role in the long time behaviour of the monopole density at these temperatures.

The introduction of dipolar interactions in the system changes the pair-assisted energy scale (even for an isolated monopole), due to the Coulomb interaction between the three charges in the intermediate stage. In addition to the energy cost of creating two new monopoles, we need to take into account two Coulomb interaction terms between opposite charges at nearest-neighbour distance and one term between like charges at second neighbour distance. An estimate of the resulting barrier for typical DyTiO parameters is given in the main text. The dependence on the magnetic field by contrast remains unchanged.

## Appendix D The triangular spins

Much of our interpretation of the behaviour of spin ice following a field quench from SatIce to KagIce in this paper relies on the triangular spins being almost fully polarised. We have verified this assumption explicitly in our numerical simulations (see Fig. 3 in the main text and Fig. S5).

Even at low fields, only a small fraction of triangular spins is temporarily reversed during the out of equilibrium relaxation. In all regimes considered in this manuscript (with the exception of quenches down to zero field, to be discussed next), to within % the magnetisation tracks the monopole density as , in units of Dy: the magnetisation is an excellent proxy for the monopole density.

In the following, we discuss what can be learnt from the behaviour of the triangular planes alone. First of all, for fields Tesla, no reversed triangular spin can be present in the system for any length of time ( MC step), as it is unstable to flipping and creating two monopoles: the Zeeman splitting of a triangular spin K exceeds the energy cost of monopole pair creation, K for DyTiO.

When Tesla, reversed triangular spins can be in a (long-lived) metastable state. Our simulations show that this occurs mostly towards the end of the initial decay and throughout the intermediate decay (see Fig. 3); the densities remain negligibly small even at very low fields: less than % of the triangular spins are reversed at any given time during, say, an Tesla field quench (see Fig. S5).

For Tesla, while it costs energy to flip a reversed triangular spin and create two monopoles, the system can then readily lower its energy by moving the two monopoles one step each. In this field range, the Coulomb cost to separate the monopoles is offset by the Zeeman energy gain so that the final state has lower energy than the starting one before we flipped the triangular spin. This process incurs an energy barrier of approximately K at Tesla (estimated in the Coulomb liquid approximation), which is in reasonable agreement with the value K obtained by collapsing the long time behaviour in the right panel of Fig. S5.

Naturally, misaligned triangular spins have another relaxation channel, whereby a monopole moves from one kagome plane to the next, by flipping the triangular spin without incurring an energy barrier. This process is controlled by the time scale for a monopole to come by, set by (among other items) the monopole hopping barrier within a kagome plane (see main text). For Tesla, said barrier is K, which is clearly less efficient than the former (and indeed in worse agreement with the data).

It is perhaps more interesting to investigate the triangular spin relaxation at lower fields, Tesla, yet large enough so that the equilibrium density of reversed triangular spins is negligibly small at the temperatures of interest ( Tesla).

In this case, flipping a triangular spin and separating the two resulting monopoles leads to an overall increase in the energy of the system unless the monopoles (which live on separate kagome planes) eventually meet other opposite monopoles and annihilate. This process likely incurs a large barrier due to the Coulomb energy that needs to be overcome in order to separate oppositely charged monopoles to large distances.

It would be therefore natural to expect that triangular spin relaxation for Tesla proceeds via stray monopole motion, which faces a rather low barrier to hopping within kagome planes ( K for Tesla). However, contrary to the field range Tesla, where the triangular spins relax whilst plenty of free monopoles are available in the system, we clearly see in Fig. 3 that the triangular spin relaxation for Tesla takes place in a regime where essentially all remaining monopoles are frozen into noncontractible pairs! Indeed, in order to collapse the long time decay of the triangular spin density in the left and middle panel of Fig. S5 we ought to invoke energy barriers of the order of K, which are entirely inconsistent with free monopole hopping barriers K.

In this regime, it appears that separating monopole pairs to large distances is unavoidable, whether the pair originates from the reversal of a triangular spin or from an existing noncontractible pair frozen in the system. Estimating the corresponding energy barrier in this case is a tall order, as it depends on the required separation distance (which is in turn related to the density of reversed spins and noncontractible pairs). Moreover, in order to compute the Coulomb energy cost to separate two charges to said distance, one needs to take into account possible screening effects that are expected in a Coulomb liquid Castelnovo2011 ().

Although this is beyond the scope of the present manuscript, it is intriguing to speculate that this scenario might in fact allow to directly probe the long range nature of the Coulomb interations by measuring the energy scale needed to collapse the long-time triangular spin density decay. Although less straightforward to measure than the overall magnetisation of the sample, it may be possible to access experimentally the magnetisation of the triangular planes alone using neutron scattering measurements in the non-spin-flip channel in the KagIce regime. NMR experiments can also in principle pick out the triangular sites, as they are symmetry distinct from the kagome ones in a field.

Here we limit ourselves to one further observation. If the monopole pair to be separated originates from a triangular spin flip, one expects a negative Zeeman energy contribution, whereas if the monopoles come from a noncontractible pair, the Zeeman contribution is positive. In Fig. S5 (left and middle panel) we observe a larger barrier for Tesla ( K) than for Tesla ( K), which suggests that the former process is the one responsible for the long time decay of the reversed triangular spin density.

## Appendix E Zero-field quenches to spin ice

In this section we discuss the case of quenches down to zero field, driving the system from the SatIce directly to the spin ice regime. This is experimentally relevant because the characterisation of the behaviour of the system following the quench can then be done using measurements in zero field, a scenario that is for instance more suitable to using NMR techniques.

The simulation results in Fig. S6 show that the intermediate-time regime, where the noncontractible pairs are only a small fraction of the total monopole density, is essentially absent in quenches from SatIce to spin ice.

The overall behaviour resembles closely the one following a thermal
quench Castelnovo2010 (). The two protocols share the feature
of being quenches from a monopole-rich phase to a spin ice phase with
low equilibrium monopole density. The only difference is that in one case
the high-monopole-density phase is a fully-ordered ionic crystal of monopoles
while in the other it is a disordered Coulomb liquid. ^{1}^{1}1Note that the
ionic crystal state and the disordered plasma state differ not all that much
in their long range correlations, because long range
Coulomb interactions suppress charge fluctuations in the plasma phase,
and topological kinematic constraints in spin ice strictly forbid volume
charge accumulation Castelnovo2010 ().

Once again, the long time behaviour is controlled by noncontractible pairs. However, since the quench is to spin ice rather than to kagome ice, the triangular spins are not fully polarised in the final state (indeed, there is no longer a distinction between kagome and triangular spins in zero field). As a result, the noncontractible pairs are allowed to decay by direct annihilation of their positive and negative monopoles, say, after they hop around a hexagonal loop on the lattice. This process has been thoroughly described in Ref. Castelnovo2010 () and has a system-size-independent energy barrier due to the cost of separating two Coulomb interacting monopoles up to third-neighbour distance.

We verified that the Poissonian decay model proposed in Ref. Castelnovo2010 () is in qualitative agreement with the long-time behaviour of the monopole density following a field quench to zero field. However, a detailed comparison of thermal and zero field quenches would require an extensive campaign of low temperature MC simulations which we leave for future study.

## Appendix F Finite-rate field ramps

Fast changes in the applied magnetic field are experimentally more manageable than correspondingly fast changes in temperature. For example, they are not limited by the sample heat capacity and thermal conductivity/contact issues. Moreover, spin ice samples are insulators and their inductance is negligible.

However, a sudden change in magnetic field in an experimental setting designed to keep the sample at sub-Kelvin temperatures is nonetheless a tall order. Magnetic field sweeps with a superconducting magnet (as in a typical NMR setting) usually achieve no more than Tesla/min. Opening the circuit in a solenoid would yield a faster field change, but of course this could generate large amounts of heat that the refrigerator is then unable to dispense with quickly enough.

Other techniques that promise to reach larger sweep rates include physically moving the sample relative to the magnet (a permanent magnet can be used given the relatively small fields involved, Tesla); or using magnetic field pulses instead of stepping the magnetic field.

Either way, rates large enough to approximate a quench sudden on the time scale of a single spin flip, ms, will require substantial experimental effort. It is therefore important to run field-sweep simulations and compare what part of the out-of-equilibrium phenomenology discussed in this paper persists for lower sweep rates.

In the left panel of Fig. S7 we compare the monopole density as a function of time from three different protocols: (i) an instantaneous quench from SatIce to Tesla; (ii) a constant rate sweep where the field is lowered from Tesla at time down to Tesla at time MC steps (i.e., approximately s in real time), and it is held constant thereafter; and (iii) a constant rate sweep where the field reaches its final value Tesla at time MC steps (i.e., approximately minute in real time). Similar results arise for other values of the final applied field. It is clear that the long-time regime remains accessible in field sweep measurements. As a matter of fact, the field sweep curves merge with their quench counterpart promptly after the field has reached its final value. Therefore, sufficiently fast field-sweeps can still be used to prepare spin ice in a quasi-static monopole-rich phase (equivalent to a frozen, over-ionised electrolyte), where direct detection of the monopoles might be within reach of zero-field NMR measurements.

In the right panel of Fig. S7 we show the behaviour of the monopole density in field sweeps down to zero field (at the constant rate of Tesla/min, starting from Tesla). Remarkably, even for such slow field sweeps, it is only a matter of reaching a sufficiently low temperature ( K) before the system enters a long-lived metastable state where the density of monopoles remains much larger than in thermal equilibrium.

## References

- (33) Bramwell S T and Gingras M J P (2001) Spin ice state in frustrated magnetic pyrochlore materials. Science 294: 1495-1501.
- (34) Castelnovo C, Moessner R, and Sondhi S L (2012) Spin Ice, Fractionalization, and Topological Order. Annu Rev Condens Matter Phys 3: 35-55.
- (35) Siddarthan R, et al. (1999) Ising Pyrochlore Magnets: Low-Temperature Properties, “Ice Rules”, and beyond. Phys. Rev. Lett. 83: 1854-1857.
- (36) Isakov S V, Moessner R, and Sondhi S L (2005) Why Spin Ice Obeys the Ice Rules. Phys. Rev. Lett. 95: 217201-4.
- (37) Castelnovo C, Moessner R, and Sondhi S L (2008) Magnetic monopoles in spin ice. Nature 451: 42-45.
- (38) Pomaranski D, et al. (2013) Absence of Pauling’s residual entropy in thermally equilibrated DyTiO. Nat Phys 9: 353-356.
- (39) Kimura K, et al. (2013) Quantum fluctuations in spin-ice-like PrZrO. Nat Comm 4: 1934-6.
- (40) Castelnovo C, Moessner R, and Sondhi S L (2010) Thermal quenches in spin ice. Phys. Rev. Lett. 104: 107201-4.
- (41) den Hertog B C and Gingras M J P (2000) Dipolar interactions and origin of spin ice in ising pyrochlore magnets. Phys Rev Lett 84, 3430-3433.
- (42) Castelnvo C, Moessner R, Sondhi S L (2011) Debye-Hückel theory for spin ice at low temperature. Phys. Rev. B 84: 144435-14