The simplest spin glass revisited: finite-size effects of the energy landscape can modify

aging dynamics in the thermodynamic limit

Bin Li1, 2 and Yuliang Jin1, 2, 3, ∗
1
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
2
School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
3
Center for Theoretical Interdisciplinary Sciences, Wenzhou Institute,
University of Chinese Academy of Sciences, Wenzhou, Zhejiang 325001, China
(Dated: January 3, 2025)
The random energy model is one of the few glass models whose asymptotic activated aging dy-
namics are solvable. However, the existing aging theory, i.e., Bouchaud’s trap model, does not
arXiv:2501.00338v1 [cond-mat.dis-nn] 31 Dec 2024

agree with dynamical simulation results obtained in finite-sized systems. Here we show that this
discrepancy originates from non-negligible finite-size corrections in the energy barrier distributions.
The finite-size effects add a logarithmic decay term in the time-correlation aging function, which
destroys the asymptotic large-time plateau predicted by Bouchaud’s trap model in the spin glass
phase. Surprisingly, the finite-size effects also give corrections, preserved even in the thermodynamic
limit, to the value of the asymptotic plateau. It results in an unexpected dynamical transition where
weak ergodicity breaking occurs, at a temperature Td above the thermodynamic spin-glass transition
temperature Tc . Based on the barrier distributions obtained by a numerical barrier-tree method
and an expansion theory, we propose a generalized trap model to incorporate such finite-size effects.
The theoretically derived aging behavior of the generalized trap model explains the Monte-Carlo dy-
namical simulation data of random energy models with Gaussian and exponential random energies.
Our results suggest that the double limits of large system size and long time are not interchangeable
for the activated aging dynamics.

Introduction. Glassy systems, including spin which was initially proposed to describe aging dynam-
glasses [1], structural glasses [2], polymers [3], col- ics in the simplest spin glass model [21], i.e., the ran-
loidal suspensions [4], granular materials [5], active mat- dom energy model (REM) [22]. A particularly interesting
ter [6] and artificial neural networks [7], are characterized outcome from the BTM is the so-called weak ergodicity
by multiple thermodynamically metastable states (glass breaking in the spin glass phase (at a temperature T be-
basins) that are separated by energy barriers in the en- low the spin glass transition temperature Tc ), where the
ergy landscape [8], and non-equilibrium aging dynam- aging function (see the definition below) converges, in the
ics [9–11]. Aging refers to the dynamical slowing down large-time limit, to a plateau whose value is determined
of relaxation processes with an increasing “age” quanti- by the so-called arcsin law. Such asymptotic solutions of
fied by the waiting time (or the aging time) tw elapsed the BTM are supported by recent rigorous mathematical
after the system is quenched. derivations [23–26]. However, this result is inconsistent
Two mechanisms of aging have been proposed previ- with Monte-Carlo (MC) simulations of the REM (with
ously. (i) The first type, called mean-field aging, is stud- Gaussian random energies) consisting of a finite number
ied by a set of closed equations for the two-time corre- N of spins [27, 28]: the obtained aging functions do not
lation and response functions in mean-field glass mod- show any sign of convergence to the predicted plateau.
els [12–15]. Mean-field aging is related to slow descent in
This study considers the activated aging dynamics in
the (free) energy landscape after quenching, during which
the REMs with Gaussian and exponential energy distri-
paths to find lower energies become more and more scarce
butions
√ (GREM and EREM) [22, 27, 29], with Tc =
– called an “entropic effect” [9, 16–18]. Because energy
1/ 2 log 2 ≈ 0.849 [22] and Tc = 1 respectively. We
barriers are infinite in mean-field models (in the thermo-
show that the activated aging dynamics obtained by MC
dynamic limit), activated barrier-crossing processes do
simulations of the REMs can be theoretically explained
not occur in mean-field aging.
by a generalized trap model (GTM), which takes into
(ii) The second type, called activated aging, corre- account the finite-size corrections in the barrier energy
sponds to activated barrier-crossing processes. With an distribution. The aging dynamics are characterized by
increasing tw , the probability of the system being trapped the function Π(tw , tw + t) that describes the probabil-
in deeper basins increases, which results in a longer hop- ity of not leaving the glass basin between two times tw
ping time to escape the basin and thus slower relaxation and tw +t, starting from a random initial configuration at
dynamics. The fundamental theoretical model for acti- t = 0. The Π(tw , tw +t) is related to the (aging part) spin
vated aging is Bouchaud’s trap model (BTM) [19, 20], auto-correlation function C(tw , tw +t) = ⟨S(tw )S(tw +t)⟩
via C(tw , tw + t) = qEA Π(tw , tw + t), where qEA is the
Edwards-Anderson order parameter [17]. Further details
on the models and methods are provided in Appendixes
∗ yuliangjin@mail.itp.ac.cn A-C.
 2

1 1
a N =8
where τ0 ∼ O(1) is a microscopic time scale, α = 1 − x,
wx 1
0.8 N =8 0.8 N = 12 C = 1+x and A = (1−x)Γ(1−x)Γ 2 (x) , with Γ(x) a gamma
N = 12 N = 16
function. The aging behavior is dominated by the large-
Πw (tw )

Πw (tw )
0.6 0.6 N = 20
N = 16
0.4
N = 20 0.4 time regime tw > τ0 , which shows power-law conver-
gence to the asymptotic plateau, similar to the early β-
0.2 0.2
b relaxation scaling in the mode-coupling theory [33, 34].
0
100 102 104 106 108 1010 10−1 100 101 102 103 104 105 106 However, the arcsin law Eq. (2) is challenged by the
tw tw
simulation results of the finite-sized GREM obtained by
MC simulations, as shown in [27]. Up to the largest sys-
FIG. 1. Aging functions for the (a) GREM and (b)
EREM. In both models, w = 0.5 and T = 0.75. Open points
tem size and time (N ∼ 20 and tw ∼ 1010 ) that can be
are MC simulation data. The dotted line represents the BTM simulated by regular CPUs, there is no sign that the sim-
plateau (a) Hx,w ≈ 0.136 with x = T /Tc ≈ 0.883 (Tc ≈ ulation data of Πw (tw ) would converge to the predicted
0.849) and (b) Hx,w ≈ 0.303 with x = 0.75 (Tc = 1). The plateau Hx,w (see Fig. 1a). In contrast, such convergence
purple line represents the short-time behavior Πsw (tw ) = 1 − is well observed in the EREM [27] (see Fig. 1b). Thus
Ctw with (a) C ≈ 0.099 and (b) C ≈ 0.214. In (a), the the aging theory built on phenomenological BTM can-
dashed line represents the GTM plateau Hxeff ,w ≈ 0.854 with not fully explain the real dynamics of microscopic models
xeff ≈ 0.246; the MC data with t > 1 are fitted (lines) to (REMs).
ΠGTM
w (tw ) = Hxeff ,w (1 + A t−α B
w ) − N ln tw (see Eq. (7)) using
Hxeff ,w ≈ 0.854, α = 1 − xeff ≈ 0.754, A ≈ 0.068 from the
theory, and B = 0.702 from the fitting (the same B for all
N ); the brown line is the ΠGTM (tw ) with an inaccessibly large
1.0 a
w
N = 5000. 0.8

Bouchaud’s trap model (BTM). The BTM as- 0.6

Hxeff,
sumes an exponential distribution of the barrier energy T = 15Tc
∆E, 0.4 T = 185 Tc
T = 95 Tc
pBTM (∆E) = exp (−∆E/Tc ) , (1) 0.2 T = 65 Tc
T = Tc
in the spin glass phase T < Tc (the Boltzmann constant T = 34 Tc
kB = 1), independent of T . This exponential form is con- 0.0
0.0 0.25 0.5 0.75 1.0 1.25
sistent with the one-step replica breaking (1RSB) scheme W
of the organization of states in the spin glass phase of 1 1
the GREM [1]. The hopping time τ is related to ∆E b c
Hxeff ,w

Hxeff ,w
through the Arrhenius law, τ (∆E) ∼ exp(∆E/T ). It can
be shown that the Arrhenius law is consistent with the
simulation results of the REM (see Appendix D; the time 0 0
unit is N ). Based on Eq. (1) and the Arrhenius law, τ 0 1 2 3 4 0 0.5 1
T /Tc T /Tc
follows a power-law distribution, ψ(τ ) ∼ τ −(1+x) , where
x = T /Tc is the reduced temperature. Because 0 < x < 1
in the spin glass phase, the mean of τ diverges. Thus the FIG. 2. Phase diagram of the GTM. (a) The asymp-
totic plateau value of the aging function Hxeff ,w vs the weight
system takes infinite time to reach equilibrium, leading
parameter W , at different T (lines). Yellow and blue areas
to long-time aging. represent the WEB (W < 1) and frozen (W > 1) phases
Based on this model, Bouchaud showed an interest- respectively, with (inset) schematic energy landscapes. The
ing property called weak-ergodicity breaking (WEB) [12, EREM and GREM correspond to W = 0 and W ≈ 0.721 (tri-
19, 20]: for a finite, fixed tw , limt→∞ Π(tw , tw + t) = 0, angle). (b) In the GREM, the dynamical transition occurs at
but if the ratio t/tw = w is fixed, the aging function Td ≈ 3.05 ≈ 3.6Tc (star) and the thermodynamical transition
Πw (tw ) ≡ Π(tw , tw + t) converges asymptotically to a at Tc ≈ 0.849 (triangle). (c) In the EREM, the two transi-
non-zero constant, tions occur simultaneously at Td = Tc = 1.

sin(πx) ∞ du
Z
lim Πw (tw ) = ≡ Hx,w , (2) A generalized trap model (GTM). To explain the
π x
tw →∞ w u (1 + u) dynamical data in the finite-sized GREM, we propose a
which is called an arcsin law [30]. The arcsin law is proven generalized trap model. The GTM retains the assump-
to be rigorous in the asymptotic regime, taking N → tions and input of the original BTM (e.g., the Arrhenius
∞ first and then tw → ∞ [24, 31]. It can be shown law), except that the energy barrier distribution is not
theoretically that[32] (see SI Sec. S1), purely exponential anymore but a mixture of exponen-
( tial and Gaussian distributions,
BTM 1 − C tw , tw < τ0 , 
1 1

Πw (tw ) ∼ −α
(3) pGTM (∆E) = exp − ∆E − 2
(∆E − Ē) , (4)
Hx,w (1 + A tw ) , tw > τ0 , Tc 2N a
 3

0
--++++
-+--++
a effective xeff . The second effect introduces a logarithmic
+---++ BT σ 2
decay term ∼ N1 ln tw in the pre-asymptotic behavior in
-++--+ -+-+--
+----+ 40 BT Ē
--+++- -+++--
−2 TE σ 2 Eq. (3), which avoids the asymptotic plateau for a small
TE Ē

σ 2 , Ē
-++-++
N . Next we discuss these two effects in detail.
E

−4 ++++--
20 Expanding Eq. (4) gives  (neglecting  the constant
+--+--
+-+--+ ĒTc
----++
+----+
-++---
-+---+
1
term), ln pGTM (∆E) ∼ − Tc 1 − N a ∆E − 2N1 a ∆E 2 .
d The linear coefficient suggests an effective x parameter,
∆E
+----- +-+-++ 10 20 30 40
+--+--
 
+-+--+ N kTc
xeff = 1 − x = (1 − W )x, (5)
b e a
ln(p(∆E)) + ∆E/Tc
10−1 10
pGREM (∆E)

where x = T /Tc is the same as in the BTM. Note that

10−4
e−∆E/Tc
the correction −kTc /a is independent of N , and thus this
0
N =8 modification does not disappear even in the thermody-
10 −7 N = 16 namic limit. Correspondingly, the asymptotic plateau
pBT
GREM (∆E) N = 24
−10 Hxeff ,w is modified with x replaced by xeff in Eq. (2).
pGTM (∆E) N = 32
10 −10 In microscopic models (such as the GREM and
0 10 20 0 0.5 1 1.5 2
∆E ∆E/N EREM), the parameters k, a, and the transition temper-
ature Tc , are model-dependent. Their combination, the
c f
ln(p(∆E)) + ∆E/Tc

10−1 10 weight parameter W = kTc /a, weights the contributions

−∆E/Tc
e of linear and Gaussian terms in Eq. (4). The Hxeff ,w ver-
pEREM (∆E)

10−3 sus W phase diagram is presented in Fig. 2 (without loss

0
N =8 of generality, we set w = 1/2). Two phases can be iden-
−5
10 N = 16 tified. (i) A WEB phase with W < 1. For a given W ,
pBT
EREM (∆E) N = 32
10−7
−10
N = 128
the WEB occurs at a temperature T = Td > Tc , where,
pBTM (∆E)
0 5 10 15 0 0.2 0.4 0.6 0.8 1 Tc
∆E ∆E/N Td = , (6)
1−W

FIG. 3. Barrier energy distributions in the GREM and

such that xeff = (1 − W )Td /Tc = 1. To distinguish from
EREM. (a) A sub-tree formed by local minima and saddle the thermodynamic spin glass transition temperature Tc ,
points for the N = 6 GREM obtained by the barrier-tree we call Td a dynamical transition temperature. When
method, with corresponding spin configurations and energies W → 0, the original BTM is recovered. (ii) A frozen
(E-axis). The blue nodes are non-minimum-non-saddle con- phase with W > 1. Here Hxeff ,w = 1 independent of T ,
figurations in the basins. The barrier energy ∆E is the dif- suggesting completely frozen dynamics. It means that
ference between the local minimum energy and saddle point when the Gaussian term dominates (xeff ≤ 0) in Eq. (4),
energy. (b) Distribution of the barrier energy pBTGREM (∆E) by the effective temperature Teff = xeff Tc ≤ 0, and trap
the barrier-tree method (points) and pGTM (∆E) (Eq. 4) with escaping is impossible.
a = 1.12 and k = 0.96 (line), for the N = 20 GREM. (c) The second finite-size effect is caused by the quadratic
pBT
EREM (∆E) and pBTM (∆E) (Eq. 1) for the N = 20 EREM.
term ∼ ∆E 2 /N in Eq (4), which is O(1/N ). This term
(e-f) ln p(∆E) + ∆E/Tc for the (e) GREM and (f) EREM,
where open points are obtained by the barrier-tree algorithm adds a logarithmic decay term (higher-order corrections
(the data do not cover the full range due to limited statistics), are neglected) to the pre-asymptotic behavior in Eq. (3),
solid-filled points by the tree-expansion theory, and lines in as derived in SI Sec. S1D [32]:
(e) are fitting to the Gaussian function. The variance σ 2 and (
mean Ē obtained by the Gaussian fitting are plotted in (d) as GTM 1 − C tw , tw < τ0 ,
Πw (tw ) ∼ B (7)
functions of N . Linear fitting of the data in (d) gives a = 1.12 Hxeff ,w (1 + A t−α
w ) − N ln tw , tw > τ0 .
and k = 0.96 (σ 2 = aN and Ē = kN ).
where α, A, C have the same expressions as in Eq. (3)
T 2H′
xeff ,w
with x replaced by xeff , and B = − a . For small
where a is a model-dependent parameter, and Ē = kN N , the logarithmic decay dominates the overall behav-
is assumed to be extensive in general. Note that the ior of the aging function, which is the reason why the
Gaussian term follows the standard central limit theo- predicted plateau is not observable in GREM. Below we
rem. Equation (4) will be explicitly examined in REMs present such analyses.
below, but for now we take it as input to the GTM Aging dynamics of finite-sized GREM and
and derive theoretically the corresponding aging behavior EREM explained by the GTM. To explain the MC
(see details in SI Sec. S1 [32]). data in Fig. 1 using the GTM, firstly we need to deter-
There are two important finite-size effects due to the mine the parameters a and k in the barrier distribution
Gaussian term in Eq. (4). The first effect modifies the Eq. (4). This requires obtaining accurate barrier distri-
asymptotic plateau Hx,w in Eq. (2) with x replaced by an butions in finite-sized GREM and EREM. To do that,
 4

1 108 To understand the MC results for the GREM in

T = 0.75 b
0.8
T = 1.00 Fig. 1a, it is necessary to consider the GTM (with the
T = 1.50 106 106 EREM
T = 2.00 GREM
above-obtained parameters a = 1.12 and k = 0.96). Due
0.6 T = 2.50 to strong higher-order corrections in small-sized GREM,
Πw (tw )

R
104 103
T = 3.00 Eq. (7) does not agree exactly with the MC data if

R
T = 3.60
0.4 100 the theoretical coefficients A, B and C are used (see SI
102
10 15 20
0.2 N Fig. S4 [32]). Nevertheless, the MC data of different N
a 100 can be fitted to Eq. (7) with only one N -independent fit-
0 ting parameter B (we keep theoretical values of A and
10−1 101 103 105 107 109 1 2 3 4
Ctw T /Tc C). As can be seen from the data, in small-N GREM, the
1
T = 0.625 T = 0.75
aging dynamics are dominated by the logarithmic decay.
108 d
0.8
T = 0.875 T = 1.00 Unlike in the EREM, the plateau at Hxeff ,w in the GREM
T = 1.10 T = 1.20 N =8
T = 1.30 106 is only visible for much larger N , which is inaccessible by
0.6 N = 12 the current computational power (see Fig. 1a). Impor-
Πw (tw )

104
N = 16
R

N = 20 tantly, in the large-N limit, the GREM plateau should

0.4
converge to Hxeff ,w predicted by the GTM, rather than
102
0.2
Hx,w by the BTM. Thus the finite-size effects, uncon-
c 100 ventionally, remain in the thermodynamic limit. This is
0 −1 1 3 5 7 9 0.6 0.8 1 1.2 1.4
because the quadratic term in Eq. (4) introduces an O(1)
10 10 10 10 10 10
Ctw T /Tc correction, as discussed above (see Eq. 5).
An unexpected dynamical transition in the
FIG. 4. Dynamical transition in the GREM. (a) Πw (tw ) GREM. An interesting prediction from the GTM is the
are obtained for the GREM (N = 20, w = 0.5) by MC sim- dynamical transition at Td (Eq. 6) for the onset of WEB,
wxeff
ulations at rescaled times Ctw with C = 1+x eff
. The solid at which the asymptotic plateau Hxeff ,w vanishes. This
s
line is Πw (tw ) = 1 − Ctw (see Eq. 7). The dotted line is result is surprising, because it is conventionally believed
Π∗ = 0.2, and the two time scales ts,∗ ∗
w (circle) and tw (crosses) that, in REMs, the onset of WEB occurs simultane-
∗ s,∗
are marked. (b) R = tw /tw vs T for the GREM with differ-
ously with the thermodynamic spin glass transition at Tc ,
ent N . Td ≈ 3.6Tc is indicated (star). (c-d) Analogous data
for the EREM. (b-inset) R as a function of N at T = 1.17Tc . where the entropy density s vanishes [19]. According to
the GTM results (Fig. 2b,c), such a coincidence (Td = Tc )
only exists in the EREM (W = 0), but not in the GREM
(W > 0 and Td > Tc ). It should be emphasized the
we employ two independent approaches. The first one is dynamic transition at Td does not disappear even in the
an exhaustive enumeration of energy barriers using the thermodynamic limit, since the correction W = kTc /a
barrier-tree (BT) algorithm [35–38] (see Fig. 3a for an in Eq. (6) is independent of N . Using a = 1.12 and
example sub-barrier-tree and Appendix B for algorithm k = 0.96 obtained above, we derive Td ≈ 3.05 ≈ 3.6Tc in
details). The method provides the exact pBT (∆E), but is the GREM.
restricted to small systems (N ≤ 20). We also develop a
tree-expansion (TE) theory to compute analytically and Can we observe the dynamic transition in the MC sim-
approximately pTE (∆E), which works for any N (see SI ulations of the GREM? To show that, in Fig. 4a we
Sec. S2A [32]). The difference between pBT (∆E) and plot the MC data of Πw (tw ) for the GREM, with tw
wxeff
pTE (∆E) is unnoticeable, as shown in Fig 3b-f. rescaled by C = 1+x eff
to collapse the short-time be-
s
While the EREM distribution pEREM (∆E) (Fig. 3c) is havior Πw (tw ) = 1 − Ctw (note that C is not a fitting
very close to the exponential Eq. (1) as expected by the parameter). In the ergodic phase (T > Td ), the dynamics
original BTM, it is obvious that the GREM distribution should be dominated by this short-time behavior. The
pGREM (∆E) (Fig. 3b) deviates from a pure exponential deviation from it, which occurs at T < Td , characterizes
distribution. To reveal further this deviation, the ex- the WEB. For this reason, we define short-time and long-
∗ ∗
ponential part Eq. (1) is subtracted, and the difference time scales, ts,∗ s s,∗
w and tw , by Πw (tw ) = 1 − Ctw = Π
s,∗

ln p(∆E) + ∆E/Tc is plotted in Fig. 3e,f for both mod- and Πw (t∗w ) = Π∗ , where (without loss of generality)
els. The data for GREM inhFig. 3e cani be well fitted Π∗ = 0.2. The ratio R = t∗w /ts,∗
w characterizes the WEB:
Ē)2 if R = 1, then Πw (tw ) follows Πsw (tw ) without dynamical
by an Gaussian function exp − (∆E− 2σ 2 . The variance slowing down; if R → ∞, then Πw (tw ) develops a plateau
σ 2 = aN and mean Ē = kN linearly depend on N , and the weak ergodicity is broken. As shown by the data
as shown in Fig. 3d. From the linear fitting, we obtain in Fig. 4b, R ≈ 1 for T > Td . For T < Td , R > 1 in-
a = 1.12 and k = 0.96 for the GREM. In sharp con- creases with N , suggesting WEB in the thermodynamic
trast, the difference for the EREM in Fig. 3f is negligi- limit. In contrast, in the EREM, R ≈ 1 as long as T > Tc
ble. Thus for the EREM, the barrier distribution is expo- (see Fig. 4c-d). To emphasize the difference between the
nential, and the original BTM works well even for small two models, we show the N -dependence of R at a fixed
systems, which explains why the BTM plateau Hx,w is T = 1.17Tc (see Fig.b-inset). At this T , the EREM is er-
well-observed in the MC data (see Fig. 1b). godic (R ≈ 1), while the GREM shows WEB (ln R ∼ N ,
 5

meaning that R → ∞ when N → ∞). The above nu- the Sherrington-Kirkpatrick model [39], and structural
merical analyses are fully consistent with the prediction glasses.
from the GTM (Fig. 2), i.e., there is an additional dy-
namical transition at Td > Tc in the GREM, but not in
the EREM.
Discussion. This study suggests that the double lim- ACKNOWLEDGMENTS
its N → ∞ and tw → ∞ are not interchangeable. Tak-
ing N → ∞ first and then tw → ∞, as done in previous We warmly thank Hajime Yoshino and Haijun
mathematical computations [23–26], reconciles Eqs. (1) Zhou for discussions. We acknowledge financial sup-
and (4), and therefore the BTM and GTM. However, gen- port from NSFC (Grants 12161141007, 11935002 and
erally in MC simulations, long-time simulations are per- 12047503) and from Wenzhou Institute (Grant WIU-
formed for small N – to describe such settings, one should CASQD2023009). In this work access was granted to
send tw → ∞ and then N → ∞. In the future, it will be the High-Performance Computing Cluster of Institute
interesting to generalize the present analysis to spin-lass of Theoretical Physics - the Chinese Academy of Sciences.
models with a hierarchical energy landscape [1], such as

[1] M. Mézard, G. Parisi, and M. A. Virasoro, Spin glass the- Physical Review E 53, 1525 (1996).
ory and beyond: An Introduction to the Replica Method [14] A. Altieri, G. Biroli, and C. Cammarota, Dynamical
and Its Applications, Vol. 9 (World Scientific Publishing mean-field theory and aging dynamics, Journal of Physics
Company, 1987). A: Mathematical and Theoretical 53, 375006 (2020).
[2] G. Parisi, P. Urbani, and F. Zamponi, Theory of sim- [15] L. F. Cugliandolo and J. Kurchan, On the out-of-
ple glasses: exact solutions in infinite dimensions (Cam- equilibrium relaxation of the sherrington-kirkpatrick
bridge University Press, 2020). model, Journal of Physics A: Mathematical and General
[3] L. C. E. Struik et al., Physical aging in amorphous poly- 27, 5749 (1994).
mers and other materials, Vol. 106 (Citeseer, 1978). [16] J. Kurchan and L. Laloux, Phase space geometry and
[4] R. E. Courtland and E. R. Weeks, Direct visualization of slow dynamics, Journal of Physics A: Mathematical and
ageing in colloidal glasses, Journal of Physics: Condensed General 29, 1929 (1996).
Matter 15, S359 (2002). [17] J.-P. Bouchaud, L. F. Cugliandolo, J. Kurchan, and
[5] P. A. Gago and S. Boettcher, Universal features of an- M. Mézard, Out of equilibrium dynamics in spin-glasses
nealing and aging in compaction of granular piles, Pro- and other glassy systems, Spin glasses and random fields
ceedings of the National Academy of Sciences 117, 33072 12, 161 (1998).
(2020). [18] E. Agoritsas, T. Maimbourg, and F. Zamponi, Out-of-
[6] G. Janzen and L. M. Janssen, Aging in thermal active equilibrium dynamical equations of infinite-dimensional
glasses, Physical Review Research 4, L012038 (2022). particle systems i. the isotropic case, Journal of Physics
[7] M. Baity-Jesi, L. Sagun, M. Geiger, S. Spigler, G. B. A: Mathematical and Theoretical 52, 144002 (2019).
Arous, C. Cammarota, Y. LeCun, M. Wyart, and [19] J.-P. Bouchaud, Weak ergodicity breaking and aging in
G. Biroli, Comparing dynamics: Deep neural networks disordered systems, Journal de Physique I 2, 1705 (1992).
versus glassy systems, in International Conference on [20] J.-P. Bouchaud and D. S. Dean, Aging on parisi’s tree,
Machine Learning (PMLR, 2018) pp. 314–323. Journal de Physique I 5, 265 (1995).
[8] F. H. Stillinger, Energy landscapes, inherent structures, [21] D. J. Gross and M. Mézard, The simplest spin glass, Nu-
and condensed-matter phenomena (Princeton University clear Physics B 240, 431 (1984).
Press, 2015). [22] B. Derrida, Random-energy model: Limit of a family
[9] J.-P. Bouchaud and M. Mézard, Aging in glasses: of disordered models, Physical Review Letters 45, 79
Traps and mode-coupling theory, Progress of Theoreti- (1980).
cal Physics Supplement 126, 181 (1997). [23] J. Černỳ and T. Wassmer, Aging of the metropolis dy-
[10] L. F. Cugliandolo, Dynamics of glassy systems, in Slow namics on the random energy model, Probability Theory
Relaxations and nonequilibrium dynamics in condensed and Related Fields 167, 253 (2017).
matter, edited by J.-L. Barrat, M. Feigelman, J. Kurchan, [24] V. Gayrard, Aging in metropolis dynamics of the rem: a
and J. Dalibard (Springer Berlin Heidelberg, Berlin, Hei- proof, Probability Theory and Related Fields 174, 501
delberg, 2003) pp. 367–521. (2019).
[11] F. Arceri, F. P. Landes, L. Berthier, and G. Biroli, [25] V. Gayrard and L. Hartung, Dynamic phase diagram of
Glasses and aging, a statistical mechanics perspective on, the rem, in Statistical Mechanics of Classical and Disor-
in Statistical and Nonlinear Physics (Springer, 2022) pp. dered Systems: Luminy, France, August 2018 (Springer,
229–296. 2019) pp. 111–170.
[12] L. F. Cugliandolo and J. Kurchan, Analytical solution of [26] B. Derrida, P. Mottishaw, and V. Gayrard, Random en-
the off-equilibrium dynamics of a long-range spin-glass ergy models: broken replica symmetry and activated dy-
model, Physical Review Letters 71, 173 (1993). namics, in Spin Glass Theory and Far Beyond: Replica
[13] L. F. Cugliandolo and P. Le Doussal, Large time nonequi- Symmetry Breaking after 40 Years (World Scientific,
librium dynamics of a particle in a random potential, 2023) pp. 657–677.
 6

[27] M. Baity-Jesi, G. Biroli, and C. Cammarota, Activated

aging dynamics and effective trap model description in
the random energy model, Journal of Statistical Mechan-
ics: Theory and Experiment 2018, 013301 (2018).
[28] D. A. Stariolo and L. F. Cugliandolo, Activated dynamics
of the ising p-spin disordered model with finite number
of variables, Europhysics Letters 127, 16002 (2019).
[29] J.-P. Bouchaud and M. Mézard, Universality classes for
extreme-value statistics, Journal of Physics A: Mathe-
matical and General 30, 7997 (1997).
[30] G. B. Arous and J. Černỳ, Course 8 dynamics of trap
models, Mathematical statistical physics, École d’ÉtÉ de
physique des houches session LXXXIII , 331 (2006).
[31] V. Gayrard, Convergence of clock processes and aging
in metropolis dynamics of a truncated rem, in Annales
Henri Poincaré, Vol. 17 (Springer, 2016) pp. 537–614.
[32] See Supplemental Material at URL, containing aging the-
ory of the generalized trap model and a tree-expansion
theory of the finite-size random energy model.
[33] W. Götze, Complex dynamics of glass-forming liquids: A
mode-coupling theory, Vol. 143 (Oxford University Press,
USA, 2009).
[34] D. R. Reichman and P. Charbonneau, Mode-coupling
theory, Journal of Statistical Mechanics: Theory and Ex-
periment 2005, P05013 (2005).
[35] C. Flamm, W. Fontana, I. L. Hofacker, and P. Schuster,
Rna folding at elementary step resolution, Rna 6, 325
(2000).
[36] F. F. Ferreira, J. F. Fontanari, and P. F. Stadler, Land-
scape statistics of the low-autocorrelation binary string
problem, Journal of Physics A: Mathematical and Gen-
eral 33, 8635 (2000).
[37] J. F. Fontanari and P. F. Stadler, Fractal geometry of
spin-glass models, Journal of Physics A: Mathematical
and General 35, 1509 (2002).
[38] C. Flamm, I. L. Hofacker, P. F. Stadler, and M. T.
Wolfinger, Barrier trees of degenerate landscapes,
(2002).
[39] D. Sherrington and S. Kirkpatrick, Solvable model of a
spin-glass, Physical review letters 35, 1792 (1975).
[40] C. Monthus and J.-P. Bouchaud, Models of traps and
glass phenomenology, Journal of Physics A: Mathemati-
cal and General 29, 3847 (1996).
 7

End Matter -5 a
Appendix A: Random energy models. A REM
comprises 2N configurations. Each configuration consists

E
of N Ising spins, whose energy is drawn randomly from
a Gaussian probability distribution,
1 -15
exp −E 2 /2N .


ρGauss (E) = √ (8) -5
2πN b
This original version with the Gaussian distribution ∆E1
∆E2

Ẽ
ρGauss (E) is called Gaussian random energy model
(GREM). An alternative version, called Exponential ran- τ1
dom energy model (EREM) [27, 29], has been introduced -15 τ2
previously, with ρGauss (E) replaced by
103 tw 104
1 10−2 10−2
ρexp (E) = exp(E/Tc )Θ(−E), (9) c
Tc d
where Θ(x) is the Heaviside step function. In each re- 10−8 10−8
alization, the energy assignment of 2N configurations is ψ(τ ) ψ(τ )

pdf

pdf
fixed (quenched disorder) in the following static barrier- ψ(τ2 |102 ) ψ(τ2 |102 )
tree analyses and dynamical simulations. The procedure 10−14 ψ(τ2 |105 ) 10−14 ψ(τ2 |105 )
is then repeated for ∼ 20000 − 1000 realizations (for ψ(τ2 |108 ) ψ(τ2 |108 )
ψ BT (τ ) ψ BT (τ )
N = 8 − 20) to take the statistical average.
10−20 10−20
102 105 108 1011 102 105 108 1011
10 0 τ τ
N =8
10−1 N = 12 FIG. 6. Basin hopping dynamics in MC simulations.
N = 16 Data are obtained for the GREM at T = 0.75 with N = 20.
−2
10 N = 20
N (a) An example of configuration trajectory, where E(tw ) is
10 15 20 25 the energy of the configuration at tw . (b) Corresponding
10−3 0.4
g(S)

basin trajectory, where Ẽ(tw ) is the energy of the basin or

10−4 0.3 the saddle point at tw . With this, a sequence of hopping
events with the barrier energy ∆Ek and the hopping time τk
⟨S⟩

−5
10 are identified. (c) The hopping time probability distribution
0.2
function (pdf) ψ(τ ) and the conditional distribution ψ(τ2 |τ1 )
10−6
0.1
measured in MC simulations, (c) without and (d) with re-
10−7 turn hops; ψ BT (τ ) is converted from the static distribution
0 0.5 1 1.5 2 2.5 3 pBT
GREM (∆E) using the Arrhenius law.
internal entropy S

FIG. 5. Probability distribution function of the inter- saddle point is then defined as the lowest energy point
nal entropy for the GREM. The data points are g BT (S)
among all maxima (min-max).
obtained by the numerical barrier-tree method, and the solid
line is Eq. (11). (inset) Average numerical entropy ⟨S⟩ as The barrier tree is constructed recursively in the fol-
a function of N , where the dashed line is S∞ ≈ 0.212439 lowing way: (i) find all Nb local minima; (ii) connect the
(Eq. 10). two lowest local minima by a saddle point, and replace
this sub-tree with a new node whose energy is equiva-
Appendix B: Barrier-tree method. The barrier- lent to the saddle point energy (the new set has Nb − 1
tree algorithm searches for all local minima and saddle nodes); (iii) repeat (i) and (ii) until only one node is left
points in the energy landscape, and organizes them into a in the set (Nb = 1). More details about the saddle-tree
barrier tree (see Fig. 3a for an example of a sub-tree). A algorithm can be found in Refs. [35–38].
spin configuration is referred to as a local minimum if its As shown in Fig. 3a, generally a local minimum can
energy is lower than the energy of any adjacent configu- be connected to multiple saddle points. For a given lo-
ration (each configuration has N adjacent configurations cal minimum with an energy Elm , we define Esp as the
related by a single spin flip). The complete set of lo- lowest energy of its connected saddle points, and the bar-
cal minima are found by exhaustive search. To find the rier energy is given by ∆E = Esp − Elm . Our definition
saddle point between two local minima, the algorithm of ∆E is also consistent with the above min-max defini-
first searches for all possible paths (a path is a series of tion of saddle points in Appendix B. In this way, each
subsequent spin flips) between the two minima, with the local minimum is assigned to an energy barrier ∆E. The
maximum energy point identified along each path – the probability distribution of ∆E gives pBT (∆E).
 8

1011 0.4 Appendix C: Obtaining the correlation function

N
−15 −10 −5
3
GREM T = 3/4
Π(tw , tw + t) in MC simulations. Dynamical trajec-
0.3 EREM T = 3/4 tories, which are time sequences of configurations C(t),
2

τ Ψ(τ |∆E)
EREM T = 3/8
are obtained by standard single-spin flip MC simulations,

τ0
0.2 1
106 starting from random initial configurations. An example
τ

0 is provided in Fig. 6a, where the energy E(tw ) of the con-

0.1 10 15 20 figuration at tw is plotted. In order to obtain the aging
0 function Π(tw , tw + t), our task is to transform the con-
20 figuration trajectory to a basin trajectory (Fig. 6b). The
0 5 10 15 102 104 106 108
∆E τ detailed procedure is described below.
During the dynamics, the system is in the basin Bi if
FIG. 7. Arrhenius law. (a) Correspondence between C(t) ∈ Bi and it leaves the basin if C(t) ∈ / Bi anymore.
τ and ∆E in MC simulations of the GREM at T = 0.75 The duration in the basin defines the hopping time (trap-
with N = 20, where the color represents the logarithm of ping time) τ . In this way, we identify a sequence of hop-
ik ik
the joint probability ln Ψ(τ, ∆E) (see the color bar). (b) The ping events with {τk , Elm , Esp , ∆Eik }, where k = 1, 2, · · ·
conditional probability Ψ(τ |∆E = 10) multiplied by τ ; the and 1 ≤ ik ≤ Nb denotes that the k-th event is in the
solid line follows Eq. (13). (b-inset) τ0 vs N . basin ik (Fig. 6b). Note that the transient time between
basins is negligible, i.e., we assume that once the system
leaves a basin, it immediately falls into the next basin.
Each glass basin corresponds to a local minimum and Occasionally, after leaving a basin, the system jumps
a saddle point defined in the above way. A glass basin back to the same basin. Following the previous
also contains other configurations that are neither local study [27], the two events are merged in such cases. The
minima nor saddle points – they form a set of configu- effects of return jumps are revealed in Fig. 6c, d, by the
rations connected to the local minimum whose energies conditional probability p(τ2 |τ1 ), which is distribution of
are all below Esp . Numerically, we start from the given the next-event hopping time τ2 when the current-event
local minimum, and search for its direct neighbours with hopping time is τ1 . Only when the return jumps are
a single-spin flip, the neighbours of neighbours, ..., until merged, the p(τ2 |τ1 ) is independent of τ1 , which means
the configuration’s energy E is larger than Esp . In this that the two consecutive jumps are independent. The in-
way, for a given realization of the REM model with a dependence of consecutive jumps (called renewal mecha-
finite N , we find all i = 1, 2, . . . , Nb basins. Each basin nism) is a basic assumption in the trap model. In order
corresponds to a set Bi of configurations belonging to it, to be consistent with the renewal mechanism, we employ
i i
a local minimum energy Elm , a saddle point energy Esp , the no-return treatment.
i i
and a barrier energy ∆Ei = Esp − Elm . Note that many Appendix D: Verification of the Arrhenius law.
configurations do not belong to any basins. In the trap model, the hopping time τ and the barrier
The internal entropy S = ln Ω is the logarithm of the energy ∆E are related through the Arrhenius law,
number of configurations Ω in a basin. Figure 5 shows
the probability distribution function g BT (S) obtained by τ (∆E) = τ0 exp(∆E/kB T ), (12)
the barrier-tree method, which is nearly independent of
N . The inset of Fig. 5 shows that the mean internal en- where τ0 the microscopic timescale. We find that the
tropy ⟨S⟩ converges to a value S∞ ≈ 0.21 in the large-N Arrhenius law is consistent with the τ̄ (∆E) data (τ̄ is
limit. Then the average number of configurations in a the mean hopping time for the given ∆E) obtained by
basin is Ω ≈ e0.21 ≈ 1.2, i.e., each basin contains ap- MC simulations (Fig. 7a), with τ0 ∼ O(1) as shown in
proximately one configuration in average. Thus in the Fig. 7b-inset (the time unit is N ).
spin glass phase, the mean internal entropy per spin van- When N is finite, for a given ∆E, τ follows a distribu-
ishes (⟨S⟩/N → 0) in the thermodynamic limit, consis- tion
tent with thermodynamic theory of the REM [22]. Al-
ternatively, S∞ can be computed by the tree-expansion Ψ(τ |∆E) = (1 − r)τ −1 r ≈ e−τ r r, (13)
theory (see SI Sec. S2B):
where r = exp(−∆E/T ). To derive Ψ(τ |∆E), consider
∞
X a discretized time τ : if a system is trapped in a basin
S∞ ≡ g∞ (Ω) ln Ω ≈ 0.212439, (10) with a ∆E-barrier for τ steps, it should remain in the
Ω=1
basin for τ − 1 steps with a probability 1 − r at each step,
where and jumps out the basin at the final step with a proba-
Ω! Ω! bility r. Equation (13) is also verified by the simulation
g∞ (Ω) = Ω+1
− (11) data (see Fig. 7b). RNote that the mean hopping time of
Ω (Ω + 1)Ω+1 ∞
Eq. (13), τ̄ (∆E) = 0 τ Ψ(τ |∆E)dτ , consistently recov-
is the probability of basins with Ω configurations. This ers the Arrhenius law Eq. (12). In short, the Arrhenius
theoretical result is fully consistent with the data ob- law is verified by our MC data.
tained by the barrier-tree method (Fig. 5-inset).
 9

Supplementary Information

CONTENTS

Acknowledgments 5

References 5

S1. Aging theory of the generalized trap model 9

A. Outline of the general strategy 9
B. Long-time power-law convergence to the asymptotic plateau 10
C. Short-time linear decay 11
D. Logarithmic decay due to finite-size effects in the generalized trap model 12

S2. Thermodynamics of the finite-size random energy model: a tree-expansion theory 14

A. Probability distributions 14
1. The first order: Ω = 1 15
2. The second order: Ω = 2 15
3. The third order: Ω = 3 16
4. Verification with the numerical barrier-tree algorithm 16
B. Internal entropy 18
C. Number of basins 19

S1. AGING THEORY OF THE GENERALIZED TRAP MODEL

A. Outline of the general strategy

Our theoretical derivation of the behavior of the aging function Π(tw , tw + t) follows the strategy provided in
Refs. [20, 40]. If the system is trapped in a basin with a trapping time τ , then the probability of the system remaining
in the basin after a time t decays exponentially as ∼ e−t/τ . According to this, one can write,
*N + Z ∞
X b
−t/τβ
Π(tw , tw + t) = Q(tw , τβ )e = Nb dτ Q(tw , τ )e−t/τ ψ(τ ), (S1)
β=1 τ0

where Nb is the total number of basins, Q(tw , τ ) is the probability of the system being in a basin with a hopping time
τ at a waiting time tw , ψ(τ ) is the hopping time probability distribution function, and τ0 is the minimum hopping
time (the unit of time). The distribution ψ(τ ) is normalized as,
Z ∞
dτ ψ(τ ) = 1. (S2)
τ0

In the original Bouchaud’s trap model (BTM), the energy barrier distribution is exponential (see Eq. 1). Using the
Arrhenius law τ (∆E) ∼ exp(∆E/T ), the exponential barrier distribution P (∆E) ∼ exp(−∆E/Tc ) can be converted
to a power-law hopping time distribution ψ(τ ),

ψ(τ ) = xτ0x τ −(1+x) , (S3)

where x = T /Tc . Note that Eq. (S3) is normalized according to Eq. (S2).
Bouchaud et al. showed that it is convenient to work in the Laplace transform of Q(tw , τβ ),
Z ∞
Q̂(s, τ ) = L{sQ(tw , τ )} = dtw e−stw sQ(tw , τ ). (S4)
0

When the consecutive hops are independent, as verified numerically in Appendix C, the probability Q̂(s, τ ) is [20, 40],
sτ
sτ +1
Q̂(s, τ ) = D E, (S5)
sτ
CNb sτ +1
 10

where C is a normalization constant. If one considers tw as an exponentially distributed random variable, as suggested
by Eq. (S4), then Eq. (S5) can be understood as Q̂(s, τ ) ∼ τ /(τ + tw ), where the mean waiting time tw ∼ 1/s.
It suggests that the probability of being in a basin with a hopping time τ increases with τ and decreases with the
waiting time tw . Let us consider several limiting cases. For any finite τ , Q̂(s, τ ) → 1 when tw → 0 (the probability
of escaping from any basin is zero without waiting), and Q̂(s, τ ) → 0 when tw → ∞ (all basins can be escaped from
after long-time waiting). For any finite tw , Q̂(s, τ ) → 0 when τ → 0 (the probability to stay in the basin with τ = 0
is zero), and Q̂(s, τ ) → 1 when τ → ∞ (the basin can not be escaped from if its τ is infinite).
With Eqs. (S1), (S3), (S4) and (S5), the aging function Π(tw , tw + t) can be obtained. For the generalized trap
model (GTM), Eq. (S3) needs to be modified in order to include the finite-size correction, but the other equations
can be kept. In the following analyses, the results in Sections. S1 B and S1 C are general for long-time and short-
time dynamics in the BTM and GTM, while the logarithmic decay discussed in Section S1 D is uniquely due to the
finite-size effects in the GTM.

B. Long-time power-law convergence to the asymptotic plateau

Here we give a theoretical derivation of the second line (tw > τ0 ) in Eq. (3). In Eq. (S5), the normalization constant
C = 1 in the large-time limit tw → ∞. With C = 1, the plateau Hx,w of the aging function can be determined as
given by the arcsin law Eq. (2). In order to obtain the behavior of how the aging function converges to the plateau,
it is crucial to consider the next-order 1/tw -correction to C, as shown below. Note that in this section, we use the
BTM distribution Eq. (S3) without finite-size corrections. With Eq. (S3), the average in the denominator of Eq. (S5)
becomes,
  Z ∞
sτ sτ 1
= dτ ψ(τ ) = 2 F 1 (1, x; 1 + x; − ), (S6)
sτ + 1 τ0 sτ + 1 sτ0

where 2 F 1 (a, b; c; z) is the ordinary hypergeometric function. In the large waiting time limit tw → ∞ (or equivalently
in the limit s → 0), 2 F 1 (1, x; 1 + x; − sτ10 ) ≈ xτ0−x s−x sin(πx)
π
. Thus an approximate expression of Q̂(s, τ ) is,

sin(πx)τ0−x sτ
Q̂(s, τ ) ≈ s−x . (S7)
CNb πx (sτ + 1)

The probability distribution Q(tw , τ ) can be obtained by performing an inverse Laplace transform to Eq. (S7) and
then applying the convolution theorem,
sin(πx)τ0−x R tw
n o ′
Q(tw , τ ) = L−1 1s Q̂(s, τ ) ≈ CNb πxΓ(x) 0
dt′ (tw − t′ )x−1 e−t /τ , (S8)

where Γ(x) is the gamma function. The normalization condition requires that,
*N + Z ∞
X b

Q(tw , τβ ) = Nb dτ Q(tw , τ )ψ(τ ) = 1, (S9)

β τ0

or,
sin(πx) R∞ Rt ′
C= × τ0 dτ 0 w dt′ τ −(1+x) (tw − t′ )x−1 e−t /τ . (S10)
πΓ(x)

Note that C depends on tw . Plugging Eqs. (S3) and (S8) into Eq. (S1), we obtain,

sin(πx) ∞
Z Z tw
′
Π(tw , tw + t) ≈ dτ dt′ τ −(1+x) (tw − t′ )x−1 e−(t +t)/τ
CπΓ(x) 0 0
Z tw −x
t + t′

sin(πx) ′
≈ dt (tw − t′ )−1 (S11)
Cπ 0 tw − t′
sin(πx) ∞
Z
= duu−x (1 + u)−1 ,
Cπ w

t+t′
where u = tw −t′ and w = t/tw .
 11

To the zeroth order of tw , Eq. (S10) gives C ≈ 1, and then Eq. (S11) recovers the arcsin law Eq. (2). Expanding
Eq. (S10) to the next order of 1/tw , we obtain,

sin(πx) ∞
 
u tw
Z
−1 −x
C =1− du(1 + u) u Γ x,
πΓ(x) 0 1 + u τ0
 x−1 (S12)
1 tw
≈1− .
(1 − x)Γ(1 − x)Γ2 (x) τ0

Combing Eqs. (S11), (S12) and (2) gives,

"  x−1 #
1 tw
ΠBTM
w (tw ) ≈ Hx,w 1 +
(1 − x)Γ(1 − x)Γ2 (x) τ0
"  −α # (S13)
tw
= Hx,w 1 + A ,
τ0

1
where α = 1 − x and A = (1−x)Γ(1−x)Γ 2 (x) . We have thus derived the large-tw form of Eq. (3). Finally, we compare

the analytic expression Eq. (S13) with the numerical results obtained by the inverse Laplace transform of the exact
expressions Eqs. (S5) and (S6), confirming the validity of Eq. (S13) in the asymptotic regime (see Fig. S1a). Note
that, without loss of generality, we fix w = 0.5 in this study.

100 10−3
x = 0.2 x = 0.81
10−1 x = 0.4 x = 0.27
x = 0.6 x = 0.09
10−2 x = 0.8 x = 0.03
(tw ) − Hx,w

10−4
(tw )

−3
10
1 − ΠBTM

10−4
w
ΠBTM

10−5
10−5
w

10−6

10−7 a 10−6 b
−5 −4 −3
1
10 10 6
10 11 16
10 10 21
10 26
10 10 10 10−2
tw /τ0 tw /τ0

FIG. S1. Comparison between the analytical approximation of the aging function for the BTM and the exact
numerical results using the inverse Laplace transform. (a) Large-tw power-law convergence to the asymptotic plateau.
(b) Small-tw linear decay. The points are exact results obtained by the numerical inverse Laplace transform of Eqs. (S5) and
(S6). The solid lines are Eq. (S13) in (a) and Eq. (S18) in (b).

C. Short-time linear decay

For short-time dynamics, similarly we begin with Eq. (S6) in the last section. In the short-time limit tw → 0
(corresponding to the limit s → ∞), 2 F 1 (1, x; 1 + x; − sτ10 ) ≈ 1, and

1 sτ
Q̂(s, τ ) ≈ . (S14)
CNb sτ + 1

The probability distribution Q(tw , τ ) is again obtained by the inverse Laplace transform:

1 −tw /τ
Q(tw , τ ) = e . (S15)
CNb
 12

Then the aging function is obtained,

Z ∞
Π(tw , tw + t) = Nb dτ Q(tw , τ )e−t/τ ψ(τ )
τ0
 −x   
x t + tw t + tw
= Γ(x) − Γ x,
C τ0 τ0 (S16)
 
 x Γ(x) − Γ x, t+tw
tw τ0
=   ,
t + tw Γ(x) − Γ x, tw
τ0

where the normalized constant is,

Z ∞
C= dτ e−tw /τ ψ(τ )
τ0
 −x    (S17)
tw tw
=x Γ(x) − Γ x, .
τ0 τ0
Expanding Eq. (S16) around tw = 0, we obtain,
   
wx tw tw
Πsw (tw ) ≈1− =1−C . (S18)
x+1 τ0 τ0
Similar to the long-time case, the analytic result Eq. (S18) can be confirmed by exact numerical calculations using
the inverse Laplace transform (see Fig. S1b).

D. Logarithmic decay due to finite-size effects in the generalized trap model

The GTM aging function is computed using the formula Eq. (S1),
Z ∞
ΠGTM
w (tw ) = dτ ψGTM (τ )QGTM (tw , τ )e−wtw /τ
τ
Z 0∞ (S19)
= dτ efGTM (tw ,ln τ ) ,
0

The GTM barrier energy distribution Eq. (4) leads to the modified trapping time distribution (via the Arrhenius law),

ψGTM (τ ) ∼ τ −xeff−1 exp −ξ ln2 (τ /τ0 )

 
(S20)
∼ τ −xeff −ξ ln(τ /τ0 )−1
 
T2 ĒTc
where ξ ≡ 2N a , and x eff = 1 − N a x. The QGTM (tw , τ ) can be computed by the inverse Laplace transform of
D E R
∞
Eq. (S5), where sτsτ+1 = τ0 dτ sτsτ+1 ψGTM (τ ).
It is hard to directly analyze Eq. (S19). However, based on a saddle-point approximation, we find that the aging
functions in the two models can be related by shifting the effective reduced temperature,
ΠGTM
w (tw ; x) ∼ ΠBTM
w (tw ; x̃eff ), (S21)
where
T2
 
ĒTc
x̃eff = 1− x+ ln tw = xeff + 2ξ ln tw . (S22)
Na Na

In other words, the aging function ΠGTM w (tw ; x) of the GTM at a temperature T = Tc x is approximately equivalent
to the aging function ΠBTM
w (t w ; x̃ eff ) of the BTM at the temperature T̃eff = Tc x̃eff with x̃eff given by Eq. (S22).
Equation (S21) is essentially due to the relation between the maximum of fGTM (tw , ln τ ) in Eq. (S19) and that of
fBTM (tw , ln τ ),
∗ ∗
ln τGTM (tw , x) ≈ ln τBTM (tw , x̃eff ), (S23)
 13

a b
15 15

(tw , xeff )

ln τGTM (tw , x)
ln τBTM
10 T = 0.21 10 T = 0.21

∗
∗

T = 0.42 T = 0.42
T = 0.64 T = 0.64
T = 0.85 T = 0.85

6 8 10 12 14 16 6 8 10 12 14 16
ln tw ln tw
16
c d
0.2

ln τ ∗ (tw )/N
14
ln τ ∗ (tw )

0.15
T = 0.21 N = 64
12 T = 0.42 N = 128
T = 0.64 N = 256
T = 0.85 0.1 N = 512
10
8 10 12 14 0.1 0.15 0.2
ln tw ln tw /N

FIG. S2. Verification of Eq. (S23). The data points are obtained by numerical evaluation of the maximum f (tw , ln τ ∗ ) of
the function f (tw , ln τ ) defined in Eq. (S19). As shown by the data, ln τ ∗ is linearly proportional to ln tw in (a) BTM and (b)
GTM, for T = 0.21, 0.42, 0.64, 0.85 (or x = 0.25, 0.5, 0.75, 1.0), and N = 64. At the same T , the BTM and GTM data do not
collapse in (a,b). However, the GTM data (dotted point-line) at T = Tc x and the BTM data (solid point-line) at T̃eff = Tc x̃eff
collapse, where x̃eff is given by Eq. (S22), for different (c) T (with N = 128 fixed) and (d) N (with T = 0.42 fixed).

where f (tw , ln τ ) is maximized at ln τ ∗ with other parameters (tw , T , etc.) fixed. Equation (S23) is consistent with
the numerically evaluated maximum point ln τ ∗ (see Fig. S2). As shown by the data, in general ln τ ∗ ∼ ln tw in both
models (Fig. S2a,b). The GTM data at x and the BTM data at the corresponding x̃eff (see Eq. S22) collapse for
different x and N (Fig. S2c,d), and thus Eq. (S23) is verified.
Equation (S22) means that, the Gaussian term in the GTM barrier distribution Eq. (4) gives rise to two effects
compared to the BTM. First, it shifts the reduced temperature effectively from x = T /Tc to an N -independent xeff =
Teff /Tc . This shift modifies the asymptotic plateau of the aging function from Hx,w to Hxeff ,w , in the thermodynamic
limit N → ∞. Second, the next-order correction adds a term N1 ln tw as x̃eff − xeff ∼ N1 ln tw . This correction
disappears in the thermodynamic limit, but brings in a ln(tw ) decay term to the aging function for a finite N , as we
show below.
With Eq. (S21), ΠGTM w (tw ) can be conveniently analyzed using the Arcsin law, with x in Eq. (S13) replaced by
x̃eff . For large N , the difference between x̃eff and xeff , x̃eff − xeff ∼ N1 ln tw , can be considered as a perturbation to
the original form. Expanding around xeff to the first order, we obtain,
"  −α #
tw
ΠGTM
w (tw ) ≈ Hx̃eff ,w 1 + A
τ0
"  −α #
T 2 Hx′ eff ,w
 
tw tw
≈ Hxeff ,w 1 + A + ln (S24)
τ0 Na τ0
"  −α #  
tw B tw
= Hxeff ,w 1 + A − ln ,
τ0 N τ0

T 2 Hx′ ,w ′
where B = − a
eff
, and Hx,w is the derivative of Hx,w with respect to x. The analytic expression of the logarithmic
 14

1 1
N = 64 N = 64
N = 128 N = 128

(tw )

(tw )
N = 256 N = 256
0.9 N = 512 0.8 N = 512

ΠGTM

ΠGTM
w

w
0.8
a 0.6 b
10−3 101 105 109 1013 10−3 101 105 109 1013
tw /τ0 tw /τ0

FIG. S3. Logarithmic decay in the GTM at(a)T = 0.64 and (b) T = 1. The data points are obtained by numerical
B
integration of Eq. (S19). The lines are the − N ln tτw0 term in Eq. (S24).

1
N =8
0.8 N = 20
Πw (tw )

0.6

0.4

0.2

100 102 104 6

108 1010
tw 10
FIG. S4. Comparison between the MC aging function (points) and Eq. (S24) (lines), at T = 0.75. No fitting
parameters are used for the lines.

T 2 Hx′ ,w
 
B tw
decay − N ln τ0 with B = − is compared to the numerical integration of Eq. (S19) in Fig. S3: the agreement
a
eff

 
is converged with increasing N . We further compare the theoretical results − NB
ln tτw0 with the MC data in Fig. S4.
While the logarithmic form is robust, the coefficient B is not exact due to strong higher-order corrections for small
N in MC simulations. In the main text Fig. 1, we treat B as a fitting parameter.

S2. THERMODYNAMICS OF THE FINITE-SIZE RANDOM ENERGY MODEL: A TREE-EXPANSION

THEORY

A. Probability distributions

Our aim is to obtain statistical properties of the energy landscape for the finite-sized REM model with N spins,
where N is sufficiently small so that the system is far from the thermodynamic limit (N → ∞). We assume that each
basin only has one nondegenerate local minimum, and that basins are independent. The probability distribution of
the local minimum Elm is,
Z ∞ N
plm
N (Elm ) = (N + 1)ρ(Elm ) dEρ(E)
Elm (S25)
= (N + 1)ρ(Elm )LN (Elm ),
 15
R∞
where L(E) ≡ E dEρ(E) is a complementary cumulative distribution function, and ρ(E) is the probability distribu-
tion function of the configuration energy E. The above expression requires that the N direct neighbors of the local
minimum, which are related to the local minimum by a single spin flip, have an energy E higher than Elm . The
factor N + 1 comes from the N + 1 choices of the local minimum among the N + 1 configurations. For the GREM,
ρ(E) = ρGauss (E) given by Eq. (8), which results in L(E) = LGauss (E), where
1 √
LGauss (E) = erfc(E/ 2N ), (S26)
2
with erfc(x) the complementary error function. For the EREM, ρ(E) = ρexp (E) as in Eq. (9), and L(E) = Lexp (E),
where
Lexp (E) = [1 − exp(E)]Θ(−E). (S27)
The probability distributions of the saddle point energy Esp and the barrier energy ∆E cannot be explicitly obtained.
They need to be computed by integrating out the other variable in the joint distribution λN (Elm , Esp ),
Z +∞
sp
pN (Esp ) = dElm λN (Elm , Esp ), (S28)
−∞

Z +∞
pN (∆E) = dEsp λN (Esp − ∆E, Esp ). (S29)
−∞

The joint distribution λN (Elm , Esp ) characterizes the probability of a basin with a local minimum energy Elm and
a saddle point energy Esp . Our key task is to compute λN (Elm , Esp ). To do that, we develop an approach named
tree-expansion theory. The starting point is to write λN (Elm , Esp ) as a summation,
N
(Ω)
X
λN (Elm , Esp ) ≡ λN (Elm , Esp ), (S30)
Ω=1

where Ω is the number of configurations in the basin. For example, Ω = 1 means that there is only a local minimum
in the basin, and Ω = 2 means that there is another configuration in the basin besides the local minimum, etc. The
summation converges quickly with the increasing Ω. We next explicitly consider the first three orders, Ω = 1, 2, 3.

1. The first order: Ω = 1

When Ω = 1, as shown in Fig. S5, the basin contains only the local minimum without any other configurations.
By definition, the saddle point should have the lowest energy among the N neighbors of the local minimum – the
probability of this condition is LN −1 (Esp ). The saddle point also has N neighbors, including the local minimum and
N − 1 other neighbors. The energies of the other N − 1 neighbors cannot be all higher than Esp – otherwise the saddle
point would be a configuration inside the basin (not a saddle point). This condition imposes a constraint given by
1 − LN −1 (Esp ). Putting these considerations together, we can write the first-order joint distribution as
(1)
λN (Elm , Esp ) = (N + 1)N ρ(Elm )ρ(Esp )LN −1 (Esp ) 1 − LN −1 (Esp ) ,
 
(S31)
where the extra factor (N +1)N comes from the permutation of choosing a local minimum from the N +1 configurations
and a saddle point from the rest of the N configurations.

2. The second order: Ω = 2

The second-order joint distribution corresponds to two graphs as shown in Fig. S6. Besides the local minimum, the
R Esp
basin contains one configuration whose energy is between Esp and Elm – the corresponding probability is Elm dEρ(E).
The permutation of choosing one saddle point, one local minimum and another configuration in the basin gives a factor
of (N + 1)N (N − 1). The two graphs in Fig. S6 coincidently have the same expression, which gives an extra factor of
two. Finally, the second-order joint distribution can be written as,
 Esp
Z
(2)
λN (Elm , Esp ) = 2(N + 1)N (N − 1)ρ(Elm )ρ(Esp )L2N −3 (Esp ) 1 − LN −2 (Esp )

dEρ(E)
Elm (S32)
2N −3 N −2
 
= 2(N + 1)N (N − 1)ρ(Elm )ρ(Esp )L (Esp ) 1 − L (Esp ) [L(Elm ) − L(Esp )]
 16

N-1

...
N-1
Saddle

...

Local minimum

FIG. S5. Graphic representation of the first-order joint distribution Eq. (S31). Each node represents a configuration,
and the height of the node represents the configuration energy E. Any pair of configurations connected by a link are related
by a single-spin flip. The saddle point energy Esp is marked by the dotted line. For Ω = 1, there is only one configuration in
the basin (E < Esp ), which is the local minimum. The number of nodes in the circled cluster is indicated (in this case, both
clusters have N − 1 nodes). For simplification, only N = 4 nodes are shown.

N-1 N-1
N-2 N-2
... ... ... ...

N-2 N-2
Saddle Saddle

... ...

Local minimum Local minimum

FIG. S6. Graphic representation of the second-order joint distribution Eq. (S32).

3. The third order: Ω = 3

To compute the third-order joint distribution, we need to consider six graphs as shown in Fig. S7. The top three
graphs correspond to the term with 1 − LN −3 (Esp ) in the following expression Eq. (S33), and the bottom three graphs
correspond to the term with 1 − LN −2 (Esp ):

(3)
λN (Elm , Esp ) = 3(N + 1)N (N − 1)ρ(Elm )ρ(Esp )L3N −6 (Esp )×
Z
N −3 N −2
    
× (N − 2) 1 − L (Esp ) + (N − 1) 1 − L (Esp ) dE1 dE2 ρ(E1 )ρ(E2 )
Esp >E1 >E2 >Elm
(S33)
3
= (N + 1)N (N − 1)ρ(Elm )ρ(Esp )L3N −6 (Esp )×
2
2
× (N − 2) 1 − LN −3 (Esp ) + (N − 1) 1 − LN −2 (Esp ) [L(Elm ) − L(Esp )] .
    

4. Verification with the numerical barrier-tree algorithm

The above theoretical results are compared with the numerical data obtained by the barrier-tree algorithm for both
GREM and EREM with N = 24 spins (Fig. S8). The expression Eq. (S25) of the local minimum energy distribution
sp
plm
N (Elm ) is tested in Fig. S8a. The tree-expansion results of the saddle point energy distribution pN (Esp ) (see Eq. S28)
and the barrier energy distribution pN (∆E) (see Eq. S29) are tested in Fig. S8b-d. It can be seen that the first-order
results are already very close to the numerical data. Note that the distribution ρ(E) in the theoretical expressions
shall be replaced by Eq. (8) and (9) in the End Matter for the corresponding REMs.
 17

N-2 N-2 N-2 N-2 N-2 N-2 N-2 N-2 N-2 N-2
... ... ...
... ... ...
N-2 ... N-2 ... ...
Saddle Saddle Saddle

... ... ...

B B
A A

Local minimum Local minimum Local minimum

N-3 N-2 N-1 N-3 N-2 N-1 N-2 N-1

... ... ... N-3
... ...
N-3 N-3 ... N-3 ...
Saddle Saddle Saddle
... ... ... ...
...

Local minimum Local minimum Local minimum

FIG. S7. Graphic representation of the third-order joint distribution Eq. (S33).

a b
GREM 0.6 GREM
0.3 EREM EREM
0.2 0.4
pdf

pdf

0.1 0.2

0.0 0.0
15 10 5 0 10 5 0
Elm Esp

10 1 GREM 10 1 EREM
3 e E/Tc 3
10 10 E/Tc
e
pdf

pdf

5 5
10 10
7 7
10 10
c d
10 9 0 10 20 10 9 0 5 10 15 20
E E

FIG. S8. Comparison between the tree-expansion theory (lines) and the barrier-tree method (points), for the
GREM and EREM with N = 24 spins. (a) Probability distribution function (pdf) of the local minimum energy. The line
represents Eq. (S25). (b) Pdf of the saddle point energy. (c,d) Pdf of the barrier energy distribution. In (b-d), the red line
represents the tree-expansion theoretical results up to the first-order (Ω = 1), and the black line represents the results up to
the third-order results (Ω = 3).
 18

0.8

0.6
Ω=1 GREM
Ω=2 EREM

gN (Ω)
0.4 Ω=3

0.2

0
5 10 15 20 25 30
N

FIG. S9. Probability gN (Ω) of having Ω configurations in a basin. Lines are tree-expansion theory results Eq. (S35), and points
are numerical data obtained by the barrier-tree algorithm.

B. Internal entropy

The probability of having Ω configurations in a basin is,

Z ∞ Z Esp
(Ω)
gN (Ω) = dEsp dElm λN (Elm , Esp ). (S34)
−∞ −∞

Using the above tree-expansion results, Eqs. (S31), (S32) and (S33), we obtain the first three orders,

Z ∞
dEsp ρ(Esp )N (N + 1)LN −1 (Esp ) 1 − LN −1 (Esp ) [1 − L(Esp )]
 
gN (1) =
−∞
Z1
= dσN (N + 1)σ N −1 (1 − σ N −1 )(1 − σ)
0
3N − 3
=
4N − 2
N →∞ 3
= ,
4
Z 1
gN (2) = dσ(N − 1)N (N + 1)σ 2N −3 (1 − σ N −2 )(1 − σ)2
0 (S35)
(N + 1)(N − 2)(19N − 12)
=
6(2N − 1)(3N − 4)(3N − 2)
N →∞ 19
= ,
108
Z 1
3N −6 N − 2 N −1
 
N −3 N −2
gN (3) = dσ(N − 1)N (N + 1)σ (1 − σ )+ (1 − σ ) (1 − σ)3 ,
0 2 2
N (N + 1)(2N − 3) 3(N − 1)N (N + 1)
= −
(3N − 5)(3N − 4)(3N − 2) 4(2N − 3)(4N − 7)(4N − 5)
N →∞ 175
= .
3456

Note that the results in Eq. (S35) are model-independent, i.e., independent of ρ(E). The numerical data obtained by
the barrier-tree algorithm confirm the tree-expansion theory Eq. (S35) (see Fig. S9).
In the thermodynamic limit (N → ∞), the probability gN (Ω) converges to an analytical expression that is universal
 19

for any Ω,
1
N!
Z
1
g∞ (Ω) = lim dσ σ ΩN − 2 Ω(Ω+1) (1 − σ N −Ω )(1 − σ)Ω
N →∞ 0 (N − Ω − 1)!
(ΩN − 12 Ω(Ω + 1))! ((Ω + 1)N − 12 Ω(Ω + 3))!
 
N !Ω!
= lim − (S36)
N →∞ (N − Ω − 1)! (ΩN − 1 Ω(Ω + 1) + 1)! ((Ω + 1)N − 12 Ω(Ω + 1) + 1)!
2
Ω! Ω!
= Ω+1 − .
Ω (Ω + 1)Ω+1

The internal entropy of the whole system in the thermodynamic limit is

∞
X
S∞ ≡ g∞ (Ω) log Ω ≈ 0.212439. (S37)
Ω=1

This result is consistent with the numerical data obtained by the barrier-tree method (see Fig. 5 in End Matter).
Note that the internal entropy per spin S∞ /N vanishes in the thermodynamic limit.

C. Number of basins

We assume that each basin contains one local minimum. Two local minima cannot be directly related via a single-
spin flip – otherwise the two basins cannot be separated. Thus, in average, we can assume that each basin occupies
N + 1 configurations, i.e., one local minimum and its N neighbours. The total number of configurations for a N -spin
REM is 2N . Therefore, the total number of basins is,

2N
Nb = . (S38)
N +1
Figure S10 shows that this simple consideration Eq. (S38) describes well the numerical data obtained by the barrier-
tree method.

GREM
10 6 EREM

104
Nb

102

100
10 20 30
N
FIG. S10. Number of basins Nb as a function of the system size N . The line represents Eq. (S38), and the points are data
obtained by the numerical barrier-tree algorithm.