Diffusion-annihilation processes in complex networks

We present a detailed analytical study of the $A+A\to\emptyset$ diffusion-annihilation process in complex networks. By means of microscopic arguments, we derive a set of rate equations for the density of $A$ particles in vertices of a given degree, valid for any generic degree distribution, and which we solve for uncorrelated networks. For homogeneous networks (with bounded fluctuations), we recover the standard mean-field solution, i.e. a particle density decreasing as the inverse of time. For heterogeneous (scale-free networks) in the infinite network size limit, we obtain instead a density decreasing as a power-law, with an exponent depending on the degree distribution. We also analyze the role of finite size effects, showing that any finite scale-free network leads to the mean-field behavior, with a prefactor depending on the network size. We check our analytical predictions with extensive numerical simulations on homogeneous networks with Poisson degree distribution and scale-free networks with different degree exponents.

We present a detailed analytical study of the A + A → ∅ diffusion-annihilation process in complex networks. By means of microscopic arguments, we derive a set of rate equations for the density of A particles in vertices of a given degree, valid for any generic degree distribution, and which we solve for uncorrelated networks. For homogeneous networks (with bounded fluctuations), we recover the standard mean-field solution, i.e. a particle density decreasing as the inverse of time. For heterogeneous (scale-free networks) in the infinite network size limit, we obtain instead a density decreasing as a power-law, with an exponent depending on the degree distribution. We also analyze the role of finite size effects, showing that any finite scale-free network leads to the mean-field behavior, with a prefactor depending on the network size. We check our analytical predictions with extensive numerical simulations on homogeneous networks with Poisson degree distribution and scale-free networks with different degree exponents.

I. INTRODUCTION
In the last few years, complex networks have become a new paradigm for complexity. The merging of graph theory together with new and classical statistical physics tools has lead to the development of a modern theory of complex network [1,2], that has found fruitful applications in domains as divers as technology (the physical Internet [3], the World-Wide Web [4], power grids [5]), biology (protein-protein interaction networks [6], metabolic networks [7], foodwebs [8,9]), social sciences (sexual contact networks [10], friendship networks [11], scientific collaboration networks [12,13]), etc.
The statistical analysis of many "real-world" networks has shown that most of these systems seem to share some typical features, the most relevant of them being the small-world property [14] and a large connectivity heterogeneity, reflected in the presence of a scale-free degree distribution [15]. The small-world property refers to the fact that, in real networks, the hop distance between two randomly chosen elements of the system is very small if compared to the total number of elements. More precisely, if ℓ is the average distance between two elements (or vertices), measured as the smallest number of connections (or edges) between any vertices, and N is the system size (number of vertices) then, usually ℓ increases logarithmically or slower with N . On the other hand, scale-free networks are characterized by a degree distribution P (k), defined as the probability that a randomly selected vertex is connected to k other vertices (has degree k), that decreases as a power-law, where γ is a characteristic degree exponent, usually in the range 2 < γ < 3. For these values of the degree exponent, the fluctuations in the degree distribution, measured by the second moment k 2 , diverge in the infinite network size limit, N → ∞, giving rise a very heterogeneous connectivity lacking any characteristic degree scale. This behavior is in opposition to the more classical homogeneous networks [16], which have a degree distribution decaying exponentially or faster and exhibit bounded degree fluctuations.
Given that complex networks are widespread in nature, it becomes a quite interesting issue the characterization of the effects that their complex topology can have on dynamical processes taking place on top of these systems. For example, it has been shown that heterogeneous networks are remarkably weak in front of targeted attacks, aimed at destroying the most connected vertices [17,18], as well as to the propagation of infective agents [19,20]. These properties, which are mainly due to the critical interplay between topology and dynamics in heterogeneous networks, are otherwise absent in their homogeneous counterparts.
Epidemic processes, chemical reactions, and many other dynamic processes, can all be modeled in terms of reaction-diffusion processes [21]. These are dynamic systems that involve particles of different "species" (A i , i = 1, . . . n) that diffuse stochastically and interact among them following a fixed set of reaction rules. One's interest is usually focused on the time evolution and steady states of the densities of the different species ρ Ai (t), and the possible presence of critical phase transitions [22]. While much is known about the behavior of reactiondiffusion processes on regular homogeneous lattices, the situation is not so well established in what respects the possible effects that a heterogeneous connectivity structure can have on them. At this respect, it is interesting the work presented in Ref. [23], in which a numerical simulation analysis of the diffusion-annihilation process A + A → ∅ [24], was performed on scale-free networks. In this reaction-diffusion process, particles of a single species A diffuse on the vertices of a network and annihilate upon contact (when two A particles fall on the same vertex).
In regular lattices of Euclidean dimension d, it is well known that the local density of A particles, ρ(x, t), is ruled by a Langevin equation [25], where η(x, t) is an uncorrelated Gaussian noise. Dynamical renormalization group arguments allow to show that the average density of A particles, ρ(t) = ρ(x, t) , behaves in the large time limit as where ρ 0 is the initial particle density, and the exponent α takes the values α = d/d c for d ≤ d c and α = 1 for d > d c , where d c = 2 is the critical dimension of this process. For d > d c one thus recovers the mean-field solution, obtained from Eq. (2) by setting the diffusion coefficient D and the noise term η(x, t) equal to zero. The numerical simulations of the A + A → ∅ diffusionannihilation process reported in Ref. [23], performed in scale-free networks with general degree exponent γ, generated using the configuration model [26,27,28,29], led the authors to conclude that the behavior in time of the average density of A particles can be approximated by Eq. (3), where the asymptotic exponent α is a decreasing function of γ and is surprisingly larger than 1 for γ < 3. The authors attributed this effect to the smallworld nature of the networks, and to the existence of hubs (vertices with a large number of connections).
In spite of the potential interest of this result, no theoretical arguments have been proposed so far to back up the numerical conclusions reached in [23]. In this paper we tackle this task, by developing a mean-field analysis of the A + A → ∅ process. This analysis, made in the continuous k approximation and inspired in previous works made for epidemic spreading [19,30,31], results in a set of differential equations for the density of A particles in the vertices of degree k, which are valid for networks with arbitrary degree distribution P (k) and two vertex correlations [32], determined by the conditional probability P (k ′ |k) that a vertex of degree k is connected to a vertex of degree k ′ [33,34]. The solution of these equations for the particular case of uncorrelated networks (in which the conditional probability P (k ′ |k) is independent of k) shows that, while homogeneous networks display a pure mean-field behavior with exponent α = 1, scale-free networks with γ < 3 in the infinite size limit exhibit instead an exponent depending on the properties of the network, i.e. α = 1/(γ − 2). Remarkably, this solution in the infinite size limit shows a crossover for any finite network to a linear behavior 1/ρ(t) ∼ t, with a slope depending on the network size. Our analytical results are confirmed by means of large scale numerical simulations for both homogeneous and heterogeneous networks.
We have organized the present paper as follows: In Sec. II we derive, from microscopic considerations, the mean-field differential equations for the A + A → ∅ diffusion-annihilation process in complex random networks with arbitrary degree distribution and two vertex correlations, quantified by means of P (k) and P (k ′ |k), respectively. We consider the case of absence of correlations, finding the density of A particles for general homogeneous networks. Sec. III is devoted to explicit results for scale-free networks, both in the infinite size limit and for finite size networks. In Sec. IV, our analytic results are compared with extensive numerical simulations of the diffusion-annihilation process running on top of homogeneous and heterogeneous (scale-free) networks. Finally, our conclusions are presented in Sec. V.

II. THE A + A → ∅ REACTION IN COMPLEX NETWORKS
Let us consider the diffusion-annihilation process A + A → ∅ on a complex network of size N which is fully defined by the adjacency matrix a ij , that takes the values a ij = 1 if vertices i and j are connected by an edge, and 0 otherwise. From a statistical point of view, the network can also be characterized by its degree distribution P (k) and its degree correlations, given by the conditional probability P (k ′ |k). Each vertex in the network can host at most one A particle, and the dynamics of the process is defined as follows: Each particle jumps at a certain rate λ to a randomly chosen nearest neighbor. If it is empty, the particle fills it, leaving the first vertex empty. If the nearest neighbor is occupied, the two particles annihilate, leaving both vertices empty.
In order to study analytically this process in a general complex network, in which vertices can show large degree fluctuations, we are forced to consider the partial densities ρ k (t), representing the density of A particles in vertices of degree k, or, in other works, the probability that a vertex of degree k contains an A particle at time t [19,30]. From these partial densities, the total density of A particles is recovered from While it is possible to obtain a rate equation for the densities ρ k (t) by means of intuitive arguments [19,31], in the following we will pursue a more microscopical approach, which can be generalized to tackle other kinds of problems. Let n i (t) be a dichotomous random variable taking values 0 or 1 whenever vertex i is empty or occupied by an A particle, respectively. Using this formulation, the state of the system at time t is completely defined by the state vector n(t) = {n 1 (t), n 2 (t), · · · , n N (t)}. Assuming that the time evolution of particles follows a Poisson process [21], the evolution of n(t) after a time increment dt can be expressed as where η(dt) and ξ(dt) are dichotomous random variables taking values and where a ij is the adjacency matrix and λ is the jumping rate that, without loss of generality, we set equal to 1. The first term in Eq. (5) stands for an event in which vertex i is occupied by a particle and, during the time interval (t, t + dt), it becomes empty, either because the particle in it decides to move to another vertex or because a particle in a nearest neighbor of i decides to jump to i, annihilating thus both particles. The second term corresponds to the case in which vertex i is empty and a particle in a neighbor vertex of i decides to move to that vertex [45]. Taking the average of Eq. (5), we obtain equation that describes the average evolution of the system, conditioned to the knowledge of its state at the previous time step. Then, after multiplying Eq. (8) by the probability to find the system at state n at time t, and summing for all possible configurations, we are led to where we have introduced the notation ρ i (t) ≡ n i (t) and The derivation presented so far is exact. To proceed further, we assume that vertices with the same degree are statistically equivalent [31]. That is, where V(k) is the set of vertices of degree k. Thus, by summing Eq. (9) for all vertices of degree k and dividing by the number of vertices with this degree, N k , we can write, after some formal manipulations, If all the vertices with the same degree are statistically equivalent, we can set [31] 1 Finally, by assuming the mean-field approximation , the rate equation for the density ρ k (t) can be written as In the case of networks with general degree correlations, the solution of Eq. (14) depends on the nature of the conditional probability P (k ′ |k) and can be a rather demanding task [35]. Therefore, in the rest of this paper we will restrict ourselves to the case of uncorrelated networks, in which the conditional probability takes the simple form P (k ′ |k) = k ′ P (k ′ )/ k [31]. For this class of networks, the rate equation Eq. (14) is simplified to the form where ρ(t) is the total density of A particles. We can obtain a differential equation for this last quantity by multiplying Eq. (15) by P (k) and summing over k, namely, where In order to solve Eq. (16), we perform a quasi-static approximation in the rate equation Eq. (15). From the mean-field solution of the A + A → ∅ process, we expect ρ(t) to be a decreasing function with a power law-like behavior. In this case, for large enough times, the time derivative of ρ(t) will be much smaller that the density proper. Extending this argument to the partial densities ρ k (t), at large times we can neglect the left-hand-side term in Eq. (15), and solve for ρ k (t) as a function of the density, obtaining Substituting this approximation into the expression for Θ(t), we get Inserting this last expression into Eq. (16), we obtain as a final equation for the density of A particles The solution of the approximate equation Eq. (20) will depend on the particular form of the degree distribution. The task becomes, however, quite simple for the class of homogeneous networks. In this case, the degree distribution decreases so quickly that all its moments are finite. So, for small ρ(t) we can perform a Taylor expansion of the right-hand-side of Eq. (20), obtaining at lowest order whose solution at large times yields being ρ 0 the initial density of A particles. This corresponds to the pure mean-field linear behavior, with a finite prefactor depending on the second moment of the degree distribution. Eq. (18) can help us to understand how this process will operate in heterogeneous networks. Indeed, for any given time t * , the partial density of vertices with degree larger than k /2ρ(t * ) is essentially constant up to time t * , that is, ρ k (t) ≃ 1/2 for k > k /2ρ(t * ) and t < t * . The reason is that vertices with high degree (hubs) are more easily reached by particles than those with small degree and, with high probability they will be surrounded by some A particle. Then, as soon as one nearby particle decides to move to the hub, both particles disappear and another nearby particle will replace the original one in the hub, keeping then the density constant, which is clearly set to 1/2. Therefore, hubs act as drains through which particles vanish. Besides, this process happens in a hierarchical way in time, that is, the partial density of small degree vertices decreases first, whereas the most connected vertices are the last to hold particles. This reasoning allows us to anticipate that the resulting dynamics will strongly depend on the number of high degree vertices of the network and, consequently, on the degree distribution.

III. SCALE-FREE NETWORKS
For heterogenous networks with a diverging second moment, as in the case of scale-free networks, we have to consider more carefully the solution of Eq. (20). Let us focus in generalized uncorrelated scale-free networks in the infinite size limit, which are completely determined by the normalized degree distribution where 2 < γ < 3 is the degree exponent, m is the minimum degree in the network, and we are approximating k as a continuous variable. The average degree is thus k = m(γ − 1)/(γ − 2). For a scale-free network in the limit N → ∞, Eq. (20) can be written as where F is the Gauss hypergeometric function [36]. Therefore, the density as function of time can be implicitly expressed in terms of the integral For very large times and small densities, we can use the asymptotic expansion of the Gauss hypergeometric function, F (1, γ − 2, γ − 1, −z) ∼ z 2−γ , z → ∞, to obtain the scaling behavior of the density with time, This same result can be derived in a more intuitive fashion starting from the argument presented at the end of the previous Section. If hubs act as drains of particles, then the rate of change in the total density of particles will be proportional to the density of hubs. More precisely, identifying the relevant hubs at time t as those vertices with degree larger than k /2ρ(t), we have that from where we obtain the same trend given by Eq. (26). As γ decreases the probability of having high degree vertices as drains in the network increases. From Eq. (26) we see as well that when this happens, the decrease in the density gets faster. Therefore, the bigger is the degree of the hubs in the network, the faster is the process of annihilation.
We can obtain an exact solution for ρ(t) in scale-free networks with degree exponent γ = 3. In this case, we have P (k) = 2m 2 k −3 and k = 2m. The differential equation for the particle density reads now whose solution is where L i is the logarithmic integral function [36]. For small ρ(t), we can exploit the asymptotic expansion of the logarithmic integral function, L i (z) ∼ z/ ln z, z → ∞, to obtain That is, in the infinite network size limit, the particle density of a scale-free network with γ = 3 follows the same decay as the mean-field solution, with a logarithmic correction.
The conclusion of this analysis is that the particle density in the A + A → ∅ process in a infinite size scale-free network with degree exponent γ follows the asymptotic form at large times with characteristic exponents The results for γ > 3 are a natural consequence of the lack of degree fluctuations ( k 2 < ∞) in this kind of scale-free networks. The analytical exponents determined above, however, are strongly affected by finite size effects. Indeed, for a power-law degree distribution, the largest weight in the sum in Eq. (20) is carried by the large k values. If the network is composed by a finite number of vertices, N , as it always happens in numerical simulations, it has a cutoff or maximum degree k c (N ), which is usually a function of the network size [37]. Thus, there exists a cross-over time t c , defined by such that, for t > t c the particle density is so small that we can approximate In this time regime, we will observe an effective linear behavior for the inverse particle density, whose effective slope depends of the network size as That is, in finite size scale-free networks, the density of A particles at very large times has the same behavior as in homogeneous networks, i.e. 1/ρ(t) is linear in t, but now with a slope that depends on the network size N . Thus, the exponents larger than one, given by Eq. (32), can only be numerically observed as transients in very large networks.

IV. NUMERICAL SIMULATIONS
As a check for the analytic predictions developed in the previous Sections, we have performed extensive numerical simulations of the diffusion-annihilation process on top of different model networks, both homogeneous and heterogeneous. The simulations are implemented using a sequential updating scheme [22]. An initial fraction ρ 0 of vertices in the networks is randomly chosen and occupied by an A particle. At time t in the simulation, a vertex is randomly chosen among the n vertices that host an A particle at that time. One of its neighbors is selected also at random. If it is empty, the particle moves and occupies it. If it contains a particle, this is annihilated with the one in the first vertex, becoming both vertices empty, and the number of occupied vertices is updated n → n − 2. In both cases, time is updated as t → t + 1/n, where n is the total number of particles at the beginning of the simulation step. In all the simulations presented in this Section, we set the initial particle density ρ 0 = 0.5.

A. Homogeneous networks
The class of homogeneous networks refers to networks with a degree distribution peaked at an average value k and decaying exponentially or faster for k ≫ k and k ≪ k . Typical examples of such networks are the Erdös and Rényi model [38] and the small-world model proposed by Watts and Strogatz (WS) [14]. We will focus in the latter in order to perform computer simulations. The WS model is defined as follows [14,39]: The starting point is a ring with N vertices, in which every vertex is symmetrically connected to its 2m nearest neighbors. Then, for every vertex, each edge connected to a clockwise neighbor is rewired to a randomly chosen vertex with probability p, and preserved with probability 1 − p. This procedure generates a graph with a degree distribution that decays faster than exponentially for large k, and average degree k = 2m. We will consider the WS model with p = 1, that is, in which all edges have been rewired. In this limit, the degree distribution of the WS network takes the form [39] Therefore, all its moments k n are finite for any value of n, and we should expect Eq. (22) to provide a good approximation for the dynamics of the A+A → ∅ process for ρ(t) → 0.
In order to check this fact, we have carried out large scale numerical simulations of the A + A → ∅ reaction on WS networks with p = 1 and minimum degree m = 2, which corresponds to an average degree k = 6. Simulations were performed on graphs of size N = 10 6 , averaging over 1000 reaction processes over 10 different realizations of the random network. In Fig. 1 we show the total density of A particles as a function of time. As we can observe, the inverse particle density follows quite precisely a linear behavior with time. A linear fit for the last two decades of t values yields a slope ≈ 1.50, which is not too different from the predicted value in Eq. (22), 2 k 2 / k 2 = 2.25.

B. Heterogeneous networks: The Barabási-Albert model
The Barabási-Albert (BA) graph was introduced as the first growing network model capable of producing as an emerging property a power-law degree distribution [15]. This model is based on the preferential attachment paradigm, a rather intuitive mechanism in which new individuals tend to develop more easily connections with already well connected individuals. The model is defined as follows: We start from a small number m 0 of vertices, and at each time step, a new vertex is introduced, with m edges that are connected to old vertices i with probability where k i is the degree of the ith vertex. After iterating this procedure a large number of times, we obtain a network composed by N vertices, minimum degree m, and average degree k = 2m, with a degree distribution  P (k) = 2m 3 k −3 and almost vanishing degree correlations [34]. In the results presented here we consider the parameters m 0 = 5 and m = 2, corresponding to an average degree k = 4. Simulations were performed on networks of size N = 10 5 and N = 10 6 , averaging over 1000 runs in 10 different network samples. A first issue that we have explored in this network is the validity of the quasi-stationary approximation made to obtain Eq. (18). If this approximation is valid, we should observe that the quantity is independent of k. We have performed this check in Fig. 2. As we can observe, the curves collapse quite nicely for t larger that t 0 ≃ 10. Therefore, we should expect that the results obtained from the approximation in Eq. (20) will be valid above this characteristic time. In Fig. 3 we represent the inverse particle density as function of time. In this plot we observe that, except for very early times, the initial time dependence is very accurately described by the analytic solution Eq. (29), which can be approximately described by a linear behavior with a logarithmic correction, 1/ρ(t) ∼ t ln t (full line). At larger times, the function 1/ρ(t) crosses over to a linear dependence (dashed line), showing the onset of finite size effects. The poor agreement at early times should be attributed to the approximation in Eq. (20), which is valid only for times larger than t 0 ≃ 10.

C. Heterogeneous networks: The uncorrelated configuration model
The construction of uncorrelated scale-free networks with arbitrary degree exponent is a nontrivial issue. Traditionally, this kind of networks were generated using the configuration model (CM) [26,27,28,29], in which, starting from a degree sequence extracted according to the desired degree distribution, edges are randomly created between vertices, respecting the preassigned degrees. It has been shown, however, that for scale-free distributions with diverging second moment, and in the absence of self-connections (a vertex joined to itself) and multiple connections (two vertices connected by more than one edge), this procedure generates in fact degree correlations [40,41]. The origin of this phenomenon can be traced back to the effects of the cut-off, or maximum expected degree, k c (N ), which must scale at most as N 1/2 in order to allow for a full random uncorrelated edge assignment [42].
To solve this question, it has been recently proposed the uncorrelated configuration model (UCM) [43], that is defined as follows: 1. Assign to each vertex i in a set of N initially disconnected vertices a degree k i , extracted from the probability distribution P (k) ∼ k −γ , and subject to the constraints m ≤ k i ≤ N 1/2 and i k i even.
2. Construct the network by randomly connecting the vertices with i k i /2 edges, respecting the preassigned degrees and avoiding multiple and selfconnections.
Using this algorithm, it is possible to create scale-free networks whose cut-off scales as k c (N ) ∼ N 1/2 for any degree exponent γ, and which are completely uncorrelated. The results presented below were obtained from simulations performed on networks generated from the UCM algorithm with minimum degree m = 2 and sizes ranging from N = 10 5 to N = 10 6 . Averages were performed over 1000 runs in 10 different network.
In Fig. 4 we plot the inverse particle density from computer simulations in networks with different degree exponent γ. From this plot we observe that, at the initial time regime, the growth of this function is faster for smaller values of γ, in agreement with the theoretical prediction Eq. (26). At larger times, on the other hand, finite size affects take over, and we observe again a linear regime (compare the slope with the reference dashed line in the Figure). For γ < 3, however, we cannot compare directly with the full analytic prediction Eq. (26), as we did for the BA model. The reason is the following: From Eq. (34), we have that the density at the crossover time scales as Then, from Eq. (26), we can estimate That is, the crossover time t c diminishes with decreasing γ. It turns out that, even for γ = 2.5, t c is so small that we cannot observe the transient regime 1/ρ(t) ∼ t 1/(γ−2) , and only the finite size effects, with its associated linear behavior, are apparent. System sizes larger than those available for our present computer resources should make visible this crossover phenomena. Finally, we have checked the slope dependence of the linear behavior given by the finite size effects. According to Eqs. (36) and (37), we should observe, for a fixed degree exponent, a slope that grows with the system size as k c (N ) 3−γ ∼ N (3−γ)/2 , where we have used the scaling of the cut-off given by the UCM algorithm. We have analyzed this point in Fig. 5. For a degree exponent γ = 2.5, we observe that the curves for increasing values of N show a correspondingly increasing of the slope in the final linear region. We have estimated this slope by performing a linear fitting to the whole curve (inset in Fig. 5). This slope increases approximately as N 1/4 , in very good agreement with our theoretical predictions.
It is noteworthy that these finite size effects, so remarkable in the simulations presented in this Section, are however not evident in the numerical work developed in Ref. [23]. The reason of this fact could be that the simulations in [23] were performed in networks generated with the CM algorithm. These networks, as we have discussed above, exhibit some degree of correlations [40,41,43], and therefore could be outside of the regime of validity of the analytical calculations we have presented here.

V. CONCLUSIONS
In this paper we have presented a detailed analytical study of the A + A → ∅ diffusion-annihilation process in complex networks [23]. For uncorrelated homogeneous networks with bounded degree fluctuations, we recover a behavior compatible with the standard mean-field solution of the process, that is, the inverse particle density grows linearly with time, 1/ρ(t) ∼ t, with a constant prefactor. In the case of uncorrelated heterogeneous networks, characterized by a scale-free degree distribution, we observe instead that, in the infinite network size limit N → ∞, the inverse particle density increases as a power law with time, 1/ρ(t) ∼ t α(γ) , with an exponent larger than one and that depends on the level of heterogeneity of the network, i.e. α(γ) = 1/(γ − 2) for γ < 3. For γ = 3, a linear growth with logarithmic corrections sets in, while for γ > 3 we recover again the mean-field solution typical of homogeneous networks. In the case of scale-free networks, we have also analyzed the effects of a finite network size in this dynamical process. For any value of γ ≤ 3, we have shown that, at large times, the inverse particle density in a finite network crosses over to a linear behavior, with a slope that is an increasing function of the network size N , i.e. 1/ρ(t) ∼ N (3−γ)/2 t. In order to check our results, we have performed extensive numerical simulations on both homogeneous and heterogeneous networks, obtaining a convincing evidence for our predictions.
An interesting conclusion that can be extracted of the present work is the very strong incidence that finite size effects can have on dynamical systems on top of scalefree networks. While these size effects have been already discussed in the context of epidemic spreading [44], the A+A → ∅ process analyzed here provides a further example, in which analytical results and numerical simulations show a striking agreement.