text:
# A Holistic Scenario of Turbulent Molecular Cloud Evolution and Control of the Star Formation Efficiency. First Tests

## 1 Introduction

Two of the main questions about molecular cloud (MC) structure and star formation concern the origin of molecular cloud turbulence and of the low efficiency of star formation (SFE). Indeed, the turbulent properties of molecular clouds are remarkably uniform, even in clouds without stellar energy sources (see, e.g., McKee 1999 and references therein). And, defined as the fraction of the mass in a MC that ends up in stars before the cloud is dispersed, the SFE is well known to be as low as a few percent with respect to the total cloud mass (see, e.g., Evans 1991 and references therein).

It is known that turbulence can provide additional support against global gravitational collapse (Chandrasekhar 1951; Bonazzola et al. 1987, 1992; Vázquez-Semadeni & Gazol 1995). In particular, numerical studies have suggested that turbulence is able to prevent the collapse when the energy injection scale is smaller than the Jeans length at typical core densities (Léorat, Passot & Pouquet 1990; Klessen, Heitsch & Mac Low 2000, hereafter KHM00). Simultaneously, supersonic, compressible turbulence is known to produce “fragmentation”, i.e., the formation of smaller-scale density substructures that can themselves possibly undergo local collapse (Sasao 1973; Hunter & Fleck 1982; Tohline, Bodenheimer & Christodolou 1987; Scalo 1987; Elmegreen 1993; Padoan 1995; Vázquez-Semadeni, Passot & Pouquet 1996; KHM00; Padoan et al. 2001; Heitsch, Mac Low & Klessen 2001). Padoan (1995) and KHM00 have proposed that turbulent fragmentation has an important role in determining the SFE. The latter authors have shown that the efficiency in numerical simulations (measured as the fraction of mass in collapsed objects) increases as the energy injection scale is increased. However, this result was obtained in simulations with an idealized random forcing applied at specific scales, and ubiquitously in time and space. Instead, the turbulence in real MCs is probably injected over a wide range of scales, but at localized sites and times (Scalo 1987; Norman & Ferrara 1996; Avila-Reese & Vázquez-Semadeni 2001). In this case, the driving scale is not uniquely defined. Also, if molecular clouds form through the collisions of larger-scale streams in the more diffuse ISM (Elmegreen 1993; Ballesteros-Paredes, Vázquez-Semadeni & Scalo 1999; Vázquez-Semadeni, Ballesteros-Paredes & Klessen 2002, hereafter Paper I), then they can be part of a turbulent cascade originating from much larger scales, and thus the driving scale can be much larger than the Jeans length in the molecular clouds (Mac Low & Ossenkopf (2001). This appears inconsistent with the requirement of a driving scale smaller than the cloud’s Jeans length in order to maintain a low efficiency.

In Paper I, we have compiled a holistic scenario in which both the turbulence in molecular clouds and the SFE are part of the life cycle of MCs: MC turbulence originates from the very process that forms the clouds, namely colliding streams in the diffuse ISM (Elmegreen 1993; Ballesteros-Paredes et al. 1999), through bending-mode instabilities (Hunter et al. 1986; Vishniac 1994; Walder & Folini 1998; Klein & Woods 1998). In this sense, MC turbulence is part of a compressible and lossy cascade which may originate at intermediate scales in the Galactic disk ($`1`$ kpc) from various instabilities (Parker 1966; Toomre 1964; Elmegreen 1991), supernova energy input (e.g., McCray & Snow 1979; MacLow 2002), and stellar winds (Norman & Silk 1980), and continues towards smaller scales, as the internal turbulence produces yet smaller-scale turbulent compressed layers (Kornreich & Scalo 2000). The compressible cascade has to be lossy because shocks are known to be dissipative by transferring a fraction of the energy involved directly to the dissipative scales (Kadomtsev & Petviashvili 1973). This compressible-cascade regime, however, must end at a “sonic” scale, denoted $`\lambda _\mathrm{s}`$, below which the turbulence internal to the new structures becomes subsonic and therefore cannot produce any further subfragmentation (Padoan 1995). This does not necessarily imply that the dissipative (or “inner”) scale of the turbulence has been reached, but only that at this point the cascade can turn into an incompressible one. Thus, in this scenario, subsonic “cores” are the natural end of the compressible part of the cascade in the ISM.

When the turbulence becomes subsonic, it ceases to be a stronger source of global support than thermal pressure, and the Jeans criterion can be applied to determine whether a “core” is gravitationally unstable and should proceed to collapse (Padoan 1995), or else rebound and merge back into its environment (Elmegreen 1993; Taylor, Morata & Williams 1996; Ballesteros-Paredes et al. 1999; Vázquez-Semadeni, Shadmehri & Ballesteros-Paredes 2002). This effectively sets an upper limit to the SFE, because, if globally the cloud is supported by its internal turbulence, only massive enough clumps among those with sizes $`l\lambda _\mathrm{s}`$ are susceptible to local collapse (Padoan 1995). In this sense, star-forming cores constitute the intersection of the set of subsonic regions with the high-mass, super-Jeans tail of the mass distribution (Padoan 1995; Ballesteros-Paredes & Vázquez-Semadeni 1997; Shadmehri, Vázquez-Semadeni & Ballesteros-Paredes 2001). In particular, Padoan (1995) gave a semi-phenomenological theory for the efficiency and mass distribution of proto-stellar cores in turbulent molecular clouds (recently extended by Padoan & Nordlund 2002). In that work, however, he considered only a power-law turbulent spectrum, neglecting the necessary existence of inner and outer scales, and the role of the energy injection scale, which was proposed to be crucial by Léorat et al. (1990) and KHM00.

In this Letter, besides putting together and advancing the above global scenario, we present (§2) a first set of numerical tests supporting it, by showing that the scale $`\lambda _\mathrm{s}`$ indeed plays a crucial role in limiting the SFE. Clearly, this is only a partial ingredient, as the mass distribution in the flow is also expected to be crucial (Padoan 1995). Our experiments also support the notion of the dual role of turbulence in providing support to large scales while promoting small-scale collapse, as we discuss in §3, where we also give a re-interpretation of the results of KHM00 in the framework of our scenario.

## 2 Numerical experiments and results

We consider a set of 12 numerical experiments similar to those reported by KHM00. We refer the reader to that paper for details (see also Klessen & Burkert 2000, 2001). Here we just mention that they are smoothed particle hydrodynamics (SPH) simulations of non-magnetic turbulent flows using 205,379 particles, with self-gravity, and random forcing applied only at wavenumbers in a narrow interval $`k_\mathrm{d}1|𝐤|k_\mathrm{d}`$, where $`k_\mathrm{d}L/\lambda _\mathrm{d}`$ is the driving wavenumber, $`L`$ is the computational box size and $`\lambda _\mathrm{d}`$ is the driving scale. The simulations are carried out in non-dimensional units, in which there are 64 thermal Jeans masses (with respect to the mean density) in the computational box, the time unit is $`t_0=t_{\mathrm{ff}}/1.5`$, where $`t_{\mathrm{ff}}`$ is the free fall time, the sound speed is $`c=0.1`$, and the mean density is $`\rho _0=1/8`$. A constant kinetic energy input rate is used, with an amplitude chosen to maintain a fixed level of turbulent kinetic energy $`E_\mathrm{t}`$. For our analysis, the SPH data are interpolated onto a cubic grid of $`128^3`$ cells, but note that the effective resolution of the SPH method in the density peaks (cores) is higher than that.

The runs are labeled mnemonically as M$`x`$K$`y`$, where $`x`$ denotes the rms Mach number and $`y`$ denotes the forcing wavenumber $`k_\mathrm{d}`$. The rms Mach number $`M_\mathrm{s}`$ and the energy injection scale $`\lambda _\mathrm{d}`$ were taken by KHM00 as the controlling parameters in their investigation of the amount of mass in collapsed objects $`M_{}`$ as a function of time for several simulations. The numerical scheme allows this mass to be measured through the usage of “sink” particles, which are particles that replace a collapsing gas clump once it becomes too dense and small to be properly resolved. A sink particle inherits the combined masses of the replaced SPH particles, as well as their linear and angular momenta, and has the ability to continue accreting further SPH particles.

In order to investigate the role of the sonic scale $`\lambda _\mathrm{s}`$, we consider the velocity dispersion-vs.-size relation, together with the mass-accretion histories of the simulations. The velocity dispersion-size relation we consider is not restricted to the densest, or highest-column density projected structures, as is customary in both oservational (Larson 1981; see Blitz 1993 for a more recent review) and numerical (Vázquez-Semadeni, Ballesteros-Paredes & Rodríguez 1997; Ostriker, Stone & Gammie 2001; Ballesteros-Paredes & Mac Low 2002) astrophysical studies. Instead, we consider a subdivision of the whole computational box in regions of size $`l`$, and measure the velocity dispersion in all such regions. The mass accretion histories are computed as the fraction of the total mass in collapsed objects (sink particles) as a function of time.

Figure 1 shows the velocity dispersion-size relation for all runs. For every run, the solid line shows the average velocity dispersion for all regions of size $`l`$; the “error” bars indicate the range of values of the velocity dispersion exhibited by all regions of size $`l`$. The dotted straight lines show an $`R^{1/2}`$ power law, and the dashed horizontal line shows the value of the sound speed. The sonic scale can be defined as that at which the solid line crosses the sound speed level, and is depicted by a vertical dotted line. Because of the large error (or, more properly, “scatter”) bars in the velocity dispersion, the sonic scale is typically defined only to within factors $`4`$. The caveat should be mentioned that the velocity dispersion-size plots we obtain are clearly not power laws, with progressively smaller scales being progressively more deficient in turbulent kinetic energy. This is in part a consequence of the spatially varying resolution of SPH, with low resolution at low density (in voids) and high resolution at high density (in cores) (Ossenkopf, Klessen, & Heitsch 2001), and in part a consequence of a power-law range not being expected at scales larger than $`\lambda _\mathrm{d}`$ for the small-driving-scale runs. In any case, the non-power-law nature of the spectrum at small scales implies that our measurements overestimate $`\lambda _\mathrm{s}`$. We discuss the implications of this in §3.

Figure 2 shows the accretion histories for all runs.It is noteworthy that the accreted mass in most runs (except the ones with the lowest accretion rates) eventually saturates at some level. This phenomenon is extremely interesting in its own right; it clearly appears to be a real phenomenon, generally seen in simulations that allow following the evolution of the system beyond the formation of individual collapse sites through the use of sink particles (KHM00; Klessen, Burkert & Bate 2001; Klessen 2001; Bate, Bonnell & Bromm 2002). We speculate that it occurs because the local collapse events decrease the mean density of the “interclump” medium, making it more difficult to form new collapsing sites from it. Concerns that this saturation level be a numerical artifact dependent of the resolution are dispelled by fig. 6 of KHM00, which shows no clear trend of the saturation level with resolution, but only statistical scatter (different random initial conditions with the same global parameters show variations of up to a factor of 2 in the final accreted masses).

Operationally, defining the SFE is not an unambiguous task, because not all runs have had time to saturate yet (assuming they all eventually do). Given that recent studies suggest that clouds have lifetimes of only a few crossing times (Ballesteros-Paredes et al. 1999; Elmegreen 2000), we choose to define the SFE as the collapsed mass fraction after two turbulent crossing times ($`t_\mathrm{c}`$).111We thank the referee for suggseting this criterion. In the units of the code, these correspond to $`t=10`$, 12.5, 6.7 and 4 for runs with rms Mach number $`M_\mathrm{s}=2`$, 3.2, 6, and 10, respectively. However, for the runs with $`M_\mathrm{s}=2`$ and 3.2, this time is longer than they have been evolved. Fortunately, the collapsed mass fractions in all these runs have saturated, and so for them we take the saturation value. Due to the statistical variation of the collapsed mass fractions mentioned above, we also put error bars of $`\pm 50`$% in the reported values of the SFE.

In order to interpret the results of the simulations, fig. 3 shows the SFE versus $`\lambda _\mathrm{d}`$ (a), $`M_\mathrm{s}`$ (b), and $`\lambda _\mathrm{s}`$ (c) for all simulations. The linear Pearson correlation coefficients for the logarithmic data in figs. 3a, 3b and 3c are respectively 0.35, $`0.73`$ and 0.75. showing that indeed the SFE correlates best with the sonic scale $`\lambda _\mathrm{s}`$. Note that this is in spite of the fact that the relationship between the SFE and $`\lambda _\mathrm{s}`$ is clearly not a power law, but instead saturates at large $`\lambda _\mathrm{s}`$. A change of variable $`\lambda _\mathrm{s}u=0.00380.012\mathrm{exp}2.39\mathrm{log}\lambda _\mathrm{s}`$ gives the fit shown by the solid curve in fig. 3c, with a much tighter correlation coefficient of 0.86. Nevertheless, the correlation with $`M_\mathrm{s}`$ is also significant. We discuss this in the next section.

## 3 Discussion and comparison with previous work

Although in the previous section we have fitted a smooth curve to the data in fig. 3c, actually it appears as if there are two functional dependences (“modes”) of the SFE on $`\lambda _\mathrm{s}`$ in this plot: for $`\lambda _\mathrm{s}0.03`$, the SFE is nearly a power-law on $`\lambda _\mathrm{s}`$, while it is roughly constant (within a factor of 2) for $`\lambda _\mathrm{s}0.03`$. All the points in the former regime correspond to cases with $`M_\mathrm{s}=6`$ or 10, which are globally supported by the turbulence, while most of the points in the latter regime correspond to $`M_\mathrm{s}=2`$ and 3.2, which are globally unsupported. Nevertheless, the large-scale driven ($`k_\mathrm{d}=2`$) high-$`M_\mathrm{s}`$ cases appear to be at the crossover between the two regimes. This demonstrates the suggested dual nature of the turbulence (Sasao 1973; Falgarone, Puget & Pérault 1992; Vázquez-Semadeni & Passot 1999; KHM00): globally supported runs can be induced to behave as unsupported ones when the bulk of the turbulent energy is at the largest scales.

It is tempting to speculate that these two “modes” correspond to the isolated and clustered observed modes of SF. In this case, our results suggest that the two would not be physically separated modes, but simply represent a continuous transition from one functional dependence of the SFE on $`\lambda _\mathrm{s}`$ to the other, as the domain transits from behaving as if globally supported to behaving as if globally unsupported. We see also that, even in the globally unsupported case, the efficiency is below unity because of the saturation described above. In this case, turbulence would only provide the seeds for the internal fragmentation of the collapsing molecular cloud region (see, e.g., Klessen et al. 1998; Klessen 2001), the efficiency being more strongly dependent on stellar feedback to destroy the cloud. Further testing is necessary in order to confirm the correspondence between the two regimes found here and the observationally identified isolated and clustered modes.

Concerning the resolution of the numerical scheme, as mentioned above, $`\lambda _\mathrm{s}`$ is overestimated at small scales by our simulations. Moreover, in the absence of numerical diffusion, one would expect the density peaks to become narrower as well, possibly increasing the number of them that go into collapse, and thus increasing the SFE. Thus, in general we expect the correlation seen in fig. 3c to be shallower if the diffusivity is decreased. Nevertheless, we do not expect significant qualitative changes in our results.

Our reinterpretation in terms of $`\lambda _\mathrm{s}`$ eliminates the apparent inconsistency implied by previous works that a low SFE would seem to require a small driving scale, yet the turbulence in molecular clouds appears to be generally driven from larger scales. Indeed, our results show that a small driving scale alone is not enough to prevent collapse if it is too weak, while driving at large scales can still decrease the efficiency if it is strong enough, as pointed out by KHM00. That is, the controlling parameter is not the driving scale nor the rms Mach number alone, but instead their combined action expressed as the energy content at the sonic scale. Previous studies did not consider large-scale-driven cases with a strong enough turbulence level to decrease the SFE.

A few final remarks: 1) The mechanism sketched here is expected to operate in addition to any other mechanism limiting the SFE based on energy feedback from the massive stars, which can disrupt the parent cloud. 2) The occurrence of a saturation of the accreted mass implies that clouds could not achieve 100% SFE even if they had arbitrarily large amounts of time available, provided the turbulence level is kept constant. This is because local collapse uses up part of the mass, rendering the remaining gas more diffuse and easy to disperse. Thus, partial dispersion may also be part of the natural course of evolution of molecular clouds. 3) In this work we have neglected the effect of the magnetic fields. Although their effect certainly must be incorporated in a full theory, the mechanism discussed here is capable of regulating the SFE without relying on their presence, depending only on the spatial redistribution of the physical quantities in a compressible turbulent flow. This is particularly relevant upon recent suggestions that molecular clouds may be super-Alfvénic (Padoan & Nordlund 1999) and that cores are in general supercritical (Nakano 1998; Hartmann, Ballesteros-Paredes & Bergin 2001; Bourke et al. 2001; Crutcher, Heiles & Troland 2002).

We thank L. Hartmann and the referee, P. Padoan, for useful suggestions. We acknowledge financial support from Conacyt grants 27752-E to E.V.-S. and I39318-E to J.B.-P., and from the Emmy Noether Program of the Deutsche Forschungsgemeinschaft (KL1358/1) and the Center for Star Formation Studies (NASA-Ames, UC Berkeley, and UC Santa Cruz), to R.S.K.


format: text

text:
# Comparison of the Total Charged-Particle Multiplicity in High-Energy Heavy Ion Collisions with 𝑒⁺𝑒⁻ and 𝑝𝑝/𝑝̄𝑝 Data

## Abstract

The PHOBOS experiment at RHIC has measured the total multiplicity of primary charged particles as a function of collision centrality in Au+Au collisions at $`\sqrt{s_{NN}}=`$ 19.6, 130 and 200 GeV. Above $`\sqrt{s}20`$ GeV, the total multiplicity per participating nucleon pair ($`N_{ch}/N_{part}/2`$) in central events scales with $`\sqrt{s}`$ in the same way as $`N_{ch}`$ in $`e^+e^{}`$ data. This is suggestive of a universal mechanism of particle production in strongly-interacting systems, controlled mainly by the amount of energy available for particle production (per participant pair for heavy ion collisions). The same effect has been observed in $`pp/\overline{p}p`$ data after correcting for the energy taken away by leading particles. An approximate independence of $`N_{ch}/N_{part}/2`$ on the number of participating nucleons is also observed, reminiscent of “wounded nucleon” scaling ($`N_{ch}N_{part}`$), but with the constant of proportionality set by the multiplicity measured in $`e^+e^{}`$ data rather than by $`pp/\overline{p}p`$ data.

Central collisions of two gold nuclei at the top energy of the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory produce thousands of charged particles. These are the largest particle multiplicities generated in man-made subatomic reactions. The hope is that these complex systems may reveal evidence of the creation and decay of a Quark-Gluon Plasma (QGP), where quarks and gluons are allowed to explore a volume larger than that of a typical hadron.

The high multiplicities in heavy ion collisions typically arise from the large number of nucleon-nucleon collisions which occur, with many of the nucleons struck several times as they pass through the oncoming nucleus. Studies of proton-nucleus collisions demonstrated that the total multiplicity ($`N_{ch}`$) does not scale proportionally to the number of binary collisions ($`N_{coll}`$) in the reaction, but rather was found to scale more closely with the number of “wounded nucleons” which participate inelastically ($`N_{part}`$) elias ; wounded . For example, the number of participants is $`N_{part}=2`$ for a proton-proton collision and $`N_{part}=(N_{coll}+1)`$ for a proton-nucleus collision. Thus, by scaling the particle yields measured in heavy ion collisions by $`N_{part}/2`$, data from heavy ion collisions may be directly compared with similar yields in elementary $`pp`$, $`\overline{p}p`$ or even the annihilation of $`e^+e^{}`$ into hadrons.

While both $`e^+e^{}`$ and $`pp/\overline{p}p`$ collisions must ultimately allow a description based on Quantum Chromodynamics (QCD), the theory of the strong interaction, the evolution of these two systems tends to be understood in different ways. The large momentum transfer to the outgoing produced quark and anti-quark in $`e^+e^{}`$ reactions allows the use of perturbative QCD (pQCD) to describe the spectrum of quarks and gluons radiated as the system fragments basics . Minimum bias collisions of hadrons are not generally thought to be amenable to such a perturbative description, since the transverse momentum exchanges involved are typically less than 1 GeV/c. Instead, phenomenological approaches (e.g. PYTHIA pythia ) are used to describe most of the (predominantly soft) particles produced in high energy $`pp`$ or $`\overline{p}p`$ collisions.

A basic connection between perturbative and non-perturbative physics has been elucidated by simultaneous measurements of the multiplicity of charged particles and the high-momentum “leading” protons in $`pp`$ collisions at the ISR. Basile et al. basile  found that the average multiplicity $`N_{ch}`$ in $`pp`$ collisions is similar to that for $`e^+e^{}`$ collisions with $`\sqrt{s}_{e^+e^{}}=\sqrt{s_{\mathrm{eff}}}`$, where $`\sqrt{s_{\mathrm{eff}}}`$ is the $`pp`$ center-of-mass energy minus the energy of the leading particles. This is interpreted as a universal mechanism of particle production controlled dominantly by the available center of mass energy basile .

In this Letter, we report results from the PHOBOS experiment on the total multiplicity of primary charged particles $`N_{ch}`$ as a function of $`N_{part}`$ for heavy ion collisions at $`\sqrt{s_{NN}}=`$ 19.6, 130 and 200 GeV, where $`\sqrt{s_{NN}}`$ is the nucleon-nucleon center-of-mass energy. Comparisons with $`pp/\overline{p}p`$ and $`e^+e^{}`$ data are made to investigate whether this universal mechanism of particle production applies in the context of heavy ion collisions.

The PHOBOS multiplicity detector consists of two arrays of silicon detectors which cover nearly the full solid angle for collision events. The “Octagon” detector surrounds the interaction region with a roughly cylindrical geometry covering $`|\eta |<3.2`$. Two sets of three “Ring” detectors are placed far forward and backward of the interaction point and surround the beam pipe, covering $`3<|\eta |<5.4`$. The methods used for measuring the multiplicity of charged particles as well as for extracting $`N_{part}`$ have been described in more detail in Refs. phobos\_cent\_200 ; phobos\_limfrag .

Using the data presented in Ref. phobos\_limfrag , Fig. 1a shows $`dN_{ch}/d\eta /N_{part}/2`$ averaged over the forward and backward hemispheres for the most central Au+Au events at $`\sqrt{s_{NN}}=`$ 200 GeV. The systematic errors (representing a 90% CL interval) depend on $`\eta `$ and are shown on the figure as a shaded band. The Au+Au data are compared with $`dN_{ch}/d\eta `$ for non-single diffractive (NSD) $`\overline{p}p`$ collisions ua5  and $`dN/dy_T`$ for $`e^+e^{}`$ collisions (with cuts applied to reject large initial-state photon radiation) ALEPH  at $`\sqrt{s}=`$ 200 GeV. The variable $`y_T`$ is the rapidity of charged particles relative to the event thrust axis, assuming all particles to have the pion mass.

It is observed that the Au+Au data are very similar in magnitude and shape to the $`e^+e^{}`$ data at the same $`\sqrt{s}`$, and similar in shape to the $`\overline{p}p`$ data (as shown in Fig. 1b), over a large range in $`\eta `$. The differences between the $`e^+e^{}`$ and Au+Au distributions shown in Fig. 1a can be partly attributed to the different kinematic variables. JETSET calculations indicate that the $`y_T`$ distribution is slightly narrower than the corresponding pseudorapidity distribution in $`e^+e^{}`$ collisions, with a higher plateau height. Yet, even without taking this into account, the difference between the distributions is no more than $`\pm 10\%`$ for $`|\eta |`$ and $`|y_T|<4`$ pythia . However, the same calculations also show that the choice in kinematic variables does not explain the difference in the forward region (above $`|\eta |=4`$), although this may not be surprising, as this region should show some residual effect of the presence of the spectator nucleons.

The similarity of the angular distributions indicates that the total yield of charged particles in $`e^+e^{}`$ and central Au+Au collisions should also be similar for the same $`\sqrt{s}`$, when the nuclear data are scaled by the number of participant pairs. To correct for the small acceptance losses in the PHOBOS apparatus (which covers $`|\eta |<5.4`$), we have used several methods inspired by the excellent agreement of the lowest energy PHOBOS data with the higher energy data when shown as a function of $`\eta ^{}=\eta y_{beam}`$ phobos\_limfrag . PHOBOS data from $`\sqrt{s_{NN}}=19.6`$ GeV for $`\eta >2.5`$, shifted by $`\mathrm{\Delta }\eta =y_{200}y_{19.6}=2.32`$ (the difference in beam rapidities between the two energies), displays the limiting behavior discussed in Ref. phobos\_limfrag . This effectively extends the rapidity coverage to $`\eta 8`$. A Woods-Saxon function for $`dN/dy`$ fit to the Au+Au data, also provides a reasonable description of the $`dN/d\eta `$ distribution, and extrapolates through the lower energy data as well. Thus, in one method, we integrate $`dN_{ch}/d\eta `$ for $`\sqrt{s_{NN}}=`$ 130 and 200 GeV for $`\eta ^{}<0`$ and then use the PHOBOS data at $`\sqrt{s_{NN}}`$ = 19.6 GeV for $`\eta ^{}>0`$. We also integrate Woods-Saxon fits, similar to that shown in Fig. 1a, for $`|\eta |<8`$. These two approaches agree within 2% for central events. For the lowest RHIC energy, we simply integrate the charged particles in the PHOBOS acceptance.

In Fig. 2a, data on $`N_{ch}`$ from $`pp`$, $`\overline{p}p`$, $`e^+e^{}`$ and central heavy ion collisions (scaled by $`N_{part}/2`$) are shown as a function of $`\sqrt{s}`$. The $`pp`$, $`\overline{p}p`$, and $`e^+e^{}`$ data and errors are taken from a compilation Groom:in  and no further corrections are applied. The errors shown are the quadratically combined statistical and systematic errors. Heavy ion data are shown for central Au+Au events at RHIC (this work), Au+Au events from E895 at the AGS ($`\sqrt{s_{NN}}=2.64.3`$ GeV) e895  and Pb+Pb events from NA49 at the SPS ($`\sqrt{s_{NN}}=`$ 8.6, 12.2 and 17.3 GeV) na49 . A PHOBOS Au+Au data point at $`\sqrt{s_{NN}}=56`$ GeV has been added by using the measured value at midrapidity phobos\_dndeta  and using the universal limiting distribution described in Ref. phobos\_limfrag  to approximate the shape of the full distribution. All of the errors shown for the heavy ion data are systematic.

Perturbative QCD calculations are able to predict the dependence of the total multiplicity in $`e^+e^{}`$ collisions as a function of $`\sqrt{s}`$, $`N_{e^+e^{}}(s)=C\alpha _s(s)^A\mathrm{exp}(\sqrt{B/\alpha _s(s)})`$, with $`A`$ and $`B`$ fully calculable within pQCDMueller:cq . The QCD scale $`\mathrm{\Lambda }_{QCD}`$ is set to 225 MeV, leaving only a constant of proportionality $`C`$ free to fit to the experimental data. A fit to the $`e^+e^{}`$ data has been made with this expression (“$`e^+e^{}`$ fit”) and has been used in Fig. 2b to see how the various systems compare with $`e^+e^{}`$ data by scaling all of the data at a given $`\sqrt{s}`$ by this function.

Fig. 2b shows that the $`pp/\overline{p}p`$ data are about 30% below $`e^+e^{}`$ over the full range of energies. However, rescaling the $`\sqrt{s}`$ of each point by a factor of $`1/2`$, $`\sqrt{s_{\mathrm{eff}}}=\sqrt{s}/2`$, brings the data into reasonable agreement with the $`e^+e^{}`$ trend, as shown by the open diamonds. This is consistent with measurements of leading protons in $`pp`$ collisions, which find $`dN/dx_F`$ (where $`x_F`$ = $`2p_z/\sqrt{s}`$ in the collider reference frame) to be approximately constant for non-diffractive events over a large range of $`\sqrt{s}`$ batista  and thus $`x_F1/2`$.

Unlike the $`pp/\overline{p}p`$ data, the heavy ion data does not follow the $`e^+e^{}`$ trend over the whole energy range. Instead, they lie below the $`pp`$ data at AGS energies, crosses through the $`pp`$ data between AGS and SPS energies, and joins smoothly with the $`e^+e^{}`$ data above the top SPS energy. Thus, at high energies, the multiplicity measured per participant pair in Au+Au collisions scales in a similar way to $`e^+e^{}`$ data at the same $`\sqrt{s}`$. If we understand the lower effective $`\sqrt{s}`$ in $`pp`$ collisions as stemming from the “leading particle effect”, where the leading proton carries off a substantial amount of the available energy, the Au+Au data suggest a substantially reduced leading particle effect in central collisions of heavy nuclei at high energy.

The alleviation of the leading particle effect might not be so surprising in central nuclear collisions. In the Glauber model, each participating nucleon is typically struck $`46`$ times on average as it passes through the oncoming gold nucleus in a central event (the exact value of $`\overline{\nu }`$ depending on the energy-dependent nucleon-nucleon inelastic cross section, $`\sigma _{NN}(s)`$). One could speculate that the multiple collisions simultaneously excite and dissociate the participating nucleons, transferring much more of the energy from the forward direction towards midrapidity than found in an average $`pp/\overline{p}p`$ collision - but ultimately limited by the total incident energy. This hypothesis should be testable in proton-nucleus collisions, by measuring particle yields as a function of $`\nu `$ as was done in Refs. e910 ; na49-pp . The data in those references suggest that pion yields, whether in the projectile region ($`y>0`$) na49-pp  or integrated over $`4\pi `$ e910  increase rapidly for $`\nu <3`$ and then much more slowly for $`\nu >3`$. However, limited experimental acceptances and theoretical uncertainties preclude making any strong conclusions regarding the relationship between the energy loss of the projectile and the total charged multiplicity.

In Fig. 3 $`N_{ch}/N_{part}/2`$ is shown for PHOBOS data at three RHIC energies as a function of $`N_{part}`$. The 90% CL systematic error on the centrality dependence of $`N_{ch}/N_{part}/2`$ is shown as a shaded band, and represents a combination of several factors, dominated by the uncertainty of the extrapolation procedure to extract $`N_{ch}`$ over the full solid angle.

It might have been expected that, in events with larger impact parameters, each participant would have fewer collisions on average and thus not be fully dissociated. However, within the systematic errors, the total yield per participant pair is approximately constant (within 10%) over the measured centrality range, $`65<N_{part}<358`$, which corresponds to $`3<\overline{\nu }<6`$, where $`\overline{\nu }`$ is the average number of collisions undergone by each oncoming nucleon. Thus, it appears that only the first few collisions have any appreciable effect on particle production. It should be noted that this simple scaling is not observed for particle yields measured in a limited pseudorapidity range near midrapidity phobos\_cent\_200 .

Proton-antiproton data exist at 200 GeV, but not for the other two RHIC energies. We use a parametrization of $`pp`$ data from Ref. heiselberg , $`N_{ch}=4.2+4.69s^{.155}`$, for 19.6 and 130 GeV. Several measurements exist in $`e^+e^{}`$ at 200 GeV, but not for the other two energies. For these we use the pQCD formula for $`N_{e^+e^{}}`$, the quality of the fit clearly indicated in Fig. 2b. Fig. 3 shows that the Au+Au data are consistent with “wounded nucleon” scaling, in that the multiplicity is proportional to $`N_{part}`$ ($`N_{ch}N_{part}`$). However, $`N_{ch}`$ clearly does not scale simply with the multiplicity measured in $`pp`$ collisions at the same energy. Rather, for a large range of impact parameter, the multiplicity scales approximately with the total multiplicity in $`e^+e^{}`$ annihilation at the same $`\sqrt{s}`$. Thus, it appears that the first few collisions per participant are sufficient to liberate as much energy for particle production as an $`e^+e^{}`$ reaction.

However, the rapid approach of $`N_{ch}/N_{part}/2`$ in central heavy ion collisions below $`\sqrt{s_{NN}}20`$ GeV towards the $`e^+e^{}`$ data complicates any simple geometric interpretation, as all of the heavy ion data compared are for a similar range of impact parameters. One feature that might point to why the particle yields at the AGS and SPS are perhaps “suppressed” relative to $`e^+e^{}`$ data (and even to $`pp`$ data at lower energies, as noted in Ref. na49 ) is the magnitude of the ratio of net baryons to pions in the system. This ratio, which scales approximately as $`N_{part}/N_{ch}`$, is $`O(50\%)`$ at AGS energies e895 , but is $`O(1\%)`$ at RHIC phenix\_lambda . In a thermal statistical approach Cleymans:2002mp , this reflects the decrease of the baryon chemical potential with increasing beam energy.

In conclusion, the PHOBOS experiment has measured the normalized charged-particle multiplicity $`N_{ch}/N_{part}/2`$ in Au+Au collisions at three RHIC energies as a function of the centrality of the collision. Above CERN SPS energies, the total multiplicity per participating nucleon pair, $`N_{ch}/N_{part}/2`$, in central events scales with $`\sqrt{s}`$ in the same way as $`e^+e^{}`$ data. This is suggestive of a universal mechanism of particle production in strongly-interacting systems, controlled mainly by the amount of energy available for particle production. This may be related to the multiple collisions suffered by each participant nucleon, which could substantially reduce the leading particle effect seen previously in $`pp`$ collisions. The weak centrality dependence for $`N_{ch}/N_{part}/2`$, reminiscent of “wounded nucleon” scaling, suggests that after the first few collisions per participant, the multiplicity per participant pair saturates near the value measured in $`e^+e^{}`$ reactions. Ultimately, the existence of simple scaling behavior with $`\sqrt{s_{\mathrm{eff}}}`$ and $`N_{part}`$ indicates stronger constraints on particle production than previously considered theoretically. Thus, these results may provide a new perspective on particle production in heavy ion collisions.

This work was partially supported by U.S. DOE grants DE-AC02-98CH10886, DE-FG02-93ER40802, DE-FC02-94ER40818, DE-FG02-94ER40865, DE-FG02-99ER41099, and W-31-109-ENG-38 as well as NSF grants 9603486, 9722606 and 0072204. The Polish groups were partially supported by KBN grant 2 PO3B 103 23. The NCU group was partially supported by NSC of Taiwan under contract NSC 89-2112-M-008-024.


text:
# Epidemic spreading in complex networks with degree correlations

## 1 Introduction

Complex networks arising in the modeling of many social, natural, and technological systems are often growing and self-organizing objects characterized by peculiar topological properties . Many empirical evidences have prompted that most of the times the resulting network’s topology exhibits emergent phenomena which cannot be explained by merely extrapolating the local properties of their constituents. Among these phenomena, two of them appear ubiquitous in growing networks. The first one concerns the small-world property , that is defined by an average path length—average distance between any pair of vertices—increasing very slowly (usually logarithmically) with the network size, $`N`$. The second one finds its manifestation in the scale-free (SF) degree distribution . This implies that the probability $`P(k)`$ that a vertex has degree $`k`$—it is connected to $`k`$ other vertices—is characterized by a power-law behavior $`P(k)k^\gamma `$, where $`2<\gamma 3`$ is a characteristic exponent.

The statistical physics approach  has been proved a very valuable tool for the understanding and modeling of these emergent phenomena in growing networks and has stimulated a more detailed topological characterization of several social and technological networks. In particular, it has been recognized that many of these networks possess non-trivial degree correlations . The use of statistical physics tools has also evidentiated several surprising results concerning dynamical processes taking place on top of complex networks. In particular, the absence of the percolation  and epidemic  thresholds in uncorrelated scale-free (SF) networks has hit the community because of its potential practical implications. The absence of the percolative threshold, indeed, prompts to an exceptional tolerance to random damages . On the other hand, the lack of any epidemic threshold makes SF networks the ideal media for the propagation of infections, bugs, or unsolicited information . While the study of uncorrelated complex networks is a fundamental step in the understanding of the physical properties of many systems, yet correlations may drastically change the obtained results, as several recent works addressing the effect of correlations in epidemic spreading have shown .

Here we want to provide a review of recent results concerning the epidemic spreading in random correlated complex networks. We will consider the susceptible-infected-susceptible (SIS) and susceptible-infected-removed (SIR) models  and we will provide an analytical description that includes two vertices degree correlations in the dynamical evolution of the infection prevalence. This will allow us to relate the presence or absence of epidemic threshold to the eigenvalue spectra of certain connectivity matrices of the networks. In particular, in the case of scale-free networks it is possible to show that for the SIS model, a SF degree distribution $`P(k)k^\gamma `$ with $`2<\gamma 3`$ in unstructured networks with any kind of degree correlations is a sufficient condition for a null epidemic threshold in the thermodynamic limit. For the SIR model, the same sufficient condition applies if the minimum possible degree of the graph is $`k_{min}2`$. The SIR model with $`k_{min}=1`$ has always a null threshold unless the SF behavior is originated only by minimum degree vertices. In other words, under very general conditions, the presence of two-point degree correlations does not alter the extreme weakness of SF networks to epidemic diffusion. The present results are derived from the divergence of the nearest neighbors average degree , which stems from the degree detailed balance condition , to be satisfied in all physical networks.

## 2 Correlated complex networks

In the following we shall consider unstructured undirected networks, in which all vertices within a given degree class can be considered statistically equivalent. Thus our results will not apply to structured networks in which a distance or time ordering can be defined; for instance, when the small-world property is not present . We will consider in particular the subset of undirected Markovian random networks , that are completely defined by the degree distribution $`P(k)`$ ant the conditional probability $`P(k^{}|k)`$ that a vertex of degree $`k`$ is connected to a vertex of degree $`k^{}`$. These two functions can have any form and are assumed to be normalized, i.e.

$$\underset{k}{}P(k)=\underset{k^{}}{}P(k^{}|k)=1.$$
(1)
The term “Markovian” refers to the fact that, in our approximation, all higher-order correlation functions can be obtained as a combination of the two fundamental functions $`P(k)`$ and $`P(k^{}|k)`$. In fact, this approximation represents a natural step towards a more complex description of real networks. In this sense, the Erdös-Rényi (ER) model  (defined starting from a set of $`N`$ vertices that are connected in pairs with an independent probability $`p`$) can be viewed as the zero-th order approximation, where the average degree is the only fixed parameter. The ER model is thus defined by the ensemble of all possible networks with a given average degree, but completely random at all other respects. The first-order approximation has been recently introduced by realizing that many real-world networked systems possess a more complex degree distribution than that predicted by the ER model (a Poisson distribution ). In this approximation the whole degree distribution, $`P(k)`$, is chosen as a constrain whereas the rest of properties are left at random . Even though this approximation represents a quantitative step forward, it only takes into account local properties and, therefore, it neglects possible correlations among different vertices, correlations that, on the other hand, are present in real networks . Thus, it is quite natural to introduce the second-order approximation as that with a fixed degree distribution, $`P(k)`$, and a two-point correlation function, $`P(k^{}|k)`$, but totally random to all other respects. As we will see, in this case the approximation must be carefully made since, due to the two-point correlation constrain, the fundamental functions $`P(k)`$ and $`P(k^{}|k)`$ must satisfy a peculiar detailed balance condition.

The degree distribution usually identifies two kinds of networks. A first class, which includes classical models of random graphs , is characterized by an exponentially bounded degree distribution. A second one refers to SF networks in which the degree distribution takes the form $`P(k)Ck^\gamma `$, usually with $`2<\gamma 3`$ . In this case the network shows a very high level of degree heterogeneity, signalled by unbounded degree fluctuations. Indeed, the second moment of the degree distribution, $`k^2`$, diverges in the thermodynamic limit $`k_c\mathrm{}`$, where $`k_c`$ is the maximum degree of the network. It is worth recalling that, in growing networks, $`k_c`$ is related to the network size $`N`$ as $`k_cN^{1/(\gamma 1)}`$ . Noticeably, it is the large degree heterogeneity of SF networks that is at the origin of their extreme weakness towards epidemic spreading.

### 2.1 Assortative and disassortative mixing

A direct study of the conditional probability $`P(k^{}|k)`$ in data from real networks usually yields results that are very noisy and difficult to interpret. In order to characterize the degree correlations, it is more useful to work with the average nearest neighbors degree (ANND) of the vertices of degree $`k`$ , defined by

$$\overline{k}_{nn}(k)\underset{k^{}}{}k^{}P(k^{}|k),$$
(2)
and to plot it as a function of the degree $`k`$. When two-point correlations are not present in the network, the conditional probability takes the form $`P^{nc}(k^{}|k)=k^{}P(k^{})/k`$, and the ANND reads $`\overline{k}_{nn}^{nc}(k)=k^2/k`$, which is independent on $`k`$. On the contrary, an explicit dependence of $`\overline{k}_{nn}(k)`$ on $`k`$ necessary implies the existence of non-trivial correlations, as often measured in real networks . For instance, in many social networks it is observed that vertices with high degree connect more preferably to highly connected vertices; a property referred to as “assortative mixing”. In this case, $`\overline{k}_{nn}(k)`$ is an increasing function of $`k`$. On the opposite side, many technological and biological networks show “disassortative mixing”; i.e. highly connected vertices are preferably connected to vertices with low degree and, consequently, $`\overline{k}_{nn}(k)`$ is a decreasing function of $`k`$. Then, the ANND provides an easy and powerful way to quantify two-point degree correlations, avoiding the fine details contained in the full conditional probability $`P(k^{}|k)`$.

### 2.2 Degree detailed balance condition

A key relation holding for all physical networks is that all edges must point from one vertex to another. This rather obvious observation turns out to have important implications since it forces the fundamental functions $`P(k)`$ and $`P(k^{}|k)`$ to satisfy the following degree detailed balance condition

$$kP(k^{}|k)P(k)=k^{}P(k|k^{})P(k^{}).$$
(3)
This condition states that the total number of edges pointing from vertices with degree $`k`$ to vertices of degree $`k^{}`$ must be equal to the total number of edges that point from vertices with degree $`k^{}`$ to vertices of degree $`k`$. This relation is extremely important since it constraints the possible form of the conditional probability $`P(k^{}|k)`$ once $`P(k)`$ is given. It may be surprising that such a detailed balance condition exists since, in fact, networks are the result of a multiplicative random process and, in principle, detailed balance conditions only holds for systems driven by additive noise . In fact, the usual detailed balance condition is the same as Eq. (3) without the prefactors $`k`$ and $`k^{}`$. These very prefactors account for the multiplicative character of the network and Eq. (3) can be viewed as a closure condition that guarantees the existence of the network. There is a simple way to derive this condition. Let $`N_k`$ be the number of vertices with degree $`k`$. Obviously, $`_kN_k=N`$ and, consequently, we can define the degree distribution as

$$P(k)=\frac{N_k}{N}.$$
(4)
The function $`P(k)`$ alone does not define completely the topology of the network, because it says nothing about how vertices are connected to each other. Thus, we need to define additionally the matrix of connections among vertices of different degrees. Let $`N_{k,k^{}}`$ be a symmetric matrix measuring the total number of edges between vertices of degree $`k`$ and vertices of degree $`k^{}`$, when $`kk^{}`$, and two times the number of self-connections, when $`k=k^{}`$. It is not difficult to realize that this matrix fulfills the identities

$`{\displaystyle \underset{k^{}}{}}N_{k,k^{}}`$ $`=`$ $`kN_k.`$ (5)
$`{\displaystyle \underset{k}{}}{\displaystyle \underset{k^{}}{}}N_{k,k^{}}`$ $`=`$ $`kN.`$ (6)
The first of this relations simply states that the number of edges emanating from all vertices of degree $`k`$ is $`kN_k`$, while the second indicates that the sum of all the vertices’s degrees is equal to two times the number of edges. The identity (6) allows us to define the joint probability

$$P(k,k^{})=\frac{N_{k,k^{}}}{kN},$$
(7)
where the symmetric function $`(2\delta _{k,k^{}})P(k,k^{})`$ is the probability that a randomly chosen edge connects two vertices of degrees $`k`$ and $`k^{}`$. The correlation coefficient computed from this joint probability has been recently used in Ref.  in order to quantify two-point degree correlations. The transition probability $`P(k^{}|k)`$, defined as the probability that an edge from a $`k`$ vertex points to a $`k^{}`$ vertex, can be easily written as

$$P(k^{}|k)=\frac{N_{k^{},k}}{kN_k}\frac{kP(k,k^{})}{kP(k)},$$
(8)
from where the detailed balance condition arises as a consequence of the symmetry of $`P(k,k^{})`$ (or $`N_{k,k^{}}`$).

From the degree detailed balance condition it is possible to derive some general exact results concerning the behavior of $`\overline{k}_{nn}(k,k_c)`$ and of $`\overline{k}_{nn}_{k_c}=_kP(k)\overline{k}_{nn}(k,k_c)`$ in SF networks . In these two functions we have now made explicit the $`k_c`$ dependence originated by the upper cut-off of the $`k`$ sum and that must be taken into account since it is a possible source of divergences in the thermodynamic limit. The results that we will derive will turn out to be fundamental in determining the epidemic spreading behavior in these networks. Let us start by multiplying by a $`k`$ factor both terms of Eq. (3) and summing over $`k^{}`$ and $`k`$. We obtain

$$k^2=\underset{k^{}}{}k^{}P(k^{})\underset{k}{}kP(k|k^{}),$$
(9)
In SF networks with $`2<\gamma <3`$ we have that the second moment of the degree distribution diverges as $`k^2k_c^{3\gamma }`$ ∗\*∗\*$``$For $`\gamma =3`$ the second moment diverges as $`k^2\mathrm{ln}k_c`$ but the argument, though more involved, is still valid.. We thus obtain from Eq. (9), using the definition (2),

$$\underset{k^{}}{}k^{}P(k^{})\overline{k}_{nn}(k^{},k_c)\frac{C}{(3\gamma )}k_c^{3\gamma },$$
(10)
where $`C`$ is the constant prefactor from the degree distribution. In the case of disassortative mixing , the function $`\overline{k}_{nn}(k^{},k_c)`$ is decreasing with $`k^{}`$ and, since $`k^{}P(k^{})`$ is an integrable function, the l.h.s. of Eq. (10) has no divergence related to the sum over $`k^{}`$. This implies that the divergence must be contained in the $`k_c`$ dependence of $`\overline{k}_{nn}(k^{},k_c)`$. In other words, the function $`\overline{k}_{nn}(k^{},k_c)\mathrm{}`$ for $`k_c\mathrm{}`$ in a non-zero measure set. In the case of assortative mixing, $`\overline{k}_{nn}(k^{},k_c)`$ is an increasing function of $`k^{}`$ and, depending on its rate of growth, there may be singularities associated to the sum over $`k^{}`$. Therefore, this case has to be analyzed in detail. Let us assume that the ANND grows as $`\overline{k}_{nn}(k^{},k_c)\alpha k^\beta `$, $`\beta >0`$, when $`k^{}\mathrm{}`$. If $`\beta <\gamma 2`$, again there is no singularity related to the sum over $`k^{}`$ and the previous argument for disassortative mixing holds. When $`\gamma 2\beta <1`$ there is a singularity coming from the sum over $`k^{}`$ of the type $`\alpha k_c^{\beta (\gamma 2)}`$. However, since Eq. (10) comes from an identity, the singularity on the l.h.s. must match both the exponent of $`k_c`$ and the prefactor on the r.h.s. In the case $`\gamma 2\beta <1`$, the singularity coming from the sum is not strong enough to match the r.h.s. of Eq. (10) since $`\beta (\gamma 2)<3\gamma `$. Thus, the function $`\overline{k}_{nn}(k^{},k_c)`$ must also diverge when $`k_c\mathrm{}`$ in a non-zero measure set. Finally, when $`\beta >1`$ the singularity associated to the sum is too strong, forcing the prefactor to scale as $`\alpha rk_c^{1\beta }`$ and the ANND as $`\overline{k}_{nn}(k^{},k_c)rk_c^{1\beta }k^\beta `$. It is easy to realize that $`r1`$, since the ANND cannot be larger than $`k_c`$. Plugging the $`\overline{k}_{nn}(k^{},k_c)`$ dependence into Eq. (10) and simplifying common factors, we obtain the identity at the level of prefactors

$$\frac{r}{2\gamma +\beta }=\frac{1}{3\gamma }.$$
(11)
Since $`\beta >1`$ and $`r<1`$, the prefactor in the l.h.s. of Eq. (11) is smaller than the one of the r.h.s. This fact implies that the tail of the distribution in the l.h.s. of Eq. (10) cannot account for the whole divergence of its r.h.s. This means that the sum is not the only source of divergences and, therefore, the ANND must diverge at some other point.

In summary, the function $`\overline{k}_{nn}(k^{},k_c)`$ must diverge when $`k_c\mathrm{}`$ in a non-zero measure set independently of the correlation behavior. The large $`k_c`$ singularity of the ANND can then be used to evaluate the quantity

$$\overline{k}_{nn}_N=\underset{k}{}P(k)\overline{k}_{nn}(k,k_c),$$
(12)
where we have explicitly considered $`k_c`$ as a growing function of the network size $`N`$. The r.h.s. of this equation is a sum of positive terms and diverges with $`k_c`$ at least as $`\overline{k}_{nn}(k,k_c)`$ both in the disassortative or assortative cases. In other words, all SF networks with $`2<\gamma 3`$ must present a $`\overline{k}_{nn}_N\mathrm{}`$ for $`N\mathrm{}`$. This statement is independent of the structure of the correlations present in the network. The quantity $`\overline{k}_{nn}_N`$ is particularly useful in model analysis and real data measurements. Degree correlation functions can be measured in several networks, but measurements are always performed in the presence of a finite $`k_c`$ that does not allow to exploit the singularity of the function $`\overline{k}_{nn}(k,k_c)`$. The most convenient way to exploit the infinite size singularity is therefore to measure the $`\overline{k}_{nn}_N`$ for increasing network sizes.

The quantity $`\overline{k}_{nn}_N`$ is very important in defining the properties of spreading processes in networks since it measures the number of individuals that can be infected in a few contagions. We shall discuss this point in relation to some specific epidemic models in the next sections.

## 3 The SIS model

As a first prototypical example for examining the properties of epidemic dynamics in SF networks we consider the susceptible-infected-susceptible (SIS) model , in which each vertex represents an individual of the population and the edges represent the physical interactions among which the infection propagates. Each individual can be either in a susceptible or infected state. Susceptible individuals become infected with probability $`\lambda `$ if at least one of the neighbors is infected. Infected vertices, on the other hand, recover and become susceptible again with probability one. A different recovery probability can be considered by a proper rescaling of $`\lambda `$ and the time. This model is conceived for representing endemic infections which do not confer permanent immunity, allowing individuals to go through the stochastic cycle susceptible $``$ infected $``$ susceptible by contracting the infection over and over again.

### 3.1 Uncorrelated homogeneous networks

In uncorrelated homogeneous networks, in which each vertex has more or less the same number of connections, $`kk`$, a general result states the existence of a finite epidemic threshold, separating an infected (endemic) phase, with a finite average density of infected individuals, from a healthy phase, in which the infection dies out exponentially fast . This is for instance the case of random networks with exponentially bounded degree distribution.

This result can be recovered by considering the dynamical evolution of the average density of infected individuals $`\rho (t)`$ (the prevalence) present in the network. The SIS model in homogeneous uncorrelated networks at a mean-field level is described by the following rate equation

$$\frac{d\rho (t)}{dt}=\rho (t)+\lambda k\rho (t)\left[1\rho (t)\right].$$
(13)
In this equation we have neglected higher order terms, since we are interested in the onset of the endemic state, close to the point $`\rho (t)0`$. Also, we have neglected correlations among vertices. That is, the probability of infection of a new vertex—the second term in Eq. (13)—is proportional to the infection rate $`\lambda `$, to the probability that a vertex is healthy, $`1\rho (t)`$, and to the probability that a edge in a healthy vertex points to an infected vertex. This last quantity, assuming the homogeneous mixing hypothesis††††††$`\mathrm{}`$The homogeneous mixing hypothesis  states that the force of the infection (the per capita rate of acquisition of the disease by the susceptible individuals) is proportional to the density of infected individuals., is approximated for homogeneous networks as $`k\rho (t)`$, i.e. proportional to the average number of connections and to the density of infected individuals. From Eq. (13) it can be proved the existence of an epidemic threshold $`\lambda _c=k^1`$ , such that $`\rho =0`$ if $`\lambda <\lambda _c`$, while $`\rho (\lambda \lambda _c)`$ if $`\lambda \lambda _c`$. In this context, it is easy to recognize that the SIS model is a generalization of the contact process model, widely studied as the paradigmatic example of an absorbing-state phase transition to a unique absorbing state .

### 3.2 Uncorrelated complex networks

For general complex networks, in which large degree fluctuations and correlations might be allowed, we must relax the homogeneous hypothesis made in writing Eq. (13) and work instead with the relative density $`\rho _k(t)`$ of infected vertices with given degree $`k`$; i.e. the probability that a vertex with $`k`$ edges is infected. Following Refs. , the rate equation for $`\rho _k(t)`$ can be written as

$$\frac{d\rho _k(t)}{dt}=\rho _k(t)+\lambda k\left[1\rho _k(t)\right]\mathrm{\Theta }_k(t).$$
(14)
In this case, the creation term is proportional to the spreading rate $`\lambda `$, the density of healthy sites $`1\rho _k(t)`$, the degree $`k`$, and the variable $`\mathrm{\Theta }_k(t)`$, that stands for the probability that an edge emanating from a vertex of degree $`k`$ points to an infected site. In the case of an uncorrelated random network, considered in Refs. , the probability that a edge points to a vertex with $`k`$ connections is equal to $`kP(k)/k`$ . This yields a $`\mathrm{\Theta }_k=\mathrm{\Theta }^{\mathrm{nc}}`$ independent of $`k`$ that reads as

$$\mathrm{\Theta }^{\mathrm{nc}}=\frac{1}{k}\underset{k^{}}{}k^{}P(k^{})\rho _k^{}(t).$$
(15)
Substituting the expression (15) into Eq. (14), it is possible to find the steady state solution where $`\mathrm{\Theta }^{\mathrm{nc}}`$ is now a function of $`\lambda `$ alone  by the following self-consistent equation:

$$\mathrm{\Theta }^{\mathrm{nc}}=\frac{1}{k}\underset{k}{}kP(k)\frac{\lambda k\mathrm{\Theta }^{\mathrm{nc}}}{1+\lambda k\mathrm{\Theta }^{\mathrm{nc}}}.$$
(16)
A non-zero stationary prevalence ($`\rho _k0`$) is obtained when the r.h.s. and the l.h.s. of Eq. (16), expressed as a function of $`\mathrm{\Theta }^{\mathrm{nc}}`$, cross in the interval $`0<\mathrm{\Theta }^{\mathrm{nc}}1`$, allowing a nontrivial solution. It is easy to realize that this corresponds to the inequality

$$\frac{d}{d\mathrm{\Theta }^{\mathrm{nc}}}\left(\frac{1}{k}\underset{k}{}kP(k)\frac{\lambda k\mathrm{\Theta }^{\mathrm{nc}}}{1+\lambda k\mathrm{\Theta }^{\mathrm{nc}}}\right)|_{\mathrm{\Theta }^{\mathrm{nc}}=0}1$$
(17)
being satisfied. The value of $`\lambda `$ yielding the equality in Eq. (17) defines the critical epidemic threshold $`\lambda _c`$, that is given for uncorrelated random networks by

$$\lambda _c^{\mathrm{nc}}=\frac{k}{k^2}.$$
(18)
In uncorrelated and infinite SF networks with $`\gamma 3`$, we therefore have $`k^2=\mathrm{}`$, and correspondingly $`\lambda _c^{\mathrm{nc}}=0`$. This is a very relevant result, signalling that the high heterogeneity of SF networks makes them extremely weak with respect to infections. These results have several implications in human and computer virus epidemiology .

### 3.3 Correlated complex networks

For a general network in which the degrees of the vertices are correlated, the above formalism is not correct, since we are not considering the effect of the degree $`k`$ into the expression for $`\mathrm{\Theta }_k`$. This effect can be taken into account, however, for Markovian networks, whose correlations are completely defined by the conditional probability $`P(k^{}|k)`$. In this case, it is easy to realize that the correct factor $`\mathrm{\Theta }_k`$ can be written as

$$\mathrm{\Theta }_k(t)=\underset{k^{}}{}P(k^{}|k)\rho _k^{}(t),$$
(19)
that is, the probability that an edge in a vertex of degree $`k`$ is pointing to an infected vertex is proportional to the probability that any edge points to a vertex with degree $`k^{}`$, times the probability that this vertex is infected, $`\rho _k^{}(t)`$, averaged over all the vertices connected to the original vertex. Eqs. (14) and (19) define together the mean-field equation describing the SIS model on Markovian complex networks,

$$\frac{d\rho _k(t)}{dt}=\rho _k(t)+\lambda k\left[1\rho _k(t)\right]\underset{k^{}}{}P(k^{}|k)\rho _k^{}(t).$$
(20)
It must be stressed that this equation is valid only if the network has no structure; i.e. the only relevant variable is the degree $`k`$. This implies that all vertices within a given degree class are statistically equivalent. This is not the case, for instance, of regular lattices or structured networks  in which a spatial ordering is constraining the connectivity among vertices. The exact solution of Eq. (20) can be difficult to find, depending on the particular form of $`P(k^{}|k)`$. However, it is possible to extract the value of the epidemic threshold by analyzing the stability of the steady-state solutions. Of course, the healthy state $`\rho _k=0`$ is one solution. For small $`\rho _k`$, we can linearize Eq. (20), getting

$$\frac{d\rho _k(t)}{dt}\underset{k^{}}{}L_{kk^{}}\rho _k^{}(t).$$
(21)
In the previous equation we have defined the Jacobian matrix $`𝐋=\{L_{kk^{}}\}`$ by

$$L_{kk^{}}=\delta _{kk^{}}+\lambda kP(k^{}|k),$$
(22)
where $`\delta _{kk^{}}`$ is the Kronecker delta symbol. The solution $`\rho _k=0`$ will be unstable if there exists at least one positive eigenvalue of the Jacobian matrix $`𝐋`$. Let us consider the connectivity matrix $`𝐂`$, defined by $`C_{kk^{}}=kP(k^{}|k)`$. Using the symmetry condition Eq. (3), it is easy to check that if $`v_k`$ is an eigenvector of $`𝐂`$, with eigenvalue $`\mathrm{\Lambda }`$, then $`P(k)v_k`$ is an eigenvector of the transposed matrix $`𝐂^T`$ with the same eigenvalue. From here it follows immediately that all the eigenvalues of $`𝐂`$ are real. Let $`\mathrm{\Lambda }_m`$ be the largest eigenvalue of $`𝐂`$. Then, the origin will be unstable whenever $`1+\lambda \mathrm{\Lambda }_m>0`$, which defines the epidemic threshold

$$\lambda _c=\frac{1}{\mathrm{\Lambda }_m},$$
(23)
above which the solution $`\rho _k=0`$ is unstable, and another non-zero solution takes over as the actual steady-state—the endemic state.

It is instructive to see how this general formalism recovers previous results , implicitly obtained for random uncorrelated networks. For any random network, in which there are no correlations among the degrees of the vertices, we have that the connectivity matrix is given by $`C_{kk^{}}^{\mathrm{nc}}=kP(k^{}/k)kk^{}P(k^{})/k`$, since the probability that an edge points to a vertex of connectivity $`k^{}`$ is proportional to $`k^{}P(k^{})`$. It is easy to check that the matrix $`\{C_{k^{}k}^{\mathrm{nc}}\}`$ has a unique eigenvalue $`\mathrm{\Lambda }_m^{\mathrm{nc}}=k^2/k`$, corresponding to the eigenvector $`v_k^{\mathrm{nc}}=k`$, from where we recover the now established result Eq. (18).

### 3.4 Correlated Scale-Free networks

The absence of an epidemic threshold in SF uncorrelated networks is an extremely important question that prompts to a possible weakness of many real-world networks. Particularly important is for the case of digital viruses spreading on the Internet  and sexually transmitted diseases diffusing on the web of sexual contacts . Both these networks show, in fact, SF properties that would imply the possibility of major epidemic outbreaks even for infections with a very low transmission rate. Immunization policies as well must be radically changed in the case that a network has a null epidemic threshold .

In view of the relevance of this framework, it is extremely important to study to which extent the presence of correlations are altering these results. The main question is therefore which conditions on the degree correlations of SF networks preserve the lack of a critical threshold. In the case of correlated networks, we have shown that the epidemic threshold is the inverse of the largest eigenvalue of the connectivity matrix $`𝐂`$. The absence of an epidemic threshold thus corresponds to a divergence of the largest eigenvalue of the connectivity matrix $`𝐂`$ in the limit of an infinite network size $`N\mathrm{}`$. In order to provide some general statement on the conditions for such a divergence we can make use of the Frobenius theorem for non-negative irreducible matrices . This theorem states the existence of the largest eigenvalue of any non-negative irreducible matrix, eigenvalue which is simple, real, positive, and has a positive eigenvector. In our case the matrix of interest is the connectivity matrix that is non-negative and irreducible. The irreducible property of the connectivity matrix is a simple consequence of the fact that all the degree classes in the network are accessible. That is, starting from the degree class $`k`$ it is always possible to find a path of edges that connects this class to any other class $`k^{}`$ of the network. If this is not the case it means that the network is built up of disconnected irreducible subnetworks and, therefore, we can apply the same line of reasoning to each subnetwork‡‡‡‡‡‡$`\mathrm{}`$Notice that being irreducible is not equivalent to being fully connected at the vertex to vertex level, but at the class to class level.. From the Frobenius theorem  it can be proved that the maximum eigenvalue, $`\mathrm{\Lambda }_m`$, of any non-negative irreducible matrix, $`A_{kk^{}}`$, satisfies the inequality

$$\mathrm{\Lambda }_m\underset{k}{\mathrm{min}}\frac{1}{\psi (k)}\underset{k^{}}{}A_{kk^{}}\psi (k^{}),$$
(24)
where $`\{\psi (k)\}`$ is any positive vector. In particular, by setting $`𝐀=𝐂^2`$ and $`\psi (k)=k`$ we obtain the inequality

$$\mathrm{\Lambda }_m^2\underset{k}{\mathrm{min}}\underset{k^{}}{}\underset{\mathrm{}}{}k^{}\mathrm{}P(\mathrm{}|k)P(k^{}|\mathrm{}).$$
(25)
This inequality relates the lower bound of the largest eigenvalue $`\mathrm{\Lambda }_m`$ to the degree correlation function and allows to find a sufficient condition for the absence of the epidemic threshold. By noting that $`_k^{}k^{}P(k^{}|\mathrm{})=\overline{k}_{nn}(\mathrm{},k_c)`$, we obtain the inequality

$$\mathrm{\Lambda }_m^2\underset{k}{\mathrm{min}}\underset{\mathrm{}}{}\mathrm{}P(\mathrm{}|k)\overline{k}_{nn}(\mathrm{},k_c).$$
(26)
The r.h.s. of this equation is a sum of positive terms, and by recalling the divergence of the ANND with $`k_c`$, we readily obtain that $`\mathrm{\Lambda }_m\mathrm{}`$ for all networks with diverging $`k^2`$ both in the disassortative or assortative cases§§§§§§$`\mathrm{\S }`$One may argue that, since we are calculating a minimum for $`k`$, if the transition probability $`P(\mathrm{}|k_0)`$ is zero at some point $`k_0`$, this minimum is zero. In this case it is possible to show that repeating the same argument with $`𝐂^3`$ instead of $`𝐂^2`$ provides us an inequality that avoids this problem, Ref. .. The divergence of $`\mathrm{\Lambda }_m`$ implies on its turn that the SIS epidemic threshold vanishes, in the thermodynamic limit, in all SF networks with assortative or disassortative mixing if the degree distribution has a diverging second moment; i.e. a SF degree distribution with exponent $`2<\gamma 3`$ is a sufficient condition for the absence of an epidemic threshold for the SIS model in unstructured networks with arbitrary two-point degree correlation function.

In physical terms, the absence of the epidemic threshold is related to the divergence of $`\overline{k}_{nn}_N`$ in SF networks. In homogeneous networks, where $`\overline{k}_{nn}_Nk`$, the epidemic spreading properties can be related to the average degree. In SF networks, however, the focus shifts to the possibility of infecting a large number of individuals in a finite number of contagions. The fact that an infected vertex has a very low degree is not very important if a hub of the network that provides connectivity to a large number of vertices is a few hops away. The infection can, in this case, very easily access a very large number of individuals in a short time. In SF networks is the ANND that takes into account more properly the level of degree fluctuations and thus rules the epidemic spreading dynamics.

It is worth stressing that the divergence of $`\overline{k}_{nn}_N`$ is ensured by the degree detailed balance condition alone, and it is a very general result holding for all SF networks with $`2<\gamma 3`$. On the contrary, the SF behavior with $`2<\gamma 3`$ is a sufficient condition for the lack of epidemic threshold only in networks with general two-point degree correlations and in absence of higher-order correlations. The reason is that the relation between the epidemic threshold and the maximum eigenvalue of the connectivity matrix only holds for these classes of networks. Higher order correlations, or the presence of an underlying metric in the network , can modify the rate equation at the basis of the SIS model and may invalidate the present discussion.

## 4 The SIR model

The susceptible-infected-removed model (SIR)  represents the other paradigmatic example of epidemic dynamics. Unlike in the SIS model, in this case infected individuals fall, after some random time, into a removed state where they cannot neither become infected again nor infect other individuals. This model tries to mimic real epidemics where individuals, after being infected, acquire permanent immunity or, in the worst case, die.

The SIR model is defined as follows. Individuals can only exist in three different states, namely, susceptible, infected, or removed. Susceptible individuals become infected with probability $`\lambda `$ if at least one of their neighbors is infected. On the other hand, infected individuals spontaneously fall in the removed state with probability $`\mu `$, which without lack of generality we set equal to unity. The main difference between both models is that whereas in the SIS, for $`\lambda >\lambda _c`$, the epidemics reaches a steady state, in the SIR the epidemics always dies and reaches eventually a state with zero density of infected individuals. The epidemic prevalence is thus defined in this case as the total number of infected individuals in the whole epidemic process.

### 4.1 Uncorrelated homogeneous networks

In a homogeneous system, the SIR model can be described in terms of the densities of susceptible, infected, and removed individuals, $`S(t)`$, $`\rho (t)`$, and $`R(t)`$, respectively, as a function of time. These three quantities are linked through the normalization condition

$$S(t)+\rho (t)+R(t)=1,$$
(27)
and they obey the following system of differential equations:

$`{\displaystyle \frac{dS}{dt}}`$ $`=`$ $`\lambda k\rho S,`$
$`{\displaystyle \frac{d\rho }{dt}}`$ $`=`$ $`\rho +\lambda k\rho S,`$ (28)
$`{\displaystyle \frac{dR}{dt}}`$ $`=`$ $`\rho .`$
These equations can be interpreted as follows: infected individuals decay into the removed class at a unity rate, while susceptible individuals become infected at a rate proportional to both the densities of infected and susceptible individuals. Here, $`\lambda `$ is the microscopic spreading (infection) rate, and $`k`$ is the number of contacts per unit time that is supposed to be constant for the whole population. In writing this last term of the equations we have assumed again, as in the case of the SIS model, the homogeneous mixing hypothesis ,

The most significant prediction of this model is the presence of a nonzero epidemic threshold $`\lambda _c`$ . If the value of $`\lambda `$ is above $`\lambda _c`$, $`\lambda >\lambda _c`$, the disease spreads and infects a finite fraction of the population. On the other hand, when $`\lambda `$ is below the threshold, $`\lambda <\lambda _c`$, the total number of infected individuals (the epidemic prevalence), $`R_{\mathrm{}}=lim_t\mathrm{}R(t)`$, is infinitesimally small in the limit of very large populations. In order to see this point, let us consider the set of equations (28). Integrating the equation for $`S(t)`$ with the initial conditions $`R(0)=0`$ and $`S(0)1`$ (i.e., assuming $`\rho (0)0`$, a very small initial concentration of infected individuals), we obtain

$$S(t)=e^{\lambda kR(t)}.$$
(29)
Combining this result with the normalization condition (27), we observe that the total number of infected individuals $`R_{\mathrm{}}`$ fulfills the following self-consistent equation:

$$R_{\mathrm{}}=1e^{\lambda kR_{\mathrm{}}}.$$
(30)
While $`R_{\mathrm{}}=0`$ is always a solution of this equation, in order to have a nonzero solution the following condition must be fulfilled:

$$\frac{d}{dR_{\mathrm{}}}\left(1e^{\lambda kR_{\mathrm{}}}\right)|_{R_{\mathrm{}}=0}1.$$
(31)
This condition is equivalent to the constraint $`\lambda \lambda _c`$, where the epidemic threshold $`\lambda _c`$ takes the value $`\lambda _c=k^1`$. Performing a Taylor expansion at $`\lambda =\lambda _c`$ it is then possible to obtain the epidemic prevalence behavior $`R_{\mathrm{}}(\lambda \lambda _c)`$ (valid above the epidemic threshold). From the point of view of the physics of non-equilibrium phase transition, it is easy to recognize that the SIR model is a generalization of the dynamical percolation model, that has been extensively studied in the context of absorbing-state phase transitions .

### 4.2 Uncorrelated complex networks

In order to take into account the heterogeneity induced by the presence of vertices with different degree, we consider the time evolution of the magnitudes $`\rho _k(t)`$, $`S_k(t)`$, and $`R_k(t)`$, which are the density of infected, susceptible, and removed vertices of degree $`k`$ at time $`t`$, respectively . These variables are connected by means of the normalization condition

$$\rho _k(t)+S_k(t)+R_k(t)=1.$$
(32)
Global quantities such as the epidemic prevalence can be expressed as an average over the various degree classes; for example, we define the total number of removed individuals at time $`t`$ by $`R(t)=_kP(k)R_k(t)`$, and the prevalence as $`R_{\mathrm{}}=lim_t\mathrm{}R(t)`$. At the mean-field level, for random uncorrelated networks, these densities satisfy the following set of coupled differential equations:

$`{\displaystyle \frac{d\rho _k(t)}{dt}}`$ $`=`$ $`\rho _k(t)+\lambda kS_k(t)\mathrm{\Theta }^{\mathrm{nc}}(t),`$ (33)
$`{\displaystyle \frac{dS_k(t)}{dt}}`$ $`=`$ $`\lambda kS_k(t)\mathrm{\Theta }^{\mathrm{nc}}(t),`$ (34)
$`{\displaystyle \frac{dR_k(t)}{dt}}`$ $`=`$ $`\rho _k(t).`$ (35)
The factor $`\mathrm{\Theta }^{\mathrm{nc}}(t)`$ represents the probability that any given edge points to an infected vertex and is capable of transmitting the disease. This quantity can be computed in a self-consistent way: The probability that an edge points to an infected vertex with degree $`k^{}`$ is proportional to $`k^{}P(k^{})`$. However, since the infected vertex under consideration received the disease through a particular edge that cannot be used for transmission anymore (since it points back to a previously infected individual) the correct probability must consider one less edge. Therefore,

$$\mathrm{\Theta }^{\mathrm{nc}}(t)=\frac{1}{k}\underset{k}{}(k1)P(k)\rho _k(t).$$
(36)
The equations (33), (34), (35), and (36), combined with the initial conditions $`R_k(0)=0`$, $`\rho _k(0)=\rho _k^0`$, and $`S_k(0)=1\rho _k^0`$, completely define the SIR model on any uncorrelated complex network with degree distribution $`P(k)`$. We will consider in particular the case of a homogeneous initial distribution of infected individuals, $`\rho _k^0=\rho ^0`$. In this case, in the limit $`\rho ^00`$, we can substitute $`\rho _k(0)0`$ and $`S_k(0)1`$. Under this approximation, Eqs. (34) and (35) can be directly integrated, yielding

$$S_k(t)=e^{\lambda k\varphi (t)},R_k(t)=_0^{\mathrm{}}\rho _k(\tau )𝑑t,$$
(37)
where we have defined the auxiliary function

$$\varphi (t)=_0^t\mathrm{\Theta }^{\mathrm{nc}}(\tau )𝑑\tau =\frac{1}{k}\underset{k}{}(k1)P(k)R_k(t).$$
(38)
In order to get a closed relation for the total density of infected individuals, it results more convenient to focus on the time evolution of the averaged magnitude $`\varphi (t)`$. To this purpose, let us compute its time derivative:

$$\frac{d\varphi (t)}{dt}=1\frac{1}{k}\varphi (t)\frac{1}{k}\underset{k}{}(k1)P(k)e^{\lambda k\varphi (t)},$$
(39)
where we have introduced the time dependence of $`S_k(t)`$ obtained in Eq. (37). Once solved Eq. (39), we can obtain the total epidemic prevalence $`R_{\mathrm{}}`$ as a function of $`\varphi _{\mathrm{}}=lim_t\mathrm{}\varphi (t)`$. Since $`R_k(\mathrm{})=1S_k(\mathrm{})`$, we have

$$R_{\mathrm{}}=\underset{k}{}P(k)\left(1e^{\lambda k\varphi _{\mathrm{}}}\right).$$
(40)
For a general $`P(k)`$ distribution, Eq. (39) cannot be generally solved in a closed form. However, we can still get useful information on the infinite time limit; i.e. at the end of the epidemics. Since we have that $`\rho _k(\mathrm{})=0`$, and consequently $`lim_t\mathrm{}\mathrm{d}\varphi (t)/\mathrm{d}t=0`$, we obtain from Eq. (39) the following self-consistent equation for $`\varphi _{\mathrm{}}`$:

$$\varphi _{\mathrm{}}=1\frac{1}{k}\frac{1}{k}\underset{k}{}(k1)P(k)e^{\lambda k\varphi _{\mathrm{}}}.$$
(41)
The value $`\varphi _{\mathrm{}}=0`$ is always a solution. In order to have a non-zero solution, the condition

$$\frac{d}{d\varphi _{\mathrm{}}}\left(1\frac{1}{k}\frac{1}{k}\underset{k}{}(k1)P(k)e^{\lambda k\varphi _{\mathrm{}}}\right)|_{\varphi _{\mathrm{}}=0}1$$
(42)
must be fulfilled. This relation implies

$$\frac{\lambda }{k}\underset{k}{}k(k1)P(k)1,$$
(43)
which defines the epidemic threshold

$$\lambda _c^{\mathrm{nc}}=\frac{k}{k^2k},$$
(44)
below which the epidemic prevalence is null, and above which it attains a finite value. It is interesting to notice that this is precisely the same value found for the percolation threshold in generalized networks . This is hardly surprising since, as it is well known , the SIR model can be mapped to a bond percolation process.

### 4.3 Correlated complex networks

In order to work out the SIR model in Markovian networks it is easier to consider the rate equations for the quantities $`N_k^I(t)`$ and $`N_k^R(t)`$, defined as the number of infected and removed individuals of degree $`k`$, present at time $`t`$, respectively. From this two quantities we can easily recover the densities $`\rho _k(t)`$ and $`R_k(t)`$ as

$$\rho _k(t)=\frac{N_k^I(t)}{N_k},R_k(t)=\frac{N_k^R(t)}{N_k},S_k(t)=1\rho _k(t)R_k(t),$$
(45)
where $`N_k=NP(k)`$ is the number of vertices with degree $`k`$. The rate equations for $`N_k^I(t)`$ and $`N_k^R(t)`$ are then given by

$`{\displaystyle \frac{dN_k^I(t)}{dt}}`$ $`=`$ $`N_k^I(t)+\lambda S_k(t)\mathrm{\Gamma }_k(t),`$ (46)
$`{\displaystyle \frac{dN_k^R(t)}{dt}}`$ $`=`$ $`N_k^I(t),`$ (47)
where we have defined the function

$$\mathrm{\Gamma }_k(t)\underset{k^{}}{}N_k^{}^I(t)(k^{}1)P(k|k^{})$$
(48)
In this case, the creation of new infected individuals—the second term in the r.h.s. of Eq. (46)—is proportional to the number of infected individuals of degree $`k^{}`$, $`N_k^{}^I(t)`$, the probability that a vertex of degree $`k`$ is susceptible, $`S_k(t)`$, and the average number of edges pointing from these infected vertices to vertices of degree $`k`$, $`(k^{}1)P(k|k^{})`$, all averaged for all the vertices of degree $`k^{}`$. This last term takes into account that one of the edges is not available for transmitting the infection, since it was used to infect the vertex considered. Dividing Eq. (46) by $`N_k`$ and making use of the detailed balance condition Eq. (3) we find the rate equations for the relative densities as

$`{\displaystyle \frac{dR_k(t)(t)}{dt}}`$ $`=`$ $`\rho _k(t),`$ (49)
$`{\displaystyle \frac{d\rho _k(t)}{dt}}`$ $`=`$ $`\rho _k(t)+\lambda kS_k(t)\mathrm{\Theta }_k(t),`$ (50)
where, for a Markovian network, the factor $`\mathrm{\Theta }_k(t)`$ takes the form

$$\mathrm{\Theta }_k(t)=\underset{k^{}}{}\rho _k^{}(t)\frac{k^{}1}{k^{}}P(k^{}|k)$$
(51)
Again, it must be stressed that no structure is allowed in the network in order for this equations to represent a valid formulation of the SIR model.

In order to extract information about the epidemic threshold, we proceed similarly to the SIS model, performing a linear stability analysis. For time $`t0`$, that corresponds to small $`\rho _k`$, $`R_k0`$ and $`S_k1`$, Eqs. (50) and (51) can be written as

$$\frac{d\rho _k(t)}{dt}\underset{k^{}}{}\stackrel{~}{L}_{kk^{}}\rho _k^{}(t),$$
(52)
where the Jacobian matrix $`\stackrel{~}{𝐋}=\{\stackrel{~}{L}_{kk^{}}\}`$ can be written as

$$\stackrel{~}{L}_{kk^{}}=\delta _{kk^{}}+\lambda \frac{k(k^{}1)}{k^{}}P(k^{}|k).$$
(53)
In order to infect a finite fraction of individuals, we need the solution $`\rho _k=0`$ to be unstable, which happens if there is at least one positive eigenvalue of the Jacobian matrix $`\stackrel{~}{𝐋}`$. Defining the matrix $`\stackrel{~}{𝐂}=\{\stackrel{~}{C}_{kk^{}}\}`$, with elements

$$\stackrel{~}{C}_{kk^{}}=\frac{k(k^{}1)}{k^{}}P(k^{}|k),$$
(54)
we know from the Frobenius theorem, since it is positive and provided that it is reducible at the degree class level, that it has a largest eigenvalue $`\stackrel{~}{\mathrm{\Lambda }}_m`$ that is real and positive. Thus, the solution $`\rho _k=0`$ of Eq. (52) is stable whenever $`1+\lambda \stackrel{~}{\mathrm{\Lambda }}_m<0`$. This relation defines the epidemic threshold for the SIR model in Markovian networks

$$\lambda =\frac{1}{\stackrel{~}{\mathrm{\Lambda }}_m}.$$
(55)
In the case of a random uncorrelated network, we have that $`\stackrel{~}{C}_{kk^{}}^{\mathrm{nc}}=k(k^{}1)P(k^{})/k`$. It can be easily seen that this matrix has a unique eigenvalue $`\stackrel{~}{\mathrm{\Lambda }}_m^{\mathrm{nc}}=k^2/k1`$, corresponding to the eigenvalue $`\stackrel{~}{v}_k^{\mathrm{nc}}=k`$, thus recovering the previous result Eq. (44) obtained for this kind of networks.

The relation between the SIR model and percolation in correlated complex networks can be closed by noticing that the relevant parameter in this last problem is the largest eigenvalue of the matrix $`𝐂^{\mathrm{perc}}=\{C_{kk^{}}^{\mathrm{perc}}\}`$, with elements $`C_{kk^{}}^{\mathrm{perc}}=(k^{}1)P(k^{}|k)`$, as pointed out in Ref. . It is easy to check that if $`v_k^{\mathrm{perc}}`$ is an eigenvector of $`𝐂^{\mathrm{perc}}`$ with eigenvalue $`\mathrm{\Lambda }`$, then $`\stackrel{~}{v}_k=kv_k^{\mathrm{perc}}`$ is an eigenvector of $`\stackrel{~}{𝐂}`$ with the same eigenvalue. Then, the eigenvalues of $`𝐂^{\mathrm{perc}}`$ and $`\stackrel{~}{𝐂}`$ coincide, yielding in this way the same description and an identical critical point.

### 4.4 Correlated Scale-Free networks

The discussion of the absence of epidemic threshold of the SIS in SF networks with any sort of degree correlations can be easily extended to the SIR model, taking again advantage of the Frobenius theorem. In this case, in the general inequality given by Eq. (24), we set $`𝐀=\stackrel{~}{𝐂}^2`$ and $`\psi (k)=k`$, obtaining

$$\stackrel{~}{\mathrm{\Lambda }}_m^2\underset{k}{\mathrm{min}}\underset{k^{}}{}\underset{\mathrm{}}{}(\mathrm{}1)P(\mathrm{}|k)(k^{}1)P(k^{}|\mathrm{}).$$
(56)
Given that $`_k^{}(k^{}1)P(k^{}|\mathrm{})=\overline{k}_{nn}(\mathrm{},k_c)1`$, the previous inequality reads

$$\stackrel{~}{\mathrm{\Lambda }}_m^2\underset{k}{\mathrm{min}}\underset{\mathrm{}}{}(\mathrm{}1)P(\mathrm{}|k)\left[\overline{k}_{nn}(\mathrm{},k_c)1\right].$$
(57)
As in the case of the SIS model, the divergence of the ANND with $`k_c`$ in the thermodynamic limit, ensures the divergence of the eigenvalue $`\stackrel{~}{\mathrm{\Lambda }}_m`$. Therefore, a SF degree distribution with diverging second moment is a sufficient condition for the absence of an epidemic threshold also for the SIR model if the minimum degree of the network is $`k_{min}2`$. The only instance in which we can have an infinite $`\overline{k}_{nn}(\mathrm{},k_c)`$ with a finite eigenvalue $`\stackrel{~}{\mathrm{\Lambda }}_m`$ is when the divergence of $`\overline{k}_{nn}`$ is accumulated in the degree $`k=1`$ and results canceled by the term $`\mathrm{}1`$ in Eq. (57). This situation happens when the SF behavior of the degree distribution is just due to vertices with a single edge that form star-like structures by connecting on a few central vertices. Explicit examples of this situation are provided in Refs. .

## 5 Conclusions

In this paper we have reviewed the analytical treatment of the epidemic SIS and SIR models in complex networks at different levels of approximation, corresponding to the different levels in which degree correlations can be taken into account. At the zero-th level, in which all the vertices are assumed to have the same degree (homogeneous networks), we observe the presence of an epidemic threshold, separating an active or endemic phase from an inactive or healthy phase, that is inversely proportional to the average degree $`k`$. At this level of approximation, both models render the same result, thus showing a high degree of universality. At the first order approximation level, in which vertices are allowed to have a different degree, drawn from a specified degree distribution $`P(k)`$, but are otherwise random, we obtain epidemic thresholds that are inversely proportional to the degree fluctuations $`k^2`$. The remarkable fact about this result is that the epidemic threshold vanishes for SF networks with characteristic exponent $`2<\gamma 3`$ in the limit of an infinitely large network. Finally, in the second order approximation level, in which degree correlations are explicitly controlled by the conditional probability $`P(k^{}|k)`$ that a vertex of degree $`k`$ is connected to a vertex of degree $`k^{}`$, our analysis yields that the epidemic threshold in the SIS and SIR models is inversely proportional to the largest eigenvalue of the connectivity matrices $`C_{kk^{}}=kP(k^{}|k)`$ and $`\stackrel{~}{C}_{kk^{}}=k(k^{}1)P(k^{}|k)/k^{}`$, respectively. In the case of the SIR model we recover the mapping with percolation at the level of correlations exclusively among nearest neighbor vertices. The analysis of the divergence of the average nearest neighbors degree $`\overline{k}_{nn}(k,k_c)`$ with the degree cut-off $`k_c`$ allows us to establish the general result that any SF degree distribution with diverging second moment is a sufficient condition for the vanishing of the epidemic threshold in the SIS model. The same sufficient condition holds in the SIR model with $`k_{min}2`$. The SIR model with $`k_{min}=1`$ always shows the absence of an epidemic threshold with the exception of the peculiar case in which the divergence of the average nearest neighbor degree is accumulated only on the nodes of minimum degree. These results have extremely important consequences, since they imply that correlations are not able to stop an epidemic outbreak in SF networks, in opposition to previous claims, and indicates that a reduction of epidemic incidence can only be obtained by means of carefully crafted immunization strategies , or trivially through finite size effects .

## Acknowledgments

This work has been partially supported by the European commission FET Open project COSIN IST-2001-33555. R.P.-S. acknowledges financial support from the Ministerio de Ciencia y Tecnología (Spain).