Complex networks and glassy dynamics: walks in the energy landscape

arXiv:1101.3490v2 [cond-mat.stat-mech] 31 Mar 2011 Complex networks and glassy dynamics: walks in the energy landscape Paolo Moretti1, Andrea Baronchelli1, Alain Barrat2,3, and Romualdo Pastor-Satorras1 1Departament de F´ısica i Enginyeria Nuclear, Universitat Polit`ecnica de Catalunya, Campus Nord B4, 08034 Barcelona, Spain 2Centre de Physique Th´eorique (CNRS UMR 6207), Luminy, 13288 Marseille Cedex 9, France 3Complex Networks Lagrange Laboratory, Institute for Scientific Interchange (ISI), Torino, Italy Abstract. We present a simple mathematical framework for the description of the dynamics of glassy systems in terms of a random walk in a complex energy landscape pictured as a network of minima. We show how to use the tools developed for the study of dynamical processes on complex networks, in order to go beyond mean-field models that consider that all minima are connected to each other. We consider several possibilities for the transition rates between minima, and show that in all cases the existence of a glassy phase depends on a delicate interplay between the network’s topology and the relationship between energy and degree of a minimum. Interestingly, the network’s degree correlations and the details of the transition rates do not play any role in the existence (nor in the value) of the transition temperature, but have an impact only on more involved properties. For Glauber or Metropolis rates in particular, we find that the low-temperature phase can be further divided into two regions with different scaling properties of the average trapping time. Overall, our results rationalize and link the empirical findings about correlations between the energy of the minima and their degree, and should stimulate further investigations on this issue. PACS numbers: 89.75.Hc,05.40.Fb,64.70.Q- Complex networks and glassy dynamics: walks in the energy landscape 2 1. Introduction In the last decade, studies about the structure and dynamics of complex networks have blossomed, thanks in particular to the versatility of the network representation, which has turned out to be adequate for systems as diverse as the Internet or social networks. A large body of knowledge about the empirical description of networked systems has thus been accumulated, together with a wealth of modeling techniques; a good level of understanding of how dynamical processes taking place on networks depend on their structure has been as well reached [1, 2, 3, 4, 5, 6]. Many network studies have been concerned with systems of interest in several scientific areas a priori remote from physics (social sciences, biology, computer science, epidemiology, ...), and they have also reached more traditional fields of statistical physics, such as the study of glassy systems, as we now describe. The many puzzles raised by the glass transition, and in particular the slow dynamics displayed by glassy systems at low temperatures, have been the subject of a large interest in the past decades [7, 8]. One of the approaches which has led to promising insights consists in the description of the dynamics of a glassy system inside its configuration space. The energy landscape of a glassy system is typically rugged, made of many local minima (metastable states), whose huge number makes it difficult to reach equilibrium. In this framework, the energy landscape is seen as a set of basins of attractions of local minima (“traps”), and the system evolves through a succession of harmonic vibrations inside traps and jumps between minima [9, 10]. This picture has stimulated the definition and study of various simplified models of dynamical evolution between traps, in order to reproduce the phenomenology of glassy dynamics [11, 12, 13, 14, 15, 16, 17]. On the other hand, several studies have focused on obtaining a better understanding of the structure of these local minima. A way to attain this goal is to perform numerical simulations of small systems, at a fixed temperature, and quenching them at regular time intervals in order to make them reach the nearest local minimum. Information is then gathered on the various local minima, and on the sizes of their basins of attraction. Various studies have investigated, among other issues, the detailed structure of the potential energy landscape, the substructure of minima, and the properties of energy barriers between minima [18, 19, 20]. Several works have also used the information on the energy landscape to study a master equation for the time evolution of the probability to be in each minimum. The considered systems range from clusters of Lennard-Jones atoms to proteins or heteropolymers [9, 10, 21, 22, 23]. An interesting property of the modeling of the energy landscape in terms of a set of traps linked by energy barriers, lies in the possibility to define and study its network representation within the context of network theory. In this representation, each local minimum is associated to a node, and a link is drawn between two nodes whenever it is possible for the system to jump between the basins of attraction of the corresponding minima. The links can then be defined as weighted and directed, as jumps between minima are not equiprobable, and may be easier in one direction than in the other. Complex networks and glassy dynamics: walks in the energy landscape 3 Networks of local minima of the energy landscape have thus been built and studied. These networks have been found to exhibit a small-world character [24]. The number of links of each node (its degree) turns out to be strongly heterogeneous, possibly with scale-free degree distributions, which have been linked to scale-free distributions of the areas of the basins of attraction [25, 26, 27]. Complex network analysis tools have also been used to investigate the structure of energy landscapes of various systems of interest, such as Lennard-Jones atoms, proteins, or spin glasses, among others [21, 22, 23, 27, 28, 29, 30, 31, 32, 33]. The energy of a minimum and its degree (i.e., the number of other minima which can be reached from this minimum) have been shown to be correlated, as well as the barriers to overcome to escape from a minimum. In particular, a logarithmic dependence of the energy of a minimum on its degree has been exhibited, as well as energy barriers increasing as a (small) power of the degree of a node [23, 25, 27]. No systematic study of these issues has however been performed, and most investigations have been limited to relatively small systems because of computational limitations. Most importantly, the investigations cited above have focused on the topology of the network of minima, conceived as a tool to characterize the energy landscape. The structure of a network has however a deep impact on the properties of the dynamical processes which take place on it [6]. It seems thus adequate to put to use the tools and techniques developed for the analysis of dynamical processes on networks to achieve a better understanding of how the energy landscape structure, represented as a network, affects the system performing a random walk in it, and how the onset of glassy dynamics can be described in this way in a general framework. In a previous paper [34], we have made a first step to fill this gap by focusing on the trap model put forward in Ref. [11]. In this paper, we generalize our approach to more involved transition rates between energy minima. We show how the heterogeneous mean-field (HMF) theory [6, 35] can be used in this context to highlight the connection between the topological properties of the network of minima and the dynamical exploration of these minima. We show in particular that the relationship between energy and degree of the minima is a crucial ingredient for the existence of a transition and the subsequent glassy phenomenology. Our results shed light on the empirically found relationship between the energy of a local minimum and its degree, and we hope that they will stimulate more systematic investigations on this issue. We have organized our paper as follows: In Section 2 we define our model of energy landscape dynamics as a random walk on a complex network. Different physical transition rates are proposed, and the corresponding numerical implementation is discussed. In Section 3 we present a theoretical analysis based on the heterogeneous mean-field approximation for dynamical processes on complex networks. This formalism is applied in Section 4, where general analytical approximate expression are presented for the main quantities characterizing the glassy transition and dynamics. These expressions are applied to the different physical transition rates considered in Section 5, where checks against numerical simulations are also presented. In Section 6 we discuss the Complex networks and glassy dynamics: walks in the energy landscape 4 Figure 1. Potential energy landscape description. Energies Eiare measured positive downwards. Energy gaps Σij are defined positive upwards. relation between energy basins and energy barriers. Finally, in Section 7 we present our conclusions. 2. Random walk models on complex energy landscapes 2.1. Definition We consider a network of Nnodes, in which each vertex icorresponds to a minimum in the energy landscape, and a link is drawn between two minima iand jif the system can jump directly from ito j. To each node iis associated the energy −Eiof the corresponding minimum (energies are defined from a reference level, in such a way that Ei>0 for all i). Moreover, an energy gap Σij is associated to the edge between vertices iand j, as depicted in Fig. 1: Σij is a symmetric function, such that the energy barrier that must be overcome to jump from vertex ito vertex jcan be written as ∆Eij =Ei+ Σij and, analogously, ∆Eji =Ej+ Σij. Obviously, we will have in general ∆Eij 6= ∆Eji. The system under investigation is pictured as a walker exploring the network through a biased random walk. The rate (probability per unit time) ri→jto go from vertex ito vertex jdepends a priori on the energy at vertices iand jand/or the energy barrier between iand jthat must be overcome. The random walk model is defined in discrete time tas follows: •At time t, the walker is in vertex i. •It chooses at random a neighbor of i, namely j. •With a probability ri→j, that depends on the energy Eiand/or on the energy barrier ∆Eij, the walker hops to vertex j. •Time is updated t→t+ 1. Complex networks and glassy dynamics: walks in the energy landscape 5 The relationship between the probabilities ri→jand the energy and energy barriers can be of different forms. In usual unbiased random walks, ri→jis a constant independent of both iand j[36]. As a first step to introduce a dependence on the nodes, a possible approach is that in which the energy barriers depend only on the local minima themselves, i.e. we consider Σij = 0. For example, in the Bouchaud trap model considered in Ref. [11], the probability to exit from a trap is just an Arrhenius law depending only on the departing trap’s depth, namely rtraps i→j=r0e−βEi,(1) where β= 1/T is the inverse temperature and r0is a constant that determines a global timescale. Other possible definitions include the Metropolis one rMetropolis i→j=r0min 1, eβ(Ej−Ei),(2) and the Glauber rate rGlauber i→j=r0 1 + e−β(Ej−Ei).(3) We note that the rates considered in the Bouchaud trap model are quite different from the Metropolis or Glauber rates. Indeed, while the former depends only on the depth of the originating trap, the latter depend also on the energy of the arriving vertex. This translates in the fact that, in the limit of zero temperature, the dynamics of the Bouchaud trap model is frozen for any Ei, while Metropolis and Glauber dynamics still allow jumps to lower energy minima [13]. Within an even more realistic representation of glassy dynamics, one can also contemplate the case Σij 6= 0, allowing for the transition rates to depend explicitly on the energy barriers between adjacent minima. As a paradigm of this choice, we propose a rate of the Arrhenius form rbarriers i→j=r0e−β∆Eij ,(4) which acts a straightforward generalization of the local transition rate (1). The case of rate (1) (local trapping) was studied in a previous publication [34]. In the following, we will consider in turn non-local rates (2), (3) and (4) and discuss the fundamental differences due to the introduction of energy barriers in the model. We emphasize that our model differs both from usual unbiased random walks, as the local energy determines the transition rates, and from mean-field trap models in which jumps between any pair of energy minima are a priori possible: Here, the system can jump only between neighboring nodes. The dynamical evolution depends therefore both on the network topology and on the energies associated to the nodes. 2.2. Numerical implementation To implement numerically the random walk, it is convenient to resort to the techniques developed for general diffusion processes on complex weighted networks [37]. The main advantage of this method is to avoid rejection steps, thus improving dramatically the computational efficiency [38, 39]. Therefore, at each simulation step the random walker Complex networks and glassy dynamics: walks in the energy landscape 6 sitting at node iselects a neighbor jwith probability ri→j/Pjri→j, where the sum in the normalizing factor is extended to all of i’s neighbors. As the walker hops on node j, the physical time is incremented by an interval ∆tdrawn from the exponential distribution P(∆t) = 1/∆texp(−∆t/∆t), where ∆t=ki/Pjri→jis the inverse of the average escape rate out of node i. In this way the simulation time is disentangled from the physical time and the latter has no impact on the simulation efficiency. No matter how much physical time a walker spends in a node, from the simulation time point of view it is always just a time step. The network substrates on which we will focus are scale-free networks with a degree distribution of the form P(k)∼k−γand 2 < γ ≤3. We will generate them with the Uncorrelated Configuration Model (UCM) [40], that allows us to tune the degree distribution to the desired form and prevents the formation of degree-degree correlations. Networks are therefore generated as follows: A number of stubs (or semilinks) extracted from the desired final degree distribution is assigned to each node. Stubs are then randomly paired to form links between nodes, with the prescription that multiple links as well as self-loops must be avoided. A minimal degree mis fixed a priori. To avoid spurious effects due to the possible presence of tree-like structures [41] it is convenient to adopt m > 2. We will choose m= 4 in all of our simulations. So far the algorithm coincides with the one of the Configuration Model [42], but the UCM introduces moreover a cutoff to the degree distribution, kc=N1/2, which avoids the formation of degree correlations by limiting the size of the hubs [40]. 3. Heterogeneous mean-field theory In order to gain analytical understanding of the role of the different transition rates in the corresponding glassy dynamics, we apply a standard heterogeneous mean-field (HMF) formalism [6, 35] . The basic tenet of HMF is the assumption that all the dynamical properties of a vertex depend only on its degree. Vertices are thus grouped into classes according to their degree, and vertices with the same degree are treated as equivalent. This approximation is consistent with previous findings that have uncovered the correlations between the energy of a local minimum and the degree of the corresponding node in the network [25]. We therefore make the assumption that there exists a relationship Ei=h(ki) where the function h(k) is a characteristic of the model. This also means that the distributions of energies ρ(E) of the system’s landscape, and the degree distribution P(k) of the corresponding network are linked through h. In the same spirit, we make the further assumption that the energy gap between minima iand jdepends only on the degrees of iand j, i.e., that it can be written as Σi,j =σ(ki, kj), where σ(k, k′) is a symmetric function of kand k′. Under the HMF approximation the dynamics will thus focus on the transitions between different degree classes. The rate to go from a vertex kto a vertex k′can be written as Wkk′=P(k′|k)r(k→k′).(5) Complex networks and glassy dynamics: walks in the energy landscape 7 The function P(k′|k), defined as the conditional probability that a vertex of degree k is connected with another vertex of degree k′[43], takes into account the topological features of the network, by gauging the probability of selecting a vertex k′as neighbor of k. The function r(k→k′) measures the rate of jumping from a vertex of degree kto a vertex of degree k′(given that they are connected by an edge), and depends on kand k′through the rates ri→jand the functions hand σ. Obviously, the rate r(k→k′) is not in general a symmetric function of kand k′. It is worth noting that, apart from a normalization, Equation (5) is simply the so-called weighted propagator describing the probability that a node in class kinteracts with a node in class k′[37]. We also note that the rates r(k→k′) depend on the inverse temperature βthrough the microscopic rates ri→j. 4. General HMF formalism In this section, we apply the HMF theory to compute different quantities relevant for the characterization of the dynamics of a random walk in a complex energy landscape represented in terms of a network of minima. 4.1. Occupation probability The description of a random walk dynamics starts from the occupation probability P(k, tw), defined as the probability for the walker to be in any node of degree kat a time tw. Its time evolution can be easily represented in terms of a master equation of the form ∂P (k, tw) ∂tw ≡˙ P(k, tw) = −X k′ Wkk′P(k, tw) + X k′ Wk′kP(k′, tw).(6) Upon describing the state at time twwith the row vector P(tw) = {P(1, tw), P(2, tw),..., P(kc, tw)}, where kcis the cutoff or largest degree in the network, Equation (6) can be rewritten in vector form as ˙ P(tw) = −P(tw)L , (7) where the matrix L, with elements Lk′k= δk′kX l Wkl −Wk′k!,(8) is a generalization of a Laplacian matrix to the case of directed weighted graphs. The matrix elements satisfy Lk′k′=X k,k6=k′ Lk′k,(9) which ensures conservation of probability and states that the columns of Lare not linearly independent. The real part of every eigenvalue of Lis non-negative [44]. As a consequence, all solutions of Eq. (7), which can be formally written as P(tw) = P(t0)e−L(tw−t0),(10) Complex networks and glassy dynamics: walks in the energy landscape 8 are stable according to Lyapunov criteria. In particular, since det L= 0, Lalways has the eigenvalue 0, which corresponds to a constant solution of the problem. At this point one can proceed in close analogy with discrete-time regular Markov chains [45]. By making the assumption that the matrix Wk′kis non-negative and irreducible (indeed it is for every choice of r(k′→k) in the following), we can prove that the 0 eigenvalue of Lhas algebraic multiplicity 1. Hence, the stationary solution of Eq. (7) is unique. 4.1.1. Steady state In order to calculate the steady solution P∞in the limit tw→ ∞, one can impose ˙ P(tw) = 0. This leads to the condition P∞L= 0,(11) so that we are left with the task of finding the left nullspace of L. Eq. (11) is a homogeneous system of algebraic linear equations. It admits non-trivial solutions since det(L) = 0. In our case, the solution to (11) can be easily found by imposing the detailed balance condition. Namely, writing Eq. (11) as X k′ [−Wkk′P∞(k) + Wk′kP∞(k′)] = 0,(12) we can obtain a solution by imposing that the terms inside the summation in Eq. (12) cancel individually, that is Wkk′P∞(k) = Wk′kP∞(k′),∀k, k′.(13) Substituting the form of Wkk′, we obtain P∞(k) P∞(k′)=Wk′k Wkk′ =P(k|k′)r(k′→k) P(k′|k)r(k→k′)=kP(k) k′P(k′) r(k′→k) r(k→k′),(14) where in the last step we have used the degree detailed balance condition kP(k)P(k′|k) = k′P(k′)P(k|k′) which simply expresses that the number of edges from a node of degree kto a node of degree k′is equal to the number of edges from a node of degree k′to a node of degree k[46]. From Eq. (14), we see that its right-hand-side must be expressible as a simple ratio of a function of kover a function of k′. A general way to obtain this is to impose a coarse-grained rate r(k′→k) taking the general form r(k′→k) = f(k′)g(k)s(k′, k).(15) In other words, we assume that the rate r(k′→k) can be written as the product of a function of k′, a function of k, and a symmetric function s(k′, k) = s(k, k′) (where kand k′need not be separable). We will see later that all the rates ri→jdefined in Sec. 2.1 (traps, Glauber, Metropolis, and energy barriers) can be written in such a form. The stationary solution is then given by P∞(k) = 1 ZkP(k)g(k)/f(k) (16) where Zis a normalizing constant determined by the condition PkP∞(k) = 1. Such a solution is unique, as proven above. Interestingly, the symmetric function s(k′, k) does not enter the steady solution, although it will play a role in affecting the transient behavior, as we will see in the next sections. Complex networks and glassy dynamics: walks in the energy landscape 9 4.1.2. Glassy phase The steady state solution found above for the occupation probability is defined if and only if the normalization constant Z=X k kP(k)g(k)/f(k).(17) is finite. When this condition is met, the random walker reaches an equilibrium state with a distribution Peq =P∞. On the other hand, whenever such condition is not met, the random walker is unable to reach a steady state, i.e. the steady solution to the rate equation does not correspond to any physical steady state in equilibrium Peq. We identify this region of the phase space with the glass phase for our random walker [11]. The functions fand gdepend on the temperature, on the precise dynamics chosen (traps, Glauber, Metropolis), and encode the relationship hbetween energy and degree of the minima. The degree distribution moreover enters explicitly the expression Z. As the various parameters of the model are changed, it is thus a priori possible to go from one phase in which Zis finite to one in which Zdiverges. In a physical system in particular, the control parameter is usually the temperature, while the topology of the network of minima and the function hare given. It is then clear from Eq. (17) that the presence or absence of a finite glass transition temperature βc, such that Z becomes infinite for β≥βc, depends on the interplay between the topology of the landscape network (as determined by P(k)) and the functions fand g. Interestingly, at this mean-field level, the existence of a transition does not depend on the network degree correlations, since the conditional probability P(k′|k) do not enter Eq. (17). Let us consider for instance a network of minima with a heavy-tailed degree distribution such as P(k)∼k−γ. A transition between a finite and an infinite Z can be observed if and only if g(k)/f(k) behaves at large kin the form ∼kawhere the exponent adepends on the temperature, and can take values smaller or larger than γ−2 depending on the temperature. Another example is given by a stretched exponential form for P(k), P(k)∼e−bka, in which case a transition is observed if and only if g(k)/f(k) is of the form eb′ka, with b′a function of the temperature (the transition is then given by b′(βc) = b). 4.1.3. Glassy dynamics In any finite system, unless the product function g(k)/f(k) exhibits some sort of singularity, the normalization constant Eq. (17) is finite and the steady state distribution P∞(k) exists, the occupation probability P(k, tw) converging to it after an equilibration time, i.e. lim tw→∞ P(k, tw) = P∞(k).(18) The corresponding thermalization of the occupation probability occurs in a way depending on the function h. Shallow energy minima are indeed explored first, while deep traps (large E) are visited at larger times [11, 13]. If his a growing function of k, as indeed found empirically [25], small degree nodes correspond to shallow minima, and deeper minima are associated to larger nodes. The evolution of P(k, tw) takes then place in a hierarchical fashion: The small degree region equilibrates first, and Complex networks and glassy dynamics: walks in the energy landscape 16 100101102103 10-6 10-3 100 tw=2 tw=3 tw=4 tw=5 tw=6 10-2 100102 10-2 100 kw P(k; tw) tw=2.5 tw=5.0 tw=7.5 tw=10.0 tw=12.5 tw=15.0 10-2 100102 k / kw 10-2 100 tw=2.5 tw=5.0 tw=7.5 tw=10.0 tw=12.5 tw=15.0 β=0.25 β=2.0 β=3.0 Figure 4. Data collapse for the time evolution of the occupation probability P(k;tw) at different temperatures (Glauber dynamics). Data refer to UCM networks with N= 106and γ= 2.5, so that βc=γ−2 = 0.5 (where we have taken E0= 1). The top panel presents data for β < βcwhile both middle and bottom panels concern the low temperature case β > βc. Accordingly, for the top panel we use kw∼t1/β w∼t4 wfor the rescaling, while for the central and the bottom ones it holds kw∼t1/βc w=t2 w(see Eq. (49)). The curves corresponding to different twcollapse well under this rescaling. Note that we use rather small values of tw, because the equilibration time, defined by kw∼kc∼N1/2, is teq ∼Nβc 2≃32 for β > βcand teq ∼Nβ 2≃6 for β < βc. Each curve is obtained by averaging over 3 ×106simulation runs. which leads to τk=1 rk ∼(kβcE0β > βc kβE0β < βc..(48) From here, using the relation τkw∼tw, we obtain kw∼(t1/(βcE0) wβ > βc t1/(βE0) wβ < βc..(49) In Fig. 4 we check the validity of Eqs. (45) and (49) by performing a data collapse analysis for different values of tw. The curves obtained for different twcollapse indeed as predicted. In the case of the Metropolis transition rates, a similar analysis yields rMetropolis k∝1 (βc−β)E0 k−βE0−β βc k−βcE0!,(50) leading to the same asymptotic behavior as in Eq. (47), and therefore to the same scaling picture as for the Glauber rate. As pointed out for the case of the traps model [34], however, in finite systems the scaling relations in Eq. (49) hold only as far as kw(tw)< kc, i.e. it exists an equilibration time teq, obtained by inverting Eq. (49), above which the system has completely relaxed and Eq. (45) is no more valid. Finally, it is worth stressing that, while for large Complex networks and glassy dynamics: walks in the energy landscape 17 temperatures the scaling exponent relating kwand twdepends on the temperature, in the low temperature phase it becomes independent, being just proportional to the transition temperature. We note that this saturation of the exponent at 1/(βcE0) is very different from the phenomenology obtained in the trap model [34], for which kw∼t1/(βE0) w. An immediate consequence is that the equilibration time strongly depends on β, as teq ∼k(βE0) c, for a system described by traps, but is given by teq ∼k(βcE0) c≪k(βE0) cfor any β > βcfor Glauber and Metropolis rates. Contrarily to results for the glass transition temperature Tcand the steady state, the glassy dynamics for barrier-mediated rates does not yield the same results as for Glauber and Metropolis rates, since rkdoes depend on the symmetric function σ(k, k′). In particular, we need here to choose a functional form for σ. We propose to use σ(k, k′) = σ0(kµ+k′µ),(51) which will be justified in Section 6. In this case we obtain rbarriers k∝r0k−βE0e−βσ0kµ.(52) In order to keep an interesting phenomenology, the constant σ0cannot be chosen arbitrarily. If σ0were independent of the system size, the escape rate would be dominated by the exponential at all temperatures. This behavior would reflect the fact that in this case the rate is suppressed in transitions involving nodes of large k, eventually generating unphysically large barriers at large kc. To prevent the system from building up infinite barriers, one can impose E0ln kc∝σ0kµ c,(53) such that the maximum barrier is always comparable to the lowest energy minimum and neither term dominates the other. As a consequence, we take σ0=ǫE0 ln kc kµ c ,(54) where ǫis now constant and size independent. Contrarily to the previous cases, kw(tw) is hard to determine, as no explicit inversion of Eqs. (21) and (52) can be provided for the range of parameters of interest in our study. A numerical evaluation of kw(tw), is reported in Figure 5. The maximum degree of equilibrated vertices kwhas an initial power law increase in time, which is reminiscent of the local trapping model. However, as larger degree nodes are equilibrated, exponential barriers come into play and the hierarchical thermalization becomes logarithmic in time. 5.4. Average escape time The average escape time, defined as the average time required by the system to escape from the vertex it occupies, can be computed in the long time limit from Eq. (22), as a function of the average trapping time τkin vertices of degree k. From the asymptotic expansions of τkin Eq. (48), valid for Glauber and Metropolis dynamics, evaluation of Eq. (22) allows us to observe that, whenever a finite Tcexists, tesc(tw→ ∞) diverges Complex networks and glassy dynamics: walks in the energy landscape 18 1001021041061081010 tw 100 101 102 103 kw β=0.4 β=0.8 β=1.2 β=1.6 Figure 5. Maximum degree of equilibrated nodes up to time twfor barrier-mediated dynamics. Curves are obtained from the numerical inversion of Eq. (21). Values of parameters are: γ= 2.5, βc= 0.5, E0=ǫ= 1, kc= 103,µ= 0.5. 100103106 100 101 N=103 N=104 100103106 tw 10-2 100 100103106 10-2 100 tesc / tesc eq N=105 N=106β=0.25 β=0.75 β=2.0 Figure 6. Rescaled average escape times (Metropolis dynamics). Data for UCM network with γ= 3.0, so that βc= 1.0 (E0= 1). The scaling forms of Eq. (55) produce a collapse of the curves concerning different systems sizes in the three regimes of high, intermediate and low temperature (top, center and bottom panels, respectively). The slightly worse collapse obtained for β= 0.75 may be due to logarithmic corrections as βis close to both βcand 2βc. Each point is averaged over 400 simulation runs (20 runs on each of 20 network realizations). at 2Tcin an infinite system, as was already observed in the case of local trapping [34]. It is noticeable that the same divergence temperature is obtained, as Eq. (22) a priori involves the network’s degree correlations and the function s. Within the continuous degree approximation, the divergence of the escape time with the system size follows Complex networks and glassy dynamics: walks in the energy landscape 19 10-1 βc1002βc101 β 100 102 104 106 108 1010 <τ> kc=1012 kc=106 kc=109 kc=103 Figure 7. Average rest time for the Glauber dynamics in a scale-free uncorrelated network with γ= 2.75, as a function of the inverse temperature, for different system sizes. Data are obtained by numerical computation of Eq. (28) (with E0= 1). Note the exponential increase with βfor βc< β < 2βc, which saturates for β > 2βcas predicted by Eq. (61). the scaling laws teq esc ≡tesc(tw→ ∞)∼ZkcP∞(k)τk∼       kβcE0 cβ > βc k(2β−βc)E0 cβc/2< β < βc const. β < βc/2 .(55) As noted in the previous paragraph for the equilibration time, we note that the scaling for β > βcdiffers from the form kβE0 cencountered for local trapping [34]. Figure 6 reports simulation data that confirm the validity of Eq. (55). As the temperature is lowered, the initial transient becomes longer, but for large enough times twthe asymptotic behavior predicted in Eq. (55) is reached, as made clear from the collapse of curves concerning different system sizes. In the case of barrier-mediated dynamics, τkdepends on the symmetric function σ(k, k′), as expressed in Eq. (52), namely τk∼kβE0eβσ0kµ.(56) Proceeding as above we obtain tesc ∼       k(1+ǫ)βE0 cβ > βc k[(2+ǫ)β−βc]E0 cβc/(2 + ǫ)< β < βc const. β < βc/(2 + ǫ) ,(57) where we recall that ǫdoes not depend on the system size. 5.5. Average rest time The HMF expression for the asymptotic average rest time, defined as the average time spent by the system in a minimum, is given by Eq. (28), namely hτi=Z/I, where the Complex networks and glassy dynamics: walks in the energy landscape 20 1031061091012 kc 100 102 104 106 <τ> ~kc βc β=2βc β=βc Figure 8. Average rest time for the Glauber dynamics in a scale-free uncorrelated network with γ= 2.5, as a function of the degree distribution cut-off and for different temperatures. Data are obtained by numerical computation of Eq. (28). hτigrows as a power-law of kc, with an exponent that grows as βincreases (going from bottom to top in the figure). The thick-gray line corresponds to β=βc, while the thick-black line corresponds to β= 2βc. For larger values of β, the power-law behavior corresponds to the predicted hτki ∼ kβcE0 c∼kγ−2 c, which no more depends on β. The kγ−2 ccurve is reported as a dashed line for reference. Values of βrepresented here are comprised between 0.25 and 3. quantities Zand I, for uncorrelated scale-free networks and a degree-energy relation h(k) = E0ln(k), take the form, in the continuous degree approximation, Z ∼ Zkck1−γ+βE0dk ∼const + k(β−βc)E0 c,(58) I ∼ Zkcdk Zkcdk′k1−γg(k)k′1−γg(k′)s(k, k′).(59) Let us first recall the case of the local trap model. Both g(k) and s(k, k′) are then constants, so that I ∼ hki2= const. Thus, the average rest time behaves as Z: it is finite for β < βc, and diverges with the system size as k(β−βc)E0 cfor β > βc, signaling the emergence of the glassy regime at low temperatures. In the case of Glauber and Metropolis dynamics (both leading to the same results), the situation is more involved, since the product g(k)g(k′)s(k, k′) entering Iis not constant. In fact, Idiverges with kcfor βE0>2(γ−2), that is, at a lower temperature given by β′ c= 2βc. The interplay of these two temperatures determines the behavior of the system for finite sizes within the low temperature phase. In particular, for the Glauber dynamics with g(k) = eβh(k)and s(k, k′) = r0/[eβh(k)+eβh(k′)], we have I ∼ const + k(β−2βc)E0 c.(60) Upon considering lower values of β,hτifirst encounters the divergence of Zat βc, which is then partially regularized by the divergence of Iat 2βc. From these results, we obtain Complex networks and glassy dynamics: walks in the energy landscape 21 the emergence of three scaling regimes for the behavior of hτias a function of the system size: hτi∞≡ hτi(tw→ ∞)∼       kβcE0 cβ > 2βc k(β−βc)E0 cβc< β < 2βc const. β < βc .(61) 10-1 βc1002βc β 100 102 104 <τ> Glauber Metropolis Barriers 10-1 100101 β 100 1012 1024 1036 <τ> Figure 9. Average rest time as a function of βfor different transition rates, for a random walker on UCM networks with γ= 2.75 (E0= 1) and N= 106. Discrete points: simulation results; continuous lines: theoretical predictions from hτi=Z/I, based on simulation parameters. The Glauber and Metropolis transition rates induce two changes of behavior at βcand at 2βc, the first being a steep increase of the average rest time hτiand the second a smoothing/saturation of this increase. No saturation of hτiis instead observed when barriers are present. The agreement with theoretical predictions is remarkable, thus corroborating the validity of the HMF assumptions. Moderate deviations are found only in the case of barriers, where exponential growth is expected to add greater fluctuations. In the inset, the difference between the two behaviors is more evident thanks to a different scale of the plot. Note that, since in the simulations E0= 1, the high temperature limit of the rest time is different for the Glauber and Metropolis dynamics, being τ(β= 0) = 2 and τ(β= 0) = 1 respectively. For the case of barriers we have chosen σ0= 10−1. Each point is obtained by averaging the rest times corresponding to the first 106hops of the random walker in each of 10 network realizations. The direct numerical computation of Eq. (28) is shown in Figs. 7 and 8, showing the validity of this analysis. In particular, the exponential increase of hτiwith βin the intermediate temperature range βc< β < 2βcis clearly apparent in Fig. 7, and Fig. 8 confirms that, for β > 2βc, the exponent in the scaling law for the system size kcdoes not depend on the temperature. While the temperature βcsignals the onset of the low temperature phase with glassy dynamics for all considered transition rates, for Glauber/Metropolis dynamics the low temperature phase can be further divided into two regions that correspond to different behaviors of the timescales with the system size. Complex networks and glassy dynamics: walks in the energy landscape 22 100103106 100 101 N=103 N=104 100103106 tw 10-3 100 100103106 10-2 100 102 <τ> / <τ>eq N=105 N=106β=0.25 β=1.5 β=3.0 Figure 10. Average rest time for Metropolis dynamics on uncorrelated scale-free networks with γ= 3.0 (E0= 1, βc= 1.0). Data for different system sizes collapse well when rescaled according to the theoretical values of Eq. (61). While the agreement is excellent both for high and low temperature (top and bottom panels, respectively), logarithmic corrections are probably present for the regime of intermediate temperatures βc< β < 2βc(central panel). In each simulation run the rest interval starting before twand ending after twis considered, and each point in the Figure is averaged over 400 simulation runs (20 runs on each of 20 network instances). Figures 9 and 10 moreover show the result of numerical simulations of random walkers on scale-free networks for Glauber and Metropolis dynamics as well as in the case of barriers, globally confirming the above discussed picture. Dynamics in the presence of barriers do not yield the same phenomenology as Glauber and Metropolis rate. In this case, we have g(k) = 1 and s(k, k′) = r0e−βσ(k,k′). Selecting σ(k, k′) = σ0(kµ+k′µ), as in Section 5.3, we are led to I ∼ "Zkcdkk1−γe−βσ0kµ#2 .(62) As for the escape time tesc, upon choosing σsize independent, the rest time hτiwill be diverging exponentially with kcat every temperature. By introducing the size dependence as in Eq. (54), instead, one can see that the Iintegral converges to a constant for large kcso that one is left with hτi ∼ (k(β−βc)E0 cβ > βc const. β < βc .(63) We therefore obtain the same picture as in the case of traps, with an exponential increase of hτias βincreases, as confirmed by numerical simulations in Fig. 9. Complex networks and glassy dynamics: walks in the energy landscape 23 6. Energy basins and energy barriers Inspired by analogies with systems governed by the Arrhenius law, we have introduced a transition rate that takes into account the energy barriers between states. Within the heterogeneous mean-field approximation, in which all variables depend only on the degree of the vertices, and choosing Ei=E0ln(ki), the transition rate we have considered becomes r(k→k′) = k−βE0e−βσ(k,k′),(64) where σ(k, k′) is a symmetric function of the degrees of the two nodes. This model represents in essence an extension of the local trap model, where non-locality enters only in the form of symmetric energy gaps Σij (see Fig. 1). The steady state has exactly the same form as the ones discussed so far, which incidentally is the same as for the local trap model. As shown in the previous Sections, the presence of barriers affect transient relaxation phenomena, but not the steady state. A different question is whether one can be more specific about the realistic functional form of the coarse-grained function σ(k, k′). In the previous section we have already introduced a definition of σ(k, k′). Here we provide the rationale behind that choice. Numerical simulations of the energy-landscape network of Lennard-Jones clusters show that the average barrier to escape from state kfollow the power law ∆Ek∼kµ, with µ > 0 [23]. In our model, such average can be computed as ∆Ek=X h P(h|k) [E0ln k+σ(h, k)] .(65) For simplicity we focus on uncorrelated networks, as simulations indeed show weak degree correlations. Under this assumption, the first term of the sum on the right-hand side of Eq. (65) will contribute as a logarithm of kand the power-law behavior of ∆Ek is possible whenever σ(h, k)∼kµ, which leads to consider the form proposed in previous sections, σ(k, k′) = σ0(kµ+k′µ),(66) where σ0has the dimensions of an energy (a discussion about the possible values of σ0is given in Section 5.3). More complicated functional forms can also be proposed, for example accounting for barriers of different signs, as long as they retain the same power-law behavior of Eq. (66) in the large klimit. It is interesting to notice that ∆Ek∼kµimplies that the average escape rate e−β∆Ekhas the form of a stretched exponential ∼exp(−βkµ), if we neglect the logarithmic correction. 7. Conclusions In this paper, we have presented a simple mathematical framework for the description of the dynamics of glassy systems in terms of a random walk in a complex energy landscape. We have shown how to incorporate into this picture the network representation of this REFERENCES 24 landscape, put forward and studied by several authors [25, 26, 27, 29, 30, 31, 32, 33], in order to go beyond simple mean-field models of random walks between traps that are all connected to each other. While our previous work had focused on the case of a landscape consisting of traps connected by a network [34], we have here generalized our study to more involved and realistic transition rates between minima, including Glauber or Metropolis rates, and the possibility of energy barriers between minima. We have shown how the interplay between the topology of the network of minima and the relationship between the energy and the degree of a minimum may determine a rich phenomenology, with the existence of two phases and of glassy dynamics at low temperature. Interestingly, the existence of these phases, and the transition temperature, do not depend on the network’s degree correlations nor on the precise form of the transition rates, but other more detailed properties do. In the case of Glauber and Metropols dynamics, the low temperature phase can be further divided into two regions with different scaling properties of the average trapping time as a function of the temperature. Overall, our results rationalize and link the empirical findings about correlations between the energy of the minima and their degree, and should stimulate further investigations on this issue. Our work has also interesting applications in terms of diffusion phenomena on complex networks, and shows that non trivial transition rates can lead to a very interesting phenomenology. Usual random walks lead to a higher probability for the random walker to be in a large degree node (∝kP(k)), with respect to the random choice of a node (∝P(k)); here, the models we have studied can lead to various stationary probabilities, such as for instance a uniform coverage which does not depend anymore on the degree. Interestingly, the biased random walks among traps that we have studied can even display a phase transition phenomenon, as either a temperature parameter or the network’s properties are changed, with the possible presence of a glassy phase with slow dynamics. Acknowledgments P. M., R.P.-S., and A. Baronchelli acknowledge financial support from the Spanish MEC, under project FIS2010-21781-C02-01, and the Junta de Andaluc´ıa, under project No. P09-FQM4682.. R.P.-S. acknowledges additional support through ICREA Academia, funded by the Generalitat de Catalunya. A. Baronchelli acknowledges support of Spanish MCI through the Juan de la Cierva program funded by the European Social Fund. References [1] Albert R and Barab´asi A L 2002 Rev. Mod. Phys. 74 47–97 [2] Dorogovtsev S N and Mendes J F F 2003 Evolution of networks: From biological nets to the Internet and WWW (Oxford: Oxford University Press) REFERENCES 25 [3] Newman M 2003 SIAM Review 45 167–256 [4] Pastor-Satorras R and Vespignani A 2004 Evolution and structure of the Internet: A statistical physics approach (Cambridge: Cambridge University Press) [5] Caldarelli G 2007 Scale-Free Networks: Complex Webs in Nature and Technology (Oxford: Oxford University Press) [6] Barrat A, Barth´elemy M and Vespignani A 2008 Dynamical Processes on Complex Networks (Cambridge: Cambridge University Press) [7] Debenedetti P and Stillinger F 2001 Nature 210 259 [8] Barrat J L, Feigelman M, Kurchan J and Dalibard J (eds) 2003 Les Houches Session LXXVII, 1-26 July, 2002 Les Houches - Ecole d’Ete de Physique Theorique, Vol. 77 (Berlin: Springer) [9] Angelani L, Parisi G, Ruocco G and Viliani G 1998 Phys. Rev. Lett. 81 4648–4651 [10] Berry R S and Breitengraser-Kunz R 1995 Phys. Rev. Lett. 74 3951–3954 [11] Bouchaud J P 1992 J. Physique I (France) 21705–1713 [12] Bouchaud J and Dean D 1995 J. Physique I (France) 5265 [13] Barrat A and M´ezard M 1995 J. Physique I (France) 5941–947 [14] Monthus C and Bouchaud J 1996 Journal of Physics A-Mathematical and General 29 3847–3869 [15] Bertin E and Bouchaud J P 2002 J. Phys. A: Math. Gen. 35 3039 [16] Bertin E and Bouchaud J P 2003 Phys. Rev. E 67 065105(R) [17] Bertin E 2003 J. Phys. A: Math. Gen. 36 10683 [18] B¨uchner S and Heuer A 2000 Phys. Rev. Lett. 84 2168 [19] de Souza V and Wales D 2009 J. Chem. Phys. 130 194508 [20] Heuer A 2008 J. Phys. Cond. Mat. 20 373101 [21] Cieplak M, Henkel M, Karbowski J and Banavar J R 1998 Phys. Rev. Lett. 80 3654–3657 [22] Bongini L, Casetti L, Livi R, Politi A and Torcini A 2009 Phys. Rev. E 79 061925 [23] Carmi S, Havlin S, Song C, Wang K and Makse H A 2009 J. Phys. A 42 105101 [24] Scala A, Amaral L A N and Barth´elemy M 2001 Europhysics Letters 55 594 [25] Doye J P K 2002 Phys. Rev. Lett. 88 238701 [26] Massen C P and Doye J P K 2005 Phys. Rev. E 71 046101 [27] Seyed-allaei H, Seyed-allaei H and Ejtehadi M R 2008 Phys. Rev. E 77 031105 [28] Doye J P K and Massen C P 2004 J. Chem. Phys. 122 084105. 14 p [29] Gfeller D, De Los Rios P, Caflisch A and Rao F 2007 Proceedings of the National Academy of Sciences 104 1817–1822 [30] Gfeller D, de Lachapelle D M, De Los Rios P, Caldarelli G and Rao F 2007 Phys. Rev. E 76 026113