Glassy dynamics on networks: local spectra and return probabilities

TL;DR

This paper models glassy dynamics as a Markov process on a random network of energy minima, analyzing the spectral properties and return probabilities to understand slow relaxation and aging phenomena.

Contribution

It introduces a spectral analysis approach using the cavity method to study local density of states and return probabilities in glassy systems with random network connectivity.

Findings

Local density of states exhibits a power law behavior at small eigenvalues.

Return probabilities have power law tails influenced by trap lifetime distributions.

Results highlight the role of network connectivity and disorder in glassy dynamics.

Abstract

The slow relaxation and aging of glassy systems can be modelled as a Markov process on a simplified rough energy landscape: energy minima where the system tends to get trapped are taken as nodes of a random network, and the dynamics are governed by the transition rates among these. In this work we consider the case of purely activated dynamics, where the transition rates only depend on the depth of the departing trap. The random connectivity and the disorder in the trap depths make it impossible to solve the model analytically, so we base our analysis on the spectrum of eigenvalues $\lambda$ of the master operator. We compute the local density of states $\rho(\lambda|\tau)$ for traps with a fixed lifetime $\tau$ by means of the cavity method. This exhibits a power law behaviour $\rho(\lambda|\tau)\sim\tau|\lambda|^T$ in the regime of small relaxation rates $|\lambda|$, which we rationalize using a simple analytical approximation. In the time domain, we find that the probabilities of return to a starting node have a power law-tail that is determined by the distribution of excursion times $F(t)\sim t^{-(T+1)}$. We show that these results arise only by the combination of finite configuration space connectivity and glassy disorder, and interpret them in a simple physical picture dominated by jumps to deep neighbouring traps.

Equations

$${\partial }_t\mathbf{p}(t)=\mathbf{Mp}(t)$$

$$M_{ij}=c_{ij}{r}_{ij}\text{mbox}fori\ne j,M_{ii}=-\sum _{j(\ne i)}{M}_{ji}$$

$${\rho }_{\tau }(\tau )=T{\tau }^{-(T+1)}$$

$$\mathbf{p}(t)=\sum _{\alpha =0}^{N-1}{e}^{{\lambda }_{\alpha }t}({\mathbf{w}}_{\alpha },\mathbf{p}(0)){\mathbf{u}}_{\alpha }$$

$$\rho (\lambda )=\frac{1}{N}\sum _{\alpha =0}^{N-1}\delta (\lambda -{\lambda }_{\alpha })$$

$${\mathbf{M}}^{\text{s}}={\mathbf{P}}_{\text{eq}}^{-1/2}\mathbf{M}{\mathbf{P}}_{\text{eq}}^{1/2}$$

$$\rho (\lambda )=\underset{\epsilon \to 0}{\lim }\frac{1}{\pi N}\sum _{i=1}^N\mathrm{Im}{G}_{ii}({\lambda }_{\epsilon })$$

$$\mathbf{G}({\lambda }_{\epsilon })=({\lambda }_{\epsilon }\mathbf{I}-{\mathbf{M}}^{\text{s}}{)}^{-1}$$

$$\rho (\lambda )=\frac{1}{N}\sum _{i=1}^N\rho (\lambda |i)$$

$$\rho (\lambda |i)=\sum _{\alpha =0}^{N-1}\delta (\lambda -{\lambda }_{\alpha })v_{\alpha ,i}^2$$

$$\rho (\lambda |\tau )=\frac{1}{N_{\tau }}\sum _{i=1}^N\delta (\tau -{\tau }_i)\rho (\lambda |i)$$

$$\rho (\lambda )=\int \text{d}\tau \rho (\lambda |\tau ){\rho }_{\tau }(\tau )$$

$$1=\int \text{d}\lambda \rho (\lambda |\tau )$$

$$G_{ii}=i\int \mathrm{d}\mathbf{x}{x}_i^{2}P(\mathbf{x})=i\int \mathrm{d}{x}_i{x}_i^{2}P(x_i)$$

$$P(y_i)=e^{-\frac{i}{2}{\lambda }_{\epsilon }\frac{y_i^2}{r_i}}\prod _{k\in {\partial }_i}\int \mathrm{d}{y}_k{e}^{-\frac{i}{2}(y_i-y_k{)}^2}{P}^{(i)}(y_k)$$

$$P^{(i)}(y_k)=e^{-\frac{i}{2}{\lambda }_{\epsilon }\frac{y_k^2}{r_k}}\prod _{l\in {\partial }_k\setminus i}\int \mathrm{d}{y}_l{e}^{-\frac{i}{2}(y_k-y_l{)}^2}{P}^{(k)}(y_l)$$

$${\omega }_k^{(i)}=i{\lambda }_{\epsilon }{\tau }_{k}c+\sum _{l\in {\partial }_k\setminus i}\frac{i{\omega }_l^{(k)}}{i+{\omega }_l^{(k)}}$$

$${\omega }_i=i{\lambda }_{\epsilon }{\tau }_{i}c+\sum _{k\in {\partial }_i}\frac{i{\omega }_k^{(i)}}{i+{\omega }_k^{(i)}},$$

$$p(\omega )=\sum _k\frac{k{p}_k}{c}\int \mathrm{d}\tau {\rho }_{\tau }(\tau )\prod _{l=1}^{k-1}\mathrm{d}{\omega }_{l}p({\omega }_l)\delta (\omega -{\Omega }_{k-1})$$

$${\Omega }_a={\Omega }_a({\lambda }_{\epsilon },\{{\omega }_l\},\tau )=i{\lambda }_{\epsilon }\tau c+\sum _{l=1}^a\frac{i{\omega }_l}{i+{\omega }_l}$$

$$\rho(\lambda)=\lim_{\varepsilon\to 0}\frac{1}{\pi}\mathrm{Re}\,\Big{\langle}\frac{\tau c}{\Omega_{k}(\lambda_{\varepsilon},\{\omega_{l}\},\tau)}\Big{\rangle}_{\{\omega_{l}\},\tau,k}$$

$$\rho(\lambda|\tau)=\lim_{\varepsilon\to 0}\frac{1}{\pi}\mathrm{Re}\,\Big{\langle}\frac{\tau c}{\Omega_{k}(\lambda_{\varepsilon},\{\omega_{l}\},\tau)}\Big{\rangle}_{\{\omega_{l}\},k}$$

$$\rho^{\mathrm{2A}}(\lambda|\tau)=\lim_{\varepsilon\to 0}\frac{1}{\pi}\text{Re}\Big{\langle}\frac{\tau c}{\Omega_{k}(\{\Omega_{c-1}(\{\bar{\omega}\},\tau_{l})\},\tau)}\Big{\rangle}_{\{\tau_{l}\},k}$$

$$\rho^{\mathrm{2A}}(\lambda|\tau)\simeq c\tau\Big{\langle}\delta(k-\sum_{l}^{k}y_{l})\Big{\rangle}_{\{y_{l}\},k}$$

$${\rho }^{2A}(\lambda |\tau )\approx \tau C(c,T)|\lambda {|}^T\text{for}|\lambda |\ll 1/\tau$$

$$P_i(t)=\sum _{\alpha }{e}^{{\lambda }_{\alpha }t}{u}_{\alpha ,i}{w}_{\alpha ,i}=\sum _{\alpha }{e}^{{\lambda }_{\alpha }t}{v}_{\alpha ,i}^2$$

$$P_i(t)=\int \mathrm{d}\lambda \sum _{\alpha }\delta (\lambda -{\lambda }_{\alpha })v_{\alpha ,i}^2{e}^{\lambda t}=\int \mathrm{d}\lambda \rho (\lambda |i)e^{\lambda t}$$

$$P_{\tau }(t)=\frac{1}{N_{\tau }}\sum _{i=1}^N\delta (\tau -{\tau }_i)P_i(t)=\int \mathrm{d}\lambda \rho (\lambda |\tau )e^{\lambda t}$$

$$S_{\tau }(t)=e^{-t/\tau }$$

$$P_{\tau }(t)\approx \tau \tilde{C}(c,T)t^{-(T+1)}\text{for}t\gg \tau$$

Full text

1 Introduction

In pursuing a better understanding of non-equilibrium glassy systems, scientists have invested much effort into characterising their complex, multidimensional potential energy landscapes in configuration space. Key properties of these energy landscapes are the number of minima, the distribution of their depths, and of the heights of barriers between them. These have been explored using both computer simulations [1, 2, 3] and theoretical approaches [4, 5, 6]. The picture that has emerged is that the energy landscape of glasses is extremely complex, consisting of (exponentially many) minima, barriers and saddles of any order. The crystalline configurations that would be occupied in equilibrium at low temperature are hidden in this maze of valleys and walls and this keeps glassy systems out of equilibrium on typical observation timescales. When a glass is prepared, for example by quenching a viscous liquid to a low enough temperature, the system is expected to start descending rapidly towards a local minimum of the energy landscape [7, 8] (though how it does so is in itself not trivial [9]). As time proceeds, the system will then slowly explore progressively lower energy minima, leading to aging effects where physical properties depend on the time since preparation of the glass [10]. The dynamics in this regime can be thought of as consisting mostly of thermal fluctuations around an energy minimum, interspersed with rare large fluctuations that allow the system to cross a barrier and reach a new energy minimum. If one ignores the small thermal fluctuations around a given minimum, i.e. within a given “basin”, and focusses on the long time exploration of the various basins, then glassy dynamics can be modelled as a Markov process on a network, with each local energy minimum represented as a network node. A complete definition of such a model requires assumptions on the network topology, i.e. the connectivity in configuration space, and the transition rates between the nodes. The latter are expressed in terms of the energies of the various energy minima and the barriers between them. The distribution of the energy minima that enters here is expected on general grounds to have an exponential tail towards the deepest minima [11, 12].

The trap model [13, 14, 15] is one of the most successful descriptions of glassy dynamics that belong to this framework, where the energy minima are thought of as traps “hanging off” a threshold level where all barriers are located; in addition, the network of traps is assumed to be fully connected so that every trap can be accessed from any other. All moves between traps then require activation to the threshold level, which is convenient to use as the zero of the energy scale, and each jump takes the system to a randomly chosen new trap. The system thus effectively forgets with each jump what trap it was in before, making the dynamics a renewal process. The transition rates only depend on the departing energy depth because activation is always to the threshold level, and directly define the inverse lifetime of any trap. These simplifications allow the model to be solved analytically and give direct access to the evaluation of time dependent quantities. In particular, aging is described by two-time correlation functions that can be found explicitly and are given by the so-called arcsine law below the glass transition temperature [16]. A variety of disordered and more complex models of glasses exhibit an emergent trap-like phenomenology and aging behaviour, as demonstrated by numerical evidence as well as analytical arguments [17, 18]. However, the presence of dynamical correlations can make it hard to access the relevant timescales via simulations, and coarse-graining the evolution into larger effective basins may be required [19]. Importantly for us, the network of traps is generically not fully connected, and the original trap model then describes only the motion between the deepest effective basins at very long times. Note that a long time reduction of correlation functions to the arcsine law has been proven explicitly for the case of regular connectivity among the traps [20]. Various works have investigated in particular the case of lattices [14, 21, 22], though this is more plausible when the dynamics is interpreted as describing movement of a particle in real space rather than of a system in a high-dimensional configuration space. In the latter case, disorder in the connectivity among traps [23] inspired models of glassy dynamics on random networks. These are impossible to solve analytically, because of the disorder in trap depths and the random connectivity among nodes. Previous studies therefore had to rely on a heterogeneous mean field approximation [24, 25], which is uncontrolled. A different approach can be taken, however, by basing the analysis on the spectral properties of the master operator. This operator is the continuous-time analogue of a Markov transition matrix, and is key in determining the dynamics of the system. In particular the spectral density or density of states (DOS) $\rho(\lambda)$ gives the spectrum of relaxation rates of the system, and the localization properties of the eigenmodes carry information about the probability flow across the network. This is the approach that we followed in our previous work [26], which was dedicated to the analysis of trap models on sparse networks. In these models the zero energy threshold level remains present but jumps among traps are only allowed along network edges, i.e. local with respect to the network, so that the renewal property is lost.

In [26] we investigated the thermodynamic limit of an infinite network of traps by means of the cavity method. This approach exploits the local tree-like structure of networks with sparse connectivity and follows in this analogous applications to the spectral analysis of symmetric random matrices; see e.g. [27, 28, 29, 30, 31] or [32, 33] for a rigorous discussion. We used two relevant limits as benchmarks: the mean field (MF) limit where the average connectivity diverges in the thermodynamic limit of infinite system size, thus giving the original Bouchaud trap model, and the infinite temperature or random walk (RW) limit, where the energy landscape no longer plays a role. Our findings confirmed the idea that the very long time dynamics is well described by the original fully connected trap model and does not depend on the topology of the network of traps: the DOS always has a small-$|\lambda|$ tail – governing the long time relaxation – with the same power law behaviour as in mean field, $\rho(\lambda)\sim|\lambda|^{T-1}$. This can be rationalized within a simple high temperature approximation. In addition to this, our results indicated a decomposition of the dynamics into three different timescales: the long time (small $|\lambda|$), network independent regime with localized eigenmodes, the short time region where eigenmodes are delocalized and dominated by the network connectivity, and an intermediate regime where the DOS is as in mean field but the eigenmodes are delocalized nonetheless.

In this work we significantly extend our analysis of the trap model on sparse networks by looking at the local DOS, $\rho(\lambda|\tau)$, which gives the contribution to the (total) DOS from all traps with a fixed average lifetime $\tau$. The high $T$ approximation scheme again proves useful for deriving an analytical approximation for $\rho(\lambda|\tau)$ in the regime relevant for the long time dynamics: we find $\rho(\lambda|\tau)\sim\tau|\lambda|^{T}$ when $|\lambda|\ll 1/\tau$. These results are then translated into the time domain to give estimates for the return probability $P_{\tau}(t)$ and the distribution $F(t)$ of excursion times, i.e. the times required by the system to return to an initial trap, leading to $P_{\tau}(t)\sim t^{-(T+1)}$ for $t\gg\tau$ and $F(t)\sim t^{-(T+1)}$. Remarkably, it is the distribution of deep minima surrounding the initial trap that determines these power laws, and they arise as a combined effect of limited connectivity and trap depth disorder: if only one of these features is present, the local DOS becomes concentrated around $-1/\tau$ implying an exponential decay of the return probability.

The paper is organised as follows: after defining the model in section 2, we present our cavity theory in section 3, describe the approximation scheme and sketch the result for the power law tail of the local DOS; the full derivation is left to appendix B. In section 4 and 5 we focus on the behaviour in the time domain by analysing, respectively, the return probability and the excursion time distribution. Finally, we summarise and discuss our results in section 6.

2 Bouchaud trap model on networks

We follow the set-up of our previous work [26]: the problem is defined by a continuous-time Markov process on a sparse network (or graph), whose nodes (or traps) represent the minima of the energy landscape where system gets trapped. These have a positive energy that represents the depth of the minimum with respect to the level zero of the energy landscape, and determines the expected lifetime of that state. Trap depths ($E>0$) are quenched random variables following the exponential distribution $\rho_{E}(E)=\theta(E)\text{exp}(-E)$. The master equation defining the Markov process is

[TABLE]

where $\mathbf{p}(t)=(p_{1}(t),\ldots,p_{N}(t))$ is the probability distribution describing the position of the system on the network. The elements of the master operator $\mathbf{M}$ are:

[TABLE]

where $r_{ij}$ are the Bouchaud transition rates $r_{ij}=e^{-\beta E_{j}}/c\equiv r_{j}$, $c$ is the average connectivity of the network, $\beta$ is the inverse temperature and $c_{ij}$ is $1$ if $i$ and $j$ are connected, and [math] otherwise. It is useful to define the quantity $\tau_{j}=(cr_{j})^{-1}=\text{exp}(\beta E_{j})$, which sets the scale of the expected waiting time ($c/k_{j})\tau_{j}$ to leave a node $j$; here $k_{j}$ is the degree of the node. The distribution of energies $E$ implies a distribution for $\tau$ given by

[TABLE]

Note that $\langle\tau\rangle$ diverges for $T\leq 1$, signalling a low $T$ regime where the dynamics gets glassy. In this work we will mostly focus on the case where the network connectivity is that of a random regular graph (RRG), i.e. where every node is connected to $c$ random nodes, with $c\geq 3$ so that the fraction of nodes outside the giant connected component of the graph vanishes in the large $N$ limit [34]. This case is the simplest and yet it exhibits the same key features as more complex network topologies. This is confirmed by the results shown at the end of section 3.1, where the cases of Erdös-Rényi and scale-free connectivities are discussed.

The assumption of a configuration space characterised by a sparse and random connectivity makes the master equation (1) impossible to solve analytically. We therefore take another route and focus on the spectral properties of $\mathbf{M}$, whose $\alpha^{\text{th}}$ eigenvalue, left and right eigenvectors we write respectively as $\lambda_{\alpha},\mathbf{w}_{\alpha}$ and $\mathbf{u}_{\alpha}$. A formal solution to (1) is then given by

[TABLE]

where $(\cdot,\cdot)$ denotes the scalar product between vectors. If the network is connected there is only one vanishing eigenvalue $\lambda_{0}=0$, and the corresponding right eigenvector represents the equilibrium distribution of the system: $\mathbf{p}_{\text{eq}}=\lim_{t\to\infty}\mathbf{p}(t)=\mathbf{u}_{0}$. All other eigenvalues have a negative real part and the contribution of the associated eigenvectors to $\mathbf{p}(t)$ is exponentially suppressed over time. We will often refer to the eigenvectors of $\mathbf{M}$ as the eigenmodes (or simply the modes) of the dynamics, so e.g. we could say that the long time behaviour of the system is governed by the slow modes, thus referring to the eigenvectors in the small $|\lambda|$ regime. The importance of the spectrum of eigenvalues for the dynamics is evident as it provides the distribution of relaxation rates of the system. For this reason, a central quantity for our analysis is the (total) density of states (DOS), defined as

[TABLE]

One could equivalently consider the spectrum of the relaxation rates $r_{\alpha}=-\lambda_{\alpha}$, which would just flip the sign of lambda. We stick to the convention in (5) for consistency with our earlier work [26]. There we discussed in some detail the features of the DOS of the trap model defined on sparse networks. In all cases that we considered the DOS showed a $|\lambda|\to 0$ power-law tail with the same exponent as found in the case of mean field connectivity, and eigenvectors exhibiting a localization transition, from delocalized fast modes to localized slow modes. We measured the degree of localization in terms of the inverse participation ratio (IPR), using a formula proposed by Bollé et al [29] to detect localization transitions in symmetric random matrices. In order to investigate the thermodynamic limit we relied on the cavity method, and used a population dynamics algorithm to solve the associated cavity equations numerically. This method links the master operator to the inverse covariance matrix of a complex Gaussian distribution and therefore requires a symmetric matrix as input, which in our case is obtained from the similarity transformation

[TABLE]

Here $\mathbf{P}_{\text{eq}}$ is a diagonal matrix with non-zero elements given by the equilibrium distribution: $(\mathbf{P}_{\text{eq}})_{ii}=(\mathbf{p}_{\text{eq}})_{i}$. This transformation preserves the eigenvalue spectrum, which is real as $\mathbf{M}^{\text{s}}$ is real and symmetric. Also, it does not affect the diagonal elements of the master operator, i.e. $(\mathbf{M}^{\text{s}})_{ii}=(\mathbf{M})_{ii}$, a fact that will turn out to be crucial for us. The eigenvectors $\mathbf{v}_{\alpha}$ of $\mathbf{M}^{\text{s}}$ are given by $\mathbf{v}_{\alpha}=\mathbf{P}_{\text{eq}}^{-1/2}\mathbf{u}_{\alpha}=\mathbf{P}_{\text{eq}}^{1/2}\mathbf{w}_{\alpha}$. These retain the same localization properties as the eigenvectors of the original system, except for finite size effects mostly appearing close to the ground state (see [26], appendix E). The symmetry of $\mathbf{M}^{\text{s}}$ is a consequence of the detailed balance condition that holds between the transition rates and the equilibrium Boltzmann distribution $\mathbf{p}_{\text{eq}}$ [35], see also [36] for a discussion in the context of Fokker-Plank evolution. Using a standard identity from random matrix theory [37], the DOS of $\mathbf{M}^{\text{s}}$ can be expressed as

[TABLE]

in terms of the resolvent

[TABLE]

Here $\mathbf{I}$ indicates the $N\times N$ identity matrix and we have used the abbreviation $\lambda_{\varepsilon}=\lambda-\operatorname*{\mathrm{i}}\varepsilon$ with $\operatorname*{\mathrm{i}}$ the imaginary unit and $\varepsilon$ small and positive. In going from (5) to (7) one replaces the delta functions in (5) with Lorentzians of width $\varepsilon$; thus $\varepsilon$ sets the numerical resolution that we have on the $\lambda$-axis when we come to evaluate quantities of interest, using in our case specifically the population dynamics algorithm outlined below.

In this work we use the cavity method to study the DOS in more detail. In particular we decompose it into a set of local DOSs, one for each node; these local DOSs are defined explicitly below. The analysis allows to probe the important effects of heterogeneity in the network, as generated by the landscape of trap depths $E_{i}$. It will also enable us to obtain insights into the dynamics directly in the time domain, as e.g. time-dependent probabilities of return to a certain trap can be easily computed using results for the local DOSs.

3 Local DOS

The $i^{\text{th}}$ term appearing in the sum on the right hand side of equation (7) is the contribution to the total DOS given by a single node. To make this explicit, we can write

[TABLE]

with $\rho(\lambda|i)=\lim_{\varepsilon\to 0}\mathrm{Im}\,G_{ii}(\lambda_{\varepsilon})/\pi$. This quantity is referred to as the (single node) local DOS. By translating the eigendecomposition of $\mathbf{M}^{\rm s}$ into one for the resolvent matrix, $\mathbf{G}(\lambda_{\varepsilon})=\sum_{\alpha}(\lambda_{\varepsilon}-\lambda_{\alpha})^{-1}\mathbf{v}_{\alpha}\mathbf{v}_{\alpha}^{\rm T}$, one sees that the local DOS can be written more explicitly as

[TABLE]

The normalization of eigenmodes implies that summing over all network nodes $i$ and dividing by $N$ gives the total DOS defined by (5), as written in (9). More generally, it is possible to decompose the total DOS according to the contribution from all traps with a given property, e.g. a fixed local timescale $\tau_{i}=\tau$. Following this idea, we define $\rho(\lambda|\tau)$ as

[TABLE]

where $N_{\tau}=\sum_{i=1}^{N}\delta(\tau-\tau_{i})$. The following relations then hold:

[TABLE]

[TABLE]

Here $\rho_{\tau}=N_{\tau}/N$ is the probability density function of $\tau$ for a given realization of the system with size $N$, which for $N\to\infty$ is self-averaging and given by the expression (3). Of course the same construction can be used to define a local DOS conditioned on generic local disorder variables; the node degree $k_{i}$ would be an obvious choice (see e.g. [38]), though we do not pursue this here.

In the next section we show how to use the cavity method to evaluate the local DOS $\rho(\lambda|\tau)$, and we derive an analytical approximation valid for the small $|\lambda|$ tail based on the high $T$ approximation scheme that we introduced in [26].

3.1 Cavity method

Here we only present a brief summary of the cavity construction and refer to our previous paper [26] for a detailed derivation of the central result, i.e. the self-consistent equation for the distribution of cavity precisions given below.

First of all, one observes that the diagonal elements of the resolvent can be expressed as

[TABLE]

where $P(\mathbf{x})\propto\text{exp}(-\operatorname*{\mathrm{i}}\mathbf{x}^{\rm T}\mathbf{G}^{-1}\mathbf{x}/2)$ is a complex Gaussian measure with covariance matrix given by $\mathbf{G}$, and $P(x_{i})$ the associated marginal distribution at node $i$. Exploiting the sparse structure of $\mathbf{G}^{-1}$, one can express the marginal distribution of $i$ in terms of the cavity distributions of its neighbouring nodes $k\in\partial_{i}$:

[TABLE]

Here $\partial_{i}$ indicates the neighbourhood of $i$, and we have used the change of variable $y_{i}=x_{i}\sqrt{r_{i}}$, which has the desired effect of confining the disorder from the transition rates $r_{i}$ to the diagonal terms. The last equation is based on the key assumption that the joint distribution of the nodes belonging to $\partial_{i}$ factorises when the central node $i$ is removed from the graph (see the illustration in Fig. 1).

This is strictly true only when the network formed by the traps and allowed transitions between them is a tree, but it also provides a valid approximation whenever the topology of the network is at least locally tree-like, as is the case of sparse networks (e.g. with random regular or scale-free connectivity) in the large $N$ limit [39]. Using the same line of reasoning as for (15) one can write for the cavity distributions the recursive relation

[TABLE]

Equation (16) is self-consistently solved by Gaussian distributions of the form $P^{(i)}(y_{k})\propto\text{exp}(-\omega^{(i)}_{k}y_{k}^{2}/2)$. This ansatz transforms Eq. (16) into an equivalent set of equations for the cavity precisions:

[TABLE]

The Gaussian nature of the cavity marginals entails that single site marginals $P(y_{i})$ are also Gaussian, with Eq. (15) implying that single site precisions $\omega_{i}$ are of the form

[TABLE]

The system of equation for the cavity precisions (17) can be solved recursively for a finite realization of the system. The marginal precisions $\omega_{i}$ are then obtained from (18) and they give the diagonal elements of the resolvent via $G_{ii}=\operatorname*{\mathrm{i}}\tau_{i}c/\omega_{i}$.

We are concerned with the thermodynamic limit of a large network. Here we can exploit that for $N\to\infty$, where due to the locally tree-like assumption loops in the network become long, the different terms in the sum in (17) become uncorrelated samples from the distribution $p(\omega)$ of cavity precisions. Requiring that the left hand side, too, is a sample from the distribution of cavity precisions, one obtains a self-consistency equation for $p(\omega)$. (Note that here and in the following we omit the superscript indicating the cavity graph in order to keep the notation simple.) For a general degree distribution $p_{k}$ this self-consistency equation reads

[TABLE]

with the abbreviation

[TABLE]

A numerical solution of (19) can be found using a population dynamics algorithm; see [40] for a detailed explanation of this method. The core idea is to take an initialised population $\mathcal{P}$ of cavity precisions $\omega$ and then update each of them using the value of $\Omega_{k-1}$ given by a sample $\tau$ and $k-1$ random elements of $\mathcal{P}$ as inputs, with the node degree $k$ sampled appropriately from the degree-weighted distribution $kp_{k}/c$. This process is then repeated until the statistics of the distribution $p(\omega)$ converge. At this point the total DOS $\rho(\lambda)$ can be evaluated as

[TABLE]

where the angle brackets $\langle\ldots\rangle_{\{\omega_{l}\},\tau,k}$ indicate averaging over the degree distribution $p_{k}$, the lifetime distribution $\rho_{\tau}(\tau)$ and the $k$ cavity precision distributions111The full expression reads $\rho(\lambda)=\lim_{\varepsilon\to 0}\frac{1}{\pi}\mathrm{Re}\sum_{k}p_{k}\int\mathrm{d}\tau\rho_{\tau}(\tau)\,\prod_{l=1}^{k}\mathrm{d}\omega_{l}\,p(\omega_{l})\,\tau c/\Omega_{k}(\lambda_{\varepsilon},\{\omega_{l}\},\tau)$ $\prod_{l=1}^{k}p(\omega_{l})$. Comparing (21) with (12) one reads off directly that the local DOS is given by an almost identical expression

[TABLE]

that differs only in the fact that $\tau$ is fixed rather than averaged over. The remaining average is, in practice, evaluated by sampling values of $k$ from the degree distribution, and sets $\{\omega_{l}\}$ of $k$ cavity precisions from the population $\mathcal{P}$ converged to equilibrium.

Figure 2-left shows the local DOS obtained by the above method on a log-log scale, for the random regular graph ensemble with $c=5$ and $T=0.8$, and for given $\tau=2$. The factor $-\lambda$ on the y-axis accounts for the transformation $\lambda\to\ln(-\lambda)$ on the x-axis, so this plot can be read as the distribution of the logarithmic relaxation rates, $\ln(-\lambda)$, with the correct normalization. The blue lines show the results obtained with the population dynamics algorithm described above, using two different values of $\varepsilon$. We note that a smaller value of $\varepsilon$ gives a better resolution on the negative lambda axis, as it should, and that straight $\varepsilon$-dependent tails appear in the small $|\lambda|$ region (below the corresponding values of $\varepsilon$), as well as for $|\lambda|>1$. These have no physical meaning: population dynamics sampling runs of finite length produce only a limited number of samples in this region, if any, and the resulting shape of the (local) DOS is strongly affected by the Lorentzians of width $\varepsilon$ that the method effectively uses to smooth the spectra. The green dashed lines were obtained from direct diagonalizations (labelled “numerics” in the plots) of instances of $\mathbf{M}^{\textbf{s}}$ with network size $N=500,1000,2000$ (light to dark green). The agreement with the population dynamics results is excellent in the regions where finite size effects are absent. The latter do show up in the small $|\lambda|$ regime and are again a consequence of limited sample sizes, with finite matrices only rarely having eigenvalues in this region222The intuition here is that the eigenvalues determine the relaxation times of the system, and exploring a smaller network must take less time, on average. More precisely, the largest trapping time $\tau_{\text{max}}$ in a finite system of size $N$ scales as $\tau_{\text{max}}\sim N^{\beta}$ [26]..

In figure 2-right we explore the dependence of the local DOS on $\tau$ as predicted by the population dynamics algorithm, for the same setting of random regular networks with $c=5$ at $T=0.8$. The grey solid line in the background is the total DOS given by (21). This plot provides two important insights about the local DOS: (i) most of its mass is peaked around a value of $-\lambda$ that scales as $1/\tau$, with this peak getting narrower as $\tau$ increases; (ii) in the small $|\lambda|$ regime well below the peak, the local DOS has a power law dependence on $|\lambda|$, with a $\tau$-independent exponent.

At this stage it is useful to compare to the simpler MF ($c\to\infty$) and RW ($T\to\infty$) limits discussed in the introduction. We discuss their local DOS in Appendix A and find in both cases a delta peak at $|\lambda^{*}|\propto 1/\tau$ for large $\tau$. (The proportionality constant is unity for the MF case, where the delta peak is the only contribution to the local DOS; in the RW case there is an additional piece to the spectrum for $|\lambda|$ of order unity.) The delta peak is the analogue of the smooth peaks visible in figure 2-right. More importantly, the local DOS turns out to be zero for $|\lambda|$ below the peak. The power law behaviour in this regime that we see in figure 2-right therefore has no analogue in either the MF or RW limits: it arises only as a combined effect of limited connectivity and trap depth disorder. Remarkably, the power law behaviour for small $|\lambda|$ is robust to changes in the network connectivity, as can be seen by comparing results for the Erdös-Rényi and scale-free network ensembles in figure 3 to those for random regular networks in 2-right. One observes that while the shape of $\rho(\lambda|\tau)$ around $\lambda=-1/\tau$ does depend on the type of network, the exponent of the power law tail for low $|\lambda|$ does not.

In the next section we present an approximation scheme that can explain the power law behaviour observed in the local DOS for $|\lambda|\ll 1/\tau$, and clarifies that this result applies to any sparse network specified by some degree distribution $p_{k}$ with finite mean. In fact, the argument that we use is insensitive to correlations among node degrees, and should therefore remain valid for networks generated e.g. by preferential attachment [39].

3.2 Approximation scheme

While the equations resulting from the cavity method do not permit closed-form solutions, we can use them to construct an approximation scheme that provides useful insights. This is done in the spirit of the single defect approximation [41, 42]: the main idea is to take into account the disorder of a certain region of interest only, e.g. a single node or a given neighbourhood in the case of a network, and assume the rest of the system to be homogeneous. In our model, a first order (or “first shell”) approximation of this kind corresponds to assuming $T=\infty$ and a $c$-regular connectivity on the cavity network, and taking only the lifetime $\tau$ and connectivity $k$ of the central node into account. With these assumptions the cavity network becomes disorder-free in the thermodynamic limit, which implies that the distribution $\rho(\omega)$ of cavity precisions becomes a delta function centred on the value of $\bar{\omega}$ that solves $\bar{\omega}=\Omega_{c-1}(\lambda_{\varepsilon},\{\bar{\omega}\},1)$. We refer to $\bar{\omega}$ as the infinite-$T$ (and $c$-regular) solution. The first order approximation is then implemented as one cavity step – involving the local $\tau$ and $k$ values – performed starting from the infinite-$T$ solution. Following this idea, the second order (or second shell) approximation consists of two cavity steps from the infinite-$T$ solution and so on. In what follows we will focus on the random regular graph ensemble for simplicity, and refer to the appendices for the demonstration of the wider applicability of the approximation.

The first order approximation was previously found to give an accurate description of the total DOS in the small $|\lambda|$ regime, where $\rho(\lambda)\sim(-\lambda)^{T-1}$ as in mean field [15, 26]. For the local DOS, on the other hand, the first order approximation is not able to provide a match to the power law behaviour shown in figure 2. In fact, this approximation coincides with the random walk limit discussed in appendix A, which yields a local DOS characterised by a single delta peak in the small $|\lambda|$ regime. This is due to the complete lack of randomness in the approximation: the attributes of the central node are fully specified (by $\tau$ and the fixed degree $k=c$) as are those of the cavity network by the infinite-$T$ solution. The second order approximation is then required if we want to include some of the original heterogeneity of the system: in the thermodynamic limit, different nodes with the same value of $\tau$ have different neighbourhoods, and averaging across these neighbourhoods gives the approximated local DOS. This can be expressed mathematically as

[TABLE]

where the superscript 2A stands for “second order approximation”. In appendix B we show that for $|\lambda|\ll 1/\tau$ Eq. (23) can be written as

[TABLE]

with $y_{l}=1/(\lambda c[\tau_{l}+(c-2)^{-1}]+c-1)$. The distribution $\rho_{y}(y_{l})$ of $y_{l}$ is strongly peaked within a region of order $|\lambda|^{T}$ around $1/(c-1)$, and it drops by a factor $|\lambda|^{T}$ if $y_{l}$ is away from this value. The average on the right hand side of (24) is then dominated by cases where all the $y_{l}$ are close to $1/(c-1)$ except for one, say $y_{1}$. The latter then has to equal $k-(k-1)/(c-1)$, which happens when the corresponding $\tau_{1}$ is of order $1/|\lambda|$. Substituting the expression for $\rho_{y}(y_{1}=k-(k-1)/(c-1))$ and multiplying by a factor $k$ (as any of the $\tau_{l}$ could be the large one) this leads to the following result:

[TABLE]

where the prefactor $C(c,T)$ depends on the degree distribution of the network and vanishes in the limit $c\to\infty$, as required to recover the mean field case where the local DOS vanishes for $|\lambda|\ll 1/\tau$. For random regular networks we find explicitly $C(c,T)=Tc(c-1)^{T-1}(c-2)^{-(T+1)}$. Equation (25) then also implies $\rho^{\mathrm{2A}}(\lambda|\tau)\to 0$ when $T\to\infty$, at least for $|\lambda|<(c-2)/(c-1)$, and so – like the first order approximation – it is consistent with the random walk limit. The predicted small-$|\lambda|$ power law of the local DOS, $\rho(\lambda|\tau)\sim|\lambda|^{T}$, matches the full population dynamics results well (see the black dashed line named “$2^{\text{nd}}$ shell approximation” in figure 2-left, and the blue dashed lines in figures 2-right and 3). Interestingly, the argument that leads to Eq. (25) implies that the observed power law arises because of deep minima surrounding the node of interest: a significant contribution to the local DOS at any given $\lambda$ comes from those traps that have at least one neighbouring minimum with expected lifetime $\gtrsim 1/|\lambda|$. As we will see, this feature of the local DOS determines the observed long time behaviour of several quantities of interest, and in particular it has important implications for the return probability discussed in the next section.

4 Return probability

We now turn our attention to the time domain. We start by using the results for the local DOS to compute the average probability of return to an initial trap, as a function of its lifetime $\tau$.

Let us call $P_{i}(t)$ the return probability to an arbitrary initial trap $i$. Physically, this gives the probability of being in trap $i$ at time $t$, given that we have started out in the same trap. Mathematically it is given by equation (4) with the initial condition $p_{j}(0)=\delta_{ij}$, yielding

[TABLE]

In the second equality we have used the relations $u_{\alpha,i}=p_{{\rm eq},i}^{1/2}\,v_{\alpha,i}$ and $w_{\alpha,i}=p_{{\rm eq},i}^{-1/2}v_{\alpha,i}$ to transform to an expression in terms of the eigenvectors of the symmetrized master operator; note that the symmetrization factors cancel. Decomposing the sum according to the eigenvalues $\lambda$ gives

[TABLE]

where $\rho(\lambda|i)$ is the single node local DOS defined in (10). We now define an average return probability over all traps with lifetime $\tau_{i}=\tau$:

[TABLE]

As the second equality shows, the return probability $P_{\tau}(t)$ can be obtained directly from the local DOS $\rho(\lambda|\tau)$, which in turn we can predict using population dynamics as explained above. We can also use (28) in reverse to deduce that the average of $\lambda$ over the local DOS $\rho(\lambda|\tau)$ is given by $-1/\tau$, in accordance with our findings about the scaling of the peaks observed in figure 2. This general result for the average of $\lambda$ can be seen by noting that $P_{i}(t)=(\exp\mathbf{M}t)_{ii}=1+M_{ii}t+O(t^{2})$. From the definition (2) of the master operator this equals $1-t/\tau_{i}+O(t^{2})$; comparing with the expansion of the r.h.s. of (28) to linear order in $t$ then gives the result.

Figure 4 shows the return probability $P_{\tau}(t)$ for the random regular graph ensemble with $c=5$, $T=0.8$ and $\tau=2$ obtained using different approaches: the blue line is computed from the population dynamics result discussed in the previous section, the green dashed lines are obtained using data from direct diagonalizations of samples of $\mathbf{M}^{\mathrm{s}}$-matrices of different sizes (light to dark green), and the red line shows the result of stochastic simulations with the Gillespie algorithm (see appendix C). From short times up to $t\sim 10^{1}$ the agreement between these three approaches is clearly very good. At larger $t$, finite size effects induce upward curvature in the direct diagonalization curves as the estimate of $P_{\tau}(t)$ approaches the Boltzmann equilibrium distribution in a finite system. The simulations, which are performed on effectively infinite networks and so do not suffer from analogous errors, start to exhibit statistical sampling fluctuations in the same range of times. (With the $10^{5}$ runs of the simulated dynamics that we use, a reasonable estimate of the return probability can only be obtained for $P_{\tau}(t)\gtrsim 10^{-4}$; we therefore do not show data beyond this point in Fig. 4.)

Looking now in more detail at the results in Figure 4, one observes that for times $t$ up to the order of the trap lifetime $\tau$, $P_{\tau}(t)$ is dominated by the “staying probability” of having never left the initial trap, which we will denote by $S_{\tau}(t)$. This quantity coincides with $P_{\tau}(t)$ for the fully-connected graph, which is the original Bouchaud model: once the initial trap has been left, the probability of coming back is $O(1/N)$ and vanishes in the thermodynamic limit. In our continuous time framework the staying probability is simply given by (see grey solid line in figure 4)

[TABLE]

The power law tail of the return probability observed beyond this, in the long time regime ($t>\tau$), is therefore clearly a network effect in the dynamics: the finite connectivity allows for returns to the initial trap even in the infinite system size limit, and as we will see this generates a power law tail provided that $T$ is finite so that the disorder in the trap depths matters. The power law predicted by the cavity theory is consistent with our simple estimate (25) of the local DOS for $|\lambda|\ll 1/\tau$, which when inserted into the integral in (28) yields

[TABLE]

The prefactor here is $\tilde{C}(c,T)=C(c,T)\,\Gamma(T+1)$, with $\Gamma(z)$ the Gamma function. This approximation agrees with the cavity theory exactly regarding the power law exponent; even the prefactor is quantitatively close (compare the blue solid line and the black dashed line in figure 4).

In the previous section we explained that the derivation of (25) implies that the small $|\lambda|$ power law tail of the local DOS originates from deep minima (with lifetimes $\gtrsim 1/|\lambda|$) surrounding the initial trap. The implication for the return probability is the following: the most likely manner in which the system can return to the initial state at $t\gg\tau$ is for it to become trapped in a neighbouring node with lifetime $\sim t$. Other possibilities of course exist, e.g. the system could come back from a trap in the second neighbour shell, but these only make a sub-dominant contribution to $P_{\tau}(t)$. This can be seen indirectly from the fact that our approximation almost overlaps with the numerically exact population dynamics curve in figure 4.

Summarizing, $P_{\tau}(t)$ exhibits an initial exponential decay typical of the mean field limit for $t\lesssim\tau$, followed by a power law behaviour with a $T$-dependent exponent for $t\gtrsim\tau$; the latter arises from deep minima surrounding the departing node. Note that the crossover point from the exponential decay to the power law regime occurs at a value of $P_{\tau}(t)$ that decreases as $\tau$ increases, which is directly related to the fact that the power law tail does not just depend on the scaled time $t/\tau$ (see (30)). Both features can be seen in figure 5, where we compare evaluations of $P_{\tau}(t)$ for different values of $\tau$: in the left plot the x-axis is scaled by $\tau$, which delivers a collapse of the exponential decay at short times, while in the right plot we have similarly scaled the y-axis to show that the prefactor of the tail is indeed proportional to $\tau$ as the approximation (30) predicts. Note that the time constant of the initial exponential decay is somewhat larger than the mean field value $\tau$. In fact for $\tau\gg 1$ this can be approximated by $\tau(c-2)/(c-1)$, as we will justify in the next section.

We note finally that the conditioning on the trap lifetime in the return probability $P_{\tau}(t)$ is essential in order to isolate the effects of finite, non-mean field network connectivity among traps. If one instead considers the probability of return to a randomly chosen initial trap $P(t)=\int\mathrm{d}\tau\,P_{\tau}(t)\rho_{\tau}(\tau)$, one finds that this is dominated by the initial exponential decay of $P_{\tau}(t)$, which is exactly the mean field result. The long time scaling $P(t)\sim t^{-T}$ that one deduces has a mean field form even for finite connectivity $c$. The decay exponent is consistent with the small $\lambda$-scaling of the total DOS [15, 26], which again is independent of $c$. This analysis tells us that the trap model on networks exhibits a long time mean field dynamics on average only, with network effects coming to the fore when considering more detailed phenomena like returns to specific initial traps as considered here.

5 Excursion times

In the previous section we discussed the dynamical properties of the system in terms of the return probability to some initial trap. We saw that the local disorder, i.e. the depth of the departing trap, determines the shape of $P_{\tau}(t)$ at short times. In contrast, the long time behaviour is always a power law in $t$, with the dependence on $\tau$ entering only via the prefactor. We suggested that this is because when $t$ is much larger than the trap lifetime $\tau$, the probability of finding the system in the original trap is dominated by the lifetimes of the neighbouring minima. In this section we consolidate this idea by excluding the time spent in the initial trap from the analysis and looking at the distribution of excursion times, i.e. the time spent by the system away from the initial trap between two visits there. To this end we decompose the return probability as

[TABLE]

where $P_{\tau}^{(n)}(t)$ indicates the probability of finding the system in the original trap at time $t$, assuming that it has returned there $n$ times in total. Note that $P_{\tau}^{(0)}(t)=S_{\tau}(t)$, while $P_{\tau}^{(1)}(t)$ can be written as

[TABLE]

Here the integrand is the probability that the system leaves the initial trap at time $t_{1}$ (we write this as $L_{\tau}(t_{1})\,\mathrm{d}t_{1}$), stays away until time $t_{2}$ ($F(t_{2}-t_{1})\,\mathrm{d}t_{2}$), then returns to the origin and remains there until $t$ ($S_{\tau}(t-t_{2})$). Reading the expression in this way implies that $F(t)$ is the distribution of excursion times, the quantity that we want to evaluate. Equation (32) is the convolution between $L_{\tau}$, $F$ and $S_{\tau}$, i.e. $P_{\tau}^{(1)}=L_{\tau}\ast F\ast S_{\tau}$, which can be generalised to $n$ returns as

[TABLE]

with $(L_{\tau}\ast F)^{(n)}$ the $n$-fold convolution between $L_{\tau}$ and $F$. Substituting (33) into (31) and taking the Laplace transform leads to

[TABLE]

which can be written more explicitly using that, from (29),

[TABLE]

Also one has $L_{\tau}(t)=-S^{\prime}_{\tau}(t)$ or in Laplace space $\hat{L}_{\tau}(s)=1-s\hat{S}_{\tau}(s)=1/(1+s\tau)$. Substituting into (34) yields

[TABLE]

This equation can be inverted to obtain an expression for the excursion time distribution in terms of the return probability,

[TABLE]

We consider first the limit $s\to 0$. The value $\hat{F}(0)=\int_{0}^{\infty}\mathrm{d}t\,F(t)$ gives the probability that an excursion lasts any finite amount of time, i.e. the probability that the system will sooner or later go back to the initial node rather than escape to infinity. On the infinite $c$-regular tree, i.e. for a random $c$-regular graph in the limit $N\to\infty$, this fixes333A simple argument to see this is the following. The probability of ever returning – call this $P_{0}$ – cannot depend on the lifetimes of the traps in the configuration space: if the system escapes to infinity it does not matter how long this will take, and so $P_{0}$ is independent of $T$ as long as $T>0$. Moreover we have $P_{0}=P_{01}$, where $P_{nm}$ is the probability to ever land on a node in the $n^{\mathrm{th}}$ neighbour shell of the initial trap, starting from the $m^{\mathrm{th}}$ shell (here the $n=0$ “shell” is the initial trap). It is immediate to see that $P_{01}=1/c+P_{02}(c-1)/c$ and also $P_{02}=P_{01}P_{12}$, which by symmetry becomes $P_{02}=P_{01}^{2}$. The resulting second order equation gives the physical solution $P_{01}=1/(c-1)$, which correctly becomes $P_{01}=1$ for a chain ($c=2$) and $P_{01}=0$ in the MF limit ($c\to\infty$). In general, the probability to reach a node that is $l$ steps away decreases exponentially, $P_{0l}=(c-1)^{-l}$, which suggests a possible explanation for why the $2^{\mathrm{nd}}$ shell approximation works so well. $\hat{F}(0)=1/(c-1)$, and so we obtain $\hat{P}_{\tau}(0)=\tau(c-2)/(c-1)$ from (37), while $\hat{S}_{\tau}(0)=\tau$. In the MF limit we have $\hat{P}_{\tau}(s)=\hat{S}_{\tau}(s)$, implying $\hat{F}_{\tau}(s)=0$ via (37), which is consistent with the $c\to\infty$ limit of $\hat{F}(0)=1/(c-1)$.

Next we analyse what (37) says about the long time behaviour of $F(t)$, by considering $\hat{F}(s)$ for small $s$. From equation (30) one obtains the following approximation for the leading singular small $s$-behaviour of the return probability in Laplace space444 One has generally $\hat{P}_{\tau}(0)-\hat{P}_{\tau}(s)=\int_{0}^{\infty}\mathrm{d}t\,P_{\tau}(t)(1-e^{-st})$. The integrand becomes dominated by large $t$ for small $s$; substituting the tail estimate (30) and integrating by parts then gives (38).

[TABLE]

Substituting this expression for $\hat{P}_{\tau}(s)$ into (37) and expanding again for small $s$ one sees that $\hat{F}(s)$ contains the same singular term:

[TABLE]

For the long time behaviour of $F(t)$ this implies the same power law that we found in the return probability:

[TABLE]

where we have used $\hat{P}_{\tau}(0)=\tau(c-2)/(c-1)$. Note that the predicted behaviour of the excursion time distribution only depends on the average connectivity and temperature, while the departing lifetime $\tau$ disappears from (40) as it did in our earlier result $\hat{F}(0)=1/(c-1)$. This is as expected: in the Bouchaud trap model, the escape time $\tau$ only contains information on the local disorder (trap depth) at the initial minimum, so once this trap has been left, the behaviour during the following excursion is independent of $\tau$.

So far our analytical reasoning was based on an approximation for the local DOS, which we converted into a return probability $P_{\tau}(t)$ and finally into the excursion distribution $F(t)$. Qualitatively, we can also alternatively argue directly from $F(t)$, by constructing a simple lower bound. The probability of an excursion taking longer than $t$, $\int_{t}^{\infty}\mathrm{d}t^{\prime}F(t^{\prime})$, is at least as large as the probability of not having left the first trap encountered during the excursion. As the depth of this trap is random, the latter probability is $\int\mathrm{d}\tau^{\prime}\,\rho_{\tau}(\tau^{\prime})e^{-t/\tau^{\prime}}$. This lower bound is just the mean field return probability $P(t)\sim t^{-T}$ discussed at the end of the previous section. Taking a derivative w.r.t. $t$ gives the estimate that $F(t)$ should decay as $t^{-(T+1)}$, exactly as we had found in (40). This then implies the analogous power law (30) in the return probability $P_{\tau}(t)$, and in turn via (28) the small $|\lambda|$ power law tail (25) we observed in the local DOS. Note that the above bound for the cumulative excursion time distribution again supports our intuitive “deep minimum in the first shell” picture: the average of $e^{-t/\tau^{\prime}}$ is dominated by traps with lifetimes $\tau^{\prime}\gtrsim t$, i.e. by the deepest minima surrounding the initial trap.

Before showing numerical results we comment briefly on the short time behaviour of $F(t)$. This is determined by the average escape rate $\int\mathrm{d}\tau^{\prime}\,\rho_{\tau}(\tau^{\prime})/\tau^{\prime}$ of the first neighbour traps, which have random depths, multiplied by the probability $1/c$ (again for the random regular graph ensemble) of making the first jump from such a neighbour back to the initial trap. Overall this yields $F(0)=T/(c\,(T+1))$. This constant translates into a $1/s$ power law in Laplace space for $s\gg 1$; such a power law also appears in $\hat{P}_{\tau}(s)$ due to $P_{\tau}(0)=1$.

Our numerical results for the excursion time distribution and the related quantities are displayed in figure 6: on the left we have the staying probability (red line), the return probability (blue line) and the excursion time distribution (green line) in Laplace space for the random regular graph ensemble with $c=5$, $T=0.8$ and considering departing traps with lifetime $\tau=2$. $\hat{P}_{\tau}(s)$ and $\hat{F}_{\tau}(s)$ are computed using the population dynamics results for the local DOS, while $\hat{S}_{\tau}(s)$ is given explicitly by Eq. (35). Note that all these quantities decrease as $1/s$ for large $s$ as expected. The inset shows the small $s$ behaviour where the horizontal lines correspond to the values $\hat{P}_{\tau}(0)$, $\hat{F}(0)$ and $\hat{S}_{\tau}(0)$ derived above. The plot on the right has the excursion time distribution $F(t)$ evaluated as the inverse Laplace transform of the data in the left plot, together with the results from direct simulations on the infinite $c$-regular tree (blue lines); we show simulations for multiple $\tau$, which produce identical results as expected. These numerical results match the approximation (40) well, see the black dashed line labelled $2^{\mathrm{nd}}$ order approximation. The inset clarifies the behaviour of $F(t)$ for $t\to 0$, with the horizontal line indicating the asymptote $F(0)$.

We stress that our results on the long time power law behaviour of the return probability and excursion time distribution are robust to changes in network topology, as long as this exhibits a locally tree-like structure. Our deep minima argument continues to apply then, with returns from the first neighbour shell surrounding the initial node giving the dominant contribution.

We return finally to the return probability $P_{\tau}(t)$ as shown in Fig 5. As discussed, this quantity initially decays exponentially in $t/\tau$, and the range where this decay is seen expands without bound as $\tau\to\infty$. The decay constant is somewhat slower than the staying probability $S_{\tau}(t)$ would suggest, however. To understand this, one can focus on the $t/\tau$ scaling of $P_{\tau}(t)$ by considering in Laplace space $\tau^{-1}\hat{P}_{\tau}(s)$ for $s=\sigma/\tau$. From (36) this is just $\tau^{-1}\hat{P}_{\tau}(\sigma/\tau)=[1+\sigma-\hat{F}(\sigma/\tau)]^{-1}$. In the large $\tau$-limit that we are interested in, $\hat{F}(\sigma/\tau)\to\hat{F}(0)=1/(c-1)$ so that $\tau^{-1}\hat{P}_{\tau}(\sigma/\tau)\to[(c-2)/(c-1)+\sigma]^{-1}$. This implies in the time domain in the same limit that

[TABLE]

which is consistent with the numerical results shown in Fig. 5-left. The exponential decay resembles that of the staying probability $S_{\tau}(t)=\exp(-t/\tau)$, but is somewhat slower except in the mean field limit $c\to\infty$. The slowing down arises from the fact that the system can return an arbitrary number of times $n$ to the initial trap, and the sum of all the contributions from $n=1,2,\ldots$ returns just conspires to produce an exponential return probability decay with a smaller decay rate. An intuitive physical explanation can be formulated as follows: when $\tau\gg 1$, excursions take negligible time compared to $\tau$. Since $\hat{F}(0)$ is the return probability, the probability of escaping in any attempt is only $1-\hat{F}(0)$. Therefore the rate of escape is $(1-\hat{F}(0))/\tau=(c-2)/(c-1)\tau^{-1}$.

6 Conclusion

In this paper we have studied a model for the evolution of glasses in configuration space: the dynamics takes place on a random network whose nodes represent the energy minima (or traps) of the system. These are characterised by the number of nearest neighbours $k$ and an average lifetime $\tau$. The latter is a quenched random variable, power law distributed, whose average diverges below the glass transition temperature $T=1$. Our focus was on the spectrum of eigenvalues $\lambda$ of the master operator governing the dynamics, and in particular on the local density of states $\rho(\lambda|\tau)$, i.e. the contribution to the spectrum of relaxation rates from all traps with a fixed lifetime $\tau$.

We employed the cavity method to exploit the tree-like structure of infinite random networks, and computed numerically the local DOS using a population dynamics algorithm. The cavity construction also allowed us to perform a simple analytical approximation that provides a very good match to the exact numerical results: the local DOS shows a small-$|\lambda|$ power law tail governing the long time dynamics, specifically $\rho(\lambda|\tau)\sim\tau|\lambda|^{T}$ for $|\lambda|\ll 1/\tau$. This result is robust to changes in the network topology as long as a locally tree-like structure is retained; in this class of networks we considered here random $c$-regular, scale-free and Erdös-Rényi graphs. The power law tail of the local DOS is associated, in the time domain, with the distribution $F(t)$ of excursion times $t$ away from some initial trap. We found $F(t)\sim t^{-(T+1)}$. This can be seen as responsible for the long time behaviour of the probability to return to traps of depth $\tau$, which is $P_{\tau}(t)\sim\tau t^{-(T+1)}$ for $t\gg\tau$.

We showed that the above dynamical properties arise as a combined effect of a sparse (loop-free) configuration space connectivity and the quenched trap depth disorder: when these features are considered separately the local DOS becomes $\delta$-shaped for small $|\lambda|$, implying an exponential decay of $F(t)$ and $P_{\tau}(t)$. In more detail, our analysis indicates that the most likely way for the system to return to the departing trap at some time $t\gtrsim\tau$ is to spend most of the interim stuck in a neighbouring trap with lifetime of order $t$ or larger. Returns from more distant minima are possible but become exponentially less probable with distance. For $t\lesssim\tau$, instead, $P_{\tau}(t)$ is dominated by the probability of having never left the original trap, $S_{\tau}(t)=e^{-t/\tau}$; in the mean field limit of infinite connectivity, $S_{\tau}(t)$ is in fact the only contribution to $P_{\tau}(t)$. Finally, the exponential shoulder of $P_{\tau}(t)$ dominates the average return probability $P(t)=\int\mathrm{d}\tau\rho_{\tau}(\tau)P_{\tau}(t)$, leading to the mean field scaling $P(t)\sim t^{-T}$, in accordance with the small $|\lambda|$ tail observed in the total DOS [26]. This tells us that the long time dynamics of the trap model on random networks is of mean field kind on average only, while the analysis of more detailed phenomena reveals the effects of the network structure.

We conclude this paper with two remarks pointing to future directions. The first one concerns the analysis of return probabilities in the Anderson model on random regular graphs put forward in [43]. In that paper the authors are able to write the return probability in terms of eigenfuction correlators evaluated at distance zero on the graph, which they can compute with a population dynamics method. The spatial correlation length of the eigenvector entries interests us too, in particular that of the slow decaying modes, and we aim at investigating this and related properties in the future using a similar approach. In the classical context the return probability to a node $i$ is given by $P_{ii}(t)=(e^{-\mathbf{M}t})_{ii}$, while in the quantum case one has $P_{ii}(t)=|\langle i|e^{-\mathrm{i}\mathcal{H}t}|i\rangle|^{2}$, where $\mathcal{H}$ is the Hamiltonian of the system. The representation that uses the eigenfunction correlators then holds true in our case for $(e^{-\mathbf{M}t})_{ii}^{2}$, which represents the probability that two independent replicas of the system starting out at the same node both return to that node at time $t$. This quantity connects directly to our second remark, on the question of the characterization of the low temperature phase of the trap model in terms of replica symmetry breaking (RSB) in trajectory space. While the 1D lattice considered by Ueda and Sasa in [44] exhibits an RSB phase for $T<1$, this is not the case for tree-like networks, because trajectories always depart from each other in such infinite-dimensional structures. However, a possible generalization of the trajectory RSB idea to random networks is to consider closed paths only: the distribution of excursion times mentioned above exhibits a diverging mean for $T<1$, a fact that we would conjecture should be associated with an RSB phase in the space of closed trajectories. Work in this direction is in progress.

7 Acknowledgements

The authors acknowledge funding by the Engineering and Physical Sciences Research Council (EPSRC) through the Centre for Doctoral Training “Cross Disciplinary Approaches to Non-Equilibrium Systems” (CANES, Grant Nr. EP/L015854/1).

Appendix A Mean field and random walk limits

In section 3 we saw that the power law behaviour of the local DOS close to the ground state ($\lambda=0$) arises as a combined effect of limited connectivity and trap depth disorder. When considering the two features separately, i.e. either taking the limit $T\to\infty$ (RW) or considering a fully connected configuration space (MF), the local DOS exhibits a $\delta$-peak at some $\lambda^{*}(\tau)$, and vanishes in the region $|\lambda|<|\lambda^{*}|$. In this appendix we discuss in more detail the local DOS of the relevant (MF and RW) limits; in particular we derive the value of $\lambda^{*}$ for the various cases.

The easiest way to obtain a fully connected network is to consider the random $c$-regular graph ensemble and then take the limit $c\to\infty$. Imposing $p_{k}=\delta_{c,k}$ in (22) we get

[TABLE]

which for $c\gg 1$ becomes

[TABLE]

where $\tilde{\omega}=\int\mathrm{d}\omega\,p(\omega)\operatorname*{\mathrm{i}}\omega/(\operatorname*{\mathrm{i}}+\omega)$. Using the self-consistency equation (19) for $p(\omega)$ we get

[TABLE]

which implies $\tilde{\omega}=\operatorname*{\mathrm{i}}$ when $c$ is large. Substituting into (43) then yields

[TABLE]

It follows that $\lambda^{*}=-1/\tau$, and the delta peak at this location is the only contribution to the local DOS. This is displayed as a vertical line (red) in figure 7-left, which also shows the mean field local DOS evaluated by averaging results obtained from direct diagonalizations of systems with different size (shades of green). Note that finite size effects are visible in the tails away from the peak. Intuitively, the local DOS in the mean field limit has to be a $\delta$-function centred in $-1/\tau$ as the return probability and the staying probability are the same in this case, and are given by the exponential function in (29).

The local DOS in the random walk limit can be analysed only if the value of $\tau$ is fixed before $T$ is sent to infinity; taking $T\to\infty$ first would make all $\tau_{i}=1$. For the random $c$-regular graph ensemble, the local DOS is given by the approximation scheme described in section 3.2 – which is exact in this case – evaluated at first order. We have:

[TABLE]

with $\tilde{\omega}=\operatorname*{\mathrm{i}}\bar{\omega}/(\operatorname*{\mathrm{i}}+\bar{\omega})$, and $\bar{\omega}$ given by the solution of the infinite temperature cavity equation

[TABLE]

The physical solution is the one with positive real part. This defines the precision of the Gaussian cavity measure (see main text before equation (18)) for an infinite random $c$-regular network in the limit $T\to\infty$, where the distribution of cavity precisions simplifies to $p(\omega)=\delta(\omega-\bar{\omega})$. Equation (46) can be written more conveniently as

[TABLE]

where $\varepsilon^{\prime}=\varepsilon\tau$ is a rescaled version of $\varepsilon$, and $\tilde{\omega}_{\mathrm{I}/\mathrm{R}}$ denote the imaginary/real part of $\tilde{\omega}$, respectively. In the limit $\varepsilon\to 0$ these read

[TABLE]

[TABLE]

Note that $\tilde{\omega}_{\mathrm{R}}$ is positive for $\lambda\in\mathcal{R}=(-1-\Delta_{c},-1+\Delta_{c})$ with $\Delta_{c}=2(c-1)^{1/2}/c$; outside of this region it vanishes.

It is instructive to consider first the case of $\tau=1$, which gives the total DOS of the random regular graph ensemble in the random walk limit. For this simple case one can directly set $\varepsilon^{\prime}=0$ to get $\rho^{\mathrm{1A}}(\lambda|1)=\tilde{\omega}_{\mathrm{R}}/[\pi(\tilde{\omega}_{\mathrm{R}}^{2}+(\lambda+\tilde{\omega}_{\mathrm{I}})^{2})]$, which when worked out explicitly is a scaled and shifted Kesten-McKay law [26].

For general $\tau>1$ one sees that the local DOS still has a contribution for $\lambda\in\mathcal{R}$ but this becomes increasingly suppressed as $\tau$ increases, scaling as $1/\tau$ for large $\tau$ from (48). This is compensated for by an additional contribution outside of $\mathcal{R}$, where again from (48) but now with $\tilde{\omega}_{\mathrm{R}}=0$ one has $\rho^{\mathrm{1A}}(\lambda|\tau)=\delta(\lambda+\tilde{\omega}_{\mathrm{I}}(\lambda)/\tau)$, which is a delta peak at a location $\lambda^{*}$ determined by $-\lambda^{*}=\tilde{\omega}_{\mathrm{I}}(\lambda^{*})/\tau$. For large $\tau$, $\lambda^{*}$ becomes small so that asymptotically $\lambda^{*}=-\tilde{\omega}_{\mathrm{I}}(0)/\tau$, hence from (50),

[TABLE]

This delta peak is the dominant contribution to the local DOS for large $\tau$, where using (28) it gives the return probability $P_{\tau}(t)=\exp(-(t/\tau)(c-2)/(c-1))$ as derived by another route in (41) in the main text.

The right plot in figure 7 shows the local DOS for the random regular graph ensemble with $c=5$, $T=0.8$ and $\tau=200$, obtained from direct diagonalizations (green lines) and using the exact cavity result (48) (blue line). The vertical dashed line corresponds to the value of $\lambda^{*}$ given by (51), while the red curve shows the case of $\tau=1$ discussed above (with support only on $\mathcal{R}$). Note the strong $\varepsilon$ dependence of the results for $\lambda\neq\lambda^{*}$ and $\lambda\notin\mathcal{R}$, which is consistent with the expectation that $\rho(\lambda|\tau)$ vanishes in these regions in the limit $\varepsilon\to 0$.

In this appendix we have only considered the case of random regular graphs, and we showed that the local DOS is composed of a continuous part with support on $\mathcal{R}$, and a $\delta$-peak whose location scales with $1/\tau$ for large $\tau$. More disordered network topologies would have the same qualitative behaviour, however. In particular, they would show a network-dependent regime of fast relaxation rates $-\lambda$, similarly to what happens in the case of finite connectivity and finite temperature discussed in section 3 (see also [26]-section 5 for results on the total DOS). A $\delta$-peak would again appear in the small $|\lambda|$ regime for large $\tau$.

Appendix B Second shell approximation

In this appendix we derive the estimate (25) for the power law behaviour of the local DOS close to the ground state ($\lambda=0$), which constitutes one of the central results of this work. We start by re-writing (23) in a more explicit form and for a general degree distribution:

[TABLE]

where

[TABLE]

and $\bar{\omega}$ is the physical solution of (47), which for small $\lambda_{\varepsilon}$ can be approximated as

[TABLE]

Substituting this approximation into the definition of $f(\tau_{l})$ leads to

[TABLE]

Since we are interested in the limit $\varepsilon\to 0$ we neglect the term in the denominator multiplying $\varepsilon^{2}$. The error is significant only when $c^{2}[\tau_{l}+(c-2)^{-1}]^{2}\gg\varepsilon^{-2}$, in which case the second term of the denominator becomes very large ($\gg 1$) and the resulting contribution to $\rho^{\mathrm{2A}}(\lambda|\tau)$ is negligible. We can therefore write the term in squared brackets in (52) as

[TABLE]

The real term on the right hand side can again be viewed as a rescaled (and still positive) $\varepsilon^{\prime}$. Once this is sent to zero we obtain

[TABLE]

with

[TABLE]

Since we are interested in the regime $|\lambda|\ll 1/\tau$, the term $\lambda c\tau$ in the $\delta$-function of equation (57) can be discarded. We observe then that the only non-zero contributions to the local DOS are given by those combinations of $\{\tau_{l}\}$ that result in $\sum_{l}y_{l}=k$. To understand when this happens let us set

[TABLE]

so that $y_{l}=(a-b\tau_{l})^{-1}$. From the distribution of lifetimes (3) we obtain

[TABLE]

for $y_{l}<0$ or $y_{l}>1/(a-b)$. Note that as $|\lambda|\to 0$, also $b\to 0$. In this limit, $\rho_{y}(y_{l})$ goes to zero everywhere except in a region of order $b$ around $y_{l}=1/(c-1)$. We can now approximate the probability distribution of $Y=\sum_{l=1}^{k}y_{l}$. This drops by a factor $b^{T}$ for each of the $y_{l}$ that is away of $1/(c-1)$, so the most likely way to realize $Y=k$ is to have $k-1$ of the $y_{l}$ equal to $1/(c-1)$, and only a single one, say $y_{1}$, equal to $k-(k-1)/(c-1)$. Note that this happens when $\tau_{1}\sim 1/|\lambda|$, i.e. if there is a single deep minimum in the first neighbouring shell of the departing trap with lifetime as large as $1/|\lambda|$. So we have $p(Y=k)\simeq k\rho_{y}(k-(k-1)/(c-1))$, where the factor $k$ arises because any of the $y_{l}$ could be the large one. Finally, from (57) (with $\lambda c\tau\to 0$) we see that

[TABLE]

which is the same as equation (25) in the main text, with $C(c,T)=T\bar{\alpha}_{k}c^{T+1}(c-1)^{1-T}$ and

[TABLE]

The simplest case of the random regular graph ensemble is obtained by imposing $p_{k}=\delta_{c,k}$ in the last equation, which leads to $C(c,T)=Tc(c-1)^{T-1}(c-2)^{-(T+1)}$.

Appendix C Simulated dynamics

The data labelled “simulations” shown in the figures 4 and 6-right have been collected by using a version of the stochastic simulation algorithm (SSA) that generates dynamically an infinite tree. The SSA became popular after Gillespie applied it to the study of chemical reactions [45], and for this reason it is also known as “Gillespie algorithm”. Its general idea is to implement the stochastic evolution of a system on a discrete state-space as follows: for each state $i$, a random waiting time $\mathrm{d}t_{i}$ is sampled from the exponential distribution $r_{i}^{\mathrm{ex}}\exp(-r_{i}^{\mathrm{ex}}t)\theta(t)$, where $r_{i}^{\mathrm{ex}}=\sum_{j(\neq i)}r_{ji}$ is the total exit rate from state $i$, and $r_{ji}$ is the transition rate from state $i$ to state $j$. Then, the next state $j$ is chosen with probability $r_{ji}/r_{i}^{\mathrm{ex}}$. This process is repeated until the total time $\sum_{i}\mathrm{d}t_{i}$ exceeds some $t_{\mathrm{max}}$ that sets the maximum running time of the simulation.

In our case, the states are represented by the nodes of the network. These have four main attributes that we track in the simulation: the energy $E_{i}$, the degree $k_{i}$, the distance (from the origin) $d_{i}$ and the (number of) visits $n_{i}$. The time spent in a given node $i$ and the next node visited $j$ are defined by the routine gillespietrap$(i)$. This works as follows:

compute the total exit rate $r_{i}^{\mathrm{ex}}=\sum_{j\in\partial i}r_{ji}=k_{i}e^{-\beta E_{i}}/c$ ($\partial i$ indicates the neighbourhood of node $i$); 2. 2.

compute the waiting time $\mathrm{d}t$ by sampling from $p_{i}(t)=r_{i}^{\mathrm{ex}}\exp(-r_{i}^{\mathrm{ex}}t)\theta(t)$; 3. 3.

select the next node $i^{\mathrm{new}}$ randomly from the $k_{i}$ neighbours: $r_{ji}/r_{i}^{\mathrm{ex}}=1/k_{i}$; 4. 4.

return $i^{\mathrm{new}}$ and $\mathrm{d}t$.

The quantities of interest, such as the current state or the distance from the origin, are measured at times defined by a time-grid with $n_{\mathrm{times}}$ values in the range $[0,t_{\mathrm{max}}]$. In order to simulate the evolution on an infinite tree, the algorithm has to create the network structure on the fly. This can be done as follows:

start from a node with $k_{0}$ leaves ($k_{0}$ is sampled from $p_{k}$), energy $E_{0}$, distance $d_{0}=0$ and visits $n_{0}=1$. Each leaf $j\in\partial i$ has a random energy sampled from $\rho_{E}(E)$, degree $k_{j}=1$, distance $d_{j}=1$ and visits $n_{j}=0$. This is the starting network configuration; 2. 2.

select the next node $j$ with gillespietrap$(i)$. If $n_{j}=0$, attach a new neighbourhood to $j$, taking into account that $j$ already has $i$ as neighbour. This is done by the routine newneighbourhood$(j)$. Then, the number of visits is set to $n_{j}=n_{j}+1$; 3. 3.

newneighbourhood$(j)$ assigns $k^{\mathrm{new}}_{j}-1$ leaves to $j$, with $k^{\mathrm{new}}_{j}$ sampled from $kp_{k}/c$, and so it replaces $k_{j}=1$ with $k^{\mathrm{new}}_{j}$. The new leaves $l\in\partial j\setminus i$ have random energies $E_{l}$ sampled from $\rho_{E}(E)$, distance $d_{l}=d_{j}+1$, degree $k_{l}=1$ and visits $n_{l}=0$.

Note that no loops, single nodes or disconnected components are created. The algorithm can also be used for running multiple copies of the dynamics in parallel, which is useful if one is interested in collecting data for a given realization of the disorder. The full structure of the implementation is explained in the following pseudo-code, where $x(m)$ indicates the position of the $m^{\mathrm{th}}$ copy of the system, $t(m)$ its time, and there are $n_{\mathrm{copies}}$ in total.

Pseudo-code: simulations on the infinite tree.

set all $x(m)=0$ where [math] is the initial trap

set all $t(m)=0$

set all $x^{\mathrm{new}}(m)=x(m)$

create the starting network configuration

$\ldots$

for ($t=0$; $t<n_{\mathrm{times}}$; $t$++) do

     $\mathrm{time}=\mathrm{time\_grid}(t)$

     for ($m=1$; $m\leq n_{\mathrm{copies}}$; $m$++) do

         while ($t(m)\leq\mathrm{time}$) do

              carry out one transition

              $x(m)=x^{\mathrm{new}}(m)$

              $x^{\mathrm{new}}(m),\,\mathrm{d}t=\mathrm{\emph{gillespietrap}}(x(m))$

              if ($(x^{\mathrm{new}}(m)).\mathrm{\emph{visits}}==0$) then

                  newneighbourhood$(x^{\mathrm{new}}(m))$

              end if

              $t(m)=t(m)+\mathrm{d}t$

         end while

         collect statistics here

         $\ldots$

     end for

end for