# Memory and rejuvenation effects in spin glasses are governed by more than one length scale

Marco Baity-Jesi, Enrico Calore, Andrea Cruz, Luis Fernandez, José Miguel Gil-Narvion, I. Gonzalez-Adalid Pemartin, Antonio Gordillo-Guerrero, David Iñiguez, Andrea Maiorano, Enzo Marinari, et al.

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

## To cite this version:

Marco Baity-Jesi, Enrico Calore, Andrea Cruz, Luis Fernandez, José Miguel Gil-Narvion, et al. Memory and rejuvenation effects in spin glasses are governed by more than one length scale. Nature Physics, 2023, 19 (7), pp.978-985. 10.1038/s41567-023-02014-6. hal-04341897

HAL Id: hal-04341897

https://inria.hal.science/hal-04341897

Submitted on 4 Apr 2024

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_2.png" /></figure>
# Memory and rejuvenation effects in spin glasses are governed by more than one length scale

Marco Baity-Jesi, Enrico Calore, Andrea Cruz, Luis Fernandez, José Miguel Gil-Narvion, I. Gonzalez-Adalid Pemartin, Antonio Gordillo-Guerrero, David Iñiguez, Andrea Maiorano, Enzo Marinari, et al.

► To cite this version:

Marco Baity-Jesi, Enrico Calore, Andrea Cruz, Luis Fernandez, José Miguel Gil-Narvion, et al. Memory and rejuvenation effects in spin glasses are governed by more than one length scale. Nature Physics, 2023, 19 (7), pp.978-985. 10.1038/s41567-023-02014-6. hal-04341897

HAL Id: hal-04341897

https://inria.hal.science/hal-04341897v1

Submitted on 4 Apr 2024

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
# Memory and rejuvenation in spin glasses: aging systems are ruled by more than one length scale

M. Baity-Jesi,¹ E. Calore,² A. Cruz,³,⁴ L.A. Fernandez,⁵,⁴ J.M. Gil-Narvion,⁴ I. Gonzalez-Adalid Pemartin,⁵

A. Gordillo-Guerrero,⁶,⁷,⁴ D. Íñiguez,⁴,⁸ A. Maiorano,⁹,⁴ E. Marinari,¹⁰ V. Martin-Mayor,⁵,⁴

J. Moreno-Gordo,⁴,³,¹¹,⁷ A. Muñoz Sudupe,⁵,⁴ D. Navarro,¹² I. Paga,¹³,* G. Parisi,¹⁰ S. Perez-Gaviro,³,⁴

F. Ricci-Tersenghi,¹⁰ J.J. Ruiz-Lorenzo,¹¹,⁷,⁴ S.F. Schifano,¹⁴ B. Seoane,⁵,⁴ A. Tarancon,³,⁴ and D. Yllanes¹⁵,⁴

(Janus Collaboration)

¹Eawag, Überlandstrasse 133, CH-8600 Dübendorf, Switzerland

²Dipartimento di Fisica e Scienze della Terra,
Università di Ferrara and INFN, Sezione di Ferrara, I-44122 Ferrara, Italy

³Departamento de Física Teórica, Universidad de Zaragoza, 50009 Zaragoza, Spain

⁴Instituto de Biocomputación y Física de Sistemas Complejos (BIFI), 50018 Zaragoza, Spain

⁵Departamento de Física Teórica, Universidad Complutense, 28040 Madrid, Spain

⁶Departamento de Ingeniería Eléctrica, Electrónica y Automática, U. de Extremadura, 10003, Cáceres, Spain

⁷Instituto de Computación Científica Avanzada (ICCAEx),
Universidad de Extremadura, 06006 Badajoz, Spain

⁸Fundación ARAID, Diputación General de Aragón, 50018 Zaragoza, Spain

⁹Dipartimento di Biotecnologie, Chimica e Farmacia, Università degli studi di Siena,
53100 Siena, Italy and INFN, Sezione di Roma 1, 00185 Rome, Italy

¹⁰Dipartimento di Fisica, Sapienza Università di Roma, and CNR-Nanotec,
Rome unit and INFN, Sezione di Roma 1, 00185 Rome, Italy

¹¹Departamento de Física, Universidad de Extremadura, 06006 Badajoz, Spain

¹²Departamento de Ingeniería, Electrónica y Comunicaciones and I3A, U. de Zaragoza, 50018 Zaragoza, Spain

¹³Dipartimento di Fisica, Sapienza Università di Roma,
and CNR-Nanotec, Rome unit, Sezione di Roma 1, 00185 Rome, Italy

¹⁴Dipartimento di Scienze Chimiche e Farmaceutiche,
Università di Ferrara e INFN Sezione di Ferrara, I-44122 Ferrara, Italy

¹⁵Chan Zuckerberg Biohub, San Francisco, CA, 94158

(Dated: August 31, 2022)

Memory and rejuvenation effects in the magnetic response of off-equilibrium spin glasses have been widely regarded as the doorway into the experimental exploration of ultrametricity and temperature chaos (maybe the most exotic features in glassy free-energy landscapes). Unfortunately, despite more than twenty years of theoretical efforts following the experimental discovery of memory and rejuvenation, these effects have thus far been impossible to simulate reliably. Yet, three recent developments convinced us to accept this challenge: first, the custom-built Janus II supercomputer makes it possible to carry out “numerical experiments” in which the very same quantities that can be measured in single crystals of CuMn are computed from the simulation, allowing for parallel analysis of the simulation/experiment data. Second, Janus II simulations have taught us how numerical and experimental length scales should be compared. Third, we have recently understood how temperature chaos materializes in aging dynamics. All three aspects have proved crucial for reliably reproducing rejuvenation and memory effects on the computer. Our analysis shows that (at least) three different length scales play a key role in aging dynamics, while essentially all theoretical analyses of the aging dynamics emphasize the presence and the crucial role of a single glassy correlation length.

The remarkable off-equilibrium behavior of glass formers at low temperatures has been described with terms such as *aging* [1] or *memory* and *rejuvenation* [2–5], which seem more suitable for living beings than for inert chunks of matter. In this context, spin glasses (which are disordered magnetic alloys, see, e.g., [6]) enjoy a privileged status. On the experimental side, their magnetic response can be studied with great accuracy using a superconducting quantum interference device (SQUID). Rejuvenation and memory (see the description below) are, furthermore, remarkably strong in spin glasses, probably because of

the large correlation length $\xi$ of the coherent spin domains. The values of $\xi$ reached in single-crystal samples [7–10] is much larger than in other glass-forming materials (for instance, the $\xi$ measured in supercooled glycerol or propylene carbonate [11] is smaller by a factor $\sim 100$). On the other hand, spin-glass theory [12] has proved applicable to distant fields that also feature rugged free-energy landscapes, such as combinatorial optimization, machine learning, biology, financial markets or social dynamics.

It is worth stressing that the main part of spin-glass experimental studies is carried out under off-equilibrium conditions [14]. In the simplest setting, the so-called zero-field-cooling (ZFC) protocol, the system is initially at equilibrium at some very high temperature. Eventually the spin glass is abruptly

<footnote>* ilaria.paga@gmail.com</footnote>
Memory and rejuvenation in spin glasses: aging systems are ruled by more than one length scale

M. Baity-Jesi,¹ E. Calore,² A. Cruz,³,⁴ L.A. Fernandez,⁵,⁴ J.M. Gil-Narvion,⁴ I. Gonzalez-Adalid Pemartin,⁵

A. Gordillo-Guerrero,⁶,⁷,⁴ D. Íñiguez,⁴,⁸ A. Maiorano,⁹,⁴ E. Marinari,¹⁰ V. Martin-Mayor,⁵,⁴

J. Moreno-Gordo,⁴,³,¹¹,⁷ A. Muñoz Sudupe,⁵,⁴ D. Navarro,¹² I. Paga,¹³,* G. Parisi,¹⁰ S. Perez-Gaviro,³,⁴

F. Ricci-Tersenghi,¹⁰ J.J. Ruiz-Lorenzo,¹¹,⁷,⁴ S.F. Schifano,¹⁴ B. Seoane,⁵,⁴ A. Tarancon,³,⁴ and D. Yllanes¹⁵,⁴

(Janus Collaboration)

¹Eawag, Überlandstrasse 133, CH-8600 Dübendorf, Switzerland

²Dipartimento di Fisica e Scienze della Terra,
Università di Ferrara and INFN, Sezione di Ferrara, I-44122 Ferrara, Italy

³Departamento de Física Teórica, Universidad de Zaragoza, 50009 Zaragoza, Spain

⁴Instituto de Biocomputación y Física de Sistemas Complejos (BIFI), 50018 Zaragoza, Spain

⁵Departamento de Física Teórica, Universidad Complutense, 28040 Madrid, Spain

⁶Departamento de Ingeniería Eléctrica, Electrónica y Automática, U. de Extremadura, 10003, Cáceres, Spain

⁷Instituto de Computación Científica Avanzada (ICCAEx),
Universidad de Extremadura, 06006 Badajoz, Spain

⁸Fundación ARAID, Diputación General de Aragón, 50018 Zaragoza, Spain

⁹Dipartimento di Biotecnologie, Chimica e Farmacia, Università degli studi di Siena,
53100 Siena, Italy and INFN, Sezione di Roma 1, 00185 Rome, Italy

¹⁰Dipartimento di Fisica, Sapienza Università di Roma, and CNR-Nanotec,
Rome unit and INFN, Sezione di Roma 1, 00185 Rome, Italy

¹¹Departamento de Física, Universidad de Extremadura, 06006 Badajoz, Spain

¹²Departamento de Ingeniería, Electrónica y Comunicaciones and I3A, U. de Zaragoza, 50018 Zaragoza, Spain

¹³Dipartimento di Fisica, Sapienza Università di Roma,
and CNR-Nanotec, Rome unit, Sezione di Roma 1, 00185 Rome, Italy

¹⁴Dipartimento di Scienze Chimiche e Farmaceutiche,
Università di Ferrara e INFN Sezione di Ferrara, I-44122 Ferrara, Italy

¹⁵Chan Zuckerberg Biohub, San Francisco, CA, 94158

(Dated: August 31, 2022)

Memory and rejuvenation effects in the magnetic response of off-equilibrium spin glasses have been widely regarded as the doorway into the experimental exploration of ultrametricity and temperature chaos (maybe the most exotic features in glassy free-energy landscapes). Unfortunately, despite more than twenty years of theoretical efforts following the experimental discovery of memory and rejuvenation, these effects have thus far been impossible to simulate reliably. Yet, three recent developments convinced us to accept this challenge: first, the custom-built Janus II supercomputer makes it possible to carry out “numerical experiments” in which the very same quantities that can be measured in single crystals of CuMn are computed from the simulation, allowing for parallel analysis of the simulation/experiment data. Second, Janus II simulations have taught us how numerical and experimental length scales should be compared. Third, we have recently understood how temperature chaos materializes in aging dynamics. All three aspects have proved crucial for reliably reproducing rejuvenation and memory effects on the computer. Our analysis shows that (at least) three different length scales play a key role in aging dynamics, while essentially all theoretical analyses of the aging dynamics emphasize the presence and the crucial role of a single glassy correlation length.

The remarkable off-equilibrium behavior of glass formers at low temperatures has been described with terms such as *aging* [1] or *memory* and *rejuvenation* [2–5], which seem more suitable for living beings than for inert chunks of matter. In this context, spin glasses (which are disordered magnetic alloys, see, e.g., [6]) enjoy a privileged status. On the experimental side, their magnetic response can be studied with great accuracy using a superconducting quantum interference device (SQUID). Rejuvenation and memory (see the description below) are, furthermore, remarkably strong in spin glasses, probably because of

the large correlation length $\xi$ of the coherent spin domains. The values of $\xi$ reached in single-crystal samples [7–10] is much larger than in other glass-forming materials (for instance, the $\xi$ measured in supercooled glycerol or propylene carbonate [11] is smaller by a factor $\sim 100$). On the other hand, spin-glass theory [12] has proved applicable to distant fields that also feature rugged free-energy landscapes, such as combinatorial optimization, machine learning, biology, financial markets or social dynamics.

It is worth stressing that the main part of spin-glass experimental studies is carried out under off-equilibrium conditions [14]. In the simplest setting, the so-called zero-field-cooling (ZFC) protocol, the system is initially at equilibrium at some very high temperature. Eventually the spin glass is abruptly

<footnote>* ilaria.paga@gmail.com</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. The zero-field-cooling (ZFC) numerical experiment measuring rejuvenation and memory. The starting random spin configuration is placed instantaneously at the working temperature and it relaxes for a time $t_w$ without a field. At time $t_w$, a magnetic field $H = 0.01$ is applied and the magnetic density, $M_{ZFC}(t, t_w; H)$, is recorded. Left panels show the relaxation function $S_{ZFC}(t, t_w; H)$, Eq. (1), for the native runs at the warmer, $T_1 = 0.9$, and colder, $T_2 = 0.5$, temperatures (both below the glass temperature $T_g = 1.102(3)$ [13]). The physically interesting peak of $S_{ZFC}(t, t_w; H)$ defines $t_H^{eff} \simeq t_w$ (in some cases there is also a peak at shorter times, see Methods). In our protocol (schematized by the green arrows), after a waiting time $t_w^{\downarrow} = 2^{31.25}$, the temperature is abruptly dropped from the initial temperature $T_1 = 0.9$ to the colder temperature $T_2 = 0.5$. Then, the system relaxes at $T_2$ for an additional time, after which the magnetic field is switched on and the function $S_{ZFC}(t, t_w; H)$ shown in the bottom-center panel is measured. Waiting times for these jump runs are reported in the legend; the rejuvenation effect is clearly visible, since $t_H^{eff} \ll t_w^{\downarrow}$ and similar to the time $t_w$ spent at $T_2$ (we use $t_w$ for the time spent at the last temperature in a given protocol). Finally, after the waiting time $t_w^{\uparrow} = 2t_w^{\downarrow} = 2^{32.25}$ [i.e., the system has spent half of its life at the initial temperature $T_1$ and half at the colder temperature $T_2$ without a field], the spin glass is suddenly heated back to $T_1$. We let the system relax for a short time, $t_w = 2^{10} \ll t_w^{\downarrow}$, after which the magnetic field is switched on. The $S_{ZFC}(t, t_w; H)$ measured after the jump back and shown in the top-right panel has a peak very similar to the one before the first jump (see top-left panel), evincing the memory of the aging at the initial temperature $T_1$, notwithstanding the rejuvenation observed when staying at the lower temperature $T_2$. In all cases, error bars are one standard deviation.</figcaption></figure>

cooled to the working temperature $T < T_g$ and relaxes for a waiting time $t_w$ ($T_g$ is the glass temperature, while $t_w$ ranges from minutes to several hours). At time $t_w$ a magnetic field $H$ is switched on and the growing magnetization $M_{ZFC}(t, t_w)$ is recorded at later times $t+t_w$. $M_{ZFC}(t, t_w)$ has turned out to have a

significant dependence on $t_w$ for as long as researchers have had the patience to wait. The relaxation rate

$$S_{ZFC}(t, t_w; H) = \frac{\frac{1}{H} \frac{dM_{ZFC}(t, t_w; H)}{d \log t} \quad (1)$$

peaks at a time $t_H^{eff}$ roughly equal to $t_w$ (see, e.g.,
<figure><img src="image_2.png" /><figcaption>FIG. 1. The zero-field-cooling (ZFC) numerical experiment measuring rejuvenation and memory. The starting random spin configuration is placed instantaneously at the working temperature and it relaxes for a time $t_w$ without a field. At time $t_w$, a magnetic field $H = 0.01$ is applied and the magnetic density, $M_{ZFC}(t, t_w; H)$, is recorded. Left panels show the relaxation function $S_{ZFC}(t, t_w; H)$, Eq. (1), for the native runs at the warmer, $T_1 = 0.9$, and colder, $T_2 = 0.5$, temperatures (both below the glass temperature $T_g = 1.102(3)$ [13]). The physically interesting peak of $S_{ZFC}(t, t_w; H)$ defines $t_H^{eff} \simeq t_w$ (in some cases there is also a peak at shorter times, see Methods). In our protocol (schematized by the green arrows), after a waiting time $t_w^{\downarrow} = 2^{31.25}$, the temperature is abruptly dropped from the initial temperature $T_1 = 0.9$ to the colder temperature $T_2 = 0.5$. Then, the system relaxes at $T_2$ for an additional time, after which the magnetic field is switched on and the function $S_{ZFC}(t, t_w; H)$ shown in the bottom-center panel is measured. Waiting times for these jump runs are reported in the legend; the rejuvenation effect is clearly visible, since $t_H^{eff} \ll t_w^{\downarrow}$ and similar to the time $t_w$ spent at $T_2$ (we use $t_w$ for the time spent at the last temperature in a given protocol). Finally, after the waiting time $t_w^{\uparrow} = 2t_w^{\downarrow} = 2^{32.25}$ [i.e., the system has spent half of its life at the initial temperature $T_1$ and half at the colder temperature $T_2$ without a field], the spin glass is suddenly heated back to $T_1$. We let the system relax for a short time, $t_w = 2^{10} \ll t_w^{\downarrow}$, after which the magnetic field is switched on. The $S_{ZFC}(t, t_w; H)$ measured after the jump back and shown in the top-right panel has a peak very similar to the one before the first jump (see top-left panel), evincing the memory of the aging at the initial temperature $T_1$, notwithstanding the rejuvenation observed when staying at the lower temperature $T_2$. In all cases, error bars are one standard deviation.</figcaption></figure>

cooled to the working temperature $T < T_g$ and relaxes for a waiting time $t_w$ ($T_g$ is the glass temperature, while $t_w$ ranges from minutes to several hours). At time $t_w$ a magnetic field $H$ is switched on and the growing magnetization $M_{ZFC}(t, t_w)$ is recorded at later times $t+t_w$. $M_{ZFC}(t, t_w)$ has turned out to have a

significant dependence on $t_w$ for as long as researchers have had the patience to wait. The relaxation rate

$$S_{ZFC}(t, t_w; H) = \frac{1}{H} \frac{dM_{ZFC}(t, t_w; H)}{d \log t} \quad (1)$$

peaks at a time $t_H^{eff}$ roughly equal to $t_w$ (see, e.g.,
Refs. [9, 10] for experimental results). The only relevant time scale that can be identified is the glass's age, namely $t_w$ (hence the term *aging*). The left panels in Fig. 1 show our results for this comparatively simple fixed-temperature protocol, which will be named *native* hereafter. The native setup is used as a standard for comparison.

## The quest for rejuvenation and memory.

An even more interesting behavior appears when temperature is made to vary with time. In fact, we shall consider here only the simplest protocol for which rejuvenation and memory have been experimentally found [15] (see our temperature-time scheme in the central part of Fig. 1). After a relaxation of duration $t_w^↓$, the temperature is lowered abruptly from the initial temperature $T_1 < T_g$ to a lower temperature $T_2$ (the choice of $T_2$ turns out to be critical, see below). The system is again let to relax at temperature $T_2$ for an additional time $t_w$, after which a magnetic field is switched on and the relaxation function $S_{ZFC}$ is measured at times $t_w^↓ + t_w + t$. Surprisingly enough, one finds that the initial relaxation at $T_1$ has been essentially forgotten: the long-time peak of $S_{ZFC}$ is found at times $t_H^{\text{eff}} \sim t_w$, which can be substantially shorter than $t_w^↓$. This is the rejuvenation effect, which was experimentally found more than 20 years ago and which we are reporting in the bottom-central panel of Fig. 1 for the first time in a simulation.

Yet, rejuvenation is not the end of the story. After a total time of $t_w^↑ = 2t_w^↓$, half of it spent at $T_1$ and half at $T_2$, the system is suddenly heated back to the original temperature $T_1$, where it is left to relax for a time $t_w \ll t_w^↓$, after which the magnetic field is switched on and the relaxation function measured. The $S_{ZFC}$ is found to peak again at time $\sim t_w^↓$, as if the excursion to temperature $T_2$ never happened (Fig. 1, top-right panel). This is the memory effect, which at first sight seems quite contradictory with the rejuvenation effect.

The physical origin of memory and rejuvenation in spin glasses has not been identified yet. Then, it is perhaps unsurprising that all past attempts to reproduce these effects in computer simulations have failed [16–21], which has even raised questions about the validity of the standard model of finite-dimensional spin glasses, the Edwards-Anderson model [22, 23]. Fortunately, the Janus II dedicated supercomputer [24] has changed this situation, attaining realistic time and length scales and allowing for the first time a thorough examination of spin-glass dynamics both in the vicinity of the critical temperature $T_g$ and in the low-temperature regime.

The spin-glass dynamics at $T < T_g$ consists in the growth of (glassy) magnetic domains of linear size $\xi(t_w)$ [25–27] (we shall later refer to this length as $\xi_{\text{micro}}$). The non-equilibrium nature of the process is evident in the growth of $\xi(t_w)$ as $t_w$ varies, which is never-ending and extremely slow. In fact, the lower the temperature, the more sluggish the growth of

$\xi(t_w)$ is, see, e.g., [7, 28]. Janus II has reached unprecedentedly large values of $\xi(t_w)$, enabling safe extrapolations from the numerical time scale of tenths of a second (when $\xi \sim 20a_0$, where $a_0$ is the typical spin-spin distance) to the experimental scale of hours [7, 28] (when $\xi \sim 200a_0$). This special-purpose computer has also made it possible to simulate [29] the experimental protocol for extracting the spin-glass coherence length from the Zeeman effect [26], thus showing the consistency between the Zeeman method and the microscopic approach. Janus II allowed us to perform computer experiments with a native (i.e., fixed-temperature) protocol and make a direct comparison of the $S_{ZFC}$ (1) obtained in the simulation with that from real experiments on a single crystal of CuMn [9, 10, 30]. The Edwards-Anderson model and CuMn turned out to be governed by the same scaling laws, where $\xi$ is the all-important scaling variable. This agreement between simulations and experiment, however, was established only for native protocols. We need to understand what happens when temperature is varied.

Experimentalists are prone to attribute the rejuvenation effect to temperature chaos (see, e.g., Ref. [15]); explanations not invoking temperature chaos have been also proposed [31, 32]). Temperature chaos [33–35] is an equilibrium notion stating that spin configurations typical from the Boltzmann distribution at temperature $T_1$ would be very atypical for temperature $T_2$, no matter how close $T_1$ and $T_2$ are (provided that $T_1, T_2 < T_g$). Temperature chaos could explain why the relaxation at temperature $T_1$ seems useless at $T_2$ (i.e., rejuvenation). Yet, even in the mean-field approximation, showing that temperature chaos is really present in equilibrium has been a real tour-de-force [36, 37]. Furthermore, extending the equilibrium concept of temperature chaos to the experimentally relevant context of off-equilibrium dynamics is a very recent achievement [38].

Dynamic temperature chaos is spatially extremely heterogeneous (see Fig. 2). To measure it we choose many spheres of linear size $R$ in random positions within the sample. We compare within each sphere spin configurations obtained at temperature $T_1$ and time $t_w^{T_1}$ with configurations from temperature $T_2$ and time $t_w^{T_2}$ (the simplifying choice $\xi(t_w^{T_1}, T_1) = \xi(t_w^{T_2}, T_2) = \xi(t_w)$ was made in [38]). The comparison is quantitative, through the computation of a correlation coefficient $X_{T_1,T_2}$, see **Methods**. Many of those spheres turn out to have very weak temperature chaos [$X_{T_1,T_2} \approx 1$]. Yet, with low probability, one picks a chaotic sphere with a significantly smaller $X_{T_1,T_2}$. In fact, the analysis in [38] identifies a crossover length scale $\xi^*(T_1, T_2)$: for $\xi(t_w) \ll \xi^*(T_1, T_2)$ chaotic spheres are *very rare* but for $\xi(t_w) \gg \xi^*(T_1, T_2)$ chaotic spheres become *fairly typical*. A scaling law was also found: $\xi^*(T_1, T_2) \propto (T_1 - T_2)^{-1/\zeta_{\text{NE}}}$ with $\zeta_{\text{NE}} = 1.19(2)$.

Our last building block comes from the experiment of Ref. [8], which identifies a minimal temperature jump $\Delta T_{\text{min}}$ in a CuMn sample. Temperature chaos in
Refs. [9, 10] for experimental results). The only relevant time scale that can be identified is the glass's age, namely $t_w$ (hence the term *aging*). The left panels in Fig. 1 show our results for this comparatively simple fixed-temperature protocol, which will be named *native* hereafter. The native setup is used as a standard for comparison.

## The quest for rejuvenation and memory.

An even more interesting behavior appears when temperature is made to vary with time. In fact, we shall consider here only the simplest protocol for which rejuvenation and memory have been experimentally found [15] (see our temperature-time scheme in the central part of Fig. 1). After a relaxation of duration $t_w^↓$, the temperature is lowered abruptly from the initial temperature $T_1 < T_g$ to a lower temperature $T_2$ (the choice of $T_2$ turns out to be critical, see below). The system is again let to relax at temperature $T_2$ for an additional time $t_w$, after which a magnetic field is switched on and the relaxation function $S_{ZFC}$ is measured at times $t_w^↓ + t_w + t$. Surprisingly enough, one finds that the initial relaxation at $T_1$ has been essentially forgotten: the long-time peak of $S_{ZFC}$ is found at times $t_H^{\text{eff}} \sim t_w$, which can be substantially shorter than $t_w^↓$. This is the rejuvenation effect, which was experimentally found more than 20 years ago and which we are reporting in the bottom-central panel of Fig. 1 for the first time in a simulation.

Yet, rejuvenation is not the end of the story. After a total time of $t_w^↑ = 2t_w^↓$, half of it spent at $T_1$ and half at $T_2$, the system is suddenly heated back to the original temperature $T_1$, where it is left to relax for a time $t_w \ll t_w^↓$, after which the magnetic field is switched on and the relaxation function measured. The $S_{ZFC}$ is found to peak again at time $\sim t_w^↓$, as if the excursion to temperature $T_2$ never happened (Fig. 1, top-right panel). This is the memory effect, which at first sight seems quite contradictory with the rejuvenation effect.

The physical origin of memory and rejuvenation in spin glasses has not been identified yet. Then, it is perhaps unsurprising that all past attempts to reproduce these effects in computer simulations have failed [16–21], which has even raised questions about the validity of the standard model of finite-dimensional spin glasses, the Edwards-Anderson model [22, 23]. Fortunately, the Janus II dedicated supercomputer [24] has changed this situation, attaining realistic time and length scales and allowing for the first time a thorough examination of spin-glass dynamics both in the vicinity of the critical temperature $T_g$ and in the low-temperature regime.

The spin-glass dynamics at $T < T_g$ consists in the growth of (glassy) magnetic domains of linear size $\xi(t_w)$ [25–27] (we shall later refer to this length as $\xi_{\text{micro}}$). The non-equilibrium nature of the process is evident in the growth of $\xi(t_w)$ as $t_w$ varies, which is never-ending and extremely slow. In fact, the lower the temperature, the more sluggish the growth of

$\xi(t_w)$ is, see, e.g., [7, 28]. Janus II has reached unprecedentedly large values of $\xi(t_w)$, enabling safe extrapolations from the numerical time scale of tenths of a second (when $\xi \sim 20a_0$, where $a_0$ is the typical spin-spin distance) to the experimental scale of hours [7, 28] (when $\xi \sim 200a_0$). This special-purpose computer has also made it possible to simulate [29] the experimental protocol for extracting the spin-glass coherence length from the Zeeman effect [26], thus showing the consistency between the Zeeman method and the microscopic approach. Janus II allowed us to perform computer experiments with a native (i.e., fixed-temperature) protocol and make a direct comparison of the $S_{ZFC}$ (1) obtained in the simulation with that from real experiments on a single crystal of CuMn [9, 10, 30]. The Edwards-Anderson model and CuMn turned out to be governed by the same scaling laws, where $\xi$ is the all-important scaling variable. This agreement between simulations and experiment, however, was established only for native protocols. We need to understand what happens when temperature is varied.

Experimentalists are prone to attribute the rejuvenation effect to temperature chaos (see, e.g., Ref. [15]); explanations not invoking temperature chaos have been also proposed [31, 32]). Temperature chaos [33–35] is an equilibrium notion stating that spin configurations typical from the Boltzmann distribution at temperature $T_1$ would be very atypical for temperature $T_2$, no matter how close $T_1$ and $T_2$ are (provided that $T_1, T_2 < T_g$). Temperature chaos could explain why the relaxation at temperature $T_1$ seems useless at $T_2$ (i.e., rejuvenation). Yet, even in the mean-field approximation, showing that temperature chaos is really present in equilibrium has been a real tour-de-force [36, 37]. Furthermore, extending the equilibrium concept of temperature chaos to the experimentally relevant context of off-equilibrium dynamics is a very recent achievement [38].

Dynamic temperature chaos is spatially extremely heterogeneous (see Fig. 2). To measure it we choose many spheres of linear size $R$ in random positions within the sample. We compare within each sphere spin configurations obtained at temperature $T_1$ and time $t_w^{T_1}$ with configurations from temperature $T_2$ and time $t_w^{T_2}$ (the simplifying choice $\xi(t_w^{T_1}, T_1) = \xi(t_w^{T_2}, T_2) = \xi(t_w)$ was made in [38]). The comparison is quantitative, through the computation of a correlation coefficient $X_{T_1,T_2}$, see Methods. Many of those spheres turn out to have very weak temperature chaos [$X_{T_1,T_2} \approx 1$]. Yet, with low probability, one picks a chaotic sphere with a significantly smaller $X_{T_1,T_2}$. In fact, the analysis in [38] identifies a crossover length scale $\xi^*(T_1, T_2)$: for $\xi(t_w) \ll \xi^*(T_1, T_2)$ chaotic spheres are very rare but for $\xi(t_w) \gg \xi^*(T_1, T_2)$ chaotic spheres become fairly typical. A scaling law was also found: $\xi^*(T_1, T_2) \propto (T_1 - T_2)^{-1/\zeta_{\text{NE}}}$ with $\zeta_{\text{NE}} = 1.19(2)$.

Our last building block comes from the experiment of Ref. [8], which identifies a minimal temperature jump $\Delta T_{\text{min}}$ in a CuMn sample. Temperature chaos in
<figure><img src="image_4.png" /><figcaption>FIG. 2. Temperature chaos is spatially heterogeneous when it is clearly present. The 8000 randomly chosen spheres in a sample of size $L = 160$ are depicted with a color code depending on $1-X$ ($X$ is the chaotic correlation parameter as computed for spheres of radius $R = 5a_0$, see Methods). For visualization purposes, spheres are represented with a radius $12(1-X)$, so that only fully chaotic spheres (i.e., $X=0$) would have the largest size. In order to avoid cluttering, we draw only spheres with $X < 0.97$. On the left, we calculate the chaotic correlation parameter $X$ between the native system at $T = 0.5$ (i.e., a fixed-temperature protocol: a completely disordered system is put at temperature $T = 0.5$ and let to evolve at this temperature for a time $t_w = 2^{31.25}$) and the jump system at the same temperature $T = 0.5$ (recall the central part of Fig. 1: the jump system has spent the first half of its life, $t_w^{\uparrow} = 2^{31.25}$, at the hot temperature $T_1 = 0.9$ and the second half, $t_w = 2^{31.25}$, at the cold temperature $T = 0.5$). Very strong chaotic heterogeneity is found. On the right, we calculate the chaotic correlation parameter $X$ between the native system at the hot temperature $T = 0.9$ and the jump-back system at the same temperature $T = 0.9$ (the jump-back system has spent a time $t_w^{\downarrow} = 2^{31.25}$ at the hot temperature $T_1 = 0.9$, a time $t_w^{\uparrow} - t_w^{\downarrow} = 2^{31.25}$ at the cold temperature $T_2 = 0.5$, and then $t_w = 2^{10}$ again at $T_1 = 0.9$ — see the temperature protocol of Fig. 1). After the cycle the system does not display chaotic heterogeneity since almost every sphere has a large correlation parameter $X$, i.e., a strong memory (more examples can be found in Supplementary Note V).</figcaption></figure>

that sample turned out to be exceedingly weak whenever $T_1 - T_2 < \Delta T_{min}$. It follows that, in a simulation, chaotic spheres will be just too rare to significantly affect the overall sample relaxation unless [39]

$$ \left. \frac{T_1 - T_2}{T_g} \right|_{\text{sim}} \approx \left. \frac{\Delta T_{\text{min}}}{T_g} \right|_{\text{CuMn}} \left[ \frac{\xi_{\text{CuMn}}(t_w)}{\xi_{\text{micro}}(t_w)} \right]^{\zeta_{\text{NE}}} , \quad (2) $$

where the subindex micro stands for the $\xi$ computed in the numerical simulation (see Methods) while $T_g$ is the glass temperature, which is different for the CuMn sample and for simulations. Plugging in typical numbers ($\Delta T_{min} = 450$ mK, $T_g = 31.5$ K, $\xi_{CuMn}(t_w) \approx 220 a_0$ and $\xi_{micro}(t_w) \approx 16.6 a_0$), we conclude from Eq. (2) that, given the correlation length reached in our simulations, a successful simulation of the rejuvenation effect should have $T_1 - T_2 > 0.32 T_g$.

In this work we have considered two temperature jumps, see Table I. The first jump, namely $T_1 = 0.9 \rightarrow T_2 = 0.5$, meets the requirement for temperature chaos expressed in Eq. (2), while the second jump, $T_1 = 0.9 \rightarrow T_2 = 0.7$ is too small. Hence we expect to find qualitative differences between the two.

## Becoming quantitative: how many controlling length scales?

Our discussion shall emphasize three different length scales, focusing on their physical interpretation and their utility to rationalize the rejuvenation and memory effects (many more details are provided in Methods [40]). Only one of these scales, named $\xi_{\text{Zeeman}}$, can be experimentally accessed nowadays (the other two lengths, however, provide invaluable microscopic information):

* $\xi_{\text{micro}}$ is the size of the (glassy) domains within the sample (is the largest length scale at which we can regard the system as ordered at time $t_w$).

* $\xi_{\text{Zeeman}}$ is obtained by counting the number of spins that react coherently to an externally applied field [26]. It provides a very direct quantification of memory and rejuvenation.

* $\zeta(t_1, t_2)$ [41–43] is obtained from the comparison of the same system at the two times $t_1 < t_2$ [$\zeta(t_1, t_2)$ is the typical size of the regions where coherent rearrangements have occurred between $t_1$ and $t_2$].

Previous analysis for native (i.e., fixed-
<figure><img src="image_3.png" /><figcaption>FIG. 2. Temperature chaos is spatially heterogeneous when it is clearly present. The 8000 randomly chosen spheres in a sample of size L = 160 are depicted with a color code depending on 1 − X (X is the chaotic correlation parameter as computed for spheres of radius R = 5 a$_{0}$ , see Methods). For visualization purposes, spheres are represented with a radius 12(1 − X ), so that only fully chaotic spheres (i.e., X = 0) would have the largest size. In order to avoid cluttering, we draw only spheres with X < 0 . 97. On the left, we calculate the chaotic correlation parameter X between the native system at T = 0 . 5 (i.e., a fixed-temperature protocol: a completely disordered system is put at temperature T = 0 . 5 and let to evolve at this temperature for a time t$_{w}$ = 2 31 . $^{25}$) and the jump system at the same temperature T = 0 . 5 (recall the central part of Fig. 1: the jump system has spent the first half of its life, t w = 2 31 . $^{25}$, at the hot temperature T$_{1}$ = 0 . 9 and the second half, t$_{w}$ = 2 31 . $^{25}$, at the cold temperature T = 0 . 5). Very strong chaotic heterogeneity is found. On the right, we calculate the chaotic correlation parameter X between the native system at the hot temperature T = 0 . 9 and the jump-back system at the same temperature T = 0 . 9 (the jump-back system has spent a time t ↓ w = 2 31 . 25 at the hot temperature T$_{1}$ = 0 . 9, a time t ↑ w − t ↓ w = 2 31 . 25 at the cold temperature T$_{2}$ = 0 . 5, and then t$_{w}$ = 2 10 again at T$_{1}$ = 0 . 9 — see the temperature protocol of Fig. 1). After the cycle the system does not display chaotic heterogeneity since almost every sphere has a large correlation parameter X , i.e., a strong memory (more examples can be found in Supplementary Note V).</figcaption></figure>

that sample turned out to be exceedingly weak when- ever T$_{1}$ − T$_{2}$ < Δ T$_{min}$ . It follows that, in a simulation, chaotic spheres will be just too rare to significantly af- fect the overall sample relaxation unless [39]

$$
\left. \frac{T_1 - T_2}{T_g} \right|_{\text{sim}} \approx \left. \frac{\Delta T_{\text{min}}}{T_g} \right|_{\text{CuMn}} \left[ \frac{\xi_{\text{CuMn}}(t_w)}{\xi_{\text{micro}}(t_w)} \right]^{\zeta_{\text{NE}}}, \quad (2)
$$

where the subindex micro stands for the ξ computed in the numerical simulation (see Methods) while T$_{g}$ is the glass temperature, which is different for the CuMn sample and for simulations. Plugging in typical numbers (Δ T$_{min}$ = 450 mK, T$_{g}$ = 31 . 5 K, ξ$_{CuMn}$ ( t$_{w}$ ) ≈ 220 a$_{0}$ and ξ$_{micro}$ ( t$_{w}$ ) ≈ 16 . 6 a$_{0}$ ), we con- clude from Eq. (2) that, given the correlation length reached in our simulations, a successful simulation of the rejuvenation effect should have T$_{1}$ − T$_{2}$ > 0 . 32 T$_{g}$ .

In this work we have considered two temperature jumps, see Table I. The first jump, namely T$_{1}$ = 0 . 9 → T$_{2}$ = 0 . 5, meets the requirement for tempera- ture chaos expressed in Eq. (2), while the second jump, T$_{1}$ = 0 . 9 → T$_{2}$ = 0 . 7 is too small. Hence we expect to find qualitative differences between the two.

Becoming quantitative: how many controlling length scales?

Our discussion shall emphasize three different length scales, focusing on their physical interpreta- tion and their utility to rationalize the rejuvenation and memory effects (many more details are provided in Methods [40]). Only one of these scales, named ξ$_{Zeeman}$ , can be experimentally accessed nowadays (the other two lengths, however, provide invaluable microscopic information):

• ξ$_{micro}$ is the size of the (glassy) domains within the sample (is the largest length scale at which we can regard the system as ordered at time t$_{w}$ ).

• ξ$_{Zeeman}$ is obtained by counting the number of spins that react coherently to an externally ap- plied field [26]. It provides a very direct quan- tification of memory and rejuvenation.

• ζ ( t$_{1}$, t$_{2}$ ) [41–43] is obtained from the comparison of the same system at the two times t$_{1}$ < t$_{2}$ [ ζ ( t$_{1}$, t$_{2}$ ) is the typical size of the regions where coherent rearrangements have occurred between t$_{1}$ and t$_{2}$ ].

Previous analysis for native ( i.e. , fixed-
temperature) protocols tell us that $\xi_{\text{Zeeman}}$ follows quite closely the behavior of the microscopic length $\xi_{\text{micro}}$ [9, 10, 26, 29]. This is what we find in the top panel of Fig. 3. There are two salient features in the time growth of either $\xi_{\text{Zeeman}}$ or $\xi_{\text{micro}}(t_w)$ at fixed temperature [7, 28]: the growth slows down as $\xi_{\text{micro}}$ increases [44] and the dynamics at lower temperatures is enormously slower [45]. In fact, see Table I and Ref. [28], at the largest temperature $T = 0.9$ it is comparatively easy to reach a large $\xi_{\text{micro}} \approx 16.6 a_0$ in a native protocol. Instead, for a similar simulation time, the native protocol at $T = 0.5$ is limited to $\xi_{\text{micro}} \approx 5.6 a_0$. It is then unsurprising that, when the temperature jumps from $T_1 = 0.9$ to $T_2 = 0.5$ or 0.7, see Fig. 3-top, the size of the glassy domains is locked to their value at jump time, namely $\xi_{\text{micro}} \approx 16.6 a_0$: the time needed for such a large domain to grow at the lower temperature $T_2$ far exceeds the scale of our simulations. The importance of this locking was also emphasized in Ref. [18].

While $\xi_{\text{micro}}$ is locked at the value it has at the jump time, the behavior of $\xi_{\text{Zeeman}}$ is different in the jump protocols. In the jump complying with Eq. (2), $T_1 = 0.9 \rightarrow T_2 = 0.5$, $\xi_{\text{Zeeman}}(t_w)$ is quite similar to the corresponding curve for the native run at $T = 0.5$. From the point of view of the response to the magnetic field, rejuvenation is almost complete for this temperature jump, because the initial relaxation at $T_1 = 0.9$ (almost) does not leave a measurable trace. Instead, for the more modest jump $T_1 = 0.9 \rightarrow T_2 = 0.7$, rejuvenation is weaker and $\xi_{\text{Zeeman}}$ is sensibly larger than in the native runs (see Supplementary Note II for more details).

Furthermore, it is also shown in Fig. 3-top that, when the system jumps back to $T_1 = 0.9$ (i.e. $T_1 = 0.9 \rightarrow T_2 = 0.5 \rightarrow T_1 = 0.9$), recall the top-right panel in Fig. 1), the response to the magnetic-field goes back to normal: $\xi_{\text{Zeeman}}$ catches up with $\xi_{\text{micro}}$ after a extremely short transient. This is another manifestation of the memory effect.

As for the third length scale, see Fig. 3-bottom, for all our jump protocols we find $\zeta \ll \xi_{\text{micro}}$, which means that the configuration right before the jump has not been substantially modified by the excursion to the low temperature $T_2 = 0.5$. Interestingly enough, for the jump-back protocol $\zeta$ gets a lot bigger than for the single jump protocols. In other words, the dynamics is unlocked when the system comes back to its original temperature $T_1$. Nevertheless, $\zeta$ is still substantially smaller than $\xi_{\text{micro}}$. The original spin configuration has, therefore, suffered only local distortions.

## Dynamic temperature chaos and rejuvenation

At this point, the elephant in the room is clear: what is the physical origin for rejuvenation and memory?

In order to answer the question, we need to compare pairs of spin configurations. One of the config-

<figure><img src="image_5.png" /><figcaption>FIG. 3. Aging dynamics is controlled by three length scales (at least). (See text and Methods for extended discussion of these three length scales). The dashed lines and empty dots are for $\xi_{\text{micro}}(t_w, T)$, the continuous lines and filled dots are for $\xi_{\text{Zeeman}}(t_w, T)$, and the continuous lines and filled triangles are for $\zeta(t_1, t_2)$. On the top, for the native protocols, $\xi_{\text{Zeeman}}(t_w)$ follows quite closely the behavior of $\xi_{\text{micro}}(t_w)$. For the jump protocol with $T_2 = 0.5$, $\xi_{\text{Zeeman}}(t_w)$ is extremely similar to the corresponding curve for the native run [this $T_2$ meets the chaos requirement expressed in Eq. (2)], which means that the system responds to an external magnetic field as if rejuvenated, $\xi_{\text{Zeeman}}^{\text{jump}}(t_w) \ll \xi_{\text{micro}}(t_w)$. When the system jumps back to $T_1 = 0.9$ (i.e. $T_1 = 0.9 \rightarrow T_2 = 0.5 \rightarrow T_1 = 0.9$), $\xi_{\text{Zeeman}}^{\text{jump}}(t_w) = \xi_{\text{micro}}(t_w)$ after a extremely short time (memory). Instead, $\xi_{\text{Zeeman}}(t_w)$ never becomes small for the jump protocols with $T_2 = 0.7$. On the bottom, the size of the regions undergoing coherent rearrangements when evolving from the initial to the final time, $\zeta(t_1, t_2)$, is much smaller than $\xi_{\text{micro}}(t_w)$ for all our jump protocols. The earliest of the two times is the very last instant that the system spent at temperature $T_1 = 0.9$ (i.e., just before the temperature jump; in the case of the jump-back protocol just before the first jump $T_1 = 0.9 \rightarrow T_2 = 0.5$). In all cases, $\zeta$ is represented as a function of $t_2 - t_1$. Error bars are one standard deviation. See Supplementary Note III for more results on $\zeta(t_1, t_2)$.</figcaption></figure>

urations will be taken from the jump protocols. The other configuration will come from the native runs at temperatures $T_2 = 0.5$ or $T_2 = 0.7$. In an attempt to make a fair comparison, we shall choose the native configurations at $T_2$ at their largest possible waiting time. In fact, the magnetic domains will be substantially smaller in the native protocol than they are in the jump protocol (at $T_2 = 0.5$, for instance, one has to compare $\xi_{\text{micro}}^{\text{native}} \approx 5.8 a_0$ with $\xi_{\text{micro}}^{\text{jump}} \approx 16.6 a_0$).

The main steps in the comparison were outlined above (for a more paused exposition see Methods and Ref. [38]). We pick at random in the sample spheres
temperature) protocols tell us that $\xi_{\text{Zeeman}}$ follows quite closely the behavior of the microscopic length $\xi_{\text{micro}}$ [9, 10, 26, 29]. This is what we find in the top panel of Fig. 3. There are two salient features in the time growth of either $\xi_{\text{Zeeman}}$ or $\xi_{\text{micro}}(t_w)$ at fixed temperature [7, 28]: the growth slows down as $\xi_{\text{micro}}$ increases [44] and the dynamics at lower temperatures is enormously slower [45]. In fact, see Table I and Ref. [28], at the largest temperature $T = 0.9$ it is comparatively easy to reach a large $\xi_{\text{micro}} \approx 16.6 a_0$ in a native protocol. Instead, for a similar simulation time, the native protocol at $T = 0.5$ is limited to $\xi_{\text{micro}} \approx 5.6 a_0$. It is then unsurprising that, when the temperature jumps from $T_1 = 0.9$ to $T_2 = 0.5$ or 0.7, see Fig. 3-top, the size of the glassy domains is locked to their value at jump time, namely $\xi_{\text{micro}} \approx 16.6 a_0$: the time needed for such a large domain to grow at the lower temperature $T_2$ far exceeds the scale of our simulations. The importance of this locking was also emphasized in Ref. [18].

While $\xi_{\text{micro}}$ is locked at the value it has at the jump time, the behavior of $\xi_{\text{Zeeman}}$ is different in the jump protocols. In the jump complying with Eq. (2), $T_1 = 0.9 \rightarrow T_2 = 0.5$, $\xi_{\text{Zeeman}}(t_w)$ is quite similar to the corresponding curve for the native run at $T = 0.5$. From the point of view of the response to the magnetic field, rejuvenation is almost complete for this temperature jump, because the initial relaxation at $T_1 = 0.9$ (almost) does not leave a measurable trace. Instead, for the more modest jump $T_1 = 0.9 \rightarrow T_2 = 0.7$, rejuvenation is weaker and $\xi_{\text{Zeeman}}$ is sensibly larger than in the native runs (see Supplementary Note II for more details).

Furthermore, it is also shown in Fig. 3-top that, when the system jumps back to $T_1 = 0.9$ (i.e. $T_1 = 0.9 \rightarrow T_2 = 0.5 \rightarrow T_1 = 0.9$), recall the top-right panel in Fig. 1), the response to the magnetic-field goes back to normal: $\xi_{\text{Zeeman}}$ catches up with $\xi_{\text{micro}}$ after a extremely short transient. This is another manifestation of the memory effect.

As for the third length scale, see Fig. 3-bottom, for all our jump protocols we find $\zeta \ll \xi_{\text{micro}}$, which means that the configuration right before the jump has not been substantially modified by the excursion to the low temperature $T_2 = 0.5$. Interestingly enough, for the jump-back protocol $\zeta$ gets a lot bigger than for the single jump protocols. In other words, the dynamics is unlocked when the system comes back to its original temperature $T_1$. Nevertheless, $\zeta$ is still substantially smaller than $\xi_{\text{micro}}$. The original spin configuration has, therefore, suffered only local distortions.

## Dynamic temperature chaos and rejuvenation

At this point, the elephant in the room is clear: what is the physical origin for rejuvenation and memory?

In order to answer the question, we need to compare pairs of spin configurations. One of the config-

<figure><img src="image_4.png" /><figcaption>FIG. 3. Aging dynamics is controlled by three length scales (at least). (See text and Methods for extended discussion of these three length scales). The dashed lines and empty dots are for $\xi_{\text{micro}}(t_w, T)$, the continuous lines and filled dots are for $\xi_{\text{Zeeman}}(t_w, T)$, and the continuous lines and filled triangles are for $\zeta(t_1, t_2)$. On the top, for the native protocols, $\xi_{\text{Zeeman}}(t_w)$ follows quite closely the behavior of $\xi_{\text{micro}}(t_w)$. For the jump protocol with $T_2 = 0.5$, $\xi_{\text{Zeeman}}(t_w)$ is extremely similar to the corresponding curve for the native run [this $T_2$ meets the chaos requirement expressed in Eq. (2)], which means that the system responds to an external magnetic field as if rejuvenated, $\xi_{\text{Zeeman}}^{\text{jump}}(t_w) \ll \xi_{\text{micro}}(t_w)$. When the system jumps back to $T_1 = 0.9$ (i.e. $T_1 = 0.9 \rightarrow T_2 = 0.5 \rightarrow T_1 = 0.9$), $\xi_{\text{Zeeman}}^{\text{jump}}(t_w) = \xi_{\text{micro}}(t_w)$ after a extremely short time (memory). Instead, $\xi_{\text{Zeeman}}(t_w)$ never becomes small for the jump protocols with $T_2 = 0.7$. On the bottom, the size of the regions undergoing coherent rearrangements when evolving from the initial to the final time, $\zeta(t_1, t_2)$, is much smaller than $\xi_{\text{micro}}(t_w)$ for all our jump protocols. The earliest of the two times is the very last instant that the system spent at temperature $T_1 = 0.9$ (i.e., just before the temperature jump; in the case of the jump-back protocol just before the first jump $T_1 = 0.9 \rightarrow T_2 = 0.5$). In all cases, $\zeta$ is represented as a function of $t_2 - t_1$. Error bars are one standard deviation. See Supplementary Note III for more results on $\zeta(t_1, t_2)$.</figcaption></figure>

urations will be taken from the jump protocols. The other configuration will come from the native runs at temperatures $T_2 = 0.5$ or $T_2 = 0.7$. In an attempt to make a fair comparison, we shall choose the native configurations at $T_2$ at their largest possible waiting time. In fact, the magnetic domains will be substantially smaller in the native protocol than they are in the jump protocol (at $T_2 = 0.5$, for instance, one has to compare $\xi_{\text{micro}}^{\text{native}} \approx 5.8 a_0$ with $\xi_{\text{micro}}^{\text{jump}} \approx 16.6 a_0$).

The main steps in the comparison were outlined above (for a more paused exposition see Methods and Ref. [38]). We pick at random in the sample spheres
<table><caption>TABLE I. Basic features of our simulations. We have simulated on the Janus II supercomputer the Edwards-Anderson model with nearest-neighbor couplings ($J = \pm 1$ with 50% probability), on simple-cubic lattices containing $160^3$ Ising spins $s = \pm 1$ (the lattice size is $L = 160 a_0$) and endowed with periodic boundary conditions. A particular set of couplings is termed *sample*. For every sample and every set of parameters, we have simulated 512 independent trajectories (i.e., 512 replicas, see Methods). This table lists the main parameters for each of our numerical simulations. Temperature-varying protocols, see the central part of Fig. 1, are named *jump protocols*. In all cases, temperature $T$ refers to the temperature at which the relaxation function in Eq. (1) is computed. All temperatures considered are in the spin-glass phase: $T < T_g = 1.102(3)$ [13]. The waiting time is the period before the magnetic field $H$ is switched on (for jump protocols, this consists of a time $t_w^↓$ at the starting temperature $T_1 = 0.9$, plus a time $t_w$ at $T_2$). The microscopic correlation length $\xi_{micro}$ is given in $a_0 = 1$ units and computed just before the magnetic field is switched on [error bars for $\xi_{micro}(t_w)$ are one standard deviation]. Finally, the longest simulation time in the presence of a field is $t_{max}$.</caption><thead><tr><th>Type</th><th>T</th><th>Waiting time</th><th>ξ<sub>micro</sub></th><th>t<sub>max</sub></th></tr></thead><tbody><tr><td>NATIVE</td><td>0.9</td><td>2<sup>31.25</sup> (= t<sub>w</sub><sup>↓</sup>)</td><td>16.63(5)</td><td>2<sup>32</sup></td></tr><tr><td>NATIVE</td><td>0.5</td><td>2<sup>10</sup></td><td>2.23926(2)</td><td>2<sup>28</sup></td></tr><tr><td>NATIVE</td><td>0.5</td><td>2<sup>15.625</sup></td><td>2.9090(4)</td><td>2<sup>28</sup></td></tr><tr><td>NATIVE</td><td>0.5</td><td>2<sup>23.5</sup></td><td>4.0865(15)</td><td>2<sup>30</sup></td></tr><tr><td>NATIVE</td><td>0.5</td><td>2<sup>31.25</sup></td><td>5.6167(4)</td><td>2<sup>32</sup></td></tr><tr><td>JUMP</td><td>0.5</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>10</sup></td><td>16.62(12)</td><td>2<sup>28</sup></td></tr><tr><td>JUMP</td><td>0.5</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>15.625</sup></td><td>16.68(12)</td><td>2<sup>28</sup></td></tr><tr><td>JUMP</td><td>0.5</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>23.5</sup></td><td>16.75(13)</td><td>2<sup>31</sup></td></tr><tr><td>JUMP</td><td>0.5</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td><td>16.81(13)</td><td>2<sup>33.5</sup></td></tr><tr><td>NATIVE</td><td>0.7</td><td>2<sup>10</sup></td><td>2.6629(4)</td><td>2<sup>28</sup></td></tr><tr><td>NATIVE</td><td>0.7</td><td>2<sup>15.625</sup></td><td>3.8230(10)</td><td>2<sup>28</sup></td></tr><tr><td>NATIVE</td><td>0.7</td><td>2<sup>23.5</sup></td><td>6.1742(4)</td><td>2<sup>28</sup></td></tr><tr><td>NATIVE</td><td>0.7</td><td>2<sup>31.25</sup></td><td>9.578(11)</td><td>2<sup>33</sup></td></tr><tr><td>JUMP</td><td>0.7</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>10</sup></td><td>16.62(12)</td><td>2<sup>28</sup></td></tr><tr><td>JUMP</td><td>0.7</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>15.625</sup></td><td>16.67(12)</td><td>2<sup>28</sup></td></tr><tr><td>JUMP</td><td>0.7</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>23.5</sup></td><td>16.76(12)</td><td>2<sup>28</sup></td></tr><tr><td>JUMP</td><td>0.7</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td><td>16.81(13)</td><td>2<sup>32</sup></td></tr></tbody></table>

of radius $R$. The results presented in this paper were obtained with $R = 5 a_0$ to make sure that the spheres will have a chance to fit within the glassy domains of the native runs (we have tried other values of $R$, finding qualitatively similar results, see Supplementary Note VI). The configurations from the two protocols are compared by computing a correlation coefficient $X$ that takes into account only the spins contained in the sphere. If $X$ is significantly smaller than unity we regard that particular sphere as *chaotic*, because typical configurations from the two protocols differ within the sphere. To be precise, we compute the probability distribution function $F(\tilde{X})$, namely the fraction of the spheres with a correlation coefficient $X < \tilde{X}$.

Our results shown in Fig. 4-bottom for the jump protocol $T_1 = 0.9 \rightarrow T_2 = 0.7$ remind us of previous studies [38]. The vast majority of the spheres have a very large correlation coefficient, and truly chaotic

spheres are found only in the tail of the distribution (probability 0.1% or smaller).

Interestingly enough, see the left panel in Fig. 2 and Fig. 4-top, the situation is radically different for the jump protocol $T_1 = 0.9 \rightarrow T_2 = 0.5$, where the spheres in percentile 10 of the distribution are as chaotic as the most chaotic spheres we could find for the jump $T_1 = 0.9 \rightarrow T_2 = 0.7$. In fact, to our knowledge, Fig. 4-top reports the strongest temperature-chaos signal ever observed in a simulation of glassy dynamics.

In order to convince ourselves that the extreme chaos is not an artifact of the disparity in domain sizes, we have tried a null experiment by simulating a model where no temperature chaos is expected, namely the link-diluted ferromagnetic Ising model (we have used the results in Ref. [46] to match as closely as possible in the diluted ferromagnet the conditions in our spin-glass simulations, with special care in matching the size of the domains, see Methods). As expected, see Fig. 4, the sphere distribution for the ferromagnet is concentrated at correlation coefficient $X \approx 1$. We conclude that the spin-glass results in Fig. 4-top are genuine evidence for dynamic temperature chaos.

It is also interesting that the distribution function in Fig. 4-top barely depends on $t_w$. This is another manifestation of the dynamic lock-down when the temperature jumps to the lower value.

The overall conclusions of this analysis are twofold. First, the requirement expressed by Eq. (2), which is based on CuMn experimental results [8], is sensible: strong temperature chaos is found only when $T_1 - T_2$ is as large as Eq. (2) demands. Second, only when temperature chaos is strong do our simulations find strong rejuvenation (recall Fig. 3-top).

### Where do we stand?

Our simulations depict a clear picture of the rejuvenation and memory effects. Provided that the temperature jump is large enough, see Eq. (2), the spin-glass state that was forming at temperature $T_1$ is completely alien at temperature $T_2$ (at least it looks like an alien when compared with the native state that grows directly at $T_2$, see Fig. 2). In fact, the response to the magnetic field (which is the quantity measured in experiments [2–5, 7–10, 15]) is not qualitatively different in the alien state and in the native state that grows from a fully disordered high-temperature state. The system just dismisses the relaxation it achieved at the higher temperature $T_1$.

Paradoxically enough, the alien state is locked at temperature $T_2$: the microscopic rearrangement at $T_2$, see Fig. 3-bottom, takes place on too small length scales to dissolve such foreign glassy domains. As a consequence, when the temperature is taken back to $T_1$, the glassy domains characteristic of $T_1$ are still there. This seems to be the physical origin of the memory effect. This reasoning is also consistent with recent experiments that find that the memory effect strongly depends on $t_w^↓$ (i.e., the time spent in the
<table><caption>TABLE I. Basic features of our simulations. We have simulated on the Janus II supercomputer the Edwards-Anderson model with nearest-neighbor couplings ($J = \pm 1$ with 50% probability), on simple-cubic lattices containing $160^3$ Ising spins $s = \pm 1$ (the lattice size is $L = 160 a_0$) and endowed with periodic boundary conditions. A particular set of couplings is termed *sample*. For every sample and every set of parameters, we have simulated 512 independent trajectories (i.e., 512 replicas, see Methods). This table lists the main parameters for each of our numerical simulations. Temperature-varying protocols, see the central part of Fig. 1, are named *jump protocols*. In all cases, temperature $T$ refers to the temperature at which the relaxation function in Eq. (1) is computed. All temperatures considered are in the spin-glass phase: $T < T_g = 1.102(3)$ [13]. The waiting time is the period before the magnetic field $H$ is switched on (for jump protocols, this consists of a time $t_w^↓$ at the starting temperature $T_1 = 0.9$, plus a time $t_w$ at $T_2$). The microscopic correlation length $\xi_{micro}$ is given in $a_0 = 1$ units and computed just before the magnetic field is switched on [error bars for $\xi_{micro}(t_w)$ are one standard deviation]. Finally, the longest simulation time in the presence of a field is $t_{max}$.</caption><thead><tr><th>Type</th><th>T</th><th>Waiting time</th><th>ξ<sub>micro</sub></th><th>t<sub>max</sub></th></tr></thead><tbody><tr><td>NATIVE</td><td>0.9</td><td>2<sup>31.25</sup> (= t<sub>w</sub><sup>↓</sup>)</td><td>16.63(5)</td><td>2<sup>32</sup></td></tr><tr><td>NATIVE</td><td>0.5</td><td>2<sup>10</sup></td><td>2.23926(2)</td><td>2<sup>28</sup></td></tr><tr><td>NATIVE</td><td>0.5</td><td>2<sup>15.625</sup></td><td>2.9090(4)</td><td>2<sup>28</sup></td></tr><tr><td>NATIVE</td><td>0.5</td><td>2<sup>23.5</sup></td><td>4.0865(15)</td><td>2<sup>30</sup></td></tr><tr><td>NATIVE</td><td>0.5</td><td>2<sup>31.25</sup></td><td>5.6167(4)</td><td>2<sup>32</sup></td></tr><tr><td>JUMP</td><td>0.5</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>10</sup></td><td>16.62(12)</td><td>2<sup>28</sup></td></tr><tr><td>JUMP</td><td>0.5</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>15.625</sup></td><td>16.68(12)</td><td>2<sup>28</sup></td></tr><tr><td>JUMP</td><td>0.5</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>23.5</sup></td><td>16.75(13)</td><td>2<sup>31</sup></td></tr><tr><td>JUMP</td><td>0.5</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td><td>16.81(13)</td><td>2<sup>33.5</sup></td></tr><tr><td>NATIVE</td><td>0.7</td><td>2<sup>10</sup></td><td>2.6629(4)</td><td>2<sup>28</sup></td></tr><tr><td>NATIVE</td><td>0.7</td><td>2<sup>15.625</sup></td><td>3.8230(10)</td><td>2<sup>28</sup></td></tr><tr><td>NATIVE</td><td>0.7</td><td>2<sup>23.5</sup></td><td>6.1742(4)</td><td>2<sup>28</sup></td></tr><tr><td>NATIVE</td><td>0.7</td><td>2<sup>31.25</sup></td><td>9.578(11)</td><td>2<sup>33</sup></td></tr><tr><td>JUMP</td><td>0.7</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>10</sup></td><td>16.62(12)</td><td>2<sup>28</sup></td></tr><tr><td>JUMP</td><td>0.7</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>15.625</sup></td><td>16.67(12)</td><td>2<sup>28</sup></td></tr><tr><td>JUMP</td><td>0.7</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>23.5</sup></td><td>16.76(12)</td><td>2<sup>28</sup></td></tr><tr><td>JUMP</td><td>0.7</td><td>t<sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td><td>16.81(13)</td><td>2<sup>32</sup></td></tr></tbody></table>

of radius $R$. The results presented in this paper were obtained with $R = 5 a_0$ to make sure that the spheres will have a chance to fit within the glassy domains of the native runs (we have tried other values of $R$, finding qualitatively similar results, see Supplementary Note VI). The configurations from the two protocols are compared by computing a correlation coefficient $X$ that takes into account only the spins contained in the sphere. If $X$ is significantly smaller than unity we regard that particular sphere as *chaotic*, because typical configurations from the two protocols differ within the sphere. To be precise, we compute the probability distribution function $F(\tilde{X})$, namely the fraction of the spheres with a correlation coefficient $X < \tilde{X}$.

Our results shown in Fig. 4-bottom for the jump protocol $T_1 = 0.9 \rightarrow T_2 = 0.7$ remind us of previous studies [38]. The vast majority of the spheres have a very large correlation coefficient, and truly chaotic

spheres are found only in the tail of the distribution (probability 0.1% or smaller).

Interestingly enough, see the left panel in Fig. 2 and Fig. 4-top, the situation is radically different for the jump protocol $T_1 = 0.9 \rightarrow T_2 = 0.5$, where the spheres in percentile 10 of the distribution are as chaotic as the most chaotic spheres we could find for the jump $T_1 = 0.9 \rightarrow T_2 = 0.7$. In fact, to our knowledge, Fig. 4-top reports the strongest temperature-chaos signal ever observed in a simulation of glassy dynamics.

In order to convince ourselves that the extreme chaos is not an artifact of the disparity in domain sizes, we have tried a null experiment by simulating a model where no temperature chaos is expected, namely the link-diluted ferromagnetic Ising model (we have used the results in Ref. [46] to match as closely as possible in the diluted ferromagnet the conditions in our spin-glass simulations, with special care in matching the size of the domains, see Methods). As expected, see Fig. 4, the sphere distribution for the ferromagnet is concentrated at correlation coefficient $X \approx 1$. We conclude that the spin-glass results in Fig. 4-top are genuine evidence for dynamic temperature chaos.

It is also interesting that the distribution function in Fig. 4-top barely depends on $t_w$. This is another manifestation of the dynamic lock-down when the temperature jumps to the lower value.

The overall conclusions of this analysis are twofold. First, the requirement expressed by Eq. (2), which is based on CuMn experimental results [8], is sensible: strong temperature chaos is found only when $T_1 - T_2$ is as large as Eq. (2) demands. Second, only when temperature chaos is strong do our simulations find strong rejuvenation (recall Fig. 3-top).

### Where do we stand?

Our simulations depict a clear picture of the rejuvenation and memory effects. Provided that the temperature jump is large enough, see Eq. (2), the spin-glass state that was forming at temperature $T_1$ is completely alien at temperature $T_2$ (at least it looks like an alien when compared with the native state that grows directly at $T_2$, see Fig. 2). In fact, the response to the magnetic field (which is the quantity measured in experiments [2–5, 7–10, 15]) is not qualitatively different in the alien state and in the native state that grows from a fully disordered high-temperature state. The system just dismisses the relaxation it achieved at the higher temperature $T_1$.

Paradoxically enough, the alien state is locked at temperature $T_2$: the microscopic rearrangement at $T_2$, see Fig. 3-bottom, takes place on too small length scales to dissolve such foreign glassy domains. As a consequence, when the temperature is taken back to $T_1$, the glassy domains characteristic of $T_1$ are still there. This seems to be the physical origin of the memory effect. This reasoning is also consistent with recent experiments that find that the memory effect strongly depends on $t_w^↓$ (i.e., the time spent in the
<figure><img src="image_6.png" /><figcaption>FIG. 4. **Strong temperature chaos correlates with full rejuvenation.** The figure shows the fraction of the spheres with radius $R = 5a_0$ that have a correlation parameter $X$ smaller than $\tilde{X}$, $F(\tilde{X})$ (see Methods; the pairs of systems for which the correlation parameter $X$ is computed are listed in Table II). In the top panel, one of the systems in the pair used to compute $\tilde{X}$ is always taken from the native protocol at $T = 0.5$ ($T = 0.7$ for the bottom panel). In the cases reported in the top panel, the partner in the pair that undergoes the temperature-jump protocol experiences strong rejuvenation. Instead, recall Fig. 3-top, rejuvenation is only partial for the cases reported in the bottom panel. Interestingly enough, small correlation parameters appear with high probability in the top panel while they are very rare events in the bottom panel. We also show a comparison with the diluted Ising Model (our null experiment, see Methods), where temperature chaos is not expected. Indeed, in the absence of temperature chaos, the probability concentrates at $X \approx 1$. In all cases, error bars are one standard deviation.</figcaption></figure>

first stay at $T_1$) [47]. Indeed, if $t_w^\downarrow$ is too small, the memory effect almost disappears. Our interpretation of this experimental finding is that the glassy domains at $T_1$ need to grow large enough as to remain mostly unaltered at the lower temperature $T_2$.

Looking back, we understand as well why rejuvenation has been so difficult to find in simulations: the correlation lengths that could be reached prior to the Janus family of supercomputers were rather limited (we are referring here to the $\xi_{\text{micro}}$ length scale).

Therefore, Eq. (2) would demand an exceedingly large temperature jump $T_1 - T_2$ if one wants to have a large fraction of chaotic spheres of the relevant size.

<table><caption>TABLE II. Identifying parameters for each of the numerical simulations appearing in Figure 4. Spin glass (SG) protocols follow the notation in Figure 1. $T$ is the final temperature in the protocol. For the diluted Ising Model (DIM; see Methods for the DIM temperature-naming convention) we write explicitly $t_w^{(1)} + t_w^{(2)}$ for jump protocols to stress that the time before the jump (at $T=0.9$), $t_w^{(1)}$, differs from the time $t_w^{(2)}$ at the final temperature. We choose $t_w^{(1)}$ such that $\xi_{\text{micro}}$ coincides for both protocols in the pairs $(A_5, A'_5)$, $(A_7, A'_7)$ and $(A_9, A'_9)$.</caption><thead><tr><th>System</th><th><i>T</i></th><th>Type</th><th>Waiting time</th></tr></thead><tbody><tr><td>A<sub>9</sub></td><td>SG</td><td>0.9</td><td>NATIVE</td><td><i>t</i><sub>w</sub><sup>↓</sup> = 2<sup>31.25</sup></td></tr><tr><td>B<sub>9</sub></td><td>SG</td><td>0.9</td><td>JUMP-BACK</td><td><i>t</i><sub>w</sub><sup>↓</sup> + <i>t</i><sub>w</sub><sup>↑</sup> + 2<sup>10</sup></td></tr><tr><td>A<sub>5</sub></td><td>SG</td><td>0.5</td><td>NATIVE</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td></tr><tr><td>B<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>10</sup></td></tr><tr><td>C<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>15.625</sup></td></tr><tr><td>D<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>23.5</sup></td></tr><tr><td>E<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td></tr><tr><td>A<sub>7</sub></td><td>SG</td><td>0.7</td><td>NATIVE</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td></tr><tr><td>B<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>10</sup></td></tr><tr><td>C<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>15.625</sup></td></tr><tr><td>D<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>23.5</sup></td></tr><tr><td>E<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td></tr><tr><td>A'<sub>5</sub></td><td>DIM</td><td>0.5</td><td>NATIVE</td><td>76</td></tr><tr><td>B'<sub>5</sub></td><td>DIM</td><td>0.5</td><td>JUMP</td><td>430 + 69</td></tr><tr><td>A'<sub>7</sub></td><td>DIM</td><td>0.7</td><td>NATIVE</td><td>197</td></tr><tr><td>B'<sub>7</sub></td><td>DIM</td><td>0.7</td><td>JUMP</td><td>430 + 165</td></tr><tr><td>A'<sub>9</sub></td><td>DIM</td><td>0.9</td><td>NATIVE</td><td>430</td></tr></tbody></table>

An open question is whether or not the only experimentally accessible coherence length, namely $\xi_{\text{Zeeman}}(t_w)$, relates to some correlation function under all circumstances. Indeed, in the case of native protocols, $\xi_{\text{Zeeman}}(t_w)$ behaves analogously to $\xi_{\text{micro}}(t_w)$, which we know how to obtain from a correlation function. However, $\xi_{\text{micro}}(t_w)$ is not a valid proxy for $\xi_{\text{Zeeman}}(t_w)$ in temperature-jump protocols.

Finally, we should also stress that the analysis of the rejuvenation and memory effects requires considering no less than three different length scales, which can be quite different from each other. Of course, one of the three, the domain size $\xi_{\text{micro}}(t_w)$, acts as a cut-off for the other lengths. Yet we have seen that $\xi_{\text{micro}}(t_w)$ is not nearly enough to describe the variety of behaviors that an aging system may present. In fact, $\xi_{\text{micro}}(t_w)$ has stayed essentially constant for all the jump simulations that we have considered here! Therefore, a useful theory of aging dynamics cannot feature just a single length scale. In this sense, we think that our work poses a new and significant question for the different theories of aging dynamics.

<footnote>[1] C. L. E. Struik, *Physical Aging in Amorphous Polymers and Other Materials* (Elsevier, Amsterdam, 1980).</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 4. **Strong temperature chaos correlates with full rejuvenation.** The figure shows the fraction of the spheres with radius $R = 5a_0$ that have a correlation parameter $X$ smaller than $\tilde{X}$, $F(\tilde{X})$ (see Methods; the pairs of systems for which the correlation parameter $X$ is computed are listed in Table II). In the top panel, one of the systems in the pair used to compute $\tilde{X}$ is always taken from the native protocol at $T = 0.5$ ($T = 0.7$ for the bottom panel). In the cases reported in the top panel, the partner in the pair that undergoes the temperature-jump protocol experiences strong rejuvenation. Instead, recall Fig. 3-top, rejuvenation is only partial for the cases reported in the bottom panel. Interestingly enough, small correlation parameters appear with high probability in the top panel while they are very rare events in the bottom panel. We also show a comparison with the diluted Ising Model (our null experiment, see Methods), where temperature chaos is not expected. Indeed, in the absence of temperature chaos, the probability concentrates at $X \approx 1$. In all cases, error bars are one standard deviation.</figcaption></figure>

first stay at $T_1$) [47]. Indeed, if $t_w^\downarrow$ is too small, the memory effect almost disappears. Our interpretation of this experimental finding is that the glassy domains at $T_1$ need to grow large enough as to remain mostly unaltered at the lower temperature $T_2$.

Looking back, we understand as well why rejuvenation has been so difficult to find in simulations: the correlation lengths that could be reached prior to the Janus family of supercomputers were rather limited (we are referring here to the $\xi_{\text{micro}}$ length scale).

Therefore, Eq. (2) would demand an exceedingly large temperature jump $T_1 - T_2$ if one wants to have a large fraction of chaotic spheres of the relevant size.

<table><caption>TABLE II. Identifying parameters for each of the numerical simulations appearing in Figure 4. Spin glass (SG) protocols follow the notation in Figure 1. $T$ is the final temperature in the protocol. For the diluted Ising Model (DIM; see Methods for the DIM temperature-naming convention) we write explicitly $t_w^{(1)} + t_w^{(2)}$ for jump protocols to stress that the time before the jump (at $T=0.9$), $t_w^{(1)}$, differs from the time $t_w^{(2)}$ at the final temperature. We choose $t_w^{(1)}$ such that $\xi_{\text{micro}}$ coincides for both protocols in the pairs $(A_5, A'_5)$, $(A_7, A'_7)$ and $(A_9, A'_9)$.</caption><thead><tr><th>System</th><th><i>T</i></th><th>Type</th><th>Waiting time</th></tr></thead><tbody><tr><td>A<sub>9</sub></td><td>SG</td><td>0.9</td><td>NATIVE</td><td><i>t</i><sub>w</sub><sup>↓</sup> = 2<sup>31.25</sup></td></tr><tr><td>B<sub>9</sub></td><td>SG</td><td>0.9</td><td>JUMP-BACK</td><td><i>t</i><sub>w</sub><sup>↓</sup> + <i>t</i><sub>w</sub><sup>↑</sup> + 2<sup>10</sup></td></tr><tr><td>A<sub>5</sub></td><td>SG</td><td>0.5</td><td>NATIVE</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td></tr><tr><td>B<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>10</sup></td></tr><tr><td>C<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>15.625</sup></td></tr><tr><td>D<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>23.5</sup></td></tr><tr><td>E<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td></tr><tr><td>A<sub>7</sub></td><td>SG</td><td>0.7</td><td>NATIVE</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td></tr><tr><td>B<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>10</sup></td></tr><tr><td>C<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>15.625</sup></td></tr><tr><td>D<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>23.5</sup></td></tr><tr><td>E<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>↓</sup> + 2<sup>31.25</sup></td></tr><tr><td>A'<sub>5</sub></td><td>DIM</td><td>0.5</td><td>NATIVE</td><td>76</td></tr><tr><td>B'<sub>5</sub></td><td>DIM</td><td>0.5</td><td>JUMP</td><td>430 + 69</td></tr><tr><td>A'<sub>7</sub></td><td>DIM</td><td>0.7</td><td>NATIVE</td><td>197</td></tr><tr><td>B'<sub>7</sub></td><td>DIM</td><td>0.7</td><td>JUMP</td><td>430 + 165</td></tr><tr><td>A'<sub>9</sub></td><td>DIM</td><td>0.9</td><td>NATIVE</td><td>430</td></tr></tbody></table>

An open question is whether or not the only experimentally accessible coherence length, namely $\xi_{\text{Zeeman}}(t_w)$, relates to some correlation function under all circumstances. Indeed, in the case of native protocols, $\xi_{\text{Zeeman}}(t_w)$ behaves analogously to $\xi_{\text{micro}}(t_w)$, which we know how to obtain from a correlation function. However, $\xi_{\text{micro}}(t_w)$ is not a valid proxy for $\xi_{\text{Zeeman}}(t_w)$ in temperature-jump protocols.

Finally, we should also stress that the analysis of the rejuvenation and memory effects requires considering no less than three different length scales, which can be quite different from each other. Of course, one of the three, the domain size $\xi_{\text{micro}}(t_w)$, acts as a cut-off for the other lengths. Yet we have seen that $\xi_{\text{micro}}(t_w)$ is not nearly enough to describe the variety of behaviors that an aging system may present. In fact, $\xi_{\text{micro}}(t_w)$ has stayed essentially constant for all the jump simulations that we have considered here! Therefore, a useful theory of aging dynamics cannot feature just a single length scale. In this sense, we think that our work poses a new and significant question for the different theories of aging dynamics.

<footnote>[1] C. L. E. Struik, *Physical Aging in Amorphous Polymers and Other Materials* (Elsevier, Amsterdam, 1980).</footnote>
[2] K. Jonason, E. Vincent, J. Hammann, J. P. Bouchaud, and P. Nordblad, Phys. Rev. Lett. 81, 3243 (1998).

[3] L. Lundgren, P. Svedlindh, and O. Beckman, J. Magn. Magn. Mater. 31–34, 1349 (1983).

[4] T. Jonsson, K. Jonason, P. E. Jönsson, and P. Nordblad, Phys. Rev. B 59, 8770 (1999).

[5] J. Hammann, E. Vincent, V. Dupuis, M. Alba, M. Ocio, and J.-P. Bouchaud, J. Phys. Soc. Jpn., Suppl A. 206 (2000).

[6] J. A. Mydosh, Spin Glasses: an Experimental Introduction (Taylor and Francis, London, 1993).

[7] Q. Zhai, V. Martin-Mayor, D. L. Schlagel, G. G. Kenning, and R. L. Orbach, Phys. Rev. B 100, 094202 (2019).

[8] Q. Zhai, R. L. Orbach, and D. L. Schlagel, Phys. Rev. B 105, 014434 (2022).

[9] Q. Zhai, I. Paga, M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pemartin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz Sudupe, D. Navarro, R. L. Orbach, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, D. L. Schlagel, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes, Phys. Rev. Lett. 125, 237202 (2020).

[10] I. Paga, Q. Zhai, M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pemartin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz-Sudupe, D. Navarro, R. L. Orbach, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, D. L. Schlagel, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes, J. Stat. Mech. 2021, 033301 (2021).

[11] S. Albert, T. Bauer, M. Michl, G. Biroli, J.-P. Bouchaud, A. Loidl, P. Lunkenheimer, R. Tourbot, C. Wiertel-Gasquet, and F. Ladieu, Science 352, 1308 (2016), arXiv:1606.04079.

[12] M. Mézard, G. Parisi, and M. Virasoro, Spin-Glass Theory and Beyond (World Scientific, Singapore, 1987).

[13] M. Baity-Jesi, R. A. Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iniguez, A. Maiorano, F. Mantovani, E. Marinari, V. Martin-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, M. Pivanti, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Phys. Rev. B 88, 224416 (2013), arXiv:1310.2910.

[14] E. Vincent, J. Hammann, M. Ocio, J.-P. Bouchaud, and L. F. Cugliandolo, in *Complex Behavior of Glassy Systems*, Lecture Notes in Physics No. 492, edited by M. Rubí and C. Pérez-Vicente (Springer, 1997).

[15] C. Djurberg, K. Jonason, and P. Nordblad, Eur. Phys. J. B 10, 15 (1999).

[16] T. Komori, H. Yoshino, and H. Takayama, Journal of the Physical Society of Japan 69, 1192 (2000).

[17] M. Picco, F. Ricci-Tersenghi, and F. Ritort, Phys. Rev. B 63, 174412 (2001).

[18] L. Berthier and J.-P. Bouchaud, Phys. Rev. B 66, 054404 (2002).

[19] H. Takayama and K. Hukushima, Journal of the Physical Society of Japan 71, 3003 (2002).

[20] A. Maiorano, E. Marinari, and F. Ricci-Tersenghi, Phys. Rev. B 72, 104411 (2005).

[21] S. Jiménez, V. Martín-Mayor, and S. Pérez-Gaviro, Phys. Rev. B 72, 054417 (2005).

[22] S. F. Edwards and P. W. Anderson, Journal of Physics F: Metal Physics 5, 965 (1975).

[23] S. F. Edwards and P. W. Anderson, J. Phys. F 6, 1927 (1976).

[24] M. Baity-Jesi, R. A. Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iniguez, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, M. Pivanti, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Comp. Phys. Comm 185, 550 (2014), arXiv:1310.1032.

[25] E. Marinari, G. Parisi, J. Ruiz-Lorenzo, and F. Ritort, Phys. Rev. Lett. 76, 843 (1996).

[26] Y. G. Joh, R. Orbach, G. G. Wood, J. Hammann, and E. Vincent, Phys. Rev. Lett. 82, 438 (1999).

[27] F. Belletti, M. Cotallo, A. Cruz, L. A. Fernandez, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, A. M. Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, D. Sciretti, A. Tarancon, R. Tripiccione, J. L. Velasco, and D. Yllanes (Janus Collaboration), Phys. Rev. Lett. 101, 157201 (2008), arXiv:0804.1471.

[28] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Phys. Rev. Lett. 120, 267203 (2018).

[29] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Phys. Rev. Lett. 118, 157202 (2017).

[30] I. Paga, Q. Zhai, M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pemartin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz-Sudupe, D. Navarro, R. L. Orbach, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, D. L. Schlagel, B. Seoane, A. Tarancon, and D. Yllanes, "Magnetic-field symmetry breaking in spin glasses," In preparation (2022).

[31] L. F. Cugliandolo and J. Kurchan, Phys. Rev. B 60, 922 (1999).

[32] L. Berthier and J.-P. Bouchaud, Phys. Rev. Lett 90, 059701 (2003).

[33] S. R. McKay, A. N. Berker, and S. Kirkpatrick, Phys. Rev. Lett. 48, 767 (1982).

[34] A. J. Bray and M. A. Moore, Phys. Rev. Lett. 58, 57 (1987).

[35] I. Kondor, J. Phys. A 22, L163 (1989).

[36] T. Rizzo and A. Crisanti, Phys. Rev. Lett. 90, 137201 (2003).
[2] K. Jonason, E. Vincent, J. Hammann, J. P. Bouchaud, and P. Nordblad, Phys. Rev. Lett. 81, 3243 (1998).

[3] L. Lundgren, P. Svedlindh, and O. Beckman, J. Magn. Magn. Mater. 31–34, 1349 (1983).

[4] T. Jonsson, K. Jonason, P. E. Jönsson, and P. Nordblad, Phys. Rev. B 59, 8770 (1999).

[5] J. Hammann, E. Vincent, V. Dupuis, M. Alba, M. Ocio, and J.-P. Bouchaud, J. Phys. Soc. Jpn., Suppl A. 206 (2000).

[6] J. A. Mydosh, Spin Glasses: an Experimental Introduction (Taylor and Francis, London, 1993).

[7] Q. Zhai, V. Martin-Mayor, D. L. Schlagel, G. G. Kenning, and R. L. Orbach, Phys. Rev. B 100, 094202 (2019).

[8] Q. Zhai, R. L. Orbach, and D. L. Schlagel, Phys. Rev. B 105, 014434 (2022).

[9] Q. Zhai, I. Paga, M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pemartin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz Sudupe, D. Navarro, R. L. Orbach, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, D. L. Schlagel, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes, Phys. Rev. Lett. 125, 237202 (2020).

[10] I. Paga, Q. Zhai, M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pemartin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz-Sudupe, D. Navarro, R. L. Orbach, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, D. L. Schlagel, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes, J. Stat. Mech. 2021, 033301 (2021).

[11] S. Albert, T. Bauer, M. Michl, G. Biroli, J.-P. Bouchaud, A. Loidl, P. Lunkenheimer, R. Tourbot, C. Wiertel-Gasquet, and F. Ladieu, Science 352, 1308 (2016), arXiv:1606.04079.

[12] M. Mézard, G. Parisi, and M. Virasoro, Spin-Glass Theory and Beyond (World Scientific, Singapore, 1987).

[13] M. Baity-Jesi, R. A. Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iniguez, A. Maiorano, F. Mantovani, E. Marinari, V. Martin-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, M. Pivanti, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Phys. Rev. B 88, 224416 (2013), arXiv:1310.2910.

[14] E. Vincent, J. Hammann, M. Ocio, J.-P. Bouchaud, and L. F. Cugliandolo, in *Complex Behavior of Glassy Systems*, Lecture Notes in Physics No. 492, edited by M. Rubí and C. Pérez-Vicente (Springer, 1997).

[15] C. Djurberg, K. Jonason, and P. Nordblad, Eur. Phys. J. B 10, 15 (1999).

[16] T. Komori, H. Yoshino, and H. Takayama, Journal of the Physical Society of Japan 69, 1192 (2000).

[17] M. Picco, F. Ricci-Tersenghi, and F. Ritort, Phys. Rev. B 63, 174412 (2001).

[18] L. Berthier and J.-P. Bouchaud, Phys. Rev. B 66, 054404 (2002).

[19] H. Takayama and K. Hukushima, Journal of the Physical Society of Japan 71, 3003 (2002).

[20] A. Maiorano, E. Marinari, and F. Ricci-Tersenghi, Phys. Rev. B 72, 104411 (2005).

[21] S. Jiménez, V. Martín-Mayor, and S. Pérez-Gaviro, Phys. Rev. B 72, 054417 (2005).

[22] S. F. Edwards and P. W. Anderson, Journal of Physics F: Metal Physics 5, 965 (1975).

[23] S. F. Edwards and P. W. Anderson, J. Phys. F 6, 1927 (1976).

[24] M. Baity-Jesi, R. A. Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iniguez, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, M. Pivanti, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Comp. Phys. Comm 185, 550 (2014), arXiv:1310.1032.

[25] E. Marinari, G. Parisi, J. Ruiz-Lorenzo, and F. Ritort, Phys. Rev. Lett. 76, 843 (1996).

[26] Y. G. Joh, R. Orbach, G. G. Wood, J. Hammann, and E. Vincent, Phys. Rev. Lett. 82, 438 (1999).

[27] F. Belletti, M. Cotallo, A. Cruz, L. A. Fernandez, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, A. M. Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, D. Sciretti, A. Tarancon, R. Tripiccione, J. L. Velasco, and D. Yllanes (Janus Collaboration), Phys. Rev. Lett. 101, 157201 (2008), arXiv:0804.1471.

[28] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Phys. Rev. Lett. 120, 267203 (2018).

[29] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Phys. Rev. Lett. 118, 157202 (2017).

[30] I. Paga, Q. Zhai, M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pemartin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz-Sudupe, D. Navarro, R. L. Orbach, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, D. L. Schlagel, B. Seoane, A. Tarancon, and D. Yllanes, "Magnetic-field symmetry breaking in spin glasses," In preparation (2022).

[31] L. F. Cugliandolo and J. Kurchan, Phys. Rev. B 60, 922 (1999).

[32] L. Berthier and J.-P. Bouchaud, Phys. Rev. Lett 90, 059701 (2003).

[33] S. R. McKay, A. N. Berker, and S. Kirkpatrick, Phys. Rev. Lett. 48, 767 (1982).

[34] A. J. Bray and M. A. Moore, Phys. Rev. Lett. 58, 57 (1987).

[35] I. Kondor, J. Phys. A 22, L163 (1989).

[36] T. Rizzo and A. Crisanti, Phys. Rev. Lett. 90, 137201 (2003).
[37] G. Parisi and T. Rizzo, Journal of Physics A: Mathematical and Theoretical 43, 235003 (2010).

[38] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pe Martin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz Sudupe, D. Navarro, I. Paga, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Taranton, R. Tripiccione, and D. Yllanes, Commun. Phys. 4, 74 (2021).

[39] We are indebted to Prof. Orbach for this observation.

[40] Our simulations are also described in the Methods section, see also Table I for crucial simulation details and Ref. [48] for useful computational tricks.

[41] F. Belletti, A. Cruz, L. A. Fernandez, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, J. Monforte, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, D. Sciretti, A. Taranton, R. Tripiccione, and D. Yllanes (Janus Collaboration), J. Stat. Phys. 135, 1121 (2009).

[42] H. E. Castillo, C. Chamon, L. F. Cugliandolo, and M. P. Kennett, Phys. Rev. Lett. 88, 237201 (2002).

[43] L. C. Jaubert, C. Chamon, L. F. Cugliandolo, and M. Picco, J. Stat. Mech. 2007, P05001 (2007).

[44] In fact, d log $t_w$/d log $\xi_{micro}$ is approximately constant when $t_w$ varies in logarithmic scale.

[45] T'd log $t_w$/d log $\xi_{micro}$ is roughly constant when different temperatures are compared.

[46] P.-E. Berche, C. Chatelain, B. Berche, and W. Janke, Eur. Phys. J. B 38, 463 (2004).

[47] Jennifer Freedberg, private communication (2022).

[48] I. Paga, From glassy bulk systems to spin-glass films: simulations meet experiments, Ph.D. thesis (2021).

## Appendix A: Methods

The layout of this note is as follows. In Sect. A1 we describe our simulations. In Sect. A2 we define some quantities characteristic of the ZFC protocol. In fact, the magnetic field plays a crucial role in the determination of the Zeeman length scale, as we explain in Sect. A3. The other two spin-glass coherence lengths, $\xi_{micro}$ and $\zeta$, are computed as explained in Sect. A4. Finally, in Sect. A5 we explain our computation of the chaotic correlation parameter.

### 1. The models simulated

We performed massive simulations on the Janus II supercomputer [24] to study the three-dimensional Edwards-Anderson (EA) model on a cubic lattice with periodic boundary conditions and size $L = 160$ (in units of the lattice constant $a_0$). The main parameters describing our simulations are provided in Tab. I.

The $N = L^3$ Ising spins, $s_x = \pm 1$, interact with their lattice nearest neighbors in presence of a magnetic field ($H$) through the Hamiltonian:

$$ \mathcal{H} = - \sum_{\langle x, y \rangle} J_{xy} s_x s_y - H \sum_x s_x, \quad (A1) $$

where the couplings are independent, identically distributed random variables: $J_{xy} = \pm 1$, with 50% probability. The couplings are chosen at simulation start, and remain fixed (quenched disorder). A particular choice of the couplings is termed a sample. In the absence of an external magnetic field $H = 0$, this model undergoes a spin-glass transition at the critical temperature $T_g = 1.102(3)$ [13].

The off-equilibrium dynamics was simulated with a Metropolis algorithm. The numerical time unit is the lattice sweep, which roughly corresponds to 1 ps of physical time.

In this work we have simulated $N_S = 4$ samples using a lattice size of $L = 160$ $a_0$. For each of these samples and for each protocol (Tab. I in the main text) we have simulated $N_R = 512$ replicas (i.e., independent simulations carried out for a given sample, following an identical protocol). We use replicas to account for the thermal noise controlling the simulation (each replica is controlled by an independent realization of the thermal noise). The average over the thermal noise will be represented as $\langle \cdots \rangle$. Only afterwards, we shall perform the average over samples, which will be indicated as $\langle \cdots \rangle$.

Some times, however, (most notably for the analysis in Sect. A3) final quantities are computed for a single sample (this is, of course, the approach followed in the laboratory). In these cases, the different samples allow us to assess to which extent our results depend on the disorder realization, see Supplementary Note I.

Besides, as a null experiment for temperature chaos, we have studied the link-diluted Ising model (DIM), also on cubic lattices of size $L = 160$ $a_0$ with periodic boundary conditions and using Metropolis dynamics. Specifically, we used the Hamiltonian in Eq. (A1) but with couplings $J_{xy} = 1$ (with 70% probability) or $J_{xy} = 0$ (with 30% probability) and magnetic field $H = 0$. Since all couplings are positive or zero, this is a ferromagnetic system without frustration, for which no temperature chaos is expected. The critical temperature for the DIM is $T_c = 3.0609(5)$ [46] (actually, this is twice the value reported in [46] due to our use of an Ising, rather than Potts, formulation). In fact, with some abuse of language, in the main text we refer to DIM temperatures as $T = 0.9$, $T = 0.7$ or $0.5$ rather than to their actual values $T = 0.9 T_c / T_g$, $T = 0.7 T_c / T_g$ or $0.5 T_c / T_g$, where $T_g$ is the critical temperature for the EA model. We follow the very same procedure, which is explained in Sect. A4a, to compute the coherence length $\xi_{micro}$ for both the spin glass and the DIM. We have chosen times for the DIM such that $\xi_{micro}$ coincides with the corresponding spin-glass value, namely $\xi_{micro} = 5.84$ (protocol A'$_5$ in Table II in the main text), $\xi_{micro} = 10.11$ (protocol A'$_7$) and $\xi_{micro} = 16.63$ (protocol A'$_9$). Of course, the necessary times are extremely shorter for the DIM than for the spin glass. Given that DIM simulations were comparatively inexpensive, we simulated 16 samples (each with 512 replicas) for this model.
[37] G. Parisi and T. Rizzo, Journal of Physics A: Mathematical and Theoretical 43, 235003 (2010).

[38] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pe Martin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz Sudupe, D. Navarro, I. Paga, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Taranton, R. Tripiccione, and D. Yllanes, Commun. Phys. 4, 74 (2021).

[39] We are indebted to Prof. Orbach for this observation.

[40] Our simulations are also described in the Methods section, see also Table I for crucial simulation details and Ref. [48] for useful computational tricks.

[41] F. Belletti, A. Cruz, L. A. Fernandez, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, J. Monforte, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, D. Sciretti, A. Taranton, R. Tripiccione, and D. Yllanes (Janus Collaboration), J. Stat. Phys. 135, 1121 (2009).

[42] H. E. Castillo, C. Chamon, L. F. Cugliandolo, and M. P. Kennett, Phys. Rev. Lett. 88, 237201 (2002).

[43] L. C. Jaubert, C. Chamon, L. F. Cugliandolo, and M. Picco, J. Stat. Mech. 2007, P05001 (2007).

[44] In fact, d log $t_w$/d log $\xi_{micro}$ is approximately constant when $t_w$ varies in logarithmic scale.

[45] T'd log $t_w$/d log $\xi_{micro}$ is roughly constant when different temperatures are compared.

[46] P.-E. Berche, C. Chatelain, B. Berche, and W. Janke, Eur. Phys. J. B 38, 463 (2004).

[47] Jennifer Freedberg, private communication (2022).

[48] I. Paga, From glassy bulk systems to spin-glass films: simulations meet experiments, Ph.D. thesis (2021).

## Appendix A: Methods

The layout of this note is as follows. In Sect. A1 we describe our simulations. In Sect. A2 we define some quantities characteristic of the ZFC protocol. In fact, the magnetic field plays a crucial role in the determination of the Zeeman length scale, as we explain in Sect. A3. The other two spin-glass coherence lengths, $\xi_{micro}$ and $\zeta$, are computed as explained in Sect. A4. Finally, in Sect. A5 we explain our computation of the chaotic correlation parameter.

### 1. The models simulated

We performed massive simulations on the Janus II supercomputer [24] to study the three-dimensional Edwards-Anderson (EA) model on a cubic lattice with periodic boundary conditions and size $L = 160$ (in units of the lattice constant $a_0$). The main parameters describing our simulations are provided in Tab. I.

The $N = L^3$ Ising spins, $s_x = \pm 1$, interact with their lattice nearest neighbors in presence of a magnetic field ($H$) through the Hamiltonian:

$$ \mathcal{H} = - \sum_{\langle x, y \rangle} J_{xy} s_x s_y - H \sum_x s_x, \quad (A1) $$

where the couplings are independent, identically distributed random variables: $J_{xy} = \pm 1$, with 50% probability. The couplings are chosen at simulation start, and remain fixed (quenched disorder). A particular choice of the couplings is termed a sample. In the absence of an external magnetic field $H = 0$, this model undergoes a spin-glass transition at the critical temperature $T_g = 1.102(3)$ [13].

The off-equilibrium dynamics was simulated with a Metropolis algorithm. The numerical time unit is the lattice sweep, which roughly corresponds to 1 ps of physical time.

In this work we have simulated $N_S = 4$ samples using a lattice size of $L = 160$ $a_0$. For each of these samples and for each protocol (Tab. I in the main text) we have simulated $N_R = 512$ replicas (i.e., independent simulations carried out for a given sample, following an identical protocol). We use replicas to account for the thermal noise controlling the simulation (each replica is controlled by an independent realization of the thermal noise). The average over the thermal noise will be represented as $\langle \cdots \rangle$. Only afterwards, we shall perform the average over samples, which will be indicated as $\langle \cdots \rangle$.

Some times, however, (most notably for the analysis in Sect. A3) final quantities are computed for a single sample (this is, of course, the approach followed in the laboratory). In these cases, the different samples allow us to assess to which extent our results depend on the disorder realization, see Supplementary Note I.

Besides, as a null experiment for temperature chaos, we have studied the link-diluted Ising model (DIM), also on cubic lattices of size $L = 160$ $a_0$ with periodic boundary conditions and using Metropolis dynamics. Specifically, we used the Hamiltonian in Eq. (A1) but with couplings $J_{xy} = 1$ (with 70% probability) or $J_{xy} = 0$ (with 30% probability) and magnetic field $H = 0$. Since all couplings are positive or zero, this is a ferromagnetic system without frustration, for which no temperature chaos is expected. The critical temperature for the DIM is $T_c = 3.0609(5)$ [46] (actually, this is twice the value reported in [46] due to our use of an Ising, rather than Potts, formulation). In fact, with some abuse of language, in the main text we refer to DIM temperatures as $T = 0.9$, $T = 0.7$ or $0.5$ rather than to their actual values $T = 0.9 T_c/T_g$, $T = 0.7 T_c/T_g$ or $0.5 T_c/T_g$, where $T_g$ is the critical temperature for the EA model. We follow the very same procedure, which is explained in Sect. A4a, to compute the coherence length $\xi_{micro}$ for both the spin glass and the DIM. We have chosen times for the DIM such that $\xi_{micro}$ coincides with the corresponding spin-glass value, namely $\xi_{micro} = 5.84$ (protocol A'5 in Table II in the main text), $\xi_{micro} = 10.11$ (protocol A'7) and $\xi_{micro} = 16.63$ (protocol A'9). Of course, the necessary times are extremely shorter for the DIM than for the spin glass. Given that DIM simulations were comparatively inexpensive, we simulated 16 samples (each with 512 replicas) for this model.
## 2. Some zero-field-cooled observables

As explained in the main text, our simulations are designed to mimic the experimental protocol named zero-field cooling (ZFC). In ZFC protocols, a sample initially in equilibrium at some very high temperature is cooled below $T_g$, always being kept at zero magnetic field. In the native protocols, the system is abruptly taken to the measuring temperature, where it is let to relax for a time $t_w$. The cooling process (always without a field) is more complex for our jump protocols, as depicted in Fig. 1 in the main text.

For both protocols, native or jump, we let the system relax for a time $t_w$ at the final, measuring temperature. Then, the external magnetic field, $H$, is switched on and we record the magnetic density

$$M_{ZFC}(t, t_w; H) = \frac{1}{N} \sum_{x} \langle s_x(t + t_w; H) \rangle, \quad (A2)$$

which grows with $t$ from its initial value $M=0$ at $t=0$. We also record the two-time autocorrelation function,

$$C_{ZFC}(t, t_w; H) = \frac{1}{N} \sum_{x} \langle s_x(t_w; 0) s_x(t + t_w; H) \rangle. \quad (A3)$$

Note that $C_{ZFC}$ is a monotonically decreasing function of time and $C_{ZFC} = 1$ at $t=0$.

## 3. Measurement of the Zeeman length through the scaling law of the effective times

The method introduced in Ref. [26] to measure the spin-glass coherence length experimentally has recently been refined. Indeed, the scaling law introduced in [9, 10] is a milestone for describing the magnetic response of a spin glass in both “lab experiments” and “numerical experiments”. We shall name $\xi_{\text{Zeeman}}$ the length scale extracted using these methods.

In experiments on a single-crystal CuMn sample, the main quantity evaluated is the relaxation function $S_{ZFC}(t, t_w; H)$, which exhibits a local maximum at time $t_H^{\text{eff}} \approx t_w$. Hence, one focuses on the $H$ dependence of $t_H^{\text{eff}}$. On the numerical side, we carry out massive numerical experiments spanning from picoseconds to tenths of a second on Janus II, from which we can also extract the $t_H^{\text{eff}}$. The numerical method proceeds as follows (see [9, 10] for a full discussion). One first changes variable by considering $S_{ZFC}$ as a function of $C(t, t_w; H)$, recall Eq. (A3), rather than time. The peak is found at some $C_{\text{peak}}(t_w)$. Finally, $t_H^{\text{eff}}$ is found by solving the equation $C(t_H^{\text{eff}}, t_w; H) = C_{\text{peak}}(t_w)$. A crucial advantage is that this equation can also be solved directly at $H=0$.

The numerical $S_{ZFC}(t, t_w; H)$, however, shows two peaks: a $t_w$-independent peak at very short times, and a second, physically interesting peak at $t \sim t_w$. Unfortunately, in fixed-temperature simulations (i.e., native protocols) with very short $t_w$, the two peaks cannot

be resolved (see, for instance, bottom-left of Fig. 1). We have not attempted to extract $\xi_{\text{Zeeman}}$ in native runs where the two peaks cannot be resolved. However, for the shortest jump protocol with $T_2 = 0.5$, namely $t_w = 2^{10}, 2^{15.625}$, we could borrow $C_{\text{peak}}$ from the jump with the largest $t_w$ (unfortunately, the same trick did not work for native runs, because important consistency checks [30] were not passed in this case).

From a phenomenological point of view, the effective time $t_H^{\text{eff}}$ can be associated with the height of the largest free-energy barrier, $\Delta_{\text{max}}$, through the usual Arrhenius law [26]

$$\Delta_{\text{max}} = k_B T (\log t_H^{\text{eff}} - \log \tau_0), \quad (A4)$$

where $\tau_0$ is a characteristic exchange time, $\tau_0 \sim \hbar/k_B T_g$. In an external magnetic field, the free-energy barriers are lowered by the Zeeman energy $E_Z$ [26]. For small magnetic field, $E_Z$ behaves as:

$$E_Z = \xi_{\text{Zeeman}}^{D-\theta/2} \chi_{\text{FC}} H^2, \quad (A5)$$

which defines $\xi_{\text{Zeeman}}$. $\chi_{\text{FC}}$ is the field-cooled magnetic susceptibility per spin, $\xi_{\text{Zeeman}}^{D-\theta/2}$ is the number of correlated spins, $D=3$ is the spatial dimension and $\theta$ is the replicon exponent [29].

We slightly depart from the previous approach by exploiting a scaling theory. We use the effective time $t_H^{\text{eff}}$ to reflect the total free-energy change at magnetic fields $H$ and $H=0^+$ [9, 10]:

$$\log \left[ \frac{t_H^{\text{eff}}}{t_H^{\text{eff}} \to 0^+} \right] = \frac{\hat{S}}{2T} \xi_{\text{micro}}^{D-\theta/2} H^2 + \xi_{\text{micro}}^{-\theta/2} \mathcal{G}(T, \xi_{\text{micro}}^{D-\theta/2} H^2), \quad (A6)$$

where $\hat{S}$ is a constant coming from the fluctuation-dissipation relations and $\mathcal{G}(x)$ is a scaling function behaving as $\mathcal{G}(x) \sim x^2$ for small $x = \xi_{\text{micro}}^{D-\theta/2} H^2$. For small-enough magnetic fields [$H \le 0.017$], we can neglect the $\mathcal{O}(H^4)$ terms in Eq. (A6):

$$\log \left[ \frac{t_H^{\text{eff}}}{t_H^{\text{eff}} \to 0^+} \right] = c_2(t_w; T) H^2, \quad (A7)$$

where we have included all the constants in the $c_2(t_w; T)$ coefficient.

Thus, fitting our data according to Eq. (A7), we can define the Zeeman coherence length $\xi_{\text{Zeeman}}$ as

$$\xi_{\text{Zeeman}}^{\text{jump}}(t_w, T_1 \to T_m) = \left[ \frac{c_2(t_w, T_1 \to T_m)}{c_2(t_w^*, T_m)} \right]^{1/(D-\theta/2)} \xi_{\text{micro}}(t_w^*; T_m), \quad (A8)$$

$$\xi_{\text{Zeeman}}^{\text{native}}(t_w, T_m) = \left[ \frac{c_2(t_w, T_m)}{c_2(t_w^*, T_m)} \right]^{1/(D-\theta/2)} \xi_{\text{micro}}(t_w^*; T_m). \quad (A9)$$

where $\xi_{\text{micro}}(t_w^*; T_m)$ plays the role of a reference length [the reference length allows us to avoid the precise determination of constants in Eq. (A7)]. The
## 2. Some zero-field-cooled observables

As explained in the main text, our simulations are designed to mimic the experimental protocol named zero-field cooling (ZFC). In ZFC protocols, a sample initially in equilibrium at some very high temperature is cooled below $T_g$, always being kept at zero magnetic field. In the native protocols, the system is abruptly taken to the measuring temperature, where it is let to relax for a time $t_w$. The cooling process (always without a field) is more complex for our jump protocols, as depicted in Fig. 1 in the main text.

For both protocols, native or jump, we let the system relax for a time $t_w$ at the final, measuring temperature. Then, the external magnetic field, $H$, is switched on and we record the magnetic density

$$M_{ZFC}(t, t_w; H) = \frac{1}{N} \sum_{x} \langle s_x(t + t_w; H) \rangle, \quad (A2)$$

which grows with $t$ from its initial value $M = 0$ at $t = 0$. We also record the two-time autocorrelation function,

$$C_{ZFC}(t, t_w; H) = \frac{1}{N} \sum_{x} \langle s_x(t_w; 0) s_x(t + t_w; H) \rangle. \quad (A3)$$

Note that $C_{ZFC}$ is a monotonically decreasing function of time and $C_{ZFC} = 1$ at $t = 0$.

## 3. Measurement of the Zeeman length through the scaling law of the effective times

The method introduced in Ref. [26] to measure the spin-glass coherence length experimentally has recently been refined. Indeed, the scaling law introduced in [9, 10] is a milestone for describing the magnetic response of a spin glass in both “lab experiments” and “numerical experiments”. We shall name $\xi_{\text{Zeeman}}$ the length scale extracted using these methods.

In experiments on a single-crystal CuMn sample, the main quantity evaluated is the relaxation function $S_{ZFC}(t, t_w; H)$, which exhibits a local maximum at time $t_H^{\text{eff}} \approx t_w$. Hence, one focuses on the $H$ dependence of $t_H^{\text{eff}}$. On the numerical side, we carry out massive numerical experiments spanning from picoseconds to tenths of a second on Janus II, from which we can also extract the $t_H^{\text{eff}}$. The numerical method proceeds as follows (see [9, 10] for a full discussion). One first changes variable by considering $S_{ZFC}$ as a function of $C(t, t_w; H)$, recall Eq. (A3), rather than time. The peak is found at some $C_{\text{peak}}(t_w)$. Finally, $t_H^{\text{eff}}$ is found by solving the equation $C(t_H^{\text{eff}}, t_w; H) = C_{\text{peak}}(t_w)$. A crucial advantage is that this equation can also be solved directly at $H = 0$.

The numerical $S_{ZFC}(t, t_w; H)$, however, shows two peaks: a $t_w$-independent peak at very short times, and a second, physically interesting peak at $t \sim t_w$. Unfortunately, in fixed-temperature simulations (i.e., native protocols) with very short $t_w$, the two peaks cannot

be resolved (see, for instance, bottom-left of Fig. 1). We have not attempted to extract $\xi_{\text{Zeeman}}$ in native runs where the two peaks cannot be resolved. However, for the shortest jump protocol with $T_2 = 0.5$, namely $t_w = 2^{10}, 2^{15.625}$, we could borrow $C_{\text{peak}}$ from the jump with the largest $t_w$ (unfortunately, the same trick did not work for native runs, because important consistency checks [30] were not passed in this case).

From a phenomenological point of view, the effective time $t_H^{\text{eff}}$ can be associated with the height of the largest free-energy barrier, $\Delta_{\text{max}}$, through the usual Arrhenius law [26]

$$\Delta_{\text{max}} = k_B T (\log t_H^{\text{eff}} - \log \tau_0), \quad (A4)$$

where $\tau_0$ is a characteristic exchange time, $\tau_0 \sim \hbar/k_B T_g$. In an external magnetic field, the free-energy barriers are lowered by the Zeeman energy $E_Z$ [26]. For small magnetic field, $E_Z$ behaves as:

$$E_Z = \xi_{\text{Zeeman}}^{D-\theta/2} \chi_{\text{FC}} H^2, \quad (A5)$$

which defines $\xi_{\text{Zeeman}}$. $\chi_{\text{FC}}$ is the field-cooled magnetic susceptibility per spin, $\xi_{\text{Zeeman}}^{D-\theta/2}$ is the number of correlated spins, $D = 3$ is the spatial dimension and $\theta$ is the replicon exponent [29].

We slightly depart from the previous approach by exploiting a scaling theory. We use the effective time $t_H^{\text{eff}}$ to reflect the total free-energy change at magnetic fields $H$ and $H = 0^+$ [9, 10]:

$$\log \left[ \frac{t_H^{\text{eff}}}{t_H^{\text{eff}} \to 0^+} \right] = \frac{\hat{S}}{2T} \xi_{\text{micro}}^{D-\theta/2} H^2 + \xi_{\text{micro}}^{-\theta/2} \mathcal{G}(T, \xi_{\text{micro}}^{D-\theta/2} H^2), \quad (A6)$$

where $\hat{S}$ is a constant coming from the fluctuation-dissipation relations and $\mathcal{G}(x)$ is a scaling function behaving as $\mathcal{G}(x) \sim x^2$ for small $x = \xi_{\text{micro}}^{D-\theta/2} H^2$. For small-enough magnetic fields [$H \le 0.017$], we can neglect the $\mathcal{O}(H^4)$ terms in Eq. (A6):

$$\log \left[ \frac{t_H^{\text{eff}}}{t_H^{\text{eff}} \to 0^+} \right] = c_2(t_w; T) H^2, \quad (A7)$$

where we have included all the constants in the $c_2(t_w; T)$ coefficient.

Thus, fitting our data according to Eq. (A7), we can define the Zeeman coherence length $\xi_{\text{Zeeman}}$ as

$$\xi_{\text{Zeeman}}^{\text{jump}}(t_w, T_1 \to T_m) = \left[ \frac{c_2(t_w, T_1 \to T_m)}{c_2(t_w^*, T_m)} \right]^{1/(D-\theta/2)} \xi_{\text{micro}}(t_w^*; T_m), \quad (A8)$$

$$\xi_{\text{Zeeman}}^{\text{native}}(t_w, T_m) = \left[ \frac{c_2(t_w, T_m)}{c_2(t_w^*, T_m)} \right]^{1/(D-\theta/2)} \xi_{\text{micro}}(t_w^*; T_m). \quad (A9)$$

where $\xi_{\text{micro}}(t_w^*; T_m)$ plays the role of a reference length [the reference length allows us to avoid the precise determination of constants in Eq. (A7)]. The
reference time $t_w^*$ is the longest available waiting time for our native runs at the measuring temperature $T_m$. For the sake of clarity, we omit in Eqs. (A8) and (A9) the explicit dependence of $\theta$ on $\xi_{\text{micro}}$ (which is dealt with as explained in Ref. [7]).

### 4. Numerical coherence lengths $\xi_{\text{micro}}$ and $\zeta$

In this paragraph, we shall consider two more length scales. One of them, $\xi_{\text{micro}}$, is computed from the correlation function for the spin-glass order parameter (hence, $\xi_{\text{micro}}$ tells us about the size of the glassy domains). The second length scale, $\zeta(t_1, t_2)$, tells us about how the system reorganizes itself when going from the earlier time $t_1$ to the later time $t_2$.

#### a. The computation of $\xi_{\text{micro}}$

For the reader's convenience, let us recall the definition of the spatial autocorrelation function that we use for computing $\xi_{\text{micro}}(t_w)$ [41]

$$C_4(\mathbf{r}, t'; T) = \overline{\langle q^{(a,b)}(\mathbf{x}, t')q^{(a,b)}(\mathbf{x} + \mathbf{r}, t') \rangle_T}, \quad (\text{A10})$$

$$q^{(a,b)}(\mathbf{x}, t') \equiv \sigma^{(a)}(\mathbf{x}, t')\sigma^{(b)}(\mathbf{x}, t'), \quad (\text{A11})$$

where $t' = t_w + t$, the indices $(a, b)$ label different real replicas and $\langle \cdots \rangle_T$ stands for the average over the thermal noise at temperature $T$.

The calculation of the correlation function is computationally costly since we have $N_R(N_R - 1)/2$ possible choices of the pair of replicas. Fortunately, it can be accelerated using the specific multispin coding methods explained in Ref. [48].

Once we have $C_4(\mathbf{r}, t'; T)$, we compute the integrals [27, 28, 41]:

$$I_k(t'; T) = \int_0^\infty d^3\mathbf{r} r^k C_4(\mathbf{r} = (\mathbf{r}, 0, 0), t'; T). \quad (\text{A12})$$

A coherence length can be computed as

$$\xi_{k,k+1}(t', T) = \frac{I_{k+1}(t', T)}{I_k(t', T)}. \quad (\text{A13})$$

We define $\xi_{\text{micro}}(t, t_w; H) = \xi_{12}(t, t_w; H)$.

#### b. The $\zeta$ length scale

This length scale was studied in details in Ref. [41] by refining earlier suggestions [42, 43].

Let us consider the thermal trajectory followed by a given replica at the two times $t_1 < t_2$. Our basic quantity will be the local correlation

$$c_x(t_1, t_2) = s_x(t_2)s_x(t_1). \quad (\text{A14})$$

Note that $c_x(t_1, t_2) = -1$ if the spin at site $\mathbf{x}$ has been flipped when going from time $t_1$ to time $t_2$ [otherwise, $c_x(t_1, t_2) = 1$]. Then, the two-time, two-site correlation function is

$$C_{2+2}(\mathbf{r}, t_1, t_2) = \frac{1}{N} \sum_{\mathbf{x}} \overline{[\langle c_{\mathbf{x}}(t_1, t_2) c_{\mathbf{x}+\mathbf{r}}(t_1, t_2) \rangle - C^2(t_1, t_2)]}, \quad (\text{A15})$$

where

$$C(t_1, t_2) = \frac{1}{N} \sum_{\mathbf{x}} \langle c_{\mathbf{x}}(t_1, t_2) \rangle. \quad (\text{A16})$$

The ideal $\zeta(t_1, t_2)$ is defined from the long-distance decay of $C_{2+2}(\mathbf{r}, t, t_w)$:

$$C_{2+2}(\mathbf{r}, t_1, t_2) \sim \frac{1}{rb} g(r/\zeta(t_1, t_2)), \quad (\text{A17})$$

where $g$ is an unknown scaling function. We bypass our lack of knowledge of $g$ exactly as we solved this problem for $\xi_{\text{micro}}$: by using integral estimators, recall Eq. (A13). Note that, by construction, $\zeta(t_1, t_2)$ tends to zero when $t_2$ approaches $t_1$. Conversely, we expect $\zeta(t_1, t_2)$ to grow with the later time $t_2$.

As for the interpretation of the length scale $\zeta$, an analogy with the theory of liquids is of help. We name a *defect* a site where $c_x(t_1, t_2) = -1$. Let $n(t_1, t_2)$ be the density of defects $[C(t_1, t_2) = 1 - 2n(t_1, t_2)]$ and let $g(\mathbf{r})$ be the pair-correlation function for defects: The conditional probability for having a defect at site $\mathbf{x}+\mathbf{r}$, given that a defect is present at site $\mathbf{x}$, is $n(t_1, t_2)g(\mathbf{r})$ (so that, at long distances, $g(\mathbf{r})$ tends to one). Given these definitions, one easily finds that

$$C_{2+2}(\mathbf{r}, t_1, t_2) = 4 \overline{n^2(t_1, t_2) [g(\mathbf{r}) - 1]}. \quad (\text{A18})$$

In other words, $\zeta$ is the length scale on which defects are correlated. Only when $\zeta(t_2, t_1) \approx \xi_{\text{micro}}(t_1)$ does the configuration at time $t_2$ start to differ *structurally* from the configuration at the earlier time $t_1$.

Finally, let us mention that a length analogous to $\zeta(t_1, t_2)$ can be obtained with the analysis tools of temperature chaos, see Supplementary Note IV.

### 5. Computation of the chaotic parameter

As we explained in the main text, our goal here is to introduce a correlation parameter that will allow us to compare two different thermal protocols. This comparison should necessarily be local in space. We adapt to that end the procedure introduced in Ref. [38].

Specifically, we select $N_{\text{sph}} = 8000$ spheres of radius $R$ randomly chosen inside the system and centered at the central points of the elementary cells of the cubic lattice. Now, let us consider two identical systems that are subjected to two different thermal protocols, which we may name protocols $A_1$ and $A_2$. Next one performs a set of independent simulations (i.e., *repli-cas*) for protocol $A_1$, and another set of independent
reference time $t_w^*$ is the longest available waiting time for our native runs at the measuring temperature $T_m$. For the sake of clarity, we omit in Eqs. (A8) and (A9) the explicit dependence of $\theta$ on $\xi_{\text{micro}}$ (which is dealt with as explained in Ref. [7]).

### 4. Numerical coherence lengths $\xi_{\text{micro}}$ and $\zeta$

In this paragraph, we shall consider two more length scales. One of them, $\xi_{\text{micro}}$, is computed from the correlation function for the spin-glass order parameter (hence, $\xi_{\text{micro}}$ tells us about the size of the glassy domains). The second length scale, $\zeta(t_1, t_2)$, tells us about how the system reorganizes itself when going from the earlier time $t_1$ to the later time $t_2$.

#### a. The computation of $\xi_{\text{micro}}$

For the reader's convenience, let us recall the definition of the spatial autocorrelation function that we use for computing $\xi_{\text{micro}}(t_w)$ [41]

$$C_4(\mathbf{r}, t'; T) = \overline{\langle q^{(a,b)}(\mathbf{x}, t')q^{(a,b)}(\mathbf{x} + \mathbf{r}, t') \rangle_T}, \quad (\text{A10})$$

$$q^{(a,b)}(\mathbf{x}, t') \equiv \sigma^{(a)}(\mathbf{x}, t') \sigma^{(b)}(\mathbf{x}, t'), \quad (\text{A11})$$

where $t' = t_w + t$, the indices $(a, b)$ label different real replicas and $\langle \cdots \rangle_T$ stands for the average over the thermal noise at temperature $T$.

The calculation of the correlation function is computationally costly since we have $N_R(N_R - 1)/2$ possible choices of the pair of replicas. Fortunately, it can be accelerated using the specific multispin coding methods explained in Ref. [48].

Once we have $C_4(\mathbf{r}, t'; T)$, we compute the integrals [27, 28, 41]:

$$I_k(t'; T) = \int_0^\infty d^3\mathbf{r} r^k C_4(\mathbf{r} = (\mathbf{r}, 0, 0), t'; T). \quad (\text{A12})$$

A coherence length can be computed as

$$\xi_{k,k+1}(t', T) = \frac{I_{k+1}(t', T)}{I_k(t', T)}. \quad (\text{A13})$$

We define $\xi_{\text{micro}}(t, t_w; H) = \xi_{12}(t, t_w; H)$.

#### b. The $\zeta$ length scale

This length scale was studied in details in Ref. [41] by refining earlier suggestions [42, 43].

Let us consider the thermal trajectory followed by a given replica at the two times $t_1 < t_2$. Our basic quantity will be the local correlation

$$c_{\mathbf{x}}(t_1, t_2) = s_{\mathbf{x}}(t_2) s_{\mathbf{x}}(t_1). \quad (\text{A14})$$

Note that $c_{\mathbf{x}}(t_1, t_2) = -1$ if the spin at site $\mathbf{x}$ has been flipped when going from time $t_1$ to time $t_2$ [otherwise, $c_{\mathbf{x}}(t_1, t_2) = 1$]. Then, the two-time, two-site correlation function is

$$C_{2+2}(\mathbf{r}, t_1, t_2) = \frac{1}{N} \sum_{\mathbf{x}} \overline{[\langle c_{\mathbf{x}}(t_1, t_2) c_{\mathbf{x}+\mathbf{r}}(t_1, t_2) \rangle - C^2(t_1, t_2)]}, \quad (\text{A15})$$

where

$$C(t_1, t_2) = \frac{1}{N} \sum_{\mathbf{x}} \langle c_{\mathbf{x}}(t_1, t_2) \rangle. \quad (\text{A16})$$

The ideal $\zeta(t_1, t_2)$ is defined from the long-distance decay of $C_{2+2}(\mathbf{r}, t, t_w)$:

$$C_{2+2}(\mathbf{r}, t_1, t_2) \sim \frac{1}{rb} g(r/\zeta(t_1, t_2)), \quad (\text{A17})$$

where $g$ is an unknown scaling function. We bypass our lack of knowledge of $g$ exactly as we solved this problem for $\xi_{\text{micro}}$: by using integral estimators, recall Eq. (A13). Note that, by construction, $\zeta(t_1, t_2)$ tends to zero when $t_2$ approaches $t_1$. Conversely, we expect $\zeta(t_1, t_2)$ to grow with the later time $t_2$.

As for the interpretation of the length scale $\zeta$, an analogy with the theory of liquids is of help. We name a defect a site where $c_{\mathbf{x}}(t_1, t_2) = -1$. Let $n(t_1, t_2)$ be the density of defects $[C(t_1, t_2) = 1 - 2n(t_1, t_2)]$ and let $g(\mathbf{r})$ be the pair-correlation function for defects: The conditional probability for having a defect at site $\mathbf{x}+\mathbf{r}$, given that a defect is present at site $\mathbf{x}$, is $n(t_1, t_2)g(\mathbf{r})$ (so that, at long distances, $g(\mathbf{r})$ tends to one). Given these definitions, one easily finds that

$$C_{2+2}(\mathbf{r}, t_1, t_2) = 4 \overline{n^2(t_1, t_2) [g(\mathbf{r}) - 1]}. \quad (\text{A18})$$

In other words, $\zeta$ is the length scale on which defects are correlated. Only when $\zeta(t_2, t_1) \approx \xi_{\text{micro}}(t_1)$ does the configuration at time $t_2$ start to differ structurally from the configuration at the earlier time $t_1$.

Finally, let us mention that a length analogous to $\zeta(t_1, t_2)$ can be obtained with the analysis tools of temperature chaos, see Supplementary Note IV.

### 5. Computation of the chaotic parameter

As we explained in the main text, our goal here is to introduce a correlation parameter that will allow us to compare two different thermal protocols. This comparison should necessarily be local in space. We adapt to that end the procedure introduced in Ref. [38].

Specifically, we select $N_{\text{sph}} = 8000$ spheres of radius $R$ randomly chosen inside the system and centered at the central points of the elementary cells of the cubic lattice. Now, let us consider two identical systems that are subjected to two different thermal protocols, which we may name protocols $A_1$ and $A_2$. Next one performs a set of independent simulations (i.e., replicas) for protocol $A_1$, and another set of independent
simulations for protocol A₂. Then, the correlation co-efficient for protocols A₁ and A₁ as computed on the k-th sphere of radius R is defined as

$$X_{A_1, A_2}^{k,R} = \frac{\langle [q_{A_1, A_2}^{k,R}]^2 \rangle_T}{\sqrt{\langle [q_{A_1, A_1}^{k,R}]^2 \rangle_T \langle [q_{A_2, A_2}^{k,R}]^2 \rangle_T}} . \quad (A19)$$

In the above expression, $q_{A_1, A_2}^{k,R}$ is the overlap between two replicas $\sigma$ and $\tau$ that have undergone thermal protocols A₁ and A₂ respectively

$$q_{A_1, A_2}^{k,R} = \frac{1}{N_r} \sum_{\mathbf{x} \in B_R^k} s_{\mathbf{x}}^{\sigma, A_1} s_{\mathbf{x}}^{\tau, A_2}, \quad (A20)$$

where $N_R$ is the number of spins within the k-th sphere $B_R^k$ of radius R.

The interpretation of the chaotic parameter is very similar to a correlation coefficient: if $X_{A_1, A_2}^{k,R} = 1$, spin configurations from thermal protocols A₁ and A₂ are completely indistinguishable inside the sphere $B_R^k$ (absence of chaos). Instead, $X_{A_1, A_2}^{k,R} = 0$ corresponds to completely different configurations, which is an extremely chaotic situation.

The reader may notice from Eq. (A19) that the computation of $X_{A_1, A_2}^{k,R}$ involves an exact thermal expectation value (which could be obtained in simulations only if one had simulated an infinite number of replicas). Unfortunately, we only have $N_R^{\max} = 512$ replicas at our disposal. Our choice has been to produce different estimates of $X_{A_1, A_2, N_R}^{k,R}$ by varying $N_R$. Specifically, our procedure has been the following:

1. For each $N_R < N_R^{\max}$ we randomly order the $N_R^{\max}$ replicas and divide them in $N_R^{\max}/N_R$ groups of $N_R$ replicas.

2. In this way, we get $N_R^{\max}/N_R$ independent estimates of $X_{A_1, A_2, N_R}^{k,R}$.

3. In order to erase the effect of the initial permutation of the $N_R^{\max}$ replicas, we repeat this procedure 10 times for all $N_R < N_R^{\max}$.

In a nutshell, for every sphere of radius R we obtain $N_{\text{thermal}}(N_R)$ estimates of $X_{A_1, A_2, N_R}^{k,R}$ where

$$N_{\text{thermal}}(N_R < N_R^{\max}) = 10 \times \frac{N_R^{\max}}{N_R}, \quad (A21)$$

or

$$N_{\text{thermal}}(N_R = N_R^{\max}) = 1. \quad (A22)$$

We average the $N_{\text{thermal}}$ estimates of $X_{A_1, A_2, N_R}^{k,R}$ for every $N_R$ and finally, in a complete analogy with Ref. [38], we compute the extrapolation of the chaotic parameter to an infinite number of replicas by means of a simple linear extrapolation

$$X_{A_1, A_2, N_R}^{k,R} = X_{A_1, A_2, \infty}^{k,R} + \frac{A_{A_1, A_2}^{k,R}}{N_R}, \quad (A23)$$

where $X_{A_1, A_2, \infty}^{k,R}$ is our best estimation of $X_{A_1, A_2}^{k,R}$. More complicated extrapolations do not seem to present advantages (see SI in [38]).

Finally, in order to explore the statistical information carried by the $N_{\text{sph}} = 8000$ spheres, we define the distribution function

$$F(\tilde{X}, A_1, A_2, R) = \text{Probability}[X_{A_1, A_2}^{k,R} < \tilde{X}]. \quad (A24)$$

Some examples of this distribution function are displayed in Fig. 4 in the main text.

## ACKNOWLEDGMENTS

We acknowledge the precious contributions and ideas of our dear late friend and collaborator Raffaele Tripiccione. The Janus project would have been impossible without Lele's technical expertise, good sense and kindness, and we dedicate this work to him.

We thank Prof. R. Orbach for discussions.

This work was partly supported by grants No. PID2020-112936GB-I00, PID2019-103939RB-I00, No. PGC2018-094684-B-C21 and No. PGC2018-094684-B-C22 funded by Ministerio de Economía y Competitividad, Agencia Estatal de Investigación and Fondo Europeo de Desarrollo Regional (FEDER) (Spain and European Union), by grants No. GR21014 and No. IB20079 (partially funded by FEDER) funded by Junta the Extremadura (Spain), and by the Atracción de Talento program (Ref. 2019-T1/TIC-12776) funded by Comunidad de Madrid and Universidad Complutense de Madrid (Spain). This project has received funding from the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant No. 694925, G. Parisi). IGAP was supported by MCIU (Spain) through FPU Grant No. FPU18/02665. JMG was supported by the Ministerio de Universidades and the European Union "NextGeneration EU/PRTR" through 2021-2023 Margarita Salas grant.

## Author contributions

D.I. and A.T. contributed to the design of the Janus II project. J.M.G.-N. and D.N. contributed Janus II/Janus simulation software. M.B.-J., E.C., A.C., L.A.F, J.M.G.-N., I.G.-A.P., A.G.-G., D.I., A.M., A.M.-S., I.P., S.P.-G., S.F.S. and A.T. contributed to Janus II hardware and software development. L.A.F., E.M., V.M.-M. and I.P. suggested undertaking this project. L.A.F., E.M., V.M.-M., I.-P., F.R.-T. and J.J. R.-L. designed the research. J.M.-G and I.P. analyzed the data. M.B.-J., L.A.F., E.M., V.M.-M., J.M.-G., I.P., G.P., B.S., J.J.R.-L., F.R.-T. and D.Y. discussed the results. L.A.F., E.M., V.M.-M., J.M.-G., I.P., J.J.R.-L., B.S., F.R.-T. and D.Y. wrote the paper.
simulations for protocol A₂. Then, the correlation co-efficient for protocols A₁ and A₁ as computed on the k-th sphere of radius R is defined as

$$X_{A_1, A_2}^{k,R} = \frac{\langle [q_{A_1, A_2}^{k,R}]^2 \rangle_T}{\sqrt{\langle [q_{A_1, A_1}^{k,R}]^2 \rangle_T \langle [q_{A_2, A_2}^{k,R}]^2 \rangle_T}} . \quad (A19)$$

In the above expression, $q_{A_1, A_2}^{k,R}$ is the overlap between two replicas $\sigma$ and $\tau$ that have undergone thermal protocols A₁ and A₂ respectively

$$q_{A_1, A_2}^{k,R} = \frac{1}{N_r} \sum_{\mathbf{x} \in B_R^k} s_{\mathbf{x}}^{\sigma, A_1} s_{\mathbf{x}}^{\tau, A_2}, \quad (A20)$$

where $N_R$ is the number of spins within the k-th sphere $B_R^k$ of radius R.

The interpretation of the chaotic parameter is very similar to a correlation coefficient: if $X_{A_1, A_2}^{k,R} = 1$, spin configurations from thermal protocols A₁ and A₂ are completely indistinguishable inside the sphere $B_R^k$ (absence of chaos). Instead, $X_{A_1, A_2}^{k,R} = 0$ corresponds to completely different configurations, which is an extremely chaotic situation.

The reader may notice from Eq. (A19) that the computation of $X_{A_1, A_2}^{k,R}$ involves an exact thermal expectation value (which could be obtained in simulations only if one had simulated an infinite number of replicas). Unfortunately, we only have $N_R^{\max} = 512$ replicas at our disposal. Our choice has been to produce different estimates of $X_{A_1, A_2, N_R}^{k,R}$ by varying $N_R$. Specifically, our procedure has been the following:

1. For each $N_R < N_R^{\max}$ we randomly order the $N_R^{\max}$ replicas and divide them in $N_R^{\max}/N_R$ groups of $N_R$ replicas.

2. In this way, we get $N_R^{\max}/N_R$ independent estimates of $X_{A_1, A_2, N_R}^{k,R}$.

3. In order to erase the effect of the initial permutation of the $N_R^{\max}$ replicas, we repeat this procedure 10 times for all $N_R < N_R^{\max}$.

In a nutshell, for every sphere of radius R we obtain $N_{\text{thermal}}(N_R)$ estimates of $X_{A_1, A_2, N_R}^{k,R}$ where

$$N_{\text{thermal}}(N_R < N_R^{\max}) = 10 \times \frac{N_R^{\max}}{N_R}, \quad (A21)$$

or

$$N_{\text{thermal}}(N_R = N_R^{\max}) = 1. \quad (A22)$$

We average the $N_{\text{thermal}}$ estimates of $X_{A_1, A_2, N_R}^{k,R}$ for every $N_R$ and finally, in a complete analogy with Ref. [38], we compute the extrapolation of the chaotic parameter to an infinite number of replicas by means of a simple linear extrapolation

$$X_{A_1, A_2, N_R}^{k,R} = X_{A_1, A_2, \infty}^{k,R} + \frac{A_{A_1, A_2}^{k,R}}{N_R}, \quad (A23)$$

where $X_{A_1, A_2, \infty}^{k,R}$ is our best estimation of $X_{A_1, A_2}^{k,R}$. More complicated extrapolations do not seem to present advantages (see SI in [38]).

Finally, in order to explore the statistical information carried by the $N_{\text{sph}} = 8000$ spheres, we define the distribution function

$$F(\tilde{X}, A_1, A_2, R) = \text{Probability}[X_{A_1, A_2}^{k,R} < \tilde{X}]. \quad (A24)$$

Some examples of this distribution function are displayed in Fig. 4 in the main text.

## ACKNOWLEDGMENTS

We acknowledge the precious contributions and ideas of our dear late friend and collaborator Raffaele Tripiccione. The Janus project would have been impossible without Lele's technical expertise, good sense and kindness, and we dedicate this work to him.

We thank Prof. R. Orbach for discussions.

This work was partly supported by grants No. PID2020-112936GB-I00, PID2019-103939RB-I00, No. PGC2018-094684-B-C21 and No. PGC2018-094684-B-C22 funded by Ministerio de Economía y Competitividad, Agencia Estatal de Investigación and Fondo Europeo de Desarrollo Regional (FEDER) (Spain and European Union), by grants No. GR21014 and No. IB20079 (partially funded by FEDER) funded by Junta the Extremadura (Spain), and by the Atracción de Talento program (Ref. 2019-T1/TIC-12776) funded by Comunidad de Madrid and Universidad Complutense de Madrid (Spain). This project has received funding from the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant No. 694925, G. Parisi). IGAP was supported by MCIU (Spain) through FPU Grant No. FPU18/02665. JMG was supported by the Ministerio de Universidades and the European Union "NextGeneration EU/PRTR" through 2021-2023 Margarita Salas grant.

## Author contributions

D.I. and A.T. contributed to the design of the Janus II project. J.M.G.-N. and D.N. contributed Janus II/Janus simulation software. M.B.-J., E.C., A.C., L.A.F., J.M.G.-N., I.G.-A.P., A.G.-G., D.I., A.M., A.M.-S., I.P., S.P.-G., S.F.S. and A.T. contributed to Janus II hardware and software development. L.A.F., E.M., V.M.-M. and I.P. suggested undertaking this project. L.A.F., E.M., V.M.-M., I.-P., F.R.-T. and J.J. R.-L. designed the research. J.M.-G and I.P. analyzed the data. M.B.-J., L.A.F., E.M., V.M.-M., J.M.-G., I.P., G.P., B.S., J.J.R.-L., F.R.-T. and D.Y. discussed the results. L.A.F., E.M., V.M.-M., J.M.-G., I.P., J.J.R.-L., B.S., F.R.-T. and D.Y. wrote the paper.
**Data availability**

The data contained in the figures of this pa-
per, accompanied by the gnuplot script files that
generate these figures, are publicly available a
https://github.com/janusII/Rejuvenation_memory.git.
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 reason-
able request
**Data availability**

The data contained in the figures of this pa-
per, accompanied by the gnuplot script files that
generate these figures, are publicly available a
https://github.com/janusII/Rejuvenation_memory.git.
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 reason-
able request
Supplementary Information. Memory and rejuvenation in spin glasses:
aging systems are ruled by more than one length-scale

M. Baity-Jesi,¹ E. Calore,² A. Cruz,³,⁴ L.A. Fernandez,⁵,⁴ J.M. Gil-Narvion,⁴ I. Gonzalez-Adalid Pemartin,⁵

A. Gordillo-Guerrero,⁶,⁷,⁴ D. Iñiguez,⁴,⁸ A. Maiorano,⁹,⁴ E. Marinari,¹⁰ V. Martin-Mayor,⁵,⁴

J. Moreno-Gordo,⁴,³,¹¹,⁷ A. Muñoz Sudupe,⁵,⁴ D. Navarro,¹² I. Paga,¹³,* G. Parisi,¹⁰ S. Perez-Gaviro,³,⁴

F. Ricci-Tersenghi,¹⁰ J.J. Ruiz-Lorenzo,¹¹,⁷,⁴ S.F. Schifano,¹⁴ B. Seoane,⁵,⁴ A. Tarancon,³,⁴ and D. Yllanes¹⁵,⁴

(Janus Collaboration)

¹Eawag, Überlandstrasse 133, CH-8600 Dübendorf, Switzerland

²Dipartimento di Fisica e Scienze della Terra,
Università di Ferrara and INFN, Sezione di Ferrara, I-44122 Ferrara, Italy

³Departamento de Física Teórica, Universidad de Zaragoza, 50009 Zaragoza, Spain

⁴Instituto de Biocomputación y Física de Sistemas Complejos (BIFI), 50018 Zaragoza, Spain

⁵Departamento de Física Teórica, Universidad Complutense, 28040 Madrid, Spain

⁶Departamento de Ingeniería Eléctrica, Electrónica y Automática, U. de Extremadura, 10003, Cáceres, Spain

⁷Instituto de Computación Científica Avanzada (ICCAEx),
Universidad de Extremadura, 06006 Badajoz, Spain

⁸Fundación ARAID, Diputación General de Aragón, 50018 Zaragoza, Spain

⁹Dipartimento di Biotecnologie, Chimica e Farmacia, Università degli studi di Siena,
53100 Siena, Italy and INFN, Sezione di Roma 1, 00185 Rome, Italy

¹⁰Dipartimento di Fisica, Sapienza Università di Roma, and CNR-Nanotec,
Rome unit and INFN, Sezione di Roma 1, 00185 Rome, Italy

¹¹Departamento de Física, Universidad de Extremadura, 06006 Badajoz, Spain

¹²Departamento de Ingeniería, Electrónica y Comunicaciones and I3A, U. de Zaragoza, 50018 Zaragoza, Spain

¹³Dipartimento di Fisica, Sapienza Università di Roma,
and CNR-Nanotec, Rome unit, Sezione di Roma 1, 00185 Rome, Italy

¹⁴Dipartimento di Scienze Chimiche e Farmaceutiche,
Università di Ferrara e INFN Sezione di Ferrara, I-44122 Ferrara, Italy

¹⁵Chan Zuckerberg Biohub, San Francisco, CA, 94158

(Dated: August 31, 2022)

SUPPLEMENTARY NOTE I. THE
MAGNETIC RESPONSE IS SAMPLE
INDEPENDENT

The main quantities used in our numerical experi-
ments are the relaxation function, S<sub>ZFC</sub>(t, t<sub>w</sub>; H) and
the time at which it peaks, the so-called effective time
t<sub>H</sub><sup>eff</sup> ~ t<sub>w</sub>. The effective time can be treated as the
bridge to connect the macroscopic feature of the sys-
tem to the microscopic one, see Methods. In exper-
iments, all these quantities are computed for a given
disorder, and this is the approach we follow as well.
However, one may wonder about the variability caused
by the particular disorder realization one is looking at.
This is the question that we address in this Supple-
mentary Note.

As a particularly relevant example, let us consider
the position of the peak of the relaxation function for
a given magnetic field *H*, namely *t*<sub>*H*</sub><sup>*eff*</sup>. We have sim-
ulated four independent samples, i.e. disorder real-
izations, that we name *S*<sub>*i*</sub> for *i* = 1, 2, 3 and 4. For
each sample, we extracted the effective time *t*<sub>*S*<sub>*i*</sub>,*H*</sub><sup>*eff*</sup> fol-
lowing the techniques introduced in Refs. [1, 2]. We
specialize to the relevant case of the *H*-dependence of
decay of the ratio log(*t*<sub>*S*<sub>*i*</sub>,*H*</sub><sup>*eff*</sup>/log(*t*<sub>*S*<sub>*i*</sub>,*H*</sub><sup>*eff*</sup>)) for *i* = 1, .., *N*<sub>S</sub>
as computed for the native protocols at our coldest

temperature ($T = 0.5$) and $t_w = t_w^\downarrow + 2^{31.25}$. As the
reader can check in Fig. 1, the data for the different
samples turned out to be statistically compatible.

<figure><img src="image_7.png" /><figcaption>FIG. 1. Logarithm of the ratio of effective times $\log(t_{S_{i,H}}^{\text{eff}}/t_{S_{i,H}\to 0+}^{\text{eff}})$ versus squared magnetic field, as computed for our four independent samples S$_{1}$, S$_{2}$, S$_{3}$ and S$_{4}$. The system undergoes a Zero-Field Cooled protocol (i.e. a native protocol) at temperature $T = 0.5$ and waiting time $t_w = t_w^\downarrow + 2^{31.25}$ (recall that $t_w = 2^{31.25}$). For each sample, we construct the relaxation function, $S_{ZFC}^S(t, t_w; H)$, and extract from it the $t_{S_{i,H}}^{\text{eff}}$ values. Errors are computed with a jackknife method applied to the 512 replicas that we simulate for every sample.</figcaption></figure>

<footnote>* ilaria.paga@gmail.com</footnote>
Supplementary Information. Memory and rejuvenation in spin glasses:
aging systems are ruled by more than one length-scale

M. Baity-Jesi,¹ E. Calore,² A. Cruz,³,⁴ L.A. Fernandez,⁵,⁴ J.M. Gil-Narvion,⁴ I. Gonzalez-Adalid Pemartin,⁵

A. Gordillo-Guerrero,⁶,⁷,⁴ D. Iñiguez,⁴,⁸ A. Maiorano,⁹,⁴ E. Marinari,¹⁰ V. Martin-Mayor,⁵,⁴

J. Moreno-Gordo,⁴,³,¹¹,⁷ A. Muñoz Sudupe,⁵,⁴ D. Navarro,¹² I. Paga,¹³,* G. Parisi,¹⁰ S. Perez-Gaviro,³,⁴

F. Ricci-Tersenghi,¹⁰ J.J. Ruiz-Lorenzo,¹¹,⁷,⁴ S.F. Schifano,¹⁴ B. Seoane,⁵,⁴ A. Tarancon,³,⁴ and D. Yllanes¹⁵,⁴

(Janus Collaboration)

¹Eawag, Überlandstrasse 133, CH-8600 Dübendorf, Switzerland

²Dipartimento di Fisica e Scienze della Terra,
Università di Ferrara and INFN, Sezione di Ferrara, I-44122 Ferrara, Italy

³Departamento de Física Teórica, Universidad de Zaragoza, 50009 Zaragoza, Spain

⁴Instituto de Biocomputación y Física de Sistemas Complejos (BIFI), 50018 Zaragoza, Spain

⁵Departamento de Física Teórica, Universidad Complutense, 28040 Madrid, Spain

⁶Departamento de Ingeniería Eléctrica, Electrónica y Automática, U. de Extremadura, 10003, Cáceres, Spain

⁷Instituto de Computación Científica Avanzada (ICCAEx),
Universidad de Extremadura, 06006 Badajoz, Spain

⁸Fundación ARAID, Diputación General de Aragón, 50018 Zaragoza, Spain

⁹Dipartimento di Biotecnologie, Chimica e Farmacia, Università degli studi di Siena,
53100 Siena, Italy and INFN, Sezione di Roma 1, 00185 Rome, Italy

¹⁰Dipartimento di Fisica, Sapienza Università di Roma, and CNR-Nanotec,
Rome unit and INFN, Sezione di Roma 1, 00185 Rome, Italy

¹¹Departamento de Física, Universidad de Extremadura, 06006 Badajoz, Spain

¹²Departamento de Ingeniería, Electrónica y Comunicaciones and I3A, U. de Zaragoza, 50018 Zaragoza, Spain

¹³Dipartimento di Fisica, Sapienza Università di Roma,
and CNR-Nanotec, Rome unit, Sezione di Roma 1, 00185 Rome, Italy

¹⁴Dipartimento di Scienze Chimiche e Farmaceutiche,
Università di Ferrara e INFN Sezione di Ferrara, I-44122 Ferrara, Italy

¹⁵Chan Zuckerberg Biohub, San Francisco, CA, 94158

(Dated: August 31, 2022)

SUPPLEMENTARY NOTE I. THE
MAGNETIC RESPONSE IS SAMPLE
INDEPENDENT

The main quantities used in our numerical experi-
ments are the relaxation function, S<sub>ZFC</sub>(t, t<sub>w</sub>; H) and
the time at which it peaks, the so-called effective time
t<sub>H</sub><sup>eff</sup> ~ t<sub>w</sub>. The effective time can be treated as the
bridge to connect the macroscopic feature of the sys-
tem to the microscopic one, see **Methods**. In exper-
iments, all these quantities are computed for a given
disorder, and this is the approach we follow as well.
However, one may wonder about the variability caused
by the particular disorder realization one is looking at.
This is the question that we address in this Supple-
mentary Note.

As a particularly relevant example, let us consider
the position of the peak of the relaxation function for
a given magnetic field *H*, namely *t*<sub>*H*</sub><sup>*eff*</sup>. We have sim-
ulated four independent samples, i.e. disorder real-
izations, that we name *S*<sub>*i*</sub> for *i* = 1, 2, 3 and 4. For
each sample, we extracted the effective time *t*<sub>*S*<sub>*i*</sub>,*H*</sub><sup>*eff*</sup> fol-
lowing the techniques introduced in Refs. [1, 2]. We
specialize to the relevant case of the *H*-dependence of
decay of the ratio log(*t*<sub>*S*<sub>*i*</sub>,*H*</sub><sup>*eff*</sup>/log(*t*<sub>*S*<sub>*i*</sub>,*H*</sub><sup>*eff*</sup>)) for *i* = 1, .., *N*<sub>S</sub>
as computed for the native protocols at our coldest

temperature ($T = 0.5$) and $t_w = t_w^\perp + 2^{31.25}$. As the
reader can check in Fig. 1, the data for the different
samples turned out to be statistically compatible.

<figure><img src="image_6.png" /><figcaption>FIG. 1. Logarithm of the ratio of effective times $\log(t_{S_{i,H}}^{\text{eff}}/t_{S_{i,H}\to 0+}^{\text{eff}})$ versus squared magnetic field, as computed for our four independent samples S$_1$, S$_2$, S$_3$ and S$_4$. The system undergoes a Zero-Field Cooled protocol (i.e. a native protocol) at temperature $T = 0.5$ and waiting time $t_w = t_w^{\perp} + 2^{31.25}$ (recall that $t_w = 2^{31.25}$). For each sample, we construct the relaxation function, $S_{ZFC}^i(t, t_w; H)$, and extract from it the $t_{S_{i,H}}^{\text{eff}}$ values. Errors are computed with a jackknife method applied to the 512 replicas that we simulate for every sample.</figcaption></figure>

<footnote>* ilaria.paga@gmail.com</footnote>
## SUPPLEMENTARY NOTE II. THE ZERO-FIELD COOLING (ZFC) NUMERICAL EXPERIMENT WITH A JUMP TO $T_2 = 0.7$

In Figure 1—main text, we have illustrated the thermal protocol in our numerical experiments for the cold temperature $T_2 = 0.5$. Now, the reader may rightly wonder about how the relaxation function looks like when the temperature jump does not meet the temperature-chaos requirement expressed by Eq. (2)—main text. The answer to this question can be found in Fig. 2, where the same thermal protocol was considered, but we simulated a more modest temperature jump from $T_1 = 0.9$ to $T_2 = 0.7$, which does not meet the chaos’ requirement.

In agreement with the results in the main text, rejuvenation in the jump to $T_2 = 0.7$ turns out to be only partial: even for the shortest $t_w$, the aging peak is visible in the form of a shoulder when jumping to $T_2 = 0.7$, see Fig. 2—central panel. Instead, the aging peak peaks is absent in the native protocol (not even a shoulder is visible for $t_w = 2^{10}$ in Fig. 2—left—bottom).

## SUPPLEMENTARY NOTE III. ON THE GROWTH OF THE COHERENCE LENGTH, $\zeta(t_1, t_2)$

As explained in the main text and in Methods, the length scale $\zeta(t_1, t_2)$ allows us to compare quantitatively the spin configurations at the two times $t_1 < t_2$. Indeed, $\zeta(t_1, t_2)$ is the length scale at which significant changes occurred. It is clear, then, that $\zeta$ is very important to rationalize the memory effect found in our numerical experiments. However, the time growth of $\zeta$ may be intricate. In this note we illustrate some features of this growth for the two main protocols considered in this work, namely the native protocol and the temperature jump protocol.

In fact, for jump protocols only and in shorter length scales, the growth of $\zeta(t_1=t_w, t_2=t_w+t)$ was studied in Ref. [4]. For $t \ll t_w$, the time growth of $\zeta$ turned out to be power law, with the same exponent found for $\xi_{micro}$. Yet, when $t$ approached $t_w$ it was found that $\zeta$ saturated at a value approximately equal to $\xi_{micro}(t_w)$ [given the way we compute it, we expect that $\zeta(t_1, t_2)$ will be upper-bounded by the smallest of $\xi_{micro}(t_1)$ and $\xi_{micro}(t_2)$]. Our aim here will be twofold. First, we shall check for jump protocols whether or not these findings to the larger length scales reached in this work. Our second goal will be presenting, for the first time we believe, results on $\zeta$ as computed for temperature-jump protocols.

Indeed, the two thermal protocols are compared on figure 3. We find rather different time-behaviours for both kinds of thermal histories:

* For native protocols, the curves for the different $t_w$ start to significantly differ for time $t \sim t_w$, when a saturation value is approached.

* For our jump protocols instead, the curves corresponding to the different waiting time are almost statistically compatible.

In other words, native protocols behave in agreement with our expectations from older simulations [4]. Indeed, Fig. 4 confirms that $\zeta(t_1 = t_w, t_2 = t_1+t)$ saturates to the $\xi_{micro}(t_w)$ value. Instead, for the jump-protocols, $\zeta(t_1, t_2)$ is significantly smaller than $\xi_{micro}$ (and, accordingly, data show no sign of saturation).

## SUPPLEMENTARY NOTE IV. RE-OBTAINING $\zeta$ WITH THE TOOLS OF TEMPERATURE CHAOS

Our goal here will be identifying the length scale $\zeta(t_1, t_2)$ discussed in the main text, in Methods and in the Supplementary Note II, using the analysis tools that we employ to discuss temperature chaos.

The main difference will be in that in the computation of $\zeta(t_1, t_2)$ the spin configurations of the single replica are compared at times $t_1$ and $t_2$. Instead, here we shall be considering two independent replicas, one taken at time $t_1$, the other at time $t_2$. The key point is that the distribution function for the chaotic correlation parameter $X$ turns out to be strongly dependent on the radius of the sphere (remember that $X$ is computed by restricting our attention to the interior points of randomly chosen spheres of radius $R$, see Methods). Now, it will be found that the distribution function extending to smaller values of $X$ (i.e. the strongest chaotic signal) is found, precisely, when $R \approx \zeta(t_1, t_2)$.

Specifically, we shall compare configurations for different times $t_w$ and for different sphere radius $R$. In order to lighten the notation we give a different label in Table I to all the protocols used to compute the chaotic correlation parameter $X$ (remember that $X$ is always computed for a pair of protocols).

In Fig. 5 we show the comparison between configurations at $T = 0.9$ just before the temperature jump against configurations at $T = 0.5$ for several times $t_w$ after the temperature jump. The analogous data for the jump to $T = 0.7$ can be seen in Fig. 6. For the shortest $t_w$ (when $\zeta$ is very small), i.e. protocols $B_5$ and $B_7$, most of the spheres have $X \approx 1$. We can only find a tail of smaller correlation parameters for spheres of radius $R = 1$. However, at later times $t_w$ after the temperature jump, corresponding to configurations $E_5$ and $E_7$, the relevant scale to observe the chaotic signal is $R = 3$ for the temperature jump to $T = 0.5$ and $R \in [3, 5]$ for the temperature jump to $T = 0.7$.

In conclusion, we find that the optimal $R$ for comparing the configurations right before the jump with configurations at time $t_w$ after the jump goes from $R \approx 1$ for very short $t_w$ to $R = 3$ (for our largest $t_w$ in the jump to 0.9) or $R \in [3, 5]$ for the temperature jump to $T = 0.7$. The comparison of these results with the values of $\zeta$ in Fig. 3 suggests that, indeed, the two
## SUPPLEMENTARY NOTE II. THE ZERO-FIELD COOLING (ZFC) NUMERICAL EXPERIMENT WITH A JUMP TO $T_2 = 0.7$

In Figure 1—main text, we have illustrated the thermal protocol in our numerical experiments for the cold temperature $T_2 = 0.5$. Now, the reader may rightly wonder about how the relaxation function looks like when the temperature jump does not meet the temperature-chaos requirement expressed by Eq. (2)—main text. The answer to this question can be found in Fig. 2, where the same thermal protocol was considered, but we simulated a more modest temperature jump from $T_1 = 0.9$ to $T_2 = 0.7$, which does not meet the chaos’ requirement.

In agreement with the results in the main text, rejuvenation in the jump to $T_2 = 0.7$ turns out to be only partial: even for the shortest $t_w$, the aging peak is visible in the form of a shoulder when jumping to $T_2 = 0.7$, see Fig. 2—central panel. Instead, the aging peak peaks is absent in the native protocol (not even a shoulder is visible for $t_w = 2^{10}$ in Fig. 2—left—bottom).

## SUPPLEMENTARY NOTE III. ON THE GROWTH OF THE COHERENCE LENGTH, $\zeta(t_1, t_2)$

As explained in the main text and in Methods, the length scale $\zeta(t_1, t_2)$ allows us to compare quantitatively the spin configurations at the two times $t_1 < t_2$. Indeed, $\zeta(t_1, t_2)$ is the length scale at which significant changes occurred. It is clear, then, that $\zeta$ is very important to rationalize the memory effect found in our numerical experiments. However, the time growth of $\zeta$ may be intricate. In this note we illustrate some features of this growth for the two main protocols considered in this work, namely the native protocol and the temperature jump protocol.

In fact, for jump protocols only and in shorter length scales, the growth of $\zeta(t_1=t_w, t_2=t_w+t)$ was studied in Ref. [4]. For $t \ll t_w$, the time growth of $\zeta$ turned out to be power law, with the same exponent found for $\xi_{micro}$. Yet, when $t$ approached $t_w$ it was found that $\zeta$ saturated at a value approximately equal to $\xi_{micro}(t_w)$ [given the way we compute it, we expect that $\zeta(t_1, t_2)$ will be upper-bounded by the smallest of $\xi_{micro}(t_1)$ and $\xi_{micro}(t_2)$]. Our aim here will be twofold. First, we shall check for jump protocols whether or not these findings to the larger length scales reached in this work. Our second goal will be presenting, for the first time we believe, results on $\zeta$ as computed for temperature-jump protocols.

Indeed, the two thermal protocols are compared on figure 3. We find rather different time-behaviours for both kinds of thermal histories:

* For native protocols, the curves for the different $t_w$ start to significantly differ for time $t \sim t_w$, when a saturation value is approached.

* For our jump protocols instead, the curves corresponding to the different waiting time are almost statistically compatible.

In other words, native protocols behave in agreement with our expectations from older simulations [4]. Indeed, Fig. 4 confirms that $\zeta(t_1 = t_w, t_2 = t_1 + t)$ saturates to the $\xi_{micro}(t_w)$ value. Instead, for the jump-protocols, $\zeta(t_1, t_2)$ is significantly smaller than $\xi_{micro}$ (and, accordingly, data show no sign of saturation).

## SUPPLEMENTARY NOTE IV. RE-OBTAINING $\zeta$ WITH THE TOOLS OF TEMPERATURE CHAOS

Our goal here will be identifying the length scale $\zeta(t_1, t_2)$ discussed in the main text, in Methods and in the Supplementary Note II, using the analysis tools that we employ to discuss temperature chaos.

The main difference will be in that in the computation of $\zeta(t_1, t_2)$ the spin configurations of the single replica are compared at times $t_1$ and $t_2$. Instead, here we shall be considering two independent replicas, one taken at time $t_1$, the other at time $t_2$. The key point is that the distribution function for the chaotic correlation parameter $X$ turns out to be strongly dependent on the radius of the sphere (remember that $X$ is computed by restricting our attention to the interior points of randomly chosen spheres of radius $R$, see Methods). Now, it will be found that the distribution function extending to smaller values of $X$ (i.e. the strongest chaotic signal) is found, precisely, when $R \approx \zeta(t_1, t_2)$.

Specifically, we shall compare configurations for different times $t_w$ and for different sphere radius $R$. In order to lighten the notation we give a different label in Table I to all the protocols used to compute the chaotic correlation parameter $X$ (remember that $X$ is always computed for a pair of protocols).

In Fig. 5 we show the comparison between configurations at $T = 0.9$ just before the temperature jump against configurations at $T = 0.5$ for several times $t_w$ after the temperature jump. The analogous data for the jump to $T = 0.7$ can be seen in Fig. 6. For the shortest $t_w$ (when $\zeta$ is very small), i.e. protocols $B_5$ and $B_7$, most of the spheres have $X \approx 1$. We can only find a tail of smaller correlation parameters for spheres of radius $R = 1$. However, at later times $t_w$ after the temperature jump, corresponding to configurations $E_5$ and $E_7$, the relevant scale to observe the chaotic signal is $R = 3$ for the temperature jump to $T = 0.5$ and $R \in [3, 5]$ for the temperature jump to $T = 0.7$.

In conclusion, we find that the optimal $R$ for comparing the configurations right before the jump with configurations at time $t_w$ after the jump goes from $R \approx 1$ for very short $t_w$ to $R = 3$ (for our largest $t_w$ in the jump to 0.9) or $R \in [3, 5]$ for the temperature jump to $T = 0.7$. The comparison of these results with the values of $\zeta$ in Fig. 3 suggests that, indeed, the two
<figure><img src="image_8.png" /></figure>

FIG. 2. The zero-field cooling (ZFC) numerical experiment at the cold temperature $T_2 = 0.7$. (This figure is the analog of Figure 1—main text, but the temperature jump is $T_1 = 0.9 \rightarrow T_2 = 0.7$). The initial random spin configuration is placed instantaneously at the working temperature and let relax for a time $t_w$ without field. At time $t_w$, a magnetic field $H = 0.01$ is applied and the magnetic density, $M_{ZFC}(t, t_w; H)$, is recorded. Left panels show the relaxation function $S_{ZFC}(t, t_w; H)$, for the native runs at the hotter, $T_1 = 0.9$, and colder, $T_2 = 0.7$, temperatures (both below the glass temperature $T_g = 1.102(3)[3]$). The physically interesting peak of the $S_{ZFC}(t, t_w; H)$ is the one at time $t_H^{\text{eff}} \approx t_w$. In our protocol (schematized by the green arrows), after a waiting time $t_w^{\downarrow} = 2^{31.25}$, the temperature is suddenly dropped from the initial temperature $T_1 = 0.9$ to the colder temperature $T_2 = 0.7$. Then, the system is let relax at $T_2$ for an additional time, after which the magnetic field is switched on and the function $S_{ZFC}(t, t_w; H)$ shown in the central-bottom panel is measured. Waiting times for these jump runs are reported in the legend. As we have explained in the Main, the cold temperature $T_2 = 0.7$ does not satisfy the chaos' requirement, and, so, the rejuvenation effect is weak (see Eq. (2) of the Main). In all cases, error bars are one standard deviation.

approaches are featuring the same length scale.

## SUPPLEMENTARY NOTE V. THE TEMPERATURE-CHAOS VIEWPOINT ON MEMORY

In this supplementary note, we show a long numerical simulation that allows us to observe memory effects from the Temperature Chaos point of view.

The studied simulation starts at temperature $T = 0.9$

where it ages for $t_w^{\downarrow} = 2^{31.25}$, then, the temperature of the simulation jumps to $T = 0.5$ where it ages for a time $t_w^{\uparrow} = 2^{31.25}$. Finally, the temperature of the thermal bath in the simulation raises up again to $T = 0.9$ where it ages for a time $t_w = 2^{31}$.

We select spheres of radius $R = 3$ and we compute the chaotic parameter restricted to those spheres, just as explained in the main text.

For the computation of the chaotic parameter, we select as a reference the configurations of the system at temperature $T = 0.9$ just before the temperature
<figure><img src="image_7.png" /></figure>

FIG. 2. The zero-field cooling (ZFC) numerical experiment at the cold temperature $T_2 = 0.7$. (This figure is the analog of Figure 1—main text, but the temperature jump is $T_1 = 0.9 \rightarrow T_2 = 0.7$). The initial random spin configuration is placed instantaneously at the working temperature and let relax for a time $t_w$ without field. At time $t_w$, a magnetic field $H = 0.01$ is applied and the magnetic density, $M_{ZFC}(t, t_w; H)$, is recorded. Left panels show the relaxation function $S_{ZFC}(t, t_w; H)$, for the native runs at the hotter, $T_1 = 0.9$, and colder, $T_2 = 0.7$, temperatures (both below the glass temperature $T_g = 1.102(3)[3]$). The physically interesting peak of the $S_{ZFC}(t, t_w; H)$ is the one at time $t_H^{\text{eff}} \approx t_w$. In our protocol (schematized by the green arrows), after a waiting time $t_w^{\downarrow} = 2^{31.25}$, the temperature is suddenly dropped from the initial temperature $T_1 = 0.9$ to the colder temperature $T_2 = 0.7$. Then, the system is let relax at $T_2$ for an additional time, after which the magnetic field is switched on and the function $S_{ZFC}(t, t_w; H)$ shown in the central-bottom panel is measured. Waiting times for these jump runs are reported in the legend. As we have explained in the Main, the cold temperature $T_2 = 0.7$ does not satisfy the chaos' requirement, and, so, the rejuvenation effect is weak (see Eq. (2) of the Main). In all cases, error bars are one standard deviation.

approaches are featuring the same length scale.

## SUPPLEMENTARY NOTE V. THE TEMPERATURE-CHAOS VIEWPOINT ON MEMORY

In this supplementary note, we show a long numerical simulation that allows us to observe memory effects from the Temperature Chaos point of view.

The studied simulation starts at temperature $T = 0.9$

where it ages for $t_w^{\downarrow} = 2^{31.25}$, then, the temperature of the simulation jumps to $T = 0.5$ where it ages for a time $t_w^{\uparrow} = 2^{31.25}$. Finally, the temperature of the thermal bath in the simulation raises up again to $T = 0.9$ where it ages for a time $t_w = 2^{31}$.

We select spheres of radius $R = 3$ and we compute the chaotic parameter restricted to those spheres, just as explained in the main text.

For the computation of the chaotic parameter, we select as a reference the configurations of the system at temperature $T = 0.9$ just before the temperature
<figure><img src="image_9.png" /><figcaption>FIG. 3. Growth of the length scale, $\zeta(t_1, t_2)$, in absence of an external magnetic field, for native protocols (triangles) and for temperature jumps (dots). In both cases, we present results for several waiting times $t_w$. For native protocols, a fully disordered system (i.e. infinite temperature) is placed at the initial time working temperature $T$ and let to evolve there ($t_1 = t_w$ and $t_2 = t + t_w$ for this protocols). Instead, in temperature jumps a fully disordered system is suddenly placed at $T = 0.9$ where it is let to relax for a time $t_w$. Next, the temperature is instantaneously lowered to $T$. For jump protocols, $t_1 = t_w + t_w$ and $t_2 = t_1 + t$. On the top, data obtained for the measuring temperature $T = 0.5$; on the bottom, data for at the hot temperature $T = 0.7$. Errors are obtained through Jack-knife propagation of the sample-to-sample fluctuations.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 4. The data for $\zeta(t_1 = t_w, t_2 = t + t_w)$ (triangles) from Fig. 3, native protocols only, are compared with the lengths $\xi_{\text{micro}}(t_2 = t + t_w)$ computed for the same simulations.</figcaption></figure>

jump to $T = 0.5$ which correspond with configuration $A_9$ in Table I, and we compare them with configurations spaced along the above described simulations, namely $B_5$, $C_5$, $D_5$, $E_5$, $B_9$, $C_9$, $D_9$, and $E_9$ in Table I.

The results of these computations are shown in Fig. 7. The initial situation, corresponding to the pure points, represents a non-chaotic situation in which configurations at $T = 0.9$ and $T = 0.5$ are similar (i.e. $X \approx 1$ for most spheres). When the color becomes lighter, the simulation at $T = 0.5$ ages, and

the configurations differ more and more from the reference configuration at $T = 0.9$. Indeed, the distribution function shows a small tail at smaller values of $X$ (however, even in the most extreme situation $x < 0.8$ only for 1% of the spheres).

When the temperature jumps back to $T = 0.9$, the distribution function evolves very quickly to a non-chaotic situation. The final situation, corresponding to the distribution function displayed with red points, is a perfect example of memory effects in Spin Glasses. The system at this point explores the same configurations that used to explore before the temperature jump, just as its excursion to $T = 0.5$ never happened.

## SUPPLEMENTARY NOTE VI. THE OPTIMAL SPHERE SIZE WHEN $\xi_{\text{micro}}$ IS TRULY DIFFERENT

In the main text (specifically, in the analysis leading to Fig. 4 in the main text) we needed to face a problem never encountered previously in numerical simulations. Specifically, we wanted to compare spin configurations from a thermal protocol producing $\xi_{\text{micro}} \approx 16.6 a_0$ with configurations with a protocol reaching $\xi_{\text{micro}} \approx 5.8 a_0$.

<table><caption>TABLE I. Main parameters for each of the numerical simulations. SG stands for Spin Glass, while DIM means Diluted Ising Model. $t_w^\perp$ corresponds to $2^{31.25}$ steps at $T = 0.9$, $t_w^\perp$ corresponds to $2^{31.25}$ steps at $T = 0.5$ as appears in the main text. For details about the thermal protocol see the main text, in particular Fig. 1, where we explicit the difference between the fixed-temperature protocol (native), the jump protocol, and the jump-back protocol.</caption><thead><tr><th></th><th>System</th><th><i>T</i></th><th>type</th><th><i>t</i><sub>w</sub></th></tr></thead><tbody><tr><td>A<sub>5</sub></td><td>SG</td><td>0.5</td><td>NATIVE</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>31.25</sup></td></tr><tr><td>B<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>10</sup></td></tr><tr><td>C<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>15.625</sup></td></tr><tr><td>D<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>23.5</sup></td></tr><tr><td>E<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>31.25</sup></td></tr><tr><td>A<sub>7</sub></td><td>SG</td><td>0.7</td><td>NATIVE</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>31.25</sup></td></tr><tr><td>B<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>10</sup></td></tr><tr><td>C<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>15.625</sup></td></tr><tr><td>D<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>23.5</sup></td></tr><tr><td>E<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>31.25</sup></td></tr><tr><td>A<sub>9</sub></td><td>SG</td><td>0.9</td><td>NATIVE</td><td><i>t</i><sub>w</sub><sup>&perp;</sup></td></tr><tr><td>B<sub>9</sub></td><td>SG</td><td>0.9</td><td>JUMP/BACK</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + <i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>10</sup></td></tr><tr><td>C<sub>9</sub></td><td>SG</td><td>0.9</td><td>JUMP/BACK</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + <i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>15.625</sup></td></tr><tr><td>D<sub>9</sub></td><td>SG</td><td>0.9</td><td>JUMP/BACK</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + <i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>23.5</sup></td></tr><tr><td>E<sub>9</sub></td><td>SG</td><td>0.9</td><td>JUMP/BACK</td><td><i>t</i><sub>w</sub><sup>&perp;</sup> + <i>t</i><sub>w</sub><sup>&perp;</sup> + 2<sup>31</sup></td></tr><tr><td>A'<sub>5</sub></td><td>DIM</td><td>0.5</td><td>NATIVE</td><td>76</td></tr><tr><td>B'<sub>5</sub></td><td>DIM</td><td>0.5</td><td>JUMP</td><td>430 + 69</td></tr><tr><td>A'<sub>7</sub></td><td>DIM</td><td>0.7</td><td>NATIVE</td><td>197</td></tr><tr><td>B'<sub>7</sub></td><td>DIM</td><td>0.7</td><td>JUMP</td><td>430 + 165</td></tr><tr><td>A'<sub>9</sub></td><td>DIM</td><td>0.9</td><td>NATIVE</td><td>430</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>FIG. 3. Growth of the length scale, $\zeta(t_1, t_2)$, in absence of an external magnetic field, for native protocols (triangles) and for temperature jumps (dots). In both cases, we present results for several waiting times $t_w$. For native protocols, a fully disordered system (i.e. infinite temperature) is placed at the initial time working temperature $T$ and let to evolve there ($t_1 = t_w$ and $t_2 = t + t_w$ for this protocols). Instead, in temperature jumps a fully disordered system is suddenly placed at $T = 0.9$ where it is let to relax for a time $t_w$. Next, the temperature is instantaneously lowered to $T$. For jump protocols, $t_1 = t_w + t_w$ and $t_2 = t_1 + t$. On the top, data obtained for the measuring temperature $T = 0.5$; on the bottom, data for at the hot temperature $T = 0.7$. Errors are obtained through Jack-knife propagation of the sample-to-sample fluctuations.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 4. The data for $\zeta(t_1 = t_w, t_2 = t + t_w)$ (triangles) from Fig. 3, native protocols only, are compared with the lengths $\xi_{\text{micro}}(t_2 = t + t_w)$ computed for the same simulations.</figcaption></figure>

jump to $T = 0.5$ which correspond with configuration $A_9$ in Table I, and we compare them with configurations spaced along the above described simulations, namely $B_5$, $C_5$, $D_5$, $E_5$, $B_9$, $C_9$, $D_9$, and $E_9$ in Table I.

The results of these computations are shown in Fig. 7. The initial situation, corresponding to the pure points, represents a non-chaotic situation in which configurations at $T = 0.9$ and $T = 0.5$ are similar (i.e. $X \approx 1$ for most spheres). When the color becomes lighter, the simulation at $T = 0.5$ ages, and

the configurations differ more and more from the reference configuration at $T = 0.9$. Indeed, the distribution function shows a small tail at smaller values of $X$ (however, even in the most extreme situation $x < 0.8$ only for 1% of the spheres).

When the temperature jumps back to $T = 0.9$, the distribution function evolves very quickly to a non-chaotic situation. The final situation, corresponding to the distribution function displayed with red points, is a perfect example of memory effects in Spin Glasses. The system at this point explores the same configurations that used to explore before the temperature jump, just as its excursion to $T = 0.5$ never happened.

## SUPPLEMENTARY NOTE VI. THE OPTIMAL SPHERE SIZE WHEN $\xi_{\text{micro}}$ IS TRULY DIFFERENT

In the main text (specifically, in the analysis leading to Fig. 4 in the main text) we needed to face a problem never encountered previously in numerical simulations. Specifically, we wanted to compare spin configurations from a thermal protocol producing $\xi_{\text{micro}} \approx 16.6 a_0$ with configurations with a protocol reaching $\xi_{\text{micro}} \approx 5.8 a_0$.

<table><caption>TABLE I. Main parameters for each of the numerical simulations. SG stands for Spin Glass, while DIM means Diluted Ising Model. $t_w^\perp$ corresponds to $2^{31.25}$ steps at $T = 0.9$, $t_w^\uparrow$ corresponds to $2^{31.25}$ steps at $T = 0.5$ as appears in the main text. For details about the thermal protocol see the main text, in particular Fig. 1, where we explicit the difference between the fixed-temperature protocol (native), the jump protocol, and the jump-back protocol.</caption><thead><tr><th></th><th>System</th><th>T</th><th>type</th><th>t<sub>w</sub></th></tr></thead><tbody><tr><td>A<sub>5</sub></td><td>SG</td><td>0.5</td><td>NATIVE</td><td>t<sub>w</sub><sup>&perp;</sup> + 2<sup>31.25</sup></td></tr><tr><td>B<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td>t<sub>w</sub><sup>&perp;</sup> + 2<sup>10</sup></td></tr><tr><td>C<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td>t<sub>w</sub><sup>&perp;</sup> + 2<sup>15.625</sup></td></tr><tr><td>D<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td>t<sub>w</sub><sup>&perp;</sup> + 2<sup>23.5</sup></td></tr><tr><td>E<sub>5</sub></td><td>SG</td><td>0.5</td><td>JUMP</td><td>t<sub>w</sub><sup>&perp;</sup> + 2<sup>31.25</sup></td></tr><tr><td>A<sub>7</sub></td><td>SG</td><td>0.7</td><td>NATIVE</td><td>t<sub>w</sub><sup>&perp;</sup> + 2<sup>31.25</sup></td></tr><tr><td>B<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td>t<sub>w</sub><sup>&perp;</sup> + 2<sup>10</sup></td></tr><tr><td>C<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td>t<sub>w</sub><sup>&perp;</sup> + 2<sup>15.625</sup></td></tr><tr><td>D<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td>t<sub>w</sub><sup>&perp;</sup> + 2<sup>23.5</sup></td></tr><tr><td>E<sub>7</sub></td><td>SG</td><td>0.7</td><td>JUMP</td><td>t<sub>w</sub><sup>&perp;</sup> + 2<sup>31.25</sup></td></tr><tr><td>A<sub>9</sub></td><td>SG</td><td>0.9</td><td>NATIVE</td><td>t<sub>w</sub><sup>&perp;</sup></td></tr><tr><td>B<sub>9</sub></td><td>SG</td><td>0.9</td><td>JUMP/BACK</td><td>t<sub>w</sub><sup>&perp;</sup> + t<sub>w</sub><sup>&perp;</sup> + 2<sup>10</sup></td></tr><tr><td>C<sub>9</sub></td><td>SG</td><td>0.9</td><td>JUMP/BACK</td><td>t<sub>w</sub><sup>&perp;</sup> + t<sub>w</sub><sup>&perp;</sup> + 2<sup>15.625</sup></td></tr><tr><td>D<sub>9</sub></td><td>SG</td><td>0.9</td><td>JUMP/BACK</td><td>t<sub>w</sub><sup>&perp;</sup> + t<sub>w</sub><sup>&perp;</sup> + 2<sup>23.5</sup></td></tr><tr><td>E<sub>9</sub></td><td>SG</td><td>0.9</td><td>JUMP/BACK</td><td>t<sub>w</sub><sup>&perp;</sup> + t<sub>w</sub><sup>&perp;</sup> + 2<sup>31</sup></td></tr><tr><td>A'<sub>5</sub></td><td>DIM</td><td>0.5</td><td>NATIVE</td><td>76</td></tr><tr><td>B'<sub>5</sub></td><td>DIM</td><td>0.5</td><td>JUMP</td><td>430 + 69</td></tr><tr><td>A'<sub>7</sub></td><td>DIM</td><td>0.7</td><td>NATIVE</td><td>197</td></tr><tr><td>B'<sub>7</sub></td><td>DIM</td><td>0.7</td><td>JUMP</td><td>430 + 165</td></tr><tr><td>A'<sub>9</sub></td><td>DIM</td><td>0.9</td><td>NATIVE</td><td>430</td></tr></tbody></table>
To address this problem, we computed the distri-
bution function $F(\tilde{X})$ for various sphere radius for
two models. One of the models is a spin glass (the
Edwards-Anderson model considered in this work),
the other is the Diluted Ising Model (DIM; this is a
ferromagnetic model where temperature chaos is ab-
sent). As explained in the main text and in **Methods**,
the conditions are carefully matched in both compu-
tations. Hence, we consider that a given $R$ is safe
to use when the DIM distribution function concen-
trates at $X \approx 1$. As Fig. 8 shows, this is the case

for $R = 3a_0$ and $5a_0$. Yet, when $R$ becomes sig-
nificantly larger than the smallest coherence length
$\xi_{\text{micro}} \approx 5.8a_0$, we start finding smaller values of the
correlation parameter $X$ even for the chaos-free DIM.
This suggests that including in the analysis spheres
larger than $\xi_{\text{micro}}$ may produce misleading results (al-
though in all cases the spin-glass distributions have a
much larger weight at small $X$). The good news are
that, for safe sphere radius $R \le \xi_{\text{micro}}$, the distribu-
tion function show a very weak dependence on $R$ (see
Fig. 8).

[1] Q. Zhai, I. Paga, M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pemartin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz Sudupe, D. Navarro, R. L. Orbach, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, D. L. Schlagel, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes, Phys. Rev. Lett. 125, 237202 (2020).

[2] I. Paga, Q. Zhai, M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pemartin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz-Sudupe, D. Navarro, R. L. Orbach, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, D. L. Schlagel, B. Seoane, A. Tarancon, R. Tripiccione, and

D. Yllanes, J. Stat. Mech. 2021, 033301 (2021).

[3] M. Baity-Jesi, R. A. Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, M. Pivanti, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Phys. Rev. B 88, 224416 (2013), arXiv:1310.2910.

[4] F. Belletti, A. Cruz, L. A. Fernandez, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, J. Monforte, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, D. Sciretti, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), J. Stat. Phys. 135, 1121 (2009).
To address this problem, we computed the distri-
bution function $F(\tilde{X})$ for various sphere radius for
two models. One of the models is a spin glass (the
Edwards-Anderson model considered in this work),
the other is the Diluted Ising Model (DIM; this is a
ferromagnetic model where temperature chaos is ab-
sent). As explained in the main text and in **Methods**,
the conditions are carefully matched in both compu-
tations. Hence, we consider that a given $R$ is safe
to use when the DIM distribution function concen-
trates at $X \approx 1$. As Fig. 8 shows, this is the case

for $R = 3a_0$ and $5a_0$. Yet, when $R$ becomes sig-
nificantly larger than the smallest coherence length
$\xi_{\text{micro}} \approx 5.8a_0$, we start finding smaller values of the
correlation parameter $X$ even for the chaos-free DIM.
This suggests that including in the analysis spheres
larger than $\xi_{\text{micro}}$ may produce misleading results (al-
though in all cases the spin-glass distributions have a
much larger weight at small $X$). The good news are
that, for safe sphere radius $R \le \xi_{\text{micro}}$, the distribu-
tion function show a very weak dependence on $R$ (see
Fig. 8).

[1] Q. Zhai, I. Paga, M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pemartin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz Sudupe, D. Navarro, R. L. Orbach, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, D. L. Schlagel, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes, Phys. Rev. Lett. 125, 237202 (2020).

[2] I. Paga, Q. Zhai, M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, I. Gonzalez-Adalid Pemartin, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. Muñoz-Sudupe, D. Navarro, R. L. Orbach, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, D. L. Schlagel, B. Seoane, A. Tarancon, R. Tripiccione, and

D. Yllanes, J. Stat. Mech. 2021, 033301 (2021).

[3] M. Baity-Jesi, R. A. Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iñiguez, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, M. Pivanti, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Phys. Rev. B 88, 224416 (2013), arXiv:1310.2910.

[4] F. Belletti, A. Cruz, L. A. Fernandez, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, J. Monforte, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, D. Sciretti, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), J. Stat. Phys. 135, 1121 (2009).
<figure><img src="image_11.png" /><figcaption>FIG. 5. Fraction of the spheres with a chaotic parameter $X < \tilde{X}$, as obtained for different pairs of protocols and different radius of the studied spherical regions. The compared configurations belong to the temperature jump from $T = 0.9$ to $T = 0.5$. Every panel shows in the title the used pair of configurations according to the notation introduced in Table I.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 6. Fraction of the spheres with a chaotic parameter $X < \tilde{X}$, as obtained for different pairs of protocols and different radius of the studied spherical regions. The compared configurations belong to the temperature jump from $T = 0.9$ to $T = 0.7$. Every panel shows in the title the used pair of configurations according to the notation introduced in Table I.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 5. Fraction of the spheres with a chaotic parameter $X < \tilde{X}$, as obtained for different pairs of protocols and different radius of the studied spherical regions. The compared configurations belong to the temperature jump from $T = 0.9$ to $T = 0.5$. Every panel shows in the title the used pair of configurations according to the notation introduced in Table I.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 6. Fraction of the spheres with a chaotic parameter $X < \tilde{X}$, as obtained for different pairs of protocols and different radius of the studied spherical regions. The compared configurations belong to the temperature jump from $T = 0.9$ to $T = 0.7$. Every panel shows in the title the used pair of configurations according to the notation introduced in Table I.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 7. Distribution function of the chaotic parameter computed in spheres of radius $R = 3$. We simulate a very long simulation experimenting two temperature jumps and we compare the last configurations before the temperature jump from $T = 0.9$ to $T = 0.5$ (i.e. configurations $A_9$ in I) with several configurations along the simulation (i.e. $B_9, C_9, D_9$, and $E_9$ in Table I).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 8. Fraction of the spheres with a chaotic parameter $X < \tilde{X}$, as obtained for different pairs of protocols and different radius of the studied spherical regions. The compared configurations are exactly the same that appear in the top part of Fig. 4 in the main text (for the reader's convenience we give as well in Table I the protocols labelling). We include in the plots a comparison with the Diluted Ising Model (DIM) where no temperature chaos is expected.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 7. Distribution function of the chaotic parameter computed in spheres of radius $R = 3$. We simulate a very long simulation experimenting two temperature jumps and we compare the last configurations before the temperature jump from $T = 0.9$ to $T = 0.5$ (i.e. configurations $A_9$ in I) with several configurations along the simulation (i.e. $B_9, C_9, D_9$, and $E_9$ in Table I).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 8. Fraction of the spheres with a chaotic parameter $X < \tilde{X}$, as obtained for different pairs of protocols and different radius of the studied spherical regions. The compared configurations are exactly the same that appear in the top part of Fig. 4 in the main text (for the reader's convenience we give as well in Table I the protocols labelling). We include in the plots a comparison with the Diluted Ising Model (DIM) where no temperature chaos is expected.</figcaption></figure>