# Multifractality in spin glasses

Marco Baity-Jesi, Enrico Calore, Andrés Cruz, Luis Antonio Fernández, José

Miguel Gil-Narvión, Isidoro González-Adalid Pemartín, Antonio

Gordillo-Guerrero, David Íñiguez, Andrea Maiorano, Enzo Marinari, et al.

► To cite this version:

Marco Baity-Jesi, Enrico Calore, Andrés Cruz, Luis Antonio Fernández, José Miguel Gil-Narvión, et al. Multifractality in spin glasses. Proceedings of the National Academy of Sciences of the United States of America, In press, 121 (2), 10.1073/pnas.2312880120. hal-04344056

HAL Id: hal-04344056

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

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>
# Multifractality in spin glasses

Marco Baity-Jesi, Enrico Calore, Andrés Cruz, Luis Antonio Fernández, José

Miguel Gil-Narvión, Isidoro González-Adalid Pemartín, Antonio

Gordillo-Guerrero, David Íñiguez, Andrea Maiorano, Enzo Marinari, et al.

► To cite this version:

Marco Baity-Jesi, Enrico Calore, Andrés Cruz, Luis Antonio Fernández, José Miguel Gil-Narvión, et al. Multifractality in spin glasses. Proceedings of the National Academy of Sciences of the United States of America, In press, 121 (2), pp.e2312880120. 10.1073/pnas.2312880120. hal-04344056

HAL Id: hal-04344056

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

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
Multifractality in spin glasses

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, Via Giuseppe Saragat, 1, 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,
3100 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

¹²Institute of Nanotechnology, Consiglio Nazionale delle Ricerche
(CNR-NANOTEC), Piazzale Aldo Moro 5, I-00185 Rome, Italy

¹³Dipartimento di Scienze dell'Ambiente e della Prevenzione Università
di Ferrara e INFN Sezione di Ferrara, I-44122 Ferrara, Italy

¹⁴Université Paris-Saclay, CNRS, INRIA Tau team, LISN, 91190, Gif-sur-Yvette, France

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

(Dated: January 23, 2024)

We unveil the multifractal behavior of Ising spin glasses in their low-temperature phase. Using the Janus II custom-built supercomputer, the spin-glass correlation function is studied locally. Dramatic fluctuations are found when pairs of sites at the same distance are compared. The scaling of these fluctuations, as the spin-glass coherence length grows with time, is characterized through the computation of the singularity spectrum and its corresponding Legendre transform. A comparatively small number of site pairs controls the average correlation that governs the response to a magnetic field. We explain how this scenario of dramatic fluctuations (at length scales smaller than the coherence length) can be reconciled with the smooth, self-averaging behavior that has long been considered to describe spin-glass dynamics.

The notion of multifractality [1, 2] refers to situations where many different fractal behaviors coexist within the same system. A major role is played in this context by scale symmetry, see, e.g., [3–5]: in many situations in physics, chemistry and beyond, apparently random objects look the same when the observation scale is changed. The scale change is often quantitatively characterized through a number, the fractal dimension. Multifractals (as opposed to fractals) are systems that need many fractal dimensions to get their scaling properties fully characterized.

Some of the first examples of multifractal behavior appeared in physics, in the contexts of turbulence [6], Anderson localization [7] and diffusion-limited aggregates [8]. A unifying language was soon introduced in a study of chaotic dynamics [9, 10]. The concept has gained popularity as the list of systems exhibiting some form of multifractality has steadily grown. To name only a few, let us recall surface growth [11], human heartbeat

dynamics [12], mating copepods [13, 14], rainfall [15] or the analysis of financial time series [16].

Here we add a (perhaps) surprising member to the list: the off-equilibrium dynamics of spin-glass systems [17, 18]. These disordered magnetic alloys have long been regarded as a paradigmatic toy model for the study of glassiness, optimization, biology, financial markets or social dynamics. It is surprising that such a prominent feature as multifractality has gone unnoticed for such a well-studied model.

The explanation for the above paradox rests on the finite coherence length $\xi(t_w)$ that develops when a spin glass, initially at some very high temperature, is suddenly cooled below the critical temperature $T_c$, and let to relax for a waiting time $t_w$—most experimental work on spin glasses is carried out under non-equilibrium conditions [19]. As $t_w$ increases, glassy domains of growing size $\xi(t_w)$ develop, see Fig. 1. The growth of $\xi(t_w)$ is sluggish for a spin glass, reaching only $\xi \sim 200$ lattice
Multifractality in spin glasses

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, Via Giuseppe Saragat, 1, 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,
3100 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

¹²Institute of Nanotechnology, Consiglio Nazionale delle Ricerche
(CNR-NANOTEC), Piazzale Aldo Moro 5, I-00185 Rome, Italy

¹³Dipartimento di Scienze dell'Ambiente e della Prevenzione Università
di Ferrara e INFN Sezione di Ferrara, I-44122 Ferrara, Italy

¹⁴Université Paris-Saclay, CNRS, INRIA Tau team, LISN, 91190, Gif-sur-Yvette, France

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

(Dated: January 23, 2024)

We unveil the multifractal behavior of Ising spin glasses in their low-temperature phase. Using the Janus II custom-built supercomputer, the spin-glass correlation function is studied locally. Dramatic fluctuations are found when pairs of sites at the same distance are compared. The scaling of these fluctuations, as the spin-glass coherence length grows with time, is characterized through the computation of the singularity spectrum and its corresponding Legendre transform. A comparatively small number of site pairs controls the average correlation that governs the response to a magnetic field. We explain how this scenario of dramatic fluctuations (at length scales smaller than the coherence length) can be reconciled with the smooth, self-averaging behavior that has long been considered to describe spin-glass dynamics.

The notion of multifractality [1, 2] refers to situations where many different fractal behaviors coexist within the same system. A major role is played in this context by scale symmetry, see, e.g., [3–5]: in many situations in physics, chemistry and beyond, apparently random objects look the same when the observation scale is changed. The scale change is often quantitatively characterized through a number, the fractal dimension. Multifractals (as opposed to fractals) are systems that need many fractal dimensions to get their scaling properties fully characterized.

Some of the first examples of multifractal behavior appeared in physics, in the contexts of turbulence [6], Anderson localization [7] and diffusion-limited aggregates [8]. A unifying language was soon introduced in a study of chaotic dynamics [9, 10]. The concept has gained popularity as the list of systems exhibiting some form of multifractality has steadily grown. To name only a few, let us recall surface growth [11], human heartbeat

dynamics [12], mating copepods [13, 14], rainfall [15] or the analysis of financial time series [16].

Here we add a (perhaps) surprising member to the list: the off-equilibrium dynamics of spin-glass systems [17, 18]. These disordered magnetic alloys have long been regarded as a paradigmatic toy model for the study of glassiness, optimization, biology, financial markets or social dynamics. It is surprising that such a prominent feature as multifractality has gone unnoticed for such a well-studied model.

The explanation for the above paradox rests on the finite coherence length $\xi(t_w)$ that develops when a spin glass, initially at some very high temperature, is suddenly cooled below the critical temperature $T_c$, and let to relax for a waiting time $t_w$—most experimental work on spin glasses is carried out under non-equilibrium conditions [19]. As $t_w$ increases, glassy domains of growing size $\xi(t_w)$ develop, see Fig. 1. The growth of $\xi(t_w)$ is sluggish for a spin glass, reaching only $\xi \sim 200$ lattice
spacings for $t_w \sim 1$ hour [20, 21]. Now, when one measures the magnetic response to an external field, which is the main experimental probe of spin-glass dynamics, an average over the whole sample is carried out. Since the sample is effectively composed of many independent domains of linear size $\sim \xi(t_w)$, the central limit theorem eliminates from the average response the large fluctuations that could ultimately cause multifractal behavior. With few exceptions (see below), most numerical work has emphasized the space-averaged correlation function in Fig. 1. Besides, see **Methods**, studying correlations without spatial averages is very demanding computationally.

It follows from the above considerations that multifractal behavior in spin glasses should be investigated in large statistical deviations that occur at a length scale smaller than (or comparable to) $\xi(t_w)$, definitively not the standard framework either for experiments[see 21–23, for instance] or for simulations [24–26]. There is, however, an important exception. Recently, progress has been achieved [27] in the theoretical interpretation of the experimental rejuvenation and memory effects in spin glasses [28]. Crucial for this achievement was the study of temperature chaos in the off-equilibrium dynamics at the $\xi(t_w)$ length scale [29], through numerical simulations using the Janus II dedicated supercomputer [30]. As we shall show below, the consideration of fluctuations at the $\xi(t_w)$ length scale still holds surprises.

Specifically, we shall consider the spin-glass correlation function, see **Methods** for definitions. The space-averaged correlation function is a well-known quantity and the basis for the computation of $\xi(t_w)$ explained in Fig. 1. We shall depart from the standard approach, however, by avoiding the spatial average. We shall compute the correlation function for a pair of sites $\mathbf{x}$ and $\mathbf{y}$, and consider the statistical fluctuations induced by varying $\mathbf{x}$ while fixing $\mathbf{r} = \mathbf{y} - \mathbf{x}$.

The reader may argue that it is difficult to find large statistical fluctuations in a mathematical object bounded between 0 and 1, such as the spin-glass correlation function. A moment of thought will reveal that large fluctuations are only possible if the average of the correlation function goes to zero as $\xi(t_w)$ grows, so that the correlation function at a given site can get large if measured in units of the averaged correlation.

Indeed, spin glasses are peculiar among systems with domain-growth off-equilibrium dynamics. Fig. 2 compares the space-averaged correlation function at distance $r = \xi(t_w)$ for two Ising systems in space dimension $D=3$, the link-diluted ferromagnet and the spin glass. In the ferromagnet, the correlation function goes to a constant value (the squared spontaneous magnetization) as $\xi$ grows. Hence, large deviations and multifractality are possible for the ferromagnet only at $T_c$, where the spontaneous magnetization vanishes [32]. In the spin glass, instead, the correlation function scales as $1/r^\theta$ for

<figure><img src="image_2.png" /><figcaption>Figure 1. The correlation function, (10), as computed for the three-dimensional Ising diluted ferromagnet (DIL) and for the Ising spin glass (EA), versus distance $r$. Data were obtained in systems of linear size $L=160$ with coherence length $\xi(t_w) = 20$ (dashed vertical line) at temperature $T=0.9$—recall that $T_c \approx 1.1$ for EA [31]. As explained in **Methods**, the coherence length is computed from the integral $I_2 = \int_0^\infty r^2 C_4^{\text{av}}(r) dr$ (the integrand is shown in the Supporting Information).</figcaption></figure>

distances up to $r \sim \xi$ [33]. Hence, unlike the diluted ferromagnet, the spin glass can accommodate large fluctuations for all $T < T_c$. This is why here we decide to focus on the spin glass.

## RESULTS

The first indication of large deviations in the statistics of the spin-glass correlation function $C_4$ is shown in Fig. 3, where we select the distance $r = \xi(t_w)$. The ratio of the second moment of $C_4$, $\overline{C_4^2}$, to the first moment squared, $\overline{C_4^2}$, nicely follows a power law as a function of $\overline{C_4}$ [this type of analysis was pioneered in 34]. If continued to $\overline{C_4} \to 0$ (i.e., as $\xi(t_w)$ grows, see Fig. 2), this power law implies that the orders of magnitude of $\overline{C_4^2}$ and $\overline{C_4^2}$ differ in the large-$\xi(t_w)$ scaling limit. This behavior is not reminiscent of a monofractal, which in the scaling limit is characterized by a single quantity (say, $\overline{C_4}$).

We also note from Fig. 3 that all our data with $T < T_c$ follow the same scaling curve, which slightly differs from its counterpart at the critical point. This is not completely unexpected, because the $\epsilon$-expansion tells us that the average $C_4$ at $T_c$ decays as a power law with distance with an exponent [35] that is twice as large as the exponent for $T < T_c$ [36]. In fact, we lack an explanation for the similarity of the two exponents that can be observed in Fig. 3. From now on, our analysis will focus on our data at $T=0.9$, namely the temperature in the spin-glass phase where we are able to reach the largest $\xi(t_w)$.

A picture of the physical situation is presented in Fig. 4. We may expect a different behavior for the aver-
spacings for $t_w \sim 1$ hour [20, 21]. Now, when one measures the magnetic response to an external field, which is the main experimental probe of spin-glass dynamics, an average over the whole sample is carried out. Since the sample is effectively composed of many independent domains of linear size $\sim \xi(t_w)$, the central limit theorem eliminates from the average response the large fluctuations that could ultimately cause multifractal behavior. With few exceptions (see below), most numerical work has emphasized the space-averaged correlation function in Fig. 1. Besides, see **Methods**, studying correlations without spatial averages is very demanding computationally.

It follows from the above considerations that multifractal behavior in spin glasses should be investigated in large statistical deviations that occur at a length scale smaller than (or comparable to) $\xi(t_w)$, definitively not the standard framework either for experiments[see 21–23, for instance] or for simulations [24–26]. There is, however, an important exception. Recently, progress has been achieved [27] in the theoretical interpretation of the experimental rejuvenation and memory effects in spin glasses [28]. Crucial for this achievement was the study of temperature chaos in the off-equilibrium dynamics at the $\xi(t_w)$ length scale [29], through numerical simulations using the Janus II dedicated supercomputer [30]. As we shall show below, the consideration of fluctuations at the $\xi(t_w)$ length scale still holds surprises.

Specifically, we shall consider the spin-glass correlation function, see **Methods** for definitions. The space-averaged correlation function is a well-known quantity and the basis for the computation of $\xi(t_w)$ explained in Fig. 1. We shall depart from the standard approach, however, by avoiding the spatial average. We shall compute the correlation function for a pair of sites $\mathbf{x}$ and $\mathbf{y}$, and consider the statistical fluctuations induced by varying $\mathbf{x}$ while fixing $\mathbf{r} = \mathbf{y} - \mathbf{x}$.

The reader may argue that it is difficult to find large statistical fluctuations in a mathematical object bounded between 0 and 1, such as the spin-glass correlation function. A moment of thought will reveal that large fluctuations are only possible if the average of the correlation function goes to zero as $\xi(t_w)$ grows, so that the correlation function at a given site can get large if measured in units of the averaged correlation.

Indeed, spin glasses are peculiar among systems with domain-growth off-equilibrium dynamics. Fig. 2 compares the space-averaged correlation function at distance $r = \xi(t_w)$ for two Ising systems in space dimension $D=3$, the link-diluted ferromagnet and the spin glass. In the ferromagnet, the correlation function goes to a constant value (the squared spontaneous magnetization) as $\xi$ grows. Hence, large deviations and multifractality are possible for the ferromagnet only at $T_c$, where the spontaneous magnetization vanishes [32]. In the spin glass, instead, the correlation function scales as $1/r^\theta$ for

<figure><img src="image_2.png" /><figcaption>Figure 1. The correlation function, (10), as computed for the three-dimensional Ising diluted ferromagnet (DIL) and for the Ising spin glass (EA), versus distance $r$. Data were obtained in systems of linear size $L=160$ with coherence length $\xi(t_w) = 20$ (dashed vertical line) at temperature $T=0.9$—recall that $T_c \approx 1.1$ for EA [31]. As explained in **Methods**, the coherence length is computed from the integral $I_2 = \int_0^\infty r^2 C_4^{\text{av}}(r) dr$ (the integrand is shown in the Supporting Information).</figcaption></figure>

distances up to $r \sim \xi$ [33]. Hence, unlike the diluted ferromagnet, the spin glass can accommodate large fluctuations for all $T < T_c$. This is why here we decide to focus on the spin glass.

## RESULTS

The first indication of large deviations in the statistics of the spin-glass correlation function $C_4$ is shown in Fig. 3, where we select the distance $r = \xi(t_w)$. The ratio of the second moment of $C_4$, $\overline{C_4^2}$, to the first moment squared, $\overline{C_4^2}$, nicely follows a power law as a function of $\overline{C_4}$ [this type of analysis was pioneered in 34]. If continued to $\overline{C_4} \to 0$ (i.e., as $\xi(t_w)$ grows, see Fig. 2), this power law implies that the orders of magnitude of $\overline{C_4^2}$ and $\overline{C_4^2}$ differ in the large-$\xi(t_w)$ scaling limit. This behavior is not reminiscent of a monofractal, which in the scaling limit is characterized by a single quantity (say, $\overline{C_4}$).

We also note from Fig. 3 that all our data with $T < T_c$ follow the same scaling curve, which slightly differs from its counterpart at the critical point. This is not completely unexpected, because the $\epsilon$-expansion tells us that the average $C_4$ at $T_c$ decays as a power law with distance with an exponent [35] that is twice as large as the exponent for $T < T_c$ [36]. In fact, we lack an explanation for the similarity of the two exponents that can be observed in Fig. 3. From now on, our analysis will focus on our data at $T=0.9$, namely the temperature in the spin-glass phase where we are able to reach the largest $\xi(t_w)$.

A picture of the physical situation is presented in Fig. 4. We may expect a different behavior for the aver-
<figure><img src="image_3.png" /><figcaption>Figure 2. Correlation function $C_4^{\text{av}}(r = \xi(t_w))$, see (10) versus the coherence length $\xi(t_w)$, as computed for DIL (top) and for EA (bottom) at temperatures $T$, $\bar{T} = 0.9, 0.8$ and $0.7$ (see Methods for a complete definition of $\bar{T}$). Error bars are smaller than the point size. The dashed line is our fit to (2), with $q = 1$, for EA at $T = 0.9$ (to avoid scaling corrections, we fit in the range $\xi(t_w) \in [10, 20]$, see SI for further information). Note that, while the DIL $C_4^{\text{av}}(r = \xi(t_w))$ tends to a T-dependent positive limit for large coherence length (which excludes multiscaling at $T < T_c$), the spin-glass correlation functions steadily decrease with $\xi(t_w)$.</figcaption></figure>

age and the local correlation function when distances up to $r \sim \xi(t_w)$ are considered ($\theta(T = 0.9) \approx 0.4$ [26]) [37]:

$$C_4^{\text{av}}(r, t_w) \sim \frac{1}{r^{\theta}}, \quad C_4(\mathbf{x}, \mathbf{x} + \mathbf{r}; t_w) \sim \frac{1}{r^{\theta} M(\mathbf{x}, \mathbf{r}, t_w)}. \quad (1)$$

As the reader can check from Fig. 4, the order-of-magnitude modulating factor $M(\mathbf{x}, \mathbf{r}, t_w)$ varies by a factor of 16, which indicates that there are site pairs $(\mathbf{x}, \mathbf{x} + \mathbf{r})$ a lot more—or a lot less—correlated than the average. In fact, see Fig. 5, the median correlation function at distance $r = \xi(t_w)$, scales as $[C_4^{\text{av}}]^a$, with $a \approx 1.5$. In other words, the typical correlation function is a lot smaller than the average value.

In order to make the above qualitative description quantitative, we consider the moments of the probability distribution of $C_4$ at distance $r = \xi(t_w)$. The q-th moment turns out to follow a scaling law

$$\bar{C}_4^q \sim \frac{A_q}{\xi \tau(q)}. \quad (2)$$

<figure><img src="image_4.png" /><figcaption>Figure 3. Ratio of the second moment of the spin-glass correlation function $C_4$ computed at $r = \xi(t_w)$, $\bar{C}_4^2$, to the squared first moment, $\bar{C}_4^2$, as a function of $\bar{C}_4$. We show the data for all temperatures considered in this work. $\bar{C}_4$ tends to zero as the coherence length $\xi(t_w)$ gets large, recall Fig. 2. Note that $\bar{C}_4^2/\bar{C}_4^2$ scales with $\bar{C}_4$ as a power law, which indicates that in the scaling limit (i.e., $\xi(t_w) \to \infty$ or $\bar{C}_4 \to 0$) the order of magnitude of $\bar{C}_4^2$ is larger than the one of $\bar{C}_4^2$. Data in the glassy phase, $T < T_c$, roughly follow the same scaling curve. At the critical point there is still a power type relation with a slightly different exponent. Error bars are smaller than the points size. The same data are shown as a function of $\xi(t_w)$ in the SI.</figcaption></figure>

Fig. 6 shows the $\tau(q)$ function, which significantly differs from the monofractal behavior $\tau^{\text{mono}}(q) = q \tau^{\text{mono}}(1)$. It is this departure from linear behavior that justifies using the term multifractal to describe spin-glass dynamics [see, e.g., 9].

For large moments, $\tau(q)$ seems to grow as $\log q$ (see the inset in Fig. 6). The origin of this logarithmic growth seems to be in the behavior of the probability distribution function $P(C_4)$ near $C_4 = 1$. As shown in the Supporting Information (SI), the numerical data are consistent with $P(C_4) \propto (1 - C_4)^{B(\xi)}$ for $C_4$ close to 1, with an exponent that grows as $B(\xi(t_w)) \sim \log \xi(t_w)$. This behavior of the correlation function would explain the observed logarithmic growth of $\tau(q)$. However, just to be on the safe side, we have tried two different functional forms to fit the numerical data in Fig. 6:

$$\tau_1(q) = mq \frac{1 + c_1 q}{1 + c_2 q}, \quad \tau_2(q) = mq \frac{1 + d_1 q \log q}{(1 + d_2 q)^2}. \quad (3)$$

Both $\tau_1$ and $\tau_2$ have the same derivative $m$ at $q=0$. We do not treat $m$ as a fitting parameter. Rather we take it from the scaling of the median of the distribution $P(C_4)$ with $\xi$ (see SI). Although both $\tau_1(q)$ and $\tau_2(q)$ make an excellent job at fitting our data (see again SI), only $\tau_2(q)$ displays the logarithmic growth with $q$, at large $q$, that we find more plausible.
<figure><img src="image_3.png" /><figcaption>Figure 2. Correlation function $C_4^{\text{av}}(r = \xi(t_w))$, see (10) versus the coherence length $\xi(t_w)$, as computed for DIL (top) and for EA (bottom) at temperatures $T$, $\bar{T} = 0.9, 0.8$ and $0.7$ (see Methods for a complete definition of $\bar{T}$). Error bars are smaller than the point size. The dashed line is our fit to (2), with $q = 1$, for EA at $T = 0.9$ (to avoid scaling corrections, we fit in the range $\xi(t_w) \in [10, 20]$, see SI for further information). Note that, while the DIL $C_4^{\text{av}}(r = \xi(t_w))$ tends to a $T$-dependent positive limit for large coherence length (which excludes multiscaling at $T < T_c$), the spin-glass correlation functions steadily decrease with $\xi(t_w)$.</figcaption></figure>

age and the local correlation function when distances up to $r \sim \xi(t_w)$ are considered ($\theta(T = 0.9) \approx 0.4$ [26]) [37]:

$$C_4^{\text{av}}(r, t_w) \sim \frac{1}{r^{\theta}}, \quad C_4(\mathbf{x}, \mathbf{x} + \mathbf{r}; t_w) \sim \frac{1}{r^{\theta} M(\mathbf{x}, \mathbf{r}, t_w)}. \quad (1)$$

As the reader can check from Fig. 4, the order-of-magnitude modulating factor $M(\mathbf{x}, \mathbf{r}, t_w)$ varies by a factor of 16, which indicates that there are site pairs $(\mathbf{x}, \mathbf{x} + \mathbf{r})$ a lot more—or a lot less—correlated than the average. In fact, see Fig. 5, the median correlation function at distance $r = \xi(t_w)$, scales as $[C_4^{\text{av}}]^a$, with $a \approx 1.5$. In other words, the typical correlation function is a lot smaller than the average value.

In order to make the above qualitative description quantitative, we consider the moments of the probability distribution of $C_4$ at distance $r = \xi(t_w)$. The $q$-th moment turns out to follow a scaling law

$$\bar{C}_4^q \sim \frac{A_q}{\xi \tau(q)}. \quad (2)$$

<figure><img src="image_4.png" /><figcaption>Figure 3. Ratio of the second moment of the spin-glass correlation function $C_4$ computed at $r = \xi(t_w)$, $\bar{C}_4^2$, to the squared first moment, $\bar{C}_4^2$, as a function of $\bar{C}_4$. We show the data for all temperatures considered in this work. $\bar{C}_4$ tends to zero as the coherence length $\xi(t_w)$ gets large, recall Fig. 2. Note that $\bar{C}_4^2/\bar{C}_4^2$ scales with $\bar{C}_4$ as a power law, which indicates that in the scaling limit (i.e., $\xi(t_w) \to \infty$ or $\bar{C}_4 \to 0$) the order of magnitude of $\bar{C}_4^2$ is larger than the one of $\bar{C}_4^2$. Data in the glassy phase, $T < T_c$, roughly follow the same scaling curve. At the critical point there is still a power type relation with a slightly different exponent. Error bars are smaller than the points size. The same data are shown as a function of $\xi(t_w)$ in the SI.</figcaption></figure>

Fig. 6 shows the $\tau(q)$ function, which significantly differs from the monofractal behavior $\tau^{\text{mono}}(q) = q \tau^{\text{mono}}(1)$. It is this departure from linear behavior that justifies using the term multifractal to describe spin-glass dynamics [see, e.g., 9].

For large moments, $\tau(q)$ seems to grow as $\log q$ (see the inset in Fig. 6). The origin of this logarithmic growth seems to be in the behavior of the probability distribution function $P(C_4)$ near $C_4 = 1$. As shown in the Supporting Information (SI), the numerical data are consistent with $P(C_4) \propto (1 - C_4)^{B(\xi)}$ for $C_4$ close to 1, with an exponent that grows as $B(\xi(t_w)) \sim \log \xi(t_w)$. This behavior of the correlation function would explain the observed logarithmic growth of $\tau(q)$. However, just to be on the safe side, we have tried two different functional forms to fit the numerical data in Fig. 6:

$$\tau_1(q) = mq \frac{1 + c_1 q}{1 + c_2 q}, \quad \tau_2(q) = mq \frac{1 + d_1 q \log q}{(1 + d_2 q)^2}. \quad (3)$$

Both $\tau_1$ and $\tau_2$ have the same derivative $m$ at $q=0$. We do not treat $m$ as a fitting parameter. Rather we take it from the scaling of the median of the distribution $P(C_4)$ with $\xi$ (see SI). Although both $\tau_1(q)$ and $\tau_2(q)$ make an excellent job at fitting our data (see again SI), only $\tau_2(q)$ displays the logarithmic growth with $q$, at large $q$, that we find more plausible.
<figure><img src="image_5.png" /><figcaption>Figure 4. Grayscale representation of the order-of-magnitude modulating factor $M(\mathbf{x}, \mathbf{r}, t_w)$, see (1), computed for site $\mathbf{x} = (64, 64, 64)$ of a sample with coherence length $\xi(t_w) = 20$, at $T = 0.9$, with an $N_R = 512$ estimator (see Methods). We show results for displacement vectors $\mathbf{r} = (r_x, r_y, r_z)$ in a cube $-40 \le r_x, r_y, r_z \le 40$. The top-left panel depicts the three visible faces of the cube, while the other three panels show sections at $r_z = -20, 0, 20$, respectively. Our color code is darker the smaller $M(\mathbf{x}, \mathbf{r}, t_w)$ (hence, the more slowly correlations decay with distance). For ease of representation, we have chosen a color code linear between the minimal value of $M(\mathbf{x}, \mathbf{r}, t_w)$ and 2.5. Displacements $\mathbf{r}$ with $M(\mathbf{x}, \mathbf{r}, t_w) > 2.5$ are depicted as if $M(\mathbf{x}, \mathbf{r}, t_w) = 2.5$. See SI for more examples of this modulating factor.</figcaption></figure>

DISCUSSION

Following Ref. [9], we shall discuss our results in terms of a different stochastic variable, $\alpha = \log C_4(r = \xi(t_w))/\log[1/\xi(t_w)]$, so that (we drop the argument in $\xi$ for the sake of shortness)

$$C_4 = \frac{1}{\xi\alpha}, \quad P(C_4) \frac{dC_4}{d\alpha} \sim \xi^{f(\alpha)}. \quad (4)$$

(4) defines the large-deviations function $f(\alpha)$. Then, we find for the moments of $C_4$

$$\overline{C_4^q} = \int_0^1 dC_4 P(C_4) C_4^q \sim \int_0^\infty d\alpha e^{\log(\xi)[f(\alpha) - q\alpha]}. \quad (5)$$

For large $\xi$, the above integral is dominated by the maximum of $[f(\alpha) - q\alpha]$ at some value $\alpha = \alpha^*$:

$$\overline{C_4^q} \sim \frac{1}{\xi - [f(\alpha^*) - q\alpha^*]} \quad (6)$$

Comparing with (2), we realize that $f(\alpha)$ is just (minus) the Legendre transform of the singularity spectrum $\tau(q)$:

$$f(\alpha) = - \max_q [\tau(q) - q\alpha]. \quad (7)$$

We show $f(\alpha)$ in Fig. 7, as computed from our fitting ansätze $\tau_1(q)$ and $\tau_2(q)$, in (3). In the range of Fig. 6 — since $\alpha(q) = \tau'(q)$— the results from the two ansätze can hardly be distinguished. The two, however, differ in that
<figure><img src="image_5.png" /><figcaption>Figure 4. Grayscale representation of the order-of-magnitude modulating factor $M(\mathbf{x}, \mathbf{r}, t_w)$, see (1), computed for site $\mathbf{x} = (64, 64, 64)$ of a sample with coherence length $\xi(t_w) = 20$, at $T = 0.9$, with an $N_R = 512$ estimator (see Methods). We show results for displacement vectors $\mathbf{r} = (r_x, r_y, r_z)$ in a cube $-40 \le r_x, r_y, r_z \le 40$. The top-left panel depicts the three visible faces of the cube, while the other three panels show sections at $r_z = -20, 0, 20$, respectively. Our color code is darker the smaller $M(\mathbf{x}, \mathbf{r}, t_w)$ (hence, the more slowly correlations decay with distance). For ease of representation, we have chosen a color code linear between the minimal value of $M(\mathbf{x}, \mathbf{r}, t_w)$ and 2.5. Displacements $\mathbf{r}$ with $M(\mathbf{x}, \mathbf{r}, t_w) > 2.5$ are depicted as if $M(\mathbf{x}, \mathbf{r}, t_w) = 2.5$. See SI for more examples of this modulating factor.</figcaption></figure>

DISCUSSION

Following Ref. [9], we shall discuss our results in terms of a different stochastic variable, $\alpha = \log C_4(r = \xi(t_w))/\log[1/\xi(t_w)]$, so that (we drop the argument in $\xi$ for the sake of shortness)

$$C_4 = \frac{1}{\xi\alpha}, \quad P(C_4) \frac{dC_4}{d\alpha} \sim \xi^{f(\alpha)}. \quad (4)$$

(4) defines the large-deviations function $f(\alpha)$. Then, we find for the moments of $C_4$

$$\overline{C_4^q} = \int_0^1 dC_4 P(C_4) C_4^q \sim \int_0^\infty d\alpha e^{\log(\xi)[f(\alpha) - q\alpha]}. \quad (5)$$

For large $\xi$, the above integral is dominated by the maximum of $[f(\alpha) - q\alpha]$ at some value $\alpha = \alpha^*$:

$$\overline{C_4^q} \sim \frac{1}{\xi - [f(\alpha^*) - q\alpha^*]} \quad (6)$$

Comparing with (2), we realize that $f(\alpha)$ is just (minus) the Legendre transform of the singularity spectrum $\tau(q)$:

$$f(\alpha) = - \max_q [\tau(q) - q\alpha]. \quad (7)$$

We show $f(\alpha)$ in Fig. 7, as computed from our fitting ansätze $\tau_1(q)$ and $\tau_2(q)$, in (3). In the range of Fig. 6 — since $\alpha(q) = \tau'(q)$— the results from the two ansätze can hardly be distinguished. The two, however, differ in that
the range of $\alpha$ for $\tau_2(q)$ goes all the way down to $\alpha = 0$ (because $\alpha_2(q) = d\tau_2/dq \sim 1/q$). Indeed, if $\tau(q)$ goes as $\log(q)$ for large $q$, then the large-deviations function goes as $f(\alpha \to 0) \sim \log(\alpha)$.

Let us recapitulate: the probability of finding a site $\mathbf{x}$ with $C_4(\mathbf{x}, \mathbf{x} + \mathbf{r})$ scaling as $1/r^\alpha$ for $r \sim \xi$ goes in the scaling limit as $\xi^{f(\alpha)}$. There are, hence, a lot more sites displaying the median scaling exponent $\alpha \approx 0.65$ than there are for the average scaling $\alpha \approx 0.4$ (because $f(0.65) > f(0.4)$, recall Fig. 5). The larger $\xi(t_w)$ grows, the more pronounced this difference is. Thus, the expression “silent majority” [38] could be aptly employed to describe spin-glass dynamics: the central limit theorem ensures that it is the (somewhat exceptional) average value the one that can be measured on length scales larger than $\xi(t_w)$ (hence, in experiments). The experimental-scale dynamics is, however, not completely blind to these short-scale fluctuations. Indeed, temperature chaos [29] —and, hence, rejuvenation [27], which is certainly experimentally observable [see, e.g., 28]— is ruled by statistical fluctuations at the scale of $r$ smaller than, or similar to, $\xi(t_w)$.

Our data show that varying $T$ simply changes $\tau(q)$ by an essentially constant factor $[e.g., \tau(q, T_c) \approx 1.5\tau(q, T = 0.9), \text{see SI}]$. Furthermore, Fig. 3 makes us confident that, taking $\bar{C}_4$ as scaling variable instead of $\xi(t_w)$, the overall picture is essentially temperature independent for $T < T_c$.

Whether or not multifractal behavior is also present in equilibrium correlation functions in the spin-glass phase stands out as an interesting open question. Statics-dynamics equivalence [25, 39–41] suggests that the answer will be positive.

As a final remark, let us stress that ongoing efforts to build a mathematically rigorous theory of non-equilibrium spin-glass dynamics through the concept of the maturation metastate (see 42 and references therein) should take into account the extreme spatial heterogeneity unveiled in this work.

# METHODS

## Model and simulations

We focus on the Edwards-Anderson model (EA) in a simple cubic lattice with linear size $L = 160$ and periodic boundary conditions. Our $S_\mathbf{x} = \pm 1$ spin, placed at the lattice sites, interact with their nearest neighbors through the Hamiltonian:

$$ \mathcal{H} = - \sum_{\langle \mathbf{x}, \mathbf{y} \rangle} J_{\mathbf{x}, \mathbf{y}} S_{\mathbf{x}} S_{\mathbf{y}}. \qquad (8) $$

The coupling constants $J_{\mathbf{x},\mathbf{y}}$ are independent random variables ($J_{\mathbf{x},\mathbf{y}} = \pm 1$ with equal probability), fixed once

and for all at the beginning of the simulation (this is named quenched disorder). A realization of the couplings is called a *sample*. We shall use 16 samples in this work. In general, errors will be computed with a jack-knife method over the samples [see, for instance, 43, 44]. We have also considered a diluted Ising model (see below), as a baseline model displaying domain-growth off-equilibrium dynamics.

We have simulated the model in (8) through a Metropolis dynamics on the Janus II supercomputer [30]. Our time unit is a full-lattice sweep, roughly equivalent to a picosecond of physical time [17]. The critical temperature for this model is $T_c = 1.1019(29)$ [31].

For each sample, we have simulated $N_R = 512$ statistically independent system copies or replicas. We denote by $\langle \cdots \rangle$ the average over thermal noise for one sample (as explained below, we obtain unbiased estimators of the thermal expectation values $\langle \cdots \rangle$ by averaging over the replicas). The subsequent average over samples is denoted by an overline ($\overline{\langle \cdots \rangle}$).

The main quantity of interest is the correlation function

$$ C_4(\mathbf{x}, \mathbf{y}, t_w) = \langle S_{\mathbf{x}}(t_w) S_{\mathbf{y}}(t_w) \rangle^2. \qquad (9) $$

Note that, for a given sample and $(\mathbf{x}, \mathbf{y}, t_w)$, $C_4(\mathbf{x}, \mathbf{y}, t_w)$ is not a stochastic variable. However, it is a stochastic variable if we regard the variations induced by the choice of couplings $J_{\mathbf{x},\mathbf{y}}$ and over the considered sites $(\mathbf{x}, \mathbf{y}, t_w)$. We term these stochastic variables $C_4$, without arguments.

As explained in the next paragraph, although $C_4(\mathbf{x}, \mathbf{y}, t_w)$ cannot be computed with a finite number of replicas, unbiased estimators of its moments can be computed. In particular, previous work has mostly focused on the average correlation function

$$ C_4^{\text{av}}(\mathbf{r}, t_w) = \frac{1}{L^3} \sum_{\mathbf{x}} \overline{C_4(\mathbf{x}, \mathbf{y} = \mathbf{x} + \mathbf{r}, t_w)}. \qquad (10) $$

Cubic symmetry, present in averages over the samples, allows us to average over the three equivalent displacements $\mathbf{r} = (r, 0, 0)$ and permutations. We shall use the shorthand $C_4^{\text{av}}(r, t_w)$ to indicate this average over the three equivalent $\mathbf{r}$. To compute the coherence length $\xi(t_w)$ we follow [26, 45, 46] and compute the integrals

$$ I_n(t_w) = \int_0^\infty r^n C_4^{\text{av}}(r, t_w) dr. \qquad (11) $$

Then, $\xi(t_w) = I_2(t_w)/I_1(t_w)$.

As stated above, we have simulated, as a null experiment, the link-diluted Ising model (DIL). The only difference with the Hamiltonian in (8) is in the choice of the couplings: $J_{\mathbf{x},\mathbf{y}} = 1$ (with 70% probability) or $J_{\mathbf{x},\mathbf{y}} = 0$ (with 30% probability). Since all couplings are positive or zero, this is a ferromagnetic system without frustration.
the range of $\alpha$ for $\tau_2(q)$ goes all the way down to $\alpha = 0$ (because $\alpha_2(q) = d\tau_2/dq \sim 1/q$). Indeed, if $\tau(q)$ goes as $\log(q)$ for large $q$, then the large-deviations function goes as $f(\alpha \to 0) \sim \log(\alpha)$.

Let us recapitulate: the probability of finding a site $\mathbf{x}$ with $C_4(\mathbf{x}, \mathbf{x} + \mathbf{r})$ scaling as $1/r^\alpha$ for $r \sim \xi$ goes in the scaling limit as $\xi^{f(\alpha)}$. There are, hence, a lot more sites displaying the median scaling exponent $\alpha \approx 0.65$ than there are for the average scaling $\alpha \approx 0.4$ (because $f(0.65) > f(0.4)$, recall Fig. 5). The larger $\xi(t_w)$ grows, the more pronounced this difference is. Thus, the expression “silent majority” [38] could be aptly employed to describe spin-glass dynamics: the central limit theorem ensures that it is the (somewhat exceptional) average value the one that can be measured on length scales larger than $\xi(t_w)$ (hence, in experiments). The experimental-scale dynamics is, however, not completely blind to these short-scale fluctuations. Indeed, temperature chaos [29] —and, hence, rejuvenation [27], which is certainly experimentally observable [see, e.g., 28]— is ruled by statistical fluctuations at the scale of $r$ smaller than, or similar to, $\xi(t_w)$.

Our data show that varying $T$ simply changes $\tau(q)$ by an essentially constant factor $[e.g., \tau(q, T_c) \approx 1.5\tau(q, T = 0.9), \text{see SI}]$. Furthermore, Fig. 3 makes us confident that, taking $\bar{C}_4$ as scaling variable instead of $\xi(t_w)$, the overall picture is essentially temperature independent for $T < T_c$.

Whether or not multifractal behavior is also present in equilibrium correlation functions in the spin-glass phase stands out as an interesting open question. Statics-dynamics equivalence [25, 39–41] suggests that the answer will be positive.

As a final remark, let us stress that ongoing efforts to build a mathematically rigorous theory of non-equilibrium spin-glass dynamics through the concept of the maturation metastate (see 42 and references therein) should take into account the extreme spatial heterogeneity unveiled in this work.

# METHODS

## Model and simulations

We focus on the Edwards-Anderson model (EA) in a simple cubic lattice with linear size $L = 160$ and periodic boundary conditions. Our $S_\mathbf{x} = \pm 1$ spin, placed at the lattice sites, interact with their nearest neighbors through the Hamiltonian:

$$ \mathcal{H} = - \sum_{\langle \mathbf{x}, \mathbf{y} \rangle} J_{\mathbf{x}, \mathbf{y}} S_{\mathbf{x}} S_{\mathbf{y}}. \qquad (8) $$

The coupling constants $J_{\mathbf{x}, \mathbf{y}}$ are independent random variables ($J_{\mathbf{x}, \mathbf{y}} = \pm 1$ with equal probability), fixed once

and for all at the beginning of the simulation (this is named quenched disorder). A realization of the couplings is called a *sample*. We shall use 16 samples in this work. In general, errors will be computed with a jack-knife method over the samples [see, for instance, 43, 44]. We have also considered a diluted Ising model (see below), as a baseline model displaying domain-growth off-equilibrium dynamics.

We have simulated the model in (8) through a Metropolis dynamics on the Janus II supercomputer [30]. Our time unit is a full-lattice sweep, roughly equivalent to a picosecond of physical time [17]. The critical temperature for this model is $T_c = 1.1019(29)$ [31].

For each sample, we have simulated $N_R = 512$ statistically independent system copies or replicas. We denote by $\langle \cdots \rangle$ the average over thermal noise for one sample (as explained below, we obtain unbiased estimators of the thermal expectation values $\langle \cdots \rangle$ by averaging over the replicas). The subsequent average over samples is denoted by an overline ($\overline{\langle \cdots \rangle}$).

The main quantity of interest is the correlation function

$$ C_4(\mathbf{x}, \mathbf{y}, t_w) = \langle S_{\mathbf{x}}(t_w) S_{\mathbf{y}}(t_w) \rangle^2. \qquad (9) $$

Note that, for a given sample and $(\mathbf{x}, \mathbf{y}, t_w)$, $C_4(\mathbf{x}, \mathbf{y}, t_w)$ is not a stochastic variable. However, it is a stochastic variable if we regard the variations induced by the choice of couplings $J_{\mathbf{x}, \mathbf{y}}$ and over the considered sites $(\mathbf{x}, \mathbf{y}, t_w)$. We term these stochastic variables $C_4$, without arguments.

As explained in the next paragraph, although $C_4(\mathbf{x}, \mathbf{y}, t_w)$ cannot be computed with a finite number of replicas, unbiased estimators of its moments can be computed. In particular, previous work has mostly focused on the average correlation function

$$ C_4^{\text{av}}(\mathbf{r}, t_w) = \frac{1}{L^3} \sum_{\mathbf{x}} \overline{C_4(\mathbf{x}, \mathbf{y} = \mathbf{x} + \mathbf{r}, t_w)}. \qquad (10) $$

Cubic symmetry, present in averages over the samples, allows us to average over the three equivalent displacements $\mathbf{r} = (r, 0, 0)$ and permutations. We shall use the shorthand $C_4^{\text{av}}(r, t_w)$ to indicate this average over the three equivalent $\mathbf{r}$. To compute the coherence length $\xi(t_w)$ we follow [26, 45, 46] and compute the integrals

$$ I_n(t_w) = \int_0^\infty r^n C_4^{\text{av}}(r, t_w) dr. \qquad (11) $$

Then, $\xi(t_w) = I_2(t_w)/I_1(t_w)$.

As stated above, we have simulated, as a null experiment, the link-diluted Ising model (DIL). The only difference with the Hamiltonian in (8) is in the choice of the couplings: $J_{\mathbf{x}, \mathbf{y}} = 1$ (with 70% probability) or $J_{\mathbf{x}, \mathbf{y}} = 0$ (with 30% probability). Since all couplings are positive or zero, this is a ferromagnetic system without frustration.
<figure><img src="image_6.png" /><figcaption>Figure 5. Median of the distribution $P[C_4(r = \xi)]$ in units of the first moment, $\overline{C}_4(r = \xi)$, versus $\overline{C}_4(r = \xi)$, as computed for the spin glass at temperature $T = 0.9$. We show data in logarithmic scale. Therefore, the dashed line (a power-law fit with exponent ~0.5, see SI for details) appears as a straight line.</figcaption></figure>

<table><caption>Table I. Maximum $t_w$ and coherence length reached for each of our models and simulation temperatures.</caption><thead><tr><th>$T$ or $\tilde{T}$</th><th>$t_w$(EA)</th><th>$\xi_{max}$(EA)</th><th>$t_w$(DIL)</th><th>$\xi_{max}$(DIL)</th></tr></thead><tbody><tr><td>0.7</td><td>46531866276</td><td>12</td><td>498</td><td>15</td></tr><tr><td>0.8</td><td>18734780191</td><td>15</td><td>919</td><td>21</td></tr><tr><td>0.9</td><td>15172184825</td><td>20</td><td>954</td><td>23</td></tr></tbody></table>

All our simulation and analysis procedures are identical for the DIL and EA models. The critical temperature is $T_c^{DIL} = 3.0609(5)$ [47]. Actually, this is twice the value reported in [47] 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 DIL temperatures as $T = 0.9$, $\tilde{T} = 0.8$ or $0.7$ rather than to their real value $T^{DIL} = \tilde{T} (T_c^{DIL}/T_c^{EA})$.

The range of coherence length and simulation times in this study can be found in Table I.

### Unbiased estimators of powers of $C_4(\mathbf{x}, \mathbf{y}, t_w)$

Given $\mathbf{x}$ and $\mathbf{y}$, we need an unbiased estimator of $C_4^q(\mathbf{x}, \mathbf{y}, t_w) = \langle S_{\mathbf{x}}(t_w)S_{\mathbf{y}}(t_w) \rangle^{2q}$. Note that the $q = 1$ instance is needed to evaluate (10).

Should we have (at least) 2q replicas at our disposal, a tentative solution would be provided by the estimator

$$[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q^{\text{poor}} = \prod_{a=1}^{2q} S_{\mathbf{x}}^{(a)}(t_w)S_{\mathbf{y}}^{(a)}(t_w). \quad (12)$$

$[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q^{\text{poor}} = (-1)^p$, where $p$ is the number of replicas for which $S_{\mathbf{x}}^{(a)}(t_w)S_{\mathbf{y}}^{(a)}(t_w) = -1$. However, the statistical independence of the different replicas ensures for the expectation value $\langle [C_4(\mathbf{x}, \mathbf{y}, t_w)]_q^{\text{poor}} \rangle = \langle S_{\mathbf{x}}(t_w)S_{\mathbf{y}}(t_w) \rangle^{2q}$.

Nevertheless, if we have at our disposal a number of replicas $N_R \gg 2q$, as is our case, the solution in (12)

<figure><img src="image_7.png" /><figcaption>Figure 6. Scaling exponent $\tau(q)$ for the $q$-th moment $\overline{C}_4(r = \xi)^q \sim \xi^{-\tau(q)}$ computed from simulations of the Ising spin glass at $T = 0.9$ (see SI for results at $T_c$). The non-linear behavior is a strong indication of multifractality. The dashed lines are fits to the functional forms in (3) (the goodness-of-fit statistics are presented in the SI). The **inset** presents the same data as a function of $\log(q)$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Legendre transformation $f(\alpha)$ of function $\tau(q)$, see (7), as a function $\alpha = d\tau/dq$, computed from the fitting ansätze in (3). Errors on both axes have been obtained as explained in **Methods** [for further details, see 44]. Since $f(\alpha)$ is the large-deviations function of the decay rate $C_4(r = \xi) \sim r^{-\alpha}$, the data show that the majority of sites have the median decay rate of approximately 0.65, much larger than the mean decay rate $\alpha \approx 0.4$.</figcaption></figure>

is very unsatisfactory. Rather, one would like to consider all possible picks of 2q different replicas (out of the $N_R$ possible choices), compute $[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q^{\text{poor}}$ for every pick, and take the average of those products.

To achieve our goal, we have solved the following auxiliary combinatorial problem. Given a set of $N_R$ different signs $c_a = \pm 1$, $M$ of which negative, we have computed $\tilde{P}(N_R, M; S, p)$, namely the probability of getting $p$ negative signs in a pick (with uniform probability) of $S$ distinct signs. With this probability in our hands, the solution is straightforward. We just need to look at our
<figure><img src="image_6.png" /><figcaption>Figure 5. Median of the distribution $P[C_4(r = \xi)]$ in units of the first moment, $\overline{C}_4(r = \xi)$, versus $\overline{C}_4(r = \xi)$, as computed for the spin glass at temperature $T = 0.9$. We show data in logarithmic scale. Therefore, the dashed line (a power-law fit with exponent ~0.5, see SI for details) appears as a straight line.</figcaption></figure>

<table><caption>Table I. Maximum $t_w$ and coherence length reached for each of our models and simulation temperatures.</caption><thead><tr><th>T or T&#771;</th><th>t<sub>w</sub>(EA)</th><th>&xi;<sub>max</sub>(EA)</th><th>t<sub>w</sub>(DIL)</th><th>&xi;<sub>max</sub>(DIL)</th></tr></thead><tbody><tr><td>0.7</td><td>46531866276</td><td>12</td><td>498</td><td>15</td></tr><tr><td>0.8</td><td>18734780191</td><td>15</td><td>919</td><td>21</td></tr><tr><td>0.9</td><td>15172184825</td><td>20</td><td>954</td><td>23</td></tr></tbody></table>

All our simulation and analysis procedures are identical for the DIL and EA models. The critical temperature is $T_c^{DIL} = 3.0609(5)$ [47]. Actually, this is twice the value reported in [47] 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 DIL temperatures as $T = 0.9$, $\tilde{T} = 0.8$ or $0.7$ rather than to their real value $T^{DIL} = \tilde{T} (T_c^{DIL}/T_c^{EA})$.

The range of coherence length and simulation times in this study can be found in Table I.

### Unbiased estimators of powers of $C_4(\mathbf{x}, \mathbf{y}, t_w)$

Given $\mathbf{x}$ and $\mathbf{y}$, we need an unbiased estimator of $C_4^q(\mathbf{x}, \mathbf{y}, t_w) = \langle S_{\mathbf{x}}(t_w)S_{\mathbf{y}}(t_w) \rangle^{2q}$. Note that the $q = 1$ instance is needed to evaluate (10).

Should we have (at least) 2q replicas at our disposal, a tentative solution would be provided by the estimator

$$[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q^{\text{poor}} = \prod_{a=1}^{2q} S_x^{(a)}(t_w)S_y^{(a)}(t_w). \quad (12)$$

$[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q^{\text{poor}} = (-1)^p$, where $p$ is the number of replicas for which $S_x^{(a)}(t_w)S_y^{(a)}(t_w) = -1$. However, the statistical independence of the different replicas ensures for the expectation value $\langle [C_4(\mathbf{x}, \mathbf{y}, t_w)]_q^{\text{poor}} \rangle = \langle S_{\mathbf{x}}(t_w)S_{\mathbf{y}}(t_w) \rangle^{2q}$.

Nevertheless, if we have at our disposal a number of replicas $N_R \gg 2q$, as is our case, the solution in (12)

<figure><img src="image_7.png" /><figcaption>Figure 6. Scaling exponent $\tau(q)$ for the $q$-th moment $\overline{C}_4(r = \xi)^q \sim \xi^{-\tau(q)}$ computed from simulations of the Ising spin glass at $T = 0.9$ (see SI for results at $T_c$). The non-linear behavior is a strong indication of multifractality. The dashed lines are fits to the functional forms in (3) (the goodness-of-fit statistics are presented in the SI). The **inset** presents the same data as a function of $\log(q)$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Legendre transformation $f(\alpha)$ of function $\tau(q)$, see (7), as a function $\alpha = d\tau/dq$, computed from the fitting ansätze in (3). Errors on both axes have been obtained as explained in **Methods** [for further details, see 44]. Since $f(\alpha)$ is the large-deviations function of the decay rate $C_4(r = \xi) \sim r^{-\alpha}$, the data show that the majority of sites have the median decay rate of approximately 0.65, much larger than the mean decay rate $\alpha \approx 0.4$.</figcaption></figure>

is very unsatisfactory. Rather, one would like to consider all possible picks of 2q different replicas (out of the $N_R$ possible choices), compute $[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q^{\text{poor}}$ for every pick, and take the average of those products.

To achieve our goal, we have solved the following auxiliary combinatorial problem. Given a set of $N_R$ different signs $c_a = \pm 1$, $M$ of which negative, we have computed $\tilde{P}(N_R, M; S, p)$, namely the probability of getting $p$ negative signs in a pick (with uniform probability) of $S$ distinct signs. With this probability in our hands, the solution is straightforward. We just need to look at our
set $S_x^{(a)}(t_w)S_y^{(a)}(t_w)$, $a = 1, 2, \dots, N_R$, count the number $M$ of them that turn out to be negative and compute the estimator

$$[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q = G(N_R, M, q), \quad (13)$$

$$G(N_R, M, q) = \sum_{p=0}^{2q} (-1)^p \tilde{P}(N_R, M; S = 2q, p). \quad (14)$$

$[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q$ is an unbiased estimator of $\langle S_x(t_w)S_y(t_w) \rangle^{2q}$, because it is an average over all possible (poor, but unbiased) estimators in (12). Our computation of $\tilde{P}(N_R, M; S = 2q, p)$ is explained in the SI.

## The probability distribution of the correlation function

We wish to study the probability distribution function (pdf) for $\langle S_x(t_w)S_{x+r}(t_w) \rangle^2$ (periodic boundary conditions are assumed for $\mathbf{x} + \mathbf{r}$). We have only considered displacements $\mathbf{r} = (r, 0, 0)$ —and permutations— and we have chosen the measuring times in such a way that $r = \xi(t_w)$.

Note that, given the starting point $\mathbf{x}$ and the sample $\{J_{\mathbf{x},\mathbf{y}}\}$, $\langle S_{\mathbf{x}}S_{\mathbf{x}+\mathbf{r}} \rangle^2$ is not a fluctuating quantity. Hence, we are referring to the pdf as $\mathbf{x}$ and the sample vary. $\langle S_{\mathbf{x}}S_{\mathbf{x}+\mathbf{r}} \rangle^2$ can be computed exactly only in the limit $N_R \to \infty$. However, as explained in the previous paragraph, we can compute without bias its $q$-th moment provided that the number of replicas at our disposal is $N_R \ge 2q$.

The basic object we compute from our simulation is the pdf $P(M; N_R, \xi)$, namely the probability, as computed over the starting point $\mathbf{x}$ and the samples, that exactly $M$ of the $N_R$ signs $S_x^{(a)}(t_w)S_{x+r}^{(a)}(t_w)$ turn out to be $-1$ in our simulation of this specific sample. Hence, the unbiased estimator of the $q$-th moment of $\langle S_{\mathbf{x}}(t_w)S_{\mathbf{x}+\mathbf{r}(t_w)} \rangle^2$ with $r = \xi(t_w)$ is

$$\overline{C_4^q}(\xi) = \sum_{M=0}^{N_R} P(M; N_R, \xi)G(N_R, M, q), \quad (15)$$

where $G(N_R, M, q)$ was defined in (14).

Unfortunately, the median of the pdf for $\langle S_x(t_w)S_{x+r}(t_w) \rangle^2$ is more difficult to compute. Our strategy, explained in full detail in the SI, consists in computing biased estimators of the median, with bias of order $1/N_R$. Then we compute these biased estimators for a sequence $N_R' = 32, 64, 128, 256$ and $512$, and proceed to an extrapolation $N_R \to \infty$. We obtain the $P(M'; N_R', \xi)$ from their $N_R = 512$ counterpart as

$$P(M'; N_R', \xi) = \sum_{M=0}^{N_R} P(M; N_R, \xi) \times P(N_R, M; S = N_R', p = M'). \quad (16)$$

The probabilities $P(N_R, M; S, p)$ were defined in the previous subsection in this **Methods** section.

## Computation of $\tau(q)$

In order to minimize corrections to scaling, we have fitted the normalized moments as

$$\frac{\overline{C_4^q}}{\overline{C_4^q}} = \left[ \frac{A_q}{\overline{C_4}} \right]^{\beta(q)}, \quad \tau(q) = \tau(1)[q - \beta(q)]. \quad (17)$$

Fig. 3 provides an example. In order to obtain good fits, we have needed to discard (at most) one data point corresponding to the smallest $\xi(t_w)$. An advantage of this method is that we only need to consider the $\xi(t_w)$ dependence to obtain $\tau(1)$, as shown in Fig. 2. The full procedure is illustrated in the SI.

To compute errors we have followed the strategy of [44], namely carrying out all fits separately for each jackknife block (when minimizing $\chi^2$ to perform the fits, we only consider the diagonal elements of the covariance matrix). Errors in the fit parameters are obtained from the fluctuations of the jackknife blocks.

## Computation of $M(\mathbf{x}, \mathbf{r})$

The order-of-magnitude factor in (1) is computed as $\log|[C_4(\mathbf{x}, \mathbf{x} + \mathbf{r}, t_w)]_1| / \log C_4^{\text{av}}(\sqrt{r_x^2 + r_y^2 + r_z^2})$, where $[C_4(\mathbf{x}, \mathbf{x} + \mathbf{r}, t_w)]_1$ is the $q=1$ estimator in (13) (as computed with $N_R = 512$). $C_4^{\text{av}}(r)$ is interpolated to non-integer arguments using a fit obtained from data with integer $r$ (see SI).

## ACKNOWLEDGMENTS

We thank Roberto Benzi and Luca Biferale for relevant discussions. This work was supported in part by Grants No.PID2022-136374NB-C21, PID2020-112936GB-I00, PID2019-103939RB-I00, PGC2018-094684-B-C21, PGC2018-094684-B-C22 and PID2021-125506NA-I00, funded by MCIN/AEI/10.13039/501100011033 by “ERDF A way of making Europe” and by the European Union. The research has received financial support from the Simons Foundation (grant No. 454949, G. Parisi) and ICSC – Italian Research Center on High Performance Computing, Big Data and Quantum Computing, funded by European Union – NextGenerationEU. DY was supported by the Chan Zuckerberg Biohub. IGAP was supported by the Ministerio de Ciencia, Innovación y Universidades (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
set $S_x^{(a)}(t_w)S_y^{(a)}(t_w)$, $a = 1, 2, \dots, N_R$, count the number $M$ of them that turn out to be negative and compute the estimator

$$[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q = G(N_R, M, q), \quad (13)$$

$$G(N_R, M, q) = \sum_{p=0}^{2q} (-1)^p \tilde{P}(N_R, M; S = 2q, p). \quad (14)$$

$[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q$ is an unbiased estimator of $\langle S_x(t_w)S_y(t_w) \rangle^{2q}$, because it is an average over all possible (poor, but unbiased) estimators in (12). Our computation of $\tilde{P}(N_R, M; S = 2q, p)$ is explained in the SI.

## The probability distribution of the correlation function

We wish to study the probability distribution function (pdf) for $\langle S_x(t_w)S_{x+r}(t_w) \rangle^2$ (periodic boundary conditions are assumed for $\mathbf{x} + \mathbf{r}$). We have only considered displacements $\mathbf{r} = (r, 0, 0)$ —and permutations— and we have chosen the measuring times in such a way that $r = \xi(t_w)$.

Note that, given the starting point $\mathbf{x}$ and the sample $\{J_{\mathbf{x},\mathbf{y}}\}$, $\langle S_{\mathbf{x}}S_{\mathbf{x}+\mathbf{r}} \rangle^2$ is not a fluctuating quantity. Hence, we are referring to the pdf as $\mathbf{x}$ and the sample vary. $\langle S_{\mathbf{x}}S_{\mathbf{x}+\mathbf{r}} \rangle^2$ can be computed exactly only in the limit $N_R \to \infty$. However, as explained in the previous paragraph, we can compute without bias its $q$-th moment provided that the number of replicas at our disposal is $N_R \ge 2q$.

The basic object we compute from our simulation is the pdf $P(M; N_R, \xi)$, namely the probability, as computed over the starting point $\mathbf{x}$ and the samples, that exactly $M$ of the $N_R$ signs $S_x^{(a)}(t_w)S_{x+r}^{(a)}(t_w)$ turn out to be $-1$ in our simulation of this specific sample. Hence, the unbiased estimator of the $q$-th moment of $\langle S_{\mathbf{x}}(t_w)S_{\mathbf{x}+\mathbf{r}(t_w)} \rangle^2$ with $r = \xi(t_w)$ is

$$\overline{C_4^q}(\xi) = \sum_{M=0}^{N_R} P(M; N_R, \xi)G(N_R, M, q), \quad (15)$$

where $G(N_R, M, q)$ was defined in (14).

Unfortunately, the median of the pdf for $\langle S_x(t_w)S_{x+r}(t_w) \rangle^2$ is more difficult to compute. Our strategy, explained in full detail in the SI, consists in computing biased estimators of the median, with bias of order $1/N_R$. Then we compute these biased estimators for a sequence $N_R' = 32, 64, 128, 256$ and $512$, and proceed to an extrapolation $N_R \to \infty$. We obtain the $P(M'; N_R', \xi)$ from their $N_R = 512$ counterpart as

$$P(M'; N_R', \xi) = \sum_{M=0}^{N_R} P(M; N_R, \xi) \times P(N_R, M; S = N_R', p = M'). \quad (16)$$

The probabilities $P(N_R, M; S, p)$ were defined in the previous subsection in this **Methods** section.

## Computation of $\tau(q)$

In order to minimize corrections to scaling, we have fitted the normalized moments as

$$\frac{\overline{C_4^q}}{\overline{C_4^q}} = \left[ \frac{A_q}{\overline{C_4}} \right]^{\beta(q)}, \quad \tau(q) = \tau(1)[q - \beta(q)]. \quad (17)$$

Fig. 3 provides an example. In order to obtain good fits, we have needed to discard (at most) one data point corresponding to the smallest $\xi(t_w)$. An advantage of this method is that we only need to consider the $\xi(t_w)$ dependence to obtain $\tau(1)$, as shown in Fig. 2. The full procedure is illustrated in the SI.

To compute errors we have followed the strategy of [44], namely carrying out all fits separately for each jackknife block (when minimizing $\chi^2$ to perform the fits, we only consider the diagonal elements of the covariance matrix). Errors in the fit parameters are obtained from the fluctuations of the jackknife blocks.

## Computation of $M(\mathbf{x}, \mathbf{r})$

The order-of-magnitude factor in (1) is computed as $\log|[C_4(\mathbf{x}, \mathbf{x} + \mathbf{r}, t_w)]_1| / \log C_4^{\text{av}}(\sqrt{r_x^2 + r_y^2 + r_z^2})$, where $[C_4(\mathbf{x}, \mathbf{x} + \mathbf{r}, t_w)]_1$ is the $q=1$ estimator in (13) (as computed with $N_R = 512$). $C_4^{\text{av}}(r)$ is interpolated to non-integer arguments using a fit obtained from data with integer $r$ (see SI).

## ACKNOWLEDGMENTS

We thank Roberto Benzi and Luca Biferale for relevant discussions. This work was supported in part by Grants No.PID2022-136374NB-C21, PID2020-112936GB-I00, PID2019-103939RB-I00, PGC2018-094684-B-C21, PGC2018-094684-B-C22 and PID2021-125506NA-I00, funded by MCIN/AEI/10.13039/501100011033 by "ERDF A way of making Europe" and by the European Union. The research has received financial support from the Simons Foundation (grant No. 454949, G. Parisi) and ICSC - Italian Research Center on High Performance Computing, Big Data and Quantum Computing, funded by European Union - NextGenerationEU. DY was supported by the Chan Zuckerberg Biohub. IGAP was supported by the Ministerio de Ciencia, Innovación y Universidades (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. IP was supported by LazioInnova-Regione Lazio under the program Gruppi di ricerca2020 - POR FESR Lazio 2014-2020, Project NanoProbe (Application code A0375-2020-36761).

## SUPPORTING INFORMATION

### The auxiliary combinatorial problem

### The computation of $\tilde{P}(N_R, M; S, p)$

As stated in the main text, we need to consider the following problem: given a set of $N_R$ different signs $c_a = \pm 1$, $M$ of which negative, we need to obtain $\tilde{P}(N_R, M; S, p)$, namely the probability of getting $p$ negative signs in a pick (with uniform probability) of $S$ distinct signs. In order to organize the computation, we shall consider both the sign labels and the pick orderings as distinguishable. Hence, the number of possible picks is

$$N_{\text{picks}} = \frac{N_R!}{(N_R - S)!}. \qquad (18)$$

Let us denote by $K(N_R, M; S, p)$ the number of picks of $S$ distinct signs that contain exactly $p$ negative signs. Hence,

$$\tilde{P}(N_R, M; S, p) = \frac{K(N_R, M; S, p)}{N_{\text{picks}}}. \qquad (19)$$

$K(N_R, M; S, p)$ can be obtained as a product of three factors

$$K(N_R, M; S, p) = F_1 F_2 F_3, \qquad (20)$$

where

$$\begin{align}
F_1 &= \frac{M!}{p!(M-p)!}, \nonumber \\
F_2 &= \frac{S!}{(S-p)!}, \tag{21} \\
F_3 &= \frac{(N_R - M)!}{(N_R - M - S + p)!}. \nonumber
\end{align}$$

In the above expression, whenever the factorial of a negative integer arises it should interpreted as $\infty$ (hence the corresponding factor vanishes). The meaning of the different factors is as follows:

* $F_1$ is the number of ways that we can choose $p$ tags of negative signs among $M$ possibilities.

* $F_2$ is the number of ways in which a given set of $p$ tags of negative signs can be extracted: the first tag can be obtained in the first selection, or in the second, etc. So there are $S$ possibilities for the first tag, which leaves us with $S-1$ possibilities for the second tag, and so on.

* Finally, $F_3$ is concerned with the $S-p$ positive signs that we need to complete the pick of $S$ signs. We have $N_R - M$ choices for the first tag to be chosen, $N_R - M - 1$ for the second tag, and so forth.

We have found it preferable, however, to make use of a Pascal-Tartaglia-like relation. To obtain it, it is useful to think of a pick of $S$ signs as two consecutive picks. We get $S-1$ signs from the first pick, and the last sign $c_\alpha$ is chosen only afterwards:

$$K(N_R, M; S, p) = A K(N_R, M; S-1, p-1) + B K(N_R, M; S-1, p), \qquad (22)$$

where $A$ is the number of negative signs available for the last pick (given that we obtained $p-1$ negative signs from the first pick), while $B$ is the number of positive signs available for the last pick (given that we already obtained $p$ negative signs from the first pick), specifically:

$$A = \max\{0, M - p + 1\}, \qquad (23)$$

$$B = \max\{0, N_R - (S-1) - \max\{0, M-p\}\}. \qquad (24)$$

The recursion relation for $K(N_R, M; S, p)$ instantaneously translates to a recursion relation for the probability:

$$\tilde{P}(N_R, M; S, p) = \frac{A}{N_R - S + 1} \tilde{P}(N_R, M; S - 1, p - 1) + \frac{B}{N_R - S + 1} \tilde{P}(N_R, M; S - 1, p). \qquad (25)$$

Starting the recursion from

$$\begin{align}
P(N_R, M; S = 1, p = 0) &= \frac{N_R - M}{N_R}, \nonumber \\
P(N_R, M; S = 1, p = 1) &= \frac{M}{N_R}, \tag{26}
\end{align}$$

we can simply and accurately compute $P(N_R, M; S, p)$ for whatever values of $M, S$ and $p$ we need, respecting, of course, the obvious bounds:

$$N_R \ge M, S, p, \quad M \ge p, \quad S \ge p. \qquad (27)$$

### A useful symmetry

The probability that we have computed in the previous paragraph presents a spin-flip symmetry

$$\tilde{P}(N_R, M; S, p) = \tilde{P}(N_R, N_R - M; S, S - p). \qquad (28)$$

The simplest way to show that the symmetry is present is noticing that the map of the signs set $\{c_a\}_{a=1}^{N_R}$ into $\{-c_a\}_{a=1}^{N_R}$ is bijective, transforming $M$ into $N_R - M$ and a pick of $p$ negative signs into a pick of $S-p$ negative signs.
grant. IP was supported by LazioInnova-Regione Lazio under the program Gruppi di ricerca2020 - POR FESR Lazio 2014-2020, Project NanoProbe (Application code A0375-2020-36761).

## SUPPORTING INFORMATION

### The auxiliary combinatorial problem

#### The computation of $\tilde{P}(N_R, M; S, p)$

As stated in the main text, we need to consider the following problem: given a set of $N_R$ different signs $c_a = \pm 1$, $M$ of which negative, we need to obtain $\tilde{P}(N_R, M; S, p)$, namely the probability of getting $p$ negative signs in a pick (with uniform probability) of $S$ distinct signs. In order to organize the computation, we shall consider both the sign labels and the pick orderings as distinguishable. Hence, the number of possible picks is

$$N_{\text{picks}} = \frac{N_R!}{(N_R - S)!}. \qquad (18)$$

Let us denote by $K(N_R, M; S, p)$ the number of picks of $S$ distinct signs that contain exactly $p$ negative signs. Hence,

$$\tilde{P}(N_R, M; S, p) = \frac{K(N_R, M; S, p)}{N_{\text{picks}}}. \qquad (19)$$

$K(N_R, M; S, p)$ can be obtained as a product of three factors

$$K(N_R, M; S, p) = F_1 F_2 F_3, \qquad (20)$$

where

$$\begin{align*} F_1 &= \frac{M!}{p!(M-p)!}, \\ F_2 &= \frac{S!}{(S-p)!}, \\ F_3 &= \frac{(N_R - M)!}{(N_R - M - S + p)!}. \end{align*} \qquad (21)$$

In the above expression, whenever the factorial of a negative integer arises it should be interpreted as $\infty$ (hence the corresponding factor vanishes). The meaning of the different factors is as follows:

* $F_1$ is the number of ways that we can choose $p$ tags of negative signs among $M$ possibilities.

* $F_2$ is the number of ways in which a given set of $p$ tags of negative signs can be extracted: the first tag can be obtained in the first selection, or in the second, etc. So there are $S$ possibilities for the first tag, which leaves us with $S-1$ possibilities for the second tag, and so on.

* Finally, $F_3$ is concerned with the $S-p$ positive signs that we need to complete the pick of $S$ signs. We have $N_R - M$ choices for the first tag to be chosen, $N_R - M - 1$ for the second tag, and so forth.

We have found it preferable, however, to make use of a Pascal-Tartaglia-like relation. To obtain it, it is useful to think of a pick of $S$ signs as two consecutive picks. We get $S-1$ signs from the first pick, and the last sign $c_\alpha$ is chosen only afterwards:

$$K(N_R, M; S, p) = A K(N_R, M; S-1, p-1) + B K(N_R, M; S-1, p), \qquad (22)$$

where $A$ is the number of negative signs available for the last pick (given that we obtained $p-1$ negative signs from the first pick), while $B$ is the number of positive signs available for the last pick (given that we already obtained $p$ negative signs from the first pick), specifically:

$$A = \max\{0, M - p + 1\}, \qquad (23)$$

$$B = \max\{0, N_R - (S-1) - \max\{0, M-p\}\}. \qquad (24)$$

The recursion relation for $K(N_R, M; S, p)$ instantaneously translates to a recursion relation for the probability:

$$\tilde{P}(N_R, M; S, p) = \frac{A}{N_R - S + 1} \tilde{P}(N_R, M; S - 1, p - 1) + \frac{B}{N_R - S + 1} \tilde{P}(N_R, M; S - 1, p). \qquad (25)$$

Starting the recursion from

$$\begin{align} P(N_R, M; S = 1, p = 0) &= \frac{N_R - M}{N_R}, \notag \\ P(N_R, M; S = 1, p = 1) &= \frac{M}{N_R}, \tag{26} \end{align}$$

we can simply and accurately compute $P(N_R, M; S, p)$ for whatever values of $M, S$ and $p$ we need, respecting, of course, the obvious bounds:

$$N_R \ge M, S, p, \quad M \ge p, \quad S \ge p. \qquad (27)$$

### A useful symmetry

The probability that we have computed in the previous paragraph presents a spin-flip symmetry

$$\tilde{P}(N_R, M; S, p) = \tilde{P}(N_R, N_R - M; S, S - p). \qquad (28)$$

The simplest way to show that the symmetry is present is noticing that the map of the signs set $\{c_a\}_{a=1}^{N_R}$ into $\{-c_a\}_{a=1}^{N_R}$ is bijective, transforming $M$ into $N_R - M$ and a pick of $p$ negative signs into a pick of $S-p$ negative signs.
(28) implies a symmetry in the quantity we use to estimate the moments of the correlation function, recall the main text,

$$[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q = G(N_R, M, q), \quad (29)$$

$$G(N_R, M, q) = \sum_{p=0}^{2q} (-1)^p \tilde{P}(N_R, M; S=2q, p). \quad (30)$$

Hence, the symmetry in (28) implies

$$G(N_R, M, q) = G(N_R, N_R - M, q). \quad (31)$$

This consideration is important when we aim to compute the median of the stochastic variable $C_4$ (remember that we only have in our hands estimators such as $[C_4]_{q=1}$) A moment's thought reveals that, for even $N_R$, the minimum value of $G(N_R, M, q = 1)$ is reached at $M = N_R/2$. Indeed,

$$G(N_R, M, q = 1) = \frac{(N_R - 2M)^2 - N_R}{N_R(N_R - 1)}. \quad (32)$$

Hence, one may like to use a symmetrized probability, rather than the probability discussed in the main text, $\mathcal{P}(M; N_R, \xi)$, (the probability, as computed over the starting point $\mathbf{x}$ and the samples, that exactly $M$ of the $N_R$ signs $S_x^{(a)}(t_w)S_{x+r}^{(a)}(t_w)$ turn out to be $-1$ in our simulation of this specific sample) [48]. Specifically,

$$\mathcal{P}_s \left( n = \frac{N_R}{2}; N_R, \xi \right) = \mathcal{P}(n; N_R, \xi), \quad (33)$$

$$\mathcal{P}_s \left( \frac{N_R}{2} < n \le N_R; N_R, \xi \right) = \mathcal{P}_s(n; N_R, \xi) + \mathcal{P}_s(N_R - n; N_R, \xi).$$

Hence, the moments of $C_4$ can be computed as

$$\overline{C_4^q}(\xi) = \sum_{n=N_R/2}^{N_R} \mathcal{P}_s(n; N_R, \xi) G(N_R, n, q). \quad (34)$$

It is then clear that we need to estimate medians from the symmetrized probability $\mathcal{P}_s(n; N_R, \xi)$, rather than from $\mathcal{P}(M; N_R, \xi)$.

## The computation of the medians

### Computing the medians from our biased estimators

Our analysis starts from the symmetrized probability $\mathcal{P}_s(n; N_R, \xi)$. We begin by computing the cumulative distribution

$$\mathcal{S}(n) = \sum_{\ell=N_R/2}^{n} \mathcal{P}_s(\ell; N_R, \xi), \quad (35)$$

and determine the smallest integer $n^* \ge N_R/2$ such that $\mathcal{S}(n^*) > 0.5$. Indeed, we assign positive weights $\omega_-$ and $\omega_+$ to $(n^* - 1)$ and $n^*$, in such a way that

$$\omega_- + \omega_+ = 1, \quad \omega_- \mathcal{S}(n^* - 1) + \omega_+ \mathcal{S}(n^*) = 0.5. \quad (36)$$

In other words we linearly interpolate $\mathcal{S}(n)$ in $n$, in order to obtain a cumulative exactly equal to 0.5.

However, we need to translate the integer $n^*$ into a real number, namely an estimator of the median of $C_4$. In order to do that, we need to take one step back and consider the physical meaning of the fact that exactly $M$ of the $N_R$ signs $S_x^{(a)}(t_w)S_y^{(a)}(t_w)$ turn out to be $-1$ in our simulation of a specific sample. Indeed, one immediately notices that

$$\begin{aligned} \frac{N_R - 2M}{N_R} &= \frac{1}{N_R} \sum_{a=1}^{N_R} S_x^{(a)}(t_w) S_y^{(a)}(t_w) \quad (37) \\ &= \langle S_x(t_w) S_y(t_w) \rangle + \eta \sqrt{\frac{1 - C_4}{N_R}}, \end{aligned}$$

where $C_4 = \langle S_x(t_w)S_y(t_w) \rangle^2$ and $\eta$ is a random variable that verifies

$$\langle \eta \rangle = 0, \quad \langle \eta^2 \rangle = 1. \quad (38)$$

Furthermore, $\eta$ is statistically uncorrelated with $\langle S_x(t_w)S_y(t_w) \rangle$ and, because of the central limit theorem, tends to a normally distributed stochastic variable in the limit of large $N_R$. Hence, we may consider the (very biased) estimator of $C_4$ ($\operatorname{sgn}(x) = x/|x|$):

$$\begin{aligned} [S_x(t_w)S_y(t_w)]^2 &\equiv \left(\frac{N_R - 2M}{N_R}\right)^2 \quad (39) \\ &= C_4 + 2\eta \operatorname{sgn}(\langle S_x(t_w)S_y(t_w) \rangle) \sqrt{\frac{C_4(1 - C_4)}{N_R}} \\ &\quad + \eta^2 \frac{(1 - C_4)}{N_R}. \end{aligned}$$

We do not expect problems from the term linear in $\eta$. Indeed for $N_R$ large enough, $\eta$ approaches a normal variable (so, symmetrically distributed around zero), hence it should not cause bias on the estimation of the median. The real problem comes from the term proportional to $\eta^2$, inducing departures from the true value of $C_4$ of order $1/N_R$. Hence we may consider our first estimator of the median

$$\begin{aligned} \operatorname{median}([S_x(t_w)S_y(t_w)]^2) = & \omega_- \left( \frac{N_R - 2(n^* - 1)}{N_R} \right)^2 \\ & + \omega_+ \left( \frac{N_R - 2n^*}{N_R} \right)^2. \end{aligned} \quad (40)$$

As explained above, we expect that this estimator of the median of $C_4$ will have a bias of order $1/N_R$.
(28) implies a symmetry in the quantity we use to estimate the moments of the correlation function, recall the main text,

$$[C_4(\mathbf{x}, \mathbf{y}, t_w)]_q = G(N_R, M, q), \quad (29)$$

$$G(N_R, M, q) = \sum_{p=0}^{2q} (-1)^p \tilde{P}(N_R, M; S=2q, p). \quad (30)$$

Hence, the symmetry in (28) implies

$$G(N_R, M, q) = G(N_R, N_R - M, q). \quad (31)$$

This consideration is important when we aim to compute the median of the stochastic variable $C_4$ (remember that we only have in our hands estimators such as $[C_4]_{q=1}$) A moment's thought reveals that, for even $N_R$, the minimum value of $G(N_R, M, q = 1)$ is reached at $M = N_R/2$. Indeed,

$$G(N_R, M, q = 1) = \frac{(N_R - 2M)^2 - N_R}{N_R(N_R - 1)}. \quad (32)$$

Hence, one may like to use a symmetrized probability, rather than the probability discussed in the main text, $\mathcal{P}(M; N_R, \xi)$, (the probability, as computed over the starting point $\mathbf{x}$ and the samples, that exactly $M$ of the $N_R$ signs $S_x^{(a)}(t_w)S_{x+r}^{(a)}(t_w)$ turn out to be $-1$ in our simulation of this specific sample) [48]. Specifically,

$$\mathcal{P}_s \left( n = \frac{N_R}{2}; N_R, \xi \right) = \mathcal{P}(n; N_R, \xi), \quad (33)$$

$$\mathcal{P}_s \left( \frac{N_R}{2} < n \le N_R; N_R, \xi \right) = \mathcal{P}_s(n; N_R, \xi) + \mathcal{P}_s(N_R - n; N_R, \xi).$$

Hence, the moments of $C_4$ can be computed as

$$\overline{C_4^q}(\xi) = \sum_{n=N_R/2}^{N_R} \mathcal{P}_s(n; N_R, \xi) G(N_R, n, q). \quad (34)$$

It is then clear that we need to estimate medians from the symmetrized probability $\mathcal{P}_s(n; N_R, \xi)$, rather than from $\mathcal{P}(M; N_R, \xi)$.

## The computation of the medians

### Computing the medians from our biased estimators

Our analysis starts from the symmetrized probability $\mathcal{P}_s(n; N_R, \xi)$. We begin by computing the cumulative distribution

$$\mathcal{S}(n) = \sum_{\ell=N_R/2}^{n} \mathcal{P}_s(\ell; N_R, \xi), \quad (35)$$

and determine the smallest integer $n^* \ge N_R/2$ such that $\mathcal{S}(n^*) > 0.5$. Indeed, we assign positive weights $\omega_-$ and $\omega_+$ to $(n^* - 1)$ and $n^*$, in such a way that

$$\omega_- + \omega_+ = 1, \quad \omega_- \mathcal{S}(n^* - 1) + \omega_+ \mathcal{S}(n^*) = 0.5. \quad (36)$$

In other words we linearly interpolate $\mathcal{S}(n)$ in $n$, in order to obtain a cumulative exactly equal to 0.5.

However, we need to translate the integer $n^*$ into a real number, namely an estimator of the median of $C_4$. In order to do that, we need to take one step back and consider the physical meaning of the fact that exactly $M$ of the $N_R$ signs $S_x^{(a)}(t_w)S_y^{(a)}(t_w)$ turn out to be $-1$ in our simulation of a specific sample. Indeed, one immediately notices that

$$\begin{aligned} \frac{N_R - 2M}{N_R} &= \frac{1}{N_R} \sum_{a=1}^{N_R} S_x^{(a)}(t_w) S_y^{(a)}(t_w) \quad (37) \\ &= \langle S_x(t_w) S_y(t_w) \rangle + \eta \sqrt{\frac{1 - C_4}{N_R}}, \end{aligned}$$

where $C_4 = \langle S_x(t_w)S_y(t_w) \rangle^2$ and $\eta$ is a random variable that verifies

$$\langle \eta \rangle = 0, \quad \langle \eta^2 \rangle = 1. \quad (38)$$

Furthermore, $\eta$ is statistically uncorrelated with $\langle S_x(t_w)S_y(t_w) \rangle$ and, because of the central limit theorem, tends to a normally distributed stochastic variable in the limit of large $N_R$. Hence, we may consider the (very biased) estimator of $C_4$ ($\operatorname{sgn}(x) = x/|x|$):

$$\begin{aligned} [S_x(t_w)S_y(t_w)]^2 &\equiv \left(\frac{N_R - 2M}{N_R}\right)^2 \quad (39) \\ &= C_4 + 2\eta \operatorname{sgn}(\langle S_x(t_w)S_y(t_w) \rangle) \sqrt{\frac{C_4(1 - C_4)}{N_R}} \\ &\quad + \eta^2 \frac{(1 - C_4)}{N_R}. \end{aligned}$$

We do not expect problems from the term linear in $\eta$. Indeed for $N_R$ large enough, $\eta$ approaches a normal variable (so, symmetrically distributed around zero), hence it should not cause bias on the estimation of the median. The real problem comes from the term proportional to $\eta^2$, inducing departures from the true value of $C_4$ of order $1/N_R$. Hence we may consider our first estimator of the median

$$\begin{aligned} \operatorname{median}([S_x(t_w)S_y(t_w)]^2) = & \omega_- \left( \frac{N_R - 2(n^* - 1)}{N_R} \right)^2 \\ & + \omega_+ \left( \frac{N_R - 2n^*}{N_R} \right)^2. \end{aligned} \quad (40)$$

As explained above, we expect that this estimator of the median of $C_4$ will have a bias of order $1/N_R$.
<figure><img src="image_9.png" /><figcaption>Figure 8. Biased estimators of the median, see (39) and (41), as functions of the inverse number of replicas, $N_R$, for $T = 0.9$ and $\xi(t_w) = 20$. The dashed lines are fits to quadratic functions in order to extrapolate to $N_R \to \infty$ (complete fit statistics are available in Table II).</figcaption></figure>

Alternatively, we may consider an unbiased estimator of $C_4$:

$$ 
\begin{aligned}
[C_4] &= G(N_R, M, q = 1) = \frac{(N_R - 2M)^2 - N_R}{N_R(N_R - 1)} \quad (41) \\
&= C_4 + 2\eta \operatorname{sgn}(\langle S_x(t_w)S_y(t_w) \rangle) \sqrt{\frac{N_R C_4 (1 - C_4)}{(N_R - 1)^2}} \\
&\quad + \frac{(\eta^2 - 1)}{N_R - 1}.
\end{aligned}
$$

Notice that, although the expectation value of $\eta^2 - 1$ is zero, this term is not symmetrically distributed around zero (not even in the limit $N_R \to \infty$, when $\eta$ is normally distributed). Hence the $\eta^2 - 1$ will also distort the computation of the median by a quantity of order $1/N_R$. Accordingly, we expect corrections of order $1/N_R$ for the corresponding estimator of the bias of the median of $C_4$

$$ 
\begin{aligned}
\operatorname{median}([C_4]) &= \omega_- G(N_R, n^* - 1, q = 1) \\
&\quad + \omega_+ G(N_R, n^*, q = 1).
\end{aligned}
\quad (42)
$$

The extrapolation to $N_R \to \infty$ from both biased estimators is illustrated in Fig. 8.

### The probability distribution near $C_4 = 1$

As we have shown in the main text, the scaling exponent for the q-th moment $\tau(q)$, i.e., $\overline{C_4^q} \sim 1/[\xi(t_w)]^{\tau(q)}$, goes as $\tau(q) \sim \log q$ for large $q$. This is to be expected if the probability distribution function for $C_4$ near 1 behaves as

$$ P(C_4 \to 1) \propto (1 - C_4)^{B(\xi)}, \quad (43) $$

<figure><img src="image_10.png" /><figcaption>Figure 9. Median of the distribution $P[C_4(r = \xi)]$ versus the first moment, $\overline{C_4}(r = \xi)$, for $T = 0.9$. The median is calculated through the extrapolation to $N_R \to \infty$ shown in Fig. 8. The dashed line is a fit to a power law in order to obtain the slope of $\tau(q)$ near $q = 0$ (complete fit information available in Table II).</figcaption></figure>

with an exponent $B(\xi)$ that grows logarithmically with $\xi$. Indeed, a simple saddle-point estimation yields

$$ 
\begin{aligned}
\overline{C_4^q} &= \int_0^1 dC_4 P(C_4) C_4^q \quad (44) \\
&\sim \int_0^1 dC_4 \exp [B(\log(1 - C_4) + q \log C_4)] \approx \frac{1}{q^B}.
\end{aligned}
$$

Hence, if $B(\xi) \sim A \log(\xi)$ ($A$ is an amplitude)

$$ \frac{1}{q^{A \log \xi}} = \frac{1}{\xi^{A \log q}}. \quad (45) $$

Equipped with the above intuition, we have checked the $P(C_4 \to 1)$ as computed from the $N_R = 512$ estimator in (41). We obtain a good fit to $(1 - C_4)^{B(\xi)}$, but the determination of the exponent $B(\xi)$ is quite difficult, as it depends significantly on the fitting range. This is why we have turned to a different strategy. We have considered the following expectation value

$$ \mathcal{I}(N_R) = \sum_{n=N_R/2}^{N_R} P_s(n; N_R, \xi) e^{-4(N_R - n)}. \quad (46) $$

The above sum is clearly dominated by values of $n$ near $N_R$ where, recall (39), $C_4 \approx 1 - 4(N_R - n)/N_R$. Hence, we can approximate

$$ 
\begin{aligned}
\mathcal{I}(N_R) &\approx \int_0^1 dC_4 P(C_4) e^{-N_R(1-C_4)} \quad (47) \\
&\sim \int_0^1 dC_4 (1-C_4)^{B(\xi)} e^{-N_R(1-C_4)} \sim \frac{1}{N_R^{1+B(\xi)}}.
\end{aligned}
$$

Then, our chosen strategy has been to fit our data for $\mathcal{I}(N_R)$ as a power law in $1/N_R$, see Fig. 10. Generally
<figure><img src="image_9.png" /><figcaption>Figure 8. Biased estimators of the median, see (39) and (41), as functions of the inverse number of replicas, $N_R$, for $T = 0.9$ and $\xi(t_w) = 20$. The dashed lines are fits to quadratic functions in order to extrapolate to $N_R \to \infty$ (complete fit statistics are available in Table II).</figcaption></figure>

Alternatively, we may consider an unbiased estimator of $C_4$:

$$ 
\begin{aligned}
[C_4] &= G(N_R, M, q = 1) = \frac{(N_R - 2M)^2 - N_R}{N_R(N_R - 1)} \quad (41) \\
&= C_4 + 2\eta \operatorname{sgn}(\langle S_x(t_w)S_y(t_w) \rangle) \sqrt{\frac{N_R C_4 (1 - C_4)}{(N_R - 1)^2}} \\
&\quad + \frac{(\eta^2 - 1)}{N_R - 1}.
\end{aligned}
$$

Notice that, although the expectation value of $\eta^2 - 1$ is zero, this term is not symmetrically distributed around zero (not even in the limit $N_R \to \infty$, when $\eta$ is normally distributed). Hence the $\eta^2 - 1$ will also distort the computation of the median by a quantity of order $1/N_R$. Accordingly, we expect corrections of order $1/N_R$ for the corresponding estimator of the bias of the median of $C_4$

$$ 
\begin{aligned}
\operatorname{median}([C_4]) &= \omega_- G(N_R, n^* - 1, q = 1) \\
&\quad + \omega_+ G(N_R, n^*, q = 1).
\end{aligned}
\quad (42)
$$

The extrapolation to $N_R \to \infty$ from both biased estimators is illustrated in Fig. 8.

### The probability distribution near $C_4 = 1$

As we have shown in the main text, the scaling exponent for the q-th moment $\tau(q)$, i.e., $\overline{C_4^q} \sim 1/[\xi(t_w)]^{\tau(q)}$, goes as $\tau(q) \sim \log q$ for large $q$. This is to be expected if the probability distribution function for $C_4$ near 1 behaves as

$$ P(C_4 \to 1) \propto (1 - C_4)^{B(\xi)}, \quad (43) $$

<figure><img src="image_10.png" /><figcaption>Figure 9. Median of the distribution $P[C_4(r = \xi)]$ versus the first moment, $\overline{C_4}(r = \xi)$, for $T = 0.9$. The median is calculated through the extrapolation to $N_R \to \infty$ shown in Fig. 8. The dashed line is a fit to a power law in order to obtain the slope of $\tau(q)$ near $q = 0$ (complete fit information available in Table II).</figcaption></figure>

with an exponent $B(\xi)$ that grows logarithmically with $\xi$. Indeed, a simple saddle-point estimation yields

$$ 
\begin{aligned}
\overline{C_4^q} &= \int_0^1 dC_4 P(C_4) C_4^q \quad (44) \\
&\sim \int_0^1 dC_4 \exp [B(\log(1 - C_4) + q \log C_4)] \approx \frac{1}{q^B}.
\end{aligned}
$$

Hence, if $B(\xi) \sim A \log(\xi)$ ($A$ is an amplitude)

$$ \frac{1}{q^{A \log \xi}} = \frac{1}{\xi^{A \log q}}. \quad (45) $$

Equipped with the above intuition, we have checked the $P(C_4 \to 1)$ as computed from the $N_R = 512$ estimator in (41). We obtain a good fit to $(1 - C_4)^{B(\xi)}$, but the determination of the exponent $B(\xi)$ is quite difficult, as it depends significantly on the fitting range. This is why we have turned to a different strategy. We have considered the following expectation value

$$ \mathcal{I}(N_R) = \sum_{n=N_R/2}^{N_R} \mathcal{P}_s(n; N_R, \xi) e^{-4(N_R-n)}. \quad (46) $$

The above sum is clearly dominated by values of $n$ near $N_R$ where, recall (39), $C_4 \approx 1 - 4(N_R - n)/N_R$. Hence, we can approximate

$$ 
\begin{aligned}
\mathcal{I}(N_R) &\approx \int_0^1 dC_4 P(C_4) e^{-N_R(1-C_4)} \quad (47) \\
&\sim \int_0^1 dC_4 (1-C_4)^{B(\xi)} e^{-N_R(1-C_4)} \sim \frac{1}{N_R^{1+B(\xi)}}.
\end{aligned}
$$

Then, our chosen strategy has been to fit our data for $\mathcal{I}(N_R)$ as a power law in $1/N_R$, see Fig. 10. Generally
<figure><img src="image_11.png" /><figcaption>Figure 10. Expectation value $\mathcal{I}(N_R)$, see (46), versus the inverse of the number of replicas $N_R$, at $T = 0.9$ and $\xi(t_w) = 20$. The dashed line is a fit to (47) with $N_R \in [64, 512]$ (full details of the fit available in Table II).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Exponent $B(\xi)$ for $T = 0.9$. The dashed line is a fit to $a + b \log(\xi)$ for $\xi \ge 10$ (see Table II for complete fit statistics). These results make more plausible a logarithmic behavior of $\tau(q)$ for large $q$.</figcaption></figure>

speaking, we obtain fair fits (although for $\xi < 10$ we had to discard the $N_R = 32$ data from the fits). The resulting exponents $B(\xi)$ are shown in Fig. 11. As it can be checked, the hypothesis $B(\xi) \propto \log(\xi)$ is tenable, particularly for $\xi > 10$.

### The $\tau(q)$ at $T_c$

We have fitted our data for $\tau(q)$ following the same functional forms that we used at $T = 0.9$, namely:

$$ \tau_1(q) = mq \frac{1 + c_1 q}{1 + c_2 q}, \quad \tau_2(q) = mq \frac{1 + d_1 q \log q}{(1 + d_2 q)^2}. \quad (48) $$

Both $\tau_1$ and $\tau_2$ have the same derivative at $q=0$, namely $m$. We do not regard $m$ as a fitting parameter. Rather we

<figure><img src="image_13.png" /><figcaption>Figure 12. The scaling exponent for the $q$-th moment $\tau(q)$, i.e., $\overline{C_4^q} \sim 1/[\xi(t_w)]^{\tau(q)}$ versus $q$, as computed for our simulations of the Ising spin glass at $T_c = 1.1019(29)$ [31] and $T = 0.9$. The dashed lines are fits to the functional forms in (48). Full fit results are shown in Table II.</figcaption></figure>

take it from the scaling of the median of the distribution $P(C_4)$ with $\xi$. As usual, we first fitted the median as a function of $\overline{C_4}$, see Fig. 9, and then translated it to a power law in $\xi$.

Fig. 12 shows our estimations and fits for $\tau(q)$ for $T = 0.9$ and $T = T_c$. Although the functional behaviors of both curves are similar, the clear difference reveals the significantly different nature of the multifractal behaviors at the critical points and in the spin-glass phase.

### Summary of the statistical information fits

As stated in the main text, and in the above sections, in order to compute and interpolate $\tau(q)$, and its errors [see 44], we need to make a huge number of fits. For example, to obtain the data shown in Fig. 9, we need to extrapolate to $N_R \to \infty$ the median of $C_4$, see Fig. 8, 16 times at each $\xi(t_w)$ in order to estimate errors. This process would be impossible without some automation. To achieve this, we use a nonlinear least-squares method using the Levenberg-Marquardt solver [49] routine available at the GNU Scientific Library. However, in order to calibrate the process and obtain a fair fits, we perform manually all the fits shown in the present work. The complete statistical information of the fits is summarized in Table II.

Fig. 13, which shows the evolution of the 10th and 20th moments of $C_4$ as function of the coherence length, includes a fit to a power law $a_q \xi^{-\tau(q)}$, where the only parameter is the amplitude constant, $a_q$, and the respective value of $\tau(q)$ is the result of our analysis. These fits, included in Table II, are a good confirmation of our analysis.
<figure><img src="image_11.png" /><figcaption>Figure 10. Expectation value $\mathcal{I}(N_R)$, see (46), versus the inverse of the number of replicas $N_R$, at $T = 0.9$ and $\xi(t_w) = 20$. The dashed line is a fit to (47) with $N_R \in [64, 512]$ (full details of the fit available in Table II).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Exponent $B(\xi)$ for $T = 0.9$. The dashed line is a fit to $a + b \log(\xi)$ for $\xi \ge 10$ (see Table II for complete fit statistics). These results make more plausible a logarithmic behavior of $\tau(q)$ for large $q$.</figcaption></figure>

Speaking, we obtain fair fits (although for $\xi < 10$ we had to discard the $N_R = 32$ data from the fits). The resulting exponents $B(\xi)$ are shown in Fig. 11. As it can be checked, the hypothesis $B(\xi) \propto \log(\xi)$ is tenable, particularly for $\xi > 10$.

### The $\tau(q)$ at $T_c$

We have fitted our data for $\tau(q)$ following the same functional forms that we used at $T = 0.9$, namely:

$$ \tau_1(q) = mq \frac{1 + c_1 q}{1 + c_2 q}, \quad \tau_2(q) = mq \frac{1 + d_1 q \log q}{(1 + d_2 q)^2}. \quad (48) $$

Both $\tau_1$ and $\tau_2$ have the same derivative at $q=0$, namely $m$. We do not regard $m$ as a fitting parameter. Rather we

<figure><img src="image_13.png" /><figcaption>Figure 12. The scaling exponent for the $q$-th moment $\tau(q)$, i.e., $C_4^q \sim 1/[\xi(t_w)]^{\tau(q)}$ versus $q$, as computed for our simulations of the Ising spin glass at $T_c = 1.1019(29)$ [31] and $T = 0.9$. The dashed lines are fits to the functional forms in (48). Full fit results are shown in Table II.</figcaption></figure>

take it from the scaling of the median of the distribution $P(C_4)$ with $\xi$. As usual, we first fitted the median as a function of $C_4$, see Fig. 9, and then translated it to a power law in $\xi$.

Fig. 12 shows our estimations and fits for $\tau(q)$ for $T = 0.9$ and $T = T_c$. Although the functional behaviors of both curves are similar, the clear difference reveals the significantly different nature of the multifractal behaviors at the critical points and in the spin-glass phase.

### Summary of the statistical information fits

As stated in the main text, and in the above sections, in order to compute and interpolate $\tau(q)$, and its errors [see 44], we need to make a huge number of fits. For example, to obtain the data shown in Fig. 9, we need to extrapolate to $N_R \to \infty$ the median of $C_4$, see Fig. 8, 16 times at each $\xi(t_w)$ in order to estimate errors. This process would be impossible without some automation. To achieve this, we use a nonlinear least-squares method using the Levenberg-Marquardt solver [49] routine available at the GNU Scientific Library. However, in order to calibrate the process and obtain a fair fits, we perform manually all the fits shown in the present work. The complete statistical information of the fits is summarized in Table II.

Fig. 13, which shows the evolution of the 10th and 20th moments of $C_4$ as function of the coherence length, includes a fit to a power law $a_q \xi^{-\tau(q)}$, where the only parameter is the amplitude constant, $a_q$, and the respective value of $\tau(q)$ is the result of our analysis. These fits, included in Table II, are a good confirmation of our analysis.
Interpolations necessary to Fig. 4 in the main text

Finally, let us recall that Fig. 4 in the main text requires the evaluation of $C_4^{av}(r, t_w)$ at non-integer values of its argument. Here, we shall address only the interpolation for our data for the spin glass at $T = 0.9$ and $t_w$ such that $\xi(t_w) = 20$.

In order to interpolate our data, obtained for integer $r$, we have performed a fit that should account for both the short- and long-distance behavior of the correlations:

$$ 
\begin{aligned}
C_4^{av}(r, t_w) &= F_{sd}(r) + F_{ld}(r), \\
F_{sd}(r) &= (b_0 + b_1 r + b_2 r^2)e^{-(r/2.6)^4}, \quad (49) \\
F_{ld}(r) &= \frac{A}{r^\alpha}e^{-(r/\xi_{exp})^\beta}.
\end{aligned}
$$

While the functional form for $F_{ld}(r)$ is well known and physically motivated [26, 45, 46], we regard $F_{sd}(r)$ as a purely ad-hoc fitting device. The above functional form fits our integer-$r$ data (within errors) for all $r \ge 1$. Of course, the fitting function can be evaluated whether the argument is integer or not. The fit's figure of merit is $\chi^2/\text{dof} = 9.18/71$ (dof stands for 'degrees of freedom'), where we have considered only the diagonal part of the covariance matrix. This partly explains the exceedingly small value of $\chi^2$ that we obtained in the fit (the departure of the $p$-value from one is $\sim 10^{-18}$). In order to compute $\chi^2$, we considered as well the first image at $L-r$ [50], i.e., we compared the numerical data to $f_{sd}(r) + f_{sd}(L-r) + f_{ld}(r) + f_{ld}(L-r)$.

We conclude this paragraph with the fit parameters (we report many digits for the sake of reproducibility). For the short-distance piece we have:

$$ 
\begin{aligned}
b_0 &= 0.0453360389027432, \\
b_1 &= -0.0264096080489446, \quad (50) \\
b_2 &= 0.00558776268187863.
\end{aligned}
$$

The parameters of the long-term decay are:

$$ 
\begin{aligned}
\alpha &= 0.45829, \quad \beta = 1.41217, \quad (51) \\
A &= 0.551495, \quad \xi_{exp} = 20.5207.
\end{aligned}
$$

Comparisons between the diluted ferromagnet and the spin glass

As explained in the main text, both the diluted ferromagnetic Ising model and the spin glass display domain-growth dynamics in their low-temperature phases. In both cases, the size of the domains can be characterized through a coherence length $\xi(t_w)$, see Fig. 14. However, while the spin-glass correlation function at distances $r = \xi$ shows strong statistical fluctuations, see Fig. 15—bottom, this is not the case for the diluted ferromagnet Fig. 15—top. Indeed, for the diluted ferromagnet we see

<figure><img src="image_14.png" /><figcaption>Figure 13. The 10th and 20th moments of $C_4$ as functions of the coherence length $\xi$ at $T = 0.9$. Dashed lines are fits to $a\xi^{-\tau(q)}$ where the only fitting parameter is the proportionality constant $a$ (fit details are available in Table II). The values of $\tau(q = 10)$, and $\tau(q = 20)$ are our final estimates presented in Fig. 12. These fits can be considered as a test of our analysis, showing that our final estimate of $\tau(q)$ is accurate.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14. Correlation function $C_4^{av}(r, t_w)$ multiplied by $r^2$, as computed for the three-dimensional Ising diluted ferromagnet (DIL) and for the Ising spin glass (EA), versus distance $r$. Data were obtained in systems of linear size $L = 160$ with coherence length $\xi(t_w) = 20$ (dashed vertical line) at temperature $T = 0.9$—recall that $T_c \approx 1.1$ for the spin glass [31]. As explained in Methods in the main text, the coherence length is computed from the integral $I_2 = \int_0^\infty r^2 C_4^{av}(r, t_w) dr$. For both models, $r^2 C_4^{av}(r)$ peaks near $\xi(t_w)$.</figcaption></figure>

from the figure that $\overline{C_4^2}/\overline{C_4^2}^2$ quickly reaches a finite limit as $\xi(t_w)$ grows. Instead, a steady increase of the ratio $\overline{C_4^2}/\overline{C_4^2}^2$ is observed for the spin glass.

[1] U. Frisch and G. Parisi, in *Turbulence and predictability in geophysical fluid dynamics and climate dynamics* (1983 International School of Physics "Enrico Fermi", Varenna), edited by M. Ghil, R. Benzi, and G. Parisi (North-Holland, Amsterdam, 1985).
Interpolations necessary to Fig. 4 in the main text

Finally, let us recall that Fig. 4 in the main text requires the evaluation of $C_4^{av}(r, t_w)$ at non-integer values of its argument. Here, we shall address only the interpolation for our data for the spin glass at $T = 0.9$ and $t_w$ such that $\xi(t_w) = 20$.

In order to interpolate our data, obtained for integer $r$, we have performed a fit that should account for both the short- and long-distance behavior of the correlations:

$$ 
\begin{align} 
C_4^{av}(r, t_w) &= F_{sd}(r) + F_{ld}(r), \\ 
F_{sd}(r) &= (b_0 + b_1 r + b_2 r^2)e^{-(r/2.6)^4}, \tag{49} \\ 
F_{ld}(r) &= \frac{A}{r^\alpha}e^{-(r/\xi_{exp})^\beta}. 
\end{align} 
$$

While the functional form for $F_{ld}(r)$ is well known and physically motivated [26, 45, 46], we regard $F_{sd}(r)$ as a purely ad-hoc fitting device. The above functional form fits our integer-$r$ data (within errors) for all $r \ge 1$. Of course, the fitting function can be evaluated whether the argument is integer or not. The fit's figure of merit is $\chi^2/\text{dof} = 9.18/71$ (dof stands for 'degrees of freedom'), where we have considered only the diagonal part of the covariance matrix. This partly explains the exceedingly small value of $\chi^2$ that we obtained in the fit (the departure of the $p$-value from one is $\sim 10^{-18}$). In order to compute $\chi^2$, we considered as well the first image at $L-r$ [50], i.e., we compared the numerical data to $f_{sd}(r) + f_{sd}(L-r) + f_{ld}(r) + f_{ld}(L-r)$.

We conclude this paragraph with the fit parameters (we report many digits for the sake of reproducibility). For the short-distance piece we have:

$$ 
\begin{align} 
b_0 &= 0.0453360389027432, \\ 
b_1 &= -0.0264096080489446, \tag{50} \\ 
b_2 &= 0.00558776268187863. 
\end{align} 
$$

The parameters of the long-term decay are:

$$ 
\begin{align} 
\alpha &= 0.45829, \quad \beta = 1.41217, \tag{51} \\ 
A &= 0.551495, \quad \xi_{exp} = 20.5207. 
\end{align} 
$$

Comparisons between the diluted ferromagnet and the spin glass

As explained in the main text, both the diluted ferromagnetic Ising model and the spin glass display domain-growth dynamics in their low-temperature phases. In both cases, the size of the domains can be characterized through a coherence length $\xi(t_w)$, see Fig. 14. However, while the spin-glass correlation function at distances $r = \xi$ shows strong statistical fluctuations, see Fig. 15—bottom, this is not the case for the diluted ferromagnet Fig. 15—top. Indeed, for the diluted ferromagnet we see

<figure><img src="image_14.png" /><figcaption>Figure 13. The 10th and 20th moments of $C_4$ as functions of the coherence length $\xi$ at $T = 0.9$. Dashed lines are fits to $a\xi^{-\tau(q)}$ where the only fitting parameter is the proportionality constant $a$ (fit details are available in Table II). The values of $\tau(q = 10)$, and $\tau(q = 20)$ are our final estimates presented in Fig. 12. These fits can be considered as a test of our analysis, showing that our final estimate of $\tau(q)$ is accurate.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14. Correlation function $C_4^{av}(r, t_w)$ multiplied by $r^2$, as computed for the three-dimensional Ising diluted ferromagnet (DIL) and for the Ising spin glass (EA), versus distance $r$. Data were obtained in systems of linear size $L = 160$ with coherence length $\xi(t_w) = 20$ (dashed vertical line) at temperature $T = 0.9$—recall that $T_c \approx 1.1$ for the spin glass [31]. As explained in Methods in the main text, the coherence length is computed from the integral $I_2 = \int_0^\infty r^2 C_4^{av}(r, t_w) dr$. For both models, $r^2 C_4^{av}(r)$ peaks near $\xi(t_w)$.</figcaption></figure>

from the figure that $\overline{C_4^2}/\overline{C_4^2}^2$ quickly reaches a finite limit as $\xi(t_w)$ grows. Instead, a steady increase of the ratio $\overline{C_4^2}/\overline{C_4^2}^2$ is observed for the spin glass.

[1] U. Frisch and G. Parisi, in *Turbulence and predictability in geophysical fluid dynamics and climate dynamics* (1983 International School of Physics "Enrico Fermi", Varenna), edited by M. Ghil, R. Benzi, and G. Parisi (North-Holland, Amsterdam, 1985).
<table><caption>Table II. Summary of the statistical information of all the fits reported in our work. The fits have been done with GNUplot. We have used this data in order to calibrate our analysis, which has been automated to estimate the errors following the strategy of [44].</caption><thead><tr><th>Identifier</th><th>Functional Form</th><th>Fitting Range</th><th>χ<sup>2</sup>/dof</th><th>Parameters</th></tr></thead><tbody><tr><td>[C<sub>4</sub>], Fig. 8</td><td rowspan="2">f(x) = a + bx + cx<sup>2</sup></td><td rowspan="2">N<sub>R</sub> ∈ [32, 512]</td><td>0.872768/2</td><td>a = 0.01492(12), b = -0.52(2), c = -5.0(6)</td></tr><tr><td>[S<sub>x</sub>S<sub>x+r</sub>]<sup>2</sup>, Fig. 8</td><td>0.946227/2</td><td>a = 0.01492(12), b = 0.46(2), c = -4.3(6)</td></tr><tr><td>I(N<sub>R</sub>), Fig. 10</td><td>f(x) = ax<sup>(b+1)</sup></td><td></td><td>3.70448/3</td><td>a = 14(2), b = 2.02(4)</td></tr><tr><td>B(ξ), Fig. 11</td><td>f(x) = a + b log x</td><td></td><td>5.95738/9</td><td>a = 0.78(11), b = 0.39(4)</td></tr><tr><td>C̄<sub>4</sub><sup>q=10</sup>, Fig. 13</td><td rowspan="2">f(x) = ax<sup>-τ(q)</sup></td><td rowspan="2">ξ ∈ [10, 20]</td><td>5.6406/10</td><td>a = 0.00513(2)</td></tr><tr><td>C̄<sub>4</sub><sup>q=20</sup>, Fig. 13</td><td>7.64613/10</td><td>a = 0.001337(11)</td></tr><tr><td>Median of C<sub>4</sub>, Fig. 9</td><td rowspan="3">f(x) = ax<sup>b</sup></td><td>C̄<sub>4</sub>(ξ ∈ [4, 20])</td><td>1.86484/15</td><td>a = 1.489(16), b = 1.575(4)</td></tr><tr><td>median[C<sub>4</sub>]/C̄<sub>4</sub>, Fig. 5</td><td>C̄<sub>4</sub>(ξ ∈ [4, 20])</td><td>8.88094/15</td><td>a = 1.488(7), b = 0.575(2)</td></tr><tr><td>C̄<sub>4</sub><sup>av</sup>(r = ξ; T = 0.9), Fig. 2</td><td>ξ ∈ [10, 20]</td><td>1.99859/9</td><td>a = 0.182(2), b = -0.411(5)</td></tr><tr><td>τ<sub>1</sub>(q) for T = 0.9, Fig. 5 and Fig. 12</td><td rowspan="2">f(x) = mx<sup>1+ax</sup>/(1+bx)</td><td rowspan="4">q ∈ [1, 20]</td><td>5.01227/18</td><td>a = 0.0129(18), b = 0.547(12)</td></tr><tr><td>τ<sub>1</sub>(q) for T = T<sub>c</sub>, Fig. 12</td><td>5.37753/18</td><td>a = 0.1718(17), b = 0.574(10)</td></tr><tr><td>τ<sub>2</sub>(q) for T = 0.9, Fig. 5 and Fig. 12</td><td rowspan="2">f(x) = mx<sup>1+ax</sup>/(1+bx)<sup>2</sup></td><td>9.33556/18</td><td>a = 0.039(2), b = 0.232(5)</td></tr><tr><td>τ<sub>2</sub>(q) for T = T<sub>c</sub>, Fig. 12</td><td>8.23336/18</td><td>a = 0.0446(2), b = 0.242(4)</td></tr></tbody></table>

[2] D. Harte, *Multifractals. Theory and applications*, 1st ed. (Chapman and Hall/CRC, New York, 2001).

[3] K. G. Wilson, *Scientific American* **241**, 158 (1979).

[4] G. Parisi, *Statistical Field Theory* (Addison-Wesley, 1988).

[5] M. Barnsley, *Fractals Everywhere*, 3rd ed. (Dover Publications Inc., Mineola, New York, 2012).

[6] R. Benzi, G. Paladin, G. Parisi, and A. Vulpiani, *Journal of Physics A: Mathematical and General* **17**, 3521 (1984).

[7] C. Castellani and L. Peliti, *Journal of physics A: mathematical and general* **19**, L429 (1986).

[8] H. E. Stanley and P. Meakin, *Nature* **335**, 405 (1988).

[9] T. C. Halsey, M. H. Jensen, L. P. Kadanoff, I. Procaccia, and B. I. Shraiman, *Phys. Rev. A* **33**, 1141 (1986).

[10] T. C. Halsey, M. H. Jensen, L. P. Kadanoff, I. Procaccia, and B. I. Shraiman, *Phys. Rev. A* **34**, 1601 (1986).

[11] A.-L. Barabási and H. E. Stanley, *Fractal Concepts in Surface Growth* (Cambridge University Press, New York, 2009).

[12] P. C. Ivanov, L. A. N. Amaral, A. L. H. Goldberger, M. G. Shlomo Rosenblum, Z. R. Struzik, and H. E. Stanley, *Nature* **399**, 461 (1999).

[13] L. Seuront and H. E. Stanley, *Proceedings of the National Academy of Sciences* **111**, 2206 (2014).

[14] A. Klopper, *Nature Physics* **10**, 183 (2014).

[15] R. Deidda, *Water Resources Research* **36**, 1779 (2000).

[16] J. Alvarez-Ramirez, J. Alvarez, and E. Rodriguez, *Energy Economics* **30**, 2645 (2008).

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

[18] P. Charbonneau, E. Marinari, M. Mézard, G. Parisi, F. Ricci-Tersenghi, G. Sicuro, and F. Zamponi, *Spin Glass Theory and Far Beyond* (WORLD SCIENTIFIC, 2023) https://www.worldscientific.com/doi/pdf/10.1142/13341.

[19] 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).

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

[21] Q. Zhai, I. Paga, 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, 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).

[22] I. Paga, Q. Zhai, 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, 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, *Journal of Statistical Mechanics: Theory and Experiment* **2021**, 033301 (2021).

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

[24] F. Belletti, M. Cotallo, A. Cruz, L. A. Fernandez, A. Gordillo, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, A. Muñoz Sudupe, D. Navarro, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, D. Sciretti, A. Tarancon, R. Tripiccione, and J. L. Velasco (Janus Collaboration), *Comp. Phys. Comm.* **178**, 208 (2008), arXiv:0704.3573.
<table><caption>Table II. Summary of the statistical information of all the fits reported in our work. The fits have been done with GNUplot. We have used this data in order to calibrate our analysis, which has been automated to estimate the errors following the strategy of [44].</caption><thead><tr><th>Identifier</th><th>Functional Form</th><th>Fitting Range</th><th>χ<sup>2</sup>/dof</th><th>Parameters</th></tr></thead><tbody><tr><td>[C<sub>4</sub>], Fig. 8</td><td rowspan="2">f(x) = a + bx + cx<sup>2</sup></td><td rowspan="2">N<sub>R</sub> ∈ [32, 512]</td><td>0.872768/2</td><td>a = 0.01492(12), b = -0.52(2), c = -5.0(6)</td></tr><tr><td>[S<sub>x</sub>S<sub>x+r</sub>]<sup>2</sup>, Fig. 8</td><td>0.946227/2</td><td>a = 0.01492(12), b = 0.46(2), c = -4.3(6)</td></tr><tr><td>I(N<sub>R</sub>), Fig. 10</td><td>f(x) = ax<sup>(b+1)</sup></td><td></td><td>3.70448/3</td><td>a = 14(2), b = 2.02(4)</td></tr><tr><td>B(ξ), Fig. 11</td><td>f(x) = a + b log x</td><td></td><td>5.95738/9</td><td>a = 0.78(11), b = 0.39(4)</td></tr><tr><td>C̄<sub>4</sub><sup>q=10</sup>, Fig. 13</td><td rowspan="2">f(x) = ax<sup>-τ(q)</sup></td><td rowspan="2">ξ ∈ [10, 20]</td><td>5.6406/10</td><td>a = 0.00513(2)</td></tr><tr><td>C̄<sub>4</sub><sup>q=20</sup>, Fig. 13</td><td>7.64613/10</td><td>a = 0.001337(11)</td></tr><tr><td>Median of C<sub>4</sub>, Fig. 9</td><td rowspan="3">f(x) = ax<sup>b</sup></td><td>C̄<sub>4</sub>(ξ ∈ [4, 20])</td><td>1.86484/15</td><td>a = 1.489(16), b = 1.575(4)</td></tr><tr><td>median[C<sub>4</sub>]/C̄<sub>4</sub>, Fig. 5</td><td>C̄<sub>4</sub>(ξ ∈ [4, 20])</td><td>8.88094/15</td><td>a = 1.488(7), b = 0.575(2)</td></tr><tr><td>C̄<sub>4</sub><sup>av</sup>(r = ξ; T = 0.9), Fig. 2</td><td>ξ ∈ [10, 20]</td><td>1.99859/9</td><td>a = 0.182(2), b = -0.411(5)</td></tr><tr><td>τ<sub>1</sub>(q) for T = 0.9, Fig. 5 and Fig. 12</td><td rowspan="2">f(x) = mx<sup>1+ax</sup> / (1+bx)</td><td rowspan="4">q ∈ [1, 20]</td><td>5.01227/18</td><td>a = 0.0129(18), b = 0.547(12)</td></tr><tr><td>τ<sub>1</sub>(q) for T = T<sub>c</sub>, Fig. 12</td><td>5.37753/18</td><td>a = 0.1718(17), b = 0.574(10)</td></tr><tr><td>τ<sub>2</sub>(q) for T = 0.9, Fig. 5 and Fig. 12</td><td rowspan="2">f(x) = mx<sup>1+ax</sup> / ((1+bx)<sup>2</sup>)</td><td>9.33556/18</td><td>a = 0.039(2), b = 0.232(5)</td></tr><tr><td>τ<sub>2</sub>(q) for T = T<sub>c</sub>, Fig. 12</td><td>8.23336/18</td><td>a = 0.0446(2), b = 0.242(4)</td></tr></tbody></table>

[2] D. Harte, *Multifractals. Theory and applications*, 1st ed. (Chapman and Hall/CRC, New York, 2001).

[3] K. G. Wilson, *Scientific American* **241**, 158 (1979).

[4] G. Parisi, *Statistical Field Theory* (Addison-Wesley, 1988).

[5] M. Barnsley, *Fractals Everywhere*, 3rd ed. (Dover Publications Inc., Mineola, New York, 2012).

[6] R. Benzi, G. Paladin, G. Parisi, and A. Vulpiani, *Journal of Physics A: Mathematical and General* **17**, 3521 (1984).

[7] C. Castellani and L. Peliti, *Journal of physics A: mathematical and general* **19**, L429 (1986).

[8] H. E. Stanley and P. Meakin, *Nature* **335**, 405 (1988).

[9] T. C. Halsey, M. H. Jensen, L. P. Kadanoff, I. Procaccia, and B. I. Shraiman, *Phys. Rev. A* **33**, 1141 (1986).

[10] T. C. Halsey, M. H. Jensen, L. P. Kadanoff, I. Procaccia, and B. I. Shraiman, *Phys. Rev. A* **34**, 1601 (1986).

[11] A.-L. Barabási and H. E. Stanley, *Fractal Concepts in Surface Growth* (Cambridge University Press, New York, 2009).

[12] P. C. Ivanov, L. A. N. Amaral, A. L. H. Goldberger, M. G. Shlomo Rosenblum, Z. R. Struzik, and H. E. Stanley, *Nature* **399**, 461 (1999).

[13] L. Seuront and H. E. Stanley, *Proceedings of the National Academy of Sciences* **111**, 2206 (2014).

[14] A. Klopper, *Nature Physics* **10**, 183 (2014).

[15] R. Deidda, *Water Resources Research* **36**, 1779 (2000).

[16] J. Alvarez-Ramirez, J. Alvarez, and E. Rodriguez, *Energy Economics* **30**, 2645 (2008).

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

[18] P. Charbonneau, E. Marinari, M. Mézard, G. Parisi, F. Ricci-Tersenghi, G. Sicuro, and F. Zamponi, *Spin Glass Theory and Far Beyond* (WORLD SCIENTIFIC, 2023) https://www.worldscientific.com/doi/pdf/10.1142/13341.

[19] 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).

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

[21] Q. Zhai, I. Paga, 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, 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).

[22] I. Paga, Q. Zhai, 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, 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, *Journal of Statistical Mechanics: Theory and Experiment* **2021**, 033301 (2021).

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

[24] F. Belletti, M. Cotallo, A. Cruz, L. A. Fernandez, A. Gordillo, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, A. Muñoz Sudupe, D. Navarro, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, D. Sciretti, A. Tarancon, R. Tripiccione, and J. L. Velasco (Janus Collaboration), *Comp. Phys. Comm.* **178**, 208 (2008), arXiv:0704.3573.
<figure><img src="image_16.png" /><figcaption>Figure 15. Ratio of the second moment of the spin-glass correlation function $C_4$ computed at $r = \xi(t_w)$, $\overline{C_4^2}$, to the squared first moment, $\overline{C_4^2}$, as a function of the coherence length $\xi(t_w)$. Data computed for the Ising diluted ferromagnet (bottom) and for Ising spin glass (top) at temperatures $T, \tilde{T} = 0.9, 0.8$ and $0.7$ (remember that $T^{\text{DIL}} = \tilde{T} (T_c^{\text{DIL}}/T_c)$). Error bars are smaller than the point size. Note that, while the diluted ferromagnet tends to one for large coherence length, the spin glass follows a power law, which indicates that in the scaling limit (i.e., $\xi(t_w) \to \infty$) the order of magnitude of $\overline{C_4^2}$ is larger than that of $\overline{C_4^2}$.</figcaption></figure>

[25] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvión, 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. Tarancón, R. Tripiccione, and D. Yllanes, Proceedings of the National Academy of Sciences **114**, 1838 (2017).

[26] 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).

[27] 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, Nature Physics **19**, 978 (2023).

[28] K. Jonason, E. Vincent, J. Hammann, J. P. Bouchaud, and P. Nordblad, Phys. Rev. Lett. **81**, 3243 (1998).

[29] 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, R. Tripiccione, and D. Yllanes, Communications physics **4**, 74 (2021), https://www.nature.com/articles/s42005-021-00565-9.pdf.

[30] 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.

[31] 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), Phys. Rev. B **88**, 224416 (2013), arXiv:1310.2910.

[32] At its critical temperature, the two-dimensional diluted ferromagnetic Potts model with more than two states presents multiscaling as well [51] —this is also the case for the diluted Ising model in $D=3$ [52, 53].

[33] The droplet picture of spin glasses [54–56] predicts $\theta = 0$, similarly to the ferromagnet. Neither simulations nor experimental data are compatible with $\theta=0$, unless one is willing to accept that the available range of $\xi(t_w)$ is too small to display the true asymptotic behavior [26].

[34] R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Masaioli, and S. Succi, Phys. Rev. E **48**, R29 (1993).

[35] J.-H. Chen and T. C. Lubensky, Phys. Rev. B **16**, 2106 (1977).

[36] C. De Dominicis and I. Kondor, Journal of Physics A: Mathematical and General **22**, L743 (1989).

[37] The correlation function behaves as $C_4^{\text{av}}(r, t_w) \sim G(r/\xi(t_w))/r^\theta$ for large $r$, where the cut-off function $G(x)$ decays faster than exponentially as $x$ grows [see, e.g., Refs. [45, 50]]. Hence, for $r \sim \xi(t_w)$ one may consider either power-law scaling in $r$ —as in (1)— or in $\xi(t_w)$ —as in (2). The analysis of scale invariance in a fractal (or multifractal) geometry typically involves power laws.

[38] M. Baity-Jesi, R. A. Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iniguez, A. Maiorano, M. F., 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, J. Stat. Mech. **2014**, P05014 (2014), arXiv:1403.2622.
<figure><img src="image_16.png" /><figcaption>Figure 15. Ratio of the second moment of the spin-glass correlation function $C_4$ computed at $r = \xi(t_w)$, $\overline{C_4^2}$, to the squared first moment, $\overline{C_4^2}$, as a function of the coherence length $\xi(t_w)$. Data computed for the Ising diluted ferromagnet (**top**) and for Ising spin glass (**bottom**) at temperatures $T$, $\tilde{T} = 0.9, 0.8$ and $0.7$ (remember that $T^{\text{DIL}} = \tilde{T} (T_c^{\text{DIL}}/T_c)$). Error bars are smaller than the point size. Note that, while the diluted ferromagnet tends to one for large coherence length, the spin glass follows a power law, which indicates that in the scaling limit (i.e., $\xi(t_w) \to \infty$) the order of magnitude of $\overline{C_4^2}$ is larger than that of $\overline{C_4^2}$.</figcaption></figure>

[25] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvión, 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. Tarancón, R. Tripiccione, and D. Yllanes, Proceedings of the National Academy of Sciences **114**, 1838 (2017).

[26] 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).

[27] 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, Nature Physics **19**, 978 (2023).

[28] K. Jonason, E. Vincent, J. Hammann, J. P. Bouchaud, and P. Nordblad, Phys. Rev. Lett. **81**, 3243 (1998).

[29] 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, R. Tripiccione, and D. Yllanes, Communications physics **4**, 74 (2021), https://www.nature.com/articles/s42005-021-00565-9.pdf.

[30] 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.

[31] 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), Phys. Rev. B **88**, 224416 (2013), arXiv:1310.2910.

[32] At its critical temperature, the two-dimensional diluted ferromagnetic Potts model with more than two states presents multiscaling as well [51] —this is also the case for the diluted Ising model in $D=3$ [52, 53].

[33] The droplet picture of spin glasses [54–56] predicts $\theta=0$, similarly to the ferromagnet. Neither simulations nor experimental data are compatible with $\theta=0$, unless one is willing to accept that the available range of $\xi(t_w)$ is too small to display the true asymptotic behavior [26].

[34] R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Masaioli, and S. Succi, Phys. Rev. E **48**, R29 (1993).

[35] J.-H. Chen and T. C. Lubensky, Phys. Rev. B **16**, 2106 (1977).

[36] C. De Dominicis and I. Kondor, Journal of Physics A: Mathematical and General **22**, L743 (1989).

[37] The correlation function behaves as $C_4^{\text{av}}(r, t_w) \sim G(r/\xi(t_w))/r^\theta$ for large $r$, where the cut-off function $G(x)$ decays faster than exponentially as $x$ grows [see, e.g., Refs. [45, 50]]. Hence, for $r \sim \xi(t_w)$ one may consider either power-law scaling in $r$ —as in (1)— or in $\xi(t_w)$ —as in (2). The analysis of scale invariance in a fractal (or multifractal) geometry typically involves power laws.

[38] M. Baity-Jesi, R. A. Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. Iniguez, A. Maiorano, M. F., 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, J. Stat. Mech. **2014**, P05014 (2014), arXiv:1403.2622.
[39] S. Franz, M. Mézard, G. Parisi, and L. Peliti, Phys. Rev. Lett. **81**, 1758 (1998).

[40] R. Alvarez Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Phys. Rev. Lett. **105**, 177202 (2010), arXiv:1003.2943.

[41] M. Wittmann and A. P. Young, Journal of Statistical Mechanics: Theory and Experiment **2016**, 013301 (2016).

[42] S. Jensen, N. Read, and A. P. Young, Phys. Rev. E **104**, 034105 (2021).

[43] D. J. Amit and V. Martín-Mayor, *Field Theory, the Renormalization Group and Critical Phenomena*, 3rd ed. (World Scientific, Singapore, 2005).

[44] D. Yllanes, *Rugged Free-Energy Landscapes in Disordered Spin Systems*, Ph.D. thesis, Universidad Complutense de Madrid (2011), arXiv:1111.0266.

[45] 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.

[46] 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), arXiv:0811.2864.

[47] P.-E. Berche, C. Chatelain, B. Berche, and W. Janke, The European Physical Journal B - Condensed Matter and Complex Systems **38**, 463 (2004).

[48] Recall that we set times such that $r = \xi(t_w)$.

[49] K. Levenberg, Quarterly of Applied Mathematics **2**, 164 (1944).

[50] L. A. Fernández, E. Marinari, V. Martín-Mayor, G. Parisi, and J. Ruiz-Lorenzo, Journal of Physics A: Mathematical and Theoretical **52**, 224002 (2019).

[51] A. W. W. Ludwig, Nucl. Phys. B **330**, 639 (1990).

[52] T. Davis and J. Cardy, Nuclear Physics B **570**, 713 (2000).

[53] E. Marinari, V. Martín-Mayor, G. Parisi, F. Ricci-Tersenghi, and J. J. Ruiz-Lorenzo, “Multiscaling in the critical site-diluted ferromagnet,” (2023), manuscript in preparation.

[54] W. L. McMillan, Phys. Rev. B **28**, 5216 (1983).

[55] A. J. Bray and M. A. Moore, Phys. Rev. Lett. **41**, 1068 (1978).

[56] D. S. Fisher and D. A. Huse, Phys. Rev. Lett. **56**, 1601 (1986).
[39] S. Franz, M. Mézard, G. Parisi, and L. Peliti, Phys. Rev. Lett. **81**, 1758 (1998).

[40] R. Alvarez Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martín-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yllanes (Janus Collaboration), Phys. Rev. Lett. **105**, 177202 (2010), arXiv:1003.2943.

[41] M. Wittmann and A. P. Young, Journal of Statistical Mechanics: Theory and Experiment **2016**, 013301 (2016).

[42] S. Jensen, N. Read, and A. P. Young, Phys. Rev. E **104**, 034105 (2021).

[43] D. J. Amit and V. Martín-Mayor, *Field Theory, the Renormalization Group and Critical Phenomena*, 3rd ed. (World Scientific, Singapore, 2005).

[44] D. Yllanes, *Rugged Free-Energy Landscapes in Disordered Spin Systems*, Ph.D. thesis, Universidad Complutense de Madrid (2011), arXiv:1111.0266.

[45] 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.

[46] 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), arXiv:0811.2864.

[47] P.-E. Berche, C. Chatelain, B. Berche, and W. Janke, The European Physical Journal B - Condensed Matter and Complex Systems **38**, 463 (2004).

[48] Recall that we set times such that $r = \xi(t_w)$.

[49] K. Levenberg, Quarterly of Applied Mathematics **2**, 164 (1944).

[50] L. A. Fernandez, E. Marinari, V. Martín-Mayor, G. Parisi, and J. Ruiz-Lorenzo, Journal of Physics A: Mathematical and Theoretical **52**, 224002 (2019).

[51] A. W. W. Ludwig, Nucl. Phys. B **330**, 639 (1990).

[52] T. Davis and J. Cardy, Nuclear Physics B **570**, 713 (2000).

[53] E. Marinari, V. Martín-Mayor, G. Parisi, F. Ricci-Tersenghi, and J. J. Ruiz-Lorenzo, “Multiscaling in the critical site-diluted ferromagnet,” (2023), manuscript in preparation.

[54] W. L. McMillan, Phys. Rev. B **28**, 5216 (1983).

[55] A. J. Bray and M. A. Moore, Phys. Rev. Lett. **41**, 1068 (1978).

[56] D. S. Fisher and D. A. Huse, Phys. Rev. Lett. **56**, 1601 (1986).