M Udagawa

L D C Jaubert

C Castelnovo

R Moessner

Out-of-equilibrium dynamics and extended textures of topological defects in spin ice

Keywords: numbers: 75.10.Hk, 75.40.Gb, 75.60.Nt

Memory effects have been observed across a wide range of geometrically frustrated magnetic materials, possibly including Pr2Ir2O7 where the spontaneous Hall effect has been observed. Frustrated magnets are also famous for the emergence of topological defects. Here we explore how the coupling between these defects can be responsible for a rich diversity of out-of-equilibrium dynamics, dominated by topological bottlenecks and multi-scale energy barriers. Our model is an extension of the spin ice model on the pyrochlore lattice, where farther-neighbour spin couplings give rise to a nearest-neighbor interaction between topological defects. This interaction can be chosen to be "unnatural" or not, i.e. attractive or repulsive between defects carrying the same topological charge. After applying a field quench, our model supports for example long-lived magnetization plateaux, and allows for the metastability of a "fragmented" spin liquid, an unconventional phase of matter where long-range order co-exists with a spin liquid. Perhaps most strikingly, the attraction between same-sign charges produces clusters of these defects at equilibrium, whose stability is due to a combination of energy and topological barriers. These clusters may take the form of a "jellyfish" spin texture, centered on a hexagonal ring with branches of arbitrary length. The ring carries a clockwise or counterclockwise circular flow of magnetisation. This emergent toroidal degrees of freedom provides a possibility for time reversal symmetry breaking with possible relevance to the spontaneous Hall effect observed in Pr2Ir2O7.

I. INTRODUCTION

In spin glasses, frustration is an essential ingredient of the glassiness 1 , responsible for nonequilibrium phenomena such as memory effects. But even in disorder-free systems, frustration, now of geometrical origin, imposes considerable constraints on the kinetics, responsible for a diverse range of unconventional dynamical behavior. The interplay between frustration and spin dynamics has been a long-standing source of questions in condensed matter and statistical mechanics [START_REF] Bouchaud | Spin Glasses and Random Fields[END_REF] . For example, the peculiar dynamical response observed in triangular-based organic systems 3 has been intensively discussed recently in terms of geometrical frustration. Indeed, while geometrically frustrated magnets lack disorder to generate a rugged energy landscape, they have topology as a new ingredient for the generation of slow dynamics. In the simplest terms, topological defects cannot spontaneously disappear. Because they carry a topological charge, they can only pair-annihilate with an oppositely charged partner. This provides connections to the physics of reaction-diffusion systems, as well as to the coarsening literature, and more generally to the study of kinetically constrained models [4][5][6] .

While such topological constraints provide hard conditions on the kinds of allowed dynamical processes, they can be supplemented by soft yet important nontopological energetic considerations, such as finite energy barriers. Adding short-range interactions between topological defects can give rise to a rich phenomenology which has not yet been thoroughly or systemati-cally explored-for a case study of quantum Hall physics and the links to different types of superconductivity, see Ref. [7]. Our work aims to add considerations of nonequilibrium dynamics to this intriguing set of phenomena.

A well-studied case in point is provided by the spin-ice model, a canonical example of strong geometrical frustration able to describe the physics of Ho 2 Ti 2 O 7 [8] and Dy 2 Ti 2 O 7 [9] materials. If one restricts the coupling to nearest neighbors, the ground state of the spin ice model has macroscopic degeneracy. Its spatial structure can be described by an emergent Coulomb gauge field arising from a divergence-free condition imposed on the spin texture. This gauge field serves as a vacuum for excited defects carrying a topological charge. Remarkably, dipolar interactions preserve the quasi-degeneracy of this extensively degenerate vacuum 10,11 and dress the topological defects with an effective magnetic charge 12 , interacting together via a magnetic Coulomb potential.

From a dynamical point of view, the divergence-free condition and interaction between monopoles impose strong microscopic kinematical constraints on the motion of spins. Among other consequences, at thermal equilibrium, the magnetic relaxation time diverges at low temperature in Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 spin-ice materials, as measured by AC-magnetic susceptibility [13][14][15] , thermal transport 16 , neutron spin echo 17 , neutron scattering 18 and muon spectroscopy 19 . The slowdown of dynamics is observed irrespective of the experimental probes, over a broad range of time scales. This spin freezing is due to the rarefaction of defects 20,21 , mediated by impurities and surface effects 22 . Indeed, defects play a fundamental role in facilitating spin dynamics 23,24 . In fact, spin ice has been shown to be a fantastic framework to engineer unconventional nonequilibrium physics 25 , whose origin can be multiple. Divergence of the intrinsic relaxation time makes experimental measurements at very low temperature challenging 26 , further complicated by a yet poorly understood spin-lattice coupling which hinders the thermalization of the samples 27 . It is also possible to induce, and even to some extent control, the out-ofequilibrium properties of spin-ice materials via thermal and field quenches where the density of defects can be made larger than at equilibrium 28,29 .

In this context of anomalous spin-ice dynamics, the metallic spin-ice material, Pr 2 Ir 2 O 7 , has recently attracted considerable attention. In this compound, the Ir 5d conduction electrons interact with a magnetic "spinice" texture formed by the localized Pr 3+ moments. The coupling between the spin-ice texture and the itinerant electrons induce anomalous scattering of conduction electrons responsible for unusual transport properties [30][31][32][33][34][35][36][37][38][39] . Among them, the so-called spontaneous Hall effect has been observed in Pr 2 Ir 2 O 7 , and is drawing considerable attention. For a finite Hall signal to be observed, timereversal symmetry of the system must be broken. However, in this compound, a finite Hall response is obtained with neither magnetic field nor spontaneous magnetizationthe usual seeds to break time-reversal symmetry.

This unusual effect is realized in an experimental setup of field quench: the system is initially placed in a magnetic field of 7 Tesla in the [111] direction. After the magnetic field is removed, magnetization relaxes to zero, while the finite Hall signal remains 31 . Possibly, this spontaneous Hall effect may be interpreted as a memory effect due to the special nature of low-energy excitations. The initial state obviously breaks time-reversal symmetry under a magnetic field, and naturally shows a finite Hall response. Accordingly, if there exist non-magnetic excitations which somehow encode the information of timereversal symmetry breaking, then the population of these excitations reflect the broken time-reversal symmetry of the initial state, and the Hall response might persist for a period after the field quench. Moreover, if the excitation has a long lifetime, and do not relax within observable time scales, even the "steady" state may retain a finite Hall response.

Motivated by these questions, we shall study the outof-equilibrium properties of a variant of the spin-ice model, induced by a field quench.

II. MODEL AND SUMMARY OF MAIN RESULTS

We study the J 1 -J 2 -J 3 spin-ice model, which extends the nearest-neighbor spin-ice model, by adding secondneighbor (J 2 ) and third-neighbor (J 3 ) couplings [Fig. 2 (a)]. In this work, we focus on the case J 2 = J 3 (≡ J),

FIG. 1. (color online)

. A spin configuration of the jellyfish structure composed of positive charges. The jellyfish structure has a central ring composed of minimally 6 tetrahedra from which branches of same-sign charges emerge. The loop of spins running along the central ring possesses a chirality whose value (clockwise or counterclockwise) is independent of the sign of the charges. In the branches however, there is an unidirectional flux of magnetization which is imposed by the sign of the charges. In this example, the magnetization flux goes away from the central ring, as can be seen from following the minority spins of each tetrahedron (shown in green).

setting J 1 = 1. This model can elegantly be rewritten in terms of charge degrees of freedom 41 . The spinice degeneracy is preserved and remains in the ground state of the model for an extended region of couplings J ∈ [-1/2 : 1/4]. While the ground state belongs to the vacuum sector of a Coulomb phase, its excitations are described as positive and negative charges, which are called monopoles in the context of dipolar spin ice. Interestingly, an emergent part of the interaction between charges for J = 0 is short-ranged, and its sign is tunable by changing that of J. Quasi-particles with the same topological charge attract (repel) for J > 0(< 0).

Given these basic properties of the model, we first summarize the main results of our analysis. Firstly, let us focus on the region of J > 0. This region is unusual in the sense that "like" charges, i.e. quasi-particles having the same charge, attract, in contrast to, e.g., monopoles in dipolar spin ice. This property immediately leads to the existence of collective excitations composed of like charges, which we pictorially name "jellyfish". A jellyfish excitation is an extended structure consisting of a central ring with an arbitrary numbers of branches attached. Its excitation energy decreases with increasing interaction J, vanishing at the critical point, J = 1/4. In addition to the energetic stability mentioned above, the jellyfish excitation enjoys "kinematic stability", i.e. the clustering of same charges decreases the opportunities of pair-annihilation, which occurs only between opposite charges. The behavior of these charges is the subject of our analysis of their stochastic dynamics.

The dynamical bottleneck caused by the jellyfish excitations leads to the possibility of interesting memory effects, including an emergent chiral degree of freedom which carries zero magnetization. Indeed, the ring part of jellyfish carries a clockwise or counterclockwise circular flow of magnetization, which gives it a well-defined toroidal moment. The stability of jellyfish excitations implies a slow relaxation of these toroidal moments. Accordingly, once the system is subject to a perturbation which breaks time-reversal symmetry, e.g. a magnetic field, this symmetry breaking can in principle persist for a long time, even after relaxation of the magnetization. In the context of a magnetic field quench, this suggests that a signal of broken time-reversal symmetry can be detected a long time after the removal of the magnetic field, e.g. through the Hall response. This provides a possible scenario for the mysterious spontaneous Hall effect observed in Pr 2 Ir 2 O 7 .

As well as in the dynamical properties, the jellyfish excitations leave their fingerprints in thermodynamic quantities, which we reveal by combining Monte Carlo simulations and the analytical Husimi-cactus approximation. Approaching the critical point, J = 1/4, the charge density exhibits a non-monotonic temperature dependence, reflecting the softening of jellyfish excitations. An immediate consequence of the softening is found in the entropy. At J = 1/4, the zero-energy jellyfish add contribution to the ground state degeneracy of the spin ice manifold, enhancing the residual value of the entropy. The crossover behavior at the energy scale of jellyfish excitation gives a clear signature in the magnetic susceptibility, χ zz . Near J = 1/4, χ zz T approaches its high temperature value as the population of jellyfish excitation grows, instead of showing a monotonic Curie-law crossover as in ideal spin ice 42 . Moreover, in addition to these thermodynamic quantities, the jellyfish affects the magnetic structure factor considerably. Instead of pinch point singularities, which serve as an icon of the vacuum of the Coulomb phase, a "half-moon" pattern appears in the magnetic structure factor, reminiscent of the dynamical spin structure factor of the classical kagome antiferromagnet 43,44 . The detection of this pattern through equal-time neutron scattering can serve as a signature of jellyfish structures.

While the region of J > 0 exhibits rich behavior both dynamically and thermodynamically, fertile nonequilibrium behavior is also observed for J < 0, where opposite charges attract as is "normal". For small negative J, the magnetization shows markedly slow relaxation, compared to the charge density. The magnetization takes a constant value over a wide time range. The formation of this magnetization plateau is associ-ated with the exhaustion of charges, and the time scale of its termination can be understood from the "vacuum creation" of pairs of charges.

Approaching J = -1/2, the dynamics changes drastically and becomes dominated by the influence of double charges [Fig. 2]. As a result, charge relaxation becomes extremely slow, while the magnetization decays rather quickly. The system finally forms the co-called fragmented Coulomb spin liquid (FCSL) 45,46 . In the FCSL, charges are long-range ordered but the spin texture remains disordered, extensively degenerate and described by a Coulomb gauge theory. We shall here make the following distinction. A configuration in the Coulomb phase of spin ice is entirely covered with two-in twoout tetrahedra, while a FCSL configuration is alternatively covered by 3-in 1-out and 3-out 1-in tetrahedra: the four nearest-neighbours of a 3-in 1-out tetrahedron are 3-out 1-in tetrahedra, and vice-versa. This phase has recently gained a certain exposure after its experimental observation in several materials: Nd 2 Zr 2 O 7 pyrochlore 47,48 , Dy 3 Mg 2 Sb 3 O 14 kagome 49 , as well as in artificial lattices 50 . In the latter case, the FCSL might also appear in a crystallite form [51][52][53] of an intermediate lowtemperature phase of the kagome dipolar spin ice 54,55 , whose geometry prevents the existence of a charge-free Coulomb phase 56 . Here, the out-of-equilibrium properties of our model offer a working example to stabilize the FCSL phase over long time scales in a spin-ice framework with classical Ising spins on the pyrochlore lattice, a noticeably difficult endeavor at equilibrium 45,46,57,58 .

The nonequilibrium magnetic-field quench proposed here provides a promising tool to realize this state and demonstrates the possibility of engineering a macroscopic state via nonequilibrium techniques 29 .

The organization of this paper is as follows. In section III, we present the general form of the spin model with first, second and third nearest-neighbor interactions. We rewrite this model with charge degrees of freedom, which we name nearest-neighbor dumbbell model, for a specific line of parameters. The corresponding phase diagram is discussed in section IV. The nonequilibrium properties of the model for a field quench are extensively investigated in section V. These results identify a promising point, J = 1/4, on the phase diagram with an enhanced emergent degeneracy. We present a thorough analysis of the equilibrium properties around this point in section VI. Consequences and future directions of research are discussed in the conclusions.

III. THE MODEL

A. J1 -J2 -J3 spin ice model
We shall consider a set of Ising moments at each site i of a pyrochlore lattice, S i = η i d i where d i defines the local easy axis, and η i = ±1 is the Ising variable. The four sublattices of the pyrochlore lattice are defined in Fig. 2, with the following easy-axes:

J1 J3 J2 Q = 0 2-in 2-out Q = -2 3-in 1-out Q = +2 1-in 3-out A D ÒupwardÓ (a) (b) Q = -4 all-in Q = +4 all-out (c)
d {A,B,C,D} = 1 √ 3      1 1 1   ,   1 -1 -1   ,   -1 1 -1   ,   -1 -1 1      (1)
By definition, the vectors d i point outward (resp. inward) the so-called "upward" (resp. "downward") tetrahedra [see Fig. 2.(b)]. We define the Hamiltonian of the J 1 -J 2 -J 3 spin ice model as

H = J1 n.n. S i • S j + J2 2nd. S i • S j + J3 3rd. S i • S j , = J 1 n.n. η i η j + J 2 2nd. η i η j + J 3 3rd. η i η j , (2)
including successively the first, second and third nearestneighbor coupling [see Fig. 2]. Here, we do not consider the coupling right across a hexagon. For J 2 = J 3 = 0, the positive (J 1 > 0, or equivalently, J1 = -J 1 /3 < 0 [59]) nearest-neighbor interaction leads to the spin-ice ground state with the 2-in 2-out Coulomb phase dominating the physics below a temperature of order J 1 . Hereafter, we set J 1 = 1, as a unit of energy and temperature. Meanwhile the second and third terms are expected to considerably alter the qualitative nature of the ground state. However, as we shall discuss in this paper, if J 2 = J 3 ≡ J, the spin-ice ground state is preserved for a broad range of values of J. One of our motivations to consider the J 1 -J 2 -J 3 spinice model comes from the growing interest in pyrochlores beyond their pristine nearest-neighbour models, e.g. in Dy 2 Ti 2 O 7 spin-ice materials [60-63], Tb 2 Ti 2 O 7 [64], Tb 2 Sn 2 O 7 [65], Gd 2 Sn 2 O 7 [66]. Furthermore, the Hamiltonian of Eq. ( 2) is also an effective model for spin degrees of freedom in pyrochlore conductors. In a group of pyrochlore oxides, such as Pr 2 Ir 2 O 7 , the rare-earth moments behave as Ising moments and interact with each other through the RKKY interaction mediated by itinerant electrons originated from, e.g., d orbitals of transition metal ions. Usually, the RKKY interaction is longranged, however, its fast-decaying and oscillating nature sometimes makes it possible to approximate it with a short-range model 67 . While the quantitative description of realistic pyrochlore compounds would certainly require an elaborate model construction, we would like to consider the Hamiltonian of Eq. ( 2) as a minimal model to capture the essential features of these compounds, both in and out of equilibrium.

B. Charge representation: J2

= J3 = J
As mentioned in the previous section, the J 1 -J 2 -J 3 spin-ice model completely preserves the two-in two-out degeneracy, as long as J 2 = J 3 is satisfied. To see this, it is convenient to rewrite the Hamiltonian (2) in terms of charge degrees of freedom, Q p , as was first done by Ishizuka et al. 41 . For a tetrahedron p composed of the sites, p1, p2, p3 and p4, the charge is defined as

Q p = ζ p (η p1 + η p2 + η p3 + η p4 ), (3)
where ζ p = +1(-1), if p corresponds to an upward (downward) tetrahedron. Please note that the sign convention used here is opposite to the one of magnetic monopoles in the dumbbell model 12 . This is an arbitrary choice, possible because while the sign of magnetic monopoles is a physical quantity imposed by the inherent nature of magnetic dipolar interactions, the objects considered here are topological charges whose convention is for us to select. If the configuration of tetrahedron p is two-in two-out, then Q p = 0, while Q p = ±2 or ±4 for single and double charges respectively [see Fig. 2]. Then, if we rewrite Hamiltonian (2) as a function of charge degrees of freedom for J 2 = J 3 , we obtain the nearestneighbour dumbbell model:

H = 1 2 -J p Q 2 p -J p,q Q p Q q . ( 4
)
The details of the derivation are given in appendix A. The two terms in the Hamiltonian (4) allow for simple interpretations: the creation cost and the nearest-neighbor interaction between charges. The two-in two-out degeneracy is trivially preserved (Q p = 0 everywhere), and, as we shall see, even remains the ground state of Hamiltonian (4) for small values of |J|. In other words, in order to lift the degeneracy of spin ice with a short-ranged perturbation, the relevant energy scale is

|J 2 -J 3 |, rather than |J 2 | or |J 3 |.
The main difference with the standard dumbbell model of magnetic monopoles in spin ice 12 is that the interaction between charges is now nearest-neighbour and can be either attractive or repulsive between same-sign charges. Also, our model is a one-parameter problem, J, which means that the creation cost of charges is directly linked to the strength of the interaction.

IV. PHASE DIAGRAM AT EQUILIBRIUM

Before going into the analysis of the dynamics, we shall discuss the equilibrium phase diagram of this model. For J 2 = J 3 = J, the phase diagram is divided into three regions, as shown in Fig. 3. At J = 0, the model is precisely the nearest-neighbor spin ice model, and its ground state is the well-known Coulomb phase where all tetrahedra are in the two-in two-out configurations [see e.g. the pedagogical review by Chris Henley [68]]. Remarkably, the nearest-neighbor and the present spin ice models are projectively equivalent 11 : the extensive degeneracy of the Coulomb phase is preserved for finite J, and it remains ground state of the model for a wide region of parameters, -1/2 < J < 1/4. 4) as a function of J. This corresponds to the J1 -J2 -J3 spin-ice model for J2 = J3 = J, and J1 = 1. The Coulomb phase is stabilized over a broad region (-1/2 < J < 1/4). For J < -1/2, it gives way to the antiferromagnetic phase made of alternating all-in and all-out tetrahedra (double charges). At J = 1/4, the jellyfish point discussed in section VI appears. For J > 1/4, the ground state ensemble is sub-extensively degenerate; spin configurations are made of one methane-like cluster per cubic unit cell [see Fig. To understand the robustness of the Coulomb phase, we rewrite the charge Hamiltonian given in equation ( 4) as

H = 1 2 + J p Q 2 p - J 2 p,q (Q p + Q q ) 2 . (5)
Since both sums are non-negative, one can immediately conclude that for -1/2 < J ≤ 0 the ground-state configurations satisfy Q p = 0 for all tetrahedra p, which make both terms exactly zero. This condition ensures two-in two-out condition for all the tetrahedra, i.e., the system is in the Coulomb phase in this region. Meanwhile, for J < -1/2, the first term of Eq. ( 5) becomes negative and is now minimized if |Q p | is maximized, i.e., Q p = ±4 for all p. At the same time, the second term is minimized if Q p + Q q = 0. These two conditions are simultaneously achieved, if and only if Q p = +4(-4) for all upward (resp. downward) tetrahedra. In terms of spins, this corresponds to all-in / allout configurations, which can be seen as a double charge crystal (zinc-blende structure).

Right at J = -1/2, the creation cost of charges in Hamiltonian (5) vanishes, making the double charge crystal degenerate with the Coulomb phase. Moreover, charge crystals made of alternating charges with Q p = +2 and -2 on upward and downward tetrahedra respectively also minimize the energy [see an example of spin configuration in Fig. 14]. This corresponds to the fragmented Coulomb spin liquid 45 where long-range charge order co-exists with extensive spin degeneracy.

As for the other side of the phase diagram, J > 0, the limit of stability of the Coulomb phase is more subtle. If one rewrites the Hamiltonian(4) as

H = 1 2 -3J p Q 2 p + J 2 p,q (Q p -Q q ) 2 , (6)
it guarantees that the ground state satisfies Q p = 0 for all tetrahedra p, at least up to J = 1/6. The interaction term [see Eq. ( 4)] clearly favors same-sign nearest neighbours, a condition that causes topological problems. First of all, charge neutrality must be preserved over the entire system: there must be as many positive as negative charges, which means that not all bonds can satisfy the condition (Q p -Q q ) = 0 for nearest neighbours p and q. Furthermore, when one tries to cover the pyrochlore lattice with only, say, positive charges, one quickly encounters steric problems imposing the inclusion of vacuum or negative charges, even at a short length scale. According to equation (4), it is relatively easy to see that there are only two kinds of spin clusters that minimize the interaction energy [Fig. 17]:

• a double charge surrounded by four same-sign single charges,

• a closed ring of same-sign single charges (this will give rise to the jellyfish structure discussed in detail in section VI).

In both cases, one gets locally for the tetrahedra which are part of these minimal clusters p∈cluster

Q 2 p = p,q ∈cluster Q p Q q . (7)
Furthermore, it is possible to attach same-sign single charges to these minimal clusters, but no more double charges. Indeed if one takes the example of negative charges, then all spins on the outskirt of the cluster point outwards, which strictly forbids to attach any more double-negative charges. The addition of samesign single charges to these minimal clusters also preserves the equality of equation (7). But as we have just discussed, not all p, q bonds can satisfy the condition (Q p -Q q ) = 0. It means that these clusters cannot cover the entire system. When summed over all tetrahedra, one gets

p Q 2 p ≥ p,q Q p Q q . ( 8
)
According to Hamiltonian (4), it means that the total energy E of any spin configuration for J > 0 has a strict lower bound:

E ≥ 1 2 -2J p Q 2 p . (9)
For 0 < J < 1/4, the two-in two-out Coulomb phase thus remains the ground state of the system. Since we have shown the existence of configurations where the inequality of equation ( 9) is a strict equality, the Coulomb phase is not stable anymore for J > 1/4. As for the physics at J = 1/4, this deserves a dedicated discussion in section VItxtr, where we also briefly discuss the nature of ground state for J > 1/4.

V. NONEQUILIBRIUM DYNAMICS VIA FIELD QUENCHES

In this section, we investigate the stochastic dynamics of the model (2). To characterize the system, we focus on the single charge density ρ 1 , and magnetization along [111] axis. ρ 1 is defined as the sum of positive and negative single charges, divided by the number of tetrahedra, so that the maximal value of ρ 1 is normalized as 1. M is also normalized to be 1, if the system is in a saturation field in [111] direction.

A. Waiting time Monte Carlo method

Since the J 1 -J 2 -J 3 spin ice model is composed of classical degrees of freedom, the rule of time evolution is not determined a priori. Here, we assume a stochastic dynamics defined as follows. Suppose that the spin configuration Ω α ≡ {S i } is realized at time t, with probability P (Ω α ). P (Ω α ) evolves with time by the following stochastic equation: (10) where τ 0 gives the unit of time, and W (Ω α → Ω β ) is the transition rate from the state Ω α to Ω β . Here, we assume only a single spin flip process, i.e., W (Ω α → Ω β ) is finite, if and only if the state Ω β can be obtained from Ω α by flipping a single spin. For the two configurations, Ω α and Ω β , we assume a transition rate of thermal bath type,

d dt P (Ω α ) = 1 τ 0 α =β [P (Ω β )W (Ω β → Ω α ) -P (Ω α )W (Ω α → Ω β )],
W (Ω α → Ω β ) = exp(-βE(Ω β )) exp(-βE(Ω α )) + exp(-βE(Ω β )) . (11)
In order to solve the stochastic equation ( 10), we resort to the waiting time Monte Carlo method. The procedure of this numerics is divided into several steps. Suppose, the system is at n-th step with time, t n , then the calculation goes as follows: (i) Make a table of transition probabilities: p i for all possible Ising variable, η i , to flip, according to the transition rate (11). (ii) determine which η i to flip, according to the probability p i /p tot (p tot = i p i ), by generating a random number Firstly, we review the results for nearest-neighbor spin ice: J 2 = J 3 = 0 in Fig. 4. Here, the temperature is initially set at T = 10, then suddenly quenched to T = 0 at time t = 0. J 1 = 1 as reference energy. In this setting at t = 0, an initial charge density ρ i is thermally excited: ρ i = 0.487. In Fig. 4, we plot the time dependence of the charge density. After the temperature quench, charges start to annihilate in pairs, and their density decreases monotonically. Their density time dependence can be described by a simple reaction-type equation 28 .

r 1 ∈ [0, 1]. (iii) obtain the time τ taken to make the flip by τ = -1 ptot log(r 2 ) (r 2 ∈ [0, 1]). (iv) update the time t n+1 = t n + τ . (v) repeat (i)-(iv)

Numerical setup for field quenches

But temperature is not the only variable that can be used for a quench. In an anisotropic system such as spin ice, field quenches offer an alternative route for nonequilibrium phenomena, which are at the core of this paper.

ρ1(t) = ρi 1 + Kρit , where K = 9 √ 3 16
in agreement with

Ref. [28]. The initial density is ρi = 0.487.

For this purpose, we supplement the Hamiltonian (2) with the Zeeman term,

H Z = -H • i S i . (12)
In a field quench, the system is initially set under the magnetic field H [111], i.e. parallel to one of the easy axes, with |H| = 100, which can be practically regarded as a saturation field. Then, the field is removed suddenly at time t = 0. At the initial stage, all the spins are aligned in the [111] direction, and all the tetrahedra are occupied with charges accordingly , i.e. the system is in an "antiferromagnetic" single charge crystal [see configuration in Fig. 5].

As shown in Fig. 6, the charge density decreases monotonically after the field quench, within a time scale of O(1). The decrease of charge density is dominated by a simple pair-annihilation process of charges. This process accompanies the reduction of energy 4J 1 , and takes place within the same order of time scale. From now on, we shall consider the stochastic dynamics of Hamiltonian (4) (J 2 = J 3 ≡ J) after a field quench of a saturated [111] field: from very large |H| = 100 to H = 0. As presented previously, it is important that the initial configuration at t = 0 is made of alternating positive and negative charges on the diamond lattice formed by the centers of the tetrahedra [see configuration in Fig. 5]. Thus not only is the magnetization fully saturated in the [111] direction; the charge degrees-of-freedom are also perfectly aligned. The multiple energy barriers associated with the most relevant dynamical processes are summarized in appendix D.

1. -0.20 < J < 0.00

We start with the negative region, J < 0. The simple comparison of Figs. 6 and 7 for J = 0 and J = -0.1 respectively shows that the introduction of J qualitatively alters the dynamical properties of the system. In Fig. 7, the total charge density and magnetization of the system are plotted together with the sublattice magnetization decomposed into triangular and kagome planes [see Fig. 5 for definitions]. The charge density decays monotonically and quickly reaches practically zero on a time scale of O(1), in a similar way to the case of J = 0. But in contrast to non-interacting spin ice, the magnetization shows a nontrivial time dependence with a wide plateau region. The dynamical process can be qualitatively divided into four time domains.

The time evolution starts with a nucleation process. Since the pyrochlore lattice is decomposed into triangular and kagome planes [Fig. 5], the initial event should either be a spin flip on a triangular plane or a kagome plane. The former process is unlikely to occur, since the spin flip on triangular plane gives rise to double charges Q = ±4 whose energy cost is ∆ 5 = 12 -12|J| > 0 [see appendix D for details]. In contrast, the spin flip on the kagome plane lowers the energy, ∆ 1 = -4 + 20|J| < 0, which is why this process occurs within the time scale of order O (1).

Consequently in the time domain I of Fig. 7, spins flip almost randomly but only on kagome layers, creating pairs of vacuum tetrahedra -two-in two-out states. This is confirmed by Fig. 7 where only the magnetization on kagome planes decreases, while that on triangular planes remains constant. However, these spin flips are not completely random. Indeed, one could have thought that the created pairs of vacuum tetrahedra can dissociate and diffuse within their kagome planes. However, vacuum tetrahedra are confined for J < 0 because their diffusion brings same-sign charges next to each other, a process whose energy cost is roughly proportional to the distance between the diffused pair of vacuum tetrahedra [see Fig. 8]. This process is thus prohibited at the beginning of the relaxation (time domain I), as opposed to the J = 0 case.

The spin flips on kagome layers strongly suppress the density of charges in the system [see Fig. 7]. Once the charges become dilute enough, they can move freely, without contact with other charges. This is the beginning of the second time domain where, in particular, charges can hop between adjacent kagome layers by flipping spins on triangular planes, resulting in the decrease of the triangular sublattice magnetization.

However, in addition to the diffusion of charges, their density continues to decrease, which eventually causes an exhaustion problem. The system may not have had the time to fully decorrelate before the exhaustion of charges, leaving a long time plateau with finite magnetization and no charges. This is the third time domain of Fig. 7. Please note that this time plateau may split into two parts for very small temperatures [see discussion in appendix E].

The system then remains frozen until the creation and dissociation of charge pairs spontaneously occur, on a time scale of the order of exp[(∆ 3 + ∆ 4 )/T ], where ∆ 3 = 4 + 4|J| and ∆ 4 = 4|J|. The dissociated charges propagate through the system and bring the magnetization to its equilibrium vanishingly small value, defining the fourth and final time domain. In Fig. 9, the time t p at which the charge density finally decays to zero is plotted as a function of 1/T . t p clearly follows the predicted scaling law (13) This confirms that the lifetime of the magnetization time plateau is determined by the pair creation and dissociation process of charges.

t p ∝ exp[(∆ 3 + ∆ 4 )/T ] = exp[(4 + 8|J|)/T ].
To summarize so far, the nonequilibrium relaxation dynamics for -0.2 < J < 0.0 displays the rich multiscale behaviour, well understood from microscopic processes. However, in light of the possibility of spontaneous Hall effect where we are looking for broken time-reversal symmetry without magnetization, this region is not particularly relevant since magnetization persists on a much longer time scale than the charge density.

= >

FIG. 8. (color online). Illustration of the relaxation dynamics for -0.2 < J < 0.0. (a) The initial process that occurs in the saturated configuration is the spin flip on the kagome plane, which leads to the pair annihilation of charges. While the diffusion of vacuum tetrahedra is energetically unfavorable, successive pair annihilations are possible, making the charge density more and more dilute. After enough dilution, the charges start to spread through the system, as highlighted by red circles. (b) Before the magnetization relaxes to zero, charges are exhausted. The relaxation of magnetization is then only possible via pair-creation and dissociation of charges. Spin flips are highlighted by yellow circles when necessary.

2. -0.25 < J < -0.20

The nature of dynamics qualitatively changes at J = -0.20. Both magnetization and single charge density remain saturated over a time scale growing exponentially with 1/T after the magnetic field is switched off at t = 0 [Fig. 10]. Their decay is then simultaneous and visible in both triangular and kagome layers. After this decay, a wide plateau follows and the dynamics becomes similar to the region for smaller |J|.

The persistence of the initial state can be attributed to the dynamical bottleneck of the initial process. As discussed in the section V C 1, the initial process is given by the spin flip on the kagome plane, which costs ∆ 1 = -4 + 20|J|. Consequently, for J < -0.20, this process requires a finite energy and takes an exponentially long time at low temperatures.

However, once a pair of charges is created, the next process becomes immediately available, since it costs no more energy to create a pair of vacuum tetrahedra adjacent to the nucleation pair: ∆ 2 = -4 + 16|J| < 0. Consequently the first spin flip serves as a nucleation seed for an avalanche effect of charge annihilation, until the charges completely disappear from the system. After the exhaustion of charges, the system reaches a plateau with residual magnetization before finally relaxing thanks to the creation and dissociation of charges, as observed for -0.2 < J < 0.0. 3. -0.50 < J < -0.25 a. Dynamics In this parameter range, both the nucleation (∆ 1 = -4 + 20|J| > 0) and the proliferation (∆ 2 = -4 + 16|J| > 0) processes cost a finite energy, preventing the avalanche effect that takes place for -0.25 < J. As an alternative mechanism of relaxation, one needs to consider the dynamics along the [111] direction, perpendicular to the kagome planes. The first two steps correspond to the creation of a vacuum tetrahedron and a double charge [see Fig . 11]. The energy of this two-step process is ∆ 1 + ∆ 9 = ∆ 5 + ∆ 6 = 4 + 8|J| and is independent of whether the nucleation occurs in a kagome or a triangular plane [see the appendix D for details]. When J -3/7 ≈ -0.43, this two-step mechanism along the [111] direction remains more expensive in energy than the two-step proliferation in kagome planes previously considered, ∆ 1 + ∆ 2 4 + 8|J|. But the dynamics after these two initial steps is mediated by coexisting energy barriers with similar magnitudes, and the system is subject to a rather complex dynamical competition for -3/7 < J < -0.25 which varies from short to long length scales.

On the other hand for -0.5 < J < -3/7, the mechanism along the [111] direction becomes the most favorable source of relaxation, both for their two-step nucleation

(∆ 1 + ∆ 9 = ∆ 5 + ∆ 6 < ∆ 1 + ∆ 2 ) and their propagation (-∆ 7 = ∆ 8 = 8 -16|J| < ∆ 2 ) [see Fig. 11].
This scenario is confirmed by simulations in Figs. 12 and 13. Most appreciably, the total magnetization drops to zero, while i) the charge density keeps its saturated value for several orders of magnitude in time and ii) the sublattice magnetizations on triangular and kagome layers both take finite values. The extinction time of the total magnetization is in quantitative agreement with an Arrhenius law exp[-(∆ 5 + ∆ 6 )/T ] = exp[-(4 + 8|J|)/T ] at low temperature, reflecting the magnitude of the energy barrier of the nucleation process.

As for the plateau of Fig. 12 after the extinction of the total magnetization, it comes from the fact that, as opposed to the proliferation of vacuum tetrahedra observed for -0.25 J, the mechanism along the [111] direction is mediated by the propagation of pairs of local excitations: vacuum tetrahedra and double charges 46 . Both excitations can be broken into two single charges: one which is fixed and respects the initial charge order, remnant of the field quench, and the second one which can move along the [111] direction. As a result, the density of single charges varies only very weakly, in order to accommodate a vanishingly small density of local excitations. For the parameter region -0.5 < J < -3/7, the propagation along the [111] direction costs zero or little energy -∆ 7 = ∆ 8 = 8 -16|J| ∈ [0; 1.14] when compared to the nucleation process ∆ 1 + ∆ 9 = ∆ 5 + ∆ 6 = 4 + 8|J| ∈ [7. 43; 8]. As a consequence the system does not relax into the vacuum ground state made of two-in two-out tetrahedra [see Fig. 3], but conserves the initial long-range charge order illustrated in Fig. 14.

b. Fragmented Coulomb spin liquid

The FCSL phase is partially ordered due to a broken Z 2 symmetry of the charge degrees-of-freedom: every positive charge has four negative nearest-neighbour charges, and viceversa [see Fig. 14]. This leaves a Z 2 degeneracy of the charge order, as for the all-in / all-out phase. But as opposed to the all-in / all-out order, the magnetic "crys- tal" of the FCSL phase is made of single charges, which retains an extensive spin degeneracy. It also supports an emergent Coulomb gauge field due to a local divergence- free constraint on the magnetization flux 45 . When averaged over its Gibbs ensemble, the FCSL phase does not bear any finite total magnetization. However, the broken Z 2 symmetry of the charge order gives rise to finite sublattice magnetizations along the [111] direction.

FIG. 14. (color online

). An example configuration of the fragmented Coulomb spin liquid (FCSL) 45,46 . All upward (downward) tetrahedra are occupied by positive (negative) charges.

Please note that this broken symmetry is not spontaneous, but selected by the initial [111] magnetic field.

As illustrated on Fig. 5, a spin on a triangular layer always has a negative charge above, and a positive one below, resulting in an averaged triangular magnetization m triangular = -1/4. On the other hand, a spin on a kagome layer always has a positive charge above, and a negative one below, resulting in an averaged kagome magnetization m kagome = +1/4. These values correspond to the plateaux observed in Fig. 12 for t > 10 4 , and imply that a substantial part of the FCSL Gibbs ensemble is visited during the relaxation.

In the present nearest-neighbour dumbbell model, the FCSL phase is part of the ground-state ensemble only for J = -0.5, degenerate with the all-in/all-out order and the two-in two-out spin-ice ground state [see Fig. 3]. It disappears however at finite temperature in favor of the two-in two-out spin-ice regime because of their entropy difference. The present setup thus offers a rare mechanism able to stabilize a metastable form of the FCSL phase in an Ising pyrochlore model, expanding the versatility of nonequilibrium physics in spin ice 25 . This out-of-equilibrium scenario contributes to the current effort to observe fragmentation in frustrated magnetism [47][48][49][50] .

c. Hall response

The final question regarding the FCSL phase is its relevance to the Hall effect. For this parameter range, the relaxation of the charge density is much slower than that of the magnetization [see Fig. 12]. In this sense, the situation may be similar to the experimental condition of the spontaneous Hall effect observed in Pr 2 Ir 2 O 7 [32]. However, based on symmetry arguments, the Hall conductivity is exactly zero in the FCSL phase. Quite generally, provided the spin configurations on upward and downward tetrahedra can be transformed into each other through time-reversal, then the Hall conductivity becomes zero, if spatial homogene-ity is preserved on average [see appendix F]. This means that a finite Hall response is not expected to occur for negative J, where positive charges tend to sit next to negative ones. Now, let us consider the other side of the phase diagram, with J > 0. Charges of opposite sign repel, making the initial [111] polarized state particularly unstable. The dynamical feature after the field quench is plotted in Fig. 15, where the decrease in charge density ρ 1 and total magnetization M start very quickly. After a time scale of O(1), both ρ 1 and M develop a shoulder-like feature after which the relaxation process slows down. Nevertheless, the magnetization decays smoothly to zero, and subsequently, charges also disappear completely. Interestingly, the whole time dependence of ρ 1 and M can be scaled with a single parameter J/T [see Fig. 15].

The early stage of the relaxation is very similar to the -0.20 J 0 case of section V C 1. There is no energy barrier preventing the relaxation from the initial [111] polarized state. Vacuum tetrahedra quickly proliferate in kagome layers, diminishing the kagome magnetization in the process. Once the density of charges becomes small enough, spins in the triangular layers can finally relax without creating pairs of double charges. However, as opposed to the J < 0 scenario, nearest-neighbour pairs of positive/negative single charges can separate at no energy cost since they repel. There is thus no energy barrier for the diffusion of charges. The resulting separation of charges makes their eventual annihilation statistically more scarce and energetically unfavorable since a pair of positive/negative charges needs to get close to each other before being able to annihilate. Accordingly, in the absence of high-energy charge creation processes, the dynamics in this time domain is dominated by a single energy scale, the nearest-neighbour repulsion, proportional to J [see equation ( 4)], which explains the J/T scaling and the quasi-plateau of charge density observed in Fig. 15. Furthermore, even if charge diffusion first slows down the magnetization decay (the diffusion process is mostly uncorrelated with the direction of the spin), it ultimately allows for the complete relaxation of the system down to zero magnetization.

For larger values of J, the magnetization shows approximately the same curve as in Fig. 15, but the previous J/T scaling breaks down for the charge density ρ on long time scales [see Fig. 16]. It implies the appearance of a new energy barrier which is not simply proportional to J, i.e. a dynamical mechanism involving charge creation/annihilation [see Eq. ( 4)]. This is understood by the fact that the creation cost of charges in equation ( 4) diminishes as J increases. The creation of a pair of charges can be further facilitated if it occurs in the appropriate "cage" of same-sign charges. This is why, for large enough J, the energy barrier of the nearestneighbour repulsion which was responsible for the J/T scaling can be replaced by a charge-creation barrier. Finally, when J reaches the value of 1/4, the relaxation of ρ 1 becomes non-monotonic, underlining a qualitative change of physics at J = 1/4 as discussed in section IV. At J = 1/4, the nearest-neighbour dumbbell model becomes

H = 1 4 p Q 2 p - 1 4 p,q Q p Q q . (14)
Hence, the energy cost for creating a charge is exactly balanced with the energy gain from placing a pair of same-sign charges next to each other. As briefly introduced in section IV, there are two minimal clusters of spins respecting this balance. The first one is a closed ring of charges, where the number of charges is equal to the number of bonds between them. The ring is made of at least six tetrahedra. The second type of cluster is reminiscent of a methane molecule, composed of a double charge next to four single charges [see Fig. 17]. These two clusters are actually the seeds of extended zero-energy textures. Attaching same-sign charges to these seeds does not cost energy since it adds the same number of charges and bonds. As illustrated in Fig. 1, it is possible to attach branches to the central ring, which we pictorially refer to as a jellyfish structure. Consequently, for J = 1/4, the two-in two-out spin-ice ground state receives additional macroscopic degeneracy coming from these composite forms. The residual entropy increases to S/k B ≈ 0.265, compared with that of nearestneighbor spin ice ≈ 0.2058 [Fig. 18(b)]. Away from the J = 1/4 point, these structures become excitations with energy:

E jellyfish = 4(N branch + 6) 1 2 -2J , (15)
E methane = 4(N branch + 8) 1 2 -2J , (16)
where N branch is the number of charges in the branch.

If the charge description is elegant, one should not forget the underlying spin texture. As illustrated in Fig. 1, each tetrahedron in a branch hosts a charge, and therefore has three majority spins and one minority spin.

Each tetrahedron is connected to its "parent" tetrahedron strictly and uniquely via the minority spin, which means that the branches can bifurcate or trifurcate, but never form closed loops nor connect to other jellyfish or methane-like structures. Such connection is energetically unfavorable for opposite-sign structures and topologically forbidden for same-sign structures. Moreover, the sequence of spins along the backbone of any branch form an oriented path. Depending on the sign of the charges, this oriented path flows either towards or away from the central seed (the ring or the double charge).

For the jellyfish structure, there is an additional degree-of-freedom present in the central ring, which carries a circular magnetization flow with two choices of chirality: clockwise or counterclockwise. Due to this circular magnetization flow, the ring acquires a toroidal moment, T ring , defined by

T ring = i r i × S i , (17)
where, r i is the relative coordinate of site i, measured from the center of the ring. Consequently, T flips its sign when the magnetization flow is reversed. It is an important property which allows the local encoding of time-reversal symmetry in a zero-energy texture, independently of the sign of the charges it belongs to, i.e. independently of global time-reversal symmetry. This emergent toroidal moment takes an enhanced flavor in the context of field-quench dynamics discussed in the previous section. If a coupling between spin and spatial degrees of freedom is taken into account -e.g. by adding spin-orbit coupling in our modelthe field induced time-reversal symmetry breaking of the initial state can be potentially encoded in the ring of the jellyfish structures. In other words, toroidal moments may persist long after the magnetic field is quenched. It is tempting to associate this possibility with the spontaneous Hall effect observed in Pr 2 Ir 2 O 7 [32]. Namely the residual Hall response may be ascribed to the combined effects of slow relaxation and toroidal ordering due to the jellyfish clusters. Indeed, the decay of a jellyfish cluster incurs an activation energy barrier of the order J and becomes exponentially slow in J/T at low temperatures (see Appendix G for a discussion of the possible decay processes and their relative barriers).

B. Crossover into the jellyfish regime

We have studied the emergence of the jellyfish structures by the means of Monte Carlo simulations and Husimi-cactus approximation at equilibrium. Technical details of the methods are discussed in appendices B and C.

Let us first define the variable (color online). Absence of long-range order for J = 1/4, probed by Monte Carlo simulations (black circles) and Husimi-cactus calculations (red dashed lines). (a) Magnetic susceptibility along the z-axis times temperature. We observe a Curie law crossover without any singularities (i.e. without any phase transition), characteristic of a classical spin liquid at low temperature 42 . (b) Entropy of the jellyfish model compared to the nearest-neighbor spin-ice model (yellow crosses). The black dashed lines show the lowtemperature limits from Monte Carlo and confirm the higher entropy of the jellyfish ground state. (c) Specific heat. (d) Quasi-order parameter defined in eq. 18, showing short-range ordering with a finite fraction of same-sign nearest-neighbour monopoles. (e) Density of tetrahedra configurations. All observables are normalized per number of spins, except for panel (e) which is normalized per number of tetrahedra. The temperature axis is on a logarithmic scale for panels (a,b).

M Q.Q = 1 N p,q Q p .Q q 4 . (18)
where N is the total number of pyrochlore sites. M Q.Q is a spatially averaged measurement of the nearestneighbour ordering of charges. In the all-in/all-out antiferromagnetic order and in the FCSL phase, it takes the values M Q.Q = -4 and -1 respectively. Please note that M Q.Q is not a proper order parameter of these phases, since M Q.Q does not differentiate between the configurations with broken Z 2 symmetry. But it is a convenient probe of the evolution of nearest-neighbour correlation and it directly couples to the interaction term of equation (14).

M Q.Q takes a finite positive value as expected for J = 1/4 > 0, and reaches a plateau for T < 0.2 [see Fig . 18.(d)]. This plateau is a signature of the jellyfish regime and it is also revealed by a corresponding plateau in the charge density ρ in Fig. 18.(e), as confirmed by the quantitative agreement between Monte Carlo simulations and Husimi-cactus calculations [see the dashed lines in Fig. 18 for the Husimi-cactus results].

The jellyfish regime is characterized by the following observables as T → 0 + :

ρ 0 = 0.670 ± 0.005, ( 19
) ρ 1 = 0.324 ± 0.005, ( 20
) ρ 2 = 0.005 ± 0.001, ( 21
) M Q.Q = 0.171 ± 0.005, ( 22
)
χ zz T = 0.333 ± 0.003. ( 23
)
Let us check if these values are consistent with our understanding of section VI A. Let N and N t be the number of spins and tetrahedra respectively, with N = 2N t . Then there are N t ρ 2 double charges, each of them surrounded by four single charges [see section VI A]. Thus 4N t ρ 2 single charges are directly linked to a double charge, while N t (ρ 1 -4ρ 2 ) are linked to a single charge in the branches of jellyfish or methane structures. One can thus estimate

M est Q.Q = N t [(ρ 1 -4ρ 2 ) + 2 (4ρ 2 )] N = ρ 1 + 4ρ 2 2 ≈ 0.172.
which quantitatively matches the measured value of equation (22). Please note that the above formula is also valid for Husimi calculations, but with somewhat lower values of M Q.Q = 0.165, ρ 1 = 0.316 and ρ 2 = 0.00363. For J = 1/4, the ground-state ensemble does not break any symmetry of the Hamiltonian, and includes all of the two-in two-out configurations (traditional spinice ground states) and all of the configurations with the structures discussed in section VI A (jellyfish and methane-like). The entropy of the J = 1/4 ground state is thus extensively large and higher than the one of spin ice [Fig. 18.(b)]. The absence of symmetry breaking and the very high entropy justify why the passage from the high-temperature paramagnet to the jellyfish regime is a crossover and not a phase transition, as confirmed by the specific heat C h and the magnetic susceptibility χ zz of Fig. 18.(a,c).

The variance of the magnetization ∆M 2 z ≡ M 2 z -M z 2 ≡ N χ zz T directly measures the build-up of correlations in classical spin liquids, visible in the Curie-law crossover of Fig. 18.(a) 42 . The fact that it reaches the value of 1/3 (within error bars) as T → 0 + , i.e. as for a standard paramagnet, is remarkable and consistent with the fact that the long-range dipolar correlations of the two-in two-out phase 71 are cut off by the jellyfish and methane-like structures.

In the absence of a critical point, the Husimi-cactus calculations provide a very good approximation of the real system. This is especially visible in Fig. 18.(e), where the temperature dependence of the charge density is essentially indistinguishable from simulations. There is however a small but clear difference for the second-order cumulant of the energy and magnetisation [see 18.(a,c)]. This small difference is due to the absence of closed loops of size 6 and higher in the Husimi-cactus approximation, which prevents the emergence of jellyfish structures. C. On the way to the jellyfish: J 1/4

The increase of single-charge density in Fig. 16 can now be understood by the emergence of charge structures precisely at J = 1/4. However, even though these structures are excitations for J < 1/4, they visibly play a role in the relaxation just below 1/4 [see the J = 0.24 curve in Fig. 16]. We have confirmed this influence by Monte Carlo simulations in Fig. 19, where the low-temperature single-charge density ρ 1 tends to zero with increasing Monte Carlo time. Despite the use of parallel tempering and worm algorithm [see appendix B], and a very slow decrease of the temperature during the preliminary Monte Carlo equilibration -made of up to 10 million Monte-Carlo steps (MCs) -simulations require 100 million MCs to reach the equilibrium vanishing value of ρ 1 . Such slow dynamics in absence of disorder is well-known in spin ice 13,14,20,25,27,28,[72][73][START_REF] Levis | Pierre et Marie Curie[END_REF] , even without any thermal showing partial ordering for J = 1/4, but not J < 1/4. As J increases, MQ.Q remains strictly positive to lower and lower temperature, showing the growing influence of the jellyfish phase. The agreement between analytics and simulations is quantitative for J = {0, 0.22} and semi-quantitative for J = 0.24, where the influence of jellyfish structures at finite temperature comes into play -especially the ring of single charges.

or field quench, but relatively rare for a phase with a finite density of charges, which are usually the source of dynamics [see e.g. the third time domain of Fig. 7].

Both single and double charges completely vanish as T → 0 for 0 J < 1/4 [see Fig. 21], but with a pronounced shoulder at finite temperature as one gets closer to the 1/4 value. This shoulder marks the change from the low-temperature two-in two-out regime, and the entropically induced apparition of jellyfish structures at intermediate temperatures. The crossover temperature is characterized by the energy scale E jellyfish of Eq. ( 15). This behaviour is most visible in the Curie-law crossover of 2/3 value characteristic of spin ice 42 , but delayed by an order of magnitude in temperature between J = 0 (yellow curve) and J = 0.24 (red curve). The same abrupt crossover is visible in

M Q.Q [see Fig. 20.(b)].
The presence of a crossover instead of a phase transition is confirmed by Husimi-cactus calculations, but could not be completely ruled out by simulations for 0.24 < J < 0.25. On a side note, the small negative value of M Q.Q for J = 0 (nearest-neighbor spin-ice model) is a consequence of the entropic Coulomb interactions between topological defects in spin ice 46,68,[START_REF] Castelnovo | [END_REF] . The magnitude of entropic interaction is proportional to T , and accordingly, the negative M Q.Q is conspicuous in the higher temperature region.

D. The ground state at J > 1/4

We also briefly discuss the nature of ground state for J > 1/4. As for J > 1/4, the energy is minimized by

-2 0 2 [hh0] -4 -2 0 2 4 [00k] T = 0.1 -2 0 2 [hh0] T = 0.3 -2 0 2 [hh0] 0.5 1 1.5 T = 1 -2 0 2 [hh0] -4 -2 0 2 4 [00k] T = 0.1 -2 0 2 [hh0] T = 0.3 -2 0 2
[hh0] 0.5 maximizing the number of bonds linking two neighboring same-sign charges. The candidate of lowest-energy configurations found by our Monte Carlo simulations is made of methane-like clusters [see Fig. 24], with exactly one such cluster per cubic unit cell. The clusters are made of 5 charges -a central double charge surrounded by four single charges -and are separated from each other by vacuum tetrahedra (two-in two-out). The densities of vacuum tetrahedra, single and double charges are thus respectively 3/8, 1/2 and 1/8. The ground-state energy normalized per number of spins is 1-4J. The intervening vacuum tetrahedra enforce a non-trivial frustrated constraint on the sign of these clusters. They impose that the four methane-like clusters connected to a given vacuum tetrahedron have to be locally charge neutral; two clusters are made of positive charges, the other two are made of negative charges. This imposes a strong constraint on the ground-state ensemble, whose degeneracy is sub-extensive. Interestingly, this ground state made of methane-like clusters is reminiscent of the 32-sublattice order that has been observed in a spin-ice type Kondo lattice 76 .

E. Neutron-scattering signature

Our final question concerns the experimental detection of the jellyfish regime. A possible fingerprint would be the Curie-law crossover of Fig. 18.(a). But since χ zz T reaches a value of 1/3 identical to a standard paramagnet, the real signature lies in the bump of the crossover which unfortunately could be modified by perturbations pertinent at intermediate or high temperatures such as finite single-ion anisotropy.

On the other hand, neutron scattering should be able to detect the fingerprints of the jellyfish structures. Neutron scattering measures the structure factor of the material

S(q) = 1 N i S ⊥ i e iq•ri 2 , (24)
where S ⊥ i are the spin components orthogonal to the vector q and r i is the position of site i. The structure factor obtained from neutron scattering for J = 1/4 is shown in Fig. 22 (top panels). At low temperature, distinctive half-moon features of scattering appear. These features are relatively broad and surroundi.e. do not sit onthe Brillouin zone centers. They persist at higher temperatures, albeit more diffuse, and are confirmed by a phenomenological analysis taking into account only the contributions from the ring part of the jellyfish [see Fig. 23]. These half-moon scattering motifs are very exotic for a spin-ice model and can thus serve as a solid signature of the jellyfish regime in experiments.

When stepping away from the J = 1/4 model, one recovers the characteristic pinch points of the spin-ice Coulomb phase at very low temperature [see Fig. 22 (bottom panels) for J = 0.22], as predicted from the discussion of section VI C. At higher temperatures (T = 1), the higher entropy of the jellyfish regime wins over the Coulomb phase, and one recovers the half-moon fea- tures of J = 1/4. In particular this broad scattering precisely replaces the pinch points, consistent with the correlated density of charges in the jellyfish regime 77,78 . Intriguingly, this motif is reminiscent of the crescentmoon pattern observed in the dynamical structure factor of the kagome antiferromagnet with classical Heisenberg spins 43,44 . In both cases, these "moon" patterns surround the Brillouin-zone centers where pinch points take place in a Coulomb phase. Hence these patterns signal the development of a characteristic length scale away from the Coulomb-phase physics. In kagome, these motifs are due to magnon excitations at finite energy 43,44 -which is why they are seen in the dynamical structure factor -while here they are due to the presence of a finite density of defects in the ground state -which is why they are seen in the equal-time structure factor down to very low temperature [Fig. 22].

VII. CONCLUSION

Motivated by the coupling between itinerant electrons and localized spins in spin-ice materials, we have investigated the physics of Ising spins on the pyrochlore lattice with first, second-and third-neighbour interactions [see equation (2)]. We have focused on the line where J 2 = J 3 = J, where our model can be rewritten in terms of charge degrees-of-freedom, giving rise to the nearest-neighbour dumbbell model [see equation (4)].

The two-in two-out Coulomb phase remains the ground state of the model over a broad range of parameters (-0.5 < J < 0.25). By varying J within this window, the relaxation processes after a field quench in the [111] direction display a rich diversity of nonequilibrium phenomena, with e.g. glassy behavior in absence of disorder due to multi-scale energy barriers, which can be quantitatively understood by a microscopic approach.

For small values of J (-0.25 < J < 0), the magnetization decay shows markedly slow relaxation when compared to the charge density. It gives rise to (multiple) magnetization plateaux over several orders of magnitude in time due to the exhaustion of charges. When approaching from J = -0.25 to the phase-diagram boundary J = -1/2, the dynamics changes drastically and is dominated by the diffusion of double charges and vacuum tetrahedra. As a result, the initial charge order imposed by the field quench at t = 0 persists while the magnetization quickly vanishes. The resulting metastable phase is known as the fragmented Coulomb spin liquid (FCSL) 45,46 , where long-range (charge) order co-exists with a (Coulomb) spin liquid. This physics of fragmentation has been observed in a recent series of experiments 47,49,50 but remains elusive in a standard spin-ice model 45,46,57,58 . By studying a variant of spin ice with further nearest-neighbors, our work provides a working example on how out-of-equilibrium phenomena can be used to engineer unconventional properties in frustrated magnetism.

As for positive J, the dynamics develop an interesting J/T scaling law due to the repulsion between oppositesign charges, which hinders their ability to annihilate pairwise. As approaching J = 1/4, this scaling breaks down as the system acquires a finite charge density in thermodynamic equilibrium, due to the high population of composite excitations which we have pictorially termed methane and jellyfish, as discussed in Sec. VI A. These structures, once created, decay slowly because of the high energy barriers accompanied with their collapse [see Sec. G].

While the jellyfish structure is non-magnetic, the central ring carries an emergent chiral degree of freedom, characterized by a toroidal moment. Accordingly, the jellyfish structure is able to encode the information of time-reversal symmetry breaking with this toroidal degree of freedom. This feature, combined with the long lifetime, implies a potential relevance to the spontaneous Hall effect observed in Pr 2 Ir 2 O 7 .

Aiming at the detection of jellyfish, we combined numerical/phenomenological analyses to show that the "half-moon" structure in the S(q) is a clear signature of the jellyfish structure. Experimentally, the detection of various types of composite excitations has been carried out for frustrated magnets, such as Zn/Mg/HgCr 2 O 4 and Tb 2 Ti 2 O 7 [79][80][81][82][83][84] . Our analysis indicates that the neutron scattering is a promising experimental probe to detect the jellyfish. In fact, though the microscopic determination of J 2 and J 3 is unavailable for Pr 2 Ir 2 O 7 , once the jellyfish structure is detected through S(q), it will impose a strong constraint on the values of J 2 and J 3 , and serve as a smoking gun evidence for this compound located close to the "jellyfish point", J = 1/4. Accordingly, it will lead to clarifying the origin of spontaneous Hall effect in Pr 2 Ir 2 O 7 .

The natural next step of this analysis would be to investigate what happens for J > 1/4. In this region the two-in two-out Coulomb phase is not part of the ground state anymore, and the ordering mechanism is completely dominated by the attraction between same-sign charges, inducing strong steric and kinetic constraints. Our simulations indicate an increase in single and double charges below a first order transition, but we found that this phase is very difficult to thermalize.

Beyond the original motivation of RKKY interactions in metallic systems, our model directly applies to insulating pyrochlores with farther-neighbour exchange. Such interactions might be present in spin-ice materials [61][62][63] and be of importance for quantum spin ices where virtual crystal field excitations can induce coupling beyond nearest-neighbors 64,65 . The J 2 = J 3 condition is an arguably strong constraint, but very useful by its simplicity in order to serve as a basis for perturbative approaches.

Our paper also provides a working example of what happens when it is possible to tune the interactions between topological defects away from their natural setting. Indeed, spin ice has often been described as a magnetic analogue of an electrolyte [START_REF] Castelnovo | [END_REF][85][86][87][88] , a picture which relies on the Coulomb attraction between opposite-sign charges. The possibility to reverse the sign of this interaction, even just at the nearest-neighbour level, drastically changes the physics of the problem. Clusters of quasiparticles can be stabilized because i) charges of same signs cannot annihilate by pair and ii) charges of opposite signs -which could annihilate -are repelled from each other. Such questions are certainly not limited to spin ice, but extend naturally to phases prone to emergent quasi-particles such as quantum spin liquids, topological phases and artificial gauge fields.

Last but not least, the apparition of the FCSL phase demonstrates the promises of engineering macroscopic states via nonequilibrium techniques 29 . We hope our work will further motivate the exploration of nonequilibrium phenomena in geometrically frustrated magnets, ranging from topological ergodicity breaking 89 to glassiness without disorder 20,28 , and with connections to kagome systems 90,91 including their realization in artificial spin ice [51][52][53][54][55] and the 16-vertex model [72][73][START_REF] Levis | Pierre et Marie Curie[END_REF] .

Note: after submission of our work, a preprint by Rau & Gingras has appeared on arXiv 93 studying the J ∈ [0 : 0.25] region of our model, especially focused on the dynamics of the J = 1/4 point. By introducing a charge variable, Qp ≡ η p1 + η p2 + η p3 + η p4 , the first term can be transformed into

F G η F η F! η F" η G η G! η G" η G η F = >
1 2 p (η p1 + η p2 + η p3 + η p4 ) 2 = 1 2 p Q2 p . (A2)
The second term can be written as

J p,q (η p1 + η p2 + η p3 + η p4 -η [p,q] )× (η q1 + η q2 + η q3 + η q4 -η [p,q] ) = J p,q ( Qp -η [p,q] )( Qq -η [p,q] ) = J p,q ( Qp Qq -Qp η [p,q] -Qq η [p,q] ) + Const. = J p,q Qp Qq -J p Qp q=n.n. of p η [p,q] + Const. = J p,q Qp Qq -J p Q2 p + Const. (A3)
The equations (A1) and (A3) add up to the representation,

H = 1 2 -J p Q2 p + J p,q Qp Qq . (A4)
Finally, by introducing a new charge variable:

Q p = Qp , Q q = -Qq , (A5)
for upward (downward) tetrahedra, p (q). We end up with the form,

H = 1 2 -J p Q 2 p -J p,q Q p Q q . ( A6
)
• a two-in two-out tetrahedron: there are two equivalent outward spins to exit the tetrahedron, chosen with probability 1/2.

• a three-in one-out tetrahedron: there is one outward spin to exit the tetrahedron, chosen with probability 1.

• a one-in three-out tetrahedron: there are three equivalent outward spins to exit the tetrahedron, chosen with probability 1/3.

• an all-in tetrahedron: there are no possibilities to exit the tetrahedron. The worm is stopped and erased.

The procedure is repeated until the worm closes. Because such a worm is made of as many inward spins as outward spins for every tetrahedron it encounters, flipping all spins in the worm does not modify the position nor the sign of the charges. Since the Hamiltonian of our model can be written in terms of charge degrees-offreedom only [se equation ( 4)], flipping a worm does not modify the energy of the system, and detailed balance is thus naturally respected.

As for the simulation setup, let a Monte Carlo step be defined as N t single-spin-flip attempts and 10 worm updates, with N t , the number of tetrahedra. Monte Carlo dynamics include 1 parallel tempering every 100 Monte Carlo steps. Let t max be the total number of Monte Carlo steps in a given simulation. The thermalization setup of the simulations consist in a slow annealing from high temperature to the temperature T of measurements during t max /10 steps, followed by further thermalization at temperature T during t max /10 steps. For T 1, the temperature scale used in parallel tempering was logarithmic in order to connect the high-and low-temperature behaviors. For T 1, the temperature scale was linear in order to precisely describe the entire low-temperature regime.

Appendix C: Husimi-cactus approximation Thermodynamic properties of the J 1 -J 2 -J 3 model can be semi-analytically evaluated by Husimi-cactus approximation. A Husimi cactus is analogue to a Bethe lattice where each vertex is replaced by a cluster of spins, here made of tetrahedra [Fig. 26]. This method is equivalent to approximating the pyrochlore lattice with its loopless variant, the tetrahedron Husimi cactus (THC) as shown in Fig. 25 (b). Thermodynamic quantities, such as specific heat and magnetic susceptibility, as well as charge density can be evaluated exactly on this network, and its temperature dependence rather precisely matches the result of Monte Carlo simulation [see section VI].

For the sake of explanation, we here summarize our convention. Firstly, we introduce four sublattices and assign anisotropy axis, d α (α = 1 -4), as defined in section III A, to discuss magnetic susceptibility in this scheme.

The sublattice index of site i is denoted as α i . To discuss the tetrahedron-based quantities, we choose one tetrahedron, T 0 , and consider the array of tetrahedra, starting from T 0 . In particular, we define T N as the N -th tetrahedron from T 0 [Fig. 25(b)]. Each tetrahedron T j has 2 4 possible states ζ j (= 1 -16) corresponding to the values of spins, η, on its four sites. We assume that the state ζ has charge Q[ζ], and the spin at sublattice α is denoted as η α [ζ].

Partition function on half cactus

The basic building blocks of this approximation are the partial partition functions evaluated on half cactus, which are obtained by terminating the tetrahedron Husimi cactus at the "root site", O, as shown in Fig. 26 (a). In particular, we call the tetrahedron involving the site O, "root tetrahedron". We consider the Hamiltonian (2) on this network, and sum up the Boltzmann factors for all the configurations of {η i }, to obtain the partition function, z. In particular, we focus on the partial summation of Boltzmann factors, with a fixed configuration of η's on the root tetrahedron. As shown in Fig. 26 (b), we define z 4 as the partial sum of Boltzmann factors, provided all η's have the same sign on the root tetrahedron. Similarly, we define z 3A , z 3B and z 2 as the partial summation with the 3-in 1-out or 2-in 2-out configurations on the root tetrahedron. We distinguish z 3A and z 3B : For the former, the root site has a minority spin, while not for the latter. Taking account of the degeneracy of configurations on the root tetrahedron, we obtain

z = 2z 4 + 6z 3B + 6z 2 + 2z 3A . (C1)
The remarkable feature of the half cactus is its selfsimilarity. By removing the root tetrahedron, one can generate three separate copies of the original half cactus network. By using this feature, one can prove that the partial partition functions satisfy the following recursive relations:

             z 4 = e -16β( 1 2 -J) e -16βJ z 4 + 3e -8βJ z 3B + 3z 2 + e 8βJ z 3A 3 , z 3B = e -4β( 1 2 -J) e -8βJ z 4 + 3e -4βJ z 3B + 3z 2 + e 4βJ z 3A 2 × e 8βJ z 4 + 3e 4βJ z 3B + 3z 2 + e -4βJ z 3A , z 2 = z 4 + 3z 3B + 3z 2 + z 3A 3 , z 3A = e -4β( 1 2 -J) e -8βJ z 4 + 3e -4βJ z 3B + 3z 2 + e 4βJ z 3A 3 . (C2)
By solving these equations, we can get three ratios, z 4 /z 2 , z 3B /z 2 and z 3A /z 2 , which will serve as basic building blocks to obtain thermodynamic quantities, as we discuss below.

Occupation rate

Next, we consider the probability that a particular spin configuration is realized in a single tetrahedron. The probability that double charge, single charge, and vacuum state is realized in a certain tetrahedron, is proportional to p 2 , p 1 , and p 0 , respectively, which are given by

p 2 = e -16β( 1 2 -J) e -16βJ z 4 + 3e -8βJ z 3B + 3z 2 + e 8βJ z 3A 4 , p 1 = e -4β( 1 2 -J) e -8βJ z 4 + 3e -4βJ z 3B + 3z 2 + e 4βJ z 3A 3 , × e 8βJ z 4 + 3e 4βJ z 3B + 3z 2 + e -4βJ z 3A p 0 = z 4 + 3z 3B + 3z 2 + z 3A 4 , (C3)
in which z 4 , z 3B , z 2 and z 3A are the solution of the selfconsistent equation (C2). By considering the number of degeneracy, the probabilities for each spin configuration are given by

ρ 2 = 2p 2 N , ρ 1 = 8p 1 N , ρ 0 = 6p 0 N , (C4)
with N = 2(p 2 + 4p 1 + 3p 0 ).

Correlation function

In order to evaluate physical quantities, such as specific heat and magnetic susceptibility, one needs spin or charge correlation function. To obtain these quantities, we need the conditional probability, P N (ζ ′ |ζ): the probability that the tetrahedron T N takes the state ζ ′ , given that T 0 is in the state ζ.

With P N (ζ ′ |ζ), one can write down the spin correlation between site i and j:

η i η j = ζ,ζ ′ P N (ζ ′ |ζ)η αi [ζ]η αj [ζ ′ ]. (C5)
Here, we assumed the site i (j) belongs to the sublattice α i (α j ) of tetrahedron T 0 (T N ). And similarly, the charge correlation between the tetrahedra, T 0 and T N :

Q T0 Q TN = ζ,ζ ′ P N (ζ ′ |ζ)Q[ζ]Q[ζ ′ ], (C6)
To obtain P N (ζ ′ |ζ), one can resort to the method of transfer matrix. Suppose a 16 × 16 matrix, K, whose (ζ, ζ ′ ) component, K(ζ|ζ ′ ) means the conditional probability that the tetrahedron T j+1 takes the state ζ, given that T j is in the state ζ ′ . This probability does not depend on the index j. P N (ζ N |ζ 0 ) can be expressed as

P N (ζ N |ζ 0 ) = ζ1•••ζN-1 K(ζ N |ζ N -1 )K(ζ N -1 |ζ N -2 ) • • • K(ζ 1 |ζ 0 ) = K N ηN η0 . (C7)
With the eigenvalue of K, λ γ , and corresponding right and left eigenvectors as u γ and v γ , one can write

P N (ζ N |ζ 0 ) = 16 α=1 λ N γ u γ ζN v γ ζ0 . (C8)

Energy and related quantities

By taking the charge representation, the internal energy of the system is given by

E = H = 1 2 -J p Q 2 p -J p,q Q p Q q . (C9)
To evaluate internal energy, the thermal expectation value of square charge, Q 2 p , and nearest-neighbor charge correlation, Q p Q q are required. These quantities can be easily obtained from the equation (C6). We can obtain specific heat C by taking the numerical derivative of E/N site in terms of temperature. We also estimate M Q.Q , defined with eq. ( 18), from Q p Q q . These quantities are plotted in Fig. 18, together with the results of Monte Carlo simulation.

Magnetic susceptibility

The magnetic susceptibility, χ zz , can be obtained from the spin correlation introduced in equation (C5). By using the anisotropy axis d α associated with each sublattice, α, one can write

χ zz = 1 T j η 0 η j d z α0 d z αj . (C10)
Appendix D: Dynamical processes

For convenience, we summarize the important dynamical processes and the energy barriers, ∆ m , associated with each process. ∆ m is defined as ∆

m = E f m -E i m , where E f(i)
m is the total energy of the system after (before) the process m takes place. For small negative J: -0.2 < J < 0.0, the magnetization exhibits a wide plateau until it is terminated due to the pair creation and dissociation processes of charges.

Here, we note one additional feature appearing in the plateau region. At low enough temperatures, the plateau splits into two regions, as shown in Fig. 28. The magnetization first drops to a value M * ∼ 0.198, then M stays the same value for a certain time range. After a while, the magnetization shows partial relaxation, and exhibits the second plateau. After a long time determined from the pair creation process, the second plateau also collapses. As temperature is raised, this first plateau shrinks, i.e., the collapse of first plateau occurs earlier at higher temperature, and it finally merges with the initial drop of magnetization.

The origin of this first plateau can be attributed to the dissociation process of non-contractible charge pairs. For negative J, charges with opposite sign attract with each other in the nearest-neighbor sites, which leads to the energy barrier ∆ 4 = 4|J| for their dessociation. Consequently, once two charges are trapped in a noncontractible position, it will take quite a long time to overcome the energy barrier to escape from that position. This mechanism is similar to that of slow relaxation found in the dipolar spin ice model, where the energy barrier is attributed to the long-range dipolar interaction 28 .

As plotted in Fig. 29, the time necessary to terminate the first plateau is well fitted by an Arrehnius law: exp(∆ 4 /T ), with the value of energy barrier associated with the dissociation process.

Appendix F: Absence of Hall effect in the charge crystal

We can show that the Hall conductivity vanishes for the FCSL, on the basis of spatial symmetry of this state. To prove this, it is instructive first to consider the case of all-in/all-out crystal [Fig. 30]. Here, we define all-in (allout) crystal as the state where all the upward tetrahedra take the all-in (all-out) configuration.

Firstly, the transverse conductivity σ xy satisfies the relation

J y = σ xy E x . (F1)
Here, J y and E x are electric current in the y direction, and electric field in the x direction.

On the basis of this property, we prove that Hall conductivity vanishes in the all-in/all-out ordered state. To begin with, we note that the all-out state is obtained by reversing all the spin directions from the all-out state. In other words, T : Time-reversal operation maps the allin state to the all-out state, and vice versa. Next, we define the spatial operation R, by combining rotational and translational operations. Firstly, we consider the 90 • rotation (of spin and position) around an line parallel to z-axis, and through one of the sites with sublattice B [Fig. 30]. Then, we successively translate the system by R 0 . This combined operation, R, maps a group of βchains to those perpendicular to the former. At the same time, this operation interchanges the upward and downward tetrahedra, and accordingly, map the all-in state to the all-out state, and vice versa. From the viewpoint of J y and E x , the operation, R, acts as The same proof holds, if the exchange of upward and downward tetrahedra can be considered equivalent to time-reversal operation. In the FCSL, the upward and downward tetrahedra are occupied with 3-in 1-out and 1in 3-out configurations, or vice versa. Accordingly, while the FCSL is spin-disordered state, on a macroscopic scale, the exchange of upward and downward tetrahedra lead to (color online).

J y → -J x , E x → E y . (
The time dependence of magnetization is plotted for J = -0.10 at T = 0.008, 0.01, 0.02, 0.04, 0.06, 0.10, 0.20 and 0.40. At low temperatures, the magnetization shows the second plateau, at which M ∼ 0.198, which is slightly larger than the next plateau at M ∼ 0.164. the time-reversal conjugate state of initial state, resulting in the absence of Hall signal in FCSL.

Appendix G: Collapse of jellyfish

The jellyfish structure is quite stable. Accordingly, the toroidal moment associated with its ring part has quite a long lifetime. Here, we discuss two dominant processes which collapse the ring structure and destroy the toroidal moment, accordingly. Both processes lead to a lifetime of the order of at least ∼ exp(J 1 /T ) = exp(1/T ). Accordingly, this structure has a quite long lifetime at suf- ficiently low temperatures.

energy barriers a. double charge-vacuum creation

Since the jellyfish is composed of the single charges with the same sign, it is subject to the creation of double charge-vacuum pair from the adjacent two single charges. In particular, if it happens in the ring part, this event immediately disturbs the chiral magnetization flow. The energy cost for this pair creation, ∆ p-c , is estimated as ∆ p-c = 4(1 -J), if the jellyfish has no branches. In case the branches take the optimal configuration for the collapse, it is reduced to ∆ min p-c = 4(1 -3J), but still this lower-bound value takes ∼ 1 near J = 1/4, so this decay process hardly occurs at T ≪ 1.

b. pair annihilation of opposite charges

Pair annihilations of opposite charges is another possible process leading to the decay of jellyfish. Due to the topological constraint, the total charges in the system must be equal to zero. Accordingly, if there is a jellyfish composed of a certain numbers of positive charges, there must be the same numbers of negative charges elsewhere in the system, which is likely to form a jellyfish or methane cluster to lower their energy, near J = 1/4. Under the circumstances, to destroy the positive jellyfish, one can consider the following three-step process: (i) one monopole is detached from the negative jellyfish/methane, (ii) it migrates to the positive jellyfish, and (iii) makes a pair-annihilation with one of the positive charges composing the target positive jellyfish. The energy barrier accompanying (i) is at least 4J. (iii) also leads to the energy barrier of 4J, since it requires the opposite charges to be placed in nearest-neighbor sites. So the total energy cost for this pair-annihilation, ∆ p-a , is estimated as ∆ p-a = 8J. However, if there are finite numbers of stray charges, possibly due to thermal excitation, the step (i) is not necessary. In this case, ∆ p-a is reduced to ∆ min p-a = 4J, and it gives a comparable value to the lower bound of double charge-vacuum creation energy barrier: ∆ min p-a ∼ ∆ min p-c ∼ 1, near J = 1/4.

kinematic constraint

In addition to the energy barriers discussed above, there are substantial contributions from the kinematic constraint. Firstly, the double charge-vacuum creation costs the energy ∆ p-c ≥ ∆ min p-c = 4(1 -3J). However, the lower-bound value, ∆ min p-c is available, only if a special configuration is realized on a ring [see Fig. 31], which is statistically rare.

Secondly, the lifetime estimated from the pairannihilation process also requires substantial kinematic correction. This process involves the migration of charge, which we call step (ii) above. If the charge density is low, the charge has to migrate for quite a long time until it successfully go overs the energy barrier due to the step (iii).

These kinematic features contrast with the dynamics at J < 0, where the dynamical bottleneck is given by a single energy barrier in most cases, and the dynamics after overcoming the barrier is avalanche-like. The kine- matic constraints give a substantially large correction to the lifetime ∼ exp(1/T ) estimated purely from energetic consideration, and contributes to the stability of jellyfish structures.

Appendix H: Effects of charge rings on structure factor

To consider how the jellyfish affects the magnetic structure factor, we evaluate the contribution from the ring part of jellyfish to the spin correlation function. Namely, we calculate S ring (q) = S ring q • S ring -q , (H1)

where

S ring q = i∈ring S i e -iq•ri . ( H2
)
The ring part can be classified into four types, according to its orientation: the hexagonal ring is perpendicular to [111], [1-1-1], [-11-1] and [-1 -11]. We assume these four types of rings are equally populated, without any correlations on their positions. For each ring, we take account of the contribution from the 18 site-cluster, as shown, in Fig. 32. For example, from the ring perpendicular to [111], we obtain

FIG. 2 .

 2 FIG. 2. (color online). (a) An example of spin configuration.The spin pairs connected by interactions, J1, J2 and J3 are shown with dashed lines. J3 does not include interactions across hexagons. (b) The charge distribution corresponding to the spin configuration shown in (a). We also show here the two classes of tetrahedra, which we call "upward" and "downward", according to their orientations. (c) Examples of spin configurations are shown, corresponding to each value of tetrahedral charge. For simplicity, we refer to the tetrahedra with Q = 0, ±2 and ±4 as vacuum, positive(negative) single charge, and positive(negative) double charge, respectively.

17

 17 FIG. 3. (color online). Ground-state phase diagram of Hamiltonian (4) as a function of J. This corresponds to the J1 -J2 -J3 spin-ice model for J2 = J3 = J, and J1 = 1. The Coulomb phase is stabilized over a broad region (-1/2 < J < 1/4). For J < -1/2, it gives way to the antiferromagnetic phase made of alternating all-in and all-out tetrahedra (double charges). At J = 1/4, the jellyfish point discussed in section VI appears. For J > 1/4, the ground state ensemble is sub-extensively degenerate; spin configurations are made of one methane-like cluster per cubic unit cell [see Fig. 17.(b)].

  FIG. 3. (color online). Ground-state phase diagram of Hamiltonian (4) as a function of J. This corresponds to the J1 -J2 -J3 spin-ice model for J2 = J3 = J, and J1 = 1. The Coulomb phase is stabilized over a broad region (-1/2 < J < 1/4). For J < -1/2, it gives way to the antiferromagnetic phase made of alternating all-in and all-out tetrahedra (double charges). At J = 1/4, the jellyfish point discussed in section VI appears. For J > 1/4, the ground state ensemble is sub-extensively degenerate; spin configurations are made of one methane-like cluster per cubic unit cell [see Fig. 17.(b)].

  at sufficient times and make a sample average. In the actual calculations, we typically choose the system size: N site = 16 × 16 × 16 × 4, number of steps: N step = 1000N site , and the number of samples: N sample = 100. B. Nearest-neighbor spin-ice model 1. Temperature quench as a test

FIG. 4 .

 4 FIG. 4. (color online). Time evolution of the single charge density, ρ1, for the nearest-neighbor spin ice model (J = 0) after a temperature quench from T = 10 to T = 0 starting at t = 0. ρ1 is defined as the sum of positive and negative single charges, divided by the number of tetrahedra. For the definition of time, see section V A. The solid curve is the mean-field solution of a creation-annihilation reaction equation:ρ1(t) = ρi 1 + Kρit , where K = 9 √ 3 16in agreement with

C.

  Field quench of the nearest-neighbour dumbbell model (J = J2 = J3)

FIG. 5 .FIG. 6 .

 56 FIG. 5. (color online). The spin configuration under the [111] saturation field. All the tetrahedra are occupied with 3-in/1out or 1-in/3-out configurations. From the [111] direction, the pyrochlore lattice can be regarded as an alternate stacking of kagome and triangular lattices. This saturated configuration belongs to the extensively degenerate ensemble of fragmented Coulomb spin liquid states. However, it is important to understand that a large fraction of FCSL states have zero or small magnetisation.

FIG. 7 .

 7 FIG. 7. (color online).Field quench process at J = -0.10 and T = 0.10. The charge density and magnetization are plotted together with the sublattice magnetizations on kagome and triangular layers. The time domain is divided into four regions I, II, III and IV, as shown with dashed lines, according to the qualitative nature of the dynamics. The sublattice magnetizations, m triangular and m kagome , are measured along the[111] direction and are normalized by the total saturated magnetization, giving rise to the value of 0.5 at t = 0.

FIG. 9 .FIG. 10

 910 FIG. 9. (color online). The extinction time tp of the magnetization is plotted against the inverse of temperature T for J = -0.05, -0.10 and -0.15. tp is defined as the time at which magnetization takes half of the plateau value. The tp's follow the predicted Arrhenius law, Eq. (13) with log tp = -8.00531 + 4.4/T , -8.12707 + 4.8/T and -8.31068 + 5.2/T for J = -0.05, -0.10 and -0.15, respectively.

FIG. 11 .FIG. 12 .

 1112 FIG. 11. (color online). Illustration of the nucleation process for -0.5 < J < -3/7 where the defects propagate in the [111] direction. It allows for spin decorrelation from the initial saturated configuration, while preserving the long-range charge order [see Fig.14]. Note that the dynamics in kagome plane cannot be ignored for J ≥ -3/7. See the main text for details.

49 FIG. 13 .

 4913 FIG. 13. (color online). The extinction time tp of the total magnetization for J = -0.45 and -0.49 as a function of the inverse temperature. The low temperature behavior of tp is fitted with the Arrhenius law: ∝ exp[(4 + 8|J|)/T ], with log tp = -8.5 + 7.6/T and -8.5 + 7.92/T for J = -0.45 and -0.49, respectively.

4 .

 4 0.00 < J < 0.25

FIG. 15 . 25 FIG. 16 .

 152516 FIG. 15. (color online). The relaxation of the charge density and magnetization after the field quench are plotted for several combinations of J and T , keeping the ratio J/T = 1.25 fixed: (J, T ) = (0.05, 0.04), (0.10, 0.08), (0.15, 0.12) and (0.20, 0.16). All relaxation data quantitatively collapse onto a single curve, suggesting a scaling relation with a single parameter, J/T .

FIG. 17 .

 17 FIG. 17. (color online). Schematic figures of the minimal clusters: (a) a ring made of six tetrahedra for the jellyfish and (b) methane-like structure.

  FIG. 18.(color online). Absence of long-range order for J = 1/4, probed by Monte Carlo simulations (black circles) and Husimi-cactus calculations (red dashed lines). (a) Magnetic susceptibility along the z-axis times temperature. We observe a Curie law crossover without any singularities (i.e. without any phase transition), characteristic of a classical spin liquid at low temperature 42 . (b) Entropy of the jellyfish model compared to the nearest-neighbor spin-ice model (yellow crosses). The black dashed lines show the lowtemperature limits from Monte Carlo and confirm the higher entropy of the jellyfish ground state. (c) Specific heat. (d) Quasi-order parameter defined in eq. 18, showing short-range ordering with a finite fraction of same-sign nearest-neighbour monopoles. (e) Density of tetrahedra configurations. All observables are normalized per number of spins, except for panel (e) which is normalized per number of tetrahedra. The temperature axis is on a logarithmic scale for panels (a,b).

1 FIG. 19 .

 119 FIG. 19. (color online). Temperature dependence of the density of single charges probed by Monte Carlo simulations for J = 0.24, just below the jellyfish phase. Each curve corresponds to a different simulation time, from 10 5 to 10 8 Monte Carlo steps. This shows that even if the thermalization process is very slow, the density of charges eventually vanishes for J < 1/4. See appendix B for details on the simulations.

FIG. 20

 20 FIG. 20. (color online). Temperature dependence of the thermodynamics probed by Monte Carlo simulations for J ∈ {0, 0.22, 0.24, 0.25} (yellow, orange, red and black respectively). The dashed lines are analytical Husimi calculations. (a) Magnetic susceptibility along the z-axis times temperature. The horizontal temperature axis is on a logarithmic scale. We observe two different kinds of Curie law crossover. For J < 1/4, χzzT ultimately goes to ≈ 2/3, characteristic of the two-in two-out Coulomb phase42 . For J = 1/4, χzzT shows a clear downturn towards ≈ 1/3, signature of a classical spin liquid different from the Coulomb phase. (b) Order parameter MQ.Q defined in eq. 18, showing partial ordering for J = 1/4, but not J < 1/4. As J increases, MQ.Q remains strictly positive to lower and lower temperature, showing the growing influence of the jellyfish phase. The agreement between analytics and simulations is quantitative for J = {0, 0.22} and semi-quantitative for J = 0.24, where the influence of jellyfish structures at finite temperature comes into play -especially the ring of single charges.

2 FIG. 21 .

 221 FIG. 21. (color online). Temperature dependence of the densities of vacuum tetrahedra (a), single (b) and double (c) charges probed by Monte Carlo simulations for J ∈ {0, 0.22, 0.24, 0.25} (yellow, orange, red and black respectively) and Husimi calculations (dashed lines).

1 1. 5 T = 1 FIG. 22 .

 15122 FIG. 22. (color online). Structure factor S(q) in the [hhk] plane as measured by neutron scattering obtained from Monte Carlo simulations for J = 0.25 (top) and 0.22 (bottom) at T = {0.1, 0.3, 1}. Top: For J = 1/4, the pattern of the structure factor is relatively similar from T = 1 to 0.1, showing sharper intensities around the half-moon features upon cooling. Bottom: For J = 0.22 on the other hand, one moves continuously from the low-temperature spin-ice regime with sharp pinch points77 to the entropically favored jellyfish regime at T = 1 [see section VI C]. The half-moon features characteristic of the jellyfish regime sit at the same wavevectors as the pinch points in spin ice. The color scale is the same for all panels and starts at finite intensity to account for the background scattering.

  FIG. 23. (color online). Contribution from the charge rings to the spin flip part of the structure factor. We assume that all rings are made of six tetrahedra and that the four possible orientations of hexagonal rings are equally populated [see appendix H]. Remnants of "half-moon" structure can be clearly seen around the Brillouin zone centers, consistent with the results of Monte Carlo simulation [Fig. 22].

FIG. 24

 24 FIG. 24. (color online). Example of a ground-state configuration for J > 1/4 made of one methane-like cluster of spins per cubic unit cell. These clusters are separated from each others by vacuum tetrahedra. The vacuum tetrahedron shared by the four clusters displayed here is colored in yellow. It imposes the local charge neutrality of the methane-like clusters: two of them must be positive (blue), the other two being negative (red).

FIG. 25 .

 25 FIG. 25. (color online). (a) The tetrahedron cluster and the convention of site indices used for the derivation of nearestneighbour dumbbell model. (b) Representation of the Husimi cactus.

FIG. 26 .

 26 FIG. 26. (color online). (a) A schematic figure of half cactus, obtained by terminating the tetrahedron Husimi cactus at site O. (b)-(e) The configurations of η's on the root tetrahedron, corresponding to the partial partition function, (b) z4, (c) z3A, (d) z3B and (e) z2. ± signs correspond to η = ±1.

( 1 )( 2 )∆ 2 = - 4 ( 3 )( 4 )( 5 )∆ 9 = 8

 122434598 Creation of a "vacuum pair" from the [111] saturated state by flipping a spin in the kagome plane ∆ 1 = -4 -20J. (D1) Creation of another vacuum pair next to the initial pair created in process (1) Creation of a dipole (= a pair of charges) in the Coulomb phase ∆ 3 = 4 -4J. (D3) Separation of two charges (starting from a dipole created in (3)) ∆ 4 = -4J. (D4) Creation of a double charge pair by flipping a spin on a triangular plane from the [111] saturated state ∆ 5 = 12 + 12J. (D5) (6) spin flip just after (5): the pair of single charge and double charge becomes a pair of vacuum tetrahedron and single charge ∆ 6 = -8 -20J. (D6) (7) spin flip just after (6): the other pair of single charge and double charge becomes a pair of vacuum tetrahedron and single charge ∆ 7 = -8 -16J. (D7) (8) spin flip just after (6) and (7): vacuum-charge paper to single charge-double charge pair ∆ 8 = 8 + 16J. (D8) (9) spin flip just after process (1): creation of a doublesingle charge pair by flipping the spin in the triangular layer,

FIG. 27 .

 27 FIG. 27. (color online). Important dynamical processes. Each process is composed of a single spin flip. The flipped spin is highlighted with a yellow circle. The energy barriers associated with each processes are listed as equations (D1) to (D9).

  FIG. 28.(color online).The time dependence of magnetization is plotted for J = -0.10 at T = 0.008, 0.01, 0.02, 0.04, 0.06, 0.10, 0.20 and 0.40. At low temperatures, the magnetization shows the second plateau, at which M ∼ 0.198, which is slightly larger than the next plateau at M ∼ 0.164.

FIG. 29 .

 29 FIG. 29. (color online). The time for the termination of first plateau is plotted for J = -0.05, -0.10, -0.15 and -0.20 from bottom to top. Each data are well fitted with Arrheniustype temperature dependence: log t = -1.5 + 4J/T , as shown above.

FIG. 30 .

 30 FIG. 30. (color online). Schematic picture of the all-in crystal state. The green dashed line shows one of the β-chains, which will be transformed to the blue line after the operation RT (see the main text).

FIG. 31

 31 FIG. 31. (color online) The schematic figure of jellyfish with two branches, each of which has one charge, respectively. The magnetization flow is shown with red arrows. The reversal of spin a costs the energy ∆p-c = 4(1 -J). Meanwhile the energy cost of reversing spin b takes the lowest possible value, ∆ min p-c = 4(1 -3J).

FIG. 32

 32 FIG. 32. (color online) The schematic picture of charge ring, perpendicular to [111]. For the sites on the hexagonal ring, the sublattice indices and the relative positions from the center of the hexagon are shown.

ACKNOWLEDGMENTS

It is a pleasure to thank John Chalker for useful discussions. This work is supported by DFG via SFB 1143, Engineering and Physical Sciences Research Council (EP-SRC) Grant No. EP/G049394/1 (C.C.), the Helmholtz Virtual Institute "New States of Matter and Their Excitations," the EPSRC NetworkPlus on Emergence and Physics far from Equilibrium, the Okinawa Institute of Science and Technology Graduate University and by KAKENHI (Nos. 26400339, 24340076, 15H05852 and 15K13533).

Appendix A: Charge representation of the J1 -J2 -J3 spin ice model

Here, we derive the charge representation of the J 1 -J 2 -J 3 spin ice model at J 2 = J 3 = J. For this purpose, we rewrite the Hamiltonian (2), up to constant, in the following form:

In the first term, the summation is taken over the tetrahedra, p, while in the second term, the summation is over neighboring pairs of tetrahedra, p, q . The first term comes from the nearest-neighbor interactions which connect spins on the same tetrahedron. Meanwhile, the second term is due to the second-and third-neighbor couplings which connect spins on neighboring tetrahedra. For this rewriting, we named the sites on the neighboring tetrahedra, p and q, as p1 -p4 and q1 -q4, as shown in Fig. 25 (a). These tetrahedra share one site, which we call η [p,q] . In the convention shown in Fig. 25 (a), η [p,q] = η p4 = η q4 .

Appendix B: Monte Carlo simulations

For the characterization of the J = 1/4 model discussed in section VI [Figs. 18, 19, 20, 21 and 22], we made extensive use of classical Monte Carlo simulations with single-spin-flip dynamics, parallel tempering and worm algorithm.

The worm algorithm especially is a powerful method to decorrelate systems where single spin flip dynamics become inefficient because of a local ice rule 92 . In spin ice, this is the case for the two-in two-out Coulomb phase. In the present study however, we are especially interested in the presence of topological defects which break this ice rule. To circumvent the problem, we used an extension of the worm algorithm developed in Ref. [45].

The worm is made of a unidirectional closed chain of spins where all spins point in the same direction along the chain. It is constructed spin by spin until it closes on itself. Let us arbitrarily choose that during its construction, the worm enters a tetrahedron via a spin pointing inward. Then, the worm may enter S ring (q) = S [111] q + S [1-1-1] q + S [-11 -1] q + S [-1 -11] q . (H4) is proportional to the total contribution from the rings on structure factor.