# Violation of the fluctuation-dissipation theorem and effective temperatures in spin ice

Valentin Raban, Ludovic Berthier, Peter C.W. Holdsworth

► To cite this version:

Valentin Raban, Ludovic Berthier, Peter C.W. Holdsworth. Violation of the fluctuation-dissipation theorem and effective temperatures in spin ice. Physical Review B, 2022, 105 (13), pp.134431. 10.1103/PhysRevB.105.134431. hal-03657988

HAL Id: hal-03657988

https://hal.science/hal-03657988v1

Submitted on 12 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Violation of the fluctuation-dissipation theorem and effective temperatures in spin ice

Valentin Raban,¹,² Ludovic Berthier,²,³ and Peter C. W. Holdsworth¹

¹Université de Lyon, ENS de Lyon, CNRS, Laboratoire de Physique, F-69342 Lyon, France

²Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS, 34095 Montpellier, France

³Yusuf Hamied Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge CB2 1EW, United Kingdom

(Dated: April 13, 2022)

We present numerical tests of the fluctuation-dissipation theorem (FDT) in the dumbbell model of spin ice with parameters suitable for dysprosium titanate. The tests are made for local spin variables, magnetic monopole density, and energy. We are able to achieve local equilibrium in which the FDT is satisfied down to $T = 0.4$ K below which the system completely freezes. Non-equilibrium dynamics, together with violation of the FDT, are nonetheless observed following a thermal quench into the non-contractable monopole pair regime. Despite FDT violation, an approximate linear response regime allows for the identification of effective non-equilibrium temperatures which are different for each variable. The spin variable appears hotter than the heat reservoir, the monopole concentration responds with a lower effective temperature while the energy has a negative effective temperature. Results are discussed in the context of the monopole picture of spin ice and compared to the structure of FDT violations in other glassy materials. Prospectives for future experiments are reviewed.

# I. INTRODUCTION

Spin ice materials and models¹ share a number of characteristics with spin glasses² and other glassy materials³. They show slow dynamics at temperatures below the interaction energy scale⁴⁻¹², ergodicity breaking between field cooled (FC) and zero-field cooled (ZFC) protocols⁵ and, under certain conditions appear to show stretched exponential slowing down characteristic of supercooled liquids¹³. It is therefore interesting to borrow tools developed in the field of spin and structural glasses to analyse spin ice in order to understand better the similarities and differences between these different types of frustrated materials undergoing some kind of ergodicity breaking.

Glassy systems, once inside the ergodicity breaking time frame, typically violate the fluctuation-dissipation theorem¹⁴ (FDT) with spontaneous fluctuations and response functions no longer related to one another with a prefactor corresponding to the temperature of the thermal bath¹⁵⁻²⁰. Many glassy materials are seen to ‘age’ as they evolve through a hierarchical structure of metastable states²¹. Despite this, integrated response functions often remain proportional to their conjugate correlation functions¹⁵,¹⁶ allowing for the definition of an effective temperature¹⁷,¹⁸ replacing that of the thermal bath and associated with the explored non-equilibrium states. As a consequence the study of effective temperatures and of the so-called fluctuation-dissipation ratio has proved a rich and powerful tool to study glassy materials and their non-equilibrium aging dynamics at low temperatures¹⁷,¹⁹. FDT violations have been studied in many disordered materials with slow dynamics, such as spin glasses²²⁻²⁵, and structural glasses²⁶⁻³⁰. The same tools have been applied to scores of non-equilibrium dynamics, like coarsening³¹,³², non-equilibrium critical dynamics³³⁻³⁷, sheared complex fluids³⁸⁻⁴⁰, and active matter⁴¹,⁴². While FDT violations can of course be gener-

ically expected in non-equilibrium regimes, the specific structure of the measured FDT violations in glassy materials has helped distinguishing different types of aging dynamics¹⁷,¹⁹,²⁰.

In this paper we present a numerical study of the fluctuation dissipation theorem, its satisfaction and violation within the monopole picture of spin ice. We show that annealing below the FC-ZFC ergodicity threshold⁵, the FDT remains satisfied, a result that we argue is consistent with the absence of a hierarchical structure in phase space. However, on making a rapid quench in temperature, which traps a finite concentration of ‘non-contractable pairs’ of monopole quasi-particle excitations⁴³, we find distinct FDT violations, with different effective temperatures for energy, spin and monopole degrees of freedom. In particular, the energy shows a negative effective temperature reminiscent of kinetically constrained models of glass forming materials⁴⁴⁻⁴⁷, where aging dynamics is similarly controlled by localised defects with thermally activated dynamics.

We explicitly limit ourselves to the monopole approximation⁴⁸,⁴⁹, that is, to the dumbbell model⁵⁰ in which the magnetic moments of spin ice materials, localised on the nodes of the corner-sharing pyrochlore lattice are replaced by needles, extended to touch at the centres of the tetrahedra, forming a diamond lattice of vertices, see Fig. 1. The vertices carry monopole charges interacting via Coulomb’s law⁴⁹. The monopole picture has been hugely successful in providing a theoretical framework for such a complex frustrated system, giving a quantitative description of the thermodynamics and a good qualitative description of both the static and dynamic magnetic properties of spin ice materials⁵¹. To this picture we highlight two further simplifications which should be separated from the use of the monopole Hamiltonian. Firstly we limit ourselves to periodic boundary conditions, and secondly we model the real dynamics using
Metropolis dynamics. This choice corresponds, at the microscopic level, to modelling spin (or in our case needle) flips with a single tunnelling rate, independently of the local environment$^{52}$.

This simplest of starting points allows for the identification of the intrinsic mechanisms that can drive FDT violations. These should be contrasted with extrinsic mechanisms or corrections to this simplest of monopole pictures which, as a consequence, will be identifiable in future work. Such effects have already been identified in the modelling of observed slow dynamics at low temperature where the influence of the rapidly falling monopole density is strongly enhanced by the presence of defects, of open boundaries, and of the temperature dependence of the microscopic tunnelling process$^{10,11,53}$.

The rest of the paper is organised as follows. In Sec. II we review the dumbbell model and the numerical techniques used to study dynamical correlation and response functions. In Sec. III we review the FDT adapted to spin ice models. In Sec. IV we show numerical data revisiting the slow equilibrium dynamics with FDT satisfaction as the monopole concentration falls at low temperature. In Sec. V we turn to non-equilibrium dynamics in the non-contractable pair regime in which our numerical data shows FDT violations. Some conclusions and perspectives are given in Sec. VI.

## II. THE MONOPOLE PICTURE AND THE DUMBBELL MODEL

Spin ice materials$^{54,55}$ such as dysprosium or holmium titanate are Ising-like frustrated ferromagnets with magnetic degrees of freedom on the sites of the corner sharing pyrochlore lattice, with moments pointing along the body-centred diagonal axes that join at the centres of the tetrahedra, which themselves form a diamond lattice of vertices, see Fig. 1. Interactions, both exchange and dipolar, are on the one Kelvin scale and are captured by the dipolar spin ice model$^{56,57}$ (DSI). The lowest energy states, with spins on each tetrahedron satisfying the Bernal-Fowler ice-rules$^{58}$ of two spins in and two out in each tetrahedron, form a low-energy extensive band containing the Pauling entropy of approximately $\frac{1}{2} \ln(\frac{3}{2})$ per spin. The band width of these states is well below the energy scale of interactions, as the long-range part of the dipole interactions are almost self-screened throughout the band$^{56,59,60}$. The separation of energy scales allows for the interpretation of the magnetic moments as elements of an emergent magnetostatic field with the low-energy states constrained by a divergence-free condition$^{61,62}$. Excitations above the low-energy band are topological defects to the emergent field$^{48}$ which, being dressed by the real magnetic fields of the dipoles, carry a magnetic charge$^{49}$. These are the magnetic monopoles that have been much studied over the last decade.

The monopole approximation for spin ice is captured

<figure><img src="image_2.png" /><figcaption>FIG. 1. From spins to dumbbells. Left: Pyrochlore lattice of corner sharing tetrahedra. Tetrahedron centres form a diamond lattice. Blue and red spheres illustrate monopoles of charge $\pm Q$. Right: The point dipoles are extended to needles touching at the diamond lattice sites. The needles carry magnetic flux and are dumbbelles of charge $q = \pm m/a$. In a 2in-2out configuration (top) the vertex is charge neutral. A 3in-1out (3out-1in) configuration carries a monopole charge $Q = 2m/a$ ($-Q = -2m/a$) (centre). A 4in (4out) configuration carries a double monopole charge $2Q = 4m/a$ ($-2Q = -4m/a$) (bottom).</figcaption></figure>

by the dumbbell model$^{49,51,63}$ in which the point-like dipole moments of the DSI are replaced by magnetic needles which touch at the centres of the tetrahedra. The magnetic flux is channeled down the needles, so that they carry dumbbells of magnetic charge $\pm m/a$ where $m$ is the magnetic moment and $a$ the distance between tetrahedron centres. Configurations with two spins in and two out are by construction charge neutral and the total charge on site $i$ is $Q_i = 0, \pm Q, \pm 2Q$, with $Q = 2m/a$ the monopole charge. Also by construction in this model, the ensemble of Pauling states are exactly degenerate, so that the model violates the third law of thermodynamics. The Hamiltonian for excitations above the extensively degenerate ground state is

$$ \mathcal{H}_{\text{db}} = \frac{u}{2} \sum_{i \neq j} \left( \frac{a}{r_{ij}} \right) \hat{n}_i \hat{n}_j - \mu \sum_i \hat{n}_i^2, \quad (1) $$

where $\hat{n}_i = Q_i/Q = 0, \pm 1, \pm 2$ is a site occupation variable, $r_{ij}$ the distance between sites $i$ and $j$, $u = \frac{\mu_0 Q^2}{4\pi a}$ is the Coulomb energy scale, $\mu_0$ is the permeability of free space and $\mu < 0$ the chemical potential$^{52,64,65}$.

Throughout the paper we follow standard notation for spin ice simulations and refer to a dimensionless length $L$, measured in cubic units. Each cubic cell contains 16 spins (dumbbells) so that the number of tetrahedra (monopole sites) is $N_0 = 8L^3$. We index diamond lattice sites hosting the monopoles with $i, j$ and the pyrochlore lattice sites hosting the spins by $\kappa, \lambda$. Quantitative measures refer to the spin ice material $Dy_2Ti_2O_7$ (DTO) for which we take the diamond lattice constant to be $a = 4.36$ Å, the nearest neighbour spin distance $r_{nn} = \sqrt{\frac{2}{3}}a = 3.56$
Å, and cube length $a_c = \frac{4a}{\sqrt{3}} \approx 10.1$ Å (see Fig. 1)$^{66}$. We use the Kelvin energy scale, fixing the Boltzmann constant to unity. Taking an estimate for the magnetic moment $m = 9.87 \mu_B$$^{67}$ yields $u = 2.88$ K and $\mu = -4.35$ K.

Dynamics are simulated using a Metropolis Monte Carlo algorithm for dumbbell flips between their discrete orientations. Through the dumbbell flips the charges associated with the monopole quasi-particles move, are created or are destroyed in nearest neighbour pairs. The long range nature of the Coulomb interaction is dealt with using the Ewald summation method. Time steps are recorded as Monte Carlo moves per dumbbell (MCS). Spin configurations at each step can be recovered by replacing each dumbbell by the corresponding spin.

The emergent monopole picture can also be cast in the language of a lattice field which can be decomposed into ‘longitudinal’, ‘transverse’ and ‘harmonic’ components via a lattice Helmholtz decomposition$^{61,62,68}$. The monopoles can be reconstituted from the longitudinal part ($m$), with the leftover, defined by the constraints of fixed spin density, shared between the other two components. The divergence-free transverse component ($d$) has characteristic dipolar correlations$^{69,70}$. With periodic boundaries and in zero field, the harmonic term can be set to zero to an excellent approximation, so that for each microstate the spin at site $\kappa$ can be written

$$S_{\kappa} = S_{\kappa}^{m} + S_{\kappa}^{d}. \qquad (2)$$

Vector fields can be built at each site by multiplying the components by the local unit vectors connecting the di-amond lattice sites with appropriate sign convention$^{62}$. The closed loops of the transverse component provide the residual entropy as the temperature falls to zero while the energy (of the monopoles) is carried uniquely by the longitudinal component. In Fourier space, the two components fall parallel and perpendicular to the propagation vector $\vec{q}$ defined within the first Brillouin zone. This apparent fragmentation of the magnetic moments into two orthogonal parts generates a form of spin-charge separation in which magnetic monopoles and transverse spin components develop largely independent fluctuations and response functions.

### III. THE FLUCTUATION DISSIPATION THEOREM OUT OF EQUILIBRIUM

#### A. The fluctuation-dissipation theorem

The fluctuation dissipation theorem (FDT) relates the equilibrium fluctuations of an observable $A$ (in this case extensive), to the linear response to an applied field $f$ conjugate to $A$, such that the Hamiltonian at finite $f$ is $H(f) = H(0) - fA$. Starting an experiment at time $t = 0$ in equilibrium and applying a field perturbation at time $t_w \ge 0$, the variable $\langle A(t) \rangle$ is measured at time $t \ge t_w$, where the brackets stand for an ensemble average. This

setting provides a general form for the FDT:

$$\chi_A(t, t_w) = \frac{\partial \langle A(t) \rangle}{\partial f(t_w)} = \frac{1}{T} [C_A(t, t) - C_A(t, t_w)], \quad (3)$$

with $C_A(t, t_w) = \langle A(t)A(t_w) \rangle - \langle A(t) \rangle \langle A(t_w) \rangle$. The proportionality constant, $1/T$, between left and right hand functions illustrates the fact that the FDT is indeed explicitly derived in the conditions of thermodynamic equilibrium. In addition, in equilibrium the functions depend on a single time variable, the time difference, $t' = t - t_w$. Taking $t' \to \infty$ we find the static limit, $\chi_A = \frac{\partial \langle A \rangle}{\partial f} = \frac{1}{T} [\langle A^2 \rangle - \langle A \rangle^2]$.

Defining the normalised response and auto-correlation functions

$$\begin{aligned} \tilde{\chi}_A(t, t_w) &= \frac{\chi_A(t, t_w)}{C_A(t, t)}, \\ \tilde{C}_A(t, t_w) &= \frac{C_A(t, t_w)}{C_A(t, t)}, \end{aligned} \qquad (4)$$

the FDT in Eq. (3) takes the very simple form $T\tilde{\chi}_A = 1 - \tilde{C}_A$ and can be tested graphically by making parametric plots of $T\tilde{\chi}_A$ as a function of $1 - \tilde{C}_A$, with FDT satisfaction corresponding to a straight line of unit slope.

Starting at $t = 0$ from a non-equilibrium situation such as a thermal or field quench$^{71}$, the FDT may in principle be violated, requiring the insertion into Eq. (4) of a further parameter, $X_A(t, t_w)$, called the fluctuation-dissipation ratio, such that

$$T\tilde{\chi}_A = X_A(t, t_w)(1 - \tilde{C}_A), \qquad (5)$$

with $X_A(t, t_w) = 1$ at equilibrium. In principle, if one waits long enough, so that $t_w$ exceeds all relaxation times of the system, equilibrium behaviour and the FDT will be recovered with $X_A \to 1$. However, in glassy regimes the equilibrium relaxation times may exceed all possible observation times with the result that they appear permanently out of equilibrium, with corresponding violations of the FDT and $X_A \neq 1$.

In following such protocols for systems where time-translation invariance is not satisfied, either $t$ or $t_w$ can be varied as the control parameter, although varying $t_w$ can be extremely time consuming, as each measurement requires an independent procedure for each value of $t_w$. This problem is particularly acute when the fluctuation-dissipation $X_A(t, t_w)$ has a non-trivial time-dependence, since its mathematical definition via Eq. (5) requires using $t_w$ as the appropriate time variable$^{72}$. This problem can be fully circumvented in numerical simulations by using the no-field techniques$^{24,30,73}$ described in Appendix A, in which the linear response to the applied field can be extracted directly from derivatives of the Boltzmann weights in zero field. With this innovation, numerical simulations varying $t$ or $t_w$ become equivalent in computational effort as both correlation and response functions are evaluated in the same set of unperturbed simulations.
In certain model systems showing a thermodynamic glass transition, $X_A(t, t_w)$ can be shown to be a constant$^{16}$, different from unity, allowing for the definition of an effective temperature$^{18}$,

$$T_{\text{eff}} = \frac{T}{X_A}, \qquad (6)$$

associated with out of equilibrium fluctuations. No such definition can be made in more realistic models showing glassy behaviour, but $X_A(t, t_w)$ does appear approximately linear in many numerical simulations$^{27,30,38,39,74}$ and in some experiments$^{28,29}$.

This phenomenology has been an important development over the last three decades as it provides a test of possible universal behaviour of the non-equilibrium dynamics of systems with very long relaxation times. This has been largely triggered by the analytic study of a broad family of mean-field glassy models in their aging regime$^{15,16,75}$. A unique effective temperature shared by all physical observables is found in supercooled liquids$^{27,39,40}$, which are characterised by a clear separation of timescales with a unique slow relaxation mechanism. Kinetically constrained and non-mean-field models of glasses partly share this phenomenology$^{72}$, but reveal in addition the possibility of an absolute negative effective temperature when considering energy fluctuations and response functions$^{44-47}$. This striking result directly follows from the thermally activated nature of the microscopic aging motion in such systems, whereby increasing the temperature leads to an acceleration of the energy decrease, and thus to negative response functions$^{44}$.

By contrast, spin glasses are characterised, at least at mean-field level, by a more complex hierarchy of relaxation timescales$^{76}$, leading to a fluctuation-dissipation ratio which becomes a function rather than a simple number. This finding has direct, deep connections to the Parisi overlap distribution describing the spin glass order in equilibrium at low temperatures$^{15,77-79}$. Physically, one can associate an effective temperature to each relaxation timescale, and the emergence of multiple timescales suggests the possibility that distinct effective temperatures can also emerge.

Critical dynamics in pure ferromagnets also display remarkable universal violation of the FDT with the fluctuation-dissipation ratio taking a universal value uniquely dependent of the universality class$^{33-37}$. Coarsening and domain growth processes in ordered phases are characterised by an effective temperature which asymptotically diverges, or equivalently, by a vanishing fluctuation-dissipation ratio$^{31,32}$.

With few notable exceptions$^{25,80}$ most experiments dedicated to tests of FDT violations in aging glassy materials are not performed in the time domain but rather in Fourier space$^{28,29,81-84}$. In this case, the FDT relates the fluctuation spectral density $S(\omega)$ of a given observable at frequency $\omega$ to the out of phase linear susceptibility

$\chi''(\omega)$ at the same frequency:

$$S(\omega) = \frac{2k_B T}{\pi \omega} \chi''(\omega), \qquad (7)$$

and its generalisation introduces a frequency-dependent fluctuation-dissipation ratio as the ratio between the right and left hand sides of Eq. (7). In the case where the fluctuation-dissipation ratio is a simple number in the slow regime, then the time and frequency domain approaches yield similar results$^{18}$, with a linear relation between response and fluctuations with an effective temperature replacing the thermal bath temperature.

## B. Correlation and response functions in spin ice

We adapt the FDT formalism to three different observables in the dumbbell model of spin ice. The first is the local moment on each site, $S_\kappa$. This is an Ising-like variable whose conjugate field, $h_\kappa$, lies parallel or anti-parallel to the axis joining the two tetrahedron centres which the spin connects. The relevant correlation function

$$\tilde{C}_S(t, t_w) = \frac{\langle S_\kappa(t) S_\kappa(t_w) \rangle - \langle S_\kappa(t) \rangle \langle S_\kappa(t_w) \rangle}{\langle S_\kappa(t) S_\kappa(t) \rangle - \langle S_\kappa(t) \rangle \langle S_\kappa(t) \rangle}, \qquad (8)$$

is related to the local spin susceptibility, $\chi_S = \frac{\partial \langle S_\kappa \rangle}{\partial h_\kappa}$ and measurements are taken from a configurational average over the $N$ spins.

Secondly, we consider fluctuations of the configurational energy $E$ calculated from Eq. (1) through the response and correlation functions $\tilde{\chi}_E(t, t_w)$ and $\tilde{C}_E(t, t_w)$. In the language of a fluid in the grand ensemble, one might call this a number enthalpy as for each configuration $E = U_C - \mu N$ contains both the Coulomb energy, $U_C$ and the energy cost of creating $N = \sum_i |\hat{n}_i|$ magnetic monopoles. In principle one should also take into account the cost of creating double monopoles (tetrahedra with all spins pointing in or out) but for the temperatures used in this paper these can be safely neglected. The conjugate field for the energy is a dimensionless temperature change, $\alpha = \frac{\delta T}{T}$. The response function in the static limit is then closely related to the specific heat at constant chemical potential, $C_\mu$, which is equivalent to the magnetic specific heat in zero field, $\chi_E = \frac{\partial E}{\partial \alpha} = T C_\mu$.

Finally, we consider fluctuations in the monopole concentration, choosing as for the spins a local measure, $|\hat{n}_i|$, with conjugate field a local chemical potential, $\mu_i = \mu + \delta\mu_i$ and with the response and correlation functions $\tilde{\chi}_Q(t, t_w)$ and $\tilde{C}_Q(t, t_w)$. The response function, $\chi_Q = \frac{\partial |\hat{n}_i|}{\partial \mu_i}$, is closely related to the monopole compressibility, $\frac{\partial N}{\partial \mu}$.
# IV. DYNAMICS AT EQUILIBRIUM

## A. Low-temperature dynamics

The evolution of magnetic time scales extracted from ac susceptibility measurements⁵ on Dy₂Ti₂O₇ is largely captured by stochastic monopole dynamics in the dumbbell model⁶ and through spin dynamics in the DSI¹⁰,¹¹,⁵²,⁸⁵. On cooling, the measured time scale exhibits a plateau-like region between 10 K and 4 K before entering a regime of rapid change in which the characteristic time increases faster than an Arrhenius law, before leaving the experimental time window for temperatures below around 0.65 K⁵. The latter also corresponds to the ergodicity breaking temperature between field-cooled and zero-field cooled protocols for dc susceptibility measurements. The increasing time scale is directly related to the fall in monopole density, with the non-Arrhenius behaviour mostly due to the Coulomb interaction between monopoles and the increased screening length as temperature is reduced. As Coulomb screening falls to zero at low temperature, one might expect the time scale to approach an Arrhenius law, $\tau = \tau_0 \exp(-\mu/T)$, with $\mu = -4.35$ K⁵²,⁶⁴ and $\tau_0$ a microscopic timescale.

Experimental results differ in two details. Firstly, the characteristic time slows down beyond this limit⁴,⁵,⁷⁻¹⁰ and, secondly, rather than a single time scale, a spread of times appears at each temperature⁸⁶ and some measurements suggest stretched exponential behaviour¹⁰,¹². Possible corrections to the monopole picture to account for these differences include moving from the dumbbell model back to the spins of the DSI which reintroduces a finite band width for the Pauling states¹⁰,⁵²,⁸⁷, considering open rather than periodic boundaries which can introduce multiple time scales and stretched exponential decay, or the effects of defects or configuration-dependent microscopic tunnelling rates, both of which increase time scales at low temperature¹⁰,¹¹.

Here we show that finite-size effects, even in this simple model, also provide corrections to the above theoretical result and longer than expected time scales at low temperature. We calculate the spin-spin autocorrelation function, $\tilde{C}_S(t, 0) = \tilde{C}_S(t)$, with a configurational average made over the $N$ spins. As this is a local measure, it contains the diagonal terms of a magnetic correlation function only. Previous studies have shown that inclusion of off-diagonal terms makes no significant difference to the measured decay of correlations with time⁶,¹⁰, at least down to the magnetic ergodicity temperature, while use of the local function allows for much better statistics allowing us to extract data at lower temperature.

Data are shown in Fig. 2 for temperatures between 0.8 and 0.4 K, for a fixed system size, $L = 6$, with correlations averaged over 2000 initial equilibrium configurations. As can be seen from the logarithmic time axis, at each temperature the correlation function falls to zero at a single, well defined time scale $\tau$ and the decay is well-represented by an exponential decay, $C_S(t) \approx \exp(-(t/\tau))$.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Top: Normalised time autocorrelation function of the spin degrees of freedom at equilibrium. Temperatures are shown in Kelvin and the energy scale is set by parameters from Dy₂Ti₂O₇. Bottom: Correlation time for the spin degrees of freedom versus the inverse of the temperature for temperature between 0.8 and 0.4 K (points). The lines show Arrhenius laws with energy scales $|\mu|$ and $2|\mu|$ with $\mu = -4.35$ K.</figcaption></figure>

The value of $\tau$ increases rapidly with decreasing temperature and exceeds $10^7$ MCS for $T \le 0.4$ K. Given that the microscopic tunnelling rate is estimated to be in the millisecond range⁶, this would correspond to experimental time scales of several hours.

The evolution of $\tau$ with inverse temperature is also shown in Fig. 2. Linear behaviour would correspond to an Arrhenius law $\tau = \tau_0 \exp(\Delta E/T)$, with a slope giving the characteristic energy scale $\Delta E$. Shown also is Arrhenius behaviour for both $\Delta E = |\mu|$ and $\Delta E = 2|\mu|$. It is clear that, in this temperature range, the behaviour remains non-Arrhenius and the energy scale at the lowest temperature considerably exceeds that set by the monopole chemical potential. The reason is that, for $L = 6$ (1728 diamond lattice sites), on reducing the temperature, one rapidly enters a regime where the average number of monopoles in the system is less than two, so that for most configurations there are none. Putting $n = \frac{2}{8L^3} = \frac{4}{3} \exp(\mu/T)$, the density expected for the Coulomb fluid at low temperature⁶³, gives a temperature $T = 0.62$ K, which is also the threshold tempera-
ture of earlier studies. Below this temperature crossover, monopoles are confined by the periodic boundaries and the system must be excited through an energy scale of $2|\mu| = 8.7$ K to create monopoles that are then free to decorrelate the system. This is the crossover we are observing in Fig. 2. This finite-size effect is rather stylised to be directly relevant for experiment but it is consistent with the effect of adding defects to the DSI, which appear to generate a longer time scale in the decay of correlation functions¹⁰. It also illustrates how mosaics and grain boundaries could lead to non-universal long-time behaviour, as seen in the ensemble of experimental results.

Despite this change of regime to finite-size, or mosaic-driven long time scales, we are able to equilibrate down to temperatures somewhat below the experimental ergodic threshold, consistently with the most recent low temperature neutron scattering and specific heat measurements on Dy₂Ti₂O₇⁸⁸. This observation begs a rather subtil question concerning the emergence of a topological time scale characterising magnetic relaxation which is different from the local equilibrium time scale exposed here. At low temperature, a change in the magnetisation requires the flipping of an extended string of spins, giving ultimately a change in topological sector as the monopole concentration falls to zero⁸⁹. Although a detailed analysis is beyond the scope of the present paper, these results suggest the possibility of an ergodic time scale for such sector fluctuations that exceeds that for local equilibrium, in which local measures find their zero-field equilibrium values. Calculation of the response function for the local spin variable, rather than the bulk magnetisation, allows us to bypass this delicate question, which we leave for future work.

## B. The FDT at equilibrium

In Fig. 3, we show the parametric plot for the local spin correlation and response functions, $\tilde{\chi}_S$ and $\tilde{C}_S$ for $T = 1$ K where equilibration is easily achieved. The data is taken by varying total time $t$ with $t_w$ held fixed, with $L = 6$. The system is started at $t = 0$ in an equilibrium configuration and results are averaged over 5000 independent runs. The longest times correspond to $1 - \tilde{C}_S \approx 1$, so that spins are fully decorrelated. As expected, the FDT is satisfied with a high degree of accuracy over the whole spectrum of correlation and response. Lowering the temperature, the correlations times increase but as shown above, local equilibrium can be achieved down to $T = 0.4$ K, with corresponding satisfaction of the FDT. This is illustrated in Fig. 3, lower panel, where we show data for the energy $E$ at $T = 0.6$ K, temperature close to the magnetic ergodicity breaking and for which the FDT is also clearly satisfied.

However, the exponentially diverging time scale guarantees that equilibration, local or otherwise is excluded at lower temperatures. For example, our results suggest

that moving to $T = 0.2$ K would give a relaxation time scale, $\tau \sim \tau_0 \exp(40) \sim 10^7$ years for dysprosium titanate. Over the vast majority of this time scale a finite system would have a monopole population of zero, which ensures that a system near equilibrium would be completely frozen throughout any experimental time window. As a consequence, low temperature states with measurable dynamics must, by construction be far from equilibrium.

<figure><img src="image_4.png" /><figcaption>FIG. 3. FDT for spin ice in equilibrium. Top: Parametric plot of $T\tilde{\chi}_S$ vs. $1 - \tilde{C}_S$, for the spin-spin autocorrelation function and spin susceptibility for the dumbbell model for $T = 1$ K. The line shows the predicted FDT result. Bottom: Parametric plot of $T\tilde{\chi}_E$ vs. $1 - \tilde{C}_E$, for the energy autocorrelation function and susceptibility for the dumbbell model for $T = 0.6$ K. The line shows the predicted FDT result.</figcaption></figure>

## V. NON-EQUILIBRIUM DYNAMICS

### A. Non-contractable pairs

A class of non-equilibrium states showing measurable dynamics at ultra-low temperatures are those with a finite concentration of ‘non-contractable’ monopole pairs⁴³, which are illustrated in Fig. 4. Here, annihilating the nearest neighbour monopole pair would lead to an energy gain, $|2\mu| - u = 5.82$ K, so that far below this energy scale such defects would disappear at equilibrium. How-
<figure><img src="image_5.png" /><figcaption>FIG. 4. A non-contractable pair of charges (left). Flipping the central spin does not annihilate the charges but instead creates two double charges (right).</figcaption></figure>

ever, the constraints of spin ice are such that flipping the spin separating the monopole pair creates a pair of double monopoles rather than annihilating the charge. As a consequence, charge annihilation requires movement of the quasi-particles around an external path and passage over an energy barrier. As a result such a pair can be locked in a metastable state over a finite period of time after the system is quenched from a high temperature state (the opposite temperature protocol would not yield any interesting metastability). A detailed discussion of this process was provided in Ref. [43].

Here we confirm the key role played by the Coulomb interaction through a study of the dumbbell model. A minimal path for monopole annihilation is a hexagon of six tetrahedra with a maximum separation of third nearest neighbour on the diamond lattice. Using the parameters quoted here for $Dy_2Ti_2O_7$ the energy barrier, $\Delta E$ for the annihilation of an isolated non-contractable pair is around 1.4 K. Such a barrier would give an Arrhenius time scale, $\tau \sim \tau_0 \exp(14) \sim 10^6$ MCS at $T = 0.1$ K which suggests the existence of slow, but accessible dynamics down to this temperature scale.

Monopole-monopole interactions are the key to the pair formation. A nearest neighbour monopole interaction corresponds to second neighbour spin interactions, while the isotropic nature of the monopole picture comes explicitly from summing terms of the original DSI to infinite distance⁵⁶. Hence, even though truncating the Coulomb potential could speed up the numerics the long range interactions would still be playing a crucial role.

To access such states we performed thermal quenches⁴³ from an initial temperature $T_0 = 1$ K where the equilibrium monopole concentration is relatively high, to a target, much lower temperature $T$. In Fig. 5 we show the time evolution of $\langle n \rangle = \langle N/N \rangle$ and $\langle E \rangle$ with the time spent since the quench. Starting from the equilibrium monopole concentration at $T_0 = 1$ K, $n \sim 0.2$ per tetrahedron, $\langle n(t) \rangle$ falls rapidly, for all temperature. For $T = 0.3$ K and above, the monopole concentration falls towards its (very low) equilibrium value in a few hundreds of Monte Carlo steps. Below this crossover temperature, the quench is sufficiently violent to trap a finite concentration of non-contractable pairs, forming a concentration plateau whose lifetime increases rapidly with decreasing the quench temperature $T$. This behaviour is accurately paralleled in the energy function as, once the

<figure><img src="image_6.png" /><figcaption>FIG. 5. Monopole concentration (top) and energy (bottom) following a quench at $t = 0$ from $T_0 = 1$ K to a final temperature $T$. For a quench below 0.3 K, a plateau regime appears for both quantities.</figcaption></figure>

free particles have been annihilated, the energy depends, almost entirely on the small concentration of isolated and locally-confined charge pairs, see Fig. 5.

The apparition of two time scales following a quench means that spin ice undergoes a simple form of aging, whereby correlation functions depend on both times $t$ and $t_w$ following the quench and not just the time difference, $t - t_w$, as in equilibrium. In Fig. 6 we show the normalised local charge correlation function, $\tilde{C}_Q(t, t_w)$, following a quench to $T = 0.12$ K for different fixed values of $t_w$. In the initial regime in which free monopoles are annihilated, the correlation function depends strongly on $t_w$ but it becomes independent, for $t_w > 100$ Monte Carlo Steps, as one enters the regime of metastable pair confinement. In this regime the correlation function depends on $t-t_w$ only, as in equilibrium, although we stress that the system remains far from equilibrium. Similar behaviour is found for both energy and spin correlations, despite the fact that most spins spend the majority of their time completely frozen.

In Fig. 7 we show $\tilde{C}_Q(t, t_w)$ for different quench temperatures, for fixed $t_w = 2000$ MCS. As anticipated⁴³, the relaxations closely follow an exponential decay from which one can extract an associated time scale, $\tau_Q$. As shown in the lower panel, the temperature evolution of
<figure><img src="image_7.png" /><figcaption>FIG. 6. Normalised time autocorrelation functions of the local charge as a function of total time $t$ for different waiting time $t_w$ after a quench from $T = 1$ K to $T = 0.12$ K. The aging of the correlation is trivial in the sense that the correlation reaches its equilibrium behaviour very quickly (about a hundred MCS) compared to the equilibration time.</figcaption></figure>

$\tau_Q(T)$ closely follows an Arrhenius law, with and energy barrier $\Delta E_Q \approx 1.3$ K, close to but slightly lower than the estimated Coulomb energy barrier ($\approx 1.4$ K) for activated monopole hopping around a hexagon ring. This small difference could be explained by taking into account the interactions between the bound monopole pairs. Taking this to be a random dipole interaction of order

$$ \delta E \sim \frac{\mu_0 m^2}{\pi L^3}, \qquad (9) $$

gives an energy scale $\delta E \sim 0.1$ K. This scale, which is the correct order of magnitude could then be used as a parameter in a stochastic decay model¹³ that captures the small shift in time scales from the Coulomb estimate.

## B. FDT violations in the non-contractable pair regime

The existence of a non-equilibrium regime with measurable fluctuations allows us to test the FDT in a non-equilibrium environment for the three designated variables, $|\hat{n}_i|$, $S_\kappa$ and $E$. In each case data was collected both for fixed $t$ varying $t_w \le t$ and for fixed $t_w$ varying $t \ge t_w$. The resulting parametric plots are shown in Figs. 8 and 9. To help with interpretation, we remind the reader that a short time difference $t-t_w$ corresponds to a small value of the horizontal axis, $1-\hat{C}_A \sim 0$, and that the equilibrium FDT appears as a straight line of slope unity in this representation.

The FDT is clearly violated in all cases. Remarkably, in each case we observe that the fluctuation-dissipation ratio, $X(t, t_w)$, is approximately constant over the accessible range of correlations, allowing for the association of an effective non-equilibrium temperature, Eq. (6). However, this temperature appears different for each variable, which shows that the non-equilibrium problem does

<figure><img src="image_8.png" /><figcaption>FIG. 7. Top: Normalised time autocorrelation functions of the local charge as a function of total time $t$ for fixed waiting time $t_w = 200$ after a quench from $T = 1$ K. Bottom: Estimated relaxation times $\tau_Q$ from the data in the upper panel are well described by an Arrhenius law with energy scale $\Delta E \approx 1.3$ K.</figcaption></figure>

not simply translate to an effective equilibrium problem, with phenomenology reproduced by a minimal variational principle. On the contrary, this is a more complex situation in which each quantity has its own fluctuation amplitude without the constraints set by the second law of thermodynamics. This type of finding is consistent with results for glassy systems characterised either by well-separated time scales¹⁸,⁷⁵,⁷⁹, or systems with independent degrees of freedom⁹⁰,⁹¹.

Figure 8 shows data for charge fluctuations, calculated using both $t_w$ (upper) and $t$ (lower) as the variable parameter. The data are consistent for short times, giving approximately linear behaviour and a constant value for $X_Q$ around $X_Q \approx 2.3$. Hence the data in this region depends on the time difference $t-t_w$ only, as in equilibrium. For $1-\hat{C}_Q$ greater than 0.5 the two methods give different results, indicating a clear two time dependence. Mathematically the fluctuation-dissipation ratio is defined by using $t_w$ as the variable⁷², and this representation should therefore be preferred. Indeed, the data indicate a much simpler structure of the FDT plot in this representation with a well-defined effective temperature characterising the charge degrees of freedom over the entire time regime.

Fig. 9 upper panel shows data for the local spin vari-
<figure><img src="image_9.png" /></figure>

FIG. 8. Top: Parametric fluctuation-dissipation plot for the local charge density, $|\hat{n}_i|$ following a quench to $T = 0.1$ K (see Sec. III). The variable parameter is the waiting time $500 < t_w < t$ at fixed $t = 3.6 \times 10^5$ MCS. Results are for $L = 6$ and data are averaged over $4 \times 10^5$ initial configurations. The solid black line corresponds to equilibrium FDT with dotted lines representing constant $X$ values as shown. Bottom: As for upper panel with fixed $t_w = 500$ MCS and variable parameter $t_w < t < 8 \times 10^5$ MCS.

ables $S_κ$ with $t_w$ as the variable parameter. The small scale on both the horizontal and vertical axes reflects the inaccessible relaxation time scale for the spin degrees of freedom as most of the spins remain completely static over the simulation time scale. As a result only small corrections from $\tilde{C}_S = 1$ are observable. The parametric plot of response and correlations can nevertheless be approximated reasonably well over this correlation range by a constant fluctuation dissipation ratio, $X \approx 0.7$, with some very weak departure from linearity. Data collected varying $t$ for fixed $t_w$ (not shown) are very similar. We have also tested magnetic fluctuations through the total moment $\vec{M}(t)$. The data is much harder to collect for the reasons discussed above, but our results are compatible with those for the local spins suggesting that local and global spin degrees of freedom are characterised by similar violations of the FDT with an equal fluctuation-dissipation ratio, as also seen in pure Ising models at criticality$^{36}$.

A fluctuation ratio, $X < 1$ ($X > 1$) corresponds to an effective temperature that is higher (lower) than the equi-

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

FIG. 9. Top: Parametric fluctuation-dissipation plot for the local spin variable, $S_κ$ following a quench to $T = 0.1$ K (see Sec. III). The variable parameter is the waiting time $500 < t_w < t$ at fixed $t = 3.6 \times 10^5$ MCS. Results are for $L = 6$ and data are averaged over $4 \times 10^5$ initial configurations. The solid black line corresponds to equilibrium FDT with dotted line representing constant $X$ values as shown. Bottom: As for upper panel for the energy variable, $E$.

librium temperature. We find that the effective temperature for the charge variable $|\hat{n}_i|$ is lower than the equilibrium temperature, $T_{\text{eff}} = T/X_Q \approx 0.04$ K while that for the spin variable, $S_κ$ is higher with $T_{\text{eff}} = T/X_S \approx 0.14$ K for quenches performed at the bath temperature $T = 0.1$ K.

Given that the monopoles are objects made up of spin textures it seems surprising that the non-equilibrium response to them and to the spin degrees of freedom are different. The difference finds its origin in the emergent Helmholtz decomposition of the moments into longitudinal and transverse fragments. The two fragments interact only weakly through the constraint fixed by the total spin density$^{62}$, allowing an apparent separation of the magnetic charge made from the longitudinal component and the total spin.

In equilibrium, the monopoles are deconfined objects, free to move independently, interacting only via a Coulomb force which falls to zero at large distance. However, in the non-contractable pair regime, the monopoles are locally confined and can only diffuse via double moves in which the first separates the particles to a second
neighbour distance and the second restores the nearest
neighbour pair. The double move requires the crossing of
an energy barrier with a consequent drastic reduction in
monopole mobility. It is possible that the reduced effec-
tive temperature reflects this reduced mobility. Effective
temperatures that are lower than the bath temperature
have been observed in a few other systems as well³⁵,⁹².
However, in parallel with the reduced monopole mobility
there is an overall enhancement of spin fluctuations asso-
ciated with the enhanced monopole concentration, lead-
ing to an increased effective temperature for the complete
spin degrees of freedom. Such an increase in effective spin
temperature is compatible with the case of glasses and
spin glasses following a temperature quench²²,²⁵,²⁷,⁷⁴.

The response to non-equilibrium energy fluctuations is
at first sight even more surprising. The data shown in
Fig. 9 (lower panel) yield a negative slope, indicating a
negative effective temperature. We observe X<sub>E</sub> &asymp; -0.7,
and so T<sub>eff</sub> = T/X<sub>E</sub> &asymp; -0.14 K, again using t<sub>w</sub> as the ap-
propriate variable parameter. Varying t at fixed t<sub>w</sub> yields
qualitatively similar results. Negative temperatures have
also been observed in kinetically constrained models of
glasses showing slow dynamics<sup>44-46</sup>. Here, the negative
value for X<sub>E</sub> was associated with a dynamic evolution
via thermally activated processes involving localised de-
fects rather than by the evolution within a hierarchical
free energy landscape provided by a mean-field analy-
sis. This is indeed the case also here: non-contractable
pairs are essentially isolated, decaying via the thermally
activated process discussed above. Intuitively, increas-
ing the temperature (which is the field conjugate to the
energy) therefore increases the rate of relaxation of the
number of pairs. Consequently, after the initial decay fol-
lowing the quench, a linear increase of the applied tem-
perature at t<sub>w</sub> lowers the observed energy after a fixed
time t-t<sub>w</sub> as energy relaxes faster at higher temperature,
which then naturally gives a negative response to a posi-
tive impulse field. As a result, the fluctuation-dissipation
ratio between response and correlation becomes negative.
It is a non-trivial finding that the negative fluctuation-
dissipation ratio takes, to a good approximation, a unique
value, thus mimicking an equilibrium system with a neg-
ative absolute effective temperature.

VI. DISCUSSION

In this paper we have shown numerically that the
dumbbell model of spin ice can be forced into a non-
equilibrium regime in which the fluctuation-dissipation
theorem is violated. This is the non-contractable pair
regime⁴³ in which neutral pairs of magnetic monopoles
are locally confined by the constraints of spin ice. The
pairs relax via the passage over an energy barrier of order
one Kelvin, allowing for slow but measurable dynamics
down to 100 mK. We have measured FDT violations as-
sociated with three quantities: local spin variables, lo-
cal monopole occupation number, and total energy. De-

spite the expected FDT violation, in each case the re-
lation between fluctuations and response was found to
be approximately linear, allowing for the definition of ef-
fective non-equilibrium temperature. This very simple
structure of the FDT violation is not trivial at all and
could not have been anticipated on simple grounds. It
therefore seems remarkable that a unique quantity, the
fluctuation-dissipation ratio, is able to characterise the
time-dependent relation between response and correla-
tion functions in the far from equilibrium aging regime
of spin ice. The effective temperatures were found to be
different in the three cases. The spin variables appear to
be at a higher temperature than that of the heat reser-
voir, the monopole concentration responds at an appar-
ent lower temperature, while the energy has an effective
temperature which is negative.

The difference in response temperature for the spin
and monopole variables appears as a direct consequence
of the fractionalisation of dipolar elements into magnetic
monopoles⁴⁹. As the monopoles only eat up a small frac-
tion of the total spin density⁶¹, the monopoles and total
spin can behave quasi-independently, which is the case
here. The spin activity is enhanced with respect to equi-
librium at such low temperatures of approximately 100
mK which reflects the higher effective temperature but
the monopoles are locally confined into pairs and have
considerably reduced mobility. It is possible that this is
reflected in the reduced effective temperature. In aging
and sheared supercooled liquids characterised by a single
glassy timescale, a unique value of the effective temper-
ature is observed²⁷,³⁹, but in that case all observables
are dynamically strongly coupled. Clearly, spin ice has a
richer physical behaviour with two independent variables.

In kinetic models of glass-forming materials, the neg-
ative energy response to an increase in the temperature
is a direct result of the thermally activated relaxation
over energy barriers. In these models⁴⁴⁻⁴⁷, microscopic or
mesoscopic objects indeed relax independently via local
energy barriers. As a consequence, on heating, the en-
ergy is relaxed more easily towards its lower, equilibrium
value via the faster motion of the defects, thus produc-
ing negative response functions and negative fluctuation-
dissipation ratio. Clearly, the bound pair regime of
spin ice is another example of this physics as the bound
pairs are essentially independent, two-body metastable
states. The similarity of the results for the energy de-
gree of freedom shows that the analogy extends to the
out of equilibrium regime. It would be interesting to see
whether the measured value of the fluctuation-dissipation
ratio, X<sub>E</sub> &asymp; -0.7 can be theoretically rationalised in the
framework of the diffusion-annihilation process proposed
in Ref.⁴³, following the type of field-theoretical analysis
done for kinetic glass models⁴⁴⁻⁴⁶.

It would also be interesting to study FDT relations in
artificial spin ice⁹³⁻⁹⁶ and colloidal ice systems⁹⁷. These
mesoscopic systems have the advantage that monopoles
are directly observable in real space and that the ac-
curacy of the monopole picture can be tuned to some
extent⁹⁸,⁹⁹ In these experimental systems the initial paradigm is slightly different in that attaining equilibrium, real or effective¹⁰⁰,¹⁰¹ for any situation is a formidable challenge so that an experiment in which the FDT were satisfied would already be a very interesting result. From here, performing quenches into the non-contractable pair regime with observation of FDT violations would seem a feasible project.

The above analogies with glass-forming materials and kinetic glass models therefore confirm that despite FDT violations, the non-contractable pair regime does not have the kind of hierarchical energy landscape structure that characterises spin glasses⁴³. The non-contractable pairs do themselves carry a magnetic dipole moment so that there is a collective dipolar interaction between them. It is therefore conceivable that they could show dipolar glass behaviour on a much lower temperature scale, but the separation of energy scales here is such that a regime of this kind is beyond the scope of our work.

The bound pairs are formed following a quench from 1 K where the monopole concentration is still high, to target temperatures below 0.3 K. Above this quench temperature the pairs fail to form and the concentration falls rapidly to the equilibrium value as if the system had been annealed. For a finite size system in this temperature range most microstates of the equilibrium measure have zero monopoles so that spin fluctuations essentially stop here. On annealing we were able to equilibrate down to around 0.4 K and the FDT was satisfied down to this limit, consistently with recent specific heat measurements⁸⁸ and with preliminary magnetic noise measurement¹⁰². Hence, for the chosen variables, following an annealing protocol, we find that the system progresses from an equilibrium to a frozen state without an intermediate regime showing non-equilibrium dynamics.

There is of course the potential for the development of hierarchical dynamics through the existence of closed loops of spin flips moving the system from one constrained spin ice configuration to another through the band of Pauling states. However, in the dumbbell model, like in the nearest-neighbour spin ice model, such states are degenerate so that any hierarchy must be entropy driven. The time and length scales necessary to address such subtle questions are at present out of reach. Moving back to the dipolar spin ice model lifts this degeneracy, opening up the possibility of a loop glass and hierarchical behaviour on the energy scale of this band width; that is, around 150 mK⁶⁰. Re-opening the bandwidth takes us beyond the monopole picture¹³,⁵⁶,⁵⁹ and into a realm in which defects and realistic boundaries will also be essential ingredients¹⁰. That is, glassy behaviour could be possible on this small energy scale but our results suggest that it will be extrinsic, depending on corrections to the monopole picture and to the translational symmetry of the pyrochlore lattice.

In preliminary experiments¹⁰² measuring magnetic noise and response in the frequency domain, the FDT

is confirmed over the full frequency spectrum down to 0.6 K. There is partial confirmation down to 0.4 K but at this lowest temperature, time scales are so long that measurements are limited to the relatively high frequency, short time scale regime. These preliminary measurements therefore confirm the results from the monopole picture that exclude FDT violations in this temperature range. Corrections to this picture, although of great interest will be subtle and present considerable experimental challenges.

Far from equilibrium states in spin ice materials, with high monopole density are attainable through field assisted monopole avalanche quench protocols⁷¹,¹⁰³,¹⁰⁴ with quenches going well below the required 300 mK. Magnetic noise measurement below this temperature range is a challenging experimental problem but our results suggest that the effort could be worth it as it will offer a robust regime of non-equilibrium response with measurable FDT violations. The origin of these violations can be found in the simplest monopole description of spin ice and they should resist the many corrections to this picture offered here and in the wider literature.

## ACKNOWLEDGMENTS

It is a pleasure to thank Elsa Lhotel and Claudio Castelnovo for fruitful discussions and for collaborations on related topics. This work was supported by ANR, France, Grant No.ANR-19-CE30-0040 and by a Grant from the Simons Foundation (Grant No. 454933, L.B.).

## Appendix A: Computation of response functions without introducing a field

The usual way to determine the response function of a given observable A is to run two different kinds of simulations: the first one with the Hamiltonian without perturbation, $H_0$, and the second one with the perturbed hamiltonian, $H_0 - f A$ which defines the field f conjugate to the observable A. We wish to measure the time-integrated response function, and so the time evolution of the field f is

$$f = \begin{cases} 0 & \text{if } t < t_w \\ \varepsilon & \text{if } t \ge t_w \end{cases} \qquad (A1)$$

In both kinds of simulations, the time evolution of $A(t)$ is computed, and then an ensemble average is done over a series of independent initial configurations. The response function is eventually obtained as the finite derivative

$$\chi(t, t_w) = \frac{\partial \langle A(t) \rangle_{\varepsilon}}{\partial f} \underset{\varepsilon \to 0}{\approx} \frac{1}{\varepsilon} \left\{ \langle A(t) \rangle_{\varepsilon} - \langle A(t) \rangle_{0} \right\}_{t \ge t_w} \qquad (A2)$$

Numerically, $\varepsilon$ has to be small so that the finite derivative is a good approximation to the true derivative (linear
response regime), but it can not be too small, otherwise the signal would be smaller than the statistical noise.

Another point to mention is the fact that each $t_w$ requires two sets of simulations (with and without the conjugated field), which ultimately results in a heavy computational load. This led us to implement instead a no-field method to obtained the directly response functions from unperturbed simulations$^{24,30,73}$.

If $P(\{\sigma\} \to \{\sigma'\})$ is the probability to accept a Monte Carlo move from the configuration $\{\sigma\}$ to the configuration $\{\sigma'\}$, then

$$ \langle A(t) \rangle = \sum_{\text{ic}} A(t) P(\{\sigma(t_w)\} \to \{\sigma(t)\}) \quad (\text{A3}) $$

where the sum runs over independent initial configurations. The integrated linear response is thus given by

$$ \chi(t, t_w) = \sum_{\text{ic}} A(t) \frac{\partial}{\partial f} P(\{\sigma(t_w)\} \to \{\sigma(t)\}), \quad (\text{A4}) $$

where the derivative with respect to the field can be evaluated exactly.

Using a markov chain of single spin flip with Metropolis rates, one can express the total probability as a product over each attempted Monte Carlo move,

$$ P(\{\sigma(t_w)\} \to \{\sigma(t)\}) = \prod_{t'=t_w}^{t} W(t' \to t' + 1) \quad (\text{A5}) $$

with the probabilities $W(t \to t + 1)$ directly following from the Metropolis acceptance rate, namely:

$$ \begin{cases} 0 & \text{if } \sigma(t+1) \text{ can not be reached from } \sigma(t) \\ 1 & \text{if } \Delta\mathcal{H} < 0 \\ e^{-\beta \Delta\mathcal{H}} & \text{if } \Delta\mathcal{H} \ge 0 \text{ and the spin flip is accepted} \\ 1 - e^{-\beta \Delta\mathcal{H}} & \text{if } \Delta\mathcal{H} \ge 0 \text{ and the spin flip is rejected} \end{cases} $$

One can then express the derivative in Eq. (A4) to get

$$ \frac{\partial}{\partial f} P(\{\sigma(t_w)\} \to \{\sigma(t)\}) = P(\{\sigma(t_w)\} \to \{\sigma(t)\}) \sum_{t'=t_w}^{t} \frac{\partial}{\partial f} \ln W(t' \to t' + 1) \quad (\text{A6}) $$

so that finally

$$ \begin{aligned} \chi(t, t_w) &= \sum_{\text{ic}} A(t) \left( \sum_{t'=t_w}^{t} \frac{\partial}{\partial f} \ln W(t' \to t' + 1) \right) \times \\ & \quad \times P(\{\sigma(t_w)\} \to \{\sigma(t)\}) \\ & \equiv \langle A(t) \mathcal{R}(t_w, t) \rangle \end{aligned} \quad (\text{A7}) $$

which defines the quantity

$$ \mathcal{R}(t_w, t) = \sum_{t'=t_w}^{t} \frac{\partial}{\partial f} \ln W(t' \to t' + 1). \quad (\text{A8}) $$

If the observable is $A$, and its conjugate field is $f$ then $-\beta\mathcal{H} \to -\beta\mathcal{H}+\beta fA$, then one can express the derivatives explicitly

$$ \frac{\partial}{\partial f} \ln W(t' \to t' + 1) = $$

$$ \left\{ \begin{array}{ll} +\beta \Delta A & \text{if } \Delta\mathcal{H} \ge 0 \text{ and the spin flip is accepted} \\ -\frac{\beta \Delta A}{e^{\beta \Delta\mathcal{H}} - 1} & \text{if } \Delta\mathcal{H} \ge 0 \text{ and the spin flip is rejected} \\ 0 & \text{otherwise.} \end{array} \right. $$

The key point of the method is the evaluation of integrated response functions in Eq. (A7) which involves the product of two physical observables that can be evaluated explicitly during the course of the simulation followed by an ensemble average with no perturbing field. Therefore, linear response is valid by construction and a single set of simulations is needed as the perturbing field is never applied. In addition, several physical observables can be analysed in the course of a single simulations. The only drawback is that the statistical noise affecting Eq. (A7) can become quite large, especially if long time scales are needed and this may imply using a large number of initial configurations.

<footnote n="1">1 L. D. C. Jaubert and M. Udagawa, Spin Ice (Springer, 2021).</footnote>

<footnote n="2">2 Marc Mézard, Giorgio Parisi, and Miguel Angel Vivasoro, Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications, Vol. 9 (World Scientific Publishing Company, 1987).</footnote>

<footnote n="3">3 Ludovic Berthier and Giulio Biroli, “Theoretical perspective on the glass transition and amorphous materials,” Rev. Mod. Phys. 83, 587–645 (2011).</footnote>

<footnote n="4">4 K Matsuhira, Y Hinatsu, and T Sakakibara, “Novel dynamical magnetic properties in the spin ice compound dy2ti2o7,” Journal of Physics: Condensed Matter 13, L737–L746 (2001).</footnote>

<footnote n="5">5 J. Snyder, B.G. Ueland, J.S. Slusky, H. Karunadasa, R.J. Cava, and P. Schiffer, “Low-temperature spin freezing in the Dy2Ti2O7 spin ice,” Physical Review B 69, 064414 (2004).</footnote>

<footnote n="6">6 L. D. C. Jaubert and P. C. W. Holdsworth, “Signature</footnote>
of magnetic monopole and dirac string dynamics in spin ice,” Nature Physics 5, 258–261 (2009).

7 J. A. Quilliam, L. R. Yaraskavitch, H. A. Dabkowska, B. D. Gaulin, and J. B. Kycia, “Dynamics of the magnetic susceptibility deep in the coulomb phase of the dipolar spin ice material $ho_{2}ti_{2}O_{7}$,” Phys. Rev. B 83, 094424 (2011).

8 Kazuyuki Matsuhira, Carley Paulsen, Elsa Lhotel, Hiroi Sekine, Zenji Hiroi, and Seishi Takagi, “Spin Dynamics at Very Low Temperature in Spin Ice $Dy_{2}Ti_{2}O_{7}$,” Journal of the Physical Society of Japan 80, 123711 (2011).

9 L R Yaraskavitch, H M Revell, S Meng, K A Ross, H M L Noad, H A Dabkowska, B D Gaulin, and J B Kycia, “Spin dynamics in the frozen state of the dipolar spin ice material $Dy_{2}Ti_{2}O_{7}$,” Phys. Rev. B 85, 20410 (2012).

10 H. M. Revell, L. R. Yaraskavitch, J. D. Mason, H. M. L. Ross, K. A. AU Noad, H. A. Dabkowska, B. D. Gaulin, P. Henelius, and J. B. Kycia, “Evidence of impurity and boundary effects on magnetic monopole dynamics in spin ice,” Nat. Phys. 9, 34–37 (2012).

11 Hiroshi Takatsu, Kazuki Goto, Hiromi Otsuka, Ryuji Higashinaka, Kazuyuki Matsubayashi, Yoshiya Uwatoko, and Hiroaki Kadowaki, “Ac susceptibility of the dipolar spin ice $dy_{2}ti_{2}o_{7}$: Experiments and monte carlo simulations,” Journal of the Physical Society of Japan 82, 104710 (2013).

12 Azar B. Eyvazov, Ritika Dusad, Timothy J. S. Munsie, Hanna A. Dabkowska, Graeme M. Luke, Ethan R. Kassner, J. C. Seamus Davis, and Anna Eyal, “Common glass-forming spin-liquid state in the pyrochlore magnets $dy_{2}ti_{2}o_{7}$ and $ho_{2}ti_{2}o_{7}$,” Phys. Rev. B 98, 214430 (2018).

13 Ethan R Kassner, Azar B Eyvazov, Benjamin Pichler, Timothy JS Munsie, Hanna A Dabkowska, Graeme M Luke, and JC Séamus Davis, “Supercooled spin liquid state in the frustrated pyrochlore $dy_{2}ti_{2}o_{7}$,” Proceedings of the National Academy of Sciences 112, 8549–8554 (2015).

14 David Chandler, “Introduction to modern statistical,” Mechanics. Oxford University Press, Oxford, UK 5 (1987).

15 Jean-Philippe Bouchaud, Leticia F Cugliandolo, Jorge Kurchan, and Marc Mezard, “Out of equilibrium dynamics in spin-glasses and other glassy systems,” Spin glasses and random fields, 161–223 (1998).

16 L. F. Cugliandolo and J. Kurchan, “Analytical solution of the off-equilibrium dynamics of a long-range spin-glass model,” Phys. Rev. Lett. 71, 173–176 (1993).

17 Leticia F Cugliandolo, “The effective temperature,” Journal of Physics A: Mathematical and Theoretical 44, 483001 (2011).

18 Leticia F. Cugliandolo, Jorge Kurchan, and Luca Peliti, “Energy flow, partial equilibration, and effective temperatures in systems with slow dynamics,” Phys. Rev. E 55, 3898–3914 (1997).

19 Andrea Crisanti and Felix Ritort, “Violation of the fluctuation-dissipation theorem in glassy systems: basic notions and the numerical evidence,” Journal of Physics A: Mathematical and General 36, R181 (2003).

20 Jorge Kurchan, “In and out of equilibrium,” Nature 433, 222–225 (2005).

21 Leendert Cornelis Elisa Struik, “Physical aging in plastics and other glassy materials,” Polymer Engineering & Science 17, 165–173 (1977).

22 Enzo Marinari, Giorgio Parisi, Federico Ricci-Tersenghi, and Juan J Ruiz-Lorenzo, “Violation of the fluctuation-dissipation theorem in finite-dimensional spin glasses,” Journal of Physics A: Mathematical and General 31, 2611 (1998).

23 Alain Barrat and Ludovic Berthier, “Real-space application of the mean-field description of spin-glass dynamics,” Phys. Rev. Lett. 87, 087204 (2001).

24 F. Ricci-Tersenghi, “Measuring the fluctuation-dissipation ratio in glassy systems with no perturbing field,” Phys. Rev. E 68, 065104 (2003).

25 D. Hérisson and M. Ocio, “Fluctuation-dissipation ratio of a spin glass in the aging regime,” Phys. Rev. Lett. 88, 257202 (2002).

26 Giorgio Parisi, “Off-equilibrium fluctuation-dissipation relation in fragile glasses,” Phys. Rev. Lett. 79, 3660–3663 (1997).

27 J-L Barrat and Walter Kob, “Fluctuation-dissipation ratio in an aging lennard-jones glass,” EPL (Europhysics Letters) 46, 637 (1999).

28 Tomás S. Grigera and N. E. Israeloff, “Observation of fluctuation-dissipation-theorem violations in a structural glass,” Phys. Rev. Lett. 83, 5038–5041 (1999).

29 Ludovic Bellon, Sergio Ciliberto, and Claude Laroche, “Violation of the fluctuation-dissipation relation during the formation of a colloidal glass,” EPL (Europhysics Letters) 53, 511 (2001).

30 Ludovic Berthier, “Efficient measurement of linear susceptibilities in molecular simulations: Application to aging supercooled liquids,” Phys. Rev. Lett. 98, 220601 (2007).

31 A. Barrat, “Monte carlo simulations of the violation of the fluctuation-dissipation theorem in domain growth processes,” Phys. Rev. E 57, 3629–3632 (1998).

32 Ludovic Berthier, Jean-Louis Barrat, and Jorge Kurchan, “Response function of coarsening systems,” The European Physical Journal B-Condensed Matter and Complex Systems 11, 635–641 (1999).

33 C Godreche and JM Luck, “Response of non-equilibrium systems at criticality: exact results for the glauber-ising chain,” Journal of Physics A: Mathematical and General 33, 1151 (2000).

34 Malte Henkel, Michel Pleimling, Claude Godrèche, and Jean-Marc Luck, “Aging, phase ordering, and conformal invariance,” Phys. Rev. Lett. 87, 265701 (2001).

35 Ludovic Berthier, Peter CW Holdsworth, and Mauro Sellitto, “Nonequilibrium critical dynamics of the two-dimensional xy model,” Journal of Physics A: Mathematical and General 34, 1805 (2001).

36 Peter Mayer, Ludovic Berthier, Juan P. Garrahan, and Peter Sollich, “Fluctuation-dissipation relations in the nonequilibrium critical dynamics of ising models,” Phys. Rev. E 68, 016116 (2003).

37 Pasquale Calabrese and Andrea Gambassi, “Ageing properties of critical systems,” Journal of Physics A: Mathematical and General 38, R133 (2005).

38 Jean-Louis Barrat and Ludovic Berthier, “Fluctuation-dissipation relation in a sheared fluid,” Phys. Rev. E 63, 012503 (2000).

39 Ludovic Berthier and Jean-Louis Barrat, “Nonequilibrium dynamics and fluctuation-dissipation relation in a sheared fluid,” The Journal of Chemical Physics 116, 6228–6242 (2002).

40 Ludovic Berthier and Jean-Louis Barrat, “Shearing a glassy material: Numerical tests of nonequilibrium mode-
coupling approaches and experimental proposals,” Phys. Rev. Lett. 89, 095702 (2002).

41 Davide Loi, Stefano Mossa, and Leticia F. Cugliandolo, “Effective temperature of active matter,” Phys. Rev. E 77, 051111 (2008).

42 Demian Levis and Ludovic Berthier, “From single-particle to collective effective temperatures in an active fluid of self-propelled particles,” EPL (Europhysics Letters) 111, 60006 (2015).

43 C. Castelnovo, R. Moessner, and S. L. Sondhi, “Thermal quenches in spin ice,” Phys. Rev. Lett. 104, 107201 (2010).

44 Peter Mayer, Sébastien Léonard, Ludovic Berthier, Juan P. Garrahan, and Peter Sollich, “Activated aging dynamics and negative fluctuation-dissipation ratios,” Phys. Rev. Lett. 96, 030602 (2006).

45 Sébastien Léonard, Peter Mayer, Peter Sollich, Ludovic Berthier, and Juan P Garrahan, “Non-equilibrium dynamics of spin facilitated glass models,” Journal of statistical mechanics: theory and experiment 2007, P07017 (2007).

46 Robert L Jack, Ludovic Berthier, and Juan P Garrahan, “Fluctuation-dissipation relations in plaquette spin systems with multi-stage relaxation,” Journal of Statistical Mechanics: Theory and Experiment 2006, P12005 (2006).

47 A. Garriga, I. Pagonabarraga, and F. Ritort, “Negative fluctuation-dissipation ratios in the backgammon model,” Phys. Rev. E 79, 041122 (2009).

48 I. A. Ryzhkin, “Magnetic relaxation in rare-earth oxide pyrochlores,” Journal of Experimental and Theoretical Physics 101, 481–486 (2005).

49 C Castelnovo, R Moessner, and SL Sondhi, “Magnetic monopoles in spin ice.” Nature 451, 42–45 (2008).

50 G. Möller and R. Moessner, “Artificial square ice and related dipolar nanoarrays,” Phys. Rev. Lett. 96, 237202 (2006).

51 C Castelnovo and PCW Holdsworth, “Modelling of classical spin ice: Coulomb gas description of thermodynamic and dynamic properties,” in Spin Ice (Springer, 2021) pp. 143–188.

52 L D C Jaubert and P C W Holdsworth, “Magnetic monopole dynamics in spin ice,” Journal of Physics: Condensed Matter 23, 164222 (2011).

53 Bruno Tomasello, Claudio Castelnovo, Roderich Moessner, and Jorge Quintanilla, “Correlated quantum tunneling of monopoles in spin ice,” Phys. Rev. Lett. 123, 067204 (2019).

54 MJ Harris, ST Bramwell, DF McMorrow, Th Zeiske, and KW Godfrey, “Geometrical frustration in the ferromagnetic pyrochlore ho 2 ti 2 o 7,” Physical Review Letters 79, 2554 (1997).

55 Steven T Bramwell and Michel JP Gingras, “Spin ice state in frustrated magnetic pyrochlore materials,” Science 294, 1495–1501 (2001).

56 Byron C. den Hertog and Michel J. P. Gingras, “Dipolar interactions and origin of spin ice in ising pyrochlore magnets,” Phys. Rev. Lett. 84, 3430–3433 (2000).

57 Roger G Melko and Michel JP Gingras, “Monte carlo studies of the dipolar spin ice model,” Journal of Physics: Condensed Matter 16, R1277 (2004).

58 J. D. Bernal and R. H. Fowler, “A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions,” The Journal of Chemical Physics 1, 515–548 (1933).

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

60 Roger G. Melko, Byron C. den Hertog, and Michel J. P. Gingras, “Long-range order at low temperatures in dipolar spin ice,” Phys. Rev. Lett. 87, 067203 (2001).

61 M. E. Brooks-Bartlett, S. T. Banks, L. D. C. Jaubert, A. Harman-Clarke, and P. C. W. Holdsworth, “Magnetic-moment fragmentation and monopole crystallization,” Phys. Rev. X 4, 011007 (2014).

62 Elsa Lhotel, Ludovic DC Jaubert, and Peter CW Holdsworth, “Fragmentation in frustrated magnets: A review,” Journal of Low Temperature Physics 1–28 (2020).

63 V. Kaiser, J. Bloxsom, L. Bovo, S. T. Bramwell, P. C. W. Holdsworth, and R. Moessner, “Emergent electrochemistry in spin ice: Debye-hückel theory and beyond,” Phys. Rev. B 98, 144413 (2018).

64 C. Castelnovo, R. Moessner, and S. L. Sondhi, “Debye-Hückel theory for spin ice at low temperature,” Physical Review B 84, 144435 (2011).

65 V. Raban, C. T. Suen, L. Berthier, and P. C. W. Holdsworth, “Multiple symmetry sustaining phase transitions in spin ice,” Phys. Rev. B 99, 224425 (2019).

66 H. D. Zhou, J. G. Cheng, A. M. Hallas, C. R. Wiebe, G. Li, L. Balicas, J. S. Zhou, J. B. Goodenough, J. S. Gardner, and E. S. Choi, “Chemical pressure effects on pyrochlore spin ice,” Phys. Rev. Lett. 108, 207206 (2012).

67 Taras Yavors’kii, Tom Fennell, Michel J. P. Gingras, and Steven T. Bramwell, “$dy_2ti_{2o7}$,” Phys. Rev. Lett. 101, 037204 (2008).

68 Steven T Bramwell, “Harmonic phase in polar liquids and spin ice,” Nature communications 8, 1–9 (2017).

69 David A. Huse, Werner Krauth, R. Moessner, and S. L. Sondhi, “Coulomb and liquid dimer models in three dimensions,” Phys. Rev. Lett. 91, 167004 (2003).

70 S V Isakov, K Gregor, R Moessner, and S L Sondhi, “Dipolar spin correlations in classical pyrochlore magnets,” Physical Review Letters 93, 167204 (2004).

71 C. Paulsen, M. J. Jackson, E. Lhotel, B. Canals, D. Prabhakaran, K. Matsuhira, S. R. Giblin, and S. T. Bramwell, “Far-from-equilibrium monopole dynamics in spin ice,” Nature Physics 10, 135–139 (2014).

72 Peter Sollich, Suzanne Fielding, and Peter Mayer, “Fluctuation-dissipation relations and effective temperatures in simple non-mean field systems,” Journal of Physics: Condensed Matter 14, 1683 (2002).

73 Christophe Chatelain, “A far-from-equilibrium fluctuation-dissipation relation for an ising-glauber-like model,” Journal of Physics A: Mathematical and General 36, 10739 (2003).

74 Francesco Sciortino and Piero Tartaglia, “Extension of the fluctuation-dissipation theorem to the physical aging of a model glass-forming liquid,” Phys. Rev. Lett. 86, 107–110 (2001).

75 Leticia Cugliandolo and Jorge Kurchan, “Thermal properties of slow dynamics,” Physica A: Statistical Mechanics and its Applications 263, 242–251 (1999).

76 Leticia F Cugliandolo and Jorge Kurchan, “On the out-of-equilibrium relaxation of the sherrington-kirkpatrick model,” Journal of Physics A: Mathematical and General 27, 5749 (1994).

77 Silvio Franz, Marc Mézard, Giorgio Parisi, and Luca Peliti, “Measuring equilibrium properties in aging sys-
tems,” Phys. Rev. Lett. 81, 1758–1761 (1998).

78 Ludovic Berthier, Jean-Louis Barrat, and Jorge Kurchan, “Dynamic ultrametricity in spin glasses,” Phys. Rev. E 63, 016105 (2000).

79 Jorge Kurchan, “Time-reparametrization invariances, multithermalization and the parisi scheme,” arXiv preprint arXiv:2101.12702 (2021).

80 Didier Hérisson and Miguel Ocio, “Off-equilibrium fluctuation-dissipation relation in a spin glass,” The European Physical Journal B-Condensed Matter and Complex Systems 40, 283–294 (2004).

81 Ph Refregier, M Alba, J Hammann, and M Ocio, “Dynamic behaviour of the insulating spin glass csnifef6,” Journal of Physics C: Solid State Physics 20, 5545 (1987).

82 Ludovic Bellon and Sergio Ciliberto, “Experimental study of the fluctuation dissipation relation during an aging process,” Physica D: Nonlinear Phenomena 168, 325–335 (2002).

83 Lionel Buisson, Ludovic Bellon, and Sergio Ciliberto, “Intermittency in ageing,” Journal of Physics: Condensed Matter 15, S1163 (2003).

84 Bérengère Abou and Francois Gallet, “Probing a nonequilibrium einstein relation in an aging colloidal glass,” Phys. Rev. Lett. 93, 160603 (2004).

85 Ritika Dusad, Franziska KK Kirschner, Jesse C Hoke, Benjamin R Roberts, Anna Eyal, Felix Flicker, Graeme M Luke, Stephen J Blundell, and JC Séamus Davis, “Magnetic monopole noise,” Nature 571, 234–239 (2019).

86 L. Bovo, J.A. Bloxsom, D. Prabhakaran, G. Aeppli, and S.T. Bramwell, “Brownian motion and quantum dynamics of magnetic monopoles in spin ice,” Nature Communications 4, 1535 (2013).

87 Stephen J. Blundell, “Monopoles, magnetic density, and the stray field from spin ice,” Physical Review Letters 108, 147601 (2012).

88 S. R. Giblin, M. Twengström, L. Bovo, M. Ruminy, M. Bartkowiak, P. Manuel, J. C. Andresen, D. Prabhakaran, G. Balakrishnan, E. Pomjakushina, C. Paulsen, E. Lhotel, L. Keller, M. Frontzek, S. C. Capelli, O. Zaharko, P. A. McClarty, S. T. Bramwell, P. Henelius, and T. Fennell, “Pauling entropy, metastability, and equilibrium in $dy_2ti_{20}7$ spin ice,” Phys. Rev. Lett. 121, 067202 (2018).

89 L.D.C. Jaubert, M.J. Harris, T. Fennell, R.G. Melko, S.T. Bramwell, and P.C.W. Holdsworth, “Topological-Sector fluctuations and Curie-Law crossover in spin ice,” Physical Review X 3, 011014 (2013).

90 Andrea Saverio Ninarello, Nicoletta Gnan, and Francesco Sciortino, “Observable-dependence of the effective temperature in off-equilibrium diatomic molecular liquids,” The Journal of Chemical Physics 141, 194507 (2014).

91 John Russo and Francesco Sciortino, “How do self-assembling polymers and gels age compared to glasses?” Phys. Rev. Lett. 104, 195701 (2010).

92 Nicoletta Gnan, Claudio Maggi, Giorgio Parisi, and Francesco Sciortino, “Generalized fluctuation-dissipation relation and effective temperature upon heating a deeply supercooled liquid,” Phys. Rev. Lett. 110, 035701 (2013).

93 Cristiano Nisoli, Roderich Moessner, and Peter Schiffer, “Colloquium: Artificial spin ice: Designing and imaging magnetic frustration,” Rev. Mod. Phys. 85, 1473–1490 (2013).

94 Sandra H. Skjærvø, Christopher H. Marrows, Robert L. Stamps, and Laura J. Heyderman, “Advances in artificial spin ice,” Nature Reviews Physics 2, 13–28 (2020).

95 Demian Levis and Leticia F Cugliandolo, “Out-of-equilibrium dynamics in the bidimensional spin-ice model,” EPL (Europhysics Letters) 97, 30002 (2012).

96 Demian Levis and Leticia F Cugliandolo, “Defects dynamics following thermal quenches in square spin ice,” Physical Review B 87, 214302 (2013).

97 A. Libál, A. del Campo, C. Nisoli, C. Reichhardt, and C. J. O. Reichhardt, “Quenched dynamics of artificial colloidal spin ice,” Phys. Rev. Research 2, 033433 (2020).

98 Yann Perrin, Benjamin Canals, and Nicolas Rouge maille, “Extensive degeneracy, coulomb phase and magnetic monopoles in artificial square ice,” Nature 540, 410–413 (2016).

99 Erik Östman, Henry Stopfel, Ioan-Augustin Chioar, Unnar B. Arnalds, Aaron Stein, Vassilios Kapaklis, and Björgvin Hjörvarsson, “Interaction modifiers in artificial spin ices,” Nature Physics 14, 375–379 (2018).

100 Gavin M. Macauley, Gary W. Paterson, Yue Li, Rair Macêdo, Stephen McVitie, and Robert L. Stamps, “Tuning magnetic order with geometry: Thermalization and defects in two-dimensional artificial spin ices,” Phys. Rev. B 101, 144403 (2020).

101 Oles Sendetskyi, Valerio Scagnoli, Naëmi Leo, Luca Anghinolfi, Aurora Alberca, Jan Lüning, Urs Staub, Peter Michael Derlet, and Laura Jane Heyderman, “Continuous magnetic phase transition in artificial square ice,” Phys. Rev. B 99, 214430 (2019).

102 Vadim Cathelin, *Les glaces de spins : monopoles magnétiques, relation fluctuation-dissipation et fragmentation*, Ph.D. Thesis, Université Grenoble Alpes (2020).

103 C. Paulsen, S. R. Giblin, E. Lhotel, D. Prabhakaran, G. Balakrishnan, K. Matsuhira, and S. T. Bramwell, “Experimental signature of the attractive coulomb force between positive and negative magnetic monopoles in spin ice,” Nature Physics 12, 661–666 (2016).

104 C Paulsen, SR Giblin, Elsa Lhotel, D Prabhakaran, K Matsuhira, G Balakrishnan, and ST Bramwell, “Nuclear spin assisted quantum tunnelling of magnetic monopoles in spin ice,” Nature communications 10, 1–8 (2019).