From entropic to energetic barriers in glassy dynamics: The Barrat-Mézard trap model on sparse networks
FFrom entropic to energetic barriers in glassy dynamics: The Barrat–M´ezard trapmodel on sparse networks
Diego Tapias, ∗ Eva Paprotzki, and Peter Sollich
1, 2, † Institut f¨ur Theoretische Physik, University of G¨ottingen,Friedrich-Hund-Platz 1, 37077 G¨ottingen, Germany Department of Mathematics, King’s College London, London WC2R 2LS, UK (Dated: May 12, 2020)Trap models describe glassy dynamics as a stochastic process on a network of configurationsrepresenting local energy minima. We study within this class the paradigmatic Barrat–M´ezardmodel, which has Glauber transition rates. Our focus is on the effects of the network connectivity ,where we go beyond the usual mean field (fully connected) approximation and consider sparsenetworks, specifically random regular graphs. We obtain the spectral density of relaxation rates ofthe master operator using the cavity method, revealing very rich behaviour as a function of networkconnectivity c and temperature T . We trace this back to a crossover from initially entropic barriers,resulting from a paucity of downhill directions, to energy barriers that govern the escape from localminima at long times. The insights gained are used to rationalize the relaxation of the energy aftera quench from high T , as well as the corresponding correlation and persistence functions. I. INTRODUCTION
Trap models have been widely used to model glassydynamics [1–15]. They abstract the full dynamics of su-percooled liquids, structural glasses and amorphous sys-tems more generally, into motion near local potential en-ergy minima (inherent structures), interspersed with rel-atively rare transitions between minima. Each minimumwith its surrounding potential energy basin is identifiedas a discrete “trap”, and the possible transitions definea stochastic dynamics on a network of these traps. Inmore sophisticated versions of the approach, traps aredefined by clustering nearby basins with rapid transitionsbetween them into so-called meta-basins [8, 9, 16].The most widely studied trap models are those inwhich the set of traps is fully connected, i.e. the meanfield trap models. The full network connectivity allowsthese models to be studied analytically using differentmethods [1–5, 10, 17, 18], yielding predictions both forthe relaxation of single-time quantities like the aver-age energy and for two–time correlations and responses.Much less is known for the more realistic case where thedegree of each node is small compared to the number oftraps. Our work in this paper is designed to fill this gap,and indeed we will find that sparse network connectivityhas very rich physical effects that have no analogue inmean-field trap models.In addition to the network connectivity as a crucial in-gredient in our analysis, a trap model is defined by the energies of the traps, and the transition rates betweentraps. These should satisfy detailed balance with respectto the trap energies, to ensure that the eventual equilib-rium that a system of finite size would reach is a Boltz-mann distribution. These aspects all remain fixed during ∗ [email protected] † [email protected] the time evolution, so for theoretical analysis it is conve-nient to think of the network structure and trap energies– which determine the transition rates – as quenched dis-order drawn from appropriate distributions [19–21].Arguably the simplest trap model is the Bouchaudmodel, which has Arrhenius transition rates that pre-sume activation to some common threshold energy leveland depend only on the departing trap. For this modelrecent studies [19, 22] have considered in detail the ef-fect of the structure of the network of traps, with a focuson the spectral and time domain properties of the mas-ter operator that defines the dynamics. In this paper weuse a similar set of tools to analyze the Barrat–M´ezard(BM) trap model on Random Regular Graphs (RRG).The BM model uses Glauber transition rates, which areessentially constant for transitions that lower the energybut of Arrhenius form for energy increases. This is an at-tractive feature that allows the model to encode both con-ventional energy barriers that need to be surmounted bythermal activation and entropic barriers that arise whenthe number of possible downhill transitions from typicaltraps becomes small [7, 23]. Thus the model includes thetwo main mechanisms that are responsible for the slowdynamics in glasses.The spectral properties we concentrate on are the re-laxation rates of the stochastic system, i.e. the negativeeigenvalues of the master operator. This operator is alarge random matrix due to the disorder in network struc-ture and trap depths and so we import tools from randommatrix theory for our analysis, specifically the CavityMethod [24–26]. From the distribution of the relaxationrates one can predict qualitatively how the system willrelax in the time domain as illustrated by recent studiesfor a range of different physical systems [27, 28].In Sec. II we define the BM model, the RRG networkstructures we consider, and the associated master oper-ator. The cavity method we use to work out its spec-tral properties is set out in Sec. III. We present results(Sec. IV) first for the mean-field limit of large connec- a r X i v : . [ c ond - m a t . d i s - nn ] M a y tivity, where we show that the relaxation rate spectrumis given by the distribution of the escape rates from theindividual traps. This is followed by the analysis for thegeneral case of sparse connectivity and arbitrary tem-perature, with a focus on temperatures below the glasstransition. In Sec. V we translate the results for the re-laxation rate spectrum into the time domain to analyzethe behavior of a one-time observable – specifically themean energy – and two-time correlation and persistencefunctions. We conclude in Sec. VI with a discussion ofthe results and perspectives for future work. II. BARRAT–M´EZARD TRAP MODELA. Barrat–M´ezard dynamics
The Barrat–M´ezard model considers dynamics on a setof traps that we label by i = 1 , . . . , N . For traps i and j that are connected, i.e. are neighbours in the network oftraps, the transition rate is of Glauber form, W ji := W i → j = 1 c
11 + exp( − β ( E j − E i )) , (1)where the normalization by c , the average number ofneighbours, ensures that the total escape rate from atypical trap is of order unity. As usual, β = 1 /T isthe inverse temperature ( k B = 1). E i and E j are thedepths, i.e. the negative energies, of the traps; we setthe origin of the energy scale so that E i > i . Two limiting cases can be obtained from (1): for E j − E i (cid:29) T , W ji → /c , which means that transi-tions to significantly deeper traps take place at a ratethat is independent of temperature.This is what allowsentropic effects to occur in the dynamics, especially atlow T where the system goes downhill in the energy land-scape [3, 7, 23]. In the opposite limit E i − E j (cid:29) T onehas W ji ∝ exp( − β ( E i − E j )); this Arrhenius form is thesignature of the thermally activated nature of transitionsto higher energy traps [13]. B. Energy distribution
The distribution of trap depths is usually taken to beeither exponential or Gaussian. The latter choice leads toa REM-like model [14] in which the energy may be posi-tive or negative. The exponential distribution, which wasused in the original mean field Bouchaud trap model [1, 4]and much of the literature since, is more physical if wethink of the traps as representing local energy minimaand it is justified from extreme–value statistics [29]. Thethreshold at depth E = 0 can then be understood asthe top of this energy landscape of local minima. In thefollowing we choose the width of the trap depth distri-bution as our energy scale, which in dimensionless unitsis then just ρ E ( E ) = exp( − E ). The Glauber transition rates of the BM model obey detailed balance with re-spect to the Boltzmann distribution , as can be checkedexplicitly from (1), so if the system equilibrates thenthe probability of finding it in a trap of depth E is ∝ ρ E ( E ) e βE = e ( β − E . This distribution becomes un-normalizable at the temperature T g = 1, which definesthe glass transition. For T < T g the system ages towardsdeeper and deeper traps and it is this regime of glassydynamics that we will mostly focus on. C. Network structure
To represent the potentially sparse connectivity of thenetwork of traps we model this network as one of thesimplest choices, a random regular graph (RRG) withconnectivity c . This means explicitly that the network israndomly selected among all networks where every node(trap) has exactly c neighbours. Consistent with the no-tation above, c is then also the average number of neigh-bours. To ensure sparse connectivity we will keep c finitewhile taking the thermodynamic limit N → ∞ , so thatalways c (cid:28) N .The RRG contains the essential features of sparse ran-dom networks by confining all dynamical transitions tothe local environment of a node. As appropriate for aconfiguration space model it is also infinite-dimensionalin the sense that the number of nodes grows exponen-tially with distance. In this work we focus on connectivity c ≥
3. This condition ensures that the fraction of nodesoutside the giant connected component vanishes for large N [19], which implies that the network is connected asexpected on physical grounds.For our analysis it will be important that, because oftheir sparse connectivity, RRGs become locally tree-likefor N → ∞ . This is a consequence of the fact that typicalloop lengths are ∼ ln( N ) and so diverge in the limit [30].The locally tree-like structure is what enables us to usethe cavity method [24, 25] to obtain the spectral proper-ties, in a way that becomes exact for N → ∞ .We mention in passing that an RRG is effectively aBethe lattice with no boundaries [31] and thus differentfrom a Cayley tree. This technical detail is relevant inthe context of localization of eigenvectors on Bethe lat-tices [31–34]. D. Master equation
We can now write down the master equation governingthe time evolution of the vector p ( t ) = ( p ( t ) , . . . , p N ( t ))of probabilities of being in any of the N traps. If A ij is the adjacency matrix of the network of traps, with A ij = 1 if two traps are connected and = 0 otherwise,then the master equation readsd p ( t )d t = M p ( t ) . (2)The master operator M has elements M ji = A ji W i → j , M ii = − (cid:88) j (cid:54) = i M ji (3)As can be seen from their definition, the (negative) diago-nal are the rates of escape − M ii = Γ i from the individualtraps. We recall that the transition rates W i → j are givenby (1) and the trap depths are sampled independentlyfrom the exponential distribution ρ E defined in II B. III. SPECTRAL ANALYSIS
The general solution of the master equation (2) can bewritten as [19] p ( t ) = N − (cid:88) α =0 e λ α t ( L α , p (0)) R α , (4)where the L α and R α are the left and right eigenvectorsof M , respectively, and the λ α are the correspondingeigenvalues. These are non-positive and we order themin the following so that 0 = λ ≤ − λ ≤ − λ ≤ . . . ≤− λ N − ; note that the connectedness of the network en-sures that there is only a single zero eigenvalue, whichdescribes equilibrium [35].Eq. (4) shows that the dynamics is a superposition ofexponentials with relaxation rates − λ α . The spectrumof these relaxation rates therefore play a key role in thephysical behaviour. It is described in the thermodynamiclimit by the eigenvalue distribution or spectral density ρ ( λ ) = lim N →∞ N (cid:42) N − (cid:88) α =0 δ ( λ − λ α ) (cid:43) , (5)The average here is over the graph realizations (struc-tural disorder) and the local trap depths (energetic dis-order), though for large N this is essentially equivalent toconsidering a single realization because of self-averaging. A. Methods
In order to develop a theory for the spectrum we firstsymmetrize the master operator. Because of detailed bal-ance, this can be achieved using a similarity transforma-tion with the equilibrium distribution, M s = P − / M P / (6)Here P eq is a diagonal matrix with entries ( P eq ) ii = p eq i and the p eq i ∝ e βE i are the equilibrium occupations ofthe traps. Such a similarity transformation preserves thespectrum, but the fact that M s is symmetric is impor-tant for setting up the cavity theory [19]. In the same spirit we write the transition rates (1) aseffective Bouchaud rates, which are of Arrhenius form,times a symmetric function: W ji = e β ( E i + E j ) / β ( E i − E j ) /
2) e − βE i c =: K ( E i , E j ) e − βE i c . (7)The spectral density can in general be deduced fromthe properties of the resolvent G matrix, which is definedas [24, 36, 37] G ( λ − i(cid:15) ) = (( λ − i(cid:15) ) I − M s ) − , (8)by the relation ρ M ( λ ) = lim (cid:15) → πN N (cid:88) j =1 Im G jj , (9)where in the equation (8), I denotes the identity matrix.As indicated by the superscript, the spectrum here isinitially for a fixed master operator and finite N , while (cid:15) is a regularizer that broadens the delta-distributions inthe definition of the spectrum into Lorentzians of width (cid:15) . The cavity method starts from the fact that the resol-vent entries G jj can be thought of as the variances ofan N -dimensional Gaussian distribution with exponent − x T G − x /
2. This exponent has nonzero cross terms x i x j only for neighbouring traps i and j , so inherits thestructure of the network of traps. If the network is a tree,the variance of each x j can then be worked out by what isknown in the machine learning literature as Belief Prop-agation [25]: the “marginal” variance 1 /ω j at any nodecan be expressed in terms of the “cavity variances” 1 /ω ( j ) k of the neighbours. “Cavity” here means that 1 /ω ( j ) k is thevariance at node k if node j had been removed from thenetwork, i.e. a cavity created; ω ( j ) k can therefore also becalled a cavity precision or cavity Green’s functions [31].Intuitively, ω ( j ) k is a “message” that node k sends to node j and that encodes all properties of the part of the net-work that feeds into j via k .In practice the equations take a simpler form if onerescales the variables x in the Gaussian probability dis-tribution so that the resolvent elements are expressed interms of the marginal variances as G jj = i e βE j cω j (10)The relation between ω j and the cavity precisions is then ω j = i ( λ − i(cid:15) )e βE j c + (cid:88) k ∈ ∂j iK ( E j , E k ) ω ( j ) k iK ( E j , E k ) + ω ( j ) k , (11)where ∂j indicates the set of neighbours of j . The cav-ity precisions themselves can be obtained from a set ofequations that is almost identical except that node j isremoved from the sum: ω ( j ) k = i ( λ − i(cid:15) )e βE k c + (cid:88) l ∈ ∂k \ j iK ( E k , E l ) ω ( k ) l iK ( E k , E l ) + ω ( k ) l . (12)Once these cavity equations have been solved, the spec-trum for a given master operator M can be obtainedfrom equation (9) with (11) inserted, i.e. ρ M ( λ ) = lim (cid:15) → πN N (cid:88) j =1 Re (cid:0) e βE j c/ω j (cid:1) . (13)In line with the fact that the function K is the factorby which BM transition rates differ from those in theBouchaud model (see (7) ), Eqs. (11, 12) reduce to thosefor the Bouchaud case [19] if we replace K by 1.The cavity equations (12) are in principle approximateon any finite random regular graph because of the pres-ence of loops, but in the thermodynamic limit N → ∞ they are expected to become exact as loop lengths di-verge. In that limit one also sees that, for given E k ,the sum on the r.h.s. of (12) consists of c − distribution ζ ( ω, E ) ofcavity precisions. We emphasize that as each ω ( j ) k ≡ ω depends on the energy E k ≡ E of the sending trap , sucha joint distribution is required, in contrast to the sim-pler Bouchaud case [19]. The cavity equations (12) thenbecome a self-consistent equation for ζ : ζ ( ω, E ) = ρ E ( E ) (cid:90) c − (cid:89) j =1 d E j d ω j δ ( ω − Ω c − ) ζ ( { ω j , E j } ) , (14)with the abbreviationΩ a ( { ω l , E l } , E ) = iλ (cid:15) e βE c + a (cid:88) l =1 iK ( E, E l ) ω l iK ( E, E l ) + ω l (15)and λ (cid:15) ≡ λ − i(cid:15) . Eq. (14) can then be solved using a Pop-ulation Dynamics algorithm [31, 38]. In this method asadapted to our case one starts from a finite populationof ( ω, E )-pairs, initialized from some in principle arbi-trary initial distribution. The population is then updatediteratively by picking c − ω l , E l )and a new energy E from the exponential distribution ρ E ( · ), calculating Ω c − ( { ω l , E l } , E ) and replacing a ran-dom population member by (Ω c − , E ) [39]. After thepopulation has equilibrated, we then calculate the spec-tral density (9) from the population analogue of (13), as ρ ( λ ) = lim (cid:15) → π Re (cid:28) e βE c Ω c ( { ω l , E l } , E ) (cid:29) ( { ω l ,E l } ,E ) . (16)For all the results presented below, a population of 10 pairs { ( ω, E ) } was used unless otherwise specified. For every given λ , 10 iterations were performed to reachequilibrium, with (cid:15) set to (cid:15) = 10 − . The populationwas then evolved for a further 10 steps and at each up-date step a sample (Ω c , E ) was drawn for the calculationof the average (16). We evaluate this average simulta-neously for different values of (cid:15) , usually { − , − } , toassess the (cid:15) -dependence of the results. For a more de-tailed discussion of the role of (cid:15) we refer to [19]. IV. RESULTS: SPECTRAL PROPERTIESA. Escape rates versus relaxation rates
We first consider the mean-field case c → ∞ , where therelaxation rate spectrum has previously been worked outfor the Bouchaud [20] but not the BM trap model (withthe exception of [5] where the case T = 0 is analyzed).We next show that in this limit the relaxation rates be-come the (negative) rates Γ i = (cid:80) j (cid:54) = i W ji , of escape fromindividual traps, with finite c -corrections scaling as 1 /c : λ i = − Γ i + O (1 /c ) (17)To see this, note from Eqs. (6,7) that one can decomposethe symmetrized master operator M s as M s = M (0) + 1 c M (1) (18)where M (0) is a diagonal matrix with elements M (0) ij = − δ ij c (cid:88) k (cid:54) = i A ki − β ( E k − E i ) = − δ ij Γ i , (19)which are nothing other than the negative escape rates.The off-diagonal terms are collected in M (1) , which haselements (recall that A ij = 0 if i = j ) M (1) ij = A ij β ( E i − E j ) / . (20)We now exploit the fact that the escape rates in the“baseline” operator M (0) remain of order unity for c →∞ , while the contribution from M (1) in (18) scales with1 /c . Using standard perturbation theory [40, 41], theeigenvalues of M s can then be expanded as λ i = − Γ i + 1 c λ (1) i + 1 c λ (2) i + . . . , (21)with λ (1) i = M (1) ii = 0 , λ (2) i = (cid:88) j (cid:54) = i ( M (1) ij ) − Γ i + Γ j , (22)The sum defining λ (2) i has c entries whose size is indepen-dent of c , hence λ (2) i = O ( c ) and we get the 1 /c scaling ofthe corrections announced in (17). This implies in par-ticular that λ i → − Γ i in the mean field limit c → ∞ , sothat the escape rate distribution directly determines thespectrum of relaxation rates.The above discussion implies that there is in fact a sec-ond limit where relaxation rates and escape rates becomeidentical, namely T → c . This is simplybecause in that limit β → ∞ and so M (1) →
0. We willsee below that the distribution of (relaxation or escape)rates then degenerates into a sum of c + 1 delta peaks. B. Escape rate distribution
Given the result of the previous subsection, we nextconsider the distribution of escape rates. For the sake ofgenerality we do this for arbitrary c to start with. Theescape rate distribution can then formally be written as ρ Γ (Γ) = (cid:104) δ (Γ − ˆΓ) (cid:105) { E,E ,...,E c } (23)withˆΓ( { E , . . . , E c } , E ) = 1 c c (cid:88) j =1
11 + e − β ( E j − E ) , (24)The average in (23) is over the exponential distributionfrom which each of the c + 1 trap depths is independentlydrawn. The distribution ρ (Γ) can then be constructednumerically simply by the appropriate sampling. Ana-lytically, the distribution can be obtained explicitly forthe (unphysical) case c = 1, with the result (see detailsin Appendix A): ρ Γ (Γ) = T (cid:26) Γ T − (1 − Γ) − T − for Γ < / − T − (1 − Γ) T − for Γ > / . (25)This expression shows that for Γ →
0, the escape ratedensity diverges as Γ T − . As demonstrated in sec-tion V C below (cf. Eq. (46)), this divergence in factcontrols the small Γ-behaviour of ρ Γ (Γ) for any finite c .Remarkably, the exponent T − c has activated features when we lookat low rates, i.e. deep traps.By way of preparation for the discussion of the relax-ation rate spectra, we mention that for c ≥ /c between these two ex-treme values (see also Appendix D). For the case c = 2one finds, for example, that the relaxation rate densityat Γ = 1 / ρ Γ (Γ = 1 / ∼ (cid:90) / (1 − R ) − T − R − T +21 d R , (26) - 1.0 - 0.5 SamplingTheory
FIG. 1: Escape rate distribution obtained by analyticalcalculation (see Appendix A) and direct evaluation bynumerical sampling of (23) for c = 2 and T = 0 . T < /
3. This gives some intu-itive justification for the fact that similar divergences at − λ = i/c , i = 1 , . . . , c − c = 2 and T = 0 . , / c → ∞ ; we defer a quantitative evaluationof this mean field limit to the next section. In the othercase where the two distributions coincide, i.e. T → i = 1 c (cid:88) k (cid:54) = i A ki Θ( E k − E i ) (27)where Θ( · ) is the Heaviside step function. Each Γ i is thenjust the fraction of deeper ( E k > E i ) traps among theneighbours of traps i . This fraction can take the values i/c with i ∈ { , , . . . , c } . If we order the c + 1 depthsof a given trap and its c neighbours into a descendinglist, then a relaxation rate of i/c results precisely whenthe central trap is at position i + 1 in this list, becauseit then has i lower-lying (deeper) neighbours. But bypermutation symmetry the central trap is equally likelyto be in any position in the sorted list of c +1 trap depthsso that ρ Γ (Γ) = 1 c + 1 c (cid:88) i =0 δ (Γ − i/c ) (28)The above argument for the uniform prefactor 1 / ( c +1) ofthe delta peaks can be confirmed by explicit calculationas sketched in Appendix B. () T = 0.2T = 0.35T = 0.5T = 0.8RW
FIG. 2: T -dependence of the spectral density for c = 5,evaluated from the cavity theory (see Eq. (16)) withsmoothing parameter (cid:15) = 10 − . Also shown is thelimiting behaviour for T → ∞ as given by the randomwalk spectrum (Kesten–McKay law (29)). C. Relaxation rate distribution
We now turn to the relaxation rate spectrum, whichencodes information about collective relaxation modes.In contrast to the escape rate distribution, it thereforecannot generally be found from just local information.The spectrum depends on two key parameters, temper-ature T and network connectivity c . We analyse theireffects separately, beginning with the former. Changes in T . Using the cavity method presented insection III we can obtain the variation of the spectral den-sity with temperature for given c in the sparse networkregime ( c < ∞ ). Exemplary results for c = 5 are shownin Figure 2. They display a rich structure, with multiplepeaks in the spectrum at low T , but as we now demon-strate the qualitative behaviour can be understood fromthe two extreme cases of high and low T . For T → ∞ , thetrap depths are irrelevant as the Glauber transition ratebetween any two connected traps energy becomes 1 / (2 c ).The dynamics thus becomes a random walk and the spec-tral density depends only on the network structure. Forlarge random regular graphs the resulting spectral den-sity ρ ( λ ) is given explicitly by the Kesten-McKay law forinfinite regular trees [42–45], which with our choice oftransition rates reads ρ RW ( λ ) = 2 cπ (1 − λ + 1 / ) (cid:115) c − c − (cid:18) λ + 12 (cid:19) (29)and is also plotted in Figure 2. For T →
0, on the otherhand, the escape rate analysis presented in the previoussection shows that the relaxation rate spectral densityconsists of a sum of delta peaks of equal height at − λ = i/c , cf. Eq. (28). Qualitatively, the results in Figure 2can therefore be understood as interpolating between therelatively flat Kesten-McKay spectrum at high T anda series of peaks at low T . An obvious question then concerns the shape of these peaks at low but nonzero T ;we next show that they are power law divergences, withexponents that depend on T and on the position of thepeak as we already saw for the escape rate distribution.We begin with the peak at λ = 0, which correspondsto the distribution of the slowest modes that govern theaging dynamics at T <
1. We find that the spectraldensity in this region grows as ρ ( λ ) (cid:39) κ c | λ | T − (30)with a prefactor that for large c scales as κ c ∼ c T − (31)so that overall ρ ( λ ) ∼ | cλ | T − for large c and small λ .Specifically one requires λ (cid:28) /c to see this power lawscaling, as one can show by an analysis similar to the onein Appendix D. The scaling (30) may be obtained analyt-ically from the cavity equations as shown in Appendix F.This includes the prefactor κ c and its scaling (31) forlarge c . Figure 5 shows the evaluation of this prefactoracross a range of c , and clearly confirms the theoreti-cal predictions. The scaling of ρ ( λ ) with λ , though notthe prefactor, can also be derived within a Single DefectApproximation as explained in Appendix G.To summarize thus far, for any finite c the power lawdependence of ρ ( λ ) for small λ is the same as in theBouchaud trap model on a RRG [19]. We thus concludethat also the collective relaxation modes are governed byactivated processes at long times. The physical pictureis that the dynamics of the system is then dominated byslow transitions between local energy minima, i.e. trapsthat have no “escape directions” in the form of lower-lying neighbours. We validate the predicted λ -scaling inFigure 3, for fixed T and different c .The low temperature behaviour of the other spectralpeaks, around λ ∗ i = − i/c with i = { , . . . , c } , is moreintricate. Numerical evaluation using the cavity theorysuggests that all peaks (except the one near λ ∗ c = − ρ ( λ ) ∼ | λ + i/c | − x i . The exponent x i depends on thepeak i being considered and in general also on c . Itdecreases with temperature T ; where it drops to zerothe corresponding peak turns into a maximum of finiteheight. The temperatures where this occurs all lie in therange 0 < T ≤ /
2. In Figure 4, the spectrum for c = 3at low temperature is displayed together with the powerlaw fit for the intermediate peaks. For the first peak at λ ∗ = − /c , our numerical data for the exponent x acrossa range of temperatures are consistent with x = 1 − T independently of c , while the exponents of higher peaks( x etc.) do have a nontrivial c -dependence.The existence of intermediate peaks in the spectraldensity is characteristic of the Barrat–M´ezard trap modeland has no analogue in the Bouchaud model [19]. Thesepeaks reflects the network structure and hence entropiceffects; on the other hand the broadening of them reflects () c = 5c = 10c = 20c = 100| | T MF FIG. 3: Log-log plot of spectral density for T = 0 . c , as obtained from the cavitymethod with (cid:15) = 10 − . The power law asymptote forsmall | λ | is consistent with the theoretical prediction ∼ | λ | T − (dashed line) (30). The mean field limit whichcan be obtained from the numerical inversion of (34) isalso shown (dotted line). The plateau for small | λ | ispredicted by (36). () x = 1 2 Tx = 1 3 T FIG. 4: Spectral density for c = 3 and T = 0 . (cid:15) = 10 − . At λ ∗ i = − i/c thespectrum diverges as a power law with the exponentsindicated in the legend.activation effects. They can be seen as an extension tononzero T of the delta peaks at T = 0, where the dynam-ics is driven purely by entropic effects, and then must beinterpreted as a signature of the entropic barriers thatare present in the model and dominate the dynamics for T < / Changes in c . Having so far focused on the dependenceof the spectral density on T , we now keep T fixed andconsider the variation with c . As an example, Figure 3shows results for T = 0 .
2. For small c we observe spec-tral peaks as before. These become more numerous withincreasing c but also less pronounced, merging for verylarge c into a smooth spectrum. This mean-field limitcan be obtained explicitly; to our knowledge even this isa new result. We exploit the discussion in section IV A,which showed that for c → ∞ the relaxation rate distri- c c T =0.2 T =0.35 c T c c T =0.2 T =0.35 c T FIG. 5: Log-log plot of prefactors of the spectra in thesmall rate regime against connectivity c , for twodifferent temperatures T . Top: Escape rates. Themarkers are obtained by numerical sampling of (23) inthe range Γ = 10 − to 10 − . The dot-dashed linesshow the theoretically predicted prefactor (47). Thedashed lines show the asymptotic power lawscaling (46) with c . Bottom: Relaxation rates. Themarkers are obtained from numerical evaluation of thespectrum at − λ = 10 − for T = 0 . − λ = 10 − for T = 0 .
35. The dot-dashed lines show the theoreticallypredicted prefactor (30), evaluated as explained inAppendix F. The dashed lines show the asymptoticpower law scaling (31) with c .bution is identical to the escape rate distribution. Thelatter can be obtained by noting that for large c the es-cape rate (24) becomes a deterministic function of thetrap depth:Γ( E ) = lim c →∞ c c (cid:88) j =1
11 + e − β ( E j − E ) (32)= (cid:90) ∞ d E (cid:48) ρ E ( E (cid:48) )1 + e − β ( E (cid:48) − E ) (33)= F (1 , T, T, − e βE ) , (34)with F the Gauss hypergeometric function. From thisrelation the escape rate distribution ρ (Γ) is then found bya simple variable transformation, ρ (Γ) = ρ E ( E (Γ)) | d E dΓ | .In general the inverse function E (Γ) has to be found nu-merically. But for large E , i.e. deep traps with low escaperates, the lower limit of the integral (33) can be sent to −∞ with negligible error provided we are in the glassphase, T <
1, giving (see also [11])Γ( E ) ≈ πT sin( πT ) e − E , (35)As this is proportional to ρ E ( E ), it results in a constant escape rate density and hence also relaxation rate density, ρ ( λ ) = πT sin( πT ) , λ → . (36)in other words the spectrum is flat for small λ . Of courseas we saw above, for any finite c the spectral density even-tually has to cross over to the divergence (30) for small λ . Comparing with (36) shows that the crossover pointmust scale as − λ ∼ /c and so moves towards λ = 0 asthe mean field limit of large c is approached. This physi-cally means that the connectivity of the network sets thescale for which the crossover from entropic to energeticdynamics is observed (more on this in the next section).These features can be seen qualitatively by looking backat Figure 3. V. RESULTS: TIME DOMAIN PROPERTIESA. Numerical simulation approach
For time domain properties we obtain stochastic simu-lation results directly for the thermodynamic limit N →∞ , by generating an effectively infinite tree on the fly,during the course of a stochastic simulation using theGillespie algorithm. The network construction methodis identical to that used for the Bouchaud trap modelin [22]; its key advantage is that it allows us to studytime-dependent properties without any finite size effects.The Gillespie simulation algorithm consists of repeatedevaluation the following steps:1. For a given node i compute the total exit rate Γ i = (cid:80) j ∈ ∂i W ji .2. Compute the waiting time ∆ t > p i (∆ t ) = Γ i exp( − Γ i ∆ t ).3. Select the node i new to which the transition occursrandomly among the c neighbours of i , with prob-ability W ji / Γ i for node j .4. Return i new as the new node and increment therunning simulation time by ∆ t .We store the nodes visited as a function of time formany different realizations. This allows us to see the evo-lution of the mean energy (cid:104) E ( t ) (cid:105) and also obtain two timequantities, namely the correlation C (0 , t ) and the per-sistence P (0 , t ), that we can compare with results fromthe relaxation rate and escape rate spectra. All of thesequantities are evaluated with respect to the initial statefor the dynamics; we assume that this is a uniform dis-tribution across traps, corresponding to an equilibriumdistribution at large temperature ( T → ∞ ). B. Mean energy
Probably the simplest manifestation of the aging dy-namics of trap models in the glass regime is the fact thatthe average energy decreases continually with time; inthis sense the system keeps track of its age through its t -7-6-5-4-3-2 E ( t ) t * c T= 0.2T= 0.35T= 0.5Bouchaud () | * | 1/ c FIG. 6: Mean energy as a function of time for c = 15,averaged over 10 different trajectories. Inset: spectrumof relaxation rates for the same temperatures. Thevertical line in both figures indicates the existence of acrossover regime beyond which (for large t , or small | λ | )activated behavior is found. t -6-5-4-3-2 E ( t ) c= 5c= 10c= 20c= 100BMBouchaud t / c -2-1012 E ( t ) + l n ( c ) FIG. 7: Mean energy as a function of time for T = 0 . different trajectories. The dashedlines correspond to the behavior expected from themean field Bouchaud model and the mean field Barrat–M´ezard (BM) model. Inset: Collapse to a master curvefor large c after rescaling by the crossover time t ∗ ∼ c ,see text for details.energy [46]. In Figure 6 we show what form this energydecay takes in the BM model with finite connectivity c , for three different temperatures. Beyond an initialtransient the mean energy decay is logarithmic in time,but with a clear change in slope at some crossover time.From the spectral point of view, this crossover time t ∗ can be understood in terms of a crossover relaxation rate | λ ∗ | ∼ /t ∗ below which the spectral density ρ ( λ ) followsthe activated behavior (30). This crossover is indicatedin the inset of Figure 6.The notion of a change in the dynamics between anentropic, non–activated and an activated regime is alsoapparent from Figure 7, where the dependence on c of theenergy decay is displayed together with the predicted be- - 1 - 2.5 - 2.0 - 1.5 - 1.0 - 0.5 T = 0.2T = 0.510 - 1 - 2.5 - 2.0 - 1.5 - 1.0 - 0.5 FIG. 8: Correlation and persistence function for c = 5and two different T . Colored lines: Stochasticsimulation results averaged over 10 trajectories; thepersistence is the lower curve in each pair of curves.Dash-dotted lines: Cavity theory predictions forcorrelation using (41). Dashed lines: Analyticalprediction for persistence (44).havior from the mean-field Bouchaud and Barrat–M´ezardmodels, which is [4, 7] (cid:104) E ( t ) (cid:105) ∼ (cid:40) − T ln( t ) Bouchaud , − ln( t ) Barrat M´ezard . (37)It is clear that the crossover between these two regimesshifts towards longer times as c increases. This of coursehas to be so, as for c → ∞ the mean-field BM behaviourmust be obtained for any finite t . In quantitative terms,it follows by comparison of the small λ divergence (30)of the spectrum with the mean-field plateau (36) thatthe crossover relaxation rate λ ∗ has to scale with 1 /c andtherefore t ∗ ∼ c . We confirm this in the inset of Figure 7,where we show the energy decay plotted against t/c ∼ t/t ∗ . Applying the corresponding shift − E → − E +ln c to the energy axis, which accounts for the entropic(BM) relaxation up to t ≈ t ∗ (see Eq. (37)), shows arather good collapse of the curves for different values ofthe connectivity c . C. Two-time observables
In this final section we study the correlation C (0 , t ) andpersistence P (0 , t ) (or survival probability) as functionsof the observation time t ; the initial time will always befixed at “waiting time” t w = 0. The correlation is definedas the probability of finding the system in the same trapat time 0 and t . If the system starts in trap i , then fromthe master equation (2) and its solution (4) this proba-bility is the ii -element of the propagator, ( e M t ) ii . Av-eraging over the uniform initial distribution across traps yields C (0 , t ) = 1 N N (cid:88) i =1 (e M t ) ii = 1 N Tr e M t = 1 N N − (cid:88) α =0 e λ α t , (38)In the limit of large N , Eq. (38) becomes C (0 , t ) = (cid:90) d λ ρ ( λ )e λt (39)so that the correlator can be deduced from the spec-tral density that we can already predict using the cav-ity method. An essentially equivalent connection can bemade to the trace of the resolvent: from (39), the Laplacetransform of the correlator is L [ C (0 , t )] = 1 N N − (cid:88) α =0 s − λ α = 1 N Tr ( s I − M ) − (40)As the eigenvalues of the master operator M and its sym-metrized version M s are identical, the last expression isthe normalized trace of the resolvent (8), giving C (0 , t ) = L − (cid:34) N N (cid:88) i =1 G ii ( s ) (cid:35) , (41)The resolvent trace can again be obtained using the cav-ity method and the correlator obtained by numerical in-verse Laplace transform. This is numerically more conve-nient than finding the spectrum first and then apply (39).Figure 8 compares the results with averages over 10 stochastic simulation trajectories. The agreement is ex-cellent, providing an explicit demonstration of how spec-tral information can be used to obtain non-trivial tem-poral properties of the BM trap model dynamics.We turn next to the persistence function P (0 , t ). Thisis defined as the probability of not leaving the initial trapup to time t . If the initial trap is i , this probability canbe expressed in terms of the escape rate Γ i of the trap ase − Γ i t . Averaging again over the initial traps gives P (0 , t ) = 1 N N (cid:88) i =1 e − Γ i t = (cid:90) dΓ ρ Γ (Γ)e − Γ t , (42)where the last expression applies in the limit of large N and relates the persistence to the escape rate distribu-tion, by direct analogy with (39). Inserting the expres-sions (23,24) for this distribution gives explicitly P (0 , t ) = (cid:104) e − ( t/c ) (cid:80) cj =1 / [1+exp( − β ( E j − E ))] (cid:105) E,E ,...,E c (43)= (cid:90) d E ρ E ( E ) (cid:20)(cid:90) d E ρ E ( E )e − tc / [1+exp( − β ( E − E )] (cid:21) c . (44)This integral can be evaluated numerically and the resultcompared to that calculated by averaging over stochastic0simulation trajectories, as shown in Figure 8 alongsidethe correlation results. As it should be the agreement isvery good. The long-time scaling of the persistence canalso be worked out analytically from the integral (44) asdiscussed in Appendix C, giving P (0 , t ) (cid:39) π c t − T (45)with π c ∼ c T − for large c . By computing the inverseLaplace transform in (42) we obtain that ρ Γ (Γ) (cid:39) ˜ π c Γ T − (46)where (with Γ( · ) the Euler Gamma function)˜ π c = π c / Γ( T ) (47)The comparison with (30) reveals that, for a given tem-perature, the spectra for the slowest processes (eitherof escape or relaxation) have the same scaling. Thisis clearly supported by the numerical data in Fig. 5.We remark that the correspondence between ρ ( λ ) and ρ Γ (Γ) translates in the time domain into a proportional-ity between C (0 , t ) and P (0 , t ) at long times, which againagrees with our numerics (see Fig. 8). VI. CONCLUSIONS
We have focussed in this study on the characteriza-tion of the dynamics of the Barrat–M´ezard trap model,which is a simple coarse-grained model for the dynamicsof glassy systems in configuration space. Its main featureis the combination of activated and non–activated path-ways on an energy landscape consisting of traps with anexponential distribution of depths. Key for the impor-tance of activated processes is the sparse connectivity ofthe network of traps in configuration space; we modelledthis as a random regular graph. This is more general(and realistic) than the original Barrat–M´ezard model,which is of mean-field type because its network is fullyconnected [3].Mathematically the model is a continuous time Markovchain on a random graph with microscopic transitionrates given by the Glauber form (1). The aim of thiswork was to obtain the spectral density of the associ-ated master operator, in the thermodynamic limit. Thiswas done via the Cavity Method. Our main finding isthat below the glass transition temperature the slowestrelaxation modes, which determine the long time be-havior, have a density that follows a divergent powerlaw (Eq. (30)), which is characteristic of activated pro-cesses [19, 22]. This contrasts with the mean field case ofinfinite connectivity ( c → ∞ ), where the spectral densityis constant for the slowest modes. We established this byshowing that in this limit the relaxation rates reduce tolocal escape rates. This allowed us to calculate the spec-trum exactly for the mean-field BM trap model, whichto the best of our knowledge is also a new result. At low temperatures and finite connectivity we found power lawsingularities also for intermediate relaxation rates. Weinterpreted these as the finite temperature analog of thedelta peaks that make up the spectrum for T = 0. Thesedelta peaks indicate that relaxation rates – again equalto escape rates for T → λ -divergence (30) of thespectral density for finite c with the mean-field plateauwe obtained a crossover relaxation rate | λ ∗ | ∼ /c . Phys-ically, this means that the network structure generatesa crossover between entropically dominated dynamics at t < / | λ ∗ | ∼ c to primarily activated dynamics for longertimes. Numerical simulations data for the decay of theaverage energy, which we obtained via a bespoke algo-rithm that eliminates all corrections due to finite networksize, confirmed this behaviour.The correlation and persistence functions following aquench from infinite temperature were also analyzed; weshowed that these can be predicted from the spectrumof relaxation and escape rates, respectively, in excellentagreement with numerical simulation data.In studying the BM trap model on networks withsparse connectivity, our broader aim was to make trapmodels into more accurate representations of real glassysystems. Indeed, the fact that the standard trap mod-els have a mean field character has been criticized in anumber of studies, and more realistic descriptions havebeen demanded [12, 16]. Our work is designed to fill thisgap. What is fascinating is the natural appearance ofa crossover in the glass regime between slow dynamicsthat are initially governed by entropy barriers but laterbecome dominated by activation across energy barriers,directly from the sparse network connectivity. In partic-ular the crossover emerges without the imposition of a re-stricted network size or of changes in the dynamical ruleswith trap depth [7]. Interestingly the entropic-energeticcrossover does not appear in Bouchaud trap models onsparse networks [19, 22], presumably because the acti-vated form of the local transition rates there rules outany downhill motion between traps for which entropicbarriers would be relevant.We comment briefly on the scaling of c in connection toreal-space glass physics. The number N of local energyminima and hence traps is expected to scale exponen-tially in system volume V . The typical number of con-figurations reachable from a trap i , on the other hand,will scale as V because transitions between traps will cor-respond to roughly independent local particle rearrange-ments inside finite volume elements. The overall rate oftransitions must then also scale as V ∼ c . Our O (1)transition rates would have to be scaled up by a factor c to accommodate this, and timescales scaled down by1 /c accordingly. The entropic-energetic crossover time t ∗ then becomes independent of c , so does not divergewith system size. This suggests that it is likely to beobservable experimentally or in targeted simulations.Open questions for future work concern in particular1the aging dynamics of the BM trap model on sparse net-works. While for Bouchaud trap models the aging be-haviour of two-time correlation functions C ( t w , t ) ulti-mately appears to become independent of the connec-tivity [10, 17, 47], one would not expect this for theBM model where the crossover between entropic and en-ergetic barriers introduces a separate timescale. Thecomparison to the aging scalings resulting from otherapproaches for generating such crossover timescales [7]should be particularly revealing. The wider context of trap models as dynamical systems in a disordered po-tential also opens up further directions by making con-tact with different models in the area of disordered sys-tems. The most relevant for our purposes is the Andersonmodel, in particular on sparse networks where it is thesubject of much ongoing research [31, 48]. A full under-standing of Anderson localization in this case remainsopen and the connection with trap models may lead tonew results, for instance, regarding dynamical universal-ity classes. [1] J. P. Bouchaud. Weak ergodicity breaking and aging indisordered systems. Journal de Physique I , 2(9):1705–1713, September 1992.[2] J.-P. Bouchaud and D. S. Dean. Aging on Parisi’s tree.
Journal de Physique I , 5(3):265–286, March 1995. arXiv:cond-mat/9410022.[3] A. Barrat and M. Mzard. Phase Space Diffusion and LowTemperature Aging.
Journal de Physique I , 5(8):941–947,August 1995.[4] Ccile Monthus and Jean-Philippe Bouchaud. Models oftraps and glass phenomenology.
Journal of Physics A:Mathematical and General , 29(14):3847–3869, July 1996.[5] R. Mlin and P. Butaud. Glauber Dynamics and Ageing.
Journal de Physique I , 7(5):691–710, May 1997. arXiv:cond-mat/9701112.[6] Bernd Rinn, Philipp Maass, and Jean-PhilippeBouchaud. Multiple Scaling Regimes in Simple AgingModels.
Physical Review Letters , 84(23):5403–5406,June 2000.[7] Eric M Bertin. Cross-over from entropic to thermal dy-namics in glassy models.
Journal of Physics A: Mathe-matical and General , 36(43):10683–10700, October 2003.[8] R. Aldrin Denny, David R. Reichman, and Jean-PhilippeBouchaud. Trap Models and Slow Dynamics in Super-cooled Liquids.
Physical Review Letters , 90(2):025503,January 2003.[9] B. Doliwa and A. Heuer. Energy barriers and activateddynamics in a supercooled Lennard-Jones liquid.
PhysicalReview E , 67(3):031506, March 2003.[10] Gerard Ben Arous and Jiri Cerny. Dynamics of trapmodels. arXiv:math/0603344 , March 2006.[11] Peter Sollich. Trap models with slowly decorrelating ob-servables.
Journal of Physics A: Mathematical and Gen-eral , 39(11):2573–2597, March 2006.[12] Paolo Moretti, Andrea Baronchelli, Alain Barrat, andRomualdo Pastor-Satorras. Complex networks andglassy dynamics: walks in the energy landscape.
Jour-nal of Statistical Mechanics: Theory and Experiment ,2011(03):P03032, March 2011.[13] M Baity-Jesi, G Biroli, and C Cammarota. Activated ag-ing dynamics and effective trap model description in therandom energy model.
Journal of Statistical Mechanics:Theory and Experiment , 2018(1):013301, January 2018.[14] Chiara Cammarota and Enzo Marinari. Numerical ev-idences of universal trap-like aging dynamics.
Jour-nal of Statistical Mechanics: Theory and Experiment ,2018(4):043303, April 2018.[15] Eric Woillez, Yariv Kafri, and Nir Gov. The active trapmodel. arXiv preprint arXiv:1910.02667 , 2019. [16] Yasheng Yang and Bulbul Chakraborty. Dynamics inthe metabasin space of a Lennard-Jones glass former:Connectivity and transition rates.
Physical Review E ,80(1):011501, 2009.[17] Grard Ben Arous and Ji ern. The arcsine law as a uni-versal aging scheme for trap models.
Communications onPure and Applied Mathematics , 61(3):289–329, 2008.[18] Vronique Gayrard. Aging in reversible dynamics of dis-ordered systems. I. Emergence of the arcsine law inBouchaud’s asymmetric trap model on the completegraph. arXiv:1008.3855 [math-ph] , August 2010. arXiv:1008.3855.[19] Riccardo Giuseppe Margiotta, Reimer K¨uhn, and PeterSollich. Spectral properties of the trap model on sparsenetworks.
Journal of Physics A: Mathematical and The-oretical , 51(29):294001, 2018.[20] Anton Bovier, Alessandra Faggionato, et al. Spectralcharacterization of aging: the REM-like trap model.
TheAnnals of Applied Probability , 15(3):1997–2037, 2005.[21] Takuma Akimoto, Eli Barkai, and Keiji Saito. Non-self-averaging behaviors and ergodicity in quenched trapmodels with finite system sizes.
Physical Review E ,97(5):052143, May 2018. Publisher: American PhysicalSociety.[22] Riccardo Giuseppe Margiotta, Reimer K¨uhn, and PeterSollich. Glassy dynamics on networks: local spectra andreturn probabilities.
Journal of Statistical Mechanics:Theory and Experiment , 2019(9):093304, 2019.[23] Chiara Cammarota and Enzo Marinari. Spontaneousenergy-barrier formation in entropy-driven glassy dy-namics.
Physical Review E , 92(1):010301, 2015.[24] Tim Rogers, Isaac P´erez Castillo, Reimer K¨uhn, andKoujin Takeda. Cavity approach to the spectral densityof sparse symmetric random matrices.
Physical ReviewE , 78(3):031116, 2008.[25] Marc Mezard and Andrea Montanari.
Information,Physics, and Computation . Oxford University Press,2009.[26] Fernando Lucas Metz, Izaak Neri, and D´esir´e Boll´e.Localization transition in symmetric random matrices.
Physical Review E , 82(3):031135, 2010.[27] Ariel Amir, Yuval Oreg, and Yoseph Imry. On relaxationsand aging of various glasses.
Proceedings of the NationalAcademy of Sciences , 109(6):1850–1855, 2012.[28] Yoav Lahini, Omer Gottesman, Ariel Amir, andShmuel M. Rubinstein. Nonmonotonic Aging and Mem-ory Retention in Disordered Mechanical Systems.
Phys-ical Review Letters , 118(8):085501, February 2017. Pub-lisher: American Physical Society. [29] Jean-Philippe Bouchaud and Marc Mzard. Universalityclasses for extreme-value statistics. Journal of Physics A:Mathematical and General , 30(23):7997–8015, December1997. Publisher: IOP Publishing.[30] Rka Albert and Albert-Lszl Barabsi. Statistical mechan-ics of complex networks.
Reviews of Modern Physics ,74(1):47–97, January 2002. Publisher: American Physi-cal Society.[31] Giulio Biroli and Marco Tarzia. Delocalization and er-godicity of the Anderson model on Bethe lattices. arXivpreprint arXiv:1810.07545 , 2018.[32] Andrea De Luca, BL Altshuler, VE Kravtsov, andA Scardicchio. Anderson localization on the Bethe lat-tice: Nonergodicity of extended states.
Physical ReviewLetters , 113(4):046806, 2014.[33] M Sonner, KS Tikhonov, and AD Mirlin. Multifractalityof wave functions on a Cayley tree: From root to leaves.
Physical Review B , 96(21):214204, 2017.[34] KS Tikhonov and AD Mirlin. Statistics of eigen-states near the localization transition on random regulargraphs.
Physical Review B , 99(2):024202, 2019.[35] William J Anderson.
Continuous-time Markov chains:An applications-oriented approach . Springer Science &Business Media, 1991.[36] Reimer Khn. Spectra of sparse random matrices.
Journal of Physics A: Mathematical and Theoretical ,41(29):295002, June 2008. Publisher: IOP Publishing.[37] Giacomo Livan, Marcel Novaes, and Pierpaolo Vivo.
In-troduction to Random Matrices: Theory and Practice ,volume 26. Springer, 2018.[38] Ragi Abou-Chacra, DJ Thouless, and PW Anderson. Aselfconsistent theory of localization.
Journal of PhysicsC: Solid State Physics , 6(10):1734, 1973.[39] Numerical code implemented in Julia can be foundat https://github.com/dapias/SparseBarratMezard .[40] R´egis M´elin and P Butaud. Glauber dynamics and age-ing.
Journal de Physique I , 7(5):691–710, 1997.[41] David J Griffiths and Darrell F Schroeter.
Introduction toquantum mechanics . Cambridge University Press, 2018.[42] Harry Kesten. Symmetric random walks on groups.
Transactions of the American Mathematical Society ,92(2):336–354, 1959.[43] Fan RK Chung and Fan Chung Graham.
Spectral graphtheory . Number 92. American Mathematical Soc., 1997.[44] Brendan D McKay. The expected eigenvalue distributionof a large regular graph.
Linear Algebra and its Applica-tions , 40:203–216, 1981.[45] Roland Bauerschmidt, Jiaoyang Huang, and Horng-Tzer Yau. Local Kesten–McKay law for random regu-lar graphs.
Communications in Mathematical Physics ,369(2):523–636, 2019.[46] Eric Bertin and Fran¸cois Bardou. From laser coolingto aging: a unified L´evy flight description.
AmericanJournal of Physics , 76(7):630–636, 2008.[47] Riccardo Margiotta.
Glassy Dynamics on networks. Spec-tra, return probabilities and aging . PhD thesis, King’sCollege London, 2019.[48] KS Tikhonov, AD Mirlin, and MA Skvortsov. Ander-son localization and ergodicity on random regular graphs.
Physical Review B , 94(22):220203, 2016.[49] Jean-Philippe Bouchaud and Antoine Georges. Anoma-lous diffusion in disordered media: Statistical mecha-nisms, models and physical applications.
Physics Re-ports , 195(4-5):127–293, 1990.
Appendix A: Escape rate distribution
The escape rate distribution (23) is most straightfor-ward to analyze from its Fourier transform (cid:104) e is Γ (cid:105) . As (42)shows, this is the same as the persistence function P (0 , t )with t replaced by − is . Writing the exponential energydistributions in (44) explicitly thus gives (cid:104) e is Γ (cid:105) = (cid:90) d E e − E [∆( E, s/c )] c (A1)with∆( E, s/c ) = (cid:90) d E e − E e i ( s/c ) / [1+exp( − β ( E − E )] (A2)If we make the change of variable R = 11 + e − β ( E − E ) , (A3)then d R = βR (1 − R )d E and e − ( E − E ) = [(1 − R ) /R ] T ,yielding∆( E, s/c ) = T e − E (cid:90) R min d R (1 − R ) T − R T +1 e i ( s/c ) R (A4)The lower integration limit R min = (1 + e βE ) − lies be-tween 0 and 1 /
2, so for further calculation it can be help-ful to split the integral into ∆ = ∆ + + ∆ − , where ∆ + isthe integral for R = 1 / . . . − the remainder.Now we focus on the case c = 1, where from (A1) wejust need the averages of ∆ ± ( E, s ) over the exponentiallydistributed central trap depth E . The first of these isstraightforward: (cid:90) ∞ d E e − E ∆ + ( E, s ) == T (cid:90) ∞ d E e − E (cid:90) / d R e isR (1 − R ) T − R T +1 (A5)= T (cid:90) / d R e isR (1 − R ) T − R T +1 . (A6)For the average of ∆ − we note that the lower integrationlimit R min = (1 + e βE ) − on R corresponds to a lowerlimit of E min = T ln[(1 − R ) /R ] for E at fixed R : (cid:90) ∞ d E e − E ∆ − ( E, s ) == T (cid:90) ∞ d E e − E (cid:90) / R min d R e isR (1 − R ) T − R T +1 (A7)= T (cid:90) / d R e isR (1 − R ) T − R T +1 e − E min T (cid:90) / d R e isR (1 − R ) T − R T +1 (cid:18) − RR (cid:19) − T (A9)= T (cid:90) / d R e isR (1 − R ) − T − R − T +1 . (A10)3Altogether we have now (cid:104) e is Γ (cid:105) = (A6)+ (A10). Thefactors e isR in both terms just produce Dirac deltas δ (Γ − R ) upon inverse Fourier transform, giving directlythe result (25). The symmetry of the distribution underΓ → − Γ follows intuitively from the fact that the escaperate R is transformed to 1 − R when E − E changes sign,together with the fact that this trap depth difference hasan even distribution.For the case c = 2, the result in equation (26) is derivedby computing (∆ + + ∆ − ) , inverse Fourier transformingit and evaluating the remaining integrals for Γ = 1 / ≤ / ρ Γ (Γ) = 4 T (cid:90) Γ0 d x f (2Γ − x ) g ( x ) , (A11)whereas for 1 / < Γ ≤ / ρ Γ (Γ) = 4 T (cid:18) (cid:90) / − d x f (2Γ − x ) g ( x )++ (cid:90) − / / d x f ( x ) f (2Γ − x )2 (cid:19) , (A12)and finally for the remaining interval 3 / < Γ ≤ ρ Γ (Γ) = 2 T (cid:90) − d x f ( x ) f (2Γ − x )2 , (A13)with f ( x ) = (1 − x ) T − x T +1 and g ( x ) = (1 − x ) − T − x − T +1 . Theabove expression produces the plot in Fig. 1. Appendix B: Spectra for T = 0 We used a permutation symmetry argument in themain text to justify why in the T = 0 spectrum (28)each delta peak has the same prefactor. This can alsobe seen explicitly as follows. The prefactor a k for each δ (Γ − k/c ) in equation (28) is the probability of having k traps among the c neighbours of a given trap that lielower, i.e. have a greater depth, and c − k traps that liehigher. Calling the depth of the given central trap E ,this gives a k = (cid:18) ck (cid:19) (cid:90) d E ρ E ( E ) P ( . . . , E k > E, E k +1 < E, . . . )(B1)= (cid:18) ck (cid:19) (cid:90) d E ρ E ( E ) (cid:18) (cid:90) ∞ E d E (cid:48) ρ E ( E (cid:48) ) (cid:19) k ×× (cid:18) (cid:90) E d E (cid:48) ρ E ( E (cid:48) ) (cid:19) c − k (B2)= (cid:18) ck (cid:19) ( c − k )! k !( c + 1)! = 1 c + 1 , (B3) where in the initial integrand we assumed a specific order-ing of the lower and higher neighbours and compensatedfor this by the binomial coefficient prefactor. The inte-gral in (B2) with the variable change q = (cid:82) E d E (cid:48) ρ E ( E (cid:48) )evaluates to a Beta function as used in the line below. Appendix C: Scaling of persistence function
We derive here the large t -scaling of the persistencefunction, which from Eq. (44) is: P (0 , t ) = (cid:90) d E e − E (cid:18)(cid:90) d E e − E e − ( t/c ) / (1+exp( − β ( E − E )) (cid:19) c . (C1)For t/c (cid:29)
1, the denominator 1 + e − β ( E − E ) in the ex-ponent must be large to get a significant contributionand can therefore be approximated by e − β ( E − E ) . Thecentral trap depth E then appears in the combination( t/c )e − βE . This suggests a change of integration vari-able to ω = ( t/c ) T e − E . With a similar transformation q = e − E /ω for E , Eq. (C1) becomes P (0 , t ) = (cid:16) ct (cid:17) T (cid:90) ( t/c ) T d ω (cid:34) ω (cid:90) /ω d q e − / ( c/t + q β ) (cid:35) c (C2)This is still exact but can be simplified for large times,specifically t/c (cid:29)
1. The upper boundary in the outerintegral can then be replaced by ∞ and the c/t in theintegrand can be neglected, giving the asymptotic scaling P (0 , t ) (cid:39) π c t − T (C3)The prefactor π c = c T (cid:90) ∞ d ω (cid:32) ω (cid:90) /ω d q e − q − β (cid:33) c (C4)is, at fixed T , just a function of the connectivity c . Tounderstand its scaling for large c , note that the innerintegral is taken to the power c and so contributes onlyin the small ω -region where it is 1 − O (1 /c ). We thereforewrite it as 1 − ω (cid:90) /ω d q (cid:16) − e − q − β (cid:17) (C5)For small ω the upper integration boundary can again betaken to ∞ , giving 1 − d T ω up to higher order corrections,with d T = (cid:90) ∞ d q (cid:16) − e − q − β (cid:17) = Γ(1 − T ) (C6)where Γ( · ) is the Euler Gamma function. The prefac-tor (C4) of the persistence for large c is then π c = c T (cid:90) ∞ d ω (1 − d T ω ) c = c T (cid:90) ∞ d ω e − cd T ω = d − T c T − (C7)4giving overall for large c and large times t ≥ cP (0 , t ) (cid:39) d − T c T − t − T (C8)which is the scaling announced in (45) in the main text.We note that the arguments above can be extended tounderstand the entire large c -scaling of the persistencearound the entropic-energetic crossover, where ˜ t = t/c is of order unity. Anticipating that again small ω willdominate, we rescale ω = ˜ ω/c in (C2) to get P (0 , t ) = ˜ t − T c (cid:90) c ˜ t T d˜ ω (cid:40) − ˜ ωc (cid:90) c/ ˜ ω d q (cid:104) − e − / (˜ t − + q β ) (cid:105)(cid:41) c (C9)For large c the upper integration boundaries again tendto ∞ , while the integrand { . . . } c becomes an exponential,leading to P (0 , t ) = ˜ t − T c (cid:90) ∞ d˜ ω exp (cid:18) − ˜ ω (cid:90) ∞ d q (cid:104) − e − / (˜ t − + q β ) (cid:105)(cid:19) (C10)= ˜ t − T c (cid:18)(cid:90) ∞ d q (cid:104) − e − / (˜ t − + q β ) (cid:105)(cid:19) − . (C11)This shows that for large c , cP (0 , t ) does indeed becomea function only of the time ˜ t scaled to the crossover time t ∗ ∼ c , so it is convenient to introduce the scaled persis-tence ˜ P (0 , ˜ t ) ≡ cP (0 , ˜ t/c ) . (C12)For large ˜ t , Eq. (C11) directly retrieves the scaling (C8).For small ˜ t , on the other hand, the exponential can belinearized so that the q -integral becomes (cid:90) ∞ d q t − + q β = πT sin( πT ) ˜ t − T (C13)and hence ˜ P (0 , ˜ t ) = sin( πT ) πT t , (C14)The unscaled persistence is then P (0 , t ) = ˜ P (0 , t/c ) c = sin( πT ) πT t (C15)This exhibits the expected 1 /t decay in the entropicallydominated regime, where well before the crossover theconnectivity c is irrelevant as long as it is large enough. Appendix D: Scaling of escape rate distribution forlarge c We show in this appendix what the scaling (C12) ofthe persistence function P (0 , t ) = c − ˜ P (0 , ˜ t ) , ˜ t = t/c (D1) for large c implies for the escape rate distribution ρ Γ (Γ).Using the relation (42), the rescaled persistence can bewritten as˜ P (0 , ˜ t ) = c (cid:90) dΓ ρ Γ (Γ)e − Γ c ˜ t = (cid:90) d˜Γ ρ Γ (˜Γ /c )e − ˜Γ˜ t (D2)where ˜Γ = c Γ. For this to have a limit for large c requiresthat also the rescaled relaxation spectrum˜ ρ (˜Γ) = ρ Γ (˜Γ /c ) (D3)must become independent of c . This master curve for thelarge c -relaxation rate spectrum can expressed in termsof an infinite series using the following steps, startingfrom (C11). First, transform q → ˜ q = q ˜ t T . Second,introduce u = 1 / (1 + ˜ q β ). This leads to˜ P (0 , ˜ t ) = (cid:18) T (cid:90) d u (1 − u ) T − u − T − (1 − e − ˜ tu ) (cid:19) − . (D4)We now extract the dominant large ˜ t -term by decompos-ing the integral into three parts:˜ P (0 , ˜ t ) − = T (cid:90) ∞ d u u − T − (1 − e − ˜ tu ) − T (cid:90) d u u − T − [1 − (1 − u ) T − ](1 − e − ˜ tu ) − T (cid:90) ∞ d u u − T − (1 − e − ˜ tu ) . (D5)The first integral can be computed analytically and givesa pure power law as intended, while the last two termscan be combined into a single integral:˜ P (0 , ˜ t ) − = Γ(1 − T )˜ t T + (cid:90) ∞ d u π ( u )(1 − e − ˜ tu ) , (D6)where Γ( · ) is the Euler Γ-function and π ( u ) is defined as π ( u ) = − T u − T − (cid:40) − (1 − u ) T − u < u ≥ . (D7)One can check that (cid:82) ∞ d uπ ( u ) = 0, giving the furthersimplification˜ P (0 , ˜ t ) − = Γ(1 − T )˜ t T (cid:18) − ˜ t − T Γ(1 − T ) (cid:90) ∞ d u π ( u )e − ˜ tu (cid:19) , (D8)For large ˜ t , the inverse can now be expanded into a geo-metric series˜ P (0 , ˜ t ) = ˜ t − T Γ(1 − T ) (cid:18) − ˜ t − T Γ(1 − T ) (cid:90) ∞ d u π ( u )e − ˜ tu (cid:19) − (D9)= ˜ t − T Γ(1 − T ) (cid:88) n ≥ ˜ t − nT Γ n (1 − T ) (cid:90) ∞ d u π ∗ n ( u )e − ˜ tu (D10)5 () c = 50c=100Theory FIG. 9: Relaxation rate spectrum plotted against thescaled rate ˜Γ = c Γ for T = 0 .
2. Numerical data for twofinite c are shown against the theoretical master curvefor c → ∞ , as predicted by equation (D11). The serieswas truncated beyond n = 3, which makes it inaccuratefor ˜Γ ≥ ∗ n denotes the n -th convolution.This expression can now be conveniently inverse Laplacetransformed (see (D2,D3)) to get the scaled relaxationrate spectrum˜ ρ (˜Γ) = ˜Γ T − Γ( T )Γ(1 − T ) + (cid:88) n ≥ (cid:82) ˜Γ0 d u π ∗ n ( u )(˜Γ − u ) ( n +1) T − Γ(( n + 1) T )Γ n +1 (1 − T ) . (D11)The first few terms of this series are straightforward toevaluate numerically. In Figure 9 we compare the result-ing prediction for the master curve with numerical datafor two different (large) connectivities, finding very goodagreement. The integral associated with n = 1 in (D11)controls the behavior at ˜Γ = 1; it is explicitly given by (cid:90) d u u − T − (1 − (1 − u )) T − (˜Γ − u ) T − (D12)It is remarkable that this diverges for T < /
3, which isexactly what we found in the case c = 2 for the corre-sponding peak at Γ = 1 / c , and the same may be true for the peaksat ˜Γ = 2 , , . . . – of course only for large enough c as˜Γ ≤ c generally. Appendix E: Cavity equations: Mean field and zerotemperature limit
In this section we show that the results for the meanfield limit (36) and zero temperature (the analogueof (28)) can be obtained via the cavity equations (11,12).We start by rewriting the equations in a way thatmakes them simpler to analyse. We start by divid-ing both sides of (12) by the factor e βE k c and cor-respondingly define rescaled cavity precisions ˜ ω ( j ) k = ω ( j ) k e − βE k /c :˜ ω ( j ) k = i ( λ − i(cid:15) ) + (cid:88) l ∈ ∂k \ j i [ K ( E k , E l )e − βE k /c ] ˜ ω ( k ) l i [ K ( E k , E l )e − βE l /c ] + ˜ ω ( k ) l (E1)From (7) the combinations in square brackets just givetransition rates:˜ ω ( j ) k = i ( λ − i(cid:15) ) + (cid:88) l ∈ ∂k \ j iW lk ˜ ω ( k ) l iW kl + ˜ ω ( k ) l (E2)The equations for the scaled marginal precisions followin the same way from (11), giving˜ ω j = i ( λ − i(cid:15) ) + (cid:88) k ∈ ∂j iW kj ˜ ω ( j ) k iW jk + ˜ ω ( j ) k (E3)Bearing in mind that the transition rates scale as 1 /c ,the sum in (E2) and hence the typical cavity precisionis O (1). The transition rates in the denominators of(E2,E3) can thus be neglected for large c and (E3) sim-plifies to ˜ ω j = i ( λ − i(cid:15) ) + (cid:88) k ∈ ∂j iW kj , (E4)= i ( λ − i(cid:15) ) + i Γ j . (E5)The spectral density becomes (cf. Eqns (13), (16)): ρ ( λ ) = lim (cid:15) → πN N (cid:88) j =1 Re(1 / ˜ ω j ) (E6)= lim (cid:15) → πN N (cid:88) j =1 Re (cid:18) i ( λ + Γ j ) + (cid:15) (cid:19) (E7)= (cid:104) δ ( λ + Γ) (cid:105) . (E8)Thus in the mean field limit the distribution of escaperates becomes equal to the distribution of relaxationrates.For the case T = 0, on the other hand, the cavityequation (E3) becomes˜ ω j = i ( λ − i(cid:15) ) + 1 c (cid:88) k ∈ ∂j i Θ( E k − E j ) , (E9)giving for the spectral density ρ ( λ ) = lim (cid:15) → πN N (cid:88) j =1 Re (cid:32) i [ λ + (cid:80) k ∈ ∂j Θ( E k − E j ) /c ] + (cid:15) (cid:33) (E10)= (cid:42) δ λ + (cid:88) k ∈ ∂j Θ( E k − E j ) c (cid:43) (E11)= 1 c + 1 c (cid:88) k =0 δ (cid:18) λ + kc (cid:19) , (E12)where in the last line we have used the result from Ap-pendix B.6 Appendix F: Cavity equations: Small λ -limit andlarge c -limit In this section we will obtain the small λ limit (30) forthe spectral density from the cavity equations (11,12).We will do this first for finite c and then show how thelarge c -behaviour of the prefactor (31) can be extracted.We start by writing (E2) explicitly as˜ ω ( j ) k = i ( λ − i(cid:15) ) + i e − βE k (cid:88) l ∈ ∂k \ j e βE l ˜ ω ( k ) l i + c ˜ ω ( k ) l (1 + e − β ( E k − E l ) ) . (F1)In the population picture that one obtains for N → ∞ ,the analogous relation for the marginal precisions reads˜Ω c = i ( λ − i(cid:15) ) + i e − βE c (cid:88) k =1 e βE k ˜ ω k i + c ˜ ω k (1 + e − β ( E − E k ) )(F2)and these marginal precisions feed into the spectral den-sity (E6), whose population form is ρ ( λ ) = lim (cid:15) → π Re (cid:28) c ( { ˜ ω l , E l } , E ) (cid:29) ( { ˜ ω l ,E l } ,E ) , (F3)Now in the limit of small λ we expect the solution ofthe cavity equations to produce cavity precisions ˜ ω ( l ) k thatare purely imaginary, up to a real part of O ( (cid:15) ) [19]. Thisallows us to simplify the expression (F3) for the spectraldensity as follows: ρ ( λ ) = (F4)= lim (cid:15) → π Re (cid:42) O ( (cid:15) ) + i (cid:16) λ + e − βE (cid:80) ck =1 e βEk ˜ ω k i + c ˜ ω k (1+e − β ( E − Ek ) ) (cid:17) (cid:43) (F5)= (cid:42) δ (cid:32) λ + e − βE c (cid:88) k =1 e βE k ˜ ω k i + c ˜ ω k (1 + e − β ( E − E k ) ) (cid:33)(cid:43) . (F6)For small λ one sees that contributions to the δ -functioncome from large E , which makes sense as slow relax-ation rates should be associated with activation from thedeepest traps in the landscape (with E − ln c (cid:29)
1) thatare typically surrounded by higher neighbours. In thisregime we can drop the exponential from the denomina-tor in (F6) to get ρ ( λ ) ≈ (cid:42) δ (cid:32) λ + e − βE c (cid:88) k =1 e βE k ˜ ω k i + c ˜ ω k (cid:33)(cid:43) (F7)= (cid:90) d E (cid:90) d ξ c ρ E ( E ) ρ ξ ( ξ c ) δ ( λ + e − βE ξ c ) (F8)with ξ c = c (cid:88) k =1 e βE k ˜ ω k i + c ˜ ω k . (F9) With this, the original average in (F7) over E and thepairs ( ω k , E k ) has been translated into the average over E and the effective variable ξ c , which has an E -independentdistribution ρ ξ ( ξ c ). By evaluating the integral over theexponential energy distribution ρ E ( E ) we then arrive at ρ ( λ ) = κ c | λ | T − , κ c = T (cid:90) ∞ d ξ c ρ ξ ( ξ c ) ξ − Tc . (F10)Up to here we have accomplished our first aim, i.e.to derive the power law dependence (30) of the spectraldensity for small λ from the cavity equations. The pref-actor κ c still has to be found numerically from the λ → c analytically. Taking the limit of large c in (F9), the i in the denominator can be neglected and one has ξ ≈ c c (cid:88) k =1 e βE k (F11)Now the τ k = e βE k with E k drawn from ρ E ( E ) = e − E have a distribution with a power law tail ∼ τ − T − andhence a divergent mean for T <
1. The sum in (F11)is therefore dominated by its largest term [49] for which E k ≈ ln c , giving ξ ∼ c − e β ln c = c β − (F12)This shows that κ c ∼ ξ − T ∼ c − T ( β − = c T − (F13)which is the result (31) announced in the main text.The above scaling argument rests on simplifying (F9)using the approximation that c ˜ ω k /i (cid:29) ω k . In fact we had shown in (E5) that (for λ →
0, and ignoring the O ( (cid:15) ) real part) the cavity pre-cisions ˜ ω k = i Γ( E k ) are the escape rates, which scale as ∼ e − E k . So for the deepest traps with E k ≈ ln c , c ˜ ω k /i is just of order unity and such traps therefore make acontribution to (F9) that is somewhat smaller than weestimated. However, one can nonetheless show that thesetraps and even rarer, deeper ones still make a contribu-tion to ξ that scales as (F12). Appendix G: Single Defect Approximation