Langevin approach for the dynamics of the contact process on annealed scale-free networks

arXiv:0810.3000v1 [cond-mat.dis-nn] 16 Oct 2008 Langevin approach for the dynamics of the contact process on annealed scale-free networks Mari´an Bogu˜n´a,1Claudio Castellano,2and Romualdo Pastor-Satorras3 1Departament de F´ısica Fonamental, Universitat de Barcelona, Mart´ı i Franqu`es 1, 08028 Barcelona, Spain 2SMC, INFM-CNR and Dipartimento di Fisica, “Sapienza” Universit`a di Roma, P.le Aldo Moro 2, I-00185 Roma, Italy 3Departament de F´ısica i Enginyeria Nuclear, Universitat Polit`ecnica de Catalunya, Campus Nord B4, 08034 Barcelona, Spain (Dated: November 13, 2018) We study the dynamics of the contact-process, one of the simplest nonequilibrium stochastic processes, taking place on a scale-free network. We consider the network topology as annealed, i.e. all links are rewired at each microscopic time step, so that no dynamical correlation can build up. This is a practical implementation of the absence of correlations assumed by mean-field approaches. We present a detailed analysis of the contact process in terms of a Langevin equation, including explicitly the effects of stochastic fluctuations in the number of particles in finite networks. This allows us to determine analytically the survival time for spreading experiments and the density of active sites in surviving runs. The fluctuations in the topological structure induce anomalous scaling effects with respect to the system size when the degree distribution has an “hard” upper bound. When the upper bound is soft, the presence of outliers with huge connectivity perturbs the picture even more, inducing an apparent shift of the critical point. In light of these findings, recent theoretical and numerical results in the literature are critically reviewed. PACS numbers: 89.75.Hc, 05.70.Jk, 05.10.Gg, 64.60.an I. INTRODUCTION The study of the effects of an heterogeneous topology on equilibrium and nonequilibrium dynamical processes has lately experienced an active interest from the statistical physics community [1]. Indeed, it has been observed in recent years that many natural and manmade systems are well characterized in terms of complex networks or graphs [2, 3], in which vertices represent elementary units in the system, while edges stand for pairwise interactions between elements. Most real networked systems can be characterized by a heterogenous complex topology, showing remarkable universal features, such as the small world property [4] and a scale-free connectivity pattern [5]. The small-world property refers to the fact that the average distance hℓibetween any two vertices—defined as the smallest number of edges on a path between one and the other—is very small, scaling logarithmically or even more slowly with the network size N[6]. This is to be compared to the power-law scaling hℓi ∼ N1/d in a d-dimensional lattice. Since the logarithm grows slower than any power-law function, even if dis very large, small-world networks can be thought of as highly compact objects of infinite dimensionality. On the other hand, scale-free (SF) networks are typically characterized by a degree distribution P(k), defined as the probability that a randomly selected vertex has degree k—is connected to kother vertices—that decreases as a power-law, P(k)∼k−γ,(1) where γis a characteristic degree exponent, usually in the range 2 < γ ≤3 [2, 3]. Dynamical processes taking place on top of complex networks arise in a wide variety of scientific and technological contexts. For example, we can mention the transmission of information packets on the Internet [7], the spreading of biological diseases on social networks or computer viruses in computer infrastructures [8, 9], etc. The interest in the study of these dynamics was triggered by the observation that the heterogeneous connectivity pattern observed in SF networks with diverging degree fluctuations can lead to very surprising outcomes, such as an extreme weakness in the face of targeted attacks aimed at destroying the most connected vertices [10, 11], or the ease of propagation of infective agents [9, 12]. These properties are due to the critical interplay between topology and dynamics in heterogeneous networks and are absent in their homogeneous counterparts. After those initial discoveries, a real avalanche of new results have been put forward, including classical equilibrium systems [13, 14, 15] and non-equilibrium processes such as epidemic spreading [9, 12], reactiondiffusion processes [16, 17, 18] and dynamics with absorbing states [19, 20]. For an extensive review of recent results we refer the reader to Ref. [1]. The analytical approach to the study of dynamical processes on complex networks is dominated by the application of heterogeneous mean-field theory [1]. Heterogeneous mean-field theory (HMF) is based in two basic assumptions: (i) the homogeneous mixing hypothesis, stating that all vertices with the same degree (within the same degree class) share the same dynamical properties; and (ii) the assumption that fluctuations are not relevant, and therefore analytical studies can be conducted within a deterministic approach. This last fact is in some sense natural, since the small-world property implies that dynamical fluctuations in a network are so close together 2 that they can be washed away in very few time steps1. HMF has proved to be extremely useful in providing a very accurate description of the behavior of most dynamical processes on complex networks [1]. On the other hand, in other instances, such as in the nonequilibrium contact process (CP) [22] a debate has arisen about the comparison between numerical simulations on SF networks and HMF predictions [19, 23, 24]. Underlying this controversy is the fact that while the HMF approach considers all relevant quantities as deterministic, and hence assumes an infinite system size, numerical simulations are performed on finite systems and thus are necessarily influenced by stochastic fluctuations due to the finite number of particles, in particular close to an absorbing state phase transition [22]. Finite size effects are very strong in networks and an appropriate theoretical framework for them is necessary to compare simulations with HMF results. For this reason, in Ref. [25] the CP was considered on the simplest network substrate, which is an annealed network, in which the quenched disorder imposed by the actual connections in the network is not considered. In this scenario, it was possible to deduce, by means of qualitative arguments, the correct size scaling of the CP in this kind of networks, in very good agreement with numerical simulations. In this paper we present a more detailed analysis of the CP in SF networks, deriving the corresponding Langevin equation describing its dynamics for the case of annealed networks. The analysis of this equation allows us to uncover the correct finite size scaling behavior of the CP in heterogeneous networks, providing the exact value of the critical exponents describing this system. Surprisingly, the critical behavior of the CP turns out to be extremely sensitive to the particular degree cutoff chosen for the construction of the network, in agreement with previous results obtained from a more phenomenological approach [25]. In particular, critical exponents depend explicitly on the way the degree cutoff diverges with the system size, if it scales sufficiently slowly. If the scaling is instead fast, additional complications arise and fluctuations of the degree distribution strongly perturb the picture. We have organized our paper as follows. In Sec. II we describe the main properties of annealed networks, which represent the simplest network substrate for a dynamical process, in which mean-field theory is supposed to be exact. We focus, in particular, on the effects of the maximum degree allowed on the network cutoff and on its fluctuations. Sec. III defines the CP on complex networks, whose mean-field analysis is reviewed in Sec. IV. In Sec. V we comment on the different approaches followed in the past to deal with finite size effects on the 1At variance with what happens in regular lattices below the critical dimension, where in particular, close to a critical point, dynamics is governed by fluctuations [21]. CP in SF networks. The general Langevin theory for the CP in networks is presented in Sec. VI, while Sec. VII focuses on the analysis of annealed networks. Sec. VIII discusses the meaning of finite size effects and finite size scaling in heterogeneous networks. In Sec. IX we present a digression to the case of annealed networks with outliers, that is, vertices with a degree much larger than the average maximum degree expected in the network. Finally, we draw our conclusions in Sec. X. Some technical questions are developed in several Appendices. II. ANNEALED SCALE-FREE NETWORKS The topological properties of any complex network are fully encoded in its adjacency matrix aij, taking the value aij = 1 if there is an edge connecting vertices iand j, and zero otherwise. In the so-called quenched networks, the values of the adjacency matrix are fixed in time. For large quenched networks, a statistical characterization in terms of the degree distribution P(k) and the degree correlations P(k′|k), defined as the conditional probability that a vertex of degree kis connected to a vertex of degree k′[26, 27], is useful as a compact way to express the essential features of the adjacency matrix2. Quenched networks are the typical output of most network models, such as the configuration model (CM) [28, 29, 30, 31], the uncorrelated configuration model [32], the class of models with hidden variables [33], linear preferential attachment models [5, 34], etc. In this case, each network must be considered as a representative of a statistical ensemble of random networks, which is characterized by the P(k) and P(k′|k) probability distributions. When a dynamical process takes place on top of such a network, one is considering the network as frozen, with respect to the characteristic time scale τDof the dynamics. In this case, in a numerical analysis of a dynamical process, one must consider the dynamics over many different quenched networks, all belonging to the network ensemble with the same statistically equivalent topological properties, and perform an ensemble average to compute the average dynamical quantities. In other instances, on the other hand, the very network is a dynamical object, changing in time over a certain time scale τN. In this case, the correct topological characterization is strictly statistical, given in terms of the degree distribution P(k) and the degree correlations P(k′|k). In the limit τN≪τD, that is, when the network connections are completely reshuffled between any two microscopic steps of the dynamics, while keeping fixed P(k) and P(k′|k), the resulting networks are called annealed [35, 36, 37]. Apart from the cases where they describe the actual evolution of real systems, an2A more detailed characterization can be made using higher order degree correlations, see Ref. [27]. 3 nealed networks are extremely important from a theoretical point of view, because mean-field predictions for dynamical processes on networks are usually obtained in this limit, via the so-called annealed network approximation [1]. In practice one replaces the adjacency matrix aij by its ensemble average ¯a(ki, kj), defining the probability that two vertices of degree kiand kjare connected. This average is given by ¯a(k, k′) = 1 NP(k) 1 NP(k′)X i∈kX j∈k′ aij ≡k′P(k|k′) NP(k),(2) where notation i∈kmeans summation for all vertices of degree k. Taking the case of uncorrelated networks, with P(k|k′) = kP(k)/hki[38], the simple form ¯a(k, k′) = kk′/Nhkiresults. From a numerical point of view, the simulation of dynamics on annealed networks implies the re-generation of the whole network every time a microscopic dynamic step is performed [37]. For uncorrelated networks this can be efficiently implemented in CP-like dynamics. In this case, an annealed network of size Nis completely defined by its degree sequence {k1,...,kN}, where the degrees ki are integer random numbers, extracted according to the degree distribution P(k), and restricted between a lower bound mand an upper bound M≤N. Degree correlations are given by P(k′|k) = k′P(k′)/hki. Thus, every time we need to find a nearest neighbor of a vertex, it is selected at random with probability k′P(k′)/hkiamong the Nvertices present in the network. Finite SF networks are additionally characterized by another parameter, the degree cutoff kc(N) [38], that is the average value of the actual maximum degree kmax in a single realization of the degree sequence: kc(N) = hkmaxi. In general, kcis a non decreasing function of the network size and, as we shall see below, the CP dynamics is very sensitive to its actual size dependence. Notice that the value of the cutoff plays a relevant role in the determination of degree correlations in finite quenched networks [39]. It is known that for the network to be closed without degree-degree correlations and no multiple edges or self-loops one must impose that degrees are smaller than the structural cutoff ∼N1/2[32]. In uncorrelated annealed networks, however, since they are by construction uncorrelated, such a restriction does not apply, and any cutoff is in principle possible. Simple considerations based on extreme value theory [39] give the probability distribution Pmax(kmax) of observing a maximum degree kmax among Ndegrees independently sampled from a distribution P(k)∼k−γand bounded by the constraint m≤k≤M. In the continuous degree approximation, the distribution of maximum degrees takes the form Pmax(kmax =k) = N(γ−1)(m1−γ−k1−γ)N−1 (m1−γ−M1−γ)Nk−γ(3) Using this expression one can compute explicitly the value of kc(N), obtaining two different behaviors, depending on whether M/m is larger or smaller than N1/(γ−1), namely kc(N) =    M, M m≪N1/(γ−1) mΓγ−2 γ−1N1/(γ−1),M m≫N1/(γ−1) , (4) where Γ(z) is the Gamma function [40]. For the network to be SF the upper bound of the degree distribution must diverge with the system size: M∼N1/ω. The parameter ω≥0 is in principle arbitrary, but its value strongly affects the nature of the actual maximum of the degree sequence. If Mdiverges not faster than N1/(γ−1) (i. e. ω≥γ−1), then kc=Mis a hard cutoff, with no degree larger than kc. For ω < γ −1 instead, kc∼N1/(γ−1) is just an average cutoff, but hkmax2igrows with the network size as N1+(3−γ)/ω ≫k2 c, indicating that fluctuations diverge. As a consequence, the maximum degree present in the degree sequence has wide fluctuations and outliers, i.e. nodes with a degree much larger than kc, may be present in the network. It is important to stress that taking ω=γ−1 is very different from setting M=∞from the beginning (ω= 0) or M=N(ω= 1), as it is usually done in the quenched configuration model [31]. In both cases the average cutoff kcscales as N1/(γ−1), but for ω=γ−1 this is a hard cutoff and fluctuations of the value of kmax around kcare bounded. The presence of outliers and large fluctuations in the maximum degree has a strong effect on the dynamics on annealed networks. In particular, as we will see, the relevant quantity characterizing the size effects on the dynamics is the second moment of the degree distribution g=hk2i hki2.(5) The fluctuating nature of this quantity from sample to sample can be assessed by looking at its standard deviation σg, that can be easily computed, given the uncorrelated nature of the degrees in annealed networks. Thus, we have the relative fluctuations σ2 g g2=1 Nhk4i hk2i2−1.(6) Assuming that hkni ∼ hkmaxn+1−γi, we have σ2 g g2∼hkmax5−γi Nhkmax3−γi2∼(N2(3−γ)(1 ω−1 (γ−1) ),for ω < γ −1 N(γ−1)/ω−1,for ω≥γ−1, (7) Thus, fluctuations vanish in the large size limit for ω≥ γ−1, while for ω < γ −1, the fluctuations of gdiverge as a power law with the network size N. In the rest of the paper, we will mainly discuss the simplest case ω≥γ−1, considering often the cases ω= 2 and ω=γ−1. The more delicate issue of the effect of outliers on the behavior of CP on SF networks will be 4 touched only in Sec. IX. Notice that the value ω= 2 has no special meaning here and it is just an example of what occurs for ω > γ −1. At odds with the case of quenched networks, the structural cutoff kc=N1/2does not play any role in annealed networks. III. THE CONTACT PROCESS ON COMPLEX NETWORKS We consider the contact process (CP) [22] on heterogeneous networks, which is defined as follows [19]. An initial fraction ρ0of vertices is randomly chosen and occupied by a particle. Dynamics evolves in continuous time by the following stochastic processes: Particles in vertices of degree kcreate offsprings into their nearest neighbors at rate λ/k, independently of the degree k′of the nearest neighbors. At the same time, particles disappear at rate µthat, without loss of generality, is set to µ= 1. From a computational point of view, the CP can be efficiently implemented by means of a sequential updating algorithm [19, 22]: At each time step t, a particle in a vertex iis chosen at random. With probability p= 1/(λ+ 1) the particle disappears. On the other hand, with probability 1 −p=λ/(λ+ 1), the particle may generate an offspring. In this case, a vertex j, nearest neighbor of i, is selected at random. If jis empty, a new particle is created of it; otherwise, nothing happens. In any case, time is updated as t→t+ [(1 + λ)n(t)]−1, where n(t) is the number of particles at the beginning of the time step. Notice that the factor (1 + λ) in the time update is due to the fact that each infected particle can perform two independent actions, either infect a neighbor (at rate λ) or become healthy again (at rate µ= 1). This factor was neglected in previous implementations of the CP in complex networks [19, 24, 25]. The results of these works remain, however, unaltered, since the factor is irrelevant for steady state properties and amounts only to a rescaling for time dependent properties. In Euclidean d-dimensional lattices, the CP undergoes a nonequilibrium phase transition [22] between an absorbing state, with zero particle density, and an active phase, with average constant density of particles, which takes place at a critical point λc. This phase transition is characterized in terms of the order parameter ρ, defined as the average density of particles in the steady state. Defining ∆ = λ−λc, we observe for ∆ <0, and in infinite lattices, an absorbing phase with ρ= 0. For ∆>0, on the other hand, the system sets in an active phase with a nonzero order parameter, obeying ρ∼∆β. Close to the critical point, the system is also characterized by diverging correlation length and time scales, namely ξ∼ |∆|−ν⊥and τ∼ |∆|−νk. The critical exponents β,ν⊥and νkcharacterize the steady state properties of the transition. It is also possible to look at the time dependent behavior at the critical point. Thus, for example, the particle density is observed to decay in time as ρ(t)∼t−θ. Different quantities can also be defined to evaluate the time properties of spreading experiments, in which the dynamics evolves starting from a single particle. In this case we can define the survival probability, S(t), as the probability that the activity lasts longer that t, finding at the critical point S(t)∼t−δ. These and other critical exponents are not independent, but are related by a set of scaling and hyperscaling relations [22]. Thus it is possible to give a full characterization of the phase transition of CP in Euclidean lattices using only three exponents, that we can take to be (without lack of generality) β,ν⊥and νk. Below the critical dimension dc= 4, the exponents are nontrivial, and depend explicitly on d. For d > dc, the exponents take the classical MF values β=νk= 1, ν⊥= 1/2. IV. HETEROGENEOUS MEAN-FIELD THEORY FOR THE CP Heterogeneous mean-field theory (HMF) is the basic starting point to obtain an analytical understanding of the behavior of any dynamical process on a complex network [1]. In order to take into account the possible fluctuations induced by the network connectivity, the partial densities ρk(t) of occupied vertices of degree k[8, 41] are considered, from which the total density of particles is obtained as ρ(t) = Pkρk(t)P(k). In the spirit of standard mean-field theories [42], the fact that the quantities ρk(t) are, in finite networks, of stochastic nature, is neglected. Instead, deterministic rate equations are considered, taking into account the changes in time of the partial densities, due to the different steps that the evolution of the model can take. In the case of the CP, the quantities ρk(t), given by ρk(t) = nk(t) NP(k),(8) where nk(t) is the number of particles in vertices of degree k, can be interpreted equivalently as the relative densities of particles in vertices of degree k, or the probabilities that a given vertex of degree kcontains a particle. In a step of the CP dynamical evolution, the partial density ρk(t) can decrease due to the annihilation of a particle in a vertex k(with rate 1), or can increase by the generation of an offspring in a vertex k′, nearest neighbor of k(with rate λ/k′). Therefore, the rate equations for the partial densities in a network characterized by a degree distribution P(k) and degree correlations given by the conditional probability P(k′|k) can be written as [19] ∂ρk(t) ∂t =−ρk(t) + λk[1 −ρk(t)] X k′ P(k′|k)ρk′(t) k′.(9) Given Eq. (9), ρk= 0 is always a solution. The conditions for the presence of non-zero steady states can be obtained by performing a linear stability analysis [43]. Neglecting 5 higher order terms, Eq. (9) becomes ∂ρk(t) ∂t ≃X k′ Lkk′ρk′(t)≡X k′−δk,k′+λk P(k′|k) k′ρk′(t). (10) It is easy to see that the Jacobian matrix Lkk′has a unique eigenvector vk=kand a unique eigenvalue Λ = λ−1. Therefore, a nonzero steady state is only possible for Λ >0, which translates in a critical threshold for the absorbing state phase transition λc= 1,(11) independent of the degree distribution and the correlation pattern. To get more detailed information on the process, and in particular on the shape of the order parameter as a function of the rate λ, we restrict our attention to uncorrelated networks. In this case, Eq. (9) reads ∂ρk(t) ∂t =−ρk(t) + λk hki[1 −ρk(t)]ρ(t).(12) Imposing the steady state condition, ∂tρk(t) = 0, yields the nonzero solutions ρk=λkρ/hki 1 + λkρ/hki(13) where ρkis now independent of time. By combining Eq. (13) with the definition of ρ, one obtains the selfconsistent equation for the order parameter ρ, ρ=λρ hkiX k kP (k) 1 + λkρ/hki,(14) that depends on the full degree distribution. In the case of SF networks, for which the degree distribution in the continuous degree approximation is given by P(k) = (γ−1)mγ−1k−γ, with mthe minimum degree in the network, the solution will depend on the degree exponent γ. Substituting the summation by an integral in Eq. (14), we obtain in the infinite network size limit (i.e. when the degree belongs to the range [m, ∞]) the expression ρ=F1, γ −1, γ, −hki λρm,(15) where F[a, b, c, z] is the Gauss hypergeometric function [40]. To evaluate the critical behavior for small ρ, we invert this expression using the asymptotic expansion of the hypergeometric function for low densities [40], obtaining the result ρ(λ)∼(λ−1)β, with β= 1/(γ−2) for 2 < γ < 3 and β= 1 for γ > 3, presenting additional logarithmic corrections at γ= 3. Right at the critical point, λ= 1, the particle density is expected to decay as a power law of time, ρ(t)∼t−θ [22], defining a new, temporal, critical exponent. This exponent can be estimated within HMF, by considering the time evolution of the total density at λ= 1, namely, ∂ρ(t) ∂t =X k P(k)∂ρk(t) ∂t =−ρ(t) hkiX k kρk(t).(16) To close this equation, we use a quasi-static approximation [17], which can be justified in terms of an adiabatic approximation for the full Langevin theory for the CP (see Sec. VI B). In essence, we consider that, even at the critical point, where no steady-state is present, the partial densities relax to a quasi-stationary state, where they take the form given by Eq. (13). In this case, for SF networks in the continuous degree approximation, Eq. (16) will read ∂ρ(t) ∂t ≃ −ρ(t)F1, γ −2, γ −1,−hki ρ(t)m,(17) Using the asymptotic approximation for the hypergeometric function, valid for low density, we obtain a decay exponent in infinite networks given by ρ(t)∼t−θ, with θ=βfor all γ(logarithmic corrections being again present at γ= 3). V. FINITE-SIZE SCALING FOR THE CP IN COMPLEX NETWORKS The exponents obtained within HMF theory in the previous Section correspond to the thermodynamic limit of a SF network of infinite size. Checking their accuracy in numerical simulations becomes thus a nontrivial task, particularly close to the critical point, due to the effects of finite network sizes. Indeed, because of the small-world property, the number of neighbors that can be reached starting from a certain node grows exponentially or faster with the geodesic distance. This implies that, even for large networks, just a few steps are sufficient to probe the finiteness of the system. Moreover, in SF networks, local topological properties show very strong fluctuations, increasing with the size of the network. For general critical phenomena, the theory of finite-size scaling (FSS) [44] has successfully overcome this problem for processes taking place on regular lattices, allowing the detection of the signature of continuous phase transitions even in very small systems. For absorbing state phase transitions, FSS is based on the observation that, even below the critical point, the density of active sites in surviving runs ρsreaches a quasi-steady state whose average is a decreasing function of the system size, and that can be expressed as a homogeneous scaling function of both the system size and the distance to the critical point. In the case of networks, system size is replaced by the number of vertices, and the surviving density is assumed to fulfill the relation [45, 46] ρs(∆, N) = N−β/¯νf(∆N1/¯ν),(18) 6 where f(x) is a scaling function that behaves as f(x)∼ xβfor x→ ∞ and f(x)∼const for x→0. Mean-field theory for homogeneous networks predicts the exponents β= 1 and ¯ν= 2. For the case of SF, Ref. [20] proposed a phenomenological Langevin equation for the particle density, taking the form dρ(t) dt = ∆ρ(t)−bρ(t)2−dρ(t)γ−1+pρ(t)η(t),(19) where η(t) is an uncorrelated Gaussian noise. Assuming a scaling form for the surviving density given by Eq. (18), and by means of a droplet-excitation argument, the authors of [20] found that β= 1/(γ−2) and ¯ν= (γ−1)/(γ−2) for γ < 3, independent of the network cutoff, whenever kc(N)> N1/γ [23]. In Ref. [25], this issue was pursued by focusing on the FSS form of survival probability, which at the critical point, and in networks of size N, was assumed to be S(t, N) = t−δf(t/tc(N)).(20) The scaling function f(x) is constant for small values of the argument and cutoff exponentially for x≫1. tc(N) is a characteristic cutoff time that, according to standard mean-field FSS theory should scale as tc(N)∼N1/2for homogeneous networks [22]. By means of a mapping to a biased, one-dimensional random walk, the authors of [25] found δ= 1, while the characteristic time showed the form, for heterogeneous networks, tc(N)∼pN/g, where gis defined in Eq. (5), and is thus dependent on the degree cutoff. This surprising result, well confirmed by numerical simulations [25], is in strong disagreement with results of Refs. [20, 23], in which no cutoff dependence was claimed. In order to fully ascertain the correct FSS behavior of CP in SF networks, we go beyond mean-field and phenomenological theories and tackle the full problem, taking into account its implicit stochastic fluctuations (particularly important in the vicinity of a critical point) by means of a Langevin approach. This problem is considered in the next Section. VI. LANGEVIN APPROACH FOR THE CP ON NETWORKS A. Generic formalism To account for the stochastic fluctuations of the CP close to the critical point, we derive here a Langevin equation describing the concentration ρk(t) or, alternatively, the number of active sites of degree k,nk(t). Our derivation follows closely the method developed in [17]. We start by deriving exact equations for the microscopic dynamics (at the vertex level) of the process. Let σi(t) be a random binary variable taking value σi(t) = 1 if node iis occupied by a particle at time t and σi(t) = 0 otherwise. Thus, the state of the process at time tis completely determined by the state vector Σ(t) = {σ1(t), σ2(t),···, σN(t)}. Variables σi(t) can undergo only two types of transition events: 1. σi(t) = 1 →σi(t+dt) = 0: Vertex iwas occupied by a particle at time t, and the particle annihilated during the time interval [t, t +dt]. 2. σi(t) = 0 →σi(t+dt) = 1: Vertex iwas empty at time tand it received an offspring from an occupied nearest neighbor during the time interval [t, t +dt]. Assuming that the temporal occurrence of these events follows Poisson processes, the previous two events can be encoded into a single dynamical equation that describes the evolution of σi(t) after an increment of time dt as σi(t+dt) = σi(t)ζi(dt) + [1 −σi(t)]ηi(dt),(21) where ζi(dt) and ηi(dt) are dichotomous random variables taking values ζi(dt) =    0 with probability dt 1 with probability 1 −dt (22) and ηi(dt) =          1 with probability λdtX j aijσj(t)1 kj 0 with probability 1 −λdtX j aijσj(t)1 kj . (23) Eqs. (22) and (21) describe the annihilation of particles, while Eq. (23) corresponds to the creation from occupied nearest neighbors. The set of random variables {ζi(dt); i= 1,···, N}are statistically independent of each other and of the conjugate random variables {ηi(dt); i= 1,···, N}. On the other hand, variables {ηi(dt); i= 1,···, N}are not totally independent since they may involve common events inducing correlations among them. For example, imagine two empty vertices, A and B, each of degree 1, connected to the same occupied vertex C. Because of the CP dynamics, during a particle reproduction event at vertex C, the particle must choose only one of its neighbors to send the offspring. Therefore, if vertex A gets the offspring, vertex B cannot receive it and vice-versa, inducing thus correlations between the random variables ηA(dt) and ηB(dt). However, it is easy to see that these correlations are of order dt2and can be then safely neglected. In any case, this effect only exists in networks with a quenched topology. In contrast, the annealed network topology changes faster than the CP dynamics and, therefore, such correlations are absent. Equations (21), (22), and (23) describe the evolution of the state of the system at the most detailed possible level of description by specifying the precise state of each and every one of the vertices of the network. This description, 7 although exact, is not very useful to derive general properties of the system, which are better described by coarsegrained quantities. In heterogeneous random networks with given degree distribution P(k) and degree-degree correlations P(k′|k), the degree of vertices kis the most appropriate indicator of the different classes of vertices. Therefore, we consider all vertices with the same degree to be statistically equivalent. Following these ideas, let nk(t) be the number of active vertices of degree kat time t, that is, nk(t)≡X i∈k σi(t).(24) As we can see, nk(t) is the sum of a large number of random variables that are nearly statistically independent in the quenched version of the network and totally independent in its annealed version. Therefore, by invoking the central limit theorem, we expect this variable to follow a Gaussian distribution and, consequently, to follow a Langevin dynamics. To derive the specific form of this Langevin equation, we need to calculate the infinitesimal moments of nk(t), which can be done using Eqs. (21), (22), and (23). Using the results in Appendix A, we can finally write the corresponding Langevin equation for the CP on annealed networks, namely dnk(t) dt =−nk(t) + λ[1 −ρk(t)] X k′ P(k|k′)nk′(t) (25) +ξk(t)snk(t) + λ[1 −ρk(t)] X k′ P(k|k′)nk′(t), where ρk(t) = nk(t)/P (k)Nis the relative density of active vertices of degree kand {ξk(t), k = 1,···, kc}are Gaussian white noises (with zero mean and unit variance) uncorrelated among them. Eq. (25) implicitly assumes that nk(t) is a continuous variable. This approximation is reasonable as long as nk(t)≫1, which is usually the case in very large systems, when we consider steady-state properties. Equation (25) is one of the main results of this paper and is also the starting point for our subsequent analysis. As one immediately recognizes, the drift term in Eq. (25) corresponds to the standard mean-field approximation derived in Eq. (9). It is easy to see that the potential associated with this drift term has a stable minimum whenever λ > λc= 1 which does not depend on the particular correlation pattern given by P(k|k′). The position of this minimum corresponds to the steady solution in the active phase in the thermodynamic limit. The diffusion term, on the other hand, points to a process with multiplicative noise that, as we shall see, has important implications when the system is close to its critical point in finite size systems. B. Uncorrelated random networks Finding solutions of Eq. (25) for networks with general degree-degree correlations is a rather difficult task. In this paper, we focus on the simplest (but instructive) case of uncorrelated random networks with a given degree distribution P(k). For this class of networks, the transition probability takes the simple form P(k|k′) = kP(k)/hki which allows us to write Eq. (25) as dρk(t) dt =−ρk(t) + λk hki[1 −ρk(t)] ρ(t) (26) +s1 Nkρk(t) + λk hki[1 −ρk(t)] ρ(t)ξk(t), where ρ(t) = Pknk(t)/N =PkP(k)ρk(t) is the global concentration of active nodes at time tand we have divided Eq. (25) by the number of vertices of degree k, Nk=NP(k). Analogously, we can write a Langevin equation for ρ(t) as dρ(t) dt =ρ(t) ∆−λX k kP (k) hkiρk(t)!(27) +X k P(k)s1 Nkρk(t) + λk hki[1 −ρk(t)] ρ(t)ξk(t), where we have defined ∆ ≡λ−1 so that the critical point corresponds to ∆ = 0. Eq. (27) is not yet a closed equation for ρ(t) because both the drift and diffusion terms involve the partial densities ρk(t). To close it, we use an adiabatic approximation [47]. From Eq. (27) we know that close to the critical point, ∆ ≈0, ρ(t) is a slowly varying variable. This is due to the fact that the first term in the right hand side of Eq. (27) is of order higher than ρ. On the other hand, ρk(t) is a variable that relaxes exponentially fast to its quasi-equilibrium state since the lowest order in Eq. (26) is linear in ρk3. The adiabatic approximation consists in neglecting the term dρk(t)/dt in front of ρk(t) and assuming that ρkis a stochastic variable that evolves much faster than ρ(t). Thus, setting dρk(t)/dt = 0 in Eq. (26) and solving for ρk(t), we obtain ρk(t)≈λkρ(t) hki+λkρ(t)(28) +hki hki+λkρ(t)s1 Nkρk(t) + λk hki[1 −ρk(t)] ρ(t)ξk(t). 3It is worth mentioning that this separation of time scales between the partial quantities ρkand the global one ρhas also been observed in other dynamics like the A+A−→ ∅ diffusionannihilation process [17] or the voter models [48]. 8 The noise term in this equation is subdominant due to its dependence on the size of the system. Thus, replacing the dominant term in the diffusion one, we finally obtain ρk(t)≈λkρ(t) hki+λkρ(t)+s1 Nk 2λkρ(t)hki2 [hki+λkρ(t)]3ξk(t).(29) In this way, we obtain an expression for the partial densities ρkas a function of kand ρ(t) only. Replacing this expression in Eq. (27) and keeping only the first order in N−1 k, we obtain dρ(t) dt =ρ(t) ∆−λX k kP (k) hki λkρ(t) hki+λkρ(t)!(30) +X k P(k)s1 Nk 2λkρ(t)hki2 [hki+λkρ(t)]3ξk(t). Notice that the sum of statistically independent Gaussian white noises is another Gaussian white noise whose variance is the sum of the individual variances. Thus, the diffusion term in the last equation is, indeed, a Gaussian white noise. Therefore, we can finally write dρ(t) dt =ρ(t) (∆ −λΘ[ρ(t)]) + r2λρ(t) NΛ[ρ(t)]ξ(t),(31) where Θ[ρ(t)] ≡X k kP (k) hki λkρ(t) hki+λkρ(t)(32) and Λ[ρ(t)] ≡X k kP (k) hkihki3 [hki+λkρ(t)]3.(33) Eq. (31) is now a closed equation for the total density of active vertices ρwhich must be solved with an absorbing boundary at ρ= 0 and a reflecting one at ρ= 1. As we can see from Eq. (31), there is an explicit dependence on the size of the network Nin the diffusion term of the Langevin equation. This size dependence, together with the specific functional forms of Θ[ρ] and Λ[ρ] will determine the finite size behavior of the system near the critical point. VII. CP IN ANNEALED SCALE-FREE NETWORKS In this section we focus on heterogeneous networks with a power law degree distribution P(k)∼k−γwith k∈[m, M], where 2 < γ < 3, and M=N1/ω is the degree upper cutoff. In particular, we consider the case ω≥γ−1, so that the average maximum of the degree distribution kc∼M=N1/ω is a hard cutoff (Sec. II). The case ω < γ −1 will be considered in Sec. IX. Given the form of the degree distribution it is possible to evaluate explicitly the functional form of Θ[ρ], that determines the dynamical properties of the CP. From its definition, Eq. (32), it is easy to see that in the limit of small density Θ[ρ] has two different functional forms depending on whether ρis larger or smaller than the quantity hki/λkc,kcbeing the network cutoff. Thus we have (see Appendix B): Θ[ρ] =        gλρ ρ ≪hki λkcregion II C(γ)λρ hkiγ−2hki λkc≪ρ≪1 region I , (34) where g=hk2i/hki2and C(γ) = mγ−2Γ(γ−1)Γ(3 −γ).(35) We denote the regime for small ρas region II, and the regime for larger ρas region I. Analogously, we can evaluate the behavior of the function Λ as Λ[ρ] =      1 region II 1−e C(γ)λρ hkiγ−2region I ,(36) with e C(γ) = 1 2mγ−2(γ−2)(γ−1)Γ(1 + γ)Γ(1 −γ).(37) As we can see, the correction term in the region I is always very small as compared to 1. Therefore, in the rest of the paper we consider that Λ[ρ] = 1. In Fig. 1 the behavior of Θ[ρ] is evaluated by numerically performing the summation in Eq. (32). The linear behavior for small densities is very well obeyed. Instead, the scaling of the region I is not cleanly observed even for the largest network considered (N= 108). This is due to the fact that region I is surrounded by two slow crossovers, one for ρ≈ hki/kc(where the transition between region I and II takes place) and the other for ρ= 1 (where Θ becomes independent of ρ). This has the consequence that some of the theoretical predictions made using the simple approximation given in Eq. (34) are difficult to observe except for extremely large system sizes. As a consequence of the form of Θ[ρ], the behavior of the system at criticality strongly depends on the type of experiment performed to probe the absorbing transition, see Fig. 2. Indeed, experiments with stochastic trajectories exploring the region ρ≫ hki/λkcfeel a drift of the form Ψ[ρ]≃ρ"∆−C(γ)λρ hkiγ−2#type I drift.(38) Instead, any experiment such that trajectories mainly stay in the region ρ≪ hki/λkcfeels a drift term of the form Ψ[ρ] = ρ(∆ −λ2gρ) type II drift.(39) 9 10-4 10-3 10-2 10-1 100101102103104105 kc ρ /<k> 10-6 10-5 10-4 10-3 10-2 10-1 100 Θ[ρ] N=105 N=106 N=107 N=108 M=N1/(γ−1) Type II Type I ρ ρ0.5 FIG. 1: Numerical evaluation of the function Θ[ρ] as a function of kcρ/hkifor different network sizes. The degree exponent is γ= 2.5 and we use ω=γ−1. In the type II region, it is clearly visible a linear behavior, in agreement with Eq. (34). In the type I region, convergence towards the theoretical expression given by Eq. (34) is much slower. 100101102103104 t 10-8 10-6 10-4 10-2 100 ρ(t) M=N1/(γ−1) 100101102103104 t 10-8 10-6 10-4 10-2 100 M=N1/2 t-2 t-1 FIG. 2: Example of two different types of experiment to study the CP dynamics at criticality performed in an annealed network of size N= 108,γ= 2.5, and m= 2. The top (black) curve is the density decay starting from a fully active network. The bottom (blue) curve corresponds to a spreading experiment starting from a single active vertex. The left plot corresponds to the hard cutoff M=N1/(γ−1) and the right plot to M=N1/2. Grey areas depict the values of ρin the domain ρ∈[N−1,hkik−1 c]. In all cases trajectories are for a single run in an instance network. The quantity g, that will play a fundamental role in the rest of the paper, diverges with the cutoff kcas kγ−3 c for γ < 3 (and ω > γ −1). Notice that for γ > 3 the leading order is linear in both cases. It is also worth stressing that if one lets kcdiverge, regime II disappears and one is left only with regime I, that coincides with what is found using HMF on infinite networks (Sec. IV). However, in any finite network, when the density gets small it is the drift of type II that rules the dynamics. Based on the explicit expression of Θ[ρ] we now provide a qualitative and quantitative description of the three types of experiment that explore the critical properties of the CP dynamics: A) density decay at criticality, B) spreading experiments and C) surviving runs. A. Density decay at criticality Starting from a configuration full of active vertices at t= 0, the concentration of active vertices is monitored as a function of time until the trajectory is trapped at the absorbing boundary. Then an average is performed over a large number of different runs up to a time such that all runs have survived. In this case, after a initial timescale t×, the system first feels the type I drift and after a crossover time t∗, at very low concentrations, the type II one. Inserting Eq. (38) into Eq. (31), a pure drift of the type I predicts a behavior ρI(t)∼t−θ with θ= 1/(γ−2). Inserting Eq. (39) gives instead ρII (t)∼(gt)−1. The crossover between the two types of behavior occurs for a time t∗such that ρII(t∗)≃ hki/kc, i.e. t∗∼kc/(ghki)∼ hkikγ−2 c. A third time scale defines the survival time of the different runs that, as we will see in the next subsection, scales as tc∼pN/g. Fig. 2 shows simulation results for this type of experiment (top curves) in annealed networks with γ= 2.5, m= 2, for ω= 2 or ω=γ−1 for a single run starting from a fully active network. Different colors (grey and white) indicate the different regions depending on the shape of the drift term. The first thing to notice is that, in the case of ω=γ−1, the region that corresponds to the type I drift is wider as compared to the case ω= 2. Nevertheless, even in this optimal case, we do not observe cleanly the two different values of the exponent θ. For comparison purposes, in Fig. 2 we also plot functions t−2and t−1that would correspond to the pure type I and II behaviors for γ= 2.5. The exponent θapproaches but does not reach the theoretical value θ= 2 even though simulations are performed in networks of size N= 108. The situation in the case of ω= 2 is even worse because the crossover happens at shorter times and the value θ= 2 is even more difficult to observe. The same is observed in Fig. 3, where we show the same as Fig. 2 for networks with γ= 2.5, ω= 2 and different network sizes but averaging over 100 runs of the process over the same instance network. Additional information is provided by the local effective exponent of the temporal decay as a function of time (Fig. 4). The effective exponent decreases initially quite fast, but even for N= 108the crossover to regime II takes place well before the asymptotic value θ= 1/(γ−2) = 2 is reached. Eventually the effective exponent sets to a constant value, that is, quite surprisingly, close to 1.2 instead of the expected value θ= 1. The reason why we do not see convincing numerical evidence of any of the two scaling exponents expected from the theory is that the separation of time scales between 16 To derive the previous equation, we have taken into account that, since σi(t) are binary variables taking only values 0 or 1, σ2 i(t) = σi(t) and that σi(t)[1 −σi(t)] = 0; ∀t. Terms of order dt2have also been neglected. Under the annealed approximation, we replace in Eqs. (A5) and (A6) the adjacency matrix aij by its average value, Eq. (2), which allows us to carry out the sums in Eqs. (A5) and (A6) and, finally, to obtain the Langevin equation Eq. (25). APPENDIX B: CALCULATION OF Θ[ρ]IN UNCORRELATED SF NETWORKS Let us consider the definition of Θ[ρ] in Eq. (32), namely Θ[ρ] = X k kP (k) hki λkρ/hki 1 + λkρ/hki.(B1) The evaluation of this quantity in finite networks depends on the value of the particle density ρ. In particular, if ρ≪ hki/λkc, where kcis the network cutoff, then the denominator in Eq. (B1) can be approximated by unity, and we have Θ[ρ]≃X k kP (k) hki λkρ hki=hk2i hki λρ hki=gλρ, (B2) where g=hk2i/hki2. On the other hand, outside this region we must keep the full denominator in Eq. (B1). To estimate Θ in this case, we perform a continuous degree approximation, that is, Θ[ρ] = Zkc m kP (k) hki λkρ/hki 1 + λkρ/hki =F»1, γ −2, γ −1,−hki λρm –−F»1, γ −2, γ −1,−hki λρkc–„m kc«γ−2 , (B3) where F[a, b, c, z] is the Gauss hypergeometric function. Using the asymptotic expansions of the hypergeometric function for small and large arguments [40], we can estimate the value of Θ in the domain hki/λkc≪ρ≪1. Within such domain, the second term in Eq. (B3) becomes an asymptotically small constant as compared to the first term, that yields Θ[ρ]≃Γ(γ−1)Γ(3 −γ)λρm hkiγ−2 .(B4) APPENDIX C: INITIAL TIME SCALE To compute the initial time scale t×needed to reach region I starting from an arbitrary initial condition ρ0at criticality, we consider Eq. (31) with the drift term given by Eq. (38) and ∆ = 0, namely dρ(t) dt =−C(γ) hkiγ−2ρ(t)γ−1.(C1) The solution of this equation is ρ(t) = ρ2−γ 0+(γ−2)C(γ) hkiγ−2t−1/(γ−2) .(C2) The asymptotic state ρ(t)∼t−1/(γ−2), independent of the initial condition, is reached for times t, such that (γ−2)C(γ) hkiγ−2t≫ρ2−γ 0,(C3) that is, for t > t×, with t×=ρ2−γ 0hkiγ−2 (γ−2)C(γ)=ρ2−γ 0(γ−1)γ−2 (γ−2)γ−1Γ(3 −γ)Γ(γ−1), (C4) where we have used the definition of C(γ) in Eq. (35) and hki= (γ−1)m/(γ−2). APPENDIX D: SURVIVAL PROBABILITY EQUATION Using standard techniques of stochastic processes theory, we can obtain the partial differential equation satisfied by the survival probability of the CP dynamics at criticality, starting from an initial concentration ρ0, S(t|ρ0), namely [47] ∂S(t|ρ0) ∂t =−ρ0Θ[ρ0]∂S(t|ρ0) ∂ρ0 +ρ0 N ∂2S(t|ρ0) ∂ρ2 0 .(D1) This equation is the result of integrating the backwards Fokker-Planck equation in the domain ρ∈[0,1] and it should be solved with the initial condition S(t= 0|ρ0) = 1 and boundary conditions S(t|ρ0= 0) = 0 and ∂S(t|ρ0) ∂ρ0ρ0=1 = 0 (D2) that correspond to an absorbing boundary at ρ= 0 and a reflecting one at ρ= 1. The survival probability Eq. (20) can then be evaluated as S(t) = S(t|ρ0= 1/N).(D3) We first start by evaluating the exponent δ. To this end, it is only necessary to solve the problem in the thermodynamic limit N→ ∞. However, the above formulation is not the most appropriate for this purpose, since the solution must be evaluated at ρ0=N−1, that is, a value that depends on the size of the system. Therefore, we perform the change of variables n0=Nρ0(D4) where n0is the initial number of active vertices, which is eventually set to n0= 1 and, therefore, is independent 17 of the system size. Using this new variable, Eq. (D1) becomes ∂S(t|n0) ∂t =−n0Θhn0 Ni∂S(t|n0) ∂n0 +n0 ∂2S(t|n0) ∂n2 0 .(D5) Notice that now the limit N→ ∞ can be taken in Eq. (D5). In this limit, the first term in the right hand side of Eq. (D5) vanishes and the process becomes a purely diffusive one with multiplicative noise. The solution is S(t|n0) = 1 −e−n0/t ≈n0 t.(D6) Setting finally n0= 1, leads to the exponent δ= 1 for any γ. To evaluate the cutoff tc(N), we compute the second moment of the survival times, T2(n0), starting from n0 active sites. However, to compute T2we first need to compute the average surviving time T1(n0). It is easy to see that T1(n0) = R∞ 0S(t|n0). Using this result in Eq.(D5), and assuming that trajectories never feel the type I drift, yields the following differential equation for T1(n0) d2T1(n0) dn2 0−g Nn0 dT1(n0) dn0 =−1 n0 .(D7) with boundary conditions T1(0) = 0 and T′ 1(N) = 0. The solution of this problem is T1(n0) = s2N gZn0√g 2N 0 dueu2Z√gN/2 u dt te−t2.(D8) When Nis very large, the upper limit in the first integral becomes very small. Therefore, we take the limit of the integrand when uis close to zero, that is, T1(n0)≃s2N gZn0√g 2N 0 du[1+u2+···][−ln u+γ+···] (D9) which finally leads to T1(n0)≃ −n0ln n0rg 2N.(D10) Similarly, the differential equation for T2(n0) can also be obtained from Eq.(D5) as T2(n0) = 2 Z∞ 0 tS(t|n0)dt. (D11) This results in the following differential equation (involving also T1(n0)) d2T2(n0) dn2 0−g Nn0 dT2(n0) dn0 =−2T1(n0) n0 .(D12) that satisfies the same boundary conditions as T1(n0). The solution of this equation is T2(n0) = 2N gZn0√g 2N 0 dueu2Z∞ u G(t)e−t2dt (D13) where G(t) = 2 tZt 0 dueu2Z∞ u dq qe−q2(D14) In the limit of large N, this expression can be approximated as T2(n0) = n0s2N gZ∞ 0 e−t2G(t)dt, (D15) proving then Eq. (43). APPENDIX E: PROBABILITY DENSITY FUNCTION FOR SURVIVING RUNS At the critical point, the probability density p(n, t|n0) of the number of active vertices at time tgiven that the process had n0active ones at time t= 0 is ruled by a Fokker-Planck equation with a drift term Ψ(n) = −nΘ[n/N] and a diffusion coefficient D(n) = 2n, that is, ∂ ∂n hnΘhn Nip(n, t|n0)i+∂2 ∂n2[np(n, t|n0)] = ∂p(n, t|n0) ∂t . (E1) A direct substitution of Eq. (46) into Eq. (E1) leads to ∂ ∂n hnΘhn Nips(n, t|n0)i+∂2 ∂n2[nps(n, t|n0)] = ∂ps(n, t|n0) ∂t +ps(n, t|n0)dln [S(t|n0)] dt .(E2) The density ps(n, t|n0) has, by construction, a welldefined steady state, that we denote by ps(n)≡lim t≫1ps(n, t|n0),(E3) which is independent of the initial condition. By taking the limit t≫1 in Eq. (E2), we obtain ∂ ∂n hnΘhn Nips(n)i+∂2 ∂n2[nps(n)] = κps(n),(E4) 18 where κ= lim t≫1 d dt ln [S(t|n0)] (E5) Using the result given in Eq. (D6) we conclude that κ= 0, meaning that ps(n) satisfies the potential solution of the Fokker-Planck equation [47]. We can then write that ρs(0, N)∝1 NZN 1 e−φ(n,N)dn, (E6) with the effective potential φ(n, N) = ZΘhn Nidn. (E7) As in the case of the function Θ, the potential φ(n, N) takes a different functional form depending on the value of n. Direct integration of Eq. (34) gives φ(n, N) =        λg 2Nn2n≪hkiN λkc C(γ) γ−1λ hkiNγ−2nγ−1hkiN λkc≪n≪N (E8) At the critical point, λ= 1, we can use this result to write ρs(0, N)∝1 N(ZhkiN kc 1 e−φII (n,N)dn +ZN hkiN kc e−φI(n,N)dn), (E9) where subindices I and II refer to which type of potential is dominating the integral. In the limit N≫1, we can evaluate the contribution of each integral, leading to Eq. (48). [1] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev. Mod. Phys. 80, 1275 (2008). [2] R. Albert and A.-L. Barab´asi, Rev. Mod. Phys. 74, 47 (2002). [3] S. N. Dorogovtsev and J. F. F. Mendes, Evolution of networks: From biological nets to the Internet and WWW (Oxford University Press, Oxford, 2003). [4] D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998). [5] A.-L. Barab´asi and R. Albert, Science 286, 509 (1999). [6] R. Cohen and S. Havlin, Phys. Rev. Lett. 90, 058701 (2003). [7] R. Pastor-Satorras and A. Vespignani, Evolution and Structure of the Internet. A Statistical Physics Approach (Cambridge University Press, Cambridge, 2004). [8] R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001). [9] A. L. Lloyd and R. M. May, Science 292, 13161317 (2001). [10] R. Cohen, K. Erez, D. ben Avraham, and S. Havlin, Phys. Rev. Lett. 86, 3682 (2001). [11] D. S. Callaway, M. E. J. Newman, S. H. Strogatz, and D. J. Watts, Phys. Rev. Lett. 85, 5468 (2000). [12] R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001). [13] M. Leone, A. V´azquez, A. Vespignani, and R. Zecchina, Eur. Phys. J. B 28, 191 (2002). [14] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Phys. Rev. E 66, 016104 (2002). [15] S. Dorogovtsev, A. Goltsev, and J. Mendes, Eur. Phys. J. B 38, 177 (2004). [16] L. K. Gallos and P. Argyrakis, Phys. Rev. Lett. 92, 138301 (2004). [17] M. Catanzaro, M. Bogu˜n´a, and R. Pastor-Satorras, Phys. Rev. E 71, 056104 (2005). [18] V. Colizza, R. Pastor-Satorras, and A. Vespignani, Nature Physics 3, 276 (2007). [19] C. Castellano and R. Pastor-Satorras, Phys. Rev. Lett. 96, 038701 (2006). [20] H. Hong, M. Ha, and H. Park, Phys. Rev. Lett. 98, 258701 (2007). [21] N. Goldenfeld, Lecture notes on phase transitions and the renormalization group, Frontiers in Physics (AddisonWesley, Massachusetts, 1992). [22] J. Marro and R. Dickman, Nonequilibrium phase transitions in lattice models (Cambridge University Press, Cambridge, 1999). [23] M. Ha, H. Hong, and H. Park, Phys. Rev. Lett. 98, 029801 (2007). [24] C. Castellano and R. Pastor-Satorras, Phys. Rev. Lett. 98, 029802 (2007). [25] C. Castellano and R. Pastor-Satorras, Phys. Rev. Lett. 100, 148701 (2008). [26] R. Pastor-Satorras, A. V´azquez, and A. Vespignani, Phys. Rev. Lett. 87, 258701 (2001). [27] M. A. Serrano, M. Bogu˜n´a, R. Pastor-Satorras, and A. Vespignani, in Large scale structure and dynamics of complex networks: From information technology to finance and natural sciences, edited by G. Caldarelli and A. Vespignani (World Scientific, Singapore, 2007), pp. 35–66. [28] A. Bekessy, P. Bekessy, and J. Komlos, Stud. Sci. Math. Hungar. 7, 343 (1972). [29] E. A. Bender and E. R. Canfield, Journal of Combinatorial Theory A 24, 296 (1978). [30] B. Bollob´as, Eur. J. Comb. 1, 311 (1980). [31] M. Molloy and B. Reed, Random Struct. Algorithms 6, 161 (1995). [32] M. Catanzaro, M. Bogu˜n´a, and R. Pastor-Satorras, Phys. Rev. E 71, 027103 (2005). [33] M. Bogu˜n´a and R. Pastor-Satorras, Phys. Rev. E 68, 036112 (2003). [34] S. N. Dorogovtsev, J. F. F. Mendes, and A. N. Samukhin, Phys. Rev. Lett. 85, 4633 (2000). [35] S. Gil and D. Zanette, Eur. Phys. J. B 47, 265 (2005). [36] D. Stauffer and M. Sahimi, Phys. Rev. E 72, 46128 (2005). [37] S. Weber and M. Porto, Phys. Rev. E 76, 046111 (2007). [38] S. N. Dorogovtsev and J. F. F. Mendes, Advances in Physics 51, 1079 (2002). [39] M. Bogu˜n´a, R. Pastor-Satorras, and A. Vespignani, Eu- 19 ropean Physical Journal B 38, 205 (2004). [40] M. Abramowitz and I. A. Stegun, Handbook of mathematical functions. (Dover, New York, 1972). [41] M. Bogu˜n´a, R. Pastor-Satorras, and A. Vespignani, in Statistical Mechanics of Complex Networks, edited by R. Pastor-Satorras, J. M. Rub´ı, and A. D´ıaz-Guilera (Springer Verlag, Berlin, 2003), vol. 625 of Lecture Notes in Physics. [42] H. E. Stanley, Introduction to phase transitions and critical phenomena (Oxford University Press, Oxford, 1971). [43] M. Bogu˜n´a and R. Pastor-Satorras, Phys. Rev. E 66, 047104 (2002). [44] V. Privman, Finite Size Scaling and Numerical Simulation of Statistical Systems (World Scientific, Singapore, 1990). [45] D. H. Zanette, Phys. Rev. E 64, 050901 (2001). [46] P. R. A. Campos, V. M. de Oliveira, and F. G. B. Moreira, Phys. Rev. E 67, 026104 (2003). [47] G. W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences (SpringerVerlag, Berlin Heidelberg New York, 2004). [48] V. Sood and S. Redner, Phys. Rev. Lett. 94, 178701 (2005).