Spectral properties of the trap model on sparse networks
Riccardo Giuseppe Margiotta Reimer Kühn Peter Sollich

Predicting Transportation Modes of GPS Trajectories using Feature Engineering and Noise Removal
    Mohammad Etemad1
Amílcar Soares Júnior1 Stan Matwin12

    Received: date / Accepted: date
===============================================================================================






One of the simplest models for the slow relaxation and aging of glasses is the trap model by Bouchaud and others, which represents a system as a point in configuration-space hopping between local energy minima. The time evolution depends on the transition rates and the network of allowed jumps between the minima. We consider the case of sparse configuration-space connectivity given by a random graph, and study the spectral properties of the resulting master operator. We develop a general approach using the cavity method that gives access to the density of states in large systems, as well as localisation properties of the eigenvectors, which are important for the dynamics.
We illustrate how, for a system with sparse connectivity and finite temperature, the density of states and the average inverse participation ratio have attributes that arise from a non-trivial combination of the corresponding mean field (fully connected) and random walk (infinite temperature) limits. In particular, we find a range of eigenvalues for which the density of states is of mean-field form but localisation properties are not, and speculate that the corresponding eigenvectors may be concentrated on extensively many clusters of network sites.






§ INTRODUCTION


Glasses are disordered materials that do not exhibit the structural periodicity of crystals but nonetheless possess the mechanical behaviour of solids. The most common way of making a glass is by quenching, i.e. cooling a viscous liquid so rapidly that crystallisation is avoided. The resulting system is called a supercooled liquid. The quench brings the molecules of the material into a configuration where the typical time needed to rearrange them is so long that the structure of the liquid appears frozen. The system falls out of equilibrium in the sense that the relaxation time becomes of the order of the observation time window.
The resulting extremely slow evolution is called glassy dynamics, and the transition into the regime of very long relaxation times is referred to as the glass transition. Technically this phenomenon is not a real phase transition as there are no discontinuous changes in any physical property. Nevertheless, one can associate a critical temperature T_G to a certain liquid, below which the rate of change of e.g. volume due to a change in temperature is comparable to that of a solid. The value of T_G also depends on the rate at which the system is cooled: slower cooling allows the material to fall out of equilibrium at lower temperatures (allowing more time for configurational sampling). 


 Several theoretical approaches have been proposed to investigate the nature of the glass transition; the general discussion is presented in a recent review by Biroli and Berthier <cit.> (see also references therein). In spite of a sustained research effort dedicated to this problem, a full understanding of glasses has not been achieved yet. One of the most successful theories (based on a microscopic description) is the mode-coupling theory, which predicts a dynamical arrest in supercooled liquids associated with a power law divergence of the `slow' time scale <cit.>. Another important class of models in the context of glassy systems is that of spin glasses, where one generally starts from a Hamiltonian ℋ with disordered interactions and derives the thermodynamic properties of the spin system and its dynamical behaviour by averaging over the disorder <cit.>.
 
A further useful angle of attack on the glass problem focusses on the dynamics in configuration-space. The energy landscape of a glass is typically very rugged, consisting of many local minima (metastable states) separated by energetic barriers, and a global minimum (the crystalline equilibrium state) that is kinetically extremely difficult to reach. One can then think of this energy landscape as a set of basins of attraction that act as `traps' for the dynamics: during its evolution towards equilibrium, the system jumps between local minima at rates that decrease strongly with decreasing temperature. Based on this picture, several studies have been developed, focusing on various aspects of glassy dynamics in configuration-space. These range from investigations of the potential energy landscape, in particular the structure and distribution of minima and energetic barriers between them <cit.>, to simplified models that describe the evolution between traps at a more phenomenological level <cit.>. 

Interestingly, once the description of the configuration-space dynamics has been simplified to motion among traps without internal structure, it is directly related to the research field of stochastic processes on networks. The structure of the energy landscape and the relative positions of neighbouring minima define a network of allowed transitions: the system can only jump between traps that are linked within this network, i.e. between minima that are close in the configuration-space. Therefore methodology and results from network theory <cit.> can be applied to understand the phenomenology of glasses. In particular, information about the energy landscape can be used to model the time evolution of the system as a Markov process on the network of minima, with rates depending on the relevant energy barriers. Mathematically, the problem thus turns into solving a master equation for the time-dependent probability distribution that describes the position of the system in configuration-space.

One of the simplest and most successful descriptions that belong to this framework is the trap model by Bouchaud and others <cit.>. The transition rates are assumed to depend on the depth of the departing trap j only, not on the arrival trap i, and have the Arrhenius-like form


    r_ij=1/N e^-β E_j

Here β is the inverse temperature, E_j>0 is the trap depth and N is the size of the network (the number of traps). Every transition effectively involves activation to the top of the energy landscape, where all the energy barriers are located, and then falling into a new state that is chosen randomly among all the minima. The latter assumption implies that this model postulates a mean field (fully connected) network structure. It is easy to show that, for an exponential density of trap depths ρ_E(E)=e^-E, a glass transition occurs at finite temperature. More specifically, below T_G=1 the equilibrium probability distribution across trap depths becomes non-normalisable. The exponential form of the trap density of states can be motivated from several points of view, e.g. the mean-field replica theory of spin-glasses <cit.>, the random energy model <cit.>, or phenomenological arguments in the context of supercooled liquids <cit.>. Also, following an extreme value statistics argument, one might expect that deep minima of potential energy landscapes are described by the Gumbel distribution, whose tail is indeed exponential <cit.>

The simple expression for the transition rates and the fully connected network structure allow one to solve the master equation for the model described above analytically. This is simplest in the Fourier-Laplace domain, from where the behaviour in the time domain can then be extracted straightforwardly <cit.>. Trap models have been also studied in Euclidean space, where the system jumps between the nodes of a regular lattice; see for example <cit.> or the work by Ben Arous and collaborators <cit.>. Variants include branching phenomena <cit.> and walks on positive integers <cit.>, though this is less plausible when modelling configuration-space dynamics. The first extension to a trap model on a network was considered relatively recently by Baronchelli et al, who used a simple heterogeneous mean field approximation to study the dynamics. This assumes that the probability to find the system on a certain site only depends on the degree (i.e. on the number of adjacent nodes) of the site. It therefore has to postulate that the trap depth at any site is fixed fully by its degree <cit.>. Numerical results do indeed show some correlation between trap depth and degree <cit.>, though the relation between the two is far from deterministic.

In this work we extend the analysis of the trap model to dynamics on generic (random) networks with sparse inter-trap connectivity. Compared to <cit.> we develop a more flexible approach to the modelling of glassy configuration-space dynamics that allows an arbitrary (deterministic or stochastic) relation between trap depth and node degree. Within this general scenario we then consider the simplest case where trap depths are uncorrelated with degrees. 

For disordered energy landscapes with sparse connectivity a direct analytical solution of the dynamics is not possible in either frequency or time domain; we therefore tackle the problem via the spectral properties of the master operator, which are key in determining the dynamics of the system. Specifically we calculate the density of states (DOS), which gives the spectra of relaxation rates of the system, and the localisation properties of the eigenvectors, measured using the inverse participation ratio (IPR). We develop a general cavity method for this purpose, leading in the infinite system size limit to an integral equation that can be solved numerically via a population dynamics algorithm. Technically, the approach follows analogous applications of the cavity method to the spectral analysis of symmetric random matrices; see e.g. <cit.> or <cit.> for a rigorous discussion, and <cit.> for related work on heavy-tailed random matrices. Based on the DOS and IPR, we are able to obtain insights into the relevant time scales and time regimes of the system. However, we do not have access to some time-dependent objects like correlation functions, which are the main quantities of interest within the literature of trap models. Our analysis will therefore be different from that of previous works <cit.>, and limited to describing the dynamics in terms of the static quantities mentioned above.

This paper is organised as follows. In  section <ref> we define the general set-up of the problem. In section <ref> we discuss by way of background the localisation of the ground state as a function of the temperature, and summarise the known results for the mean field and random walk limits. In section <ref> we address the general case of trap model dynamics on networks with finite connectivity and at finite temperature, and we propose a simple analytical approximation for the DOS. Also, we explain how the parameter ε that appears in the evaluation of the DOS can be exploited as a detection tool for localisation transitions within the spectrum of the system. We then use these methods to extract dynamical properties of the trap model on random regular graphs, where all nodes have the same degree.
In section <ref> we extend the analysis to other network topologies including scale-free graphs. Section <ref> summarises our conclusions and outlines perspectives for future work.




§ PROBLEM SET-UP


The general setting of the problem is the following: we consider a continuous-time Markov process defined on a network of N nodes that represent the energy states accessible by the system, i.e. the minima of the potential energy landscape or simply the traps. The starting point is then given by the master equation for the probability distribution 𝐩(t)=(p_1(t),…,p_N(t)), where p_i(t) is the probability to find the system in trap i at time t:

    ∂_t 𝐩(t)=𝐌𝐩(t)

The master operator 𝐌 has the following structure:

    M_ij=c_ijr_ij     M_ii=-∑_j≠ i M_ji

where c_ij=c_ji=1 if nodes i and j are connected and c_ij=0 otherwise, also c_ii=0 (there are no self-loops), and r_ij is the transition rate from node j to node i. Note that ∑_j M_ji = 0, which ensures that probability is conserved.
We shall now specify the trap-depth distribution, the transition rates and the network topology. We assume


  * exponentially distributed energies: E∼ρ_E(E)= e^-E, E≥ 0;

  * random graph structure: the c_ij are sampled from a random graph ensemble with finite connectivity. The simplest case for our purposes is one where all those graphs have equal probability for which each node i is connected to exactly c others; the probability distribution of node degrees k_i=∑_j c_ij is then p_k=δ_c,k.  Samples that belong to this ensemble are called random regular graphs (RRG). In this work we are interested in the case of c≥ 3. This condition ensures that the fraction of nodes outside the giant component vanishes in the large system limit <cit.>, which also implies that the configuration-space is connected, therefore ergodic. We will develop our theory for  general random graph ensembles where the degree distribution is constrained to some p_k; in that case c is defined as the average degree c=∑_k k p(k).


  * Bouchaud transition rates: r_ij= e^-β E_j/c ≡ r_j. The total escape rate from node j is defined as r̂_j = ∑_i c_ij r_ij and can be written in terms of the node degree as r̂_j=k_j r_j. We will find it useful to define τ_j=(c r_j)^-1= e^β E_j. This gives the expected waiting time to exit from trap j, exactly so for regular graphs and up to a factor c/k_j in the general case. From the energy distribution we obtain that the τ_j have the distribution

    ρ_τ(τ)=Tτ^-(T+1)  with  τ∈[1,∞)

 which implies that the average waiting time ⟨τ⟩ diverges for T<1, signalling the occurrence of glassy dynamics and aging.

With the above assumptions the master operator 𝐌 is a sparse random matrix. It will be important to bear in mind that two sources of randomness come into play: the disorder in the trap depths {E_i} and in the inter-trap connectivity {c_ij}. Accordingly we also have two different notions of distance that are relevant for this model: the distance in energy, i.e. the energy difference among the minima, and the distance on the graph structure. As we will see in the following sections, these notions of distance play different roles, depending on the case being studied, with regards to their relevance for the degree of localisation of the eigenstates.


The formal solution of equation (<ref>) is given by 

    𝐩(t)=∑_αe^λ_αt(𝐰_α,𝐩(0))𝐮_α




where λ_α, 𝐰_α, 𝐮_α are respectively the eigenvalues, left eigenvectors and right eigenvectors of 𝐌, indexed by α=0,1,…,N-1, and (𝐰_α,𝐩(0)) denotes the scalar product between the left eigenvector 𝐰_α and the initial probability distribution 𝐩(0). In the following we will refer to r_α = -λ_α as the relaxation rates of the system, and write r for a generic relaxation rate. If the network is connected, there is a single vanishing eigenvalue λ_0=0. All other λ_α must then have negative real part so that the corresponding modes 𝐮_α make a contribution to 𝐩(t) that is exponentially suppressed over time. In the long-time limit only 𝐮_0 survives, which is the equilibrium Boltzmann distribution of the system associated with the ground state λ_0. The corresponding left eigenvector is 𝐰_0 = (1,… ,1). So

    lim_t→∞𝐩(t)=𝐮_0=𝐩_eq=1/Z(e^β E_1,…,e^β E_N)

Within the present formulation the energies are positive as they represent the depth of each trap, so e^β E_i is the correct Boltzmann weight for node i.

The evolution of the probability at finite t depends on the spectral properties of the master operator. In particular, slowly decaying modes govern the long-time behaviour of the system, and solving the master equation amounts to diagonalising 𝐌. This operation can be performed analytically only for a few special cases presented in section <ref>. However, information about the spectrum and the localisation properties of 𝐌 can still be obtained in the large system limit; we use the cavity method <cit.> for this purpose. 
This method links 𝐌 to the inverse covariance matrix of a Gaussian distribution, and therefore requires the symmetrised form of the master operator

    𝐌^s=𝐏_eq^-1/2𝐌𝐏_eq^1/2

or in components M^s_ij=r_i^1/2M_ijr_j^-1/2,
where we have introduced a diagonal matrix 𝐏_eq with (𝐏_eq)_ii=p_i^eq∝ r_i^-1. This transformation preserves the eigenvalue spectrum of 𝐌, implying that the associated eigenvalues are real, as 𝐌^s is real and symmetric. We note that the diagonal elements of 𝐌 remain unchanged: (𝐌^s)_ii=(𝐌)_ii. The eigenvectors 𝐯_α of 𝐌^s are given by 𝐯_α = 𝐏_eq^-1/2𝐮_α = 𝐏_eq^1/2𝐰_α. Physically, the symmetry of 𝐌^s means that the dynamics we are considering obeys detailed balance with respect to the Boltzmann steady state. 

Our study aims to predict the statistics of the eigenvalues and eigenvectors of 𝐌^s. The first quantity of interest is the density of states (DOS), i.e. the fraction of eigenvalues lying between λ and λ + dλ, defined as ρ(λ)dλ with

    ρ(λ)= 1/N∑_α=0^N-1δ(λ-λ_α)

We average the DOS over random samples, finitely sized, and assume self-averaging in the thermodynamic limit. The DOS is crucial as it defines the time scales {|λ_α^-1|} of the dynamics, or, more precisely, it gives the full spectra of relaxation rates {r_α}.  It is essential to keep this in mind as, for consistency, we will present the results in terms of the DOS throughout the paper, and occasionally remind the reader of the simple relation r_α = -λ_α.

The second key quantity that we are interested in is the degree of localisation of the eigenstates, which carries information about their ability to contribute to the transport properties of the system across the network: localised modes can only contribute to local probability-flows. The rationale behind this is clear: assuming p_i(0)=δ_ij, then if the vectors {𝐮_α, 𝐰_α} are mostly localised, only a few terms in the sum on the r.h.s. of equation (<ref>) give a significant contribution to the probability distribution 𝐩(t), which should therefore spread only slowly over time away from the initial node j. In general, we expect that the ability of the system to explore the configuration-space depends on the degree of localisation of the eigenvectors of 𝐌. To quantify this we use the inverse participation ratio (IPR) defined as

    I_q(𝐯)=∑_i=1^N v_i^2q/(∑_i=1^N v_i^2 )^q∼ N^-t_q

where 𝐯=(v_1,…,v_N) is an eigenstate; if v is normalised the denominator equals unity. The exponent t_q
defines the scaling of I_q with N. We refer to <cit.> for a general introduction to the IPR and related quantities. In what follows we concentrate on the standard IPR with q=2. We can distinguish two extreme situations: if the “mass” of the eigenstate 𝐯 is evenly spread over all the states of the system, namely each element v_i is of order 1/√(N), then the eigenstate is delocalised and I_2(v)=O(1/N), t_2=1. If instead only a few elements of 𝐯 differ from zero, the eigenstate is localised and I_2(v)=O(1), t_2=0.

Interestingly, the localisation properties of eigenstates defined on random regular graphs (and random matrices) are studied also in the context of quantum many body systems – with similar terminology and methodology – where they are linked to the problem of ergodicity and equilibration dynamics <cit.>. 

A number of studies have looked instead at the localisation of the time-dependent probability distribution of trap models on lattices. Particularly interesting is the 1D case, which exhibits dynamical localisation where localisation properties differ between the aging regime and the final Boltzmann distribution <cit.>. Flegel and Sokolov analysed this phenomenon using the spectral properties of the master operator <cit.>; they trace the non-equilibrium value of the IPR during aging back to the eigenvector statistics, while the eigenvalue statistics only make a minor contribution. Dynamical localisation is also discussed in the context of statistical mechanics of trajectories <cit.>, which represents another interesting approach to describing the glass transition in terms of configuration-space evolution.

The model we study, which is a Markov process on a random graph with Bouchaud transition rates, is described by two main parameters: the temperature T and the mean connectivity c. As depicted in figure <ref>, there are two obvious limits that can be considered: the mean field (MF) c→∞ limit, where the network structure becomes trivial and only the disorder in energy is present, i.e. there are “glassiness effects” only, and the T→∞ limit, where the trap depths become irrelevant and the system effectively performs a random walk (RW) among neighbouring traps. The point shown in the (1/c, 1/T) plane in figure <ref> represents our model with finite connectivity, at finite temperature. This general case can be thought of as a combination of the two limiting situations of mean field and random walk. In this work we illustrate how, for a system with finite c and T, quantities such as the DOS and the average IPR of eigenstates have attributes that arise by a non-trivial combination of the corresponding MF and RW limits.




§ GROUND STATE AND LIMITING CASES




 §.§ Ground state: λ=0


The eigenvector 𝐮_0 represents the equilibrium probability distribution of the system,
𝐩_eq = lim_t→∞𝐩(t). This is independent of the network structure and its statistics depend only on the energy distribution ρ_E(E). One can assess the degree of localisation of the equilibrium distribution via the IPR of either 𝐮_0 or its symmetrised analogue 𝐯_0. 
Explicitly, these are proportional to 

    𝐮_0=𝐩_eq∝ (e^β E_1, …, e^β E_N) 
        𝐯_0=𝐩_eq^s∝ (e^βE_1/2, …, e^βE_N/2)

The localisation of the ground state depends on whether the Boltzmann weights are concentrated on the deepest traps or not. Since the energies are randomly allocated to the vertices of the network, its topology will not affect the IPR of the equilibrium distribution; the distance in energy is the only relevant one here. From the definition of the IPR, we get for the ground state of the symmetrised master operator

    I_2(𝐯_0) = ∑_i e^2β E_i/(∑_i e^β E_i)^2∼N∫^N^β_1 dττ^1-T/(N∫^N^β_1 dττ^-T)^2
    ≃ N^-1  if   T>2  
     N^-2 + 2/T  if   1<T<2  
     N^0   if   T<1

where the cutoff N^β derives from the extreme value statistics of the distribution ρ_τ <cit.>: the k largest waiting times of N samples τ_N+1-k<…<τ_N fall in the range [τ_N+1-k,∞), therefore the fraction k/N is of the order of the area under ρ_τ over this range, which for k=1 gives 1/N ≃τ_N^-T = τ_max^-T.
From (<ref>), the result for the non-symmetrised version is the same except for the replacement of T by T/2.
We focus on the symmetrised case as this is the most sensible from the random matrix perspective that we use. The symmetrised eigenvectors are also the natural objects to appear in our cavity approach, which starts from a (complex) Gaussian distribution and hence requires a symmetric covariance matrix as input (see section <ref> and references therein). While the symmetrised eigenvectors do not describe the Markov process that obeys equation (<ref>), away from the ground state the symmetrisation is not expected to affect their localisation properties. In other words, a localised/delocalised symmetrised eigenvector should stay localised/delocalised also in its (either left or right) unsymmetrised form. We refer to appendix <ref> for further discussion and data showing that qualitative localization statistics for 𝐮_α, 𝐰_α and 𝐯_α are the same except in the small finite-size region of the crossover towards the ground state. We also note that the IPR of the symmetric ground state coincides with the measure of ground state localisation considered in previous works <cit.>. Finally, using symmetrised eigenvectors to calculate the IPR has an additional benefit: the characteristic temperature where the IPR ceases to be of 𝒪(1) coincides with the glass transition temperature that is known from the dynamics. Indeed, according to equation (<ref>), the ground state 𝐯_0 is localised below the glass transition T_G=1, delocalised for T>2, and has an intermediate behaviour for 1<T<2. Figure <ref> shows the exponent t_2(𝐯_0) of this prediction compared with the average t̅_2 (𝐯_0) = ⟨ -ln I_2(𝐯_0)/ln N⟩
of data taken from direct diagonalisations of the symmetrised master operator (hereafter labelled “numerics" in the plots). Note that in the limit N→∞  the IPR is of order unity for T<1, and drops to zero above, so that the intermediate temperature region (1<T<2) should also be regarded as delocalised. In the localised regime, an infinite-N calculation
shows that the O(1) value of the IPR is given explicitly by I_2(𝐯_0)=1-T for T<1 <cit.>, dropping to zero at T=1 in agreement with our result. 





 §.§ Random walk limit: T→∞


In the infinite temperature limit the dynamics is only affected by the graph topology, i.e. there are network effects only and the differences in energy depth among traps become immaterial. In this case the master operator coincides with its symmetrised form, and it simplifies to

    M_ij=c_ij/c-δ_ij

where the first term on the right hand side is the off-diagonal contribution (because c_ii=0).
Given that 𝐌=c^-1𝐀-𝐈 is directly related to the adjacency matrix 𝐀, its DOS can be deduced where that of 𝐀 is known.
For the case of a random regular graph, one obtains the DOS for N→∞ as a shifted and scaled Kesten-McKay law <cit.>:

    ρ(λ)=c/2π√(4 c-1/c^2-(λ+1)^2)/(1-(λ+1)^2)

which can alternatively be derived using e.g. the cavity construction explained below. 
For this graph ensemble all the eigenvectors are delocalised with high probability <cit.>. These results are illustrated in figure <ref>-left. 

The dynamics in this case has no glassy features, all local waiting times equal unity so that jumps occur at a constant rate, and the average distance from the initial node grows linearly with time. This is true because, at every jump, the particle has c-1 outward paths, and only one inward path pointing towards the starting node, thus the motion effectively resembles a 1D biased random walk. 



 §.§ Mean field limit: c→∞


In the infinite c limit the master operator reduces to that of the mean field (fully connected) case. There is no notion of space and only the distance in energy is relevant for the degree of localisation of the eigenstates. For such a fully connected system of size N  we have (approximating c=N-1≈ N, which is immaterial for N→∞)

    M_ij=e^-β E_j/N(1-δ_ij)-e^-β E_j(1-1/N)δ_ij

The eigenvalue equation in this case reads

    1/N∑_j≠ i e^-β E_ju_α, j - e^-β E_i(1-1/N)u_α, i=λ_α u_α,i

This can be written as 

    1/N∑_j e^-β E_ju_α, j - e^-β E_iu_α, i = λ_α u_α,i

and, as the first term is independent of i, one has

    u_α, i∝ (λ_α + e^-β E_i)^-1

Similarly, for the symmetrised case one obtains

    v_α, i∝e^-β E_i/2/λ_α + e^-β E_i

Note that for λ_α=0 these expressions recover the equilibrium distribution (<ref>) as they should. The solution (<ref>) is in agreement with <cit.> where the spectral properties of the mean field trap model are discussed extensively. From (<ref>, <ref>) we expect the IPR to be of order one for all λ≠ 0, for either their symmetrised or non-symmetrised forms. A simple argument for this localisation result in mean field is presented in appendix <ref>; we note here only that the eigenvector entries decay as a power law with energy difference to the “centre” of the eigenvector at  e^-β E_i≈-λ_α (the “centre” has to be understood here as defined on the energy axis). Returning to the eigenvalues, the condition for λ_α follows from equations (<ref>) and (<ref>) as

    ∑_je^-β E_j/N(λ_α+e^-β E_j) =1

which implies that there is an eigenvalue in each interval (-τ_i^-1, -τ_i+1^-1), assuming that the energies are ordered so that E_i<E_i+1 (recall that τ_i=e^β E_i). Therefore in the large N limit the DOS is given by

    ρ(λ)=∫dτ ρ_τ(τ) δ(λ+1/τ)

which gives

    ρ(λ)=T(-λ)^T-1

for -1<λ<0. These results are shown in figure <ref>-right. Note that the eigenvalue condition (19) and the interleaving of eigenvalues between the (negative) inverse trap lifetimes can also be seen from the fact that the mean field master operator (<ref>) is a diagonal matrix with a rank one perturbation <cit.>, as all elements in each column are the same except for those appearing on the diagonal.

Interestingly, all modes remain localised at any finite temperature, while glassiness manifests itself only for T<1, and even then only for the ground state. This circumstance has to be attributed to the slow decay of the mass of MF eigenvectors away from their localization centre, which does not impair the mobility of the particle. This agrees with the intuition that, in the absence of spatial structure, the particle is always able to reach any node of the network in finite time, as long as the average trapping time is finite, i.e. for T>1. 




§ FINITE CONNECTIVITY AND FINITE TEMPERATURE




 §.§ The cavity method


We now turn to the main contribution of our work: moving on from the two limiting scenarios discussed above, we study the general case of finite connectivity and finite temperature where the distance on the graph structure and separation between trap energies are both relevant. We do this by means of the cavity method, exploiting the fact that the random graphs we consider become locally treelike in the large N limit. The master operator in this case has the general structure (<ref>). In what follows we omit the index “s" and consider the symmetrised master operator only. The DOS of the matrix 𝐌 can be written in terms of the resolvent 𝐆(λ_ε) as

    ρ(λ)=lim_ε→ 01/π N∑_i=1^NIm G_ii(λ_ε)

where

    𝐆(λ_ε)=[λ_ε𝐈 - 𝐌]^-1

Here
λ_ε = λ - iε, with ε a small, positive quantity, while i is the imaginary unit; 𝐈 indicates the N× N identity matrix. To derive equation (<ref>) one replaces the delta distributions in (<ref>) with Lorentzians of width ε and takes the limit ε→ 0; this explains the origin of the small imaginary term in λ_ε. For a detailed description we refer to the original work of Edward and Jones <cit.>.
We define the complex Gaussian measure P(𝐱) as 

    P(𝐱)∝ e^-i/2𝐱^T 𝐆^-1𝐱= e^-i/2∑_i,j(λ_εδ_ij-M_ij) x_i x_j

with 𝐱=(x_1,… ,x_N). The diagonal entries of the resolvent are then given by the local variances

    G_jj= i∫dx_j  x_j^2 P(x_j)

where P(x_j) is the marginal distribution

    P(x_j)=∫∏_k≠ jdx_k  P(𝐱)


To make further progress we recall that the off-diagonal terms of the symmetrised master operator are M_jk=c_jkr_j^1/2r_k^1/2, while the diagonal terms are M_jj=-∑_k c_kjr_j. This gives  

    P(𝐱)∝ e^-i/2[λ_ε∑_j x_j^2 - ∑_jk c_jk (-r_j x_j^2 +  r_j^1/2 r_k^1/2 x_j x_k )]

Symmetrising r_j x_j^2 to (r_j x_j^2 + r_k x_k^2)/2
allows the term in brackets in the last sum to be written as a complete square. 
Equation (<ref>) can be further simplified by the change of variables y_j=x_jr_j^1/2, which has the benefit of confining the disorder from the transition rates r_j to the local terms:

    P(𝐲)∝ e^-i/2[λ_ε∑_j y_j^2/r_j + 1/2∑_jk c_jk (y_j-y_k)^2]
    = ∏_j e^-i/2λ_ε∑_j y_j^2/r_j∏_(jk)∈𝒢
    e^-i/2
    (y_j-y_k)^2

where the last product runs over all distinct edges of the graph 𝒢 defined by the inter-trap connectivity {c_ij}.

The core of the cavity approach is to decompose P(𝐲) into the factors involving a given node j, and the remaining factors. The latter define the cavity 
graph 𝒢^(j), where node j and all its connections have been removed from 𝒢, and a corresponding cavity distribution denoted P^(j)(·). This leads to the following equation for the marginal distribution P(y_j):

    P(y_j)=e^-i/2λ_εy_j^2/r_j∫d𝐲_∂ j  e^-i/2∑_k∈∂ j(y_j-y_k)^2 P^(j)(𝐲_∂ j)

where P^(j)(𝐲_∂ j) denotes the (complex) probability distribution of the variables {y_k} on the nodes that are neighbours of j on the graph 𝒢.
The cavity method is based on the assumption that the joint distribution P^(j)(𝐲_∂ j) factorises on 𝒢^(j) as

    P^(j)(𝐲_∂ j)=∏_k∈∂ jP^(j)(y_k)

This is exact if the original graph 𝒢 is a tree, because the cavity graph 𝒢^(j) then consists of disconnected branches. Sparse graphs do contain loops, but these have an average length of order ln(N) <cit.>. Intuitively, as the total number of nodes in the k-th coordination shell is c^k, these will typically be distinct as long as c^k≪ N. Conversely, different sub-trees rooted in the neighbourhood of j can be connected to a common site, hence producing a loop, if c^k = 𝒪 (N), or k = 𝒪 (ln (N)). In the large N limit these graphs therefore become locally treelike and the factorisation (<ref>) will again become exact: conditional on a given node, the branches rooted at that node become independent of each other. Equation (<ref>) then simplifies to
 
    P(y_j)=e^-i/2λ_εy_j^2/r_j∏_k∈∂ j∫dy_k   e^-i/2(y_j-y_k)^2P^(j)(y_k)

Similarly one can show for the marginals of the cavity distribution around node j

    P^(j)(y_k)=e^-i/2λ_εy_k^2/r_k∏_l∈∂ k∖ j∫dy_l   e^-i/2(y_l-y_k)^2P^(k)(y_l)

where ∂ k∖ j indicates the neighbourhood of node k excluding node j. As all distributions involved are zero mean Gaussians, also the marginals must be of this form, i.e. 

    P^(j)(y_k)=√(ω_k^(j)/2π)e^-1/2ω_k^(j)y_k^2    
    P(y_j)=√(ω_j/2π)e^-1/2ω_j y_j^2

We follow statistical terminology and call the ω, which are inverse variances, precisions <cit.>. Their real part must be positive in order to preserve normalisability of the corresponding Gaussians. 

In terms of the precisions and using r_j^-1 = τ_j c,  equations (<ref>), (<ref>) become 

    ω^(j)_k=iλ_ετ_k c + ∑_l∈∂ k ∖ jiω_l^(k)/i+ω_l^(k)    ω_j=iλ_ετ_j c + ∑_k∈∂ jiω_k^(j)/i+ω_k^(j)

as derived in more detail in appendix <ref>. 
The equations for the cavity precisions form a closed set {ω_k^(j)} that can be solved iteratively.
Note that for an actual tree, no iteration is required as the equations can be solved recursively by working inwards from the leaves. Once the cavity precisions are known, the marginal precisions {ω_j} can be deduced.
Finally, from (<ref>) and (<ref>) one obtains the diagonal entries of the resolvent: 

    G_jj=iτ_j c/ω_j

The factor τ_j c = r_j^-1 arises here from the transformation from y_j back to x_j.

Given a specific realisation of the disorder, i.e. for a single instance of the matrix 𝐌, we know the rates {r_i} and the connections {c_ij}, therefore (<ref>) can be solved iteratively starting from a suitable initial condition. The eigenvalue spectrum of the system is finally given by (<ref>) and (<ref>). We will refer to this procedure as the single instance cavity method.

In the large N limit (<ref>-left) turns into a self-consistent equation for the distribution p (ω)
of the cavity precisions – here to keep the notation clean we drop the superscript indicating the cavity graph. For the simplest case of a random regular graph this reads 

    p(ω)=∫dτ ρ_τ(τ)∏_l=1^c-1dω_l p(ω_l) δ(ω - Ω_c-1)

where 

    Ω_a = Ω_a({ω_l},τ)=iλ_ετ c + ∑_l=1^aiω_l/i+ω_l

The intuition here is that because p(ω) is a distribution resulting from the solution of the equations (<ref>) for the cavity precisions, updating the precision on a randomly chosen edge of the graph with the r.h.s. of (<ref>) does not change the distribution.
Technically, one  assumes here that the distribution of G_jj, and consequently p(ω), is  self-averaging in the limit N→∞. In the case of a general random graph, the only change is an additional average over the number of neighbours k of a randomly chosen edge, with the appropriate probability weight kp_k/c; see appendix <ref> for details.

 A numerical solution for p(ω) at any given λ_ε can be obtained using a population dynamics algorithm <cit.>. The basic idea is to represent the distribution p(ω) with a population of N_p cavity precisions 𝒫=(ω_1,…,ω_N_p). One starts with a certain initial condition and lets the population evolve according to the update rule given by the delta function in (<ref>). Once equilibrated, the histogram of 𝒫 should give an approximation of p(ω). In summary the algorithm works as explained in the following box:



 

Population Dynamics Algorithm




  * Start with an initial (complex) population 𝒫=(ω_1,…,ω_N_p).

  * Pick c-1 random elements {ω_l} from 𝒫 and a sample τ from ρ_τ.

  * Replace a random element of the population with Ω_c-1({ω_l},τ).

  * Repeat 2 and 3 until equilibration is reached.

 




Finally, we use (<ref>) and (<ref>) to write the DOS as an average over the distributions p(ω) and ρ_τ

    ρ(λ) = lim_ε→ 0lim_N→∞1/π N∑_i=1^NIm G_ii(λ_ε)=lim_ε→ 01/πRe⟨τ c/Ω_c({ω_l},τ)⟩_{ω_l},τ

where the {ω_l} are sampled from the population of cavity precisions converged to equilibrium.

We next discuss the relative merits of single instance cavity method versus population dynamics, and the influence of ε.
The single instance method allows us to find the spectrum of (large) sparse symmetric matrices 𝐌, under the cavity approximation of factorisation in each cavity graph. In terms of computational cost this method is in principle much faster than direct diagonalisation because one only has to find the O(N) cavity precisions, typically from a number of iterations of the cavity equations that does not grow with N.  However, one still has to store all the information on the disorder {τ_i, c_ij} and, as is true generally with the cavity technique, one obtains little information about the eigenstates. The calculation also has to be repeated across a suitably fine grid of λ-values in order to find the spectrum.

In choosing the λ-grid, one has to bear in mind that the general approach replaces the N delta-functions in (<ref>) by Lorentzians of width ε, therefore ε is the “resolution" that we have on the lambda axis. To catch all eigenvalues of a single instance, one therefore requires a grid spacing in λ of order ε or smaller. Conversely, for a fixed λ-grid, ε has to be chosen larger than the grid spacing, otherwise the chance of hitting all eigenvalues becomes too low to obtain accurate results.

In practice, we always perform the cavity iterations themselves with ε=ε_0→ 0 (specifically we set ε_0∼10^-300) so that the resulting cavity precisions are not affected by the width of the Lorentzians. The required nonzero ε (≫ε_0) is then applied only in the evaluation of the average (<ref>), i.e. in the measurement step. 
This makes it easy to explore the effect of changes in ε, without having to solve the cavity equations afresh.
From (<ref>, <ref>) one sees that using ε≠ε_0 to calculate the marginal precisions is equivalent to adding ε-ε_0 to the inverse variance of each x_j. This provides a regularization for the case where the variance calculated using ε_0 is close to imaginary because the chosen λ has hit an eigenvalue.
 
In contrast to the single instance approach, the population dynamics algorithm is designed to give the DOS of infinitely large systems. There is no need to keep track of the disorder because of self-averaging, and we only have to let the population equilibrate. The eigenvalue spectrum becomes densely populated, typically showing a continuous part referred to as the bulk. This means that we are always able to compute the DOS over this region, even with ε_0∼ 10^-300. A common feature of (sparse) random matrices is that the states covering the bulk are in fact delocalised (or extended), and localisation (Lifshitz) tails are present at the edges of the spectrum <cit.>. The values of λ where these localisation transitions occur are called mobility edges. Pure points, i.e. isolated eigenvalues <cit.>, do sometimes occur within the bulk of the spectrum, as is the case for e.g. sparse adjacency matrices with varying node degrees <cit.>. In the following sections we refer to the density of all the states of the system as the total DOS (tDOS), obtained by the population dynamics algorithm with ε small but finite, and to the density of the extended states only as the extended DOS (eDOS), obtained with ε effectively equal to zero. 



 §.§ Total DOS via population dynamics


We next present the results for the total DOS of a trap model on a random regular graph with connectivity c=5. Figure <ref>-left shows the results obtained using the population dynamics algorithm compared with data from direct diagonalisation of the master operator for finite N (labelled “numerics"). The agreement across the entire λ range is clearly very good. In figure <ref>-right we include the MF and RW-limits of the DOS for comparison; recall that the quantity r=-λ represents the relaxation rate of the system, so the plots showing -λρ(λ) vs -λ can equivalently be read as the density of ln (r), plotted against r on a logarithmic x-axis. We observe that the small |λ| tails (the slow modes governing the long-time dynamics) follow the MF trend (blue dashed lines), showing the same power law exponent asymptotically. Conversely, fast modes (large |λ|) show a non-linear DOS which originates primarily from the Kesten-McKay law (RW limit). We note here that because of the MF tail, one expects systems of finite size to have a spectral gap that scales with N as in mean field <cit.>, so the second largest eigenvalue should be bounded from above by -τ_max^-1∼-N^-β; a detailed analysis of the N-scaling of the spectral gap, however, is beyond the scope of the present work. Note that at the highest temperature T=2.5, the small |λ| tail of the DOS shows larger statistical uncertainties because of finite size effects: in direct diagonalisation, finite-sized matrices only rarely have eigenvalues in this region; similarly population dynamics sampling runs of finite length produce only a limited number of samples contributing to the slow mode regime. 


To understand the structure of the DOS in more qualitative terms, we can perform a simple (high T) analytical approximation: we take one cavity iteration at finite temperature starting from the infinite temperature solution. This means that only the local disorder is taken into account when computing the DOS, i.e. the central node receives its messages from c neighbours belonging to an infinite temperature cavity network. A similar idea, called the single defect approximation, has been used to explain localisation phenomena arising from topological disorder in random lattices <cit.>. In the large T and N limits, where all nodes become equivalent, the cavity precision distribution p(ω) becomes delta-peaked on the value ω̅ that  solves (<ref>), i.e.

    ω̅=Ω_c-1({ω̅},1)

The approximated total DOS is then evaluated as in (<ref>)

    ρ^A(λ)=lim_ε→ 01/πRe⟨τ c/Ω_c({ω̅},τ)⟩_τ

The average can be performed analytically, as detailed in appendix <ref>. Figure <ref>-left shows the resulting first order approximation ρ^A(λ) against the DOS obtained by direct diagonalisation. One observes that the approximation is in remarkably good agreement with the numerical data in the region of slow (MF-like) modes, though even in the RW-like regime it is qualitatively correct. We can iterate the scheme to obtain higher order approximations: to the second order we perform two cavity iterations at finite temperature starting from the infinite temperature solution, and so on. The second order approximation is then given by

    ρ^2A(λ)=lim_ε→ 01/πRe⟨τ c/Ω_c({Ω_c-1({ω̅},τ_l)},τ)⟩_{τ_l},τ

This average cannot be carried out analytically but is straightforward to  perform by sampling from the distribution of waiting times. Figure <ref>-right shows the first and second order approximations on a linear scale. One gets a slightly better result with the second order approximation ρ^2A in the large |λ|-region, though not yet a quantitative  match to the full cavity predictions. In the small |λ| tail we find (not shown here) that there is no significant difference between the first and second order approximations. As a final remark, we note that an infinite order approximation would give the population dynamics result: in this case, the infinite temperature solution ω̅ corresponds to a particular initial condition for 𝒫, which is lost after a large number of iterations of the approximation scheme. 





 §.§ Extended DOS and IPR


As explained above, we can measure the extended DOS (only) by evaluating the cavity predictions in the limit ε→ 0 as it is this part of the spectrum that becomes continuous in the thermodynamic limit. Comparing the extended and total DOS then allows us to locate the mobility edges of the system. Figure <ref>-left shows the total DOS and the extended DOS on a linear scale, with an inset zooming in on the localisation transition occurring on the right end of the spectrum. Figure <ref>-right displays the same plot with a logarithmic y-axis, where we have included evaluations of the total DOS for different ε values. This allows one to estimate the left end of the spectrum from the point on the λ-axis where the total DOS ceases to be ε-independent. Note that on approaching the mobility edges, the convergence of the population dynamics to its steady state becomes very slow. The  peaks in the extended DOS that are visible in the inset of figure <ref>-left are caused by this and should accordingly be ignored as unphysical. While it is not surprising to find localisation tails at the edges of the spectrum, at least from a random matrix perspective, it is remarkable that the appearance of mobility edges arises directly from the combination of two limiting cases with exclusively extended (RW) and localised (MF) eigenvectors, respectively. We can already argue that the fastest and slowest processes are governed by localised modes, with an intermediate regime where all modes are delocalised. Since we are interested in the long time dynamics, our attention will be focused on the bulk of extended states and on the slow (MF) localised modes only; we will also show that the fraction of fast localised states is relatively small compared to that of slow modes. As we will see, the mobility edge occurring on the slow end of the spectrum (of eigenvalues, or similarly relaxation rates) allows one to identify three different regimes in the time domain.


The localisation transition described above is associated with a change in the distribution p(ω): the population of cavity precisions converges to a steady state which has complex support for λ in the bulk of the spectrum, and purely imaginary support outside. 
This transition can be detected by considering the average real part of the cavity precisions, which is shown in figure <ref>-left: 
as these precisions must have non-negative real part, the vanishing of the average real part means all real parts are zero. “Zero” is to be interpreted here as of order ε_0, the value of ε used in the population dynamics; our ε_0  ∼ 10^-300 is indistinguishable from zero even on the logarithmic scale of figure <ref>-left. The figure shows the average real part of the cavity precisions and the total/extended DOS as a function of λ. As claimed above, the average real part is nonzero in the bulk of the spectrum, goes to zero exactly where the extended DOS does, and then vanishes within the localised spectrum. The approach of the average real part to zero is continuous (see inset), indicating that the transition in the structure of the distribution of cavity precisions is likewise continuous.

Note that having imaginary cavity precisions amounts to having real diagonal entries of the resolvent (via eq. (<ref>)), which in the context of field theory and Anderson localisation are related to the so called self-energies. Similarly to what we have outlined above,
the state of an electron in a disordered medium is classified as localised or delocalised depending on whether the electron's self-energy is real or complex <cit.>.

We complement the above results by measuring the average degree of localisation of the eigenvectors.
For N→∞ one cannot access the IPR of individual eigenvectors. Instead one can consider the average IPR in a small range ε around λ and then take ε to zero: 

    I̅_̅2̅(λ)=lim_ε→ 0lim_N→∞1/Nρ(λ)∑_α=0^N-1δ_ε(λ-λ_α)I_2(𝐯_α)

where δ_ε(x)=ε/[π(x^2+ε^2)] is a Lorentzian of width ε and ρ(λ) is assumed to be calculated similarly, using δ_ε(x) instead of δ(x) in the definition (<ref>). The order of the limits in the definition ensures self-averaging because the number of λ_α that contribute, which is of order Nε, becomes large.

Swapping the two limits, i.e. assuming that at any given λ at most one eigenvector contributes to the average IPR, one can relate I̅_2(λ) to the squared modulus |G_jj(λ)|^2 of the resolvent entries. Bollé et al. obtained from this a formula that allows the IPR to be evaluated within population dynamics <cit.>, and used this to study localisation transitions in Laplacian and Levy matrices.  In our notation their expression reads

    I̅_2(λ)=lim_ε→ 0ε/πρ(λ)⟨| τ c/Ω_c({ω_l},τ)|^2 ⟩_{ω_l},τ

which can be rewritten more explicitly as

    I̅_2(λ)=lim_ε→ 0⟨ε/(ε+A_r)^2+(λ+A_i)^2⟩/⟨ε+A_r/(ε+A_r)^2+(λ+A_i)^2⟩

where, to keep the notation simple, we have used

    A = A_r +iA_i = Ω_c/τ c - iλ_ε = 1/τ c∑_l=1^ciω_l/i+ω_l

with A_r and A_i respectively the real and imaginary part of A. If the cavity precisions have zero/positive real part, then A_r is zero/positive accordingly. It follows from (<ref>) that I̅_2(λ)=1 in the localised part of the spectrum, where the distribution p(ω) has purely imaginary support, and it is of order ε within the bulk, where the support of p(ω) is complex.
While we expect an average IPR of order unity, a value exactly equal to one is implausible in our case. This can be seen from the large c-limit, where we must recover the IPR of the MF eigenvectors (<ref>), for which clearly I̅_2 < 1 (see also figure <ref>-right). The discrepancy indicates that the swapping of the limits ε→ 0  and N→∞ is not in general justified. 
Nonetheless, (<ref>) remains useful as a tool for differentiating between localised and extended parts of a spectrum.

As an alternative to the treatment of Bollé et al, we suggest an approximation to the IPR that is derived by taking N→∞ at fixed ε, and therefore is suitable for use within population dynamics based on (<ref>). We leave the derivation to appendix <ref> and only give the result

    I̅_2^⋆(λ) = lim_ε→ 02ε/πρ(λ)Var[Re(τ c/Ω_c({ω_l},τ))]_{ω_l},τ

where Var(·) indicates the variance. The order of limits (N→∞ first, then ε→ 0) used ensures that there are always enough states within the λ-range of width ε where quantities are measured. 
In this regard, our approach is opposite to that of Bollé et al, where the two limits are inverted. Even so, we observe a very close agreement between the IPR estimates (<ref>) and (<ref>), as shown in figure <ref> – appendix <ref>.

Figure <ref>-right shows the total DOS, the extended DOS, and the average IPR predicted by (<ref>). As explained above, in the localised region I̅_2(λ) has a constant value of one, even where the total DOS drops to the ε-value used in the measurement step of the population dynamics algorithm, i.e. outside of the support of the spectrum. The localisation transitions, detected as the points on the λ-axis where the value of I̅_2 changes from 𝒪(ε) to 𝒪(1), occur where the extended DOS drops to 𝒪(ε_0). This confirms what we have discussed before: in the thermodynamic limit the spectrum has a continuous part of extended states, with localisation tails of pure point states occurring at the edges of the bulk.





 §.§ Finite size effects


The advantage of the population dynamics approach is that it allows us to evaluate the spectral properties of infinitely large systems at a relatively low computational cost. As described above, the distribution p(ω) is approximated by a large population 𝒫 of representative cavity precisions samples. This population converges to steady states that depend on the value of λ, and it undergoes critical transitions at the mobility edges. The location of these transitions turns out to have a non-negligible dependence on the population size N_p (see figure <ref>-left). Finite size effects in population dynamics algorithms have been discussed in the context of the Moran model <cit.> and in the evaluation of large deviation functions <cit.>. In particular in <cit.> the authors show that in their case, systematic errors in the algorithm decrease proportionally to the inverse of the population size.
In order to determine the actual position of the mobility edges occurring within our spectra we assume that

    |λ^L/R_∞-λ^L/R_c(N_p)| ∼ N_p^-a

where λ^L/R_c(N_p) is the left/right mobility edge measured using a population of size N_p, and λ^L/R_∞=lim_N_p→∞λ^L/R_c(N_p). Accordingly, we gather data for different N_p and  then fit them using 

    λ^L/R_c(N_p) = c_1 + c_2 N_p^-a

where c_1, c_2 and the exponent a are determined by minimizing the least-squares deviation. This then identifies, in particular, the extrapolated value λ^L/R_∞ =c_1.
As shown in figure <ref>-right, with the exponent a chosen in this way our data {(λ^L/R_c(N_p),N_p^-a)} do lie on a straight line to a good approximation as (<ref>) assumes. This fitting method is applied to determine the right and left mobility edges for different values of the temperature, with mean connectivity c=5. 
The exponent a is non-trivial in our case: it shows a monotonic increase with inverse temperature and typically lies between 0 and 1 (see figure <ref>, top-right), in contrast to the setting in <cit.> where a=1. 
We conjecture that the T-dependence of a is related to the fact that also the exponent of the distribution of waiting times ρ_τ(τ) varies with T; a more precise quantitative understanding of the value of a remains an open problem, however.
Figure <ref>-left shows the DOS with the extrapolated right λ_∞^R  and left λ_∞^L mobility edges (respectively on the left and right of the plot, because the x-axis shows -λ), the RW-DOS and the power-law MF-DOS for the slow modes. We note that λ_∞^R lies at a point on the λ-axis where the full DOS is already MF-like. In facts, it seems natural to describe the spectrum as composed of three main regions: the slowest modes possess MF-like features as they are localised and power-law distributed. The fastest modes are delocalised and exhibit a non-monotonic DOS that is closely related to  the Kesten-McKay law for the RW limit. Finally, the intermediate region (green shaded area) has mixed properties of the two limiting cases: here the eigenstates are delocalised (RW-like) but show a power-law distribution (MF-like). We also observe that this intermediate region becomes wider – the fraction of delocalised modes with a MF-like density of states increases – as the temperature decreases (figure <ref>-right).
 




From the extrapolated position of the mobility edges in the eigenvalue spectrum we can estimate the fraction of localised modes π_tot in our system. This is given by the integral of the total DOS over the λ-regions containing localised states, which in our case lie at the edges of the spectrum. It is in fact simpler to evaluate π_tot by working out the complement, i.e. integrating over the bulk of the DOS:

    π_tot=1-∫_λ_∞^L^λ_∞^Rdλ ρ(λ)

Here λ_∞^L and λ_∞^R are the left and right mobility edges of the system, extrapolated to infinite population size as explained above. 

We can similarly obtain the fraction of localised fast/slow modes by integrating over the λ-region at the left/right end of the spectrum. Since we have Lorentzian tails of width ε affecting the total DOS, the most accurate way of computing these fractions is to locate the left end of the spectrum λ_end^L by exploiting the ε dependence of the total DOS (see figure <ref>-right), then evaluating the fraction of localised fast modes as 

    π_L = ∫_λ_end^L^λ_∞^Ldλ ρ(λ)

The fraction of localised slow modes is finally given by π_R = π_tot - π_L. Figure <ref> shows the fractions of localised modes (left) and the right mobility edge (bottom-right) as functions of the temperature, for mean connectivity c=5. We observe that as the temperature decreases the localisation region on the right edge of the spectrum becomes narrower while the fraction of slow localised modes in this increases. Overall, the total fraction of localised eigenstates becomes larger as the temperature decreases. Nevertheless, even at T=0.5 – the lowest temperature considered here – the fraction of localised modes only amounts to around 10% of the total DOS. The majority of these localised eigenstates lie in the low |λ| tail, as can be seen from the fact that the quantities π_R and π_tot are almost overlapping on the log scale shown in figure <ref>-left. Importantly for the long-time dynamics, all the slowest modes in the system are localised, at least for the temperature regime that we have considered here.
The temperature trend for low T is consistent, at the other end, with the T→∞ limit: here we obtain a RW spectrum with only extended and no localised modes. 





§ MORE DISORDERED NETWORK TOPOLOGIES


So far we have focused on the case of random regular graph (RRG) connectivity, where the network defining the possible paths among minima in the potential energy landscape has a regular structure that becomes free of disorder in the thermodynamic limit. From a topological perspective the absence of disorder might seem as unrealistic as, say, in the n-dimensional hypercubic lattice or the complete graph with its mean field connectivity. However, the RRG does introduce the essential features of sparse random networks, i.e. it is “infinite dimensional” – the number of nodes grows exponentially with distance – and it confines all dynamical transitions to a local environment. The RRG case is also interesting as the localisation properties of the eigenvectors of the master operator are the opposite of those in the mean-field Bouchaud model, where the connectivity is infinite and all eigenmodes are power-law localised on the energy axis. 

The question that we want to address in this section is whether the RRG case possesses all the relevant features of sparsely connected energy landscapes, at least in terms of the spectral properties discussed so far, or whether more disordered network topologies add new features (see e.g. <cit.>). We therefore extend our analysis to  Erdös-Rényi (ER) and scale-free (SF) graph structures, which have been widely studied in other contexts <cit.>. 
They both have finite average degree c but are paradigmatic as graph ensembles with finite (ER) and infinite (SF) degree variance, respectively.
Further motivation comes from the fact that numerical studies on a relatively small number of Lennard-Jones interacting atoms have suggested a configuration-space connectivity of the scale-free type <cit.>, which is also the network topology assumed by Baronchelli et al <cit.>
as discussed in the introduction. The SF case may therefore represent the best candidate for modelling configuration space connectivity, though we stress that our approach is flexible and can be applied to any network topology without short loops. 

Looking at the random walk (RW) and mean-field (MF) limits, we note first that in the former case there is no simple closed form expression for the DOS, analogous to (<ref>) for regular graphs, on complex network structures: we will have to obtain results by population dynamics instead for the T→∞ limit.
The limit c→∞, on the other hand, effectively brings us back to the fully connected case so our previous results and discussion for the MF limit still apply.

Erdös-Rényi graphs <cit.> of size N are constructed by assigning an edge between any pair of vertices with probability p, so the average number of edges in the network is N(N-1)p/2 and we need (N-1)p to be finite to ensure that the resulting graphs are sparse. The probability that a given node has k neighbours then follows a binomial distribution, which in the large N limit approaches a Poisson distribution with parameter c=(N-1)p. We therefore apply our cavity method assuming p_k = e^-c c^k/k!, with ⟨ k ⟩ = Var(k)=c. Since the Poisson distribution is strongly peaked around c, the local environment of these graphs is typically subject to weak fluctuations, and the overall structure is not far from that of random regular graphs.
In the following we will assume c=5, which ensures that the fraction of nodes in the giant cluster is approximately equal to one, ignoring the effects of very small disconnected components on the spectral properties of the whole system. The population dynamics algorithm applies as explained in section <ref>, with the only difference that at each update we pick k-1 elements from the population of cavity marginals with probability p_k k/c; see appendix <ref> for further details.

The spectral features of the ER ensemble with c=5 and T=1.5 are shown in figure <ref>-left. The total DOS is displayed for evaluations involving two different values of ε, whose effect is visible on the left of the plot. The total DOS of the RW limit would have the same ε tail for small |λ| but we do not show this region as the RW DOS becomes too small to estimate reliably there. Similarly to the case of random regular graphs, the DOS is composed of three main parts: a mean field power-law tail occurs at the slow end of the spectrum, covering localised (left) and delocalised (centre) modes, while the distribution of fast modes (right) is non-monotonic and follows closely the DOS of the corresponding RW limit. In contrast to the RRG case, the connectivity disorder alone is enough to induce localisation transitions within the spectrum; the mobility edges are extrapolated by the least squares fit discussed in the previous section, and they are marked by the green (T=1.5) and red (T→∞) dashed lines in the plot. We observe that the area under the total DOS of MF localised modes at T=1.5 is much larger than that corresponding to the RW case. This means that the intrinsic localisation attributes of ER graphs are sub-dominant with respect to the effects introduced by energy disorder, which become stronger when the temperature is lowered.



The last class of networks that we address in this work is the scale-free type. In the form originally proposed <cit.>, these networks are constructed via preferential attachment: starting with a dimer of two nodes linked together, one connects a new node to the existing ones with a probability that is proportional to the number of links that they already have, repeating the process until the network has the desired size. More generally SF networks are characterised by a degree distribution p_k that decays as k^-γ. Here γ is typically in the range 2<γ<3, implying that second and higher order moments diverge. This motivates the appellative “scale-free" as, by contrast to Erdös-Rényi and random regular graphs, the degree fluctuations are infinitely large and have no intrinsic scale. Defining a SF graph ensemble by assigning equal probability to all networks with the given degree distribution, one typically finds many hubs within the network, i.e. nodes with very high degree, occurring at all (degree) scales. In order to minimise the number of disconnected sub-graphs we introduce a lower bound on the range of degrees by imposing p_0 = p_1 = 0. Also, in practice we cannot deal numerically with unbounded probability distributions, and a cutoff K_MAX has to be specified so that p_k = 0   ∀   k> K_MAX; we take K_MAX=1000. The cutoff is entirely immaterial for the ER case, where for the Poisson degree distribution with e.g. c=5 one has p_K_MAX∼ 5^1000/1000! ∼ 10^-1869. Even for SF graphs, with γ=2.5, p_K_MAX∼ 1000^-2.5∼ 10^-8 so the cutoff lies far in the tail of the distribution. It does make all moments of the distribution finite, but still retains much larger degree fluctuations than for random regular and ER graphs.

The results for the SF ensemble with γ=2.5 and T=1.5 are shown in figure <ref>-right, using the same colour scheme as for the ER plot on the left and a population of 2500 cavity precisions. A striking difference is that the spectrum is much broader than for the previous cases, by at least one order of magnitude (scaling the rates by the average connectivity is ineffective when the variance of degrees is large as here): a long tail of fast, localised modes appears (at the left end of the spectrum, which on the plot is on the right as the x-axis shows -λ). Similarly to the case of ER graphs, the area under the total DOS of MF localised modes at finite temperature (T=1.5) is much larger than that corresponding to the RW case, which demonstrates the strengthening of slow mode localisation when the temperature is lowered. Overall, in spite of some differences in the details, the DOS for SF networks has the same structure as for random regular and Erdös-Rényi graphs, with a tail of slow modes following the mean field statistics, a mixed region where the DOS is MF-like but eigenstates are delocalised, and a remaining part of the spectrum that is non-monotonic and closely related to the associated RW case. By lowering the temperature we induce a shift in the DOS towards slower modes, the range of MF-RW mixed modes becomes wider, and the fraction of slow localised modes increases.



§ CONCLUSIONS AND FUTURE PERSPECTIVES


In this paper we have considered the problem of walks on the potential energy landscape as described by the trap model of Bouchaud and others, extending previous analyses to the case of sparse inter-trap connectivity. In this scenario there are two different sources of disorder: one is associated with the topology defining the connectivity among minima, and the other one is given by the different energy depth of the traps. Accordingly there are two important notions of distance: the distance on the graph structure and the distance on the energy axis. The sparse structure of the master operator 𝐌 makes the problem impossible to solve with analytical tools, and we then approached it by means of the cavity method, which in the thermodynamic limit leads to a population dynamics algorithm. This allowed us to evaluate the eigenvalue spectrum of 𝐌, and the localisation properties of the associated eigenstates (the modes of the dynamics), which are key to understanding the dynamical behaviour of the model.

We first discussed the spectral properties of the ground state, i.e. the equilibrium distribution, focussing on how the IPR scales with system size for different temperatures; these results are independent of network structure because the transition rates obey detailed balance. 
In the bulk of the paper we considered the case of random networks with regular connectivity, where the key system parameters are the temperature T and the mean connectivity c. We discussed the limiting situation of infinite temperature, where the dynamics is a random walk (RW) and only the distance on the graph structure matters. Here the density of states (DOS) is given by a shifted and scaled Kesten-McKay law, and all eigenmodes are delocalised. In the opposite mean-field (MF) limit c→∞, where traps are distinguished only by their energy depth, the eigenstates are localised and the DOS 
follows a power law with a T-dependent exponent. We found that these features are combined in the general case of finite connectivity and finite temperature where both notions of distance are relevant: a MF-like tail of slow localised modes (governing the long time dynamics) appears, while fast modes follow the RW case, being delocalised and showing a non-monotonic DOS related to the Kesten-McKay law.
Localisation transitions appear within the spectrum at the changeover between these two behaviours; they correspond to continuous transitions in the nature of the support of p(ω), the distribution of cavity precisions that is the key quantity within the population dynamics algorithm. The location of these transitions is affected by population size, and we extrapolated them to the infinite population limit using a simple power law form.
This revealed a surprise: the combination of RW and MF features give rise to a mixed region separating the fast modes from the slow modes, where the DOS has a MF shape but eigenstates are nonetheless delocalised. 

The shape of the DOS, and particularly the power law MF-like tail of slow modes, are well captured by a simple “high temperature" approximation scheme. At first order, one cavity iteration (involving disorder) is performed starting from the infinite temperature solution found on the cavity graph (where, in the case of RRG, disorder is absent). This is similar to the “single defect approximation" <cit.>: the central vertex is the only source of randomness and this allows one to find an analytical solution for the DOS.

We observed the same overall structure in the spectra for more disordered network topologies, specifically Erdös-Rényi and scale-free graphs, though with some changes in the details particularly for the fastest modes. The broader degree distributions of these networks are sufficient to induce localisation transitions even for infinite temperature, i.e. without any effects from the trap depths.
However, the fraction of slow localised modes is greatly enhanced at finite temperature, where the corresponding DOS has a power-law MF shape. The latter feature arises in every graph ensemble considered here. Bearing in mind that the spectrum of relaxation times is simply given by the collection of inverse eigenvalues {1/|λ_α|}, the ultimate long time dynamics should then always be of mean field kind. This asymptotic independence from the network structure is in fact consistent with the way the trap model was originally designed, 
namely to describe dynamics in configuration-space on a very long time scale where deep minima are effectively fully connected by paths passing through shallow traps.
The presence of a region in the spectrum with mixed MF and RW properties suggests the existence of an intermediate time scale during which the dynamics should exhibit features of both the short-time random walk behaviour, which is network dependent, and the long time mean field evolution. How this distinction based on the spectral analysis can be quantified in the time domain remains an open problem and will have to be addressed in future. 

More generally, the lack of detailed information about the eigenvectors remains the major limitation of our approach, as it impedes a direct evaluation and classification of the ageing dynamics. Nevertheless, we believe that there is scope here for significant improvement and we hope that this work stimulates further investigations in this direction. In particular, these should include exploring the time domain and trying to characterise the three different regimes that we have highlighted above. This can be done e.g. by looking at the time-dependent probability of return to a given trap, which can be expressed as the Laplace transform of the local density of states (i.e. the contribution to the total DOS from a local node with a given trap depth). The local DOS will typically be peaked around the value of -λ that corresponds to the initial decay of the return probability, and this would allow one to focus on a single one of the three distinct regions that we have identified in the spectrum. Alongside the return probability, other time dependent quantities, such as the mean number of distinct nodes visited within some time t, will help to characterise the (non-equilibrium) dynamics and identify the system's time scales. This question can be tackled with the techniques of <cit.>, accompanied by numerical simulations of random walks on trapping networks to assess finite size and pre-asymptotic (short t) effects.

Future work should aim also to elucidate the link between the dynamics and the localisation properties of glassy dynamics on networks. Insights might come from a closer look at the structure of the eigenvectors. The IPR carries no information on the spatial distribution of the eigenvector components, nor is it able to distinguish between exponential or non-exponential localisation, either in energy or on the graph structure. For this reason it is not clear from our results how the two sources of randomness influence the localisation strength of, say, the localised eigenstates governing the long time dynamics. The region of the spectrum that we have characterised as “MF-localised” might exhibit eigenvectors that are e.g. exponentially localised within small areas of the network, rather than covering all the nodes through a power law decay with difference in trap energy (as happens for mean field connectivity). Then, for waiting times in the range of the slow localised modes (which will be the typical case at low temperature), when the system escapes from a deep trap all delocalised modes will have decayed and the motion will remain confined to the neighbourhood of the initial node, in contrast to MF dynamics. Similarly, in the “mixed” region of the spectrum the delocalised eigenvectors might have non-trivial spatial structure; one might conjecture that they should be concentrated onto an extensive number of clusters of nodes, interpolating between the localised slow modes and the delocalised fast ones. Overall, therefore, more detailed information on the distribution of the entries of the slow eigenvectors and their spatial correlations should allow a better understanding of the asymptotic dynamics. The former quantity can easily be obtained numerically via e.g. the multifractality spectrum <cit.>, while the latter could be assessed with an analogue of the radial distribution function (or similar measures) from liquid state theory <cit.>. To calculate them in the large system size limit, on the other hand, as outputs from a population dynamics algorithm, remains a technical challenge.

Finally, we point out two possible extensions of the present work: given the evidence for correlations between the depth and the number of neighbours of an energy minimum, this would be a feature worth including within our model. Such a more general setting should be amenable to an analysis similar to the one in this paper as we sketch briefly in appendix <ref>. In principle, correlations in the degree of neighbouring minima could also be taken into account. The idea, supported by previous works on simulations of L-J interacting atoms <cit.>, is that deep minima in the potential energy landscape are surrounded by many shallower ones, creating a hierarchical structure and thus inducing degree-degree and degree-energy correlations.

The second interesting model extension would be to consider alternative transition mechanisms among minima, as considered e.g. by Barrat and Mézard <cit.>, who used Glauber transition rates. These rates depend on the difference in energy between the departure and arrival nodes, and as energy-decreasing transitions are always allowed one can picture the situation on a fully-connected graph as an energy landscape made up of steps rather than traps. Here the entropy of relaxation paths is key and leads to a completely different phenomenology (see for example <cit.>), with ageing arising from entropic rather than energetic barriers. On sparse graphs, on the other hand, even Glauber dynamics will encounter energy barriers – consider a deep trap with only shallow neighbours – and so a much richer and possibly more realistic dynamics should result. Work towards analysing this case is in progress. 





§ 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). RGM gratefully acknowledges insightful discussions with Chiara Cammarota, Davide Facoetti and Aldo Glielmo.









unsrt 
 













§ LOCALISATION IN MEAN FIELD LIMIT


We give here a qualitative argument why in the mean-field limit a generic eigenvector associated with eigenvalue λ<0 will be localised. We regard λ<0 as fixed here and take N→∞, to stay well away from the groundstate. For finite N one would expect a crossover to the localisation properties of the (delocalised, for T>1) ground state as λ→ 0.

The explicit form of the eigenvector components is given in (<ref>). To estimate how these components vary across nodes i, consider a typical realisation of the trap depths {E_1, E_2, …, E_N}, arranged in ascending order such that E_i<E_i+1. The inverse trapping times τ_i^-1=exp(-β E_i), which determine the eigenvector components u_i ∝ (λ + τ_i^-1)^-1 (we drop the eigenvector label α here), are then in descending order. The largest component will occur at the node i with τ_i^-1 closest to |λ|; call this node j. The number of inverse trapping times at other nodes that lie in an interval [τ_j^-1,τ_i^-1] is typically Nρ_τ^-1(τ_j^-1)
(τ_i^-1 - τ_j^-1) where ρ_τ^-1 denotes the distribution of inverse trapping times. Abbreviating this density of states-factor as simply ρ, we can therefore write

    j-i=Nρ(τ_i^-1-τ_j^-1)

as a deterministic approximation for the values of the inverse trapping times around τ_j^-1.

Now call S_N and S_D respectively the sum on the numerator and denominator in the definition (<ref>) of I_2. Using that λ≈ -τ_j^-1, we have

    S_N=∑_i u_i^4∝∑_i (λ+τ_i^-1)^-4≃∑_i (j-i/N ρ)^-4

The last sum can be approximated as twice the integral over the positive values of m=i-j 

    S_N≃ 2 ρ^4 ∫_1^∞dm 1/(m/N)^4∝ρ^4 N^4

Similarly, for the sum in the denominator in I_2 we get

    S_D≃( 2 ρ^2 ∫_1^∞dm 1/(m/N)^2)^2∝ρ^4 N^4

Taking the ratio, it follows that I_2 = O(1) for any eigenvector 𝐮 with an eigenvalue away from zero; this statement holds at any temperature. Note that the sums or integrals defining S_N and S_D all converge at the upper end, i.e. have their mass concentrated around small m=i-j. This justifies our initial approximation of focussing on inverse trapping times close to τ_j^-1. It also implies that the above argument for the IPR of the right eigenvectors of the master operator applies equally to the eigenvectors of the symmetric master operator: these differ only by factors of τ_i^-1/2, which vary weakly (by O(1/N)) across the relevant range where m=i-j is finite.




§ THE CAVITY METHOD

In this appendix we illustrate how to derive the equations (<ref>) relating the (cavity) precisions, starting from the cavity marginal probability distributions as expressed in (<ref>), i.e. 
    P^(j)(y_k)=e^-i/2λ_εy_k^2/r_j∏_l∈∂ k∖ j∫dy_l   e^-i/2(y_l-y_k)^2P^(k)(y_l)

which is based on the factorisation  (<ref>).
As depicted in figure <ref>, this assumption works well if the graph lacks short loops, i.e. when it is locally treelike. Under this condition the correlations between the variables 𝐲_∂ j belonging to different branches become negligible when the common root is removed from the graph. Inserting the ansatz given in (<ref>), i.e.

    P^(j)(y_k)=√(ω_k^(j)/2π)e^-1/2ω_k^(j)y_k^2

we obtain for the cavity marginals

    P^(j)(y_k)∝ e^-i/2λ_εy_k^2/r_k∏_l∈∂ k∖ j∫dy_l  e^-i/2(y_l-y_k)^2 - 1/2ω^(k)_l y_l^2

Completing the square and integrating out y_l one finds

    P^(j)(y_k)    ∝ e^-i/2λ_εy_k^2/r_k∏_l∈∂ k∖ j e^-i/2y_k^2 - 1/2 (i+ω^(k)_l)^-1y_k^2

As by definition this must be proportional to 
exp(-ω_k^(j)y_k^2/2), it follows that

    ω^(j)_k=iλ_ε/r_k + ∑_l∈∂ k ∖ jiω_l^(k)/i+ω_l^(k)

which after replacing r_k=(τ_k c)^-1 is  the desired cavity equation. The calculation for the marginal precisions {ω_j} (see (<ref>)) is exactly analogous. Note that the discussion so far allows any kind of graph structure, i.e. it is independent of a specific choice for the degree distribution p_k, and it also allows correlations between degree and energies.

When going from the above considerations for a single finite-sized graph to the thermodynamic limit, one assumes that the cavity precisions {ω_k^(j)} are random variables taken from some distribution p(ω). Equation (<ref>) then turns into a self-consistent equation for p(ω). For a general degree distribution p_k and a joint distribution ρ_τ,k(τ,k) = ρ_τ|k(τ|k)p_k this reads


    p(ω)=∑_k p_k k/c∫dτ ρ_τ|k(τ|k)∏_l=1^k-1dω_l p(ω_l) δ(ω - Ω_k-1)

where c is the average degree of the network, p_k k/c is the probability that a randomly chosen edge connects the root-node to a neighbour with degree k, and

    Ω_k-1 = Ω_k-1 ({ω_l},τ)= iλ_ετ c + ∑_l=1^k-1iω_l/i+ω_l

Clearly (<ref>) reduces to the result for random regular graphs (<ref>) in the main text once we impose that ρ_τ|k(τ|k)=ρ_τ(τ) and p_k=δ_c,k. The population dynamics algorithm for the general case (<ref>) follows the same protocol as discussed in section <ref>, with the only difference that, at each update, one has to pick k randomly with weight p_k k/c, then draw k-1 elements from 𝒫 and τ from ρ_τ|k(τ|k).

We conclude this appendix with a final remark: while the change of variable y_i=x_i r_i^1/2 is not essential for single instance cavity evaluations, i.e. for fixed realisations of the disorder, this step becomes necessary in going to the thermodynamic limit. This is because otherwise correlations between cavity precisions on different branches of a cavity graph would be created by the coupling to the local disorder, and therefore the assumption of statistical independence between these cavity precisions would be violated.







§ HIGH T APPROXIMATION


In this appendix we discuss the construction of the high T approximation for the DOS, ρ^A(λ) (see (<ref>)). As mentioned in section <ref>, we take one cavity iteration at finite T, starting from the infinite temperature solution. Consequently the cavity precisions are evaluated without on-site disorder as the limit T→∞ gives τ_k=1 for all k. Also, in the thermodynamic limit the random regular graph structure becomes effectively a regular tree, and the problem of finding the cavity precisions becomes free of disorder. Figure <ref>-left shows a schematic representation of this procedure. In this non-disordered framework the distribution p(ω) is then expected to be delta peaked on some value ω̅. We have

    p(ω) = δ(ω-ω̅)     =∫dτρ_τ(τ)∏_l=1^c-1dω_l  p(ω_l)δ(ω-(iλ_εc+∑_l=1^c-1iω_l/i+ω_l))
        = ∫∏_l=1^c-1dω_l  δ(ω_l-ω̅)δ(ω-(iλ_εc+∑_l=1^c-1iω_l/i+ω_l))
        = δ(ω-(iλ_εc+(c-1)iω̅/i+ω̅))

Thus ω̅ needs to satisfy

    ω̅=iλ_εc+(c-1)iω̅/i+ω̅

and out of the two roots we need to pick the one with Re ω̅≥ 0, to which we also simply refer as ω̅. We recall that the DOS in the thermodynamic limit is obtained by averaging over the cavity precisions and waiting time distributions according to

    ρ(λ) =lim_ε→ 01/πRe⟨τ c/Ω_c({ω_l},τ)⟩_{ω_l},τ

The first order of our approximation scheme thus gives

    ρ^A(λ) =lim_ε→ 01/πRe⟨τ c/Ω_c({ω̅},τ)⟩_τ

where

    ⟨τ c/Ω_c({ω̅},τ)⟩_τ    =∫_1^∞dτ ρ_τ(τ) τ/iλ_ετ +iω̅/i+ω̅
        = T/iλ_ε∫_1^∞dτ τ^-T/τ+C(λ_ε)

with C(λ_ε)=ω̅/(λ_ε(i+ω̅)). 
The integral over τ can be done directly giving a hypergeometric function of T and C. Explicitly, one obtains

    ρ^A(λ)=Re[_2F_1(1,T;1+T |-C(λ))/iπλ]


The second order approximation consists of two cavity steps at finite T – and similarly the n^th order approximation would have n cavity steps – starting from the infinite temperature solution (see figure <ref>-right). The cavity precisions are evaluated with their on-site disorder. As a result, the average giving the DOS contains the disorder of the neighbouring environment {τ_k}, plus the local disorder of the central node τ: 

    ρ^2A(λ)=lim_ε→ 01/πRe⟨τ c/Ω_c({Ω_c-1({ω̅},τ_k)},τ)⟩_{τ_k},τ

Similarly to (<ref>) we have

    ⟨τ c/Ω_c({Ω_c-1({ω̅},τ_k)},τ)⟩_{τ_k},τ    =∫dτ ρ_τ(τ) ∏_k=1^cdτ_k ρ_τ(τ_k) τ c/iλ_ετ c+∑_k=1^ciΩ_c-1({ω̅},τ_k)/i+Ω_c-1({ω̅},τ_k)

The results of numerical evaluation of ρ^2A(λ) are discussed in section <ref>.

Since the equation for ω̅ is T-independent, the localisation transitions detected by this approximation do not depend on temperature and they always lie at the ends of the RW-limit spectrum, specifically for the RRG case we have λ^L/R=-(2c-1)/c, -1/c (from (<ref>)). This is true for any finite iteration of the approximation scheme, as can be argued inductively: the m-th iteration cavity precisions will be imaginary – indicating a localised region of the spectrum – whenever they are evaluated using an imaginary (m-1)-th iteration solution, as long as the imaginary part in λ_ε used for these evaluations is kept small enough.







§ THE INVERSE PARTICIPATION RATIO I_2^⋆


In this appendix we explain how to derive the equation (<ref>), which constitutes an alternative to the formula that Bollé et al proposed for estimating the average IPR <cit.>. We start by expressing the diagonal resolvent entries in terms of the eigenvector components

    G_jj (λ_ε) = ∑_α=1^N v_α,j^2/λ-iε + λ_α

whose imaginary part reads

    Im G_jj (λ_ε) = ∑_α=1^N ε/(λ-λ_α)^2 + ε^2  v_α,j^2

In order to simplify the notation we will omit the lambda argument in G_jj(λ_ε) and ρ(λ) where necessary. We also take N as large but finite and assume that ρ(λ) is finite, too. This can be ensured by choosing ε small but such that Nε≫ 1: for any given λ many eigenvalues then contribute to ρ(λ), which 
evaluates the DOS using Lorentzians δ_ε(λ - λ_α) of width ε. Since the statistics of the Im  G_jj are crucial in determining the value of important quantities like the DOS, we consider the associated cumulant generating function F( q) = ln⟨exp(∑_j q_j Im G_jj)⟩, 
which can be expressed as

    F( q) ≃⟨ln∏_α,j⟨ e^q_jε v^2 / [(λ-λ_α)^2 +ε^2]⟩_v ⟩_{λ_α}

Here we have made two approximations. The first is to treat the v_α,j for different α and j as independent, thus ignoring normalisation and orthogonality constraints on the eigenvectors. This is plausible as the number of constraints is much smaller than the number of variables v_α,j, producing only weak correlations. To see this, note that there are O(Nρε) eigenvectors contributing significantly to (<ref>). These have O(N× Nρε) components, while the number of orthonormality constraints between them is O(N^2ρ^2ε^2) and so smaller by a factor ε.
The second approximation in (<ref>) is that we are neglecting correlations between eigenvalues and eigenvectors. This again seems plausible given that the eigenvalues λ_α that contribute lie within a small range of O(ε) around λ where the statistics of the v_α,j should change little.

To evaluate F( q) it now remains to average over the λ_α. As the number of contributing eigenvalues is O(Nρε) and hence large, small fluctuations of the eigenvalues around their mean positions should be immaterial. We therefore approximate the λ_α as lying on a linear grid with the relevant spacing (Nρ)^-1, and in the same spirit replace the sum over α by an integral, shifting its origin so that α=0 designates the eigenvalue closest to λ. Then 
(<ref>) becomes

    F( q) ≃∑_j 
    ∫ dα ln⟨ e^q_jε v^2 / [(α/(Nρ))^2 +ε^2]⟩_v = 
    ∑_j 
    ∫ dα ln⟨ e^β_α(q_j)v^2⟩_v

where we have defined

    β_α(q) = N^2ρ^2ε q/α^2 + N^2ρ^2ε^2

Expanding in the q_j now gives a conventional cumulant expansion

    F( q) ≃
    ∑_j,n∫ dα β_α^n(q_j)/n!𝒦_v^2^n

where we denote by 𝒦^n_x the n-th cumulant of x. The remaining integral is

    ∫ dα β_α^n(q)=
    q^n N ρε^-(n-1)c_n

where c_n = ∫dx (1/(1+x^2))^n, so that

    F( q) ≃∑_j,nq_j^n/n!
    N ρε^-(n-1)c_n
    𝒦_v^2^n

Picking out the term of order q_j^n finally leads to the following correspondence between the n-th cumulants of Im G_jj and v^2

    𝒦_Im G_jj^n ≃ Nρ(λ)ε^-(n-1)c_n 𝒦_v^2^n



We can now use the above general result to relate the second cumulant of v^2 to the IPR: from the definition (<ref>) we have I_2(λ)≃⟨ v^4⟩/(N⟨ v^2⟩^2), while generally 𝒦_v^2^2=⟨ v^4⟩ - ⟨ v^2 ⟩^2. Imposing the eigenvector normalization condition ⟨ v^2 ⟩ = 1/N we obtain

    𝒦_v^2^2≃1/N^2 (NI_2(λ) - 1)

When I_2(λ) = 𝒪(1), i.e. in a localised part of the spectrum, the first term dominates for large N and we obtain 𝒦_v^2^2≃ I_2(λ)/N. Equation (<ref>) evaluated at second order then gives a formula for the IPR that is N-independent:

    I_2(λ) = lim_ε→ 0lim_N→∞2ε/πρ(λ)Var[Im G_jj]

where we have noted explicitly the order of limits involved. In the large N-limit taken, 
one expects I_2 to be self-averaging and thus approach some value I̅^⋆_2. Substituting (<ref>) into (<ref>), in the infinite-N limit we finally obtain the IPR estimate (<ref>) in the main text.

Figure <ref>-left shows the IPR estimated via (<ref>) and (<ref>) using different values of ε, at temperature T=1.5 and for average connectivity c=5; the total and extended DOS are also included. We observe that the IPR scales linearly with ε within the bulk of the spectrum, which is as expected for (<ref>). Our alternative estimate (<ref>) is directly applicable only within the localised part of the spectrum, but also turns out to be O(ε) for extended states. 
In finite systems, the IPR for extended states is O(1/N). Intuitlvely, one can therefore say that
in the population dynamics algorithm, which assumes N→∞, the “regulariser” ε effectively plays the role of the inverse system size, 1/N.

In the localised part of the spectrum, figure <ref>-left shows that both IPR estimates are of order unity, though I̅_2(λ)=1 throughout for ε→ 0 while I̅_2^⋆(λ) remains below unity as one would expect physically (an average IPR of one would require all eigenvectors to be localised onto a single node, which is not even true in the mean-field limit). We note, however, that I̅_2^⋆(λ) can be written as


    I̅_2^⋆(λ) = lim_ε→ 0 2ε⟨ x^2⟩ - ⟨ x⟩^2/⟨ x⟩

where, using the notation of (<ref>),

    x = ε+A_r/(ε+A_r)^2+(λ+A_i)^2

The second (mean squared) term in (<ref>) is irrelevant in the limit ε→ 0, because ⟨ x⟩=πρ(λ) is of order unity.
In the localised part of the spectrum A_r=0 so x simplifies to x=ε/(ε^2+(λ+A_i)^2). Here, the remaining term

    I̅_2^⋆(λ) = lim_ε→ 0 2ε⟨ x^2⟩/⟨ x⟩

makes clear that the second moment of x must be O(1/ε) and hence significantly larger than the squared mean. In fact, if the distribution of A_i approaches a smooth limit ρ(A_i) for ε→ 0, then (for small ε)

    ⟨ x⟩ = πρ(-λ),     ⟨ x^2⟩ = π/2ερ(-λ)

because both are given by integrals that are sharply peaked at A_i=-λ. Finally (<ref>) and (<ref>) give

    I̅_2^⋆(λ) = 1

This argument is confirmed by the data shown in figure <ref>-right: on the top we can clearly see that the histogram of A_i is smooth everywhere, and in particular around the value of -λ that dominates the computation. The plot at the bottom shows the convergence (green points) of I̅_2^⋆ when epsilon decreases; here we have used the same data as in the plot for ρ(A_i). The reason why the last few points (in black) drop to zero is that the amount of data collected was enough to give a smooth histogram on a scale as small as 10^-4, but not less. 

In conclusion, our estimate for the average IPR is expected to give a value of unity in the localised region of the spectrum, in the limit ε→ 0. Surprisingly this is the same result given by the Bollé et al formula, even though the latter is based on the opposite assumption of the estimate of the IPR being dominated by a single eigenvector.
 






§ NUMERICAL RESULTS FOR THE IPR




In section <ref> we discussed briefly the effect of the symmetrisation (<ref>) on the localisation properties of the eigenvectors, focussing particularly on the ground state. Recall that the right, left and symmetrised eigenvectors, respectively 𝐮_α, 𝐰_α and 𝐯_α, are related via 𝐏_eq as 𝐯_α = 𝐏_eq^-1/2𝐮_α = 𝐏_eq^1/2𝐰_α, where (𝐏_eq)_ii=p_i^eq∝τ_i = exp(β E_i) and (𝐏_eq)_ij=0. In the infinite temperature limit (RW) the symmetrisation is immaterial as 𝐏_eq reduces to the identity matrix, which implies 𝐯_α = 𝐮_α = 𝐰_α. In appendix <ref> we have discussed the localisation in the mean field limit where the system has no spatial structure and the IPR of symmetrised or non-symmetrised eigenvectors is dominated by the pole in -λ (see equations (<ref>, <ref>)); the factor τ_i^-1/2 in the numerator of the symmetrised case does not affect the value of the IPR qualitatively. Likewise, we generally expect that multiplying element-wise the eigenvectors by a smooth function of the energy will not change the qualitative behaviour of I_2(λ). This idea is confirmed by the numerical results presented in this appendix. Figure <ref> shows the IPR of left, right and symmetrised eigenvectors across the entire λ-range (except for the ground state λ = 0) for the case of random regular graphs with mean connectivity c=5 and temperature T=1.5. We observe that the different choices of eigenvectors have qualitatively the same localisation behaviour, except in the range of small r=-λ, where there is a natural crossover to the ground state value. We also observe that the IPR values in the scatterplots are mostly concentrated on their bin-wise average (bottom-right). The latter also indicates that symmetric and right eigenvectors have almost overlapping values of I_2. Figure <ref> shows the average IPR of left, right and symmetric eigenvectors across the λ-range, for different values of the system size N, together with the DOS and the extrapolated mobility edge from population dynamics (see section <ref>).  The N-dependence is as expected: the average IPR scales as 1/N in delocalised regions and it is of order 1 for localised regions. The scaling with 1/N is illustrated by the horizontal lines in the right plot. These are separated by a factor of 2 on the y-axis corresponding to the change from N=2000 to N=4000. The separation on the r-axis between the two regimes (localised and delocalised) is consistent with the extrapolated mobility edge (see particularly the right plot), though the decrease of the IPR towards values of order 1/N is slow in the delocalised regime near the mobility edge.






