Crossover from string to cluster dynamics following a field quench in spin ice

We investigate quench dynamics of spin ice after removal of a strong magnetic field along the [100] crystal direction, using Monte Carlo simulations and theoretical arguments. We show how the early-time relaxation of the magnetization can be understood in terms of nucleation and growth of strings of flipped spins, in agreement with an effective stochastic model that we introduce and solve analytically. We demonstrate a crossover at longer times to a regime dominated by approximately isotropic clusters, which we characterize in terms of their morphology, and present evidence for a percolation transition as a function of magnetization.


I. INTRODUCTION
Studying the nonequilibrium dynamics of many-body systems provides a way to explore phenomena that are not accessible at or near equilibrium.The simplest protocol, at least conceptually, is the "quench" [1][2][3], where a system is initially prepared in equilibrium and then a sudden change is made to one or more external parameters, such as temperature [3] or an applied field [4].Quenches have been realized in numerous experiments, principally in cold atomic gases [5][6][7][8][9], but also in magnetic systems [10].
These include the spin ice materials [11,12], a class of frustrated pyrochlore oxides with unusual low-temperature properties.The frustration results from a combination of the pyrochlore lattice (a network of corner-sharing tetrahedra; see Fig. 1), strong easy-axis anisotropy, and effectively ferromagnetic nearest-neighbor interactions [13,14] (see Sec. II A for details).
In the "classical" spin ices, such as Dy 2 Ti 2 O 7 [15] and Ho 2 Ti 2 O 7 [16,17], quantum effects are negligible [18] and the magnetic moments are well approximated as classical Ising degrees of freedom.(By contrast, in "quantum spin ices" [19] such as Pr 2 Sn 2 O 7 and Pr 2 Zr 2 O 7 significant tunneling occurs and a classical description is insufficient.)The spin configurations with lowest energy are those that obey the "ice rule" on each tetrahedron: two spins point outwards and two point inwards.Such states form an extensive low-energy manifold, resulting in a large residual entropy at least down to temperatures  ≃ 0.35 K [20].In this regime, the system behaves as a strongly correlated paramagnet, referred to as a "Coulomb phase" [21,22], while magnetic ordering is predicted to occur at still lower  ≃ 0.15 K [23].
The elementary excitations above the low-energy manifold are tetrahedra at which the ice rule is broken, where three spins point out and one in, or vice versa.(Tetrahedra where all four spins point out or in also exist and have still higher energy.)Such excitations, which occur at finite density for any nonzero  , are points where the local magnetization has nonzero divergence and are hence monopoles of the magnetic field  [14].These monopoles are deconfined [24], interacting through a magnetic Coulomb law, and can be manipulated by applied magnetic fields [4,25,26].In this work, we use Monte Carlo simulations to study dynamics in classical spin ice following a quench of a magnetic field applied along the [100] crystal direction.In the protocol we consider, a large field is initially applied, polarizing the spins along this direction [27], and then suddenly removed, leaving the spins to relax in zero field.
The equilibrium properties of spin ice in a [100] field are naturally described in terms of strings of flipped spins terminated by a monopole at either end [28,29].At a critical ratio of temperature to field, there is a crossover where such strings proliferate [28], which becomes a ("Kasteleyn" [30]) phase transition in the limit where the proliferation temperature is much smaller than the energy cost of a monopole [28,31,32].
The quench that we consider here effectively drives the system across this transition, starting on the low-temperature side.We show that the dynamics is initially driven by the nucleation and growth of strings and derive an effective stochastic model for these processes, which gives quantitatively accurate results at early times.At later times, we observe a crossover from string to cluster dynamics, which we characterize in terms of the size and shape of the largest cluster.We also find a percolation transition analogous to the equilibrium Kasteleyn transition, which appears to demonstrate a crossing point with system size.
The paper is organised as follows.We begin in Sec.II with a description of the model of spin ice and its dynamics.Following that, in Sec.III, we consider the relaxation of bulk properties, specifically the density of magnetic monopoles and the magnetization, following the field quench.In Sec.IV, we show how the relaxation process can be understood in terms of strings and clusters of flipped spins.The cluster percolation transition is described in Sec.V. Finally, we conclude in Sec.VI, including a discussion of relevance to experiments.

A. Hamiltonian
We study a model of spin ice with classical spins   of magnetic moment  ≃ 10  B on the sites  of a pyrochlore lattice.The spins are constrained to   =   n , where   = ±1 is an effective Ising degree of freedom and the fixed unit vector n points along the local ⟨111⟩ easy axis between the centers of the two tetrahedra to which each site belongs [33].More precisely, we define   = ±1 for the two orientations of tetrahedra  in the pyrochlore lattice, and choose n to point from the tetrahedron with   = +1 to the tetrahedron with   = −1.
The interactions between the spins are well described by the dipolar spin ice Hamiltonian [23] where  is the nearest-neighbor distance in the pyrochlore lattice,  =  0  2 ∕(4 3 ) is the dipole energy scale, and is the interaction energy of a pair of magnetic dipoles.(For simplicity, we neglect further-neighbor exchange interactions, which may become significant at lower temperatures [34,35].)For Dy 2 Ti 2 O 7 , the coefficients take values  = −3.72K and  = 1.41K [23] (we set  B = 1 throughout).Because of the Ising constraint, for every pair of nearest-neighbor sites  and ,   ⋅   = − 1  3     while  dd (   ,   , 3     [23].Taking into account both terms in  DSI , the net interaction between nearest neighbors can therefore be written as −3 eff   ⋅   = + eff     where for Dy 2 Ti 2 O 7 .This antiferromagnetic nearest-neighbor interaction (in terms of   ) is frustrated and is minimized by the "ice rules" states in which two spins point into each tetrahedron and two point out.
Rather than using the full Hamiltonian  DSI , we approximate the dipolar interactions using the dumbbell model [14], which replaces each magnetic dipole, of moment   , by a pair of magnetic charges ±∕ d at the centers of the two tetrahedra to which that spin belongs.These tetrahedron centers form a diamond lattice with nearest-neighbor distance  d = √ 3∕2 (see Fig. 1); in the following, we use  to label both a tetrahedron and the corresponding diamond site.This approximation, which differs from the DSI model by quadrupolar corrections [13], considerably reduces the computational complexity of the problem and provides an accurate approximation except at very low temperatures [36].
Within this model, the Hamiltonian is given (up to an unimportant constant) by [14] where   is the total magnetic charge on diamond site , and the on-site energy [37] can be fixed by considering the energy change due to a single spin flip [14].
The magnetic charge   is given by summing the contributions from the four dumbbells representing the four spins on tetrahedron .We write it as   = 2  ∕ d , where takes integer values 0, ±1, and ±2.(The sum is over sites  belonging to tetrahedron .)The Hamiltonian can then be rewritten as [37] where  C () = || −1 is the (magnetostatic) Coulomb potential and for Dy 2 Ti 2 O 7 .We include the dimensionless parameter , which in reality takes the value  = 1, in order to vary the strength of the dipolar interactions in our simulations.
In practice, to calculate the long-range interactions between the magnetic monopoles, we modify the Coulomb potential  C by including mirror charges, implemented using Ewald summation [38][39][40].
Within the dumbbell model, all configurations that obey the ice rules, i.e., that have   = 0 on every tetrahedron , are degenerate ground states, with energy  gs = 0.The lowestenergy excited states each have a single spin flipped relative to a ground state, or equivalently a pair of charges   = ±1 on adjacent tetrahedra.Their energy is therefore  min = 2 −  C , with the two terms coming from the on-site and Coulomb terms in Eq. ( 7), respectively.Since  min gives the activation energy for dynamics based on single spin flips, we define 2Δ =  min −  gs , giving Rather than treating , the on-site interaction, and , the relative strength of the Coulomb interaction, as independent parameters in our simulations, we choose to fix Δ = 2.8 K, corresponding to Dy 2 Ti 2 O 7 , while allowing  to vary with  according to Eq. (9).Our motivation for this is that we expect the Boltzmann weight for a monopole,  −Δ∕ , to set the principal timescale for the dynamics, and so holding this fixed while varying  allows us to isolate the effects of the long-range interactions.In the dumbbell picture, it effectively means that changing  tunes the strength of the Coulomb interaction between further-neighbor tetrahedra but not between nearest neighbors.
Within the nearest-neighbor model,  = 0, the system can be thought of as a collection of vertices (sites of the diamond lattice) of sixteen different types.Six of these satisfy the ice rule [neutral,   = 0, see Fig. 2(a)], two are all-in or all-out [  = ±2, see Fig. 2(c1,c2)], and the remaining eight are threein-one-out or three-out-one-in [  = ±1, see Fig. 2(b1,b2)].For  > 0, the configuration energy is no longer simply a sum of vertex terms.

B. Dynamics and simulation parameters
Our simulations are performed on pyrochlore lattices with  s = 16  u   u   u sites, where  ,, u are the numbers of cubic unit cells of the fcc Bravais lattice of pyrochlore.Periodic boundary conditions are applied in each direction.The total number of tetrahedra (of both orientations), equal to the number of sites of the diamond lattice, is  d = 1  2  s .We treat the dynamics using the so-called "standard model" [41] of uncorrelated spin flips with a single temperature- independent timescale  flip ∼ 3 ms [42,43].Flips are attempted at randomly chosen sites at a rate  s  −1 flip , and accepted with a Glauber probability [44][45][46] where  is the associated change in energy and  is the temperature.The time  in our numerical results is therefore effectively measured in units of  flip .We consider dynamics following an instantaneous "quench" of the applied magnetic field , which couples to the spins through a Zeeman term  Z = − ⋅ , where is the total magnetization (in units of the ionic magnetic moment ).The initial field is chosen along the [100] crystal direction, which we take as the  axis, taking the value  = (0, 0, ℎ  ), with magnitude ℎ  ≫  .As a result, all spins are aligned with the field to the maximum extent consistent with the easy-axis constraint, as illustrated in the top-left panel of Fig. 1, where [100] points upwards.(We assume that ℎ  is much smaller than the crystal-field term that enforces the easyaxis constraint, which is of order 300 K in the classical spin ice materials [12].For example, an applied field of magnitude  0 || = 1 T corresponds to || ∼ 7 K.)This configuration satisfies the ice rules and so minimizes the (dumbbell-model) Hamiltonian , since each tetrahedron has its top two spins pointing outward and bottom two spins pointing inward [47].
At  = 0, the field is instantaneously reduced to  = , and the subsequent dynamics takes place in zero field starting from the fully magnetized configuration.Since the spin-flip dynamics is ergodic and the (post-quench) Hamiltonian preserves the full symmetry of the lattice, the dynamics progresses towards thermal equilibrium at temperature  , a spin-ice state with a finite density of monopoles and zero mean magnetization.
Our interest in this work is in studying how the thermal equilibrium state is reached and how the nature of the excitations controls the dynamics into the equilibrium state.We first characterize the relaxation to equilibrium by considering bulk properties, before studying how this happens in terms of the appearance of strings and clusters of flipped spins.

A. Density of monopoles
The starting configuration, with all spins pointing upward, obeys the ice rule on every tetrahedron, and so has no monopoles.In the long-time limit, the system reaches thermal equilibrium at temperature  and so monopoles with both |  | = 1 and 2 will occur at finite density.
Our MC results for the density of monopoles (absolute number of monopoles per tetrahedron) of both types,  1 and  2 respectively, are shown in Fig. 3(a,b), for nearest-neighbor spin ice,  = 0, with  s = 128 spins.As expected, both densities increase from zero and quickly reach equilibrium values that increase with temperature.
Starting from an ice-rules configuration, a spin flip produces a pair of monopole excitations of charge ±1 on adjacent tetrahedra.We therefore expect  1 ∼  at very early time, which agrees with our results for  ≲ 0.1.A charge-±2 monopole can be created by a spin flip at a tetrahedron already containing a monopole.Their density should therefore increase as  2 ∼  2 , which is also consistent with the MC results at similar times.
At long times, the equilibrium densities are determined by the activation energy Δ, as well as the further neighbor interactions .For  = 0, the equilibrium densities can be calculated by treating the tetrahedra as independent and considering the Boltzmann weight and multiplicity of each configuration [48] where  1 =  −Δ∕ and  2 =  −4Δ∕ .These values are plotted as horizontal dashed lines in Fig. 3(a,b), and are in good agreement with the MC results.
As a final check that the system is reaching thermal equilibrium and that there is no dependence on the initial state, we compare with simulations starting from a random configuration (effectively  ≫ Δ).In this case, shown in Fig. 3(a,b) with dash-dot lines, the densities of both types of monopoles start at a large value and rapidly decrease, reaching an identical final density as with the field-quench protocol.

B. Magnetization
In Fig. 3(c), we show the time evolution of the magnetization density, which decreases from its maximum value at  = 0 to zero at long times.We define which gives the  component of the magnetization relative to its saturation value is the component of the fixed unit vector n along the  axis.The initial configuration is fully magnetized and so has   = 1, while the full symmetry of the lattice is restored in the equilibrium configuration at long times, and so   = 0.
As with the monopole densities, the relaxation is faster at higher temperatures.The dashed lines in Fig. 3(c) show exponential fits,   =  −∕ , for each temperature, and the fitted relaxation timescale  is plotted as a function of  in the left inset.The observed exponential growth of  with  −1 is consistent with previous studies using the standard model of spin ice dynamics [41,49].The exponential fits match the MC data quite well for  ≳ 3 K, but become increasingly poor as  decreases.This is compatible with relaxation dominated by isolated spin flips at higher temperatures, with collective effects becoming significant only at lower  .We discuss this early-time behavior in more detail in Sec.IV A.
The right inset of Fig. 3(c) shows the same data plotted on a double-logarithmic vertical scale, so that a decrease of the form ∼  −(∕)  , would appear as a straight line with slope −.We see no indication of stretched-exponential decay, which would correspond to  < 1, though the relaxation at the lowest temperatures is broadly compatible with a compressed exponential ( > 1).

C. Finite-size effects and long-range interactions
As a test of the sensitivity of our results to finite-size effects and long-range interactions, we show the dependence of monopole density and magnetization for various system sizes in Fig. 4 and for nonzero  in Fig. 5.In both cases, only small quantitative effects are seen.
The most noticeable effect of increasing  from zero is a decrease in the equilibrium (i.e., long-time) density of monopoles  1 , shown in the inset of Fig. 5(a).This occurs at all three temperatures shown, although for  = 1.22 K the slight decrease for small nonzero  reverses as  approaches 1. (Simulations at  > 1 [not shown] confirm that  eq 1 increases at larger  in all three cases.) The decrease of monopole density with  is, we believe, merely a consequence of the parameterization of the Coulomb interactions that we choose in our simulation: As  increases, we also increase  in order to keep the activation energy Δ fixed [see Eq. ( 9) and the following paragraph].This has the effect of slightly increasing the energy of an oppositely charged pair of monopoles on tetrahedra beyond nearest neighbors, and hence reducing their density.At larger , the longerrange interaction [second term in Eq. ( 7)] can become larger, reversing the effect.
Increasing  from zero also results in slower relaxation of the magnetization, as shown in Fig. 5(b).This effect is more pronounced at lower temperatures, and can be understood as a consequence of the reduced monopole density.

IV. STRING AND CLUSTER DYNAMICS
Up to this point, we have been considering bulk properties, in terms of which the relaxation appears to be quite conventional.We now consider the microscopic processes by which this magnetization occurs, considering the basic processes that allow demagnetization of the fully saturated configuration.
The starting configuration has all spins polarized along the field and hence no monopoles.Starting from the polarized configuration, flipping any spin from up to down produces a pair of monopoles on neighboring tetrahedra, as illustrated in Fig. 6, and so involves a large energy cost 2Δ.Once such a pair has been produced, however, another process becomes possible: by flipping a second spin on either of the tetrahedra, its monopole can be moved to a neighboring tetrahedron; see Fig. 6(c).The only energy cost associated with this second process is due to the change in the Coulomb interaction from separating the two monopoles, which is much smaller than Δ (and zero in the case  = 0).
This process can be continued, separating the two mono- poles along the  direction, and leaving behind a string of downward-pointing spins [28,31], sometimes referred to as a Dirac string [29].Importantly, since these strings are defined with respect to the fully polarized initial configuration, an isolated string cannot form a closed loop unless it spans the periodic boundaries; an open string is always aligned along the  direction and has one monopole at each end. 1  The dynamics following the quench, at least at short times, can therefore be understood in terms of two processes: a slow process whereby a single spin is flipped, producing a pair of monopoles on neighboring tetrahedra; and a fast process where subsequent spin flips separate the monopoles and produce a string of flipped spins.On longer time scales, once a significant fraction of the spins have been flipped downward, the system is no longer well described by isolated strings but rather 1 Given a single spin configuration, there is no way to uniquely define the strings.In the context of the field quench, however, we can define the strings by reference to the initial configuration [29].clusters of flipped spins.

A. Early-time dynamics: Formation of strings
To study the string-cluster crossover within our MC simulations, we identify, for each configuration, all connected sets of flipped (down) spins, where two spins are connected if they are joined by a nearest-neighbor path of down spins.In Fig. 7, we show (, ), the mean number of such sets of each size  ≤ 4 at time  in a system with  s = 1024 spins, for different temperatures  and dipolar interaction strengths .
In all cases, (, ) grows at short times before decreasing, with successively larger  values growing more slowly at first.Comparing with Fig. 4, the maximum is seen to coincide roughly with the time at which the magnetization begins to decrease significantly from its initial value of   = 1.We therefore interpret the initial increase as the regime where the strings are formed and grow independently.The decrease at later times indicates the crossover into cluster dynamics, which we address in Sec.IV B.

Single-string model
In the independent-string regime, it is possible to describe the dynamics analytically by considering the population of strings of each length  (in units of  d ).We assume that the density of down spins is low enough that we can treat each string as well isolated from all others, so that each string effectively grows in a background where all other spins point upwards.In addition, we assume that we can neglect finitesize effects, which is valid if all strings are much smaller than In both main figures and insets, solid lines show the numerical solution of the single-string model defined in Sec.IV A 1, which assumes  = 0 and is valid only at early times.[We solve the coupled equations Eqs. ( 14) and ( 15) with a finite cutoff on , which is varied to confirm convergence.]MC results, with  s = 1024 spins (4 × 4 × 4 cubic unit cells) and averaged over 500 independent runs, are shown with symbols in the main figures, where  > 0, and dashed lines in the inset, where  = 0. Note that mean number is well below 1 at early times for most , indicating that most samples have no sets at all of that size.
the system size in the  direction.We describe the  = 0 case first, and then comment briefly on the effect of long-range interactions.
The number of strings of length  > 1 changes due to growth and contraction of existing strings.In addition, we must consider creation and annihilation processes for strings of length  = 1.To describe the dynamics, we define a stochas-FIG.8. Graph representing the dynamical model for a single string.Vertices (circles) denote strings of length , with  = 0 representing a string that has shrunk to zero length and hence disappeared.Edges (arrows) are transitions with associated rates for growth  + = 2, contraction  − = 1, and annihilation  0 =  G (−2Δ), in units of the spinflip rate  −1 flip .
tic model for a single string, illustrated in Fig. 8: consider a string created at time  0 and let Pr(,  =  −  0 ) be the probability that it has length  ∈ ℤ ≥0 at time .By definition, Pr(, 0) =  ,1 , while the only stationary state is  = 0, an absorbing state representing a string that has shrunk to zero length and disappeared.Growth of a string of length  ≥ 1, producing one of length  + 1, occurs when one flips any of four up spins, two in the row immediately above the top of the string and two in the row immediately below its bottom [see Fig. 6(b)].On the other hand, contraction of a string of length  > 1 requires that one of two down spins, either the top-most or bottom-most of those comprising the string, should be flipped.For each of these processes, there is no change in energy (for  = 0), because a monopole is effectively moved from one tetrahedron to another.Either update is therefore accepted with Glauber probability  G (0) = 1  2 .In units of the overall spin-flip rate  −1 flip , growth therefore occurs with rate  + = 4 G (0) = 2 and contraction with rate  − = 2 G (0) = 1. 2nnihilation of a unit-length string is represented in the single-string model by the absorbing transition from  = 1 to  = 0.This occurs when the single down spin comprising the string is randomly chosen and flipped back upwards.This process reduces the energy by 2Δ, where Δ is the energy of a single monopole.The acceptance probability is therefore  G (−2Δ), and so the rate is simply  0 =  G (−2Δ).
This stochastic model can be solved by expanding in terms of the eigenvectors of the rate matrix, as we describe in Appendix B. For the limiting case Δ∕ → ∞, where  0 →  − , the result can be expressed in terms of modified Bessel functions, while the general result can be expressed in closed form in terms of a contour integral, Eq. (B7).
The single-string stochastic model could be modified to include the effect of long-range interactions,  ≠ 0. In this case, changing the length of the string does involve an energy change, due to the Coulomb attraction between the monopoles at its ends.For  > 2, the energy depends on the path of the string, rather than merely its length, and so one needs to distinguish all possible string shapes, with the number of states increasing exponentially with .The qualitative effect will be to reduce the rate at which short strings grow into longer ones, as seen in the numerical results in Fig. 7.

String populations
A string of length  = 1 is created when any spin is flipped from up to down, as long as all of its neighbors point upward.Since we assume that the density of down spins remains small, we approximate the number of sites where such a string can be created by the total number of up spins, The process costs energy 2Δ and hence the spin flip is accepted with probability  G (2Δ).The rate at which strings are created at time  0 is therefore  G (2Δ) ↑ ( 0 ), within this approximation.The string distribution at time , is then found by considering strings created at all previous times  0 and their probability of reaching size .
The behavior of (, ) to leading-order in  can be determined using Eqs.( 15) and (B13).To this order, we set  ↑ =  s , replacing Eq. ( 14), and so the integral gives simply to leading order in  for each .
Results obtained from the analytical model by numerical integration, with initial condition (, 0) = 0 for all  ≥ 1, are shown with solid lines in Fig. 7. (We apply a cutoff  max on the maximum length  and check that the results are insensitive to the value of  max .)At very short times, they are in good quantitative agreement with those from the MC simulations for  = 0 (with no fitting parameters).
As expected, the results of the analytical model deviate significantly around the time where the maxima are reached and the assumptions cease to apply.As the MC results in Fig. 7 show, the number of short strings decreases rapidly at larger times.
The model of independent strings necessarily fails when the magnetization falls well below its saturated value, and so the number of flipped spins can no longer be treated as small.In the terms of the number of up spins, Eq. ( 14), the magnetization density can be written as For  of order Δ, the rate of string creation is high and significant deviation from saturation magnetization occurs when a large number of strings has been created, even if most strings remain short.As Fig. 7(a) shows, even for fairly low temperature,  = 1.22 K ≃ 0.44Δ significant deviations from the single-string model occur while strings of length  = 1 remain a clear majority.
On the other hand, for  ≪ Δ, creation of strings is suppressed by  G (2Δ) ≈  −2Δ∕ , and string growth is the main process that reduces the magnetization.For intermediate times, where most strings are much longer than the lattice scale but shorter than their separation (and the system size), the single-string model is still valid.According to Eq. (B12), the mean string length is given in this regime by ∑   Pr(, ) ≈ 3 2 + 1 2 , using  G (−2Δ) ≃ 1 for  ≪ Δ. Integrating over time as in Eq. ( 15) and again approximating  ↑ =  s within the integral gives at early time.While this result gives a quantitatively good description of the data only at very short times, it provides a qualitative explanation of the deviation from exponential behavior,   =  −∕ , noted in Sec.III B.

B. Late-time dynamics: From strings to clusters
While the short-time dynamics can be understood in terms of growth of isolated strings, this picture ceases to apply at longer times when the density of strings becomes large.In fact, on the pyrochlore lattice, strings cease to be uniquely defined at higher densities of flipped spins: where all four spins on a single tetrahedron point downwards, there are two equivalent choices of pairings that define different paths for the two strings.We therefore describe the dynamics at later times in terms of clusters of flipped spins; these may be viewed as dense networks of interwoven strings, with the caveat that there is no unique way to "untangle" the strings.
The crossover from strings to clusters is illustrated in Fig. 9(a) and (b), where connected sets of flipped spins are joined by solid blue lines.Inset (a) shows a typical configuration at early time, where a single string has grown along the [100] direction (upwards), while inset (b) shows a late-time configuration containing a large cluster of flipped spins.In the latter case, the cluster includes nearly half of the spins in the lattice and is roughly isotropic, with no privileged orientation.
Using our MC algorithm, we produce a sample of independent runs at  = 0.81 K, and, at each time  in each run, identify all clusters (connected sets of flipped spins) in the configuration.The distribution of cluster volumes  {C} (i.e., the number of clusters containing  {C} flipped spins) is shown in Fig. 9 (main figure).Superimposed on this (red solid line), we show the mean volume  {C} max of the largest cluster, with the average taken over the sample of runs.
We show in Fig. 10(a) and (b) respectively the total number  {C} of connected sets of flipped spins and the mean volume of the largest set  {C} max for several  values and  s = 128 spins.
{C} (V {C} , t) To measure the spatial extent of the largest set, we define where and similarly for  {C} max and  {C} max .In Eq. ( 20), the sum runs over the  {C} d tetrahedra , with centers at   = (  ,   ,   ), that contain at least one flipped spin belonging to the largest connected set.This definition of  {C} max reduces to the root-meansquare radius (rms distance from the centroid) for a small set not spanning the system boundaries, but correctly accounts for the periodic boundary conditions, depending on   only mod- for a cluster that fills the system uniformly.
In Fig. 10(c) and (d) respectively, we plot, as functions of time,  {C} max and the aspect ratio  {C} max =  {C} max ∕ {C} max , a measure of the anisotropy of the largest set.When this forms a long isolated string, such as shown in Fig. 9(a),  > 1, while for a cluster that fills the (isotropic) system, such as Fig. 9(b),  ≃ 1.For the smallest possible set, where a single flipped spin shared by two adjacent tetrahedra,  rms =  rms ≃   ∕2, and so  = 1.
Considering the behavior of these quantities, the dynamics can be divided into three stages.The initial rise in  {C} is due to the proliferation of isolated short strings, as discussed in Sec.IV A. For all temperatures,  {C} max grows linearly with  at the earliest times and the anisotropy measure  {C} max increases rapidly from 1, consistent with linear growth of isolated strings.For small  , this growth of  {C} max continues until it reaches a value of order the linear system size.
The process of string formation is suppressed at low temperatures, since it involves the creation of monopole pairs.The peak in  {C} is therefore lower and occurs at later time for lower  .For all temperatures, this peak is reached at roughly the time when largest set begins to occupy a significant fraction of the system, with  {C} max of order 10% of  s .The subsequent decrease of  {C} can be interpreted as a consolidation process, whereby strings merge into longer strings and clusters.A pair of strings can join if the monopoles of opposite charge on their ends reach the same tetrahedron and hence annihilate.If instead a monopole enters a tetrahedron through which another string passes, the result is a single cluster that can no longer be viewed in terms of isolated strings.
In the third stage, at long times,  {C} reaches a constant value that approaches 1 at low temperature, while  {C} max saturates at 1  2  s .In this regime, a single large cluster fills the system, as illustrated in Fig. 9(b), containing approximately half of the spins.Such a cluster has linear extent  {C} max equal to its maximum value, L, and is approximately isotropic, giving  {C} max = 1, as seen at late times in Fig. 10(c) and (d).
To investigate the consequences of finite-size effects, we show results for different system sizes in Fig. 11.At early time, the number  {C} of connected sets of flipped spins is approximately proportional to system size (number of tetrahedra)  d , consistent with independent formation and growth of strings.For the higher temperature values,  {C} continues to scale with  d up to and somewhat beyond its peak, indicating that consolidation of clusters dominates over formation once they reach a certain ( -dependent, but  d -independent) density.
By contrast, for the lowest temperature [ = 0.81 K; see inset of Fig. 11(a)], the maximum value of  {C} grows more slowly than  d .This suggests that in this case finite-size effects on individual strings are already important at the start of the consolidation process.
At late time, the number of sets  {C} instead approaches 1 for all system sizes, with a single large cluster of down spins filling the system.
Results including long-range interactions are shown in Fig. 12.Comparison with the case  = 0 indicates that there is no significant change in the qualitative behavior, and indeed little quantitative change.Coulomb interactions between monomers are expected to suppress the initial growth of strings, as noted in Sec.IV A, and may also favor reconnection of long strings at the expense of cluster formation.

V. PERCOLATION
Finally, we consider the growth of clusters from the perspective of percolation theory [50,51].Roughly speaking, we consider the set of flipped spins to be percolating if there exists FIG.11.Finite-size effects on the string-cluster crossover: (a) Number  {C} of connected sets of flipped spins versus time for various system sizes (line styles; same dimensions as in Fig. 4) and temperatures (colors).In the main figure,  {C} is scaled by the system size (specifically, the number of tetrahedra  d ).The inset shows the data for  = 0.81 K without scaling, which indicates that  {C} → 1 at long time.(b) Mean volume  {C} max of largest connected set of flipped spins, scaled by the number of pyrochlore lattice sites  s .a cluster that spans the periodic boundaries in the [100] direction.Examples of this criterion are shown in Fig. 13; note that we choose to include only cases where a cluster has nontrivial winding numbers and exclude open strings.(The motivation for this choice is that it corresponds to the criterion for the equilibrium Kasteleyn transition in the limit Δ∕ → ∞ [28], and therefore provides a way to extend this to quench dynamics.) In Fig. 14(a), we show the probability  that a percolating cluster exists (i.e., the fraction of samples that contain such a cluster), as a function of time  for various system sizes  s .The same data are shown in Fig. 14(b) but plotted as a function of magnetization   .(These results are all for the case without long-range interactions,  = 0.) In both cases we see a crossover from  = 0 to 1, which becomes sharper as the system size increases.
In the limit  ∕Δ → ∞, where the spins are independent and flipped randomly, this process is identical to standard bond percolation [52,53] on the diamond lattice, with occupation probability  =  ↓ ∕ s = 1 2 (1 −   ).(Our model has spins on the sites of pyrochlore lattice, which are equivalent to the links of the diamond lattice.)We therefore expect a percolation transition at   = 0.39 [50,54,55], manifesting as a crossing in  as a function of   at (  ) c = 0.22, becoming sharper with , where an open string spans the system, but cannot be assigned a nonzero winding number.
larger system size.
For lower temperatures, we expect strings to percolate at lower  for any given system size, and hence higher   , than independent flipped spins.(As an example, the configuration shown in Fig. 9(a) has a single string, and hence a low density of flipped spins, but is nonetheless percolating.)In our simulation results, we indeed find that decreasing Δ∕ causes each curve of  to shift towards higher   .In fact, we see evidence that the crossing point characteristic of a continuous transition remains for all Δ∕ .

VI. CONCLUSIONS AND OUTLOOK
In this work, we have used MC simulations to study the dynamics of classical spin ice following a magnetic-field quench, where a strong field along the [100] crystal direction is suddenly removed.We have shown how the early-time dynamics can be understood through the formation and growth of strings of flipped spins, terminated at each end by magnetic monopoles, and presented exact results for a simple analytical model of this process.During this stage, the magnetization relaxes rapidly, though with significant deviations from a sim- ple exponential decay for  ≲ 2 K, and the monopole density increases from zero towards its equilibrium value.
At longer times, the system approaches a completely demagnetized equilibrium state.This can be interpreted as a process where the strings merge into clusters that eventually form a network extending throughout the whole system.We have also provided evidence for a percolation transition as a function of magnetization, which extends the equilibrium Kasteleyn transition [28] to the dynamics.Throughout our results, we find that finite system size and dipolar spin-spin interactions (incorporated here using the dumbbell approximation) have relatively modest effects.
This work is partly motivated by experimental realizations of field quenches in spin ice materials [10] that have been performed using Dy 2 Ti 2 O 7 in fields along both the [100] and [111] crystal directions.Our results are in general qualitative agreement with these experiments, which observed relaxation of the magnetization with a timescale increasing rapidly with decreasing temperature [10].(They also observed only minor effects from dipolar interactions, finding good agreement with nearest-neighbor simulations.)One important confounding factor for quantitative comparisons is the question of the extent to which the system remains in thermal equilibrium.Our simulations make the significant simplification of using Glauber dynamics with a fixed temperature, in effect assuming that thermal energy is always transferred rapidly enough for the system to remain locally in thermal equilibrium.In reality, this is certainly not always the case, as evidenced by occurrence of magnetic avalanches induced by local heating [27,56].Incorporating these effects into the theoretical framework developed here is a challenging problem that we defer to future work.
Our simulations use the so-called "standard model" of spin ice dynamics [42,43], where spin flips are attempted at a single fixed rate and accepted with a temperature-dependent probability.Recent work [41,57] has found that the low-temperature relaxation is better described by including two distinct flip rates, reflecting a bimodal distribution of local field strengths.
The consequences of this for the quench dynamics described here can be inferred at early times: For an isolated string, growth and contraction always involves flipping a spin in an asymmetric environment (see Fig. 6), which occurs at the faster rate.By contrast, creation and annihilation of unitlength strings occurs in a symmetric environment and hence at the slower rate.The main result at early times is therefore an effective renormalization of the "seeding" rate to a smaller value.At later times the string density is higher, making the effects more difficult to predict, and simulations incorporating the two rates are required.
Future work will study quenches where the final magnetic field strength is nonzero, including close to the Kasteleyn transition and to the low-temperature side.

FIG. 1 .
FIG. 1. Part of the pyrochlore lattice, a network of corner-sharing tetrahedra, with the [100] crystal direction shown vertically.Topleft: Initial configuration, where all spins are aligned with the magnetic field  = (0, 0, ℎ  ) (gray vertical arrow), as far as possible given the local easy-axis constraints.Bottom-right: Example configuration following the quench of the field to zero.The flipped spins form a "string" of length  = 3, terminated by monopoles (red and blue spheres) at either end.The tetrahedron centers lie on a diamond lattice with nearest-neighbor distance  d = √ 3∕2, where  is the nearestneighbor distance in the pyrochlore lattice.

15 FIG. 3 .
FIG. 3. Time evolution of density of magnetic monopoles with (a)|  | = 1 (b) |  | = 2 at different temperatures  for a system of  s = 128 spins (2 × 2 × 2 cubic unit cells), averaged over 1000 independent runs.Solid (resp.dash-dot) lines corresponds to the fully magnetized (disordered) initial configuration.At long time both densities reach their equilibrium values shown with horizontal dashed lines and obtained from Eq. (12).(c) Time evolution of the  component of the magnetization,   (), for the same  values as in panel (a).Dashed lines show fits to  −∕ for each  .The left inset shows a plot of the fitted  values, along with a fit to  ∝ exp( 0 ∕ ) with  0 = 3.5 K.The right inset shows the same data as the main panel with a double-logarithmic vertical axis, on which a stretched exponential ∼  −(∕)  would appear as a straight line with slope −.

1 FIG. 4 .
FIG.4.Finite-size effects on relaxation of bulk properties: (a) Time evolution of density of (single) magnetic monopoles  1 for various system sizes  s and temperatures  .System dimensions are 2 × 2 × 2 ( s = 128), 2×2×4 ( s = 256) and 4×4×4 ( s = 1024) cubic unit cells.Dashed horizontal lines show the equilibrium monopole density in the thermodynamic limit for each  .At the lowest  values, there is a significant finite-size effect in the equilibrium value of  1 .(b) Time evolution of magnetization   () for the same parameters.

1 FIG. 5 .
FIG. 5. Effects of long-range interactions on relaxation of bulk properties: (a) Time evolution of monopole density  1 for various strengths of long-range interactions, where  = 0 has only nearestneigbor interactions and  = 1 is the dumbbell model (see Sec. II A for details).The temperatures are  = 1.22 K (S1),  = 0.97 K (S2), and  = 0.81 K (S3), and the system size is  s = 128 in all cases.Inset: Long-time values  eq 1 of  1 , evaluated at the time corresponding to the vertical dashed line in the main figure, plotted versus .(b) Time evolution of magnetization   for the same parameters.

FIG. 6 .
FIG.6.Illustration of string formation and growth.(a) In the starting configuration, all spins point upwards, aligned with the external field applied at  < 0. (b) A spin flip produces a pair of monopoles on adjacent tetrahedra, which can be interpreted as a string of unit length,  = 1.(c) Flipping a neighboring spin in the layer either above or (as shown) below the first moves one of the two monopoles, increasing the length of the string to  = 2. From this configuration, the string could be further extended to length  = 3, for example by flipping one of the two bottom-most spins.Note that flipping the top-right spin in the bottom tetrahedron would create a double monopole of charge   = +2, with a much higher energy cost.For this reason, an isolated string always follows the  direction and cannot form a closed loop (unless it spans the system boundaries).

1 FIG. 7 .
FIG.7.Early-time dynamics: Mean number (, ) of connected sets of flipped spins of size  ≤ 4 at time , for various temperatures.In both main figures and insets, solid lines show the numerical solution of the single-string model defined in Sec.IV A 1, which assumes  = 0 and is valid only at early times.[We solve the coupled equations Eqs.(14) and (15) with a finite cutoff on , which is varied to confirm convergence.]MC results, with  s = 1024 spins (4 × 4 × 4 cubic unit cells) and averaged over 500 independent runs, are shown with symbols in the main figures, where  > 0, and dashed lines in the inset, where  = 0. Note that mean number is well below 1 at early times for most , indicating that most samples have no sets at all of that size.

FIG. 9 . 4 FIG. 10 .
FIG. 9. Crossover from string to cluster dynamics following a field quench: Distribution of sizes  {C} of connected sets of flipped spins across a sample of configurations, plotted on a color scale, as a function of time.The system contains  s = 256 spins (2 × 2 × 4 cubic unit cells) and has temperature  = 0.81 K.The red solid line shows the mean size  {C} max of the largest set, which reaches ∼  s ∕2 = 128 at late time.Insets illustrate typical configurations at early and late times, indicated with blue dots on the main figure.(A smaller system, with 2 × 2 × 2 cubic unit cells, is shown for clarity.)(a) Typical early-time configuration.A string of down (flipped) spins that spans the system is highlighted with a blue solid line; all other spins point upward.(b) Typical late-time configuration.Blue solid lines join the centers of neighboring tetrahedra containing flipped spins, which form a single large cluster filling nearly all of the system.

2 FIG. 12 .FIG. 13 .
FIG.12.Effects of long-range interactions on the string-cluster crossover: Same quantities as in Fig.10for various strengths of long-range interactions  (see Fig.5).The top and bottom rows have  = 1.22 K and 0.82 K respectively, and  s = 1024 spins (4 × 4 × 4 cubic unit cells) in both cases.Dashed vertical lines (at the same times in each panel) mark the peak of  {C} at each  .