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

We find a dynamic effect in the non-equilibrium dynamics of a spin glass that closely parallels equilibrium temperature chaos. This effect, that we name dynamic temperature chaos, is spatially heterogeneous to a large degree. The key controlling quantity is the time-growing spin-glass coherence length. Our detailed characterization of dynamic temperature chaos paves the way for the analysis of recent and forthcoming experiments. This work has been made possible thanks to the most massive simulation to date of non-equilibrium dynamics, carried out on the Janus~II custom-built supercomputer.

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 off-equilibrium 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-8, 11, 13, 14, 16-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 custombuilt 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-to-sample 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, off-equilibrium a weak chaotic effect grows with ξ(t w ) when the whole system is considered on av-erage. 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. Staticsdynamics equivalence suggests regarding these spheres as the non-equilibrium analogue of small equilibrated samples of linear size ∼ ξ(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,9,10), where one should set 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 which represents the maximum lengthscale at which configurations at temperatures T1 and T2 still look similar, see Methods section for further details. The two times tw1 and tw2 are chosen in such a way that the configurations at both temperatures have glassy-domains of the same size, namely ξ1,2(tw1, T1) = ξ1,2(tw2, T2) = ξ. The figure shows the ratio ξ T 1 T 2 1,2 /ξ as a function of ξ for two pairs of temperatures (T1, T2), recall that Tc ≈ 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 T2 − T1 the more pronounced the decrease). Error bars represent one standard deviation. 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 ): q s,r;σ,τ T1,T2 (ξ) = 1 N r x∈Bs,r s σ,T1 x (t w1 )s τ,T2 x (t w2 ) , (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 T1,T2 = 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 T1,T2 = 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 T1,T2 (ξ) < X] and on its inverse X(F, T 1 , T 2 , ξ, r). 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 Temperature chaos increases with the coherence length. The figure shows the distribution function F (X, T1, T2, ξ, r) for temperatures T1 = 0.625 and T2 = 0.9, for spheres of radius r = 4 and r = 8, as computed for various values of ξ. The distributions have been extrapolated to infinite number of replicas NRep = ∞, see Supplementary Note 1 for further details. Error bars, that represent one standard deviation, are horizontal, because we have actually extrapolated the chaotic parameter, which is its inverse function X(F, T1, T2, ξ, 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.
systems [16,17,19], while most spheres exhibit a very weak TC (X > 0.9, say), there is a fraction of spheres 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. Fig. 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 Emergence of an optimal scale to observe temperature chaos. The difference 1 − X(F, T1, T2, ξ, r) [recall that X(F, T1, T2, ξ, 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 T1 and T2, and the coherence length ξ. In this representation, the optimal size of the spheres for the observation of chaos (for given parameters F , T1, T2 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.
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 strongchaos 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 I. 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 I are very clear: ξ * grows upon increasing F or upon decreasing T 2 − T 1 . Identifying ξ * as the non-equilibrium part- 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 T1 = 0.625 and ξmax = 12.5 for T1 = 0.7). T1 and T2 represent the temperatures involved in the computation of the chaotic parameter. Unfortunately, the flatness of the peak for (T1 = 0.625, T2 = 0.7, F = 0.01) did not allow us to compute the peak's parameters.
ner 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 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.  (2). Error bars represent one standard deviation.

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 chaos-oriented 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 I). 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 time-dependent 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 where [. . .] J denotes the average over the samples. Building on this function we can define our integral estimator for the coherence length [60]: and ξ T1T2 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.

ACKNOWLEDGMENTS
We are grateful for discussions with R. Orbach and Q. Zhai. This work was partially supported by Ministerio de Economía, Industria y Competitividad (MINECO, Spain), Agencia Estatal de Investigación (AEI, Spain), and Fondo Europeo de Desarrollo Regional The thermal expectation values necessary to compute the chaotic parameter are defined in the limit of infinite replicas, so an extrapolation is necessary to avoid bias. Fortunately, with all other parameters fixed, the evolution of X(F, T 1 , T 2 , ξ, r) with N Rep is smooth (see Fig. 6 and Fig. 7). We have mainly used a linear ansatz for the extrapolation, where X NRep is shorthand for X(F, T 1 , T 2 , ξ, r; N Rep ), X ∞ = X(F, T 1 , T 2 , ξ, r; N Rep → ∞) and A is a constant. As a check for the linear ansatz in Eq. (11), we have considered two alternative functional forms: where B, C and D are amplitudes and γ is a free exponent. We perform independent fits to Eq. (11), Eq. (12) and Eq. (13) for every value of the parameters (F, T 1 , T 2 , ξ, r). We reject fits with a diagonal Errors in X ∞ are computed by performing separate fits for each jackknife block (the fitting procedure consists in minimizing the diagonal χ 2 , see [60]). As a first check, we compare the linear and quadratic extrapolations (see Fig.6 for an illustrative example). The figure shows that even for our largest N Rep , namely N Rep = 256 and N Rep = 512, we are still far from the N Rep → ∞ limit. Fortunately, the linear and the quadratic extrapolations yield compatible results in our region of interest, i.e., the tail of the distribution function. We remark that the consistency condition χ 2 /d.o.f < 1.1 is met in a larger range for the quadratic extrapolation (F < 0.9) than for the linear extrapolation (F < 0.7). However, because both coincide in the low-F range we are interested in, we have kept the simpler linear extrapolation.  (12), turn out to be equivalent for the tail of the distribution function. Continuous lines are the linear (golden curves) and quadratic (blue curves) extrapolations to NRep → ∞ for F (X, T1, T2, ξ, r) as a function of X. The data shown correspond to the case T1 = 0.7, T2 = 0.8, ξ = 11 and r = 8. The two curves shown for each extrapolation correspond to the central value plus or minus the standard error. We show horizontal error bars because we are computing the inverse distribution function X(F, T1, T2, ξ, r). We only show extrapolated data when χ 2 /d.o.f < 1.1 in the fits to Eq. 11 or to Eq. 12. For comparison, we also plot the data corresponding to NRep = 512 and NRep = 256. Inset: As in the main plot, but with the vertical axis in log scale.
Our second check in Eq. (13) seeks the natural exponent γ for the extrapolation as a fitting parameter. We have found that the consistency condition χ 2 /d.o.f < 1.1 is met for F < 0.85. Fortunately, γ turns out to be very close to the value γ = 1, with the exception of an instability in the crossing region around F ≈ 0.35, see Fig. 7.
In summary, the quadratic and the free-exponent extrapolations support our choice of Eq. 11 as the preferred form for the N Rep → ∞ extrapolation.

SUPPLEMENTARY NOTE 2: CHARACTERIZATION OF THE PEAK
The complementary of the chaotic parameter 1 − X, as a function of the sphere size, has a well-defined peak. Characterizing the peak is fundamental to the analysis because it provides information about the optimal sphere size for the study of temperature chaos and about the strength of the phenomenon.
Let us remark that, at least close to a maximum, any smooth curve is characterized by the position, height and width of the peak. In order to meaningfully compute these three parameters from our data, we fit 1 − X to Eq. (14) with z = N 1/3 r . We extract the position, width and height of the maximum from the fitted function f (z). Errors are computed with a jackknife method [60]. In this section we explain our rationale for choosing the cubic root of the number of spins contained in the sphere, N 1/3 r , rather than its radius r, to characterize the size of the spheres considered in our analysis.
The objective is fitting the peaks of 1−X to a function of the form with a, b, c, and d as fit parameters. If one uses the obvious choice of z = r, however, the fit fails. Indeed, see Fig. 8 upper panel, 1 − X is not a smooth function of r. The reason is that the number of lattice points in the spheres is not a smooth function of r either (see Fig. 8 lower panel). It is natural, therefore, to replace r with N 1/3 r as independent variable. This substitution makes Eq. (14) work down to smaller spheres. The difference between both independent variables becomes negligible for very large spheres. We enlarge the small-sphere region, where both independent variables most differ. Error bars represent one standard deviation. Lower panel: the cubic root of the volume of a sphere (blue curve) is plotted as a function of the radius of the sphere r. The golden curve is N 1/3 r , namely the cubic root of the number of lattice points contained in a sphere of radius r, centered at a node of the dual to our cubic lattice. Values of N 1/3 r corresponding to integer r are highlighted as black dots.

SUPPLEMENTARY NOTE 4: GLOBAL VERSUS LOCAL DESCRIPTION OF THE PEAKS
In the main text, we reduce the study of the scaling of temperature chaos with the coherence length ξ to the study of the peak of 1 − X against the size of the spheres. The reader may wonder whether the local fit of the peak will extend to describe the full curve. Here we present some positive evidences in this respect. to Eq. (16). For each fit, we also report the figure of merit χ 2 /d.o.f (we include in the fit data with ξ ≥ ξmin ; ξmin is set by requiring the fit's P value to be larger than 0.05).
Consider any smooth, positive function H(z), with a local maximum at z = z max . Close to this peak, Taylor's theorem implies some (trivial) universality where y = |H (zmax)| H(zmax) (z − z max ). Note that the peak's position is z max , its heigth is H(z max ) and its (inverse) width is |H (z max )|/H(z max ). In principle, there is no reason for Eq. (15) to be accurate away from the peak, but this formula suggests an alternative representation for our 1 − X curves, see Fig. 9. We note that, in this new representation, the 1 − X curves are invariant under changes of the coherence length ξ (upper panel). When considering changes in the temperatures T 1 and T 2 and the probability F , however, the curves mildly differ away from the peak (see Fig. 9 lower panel). This (approximate) independence of (T 1 , T 2 , F, ξ) is a fortunate fact because the complexity of the problem gets reduced to the study of the scaling with ξ of the three peak parameters while keeping (T 1 , T 2 , F ) constant.

SUPPLEMENTARY NOTE 5: POSITION AND WIDTH OF THE PEAKS
In this section we analyze the scaling of the peaks' position and width with the coherence length ξ.
We first focus on the peak's position, which signals the most convenient length scale for studying TC (for a given coherence length ξ, probability F and temperatures T 1 and T 2 ). Dimensional analysis suggests a linear fit as the natural ansatz to study the scaling of the peak's FIG. 9. Universality in 1 − X extends beyond the trivial Taylor's universality. The upper panel shows 1 − X in units of its peak value, for the temperatures T1 = 0.7, T2 = 1.0 and F = 0.01. Taylor's theorem implies that, using the independent variable y [see Eq. (15)], the different curves should coincide close to y = 0. However, we see that the collapse holds beyond the quadratic approximation (as evinced by the strong asymmetry of the master curve). The lower panel mixes different values of F , T1 and T2, which leads to significant discrepancies for large values of |y|. Nevertheless, the curves still collapse in a range y ∈ (−0.3, 0.5) where the asymmetry is significative. Error bars represent one standard deviation. position N 1/3 r,max with the coherence length ξ(t w ) (indeed, both quantities are lengths): Fits of the data to Eq. 16 are shown in Table II. In all cases, values of parameter b are compatible with 0 (at the two-σ level). In addition, amplitude a exhibits monotonic behavior with T 2 −T 1 and with the probability F . Hence, our naive expectation N 1/3 r,max ∝ ξ(t w ) is confirmed. The peak's width determines how delicate the selection of the right length scale is to study TC (i.e., if the peak's width becomes larger than its position, this choice is no  (17). For each fit, we also report the figure of merit χ 2 /d.o.f (we include in the fit data with ξ ≥ ξmin ; ξmin is set by requiring the fit's P value to be larger than 0.05).
We study the inverse peak's width (i.e., the curvature κ(ξ)) and propose a power law decaying with ξ(t w ) characterized by the ansatz Results are shown in Table III.
The value of A(F ) turns out to be compatible for all pairs of temperatures (T 1 , T 2 ) at fixed probability F . Furthermore, at the current precision of the data, exponent β does not exhibit any significant dependence either on the temperature pair (T 1 , T 2 ) or on the probability F .
Let us now recall the linear relation between the peak's position and the coherence length, see Eq. (16). Consider the ratio between the position of the maximum and its width, N r,max κ(ξ) ∼ ξ 1−β . Table III mildly suggest that β is slightly greater than 1, which implies that the ratio goes to zero (very slowly) in the limit of large ξ.