Temperature chaos is present in off-equilibrium spin-glass dynamics

Experiments featuring non-equilibrium glassy dynamics under temperature changes still await interpretation. There is a widespread feeling that temperature chaos (an extreme sensitivity of the glass to temperature changes) should play a major role but, up to now, this phenomenon has been investigated solely under equilibrium conditions. In fact, the very existence of a chaotic effect in the non-equilibrium dynamics is yet to be established. In this article, we tackle this problem through a large simulation of the 3D Edwards-Anderson model, carried out on the Janus II supercomputer. We find a dynamic effect that closely parallels equilibrium temperature chaos. This dynamic temperature-chaos effect is spatially heterogeneous to a large degree and turns out to be controlled by the spin-glass coherence length ξ. Indeed, an emerging length-scale ξ* rules the crossover from weak (at ξ ≪ ξ*) to strong chaos (ξ ≫ ξ*). Extrapolations of ξ* to relevant experimental conditions are provided. While temperature chaos is an equilibrium notion that denotes the extreme fragility of the glassy phase with respect to temperature changes, it remains unclear whether it is present in non-equilibrium dynamics. Here the authors use the Janus II supercomputer to prove the existence of dynamic temperature chaos, a nonequilibrium phenomenon that closely mimics equilibrium temperature chaos.

Specifically, TC means that the spin configurations that are typical from the Boltzmann weight at temperature T 1 are very atypical at temperature T 2 (no matter how close the two temperatures T 1 and T 2 are).
This equilibrum notion of TC has turned out to be remarkably elusive, even in the context of Mean-Field models (i.e., models that can be solved exactly in the Mean-Field approximation). Indeed, establishing the existence of TC for the Sherrington-Kirkpatrick model has been a real tour de force 12 . Although Sherrington-Kirkpatrick's model is the Mean-Field model of more direct relevance for this work, let us recall for completeness that TC has been investigated as well in other Mean-Field systems named p-spin models. In these models, groups of p ≥ 3 spins interact (instead, p = 2 for Sherrington-Kirkpatrick). Surprisingly enough, one finds different behaviors. On the one hand, we have a recent mathematical proof of the absence of TC in the homogeneous spherical p-spin model 37 , in agreement with a previous claim based on physical arguments 38 . On the other hand, TC should be expected when one mixes several values of p 39 , as confirmed by a quite recent mathematical analysis [40][41][42][43] . Unfortunately, the mathematically rigorous analysis of TC in offequilibrium dynamics seems out of reach for now, even in the Mean-Field context.
In order to obtain experimentally relevant results, one needs to go beyond the Mean-Field approximation and study short-range spin glasses, represented by the Edwards-Anderson model 44,45 . In this case, analytical investigations are even more difficult, but the equilibrium notion of TC that we have outlined above has been studied through numerical simulations. Yet, these equilibrium simulations have been limited to small system sizes by the severe dynamic slowing down [6][7][8]11,13,14,[16][17][18][19] .
Here we tackle the problem from a different approach by showing that a non-equilibrium TC effect is indeed present in the dynamics of a large spin-glass sample in three spatial dimensions (our simulations of the Edwards-Anderson model are carried out on the Janus II custom-built supercomputer 46 ). In a reincarnation of the statics-dynamics equivalence principle [47][48][49][50] , just as equilibrium TC is ruled by the system size, dynamic TC is found to be governed by the time-growing spin-glass coherence length ξ(t w ), where the waiting time t w is the time elapsed since the system was suddenly quenched from some very high temperature to the working temperature T. Below the critical temperature, T < T c , the spin glass is perennially out of equilibrium as evinced by the never-ending (and sluggish) growth of glassy magnetic domains of size ξ(t w ), see refs. 51,52 for instance. Now, the extreme sample-tosample variations found in small equilibrated systems 16,17,19,[53][54][55] translate into a strong spatial heterogeneity of dynamic TC. Despite such strong fluctuations, our large-scale simulations allow us to observe traces of the effect even in averages over the whole system. In close analogy with equilibrium studies 16 , however, dynamic TC can only be fully understood through a statistical analysis of the spatial heterogeneity. A crossover length ξ * emerges such that TC becomes sizeable only when ξ(t w ) > ξ * . We find that ξ * diverges when the two observation temperatures T 1 and T 2 approach. The analysis of this divergence reveals that ξ * is the non-equilibrium partner of the equilibrium chaotic length 3,56 . The large values of ξ(t w ) that we reach with Janus II allow us to perform mild extrapolations to reach the most recent experimental regime 57 .
In equilibrium, sample-averaged signals of TC become more visible when the size of the system increases 16 . Analogously, offequilibrium a weak chaotic effect grows with ξ(t w ) when the whole system is considered on average. Hoping that studying spatial heterogeneities will help us to unveil dynamic TC, we shall consider spatial regions of spherical shape and linear size~ξ(t w ), chosen randomly within a very large spin glass. Statics-dynamics equivalence suggests regarding these spheres as the nonequilibrium analog of small equilibrated samples of linear sizẽ ξ(t w ). The analogy with equilibrium studies 16,17,19 suggests that a small fraction of our spheres will display strong TC. The important question will be how this rare-event phenomenon evolves as ξ(t w ) grows (in equilibrium, the fraction of samples not displaying TC is expected to diminish exponentially with the number of spins contained in the sample 12,15 ).

Results and discussion
Model. We simulate the standard Edwards-Anderson model in a three-dimensional cubic lattice of linear size L = 160 and periodic boundary conditions. In each lattice node x, we place an Ising spin (S x = ± 1). Lattice nearest-neighbors spins interact through the Hamiltonian H = − ∑ 〈x, y〉 J xy S x S y . The couplings J xy are independent and identically distributed random variables (J xy = ± 1 with 1/2 probability), fixed when the simulation starts (quenched disorder). This model exhibits a spin-glass transition at temperature T c = 1.1019 (29) 58 . We refer to each realization of the couplings as a sample. Statistically independent simulations of a given sample are named replicas. We have considerably extended the simulation of Baity et al. 51 , by simulating N Rep = 512 replicas (rather than 256) of the same N S = 16 samples considered in 51 , in the temperature range 0.625 ≤ T ≤ 1.1.
We simulate the non-equilibrium dynamics with a Metropolis algorithm. In this way, one picosecond of physical time roughly corresponds to a full-lattice Metropolis sweep. At the initial time t w = 0 the spin configuration is fully random (i.e., we quench from infinite temperature). The subsequent growth of spin-glass domains is characterized by the spin-glass coherence length ξ(t w ). Specifically, we use the ξ 1,2 integral estimators, see refs. 48,51,59,60 for details [the main steps in the computation of ξ 1,2 are also sketched in Eqs. (8-10)), where one should set T 1 = T 2 = T].
Finally, let us briefly comment on our choices for N Rep and N S . A detailed analysis 51,61 shows that, for a given total numerical effort N S × N Rep , errors in ξ are minimized if N Rep ≫ N S . Furthermore, Supplementary Note 1 shows that having N Rep ≫ 1 is crucial as well for the main quantities considered in this work (see definitions below). Therefore, given our finite computational resources, we have chosen to limit ourselves to N S = 16. This small number of samples is partly compensated by the fact that we are working close to the experimental regime L ≫ ξ [we remark that N S = 1 in typical experiments: indeed, statics-dynamics equivalence suggests that the number of statistically independent events is proportional to N S (L/ξ) 3 ].
The local chaotic parameter. We shall compare the spin textures from temperature T 1 and waiting time t w1 with those from temperature T 2 and waiting time t w2 (we consider T 1 ≤ T 2 ≤ T c ). A fair comparison requires that the two configurations be ordered at the same lengthscale, which we ensure by imposing the condition A first investigation of TC is shown in Fig. 1. The overlap, computed over the whole sample, of two systems satisfying condition Eq. (1) is used to search for a coarse-grained chaotic effect. The resulting signal is measurable but weak. Instead, as explained in the introduction, spin configurations should be compared locally. Specifically, we consider spherical regions. We start by choosing N sph = 8000 centers for the spheres on each sample. The spheres' centers are chosen randomly, with uniform probability, on the dual lattice which, in a cubic lattice with periodic boundary conditions, is another cubic lattice of the same size, also periodic boundary condition. The nodes of the dual lattice are the centers of the elementary cells of the original lattice. The radii of the spheres are varied, but their centers are held fixed. Let B s,r be the s-th ball of radius r. Our basic observable is the overlap between replica σ (at temperature T 1 ), and replica τ ≠ σ (at temperature T 2 ): where N r is the number of spins in the ball, and t w1 and t w2 are chosen according to Eq. (1). Averages over thermal histories, indicated by 〈…〉 T , are computed by averaging over σ and τ. Next, we generalize the so-called chaotic parameter 6,16,17,20 as The extremal values of the chaotic parameter have a simple interpretation: X s;r T 1 ;T 2 ¼ 1 corresponds with a situation in which spin configurations in the ball B s,r , at temperatures T 1 and T 2 , are completely indistinguishable (absence of chaos) while X s;r T 1 ;T 2 ¼ 0 corresponds to completely different configurations (strong TC). A representative example our results is shown in Fig. 2.
Our main focus will be on the distribution function FðX; T 1 ; T 2 ; ξ; rÞ ¼ Probability ½X s;r T 1 ;T 2 ðξÞ < X and on its inverse The rare-event analysis. Representative examples of distribution functions F(X, T 1 , T 2 , ξ, r) are shown in Fig. 3. We see that, in close analogy with equilibrium systems 16,17,19 , while most spheres exhibit a very weak TC (X > 0.9, say), there is a fraction of spheres Fig. 1 Non-equilibrium temperature chaos is weak when averaging over the whole system. We compare typical spin configurations at temperature T 1 and time t w1 with configurations at T 2 and time t w2 . The comparison is carried through a global estimator of the coherence length of their overlap ξ T 1 T 2 1;2 which represents the maximum lengthscale at which configurations at temperatures T 1 and T 2 still look similar, see Methods section for further details. The two times t w1 and t w2 are chosen in such a way that the configurations at both temperatures have glassy-domains of the same size, namely ξ 1,2 (t w1 , T 1 ) = ξ 1,2 (t w2 , T 2 ) = ξ. The figure shows the ratio ξ T 1 T 2 1;2 =ξ as a function of ξ for two pairs of temperatures (T 1 , T 2 ), recall that T c ≈ 1.1 58 . Under the hypothesis of fully developed Temperature Chaos, we would expect ξ T 1 T 2 1;2 to be negligible compared to ξ. Instead, our data shows only a small decrease of ξ T 1 T 2 1;2 =ξ with growing ξ (the larger the difference T 2 − T 1 the more pronounced the decrease). Error bars represent one standard deviation. displaying smaller X (stronger chaos). Note that the probability F of finding spheres with X smaller than any prefixed value increases when ξ grows.
In order to make the above finding quantitative, we consider the (inverse) distribution function X(F, T 1 , T 2 , ξ, r). We start by fixing (T 1 , T 2 ), ξ and some small probability F, which leaves us with a function of only r. In order to obtain smoother interpolations for small radius, however, we have used N 1=3 r instead of r as our independent variable, a technical detailed discussion can be found in Supplementary Note 3. Figure 4 shows plots of 1 − X under these conditions, which exhibit well-defined peaks (see further information about the fitting function to the peaks in the Supplementary Note 2). Now, to a first approximation we can characterize any peak by its position, height and width. Fortunately, these three parameters turn out to describe the scaling with ξ of the full 1 − X curve, see Supplementary Note 4.
The physical interpretation of the peak's parameters is clear. The peak's height represents the strength of dynamic TC (the taller the peak, the larger the chaos). The peak's position indicates the optimal lengthscale for the study of TC, given the probability F, ξ and the temperatures T 1 , T 2 . The peak's width indicates how critical it is to spot this optimal lengthscale (the wider the peak, the less critical the choice). Perhaps unsurprisingly, the peak's position is found to scale linearly with ξ, while the peak's width scales as ξ β , with β slightly larger than one, see Supplementary Note 5 for further details. We shall focus here on the temperature and ξ dependence of the peak's height (i.e., the strength of chaos), which has a richer behavior.
The ξ dependence of the peak's height (for a given probability F and temperatures T 1 and T 2 ) turns out to be reasonably well described by the following ansatz: This formula describes a crossover phenomenon, ruled by a characteristic length ξ * . For ξ ≪ ξ * the peak's height grows with ξ as a power law, while for ξ ≫ ξ * the strong-chaos limit [i.e., (1 − X) → 1] is approached. However, some consistency requirements should be met before taking the crossover length ξ * seriously. Not only should the fit to Eq. (4) be of acceptable statistical quality (the fit parameters are the characteristic lengthscale ξ * and the exponent α). One would also wish exponent α to be independent of the temperatures T 1 and T 2 and of the chosen probability F. We find fair fits to Eq. (4), see Table 1. In all cases, exponent α turns out to be compatible with 2.1 at the two-σ level [except for the (F = 0.01, T 1 = 0.625, T 2 = 0.8) fit]. Under these conditions, we can interpret ξ * as a characteristic length indicating the crossover from weak to strong TC, at the probability level indicated by F. Furthermore, the relatively large value of exponent α indicates that this crossover is sharp.
The trends for the crossover length ξ * in Table 1 are very clear: ξ * grows upon increasing F or upon decreasing T 2 − T 1 . Identifying ξ * as the non-equilibrium partner of the equilibrium chaotic length ℓ c (T 1 , T 2 ) 3,56 will allow us to be more quantitative (indeed, the two lengthscales indicate the crossover between weak Error bars, that represent one standard deviation, are horizontal, because we have actually extrapolated the chaotic parameter, which is its inverse function X(F, T 1 , T 2 , ξ, r). Most of the spheres have a chaotic parameter very close to X = 1 (absence of chaos). However, if we fix our attention, for instance, on percentile 1 (i.e., F = 0.01) we see that the corresponding value of X decreases monotonically (and significantly) as ξ grows, signaling a developing chaotic effect. This trend is clear both for spheres of radius r = 4 and r = 8. Fig. 4 Emergence of an optimal scale to observe temperature chaos. The difference 1 − X(F, T 1 , T 2 , ξ, r) [recall that X(F, T 1 , T 2 , ξ, r) is the inverse of the distribution function] as a function of the cubic root N 1=3 r of the number of spins in the spheres, as computed for different values of the probability level F, the temperatures T 1 and T 2 , and the coherence length ξ. In this representation, the optimal size of the spheres for the observation of chaos (for given parameters F, T 1 , T 2 and ξ) appears as the maximum of the curves. Continuous lines are fits to a smooth interpolating function, further details can be found in Supplementary Note 2. Error bars represent one standard deviation. chaos and strong chaos). Now, the equilibrium ℓ c (T 1 , T 2 ) has been found to scale for the 3D Ising spin glass as with ζ ≈ 1.07 14 or ζ ≈ 1.07(5) 16 . These considerations suggest the following ansatz for the non-equilibrium crossover length where B(F, T 1 ) is an amplitude. We have tested Eq. (6) by computing a joint fit for four (T 1 , F) pairs as functions of T 2 − T 1 , allowing each curve to have its own amplitude but enforcing a common ζ NE (see Fig. 5). The resulting χ 2 /d.o.f. = 7.55/7 validates our ansatz, with an exponent ζ NE = 1.19(2) fairly close to the equilibrium result ζ = 1.07(5) 16 . This agreement strongly supports our physical interpretation of the crossover length. We, furthermore, find that B is only weakly dependent on T 1 . Nevertheless, the reader should be warned that it has been suggested 16 that the equilibrium exponent ζ may be different in the weak-and strong-chaos regimes.

Conclusions
We have shown that the concept of temperature chaos can be meaningfully extended to the non-equilibrium dynamics of a large spin glass. This is, precisely, the framework for rejuvenation and memory experiments [28][29][30][31] , as well as other more chaosoriented experimental work 57 . Therefore, our precise characterization of dynamical temperature chaos paves the way for the interpretation of these and forthcoming experiments. Our simulation of spin-glass dynamics doubles the numerical effort in 51 and has been carried out on the Janus-II special-purpose supercomputer.
The key quantity governing dynamic temperature chaos is the time-dependent spin-glass coherence length ξ(t w ). The very strong spatial heterogeneity of this phenomenon is quantified through a distribution function F. This probability can be thought of as the fraction of the sample that shows a chaotic response to a given degree. When comparing temperatures T 1 and T 2 , the degree of chaoticity is governed by a lengthscale ξ * (F, T 1 , T 2 ). While chaos is very weak if ξ(t w ) ≪ ξ * (F, T 1 , T 2 ), it quickly becomes strong as ξ(t w ) approaches ξ * (F, T 1 , T 2 ). We find that, when T 1 approaches T 2 , ξ * (F, T 1 , T 2 ) appears to diverge with the same critical exponent that it is found for the equilibrium chaotic length 16 .
Although we have considered in this work fairly small values of the chaotic system fraction F, a simple extrapolation, linear in log F, predicts ξ * ≈ 60 for F = 0.1 at T 1 = 0.7 and T 2 = 0.8 (our closest pair of temperatures in Table 1). A spin-glass coherence length well above 60a 0 is experimentally reachable nowadays 52,57,62,63 (a 0 is the typical spacing between spins), which makes our dynamic temperature chaos significant. Indeed, while completing this manuscript, a closely related experimental study 57 reported a value for exponent ζ NE in fairly good agreement with our result of ζ NE = 1.19 (2) in Fig. 5.
Let us conclude by commenting on possible venues for future research. Clearly, it will be important to understand in detail how dynamic temperature chaos manifests itself in non-equilibrium experiments. Simple protocols (in which temperature sharply drops from T 2 to T 1 , see, e.g., Zhai et. al. 57 ) seems more accessible to a first analysis than memory and rejuvenation experiments [28][29][30][31] . An important problem is that the correlation functions that are studied theoretically are not easily probed experimentally. Instead, experimentalists privilege the magnetization density (which is a spatial average over the whole sample). Therefore an important theoretical goal is to predict the behavior of the non-equilibrium timedependent magnetization upon a temperature drop. One may speculate that the Generalized Fluctuation-Dissipation Relations 64 might be the route connecting the correlation functions with the response to an externally applied magnetic field. Interestingly enough, these relations (that apply at fixed temperature) can be defined locally as well 65 . The resulting spatial distribution function allows the reconstruction of the global response to the magnetic field. Extending this analysis to a temperature drop may turn out to be fruitful in the future.

Methods
All the observables involved in the computation of temperature chaos depend on a pair of replicas (σ, τ). The basic quantity is the overlap field Usually, this pair of replicas are at the same temperature T. All the definitions are, however, straightforwardly extended to two temperatures. For instance, the fourpoint two-temperature spatial correlation function is  Parameters obtained in the fits to Eq. (4) of our data for the peak's height, see Fig. 4, with ξ min ≤ ξ ≤ ξ max . We also report the fits' figure of merit χ 2 /d.o.f. We chose ξ min by requiring a P value greater than 0.05 in the fits (ξ max = 9.5 for T 1 = 0.625 and ξ max = 12.5 for T 1 = 0.7). T 1 and T 2 represent the temperatures involved in the computation of the chaotic parameter. Unfortunately, the flatness of the peak for (T 1 = 0.625, T 2 = 0.7, F = 0.01) did not allow us to compute the peak's parameters. ðr; t w1 ; t w2 Þ dr; ð9Þ and ξ T 1 T 2 k;kþ1 ðt w1 ; t w2 Þ ¼ As explained in the main text, times t w1 and t w2 are fixed through the condition expressed in Eq. (1), which ensures that we are comparing spin configurations that are ordered on the same length scale. Since our t w are on a discrete grid, we solve Eq. (1) for the global overlaps through a (bi)linear interpolation.

Data availability
The data contained in the figures of this paper, accompanied by the gnuplot script files that generate these figures, are publicly available at https://github.com/ JanusCollaboration/caosdin.
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The codes that support the findings of this study are available from the corresponding author upon reasonable request.