text:
# The spread of epidemic disease on networks

## I Introduction

Many diseases spread through human populations by contact between infective individuals (those carrying the disease) and susceptible individuals (those who do not have the disease yet, but can catch it). The pattern of these disease-causing contacts forms a network. In this paper we investigate the effect of network topology on the rate and pattern of disease spread.

Most mathematical studies of disease propagation make the assumption that populations are “fully mixed,” meaning that an infective individual is equally likely to spread the disease to any other member of the population or subpopulation to which they belong Bailey75 ; AM91 ; Hethcote00 . In the limit of large population size this assumption allows one to write down nonlinear differential equations governing, for example, numbers of infective individuals as a function of time, from which solutions for quantities of interest can be derived, such as typical sizes of outbreaks and whether or not epidemics occur. (Epidemics are defined as outbreaks that affect a non-zero fraction of the population in the limit of large system size.) Epidemic behavior usually shows a phase transition with the parameters of the model—a sudden transition from a regime without epidemics to one with. This transition happens as the “reproductive ratio” $`R_0`$ of the disease, which is the fractional increase per unit time in the number of infective individuals, passes though one.

Within the class of fully mixed models much elaboration is possible, particularly concerning the effects of age structure in the population, and population turnover. The crucial element however that all such models lack is network topology. It is obvious that a given infective individual does not have equal probability of infecting all others; in the real world each individual only has contact with a small fraction of the total population, although the number of contacts that people have can vary greatly from one person to another. The fully mixed approximation is made primarily in order to allow the modeler to write down differential equations. For most diseases it is not an accurate representation of real contact patterns.

In recent years a large body of research, particularly within the statistical physics community, has addressed the topological properties of networks of various kinds, from both theoretical and empirical points of view, and studied the effects of topology on processes taking place on those networks Strogatz01 ; AB02 . Social networks WS98 ; ASBS00 ; Liljeros01 ; GN02 , technological networks ABW98 ; AJB99 ; FFF99 ; Broder00 , and biological networks Jeong00 ; Jeong01 ; FW00 ; WM00 ; MS02  have all been examined and modeled in some detail. Building on insights gained from this work, a number of authors have pursued a mathematical theory of the spread of disease on networks MN00a ; KA01 ; PV01a ; MPV02 ; WSS01 ; Warren02 . This is also the topic of the present paper, in which we show that a large class of standard epidemiological models can be solved exactly on networks using ideas drawn from percolation theory.

The outline of the paper is as follows. In Section II we introduce the models studied. In Section III we show how percolation ideas and generating function methods can be used to provide exact solutions of these models on simple networks with uncorrelated transmission probabilities. In Section IV we extend these solutions to cases in which probabilities of transmission are correlated, and in Section V to networks representing some types of structured populations. In Section VI we give our conclusions.

## II Epidemic models and percolation

The mostly widely studied class of epidemic models, and the one on which we focus in this paper, is the class of susceptible/infective/removed or SIR models. The original and simplest SIR model, first formulated (though never published) by Lowell Reed and Wade Hampton Frost in the 1920s, is as follows. A population of $`N`$ individuals is divided into three states: susceptible (S), infective (I), and removed (R). In this context “removed” means individuals who are either recovered from the disease and immune to further infection, or dead. (Some researchers consider the R to stand for “recovered” or “refractory.” Either way, the meaning is the same.) Infective individuals have contacts with randomly chosen individuals of all states at an average rate $`\beta `$ per unit time, and recover and acquire immunity (or die) at an average rate $`\gamma `$ per unit time. If those with whom infective individuals have contact are themselves in the susceptible state, then they become infected. In the limit of large $`N`$ this model is governed by the coupled nonlinear differential equations Bailey75 :

$$\frac{\mathrm{d}s}{\mathrm{d}t}=\beta is,\frac{\mathrm{d}i}{\mathrm{d}t}=\beta is\gamma i,\frac{\mathrm{d}r}{\mathrm{d}t}=\gamma i,$$
(1)
where $`s(t)`$, $`i(t)`$, and $`r(t)`$ are the fractions of the population in each of the three states, and the last equation is redundant, since necessarily $`s+i+r=1`$ at all times. This model is appropriate for a rapidly spreading disease that confers immunity on its survivors, such as influenza. In this article we will consider only diseases of this type. Diseases that are endemic because they propagate on timescales comparable to or slower than the rate of turnover of the population, or because they confer only temporary immunity, are not well represented by this model; other models have been developed for these cases Hethcote00 .

The model described above assumes that the population is fully mixed, meaning that the individuals with whom a susceptible individual has contact are chosen at random from the whole population. It also assumes that all individuals have approximately the same number of contacts in the same time, and that all contacts transmit the disease with the same probability. In real life none of these assumptions is correct, and they are all grossly inaccurate in at least some cases. In the work presented here we remove these assumptions by a series of modifications of the model.

First, as many others have done, we replace the “fully mixed” aspect with a network of connections between individuals MN00a ; KA01 ; PV01a ; WSS01 ; Warren02 ; MPV02 ; SS88 ; Longini88 ; KM96 ; BMS97 . Individuals have disease-causing contacts only along the links in this network. We distinguish here between “connections” and actual contacts. Connections between pairs of individuals predispose those individuals to disease-causing contact, but do not guarantee it. An individual’s connections are the set of people with whom the individual may have contact during the time he or she is infective—people that the individual lives with, works with, sits next to on the bus, and so forth.

We can vary the number of connections each person has with others by choosing a particular degree distribution for the network. (Recall that the degree of a vertex in a network is the number of other vertices to which it is attached.) For example, in the case of sexual contacts, which can transmit STDs, the degree distribution has been found to follow a power-law form Liljeros01 . By placing the model on a network with a power-law degree distribution we can emulate this effect in our model.

Our second modification of the model is to allow the probability of disease-causing contact between pairs of individuals who have a connection to vary, so that some pairs have higher probability of disease transmission than others.

Consider a pair of individuals who are connected, one of whom $`i`$ is infective and the other $`j`$ susceptible. Suppose that the average rate of disease-causing contacts between them is $`r_{ij}`$, and that the infective individual remains infective for a time $`\tau _i`$. Then the probability $`1T_{ij}`$ that the disease will not be transmitted from $`i`$ to $`j`$ is

$$1T_{ij}=\underset{\delta t0}{lim}(1r_{ij}\delta t)^{\tau _i/\delta t}=\mathrm{e}^{r_{ij}\tau _i},$$
(2)
and the probability of transmission is

$$T_{ij}=1\mathrm{e}^{r_{ij}\tau _i}.$$
(3)
Some models, particularly computer simulations, use discrete time-steps rather than continuous time, in which case instead of taking the limit in Eq. (2) we simply set $`\delta t=1`$, giving

$$T_{ij}=1(1r_{ij})^{\tau _i},$$
(4)
where $`\tau `$ is measured in time-steps.

In general $`r_{ij}`$ and $`\tau _i`$ will vary between individuals, so that the probability of transmission also varies. Let us assume initially that these two quantities are iid random variables chosen from some appropriate distributions $`P(r)`$ and $`P(\tau )`$. (We will relax this assumption later.) The rate $`r_{ij}`$ need not be symmetric—the probability of transmission in either direction might not be the same. In any case, $`T_{ij}`$ is in general not symmetric because of the appearance of $`\tau _i`$ in Eqs. (3) and (4).

Now here’s the trick: because $`r_{ij}`$ and $`\tau _i`$ are iid random variables, so is $`T_{ij}`$, and hence the a priori probability of transmission of the disease between two individuals is simply the average $`T`$ of $`T_{ij}`$ over the distributions $`P(r)`$ and $`P(\tau )`$, which is

$$T=T_{ij}=1_0^{\mathrm{}}drd\tau P(r)P(\tau )\mathrm{e}^{r\tau }$$
(5)
for the continuous time case or

$$T=1_0^{\mathrm{}}dr\underset{\tau =0}{\overset{\mathrm{}}{}}P(r)P(\tau )(1r)^\tau $$
(6)
for the discrete case WSS01 . We call $`T`$ the “transmissibility” of the disease. It is necessarily always in the range $`0T1`$.

Thus the fact that individual transmission probabilities vary makes no difference whatsoever; in the population as a whole the disease will propagate as if all transmission probabilities were equal to $`T`$. We demonstrate the truth of this result by explicit simulation in Section III.5. It is this result that makes our models solvable. Cases in which the variables $`r`$ and $`\tau `$ are not iid are trickier, but, as we will show, these are sometimes solvable as well.

We note further that more complex disease transmission models, such as SEIR models in which there is an infected-but-not-infective period (E), are also covered by this formalism. The transmissibility $`T_{ij}`$ is essentially just the integrated probability of transmission of the disease between two individuals. The precise temporal behavior of infectivity and other variables is unimportant. Indeed the model can be generalized to include any temporal variation in infectivity of the infective individuals, and transmission can still be represented correctly by a simple transmissibility variable $`T`$, as above.

Now imagine watching an outbreak of the disease, which starts with a single infective individual, spreading across our network. If we were to mark or “occupy” each edge in the graph across which the disease is transmitted, which happens with probability $`T`$, the ultimate size of the outbreak would be precisely the size of the cluster of vertices that can be reached from the initial vertex by traversing only occupied edges. Thus, the model is precisely equivalent to a bond percolation model with bond occupation probability $`T`$ on the graph representing the community. The connection between the spread of disease and percolation was in fact one of the original motivations for the percolation model itself FH63 , but seems to have been formulated in the manner presented here first by Grassberger Grassberger83  for the case of uniform $`r`$ and $`\tau `$, and by Warren et al. WSS01 ; Warren02  for the non-uniform case.

In the next section we show how the percolation problem can be solved on random graphs with arbitrary degree distributions, giving exact solutions for the typical size of outbreaks, presence of an epidemic, size of the epidemic (if there is one), and a number of other quantities of interest.

## III Exact solutions on networks with arbitrary degree distributions

One of the most important results to come out of empirical work on networks is the finding that the degree distributions of many networks are highly right-skewed. In other words, most vertices have only a low degree, but there are a small number whose degree is very high Price65 ; AJB99 ; ASBS00 ; AB02 . The network of sexual contacts discussed above provides one example of such a distribution Liljeros01 . It is known that the presence of highly connected vertices can have a disproportionate effect on certain properties of the network. Recent work suggests that the same may be true for disease propagation on networks PV01a ; LM01 , and so it will be important that we incorporate non-trivial degree distributions in our models. As a first illustration of our method therefore, we look at a simple class of unipartite graphs studied previously by a variety of authors BC78 ; MR95 ; MR98 ; NSW01 ; NWS02 ; CEBH00 ; CEBH01 ; CNSW00 ; DMS01a ; ALPH01 , in which the degree distribution is specified, but the graph is in other respects random.

Our graphs are simply defined. One specifies the degree distribution by giving the properly normalized probabilities $`p_k`$ that a randomly chosen vertex has degree $`k`$. A set of $`N`$ degrees $`\left\{k_i\right\}`$, also called a degree sequence, is then drawn from this distribution and each of the $`N`$ vertices in the graph is given the appropriate number $`k_i`$ of “stubs”—ends of edges emerging from it. Pairs of these stubs are then chosen at random and connected together to form complete edges. Pairing of stubs continues until none are left. (If an odd number of stubs is by chance generated, complete pairing is not possible, in which case we discard one $`k_i`$ and draw another until an even number is achieved.) This technique guarantees that the graph generated is chosen uniformly at random from the set of all graphs with the selected degree sequence.

All the results given in this section are averaged over the ensemble of possible graphs generated in this way, in the limit of large graph size.

### III.1 Generating functions

We wish then to solve for the average behavior of graphs of this type under bond percolation with bond occupation probability $`T`$. We will do this using generating function techniques Wilf94 . Following Newman et al. NSW01 , we define a generating function for the degree distribution thus:

$$G_0(x)=\underset{k=0}{\overset{\mathrm{}}{}}p_kx^k.$$
(7)
Note that $`G_0(1)=_kp_k=1`$ if $`p_k`$ is a properly normalized probability distribution.

This function encapsulates all of the information about the degree distribution. Given it, we can easily reconstruct the distribution by repeated differentiation:

$$p_k=\frac{1}{k!}\frac{\mathrm{d}^kG_0}{\mathrm{d}x^k}|_{x=0}.$$
(8)
We say that the generating function $`G_0`$ “generates” the distribution $`p_k`$. The generating function is easier to work with than the degree distribution itself because of two crucial properties:

Powers If the distribution of a property $`k`$ of an object is generated by a given generating function, then the distribution of the sum of $`k`$ over $`m`$ independent realizations of the object is generated by the $`m^{\mathrm{th}}`$ power of that generating function. For example, if we choose $`m`$ vertices at random from a large graph, then the distribution of the sum of the degrees of those vertices is generated by $`[G_0(x)]^m`$.

Moments The mean of the probability distribution generated by a generating function is given by the first derivative of the generating function, evaluated at 1. For instance, the mean degree $`z`$ of a vertex in our network is given by

$$z=k=\underset{k}{}kp_k=G_0^{}(1).$$
(9)
Higher moments of the distribution can be calculated from higher derivatives also. In general, we have

$$k^n=\underset{k}{}k^np_k=\left[\left(x\frac{\mathrm{d}}{\mathrm{d}x}\right)^nG_0(x)\right]_{x=1}.$$
(10)
A further observation that will also prove crucial is the following. While $`G_0`$ above correctly generates the distribution of degrees of randomly chosen vertices in our graph, a different generating function is needed for the distribution of the degrees of vertices reached by following a randomly chosen edge. If we follow an edge to the vertex at one of its ends, then that vertex is more likely to be of high degree than is a randomly chosen vertex, since high-degree vertices have more edges attached to them than low-degree ones. The distribution of degrees of the vertices reached by following edges is proportional to $`kp_k`$, and hence the generating function for those degrees is

$$\frac{\underset{k}{}kp_kx^k}{_kkp_k}=x\frac{G_0^{}(x)}{G_0^{}(1)}.$$
(11)
In general we will be concerned with the number of ways of leaving such a vertex *excluding* the edge we arrived along, which is the degree minus 1. To allow for this, we simply divide the function above by one power of $`x`$, thus arriving at a new generating function

$$G_1(x)=\frac{G_0^{}(x)}{G_0^{}(1)}=\frac{1}{z}G_0^{}(x),$$
(12)
where $`z`$ is the average vertex degree, as before.

In order to solve the percolation problem, we will also need generating functions $`G_0(x;T)`$ and $`G_1(x;T)`$ for the distribution of the number of *occupied* edges attached to a vertex, as a function of the transmissibility $`T`$. These are simple to derive. The probability of a vertex having exactly $`m`$ of the $`k`$ edges emerging from it occupied is given by the binomial distribution $`\left(\genfrac{}{}{0pt}{}{k}{m}\right)T^m(1T)^{km}`$, and hence the probability distribution of $`m`$ is generated by

$`G_0(x;T)`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{k=m}{\overset{\mathrm{}}{}}}p_k\left({\displaystyle \genfrac{}{}{0pt}{}{k}{m}}\right)T^m(1T)^{km}x^m`$ (13)
$`=`$ $`{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}p_k{\displaystyle \underset{m=0}{\overset{k}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{m}}\right)(xT)^m(1T)^{km}`$
$`=`$ $`{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}p_k(1T+xT)^k`$
$`=`$ $`G_0(1+(x1)T).`$
Similarly, the probability distribution of occupied edges leaving a vertex arrived at by following a randomly chosen edge is generated by

$$G_1(x;T)=G_1(1+(x1)T).$$
(14)
Note that, in our notation

$`G_0(x;1)`$ $`=`$ $`G_0(x),`$ (15a)
$`G_0(1;T)`$ $`=`$ $`G_0(1),`$ (15b)
$`G_0^{}(1;T)`$ $`=`$ $`TG_0^{}(1),`$ (15c)
and similarly for $`G_1`$. ($`G_0^{}(x;T)`$ here represents the derivative of $`G_0(x;T)`$ with respect to its first argument.)

### III.2 Outbreak size distribution

The first quantity we will work out is the distribution $`P_s(T)`$ of the sizes $`s`$ of outbreaks of the disease on our network, which is also the distribution of sizes of clusters of vertices connected together by occupied edges in the corresponding percolation model. Let $`H_0(x;T)`$ be the generating function for this distribution:

$$H_0(x;T)=\underset{s=0}{\overset{\mathrm{}}{}}P_s(T)x^s.$$
(16)
By analogy with the previous section we also define $`H_1(x;T)`$ to be the generating function for the cluster of connected vertices we reach by following a randomly chosen edge.

Now, following Ref. NSW01, , we observe that $`H_1`$ can be broken down into an additive set of contributions as follows. The cluster reached by following an edge may be:

1. a single vertex with no occupied edges attached to it, other than the one along which we passed in order to reach it;
2. a single vertex attached to any number $`m1`$ of occupied edges other than the one we reached it by, each leading to another cluster whose size distribution is also generated by $`H_1`$.

We further note that the chance that any two finite clusters that are attached to the same vertex will have an edge connecting them together directly goes as $`N^1`$ with the size $`N`$ of the graph, and hence is zero in the limit $`N\mathrm{}`$. In other words, there are no loops in our clusters; their structure is entirely tree-like.

Using these results, we can express $`H_1(x;T)`$ in a Dyson-equation-like self-consistent form thus:

$$H_1(x;T)=xG_1(H_1(x;T);T).$$
(17)
Then the size of the cluster reachable from a randomly chosen starting vertex is distributed according to

$$H_0(x;T)=xG_0(H_1(x;T);T).$$
(18)
It is straightforward to verify that for the special case $`T=1`$ of 100% transmissibility, these equations reduce to those given in Ref. NSW01,  for component size in random graphs with arbitrary degree distributions. Equations (17) and (18) provide the solution for the more general case of finite transmissibility which applies to SIR models. Once we have $`H_0(x;T)`$, we can extract the probability distribution of clusters $`P_s(T)`$ by differentiation using Eq. (8) on $`H_0`$. In most cases however it is not possible to find arbitrary derivatives of $`H_0`$ in closed form. Instead we typically evaluate them numerically. Since direct evaluation of numerical derivatives is prone to machine precision problems, we recommend evaluating the derivatives by numerical contour integration using the Cauchy formula:

$$P_s(T)=\frac{1}{s!}\frac{\mathrm{d}^sH_0}{\mathrm{d}x^s}|_{x=0}=\frac{1}{2\pi \mathrm{i}}\frac{H_0(\zeta ;T)}{\zeta ^{s+1}}d\zeta ,$$
(19)
where the integral is over the unit circle MN00b . It is possible to find the first thousand derivatives of a function without difficulty using this method NSW01 . By this method then, we can find the exact probability $`P_s`$ that a particular outbreak of our disease will infect $`s`$ people in total, as a function of the transmissibility $`T`$.

### III.3 Outbreak sizes and the epidemic transition

Although in general we must use numerical methods to find the complete distribution $`P_s`$ of outbreak sizes from Eq. (19), we can find the mean outbreak size in closed form. Using Eq. (9), we have

$$s=H_0^{}(1;T)=1+G_0^{}(1;T)H_1^{}(1;T),$$
(20)
where we have made use of the fact that the generating functions are 1 at $`x=1`$ if the distributions that they generate are properly normalized. Differentiating Eq. (17), we have

$$H_1^{}(1;T)=1+G_1^{}(1;T)H_1^{}(1;T)=\frac{1}{1G_1^{}(1;T)},$$
(21)
and hence

$$s=1+\frac{G_0^{}(1;T)}{1G_1^{}(1;T)}=1+\frac{TG_0^{}(1)}{1TG_1^{}(1)}.$$
(22)
Given Eqs. (7), (12), (13), and (14), we can then evaluate this expression to get the mean outbreak size for any value of $`T`$ and degree distribution.

We note that Eq. (22) diverges when $`TG_1^{}(1)=1`$. This point marks the onset of an epidemic; it is the point at which the typical outbreak ceases to be confined to a finite number of individuals, and expands to fill an extensive fraction of the graph. The transition takes place when $`T`$ is equal to the critical transmissibility $`T_c`$, given by

$$T_c=\frac{1}{G_1^{}(1)}=\frac{G_0^{}(1)}{G_0^{\prime \prime }(1)}=\frac{\underset{k}{}kp_k}{_kk(k1)p_k}.$$
(23)
For $`T>T_c`$, we have an epidemic, or “giant component” in the language of percolation. We can calculate the size of this epidemic as follows. Above the epidemic threshold Eq. (17) is no longer valid because the giant component is extensive and therefore can contain loops, which destroys the assumptions on which Eq. (17) was based. The equation *is* valid however if we redefine $`H_0`$ to be the generating function only for outbreaks other than epidemic outbreaks, i.e., isolated clusters of vertices that are not connected to the giant component. These however do not fill the entire graph, but only the portion of it not affected by the epidemic. Thus, above the epidemic transition, we have

$$H_0(1;T)=\underset{s}{}P_s=1S(T),$$
(24)
where $`S(T)`$ is the fraction of the population affected by the epidemic. Rearranging Eq. (24) for $`S`$ and making use of Eq. (18), we find that the size of the epidemic is

$$S(T)=1G_0(u;T),$$
(25)
where $`uH_1(1;T)`$ is the solution of the self-consistency relation

$$u=G_1(u;T).$$
(26)
Results equivalent to Eqs. (22) to (26) were given previously in a different context in Ref. CNSW00, .

Note that it is not the case, even above $`T_c`$, that all outbreaks give rise to epidemics of the disease. There are still finite outbreaks even in the epidemic regime. While this appears very natural, it stands nonetheless in contrast to the standard fully mixed models, for which all outbreaks give rise to epidemics above the epidemic transition point. In the present case, the probability of an outbreak becoming an epidemic at a given $`T`$ is simply equal to $`S(T)`$.

### III.4 Degree of infected individuals

The quantity $`u`$ defined in Eq. (26) has a simple interpretation: it is the probability that the vertex at the end of a randomly chosen edge remains uninfected during an epidemic (i.e., that it belongs to one of the finite components). The probability that a vertex becomes infected via one of its edges is thus $`v=1T+Tu`$, which is the sum of the probability $`1T`$ that the edge is unoccupied, and the probability $`Tu`$ that it is occupied but connects to an uninfected vertex. The total probability of being uninfected if a vertex has degree $`k`$ is $`v^k`$, and the probability of having degree $`k`$ given that a vertex is uninfected is $`p_kv^k/_kp_kv^k=p_kv^k/G_0(v)`$, which distribution is generated by the function $`G_0(vx)/G_0(v)`$. Differentiating and setting $`x=1`$, we then find that the average degree $`z_{\mathrm{out}}`$ of vertices outside the giant component is

$$z_{\mathrm{out}}=\frac{vG_0^{}(v)}{G_0(v)}=\frac{vG_1(v)}{G_0(v)}z=\frac{u[1T+Tu]}{1S}z.$$
(27)
Similarly the degree distribution for an infected vertex is generated by $`[G_0(x)G_0(vx)]/[1G_0(v)]`$, which gives a mean degree $`z_{\mathrm{in}}`$ for vertices in the giant component of

$$z_{\mathrm{in}}=\frac{1vG_1(v)}{1G_0(v)}z=\frac{1u[1T+Tu]}{S}z.$$
(28)
Note that $`1S=G_0(u;T)u`$, since all coefficients of $`G_0(x;T)`$ are by definition positive (because they form a probability distribution) and hence $`G_0(x;T)`$ has only positive derivatives, meaning that it is convex everywhere on the positive real line within its domain of convergence. Thus, from Eq. (27), $`z_{\mathrm{out}}z`$. Similarly, $`z_{\mathrm{in}}z`$, and hence, as we would expect, the mean degree of infected individuals is always greater than or equal to the mean degree of uninfected ones. Indeed, the probability of a vertex being infected, given that it has degree $`k`$, goes as $`1v^k=1\mathrm{e}^{k\mathrm{log}(1/v)}`$, i.e., tends exponentially to unity as degree becomes large.

### III.5 An example

Let us now look at an application of these results to a specific example of disease spreading. First of all we need to define our network of connections between individuals, which means choosing a degree distribution. Here we will consider graphs with the degree distribution

$$p_k=\{\begin{array}{cc}0\hfill & \text{for }k=0\hfill \\ Ck^\alpha \mathrm{e}^{k/\kappa }\hfill & \text{for }k1\text{.}\hfill \end{array}$$
(29)
where $`C`$, $`\alpha `$, and $`\kappa `$ are constants. In other words, the distribution is a power-law of exponent $`\alpha `$ with an exponential cutoff around degree $`\kappa `$. This distribution has been studied before by various authors ASBS00 ; NSW01 ; NWS02 ; CNSW00 . It makes a good example for a number of reasons:

1. distributions of this form are seen in a variety of real-world networks ASBS00 ; Newman01a ;
2. it includes pure power-law and pure exponential distributions, both of which are also seen in various networks Price65 ; AJB99 ; FFF99 ; ASBS00 , as special cases when $`\kappa \mathrm{}`$ or $`\alpha 0`$;
3. it is normalizable and has all moments finite for any finite $`\kappa `$;

The constant $`C`$ is fixed by the requirement of normalization, which gives $`C=[Li_\alpha (\mathrm{e}^{1/\kappa })]^1`$ and hence

$$p_k=\frac{k^\alpha \mathrm{e}^{k/\kappa }}{Li_\alpha (\mathrm{e}^{1/\kappa })}\text{for }k1\text{,}$$
(30)
where $`Li_n(x)`$ is the $`n`$th polylogarithm of $`x`$.

We also need to choose the distributions $`P(r)`$ and $`P(\tau )`$ for the transmission rate and the time spent in the infective state. For the sake of easier comparison with computer simulations we use discrete time and choose both distributions to be uniform, with $`r`$ real in the range $`0r<r_{\mathrm{max}}`$ and $`\tau `$ integer in the range $`1\tau \tau _{\mathrm{max}}`$. The transmissibility $`T`$ is then given by Eq. (6). From Eq. (30), we have

$$G_0(x)=\frac{\underset{\alpha }{Li}(x\mathrm{e}^{1/\kappa })}{Li_\alpha (\mathrm{e}^{1/\kappa })}.$$
(31)
and

$$G_1(x)=\frac{\underset{\alpha 1}{Li}(x\mathrm{e}^{1/\kappa })}{xLi_{\alpha 1}(\mathrm{e}^{1/\kappa })}.$$
(32)
Thus the epidemic transition in this model occurs at

$$T_c=\frac{\underset{\alpha 1}{Li}(\mathrm{e}^{1/\kappa })}{Li_{\alpha 2}(\mathrm{e}^{1/\kappa })Li_{\alpha 1}(\mathrm{e}^{1/\kappa })}.$$
(33)
Below this value of $`T`$ there are only small (non-epidemic) outbreaks, which have mean size

$`s`$ $`=`$ $`1+`$
$`{\displaystyle \frac{T[\underset{\alpha 1}{Li}(\mathrm{e}^{1/\kappa })]^2}{Li_\alpha (\mathrm{e}^{1/\kappa })[(T+1)Li_{\alpha 1}(\mathrm{e}^{1/\kappa })TLi_{\alpha 2}(\mathrm{e}^{1/\kappa })]}}.`$
Above it, we are in the region in which epidemics can occur, and they affect a fraction $`S`$ of the population in the limit of large graph size. We cannot solve for $`S`$ in closed form, but we can solve Eqs. (25) and (26) by numerical iteration and hence find $`S`$.

In Fig. 1 we show the results of calculations of the average outbreak size and the size of epidemics from the exact formulas, compared with explicit simulations of the SIR model on networks with the degree distribution (30). Simulations were performed on graphs of $`N=\mathrm{100\hspace{0.17em}000}`$ vertices, with $`\alpha =2`$, a typical value for networks seen in the real world, and $`\kappa =5`$, 10, and 20 (the three curves in each panel of the figure). For each pair of the parameters $`\alpha `$ and $`\kappa `$ for the network, we simulated $`\mathrm{10\hspace{0.17em}000}`$ disease outbreaks each for $`(r,\tau )`$ pairs with $`r_{\mathrm{max}}`$ from 0.1 to 1.0 in steps of $`0.1`$, and $`\tau _{\mathrm{max}}`$ from 1 to 10 in steps of 1. Fig. 1 shows all of these results on one plot as a function of the transmissibility $`T`$, calculated from Eq. (6).

The figure shows two important things. First, the points corresponding to different values of $`r_{\mathrm{max}}`$ and $`\tau _{\mathrm{max}}`$ but the same value of $`T`$ fall in the same place and the two-parameter set of results for $`r`$ and $`\tau `$ collapses onto a single curve. This indicates that the arguments leading to Eqs. (5) and (6) are correct (as also demonstrated by Warren et al. WSS01 ; Warren02 ) and that the statistical properties of the disease outbreaks really do depend only on the transmissibility $`T`$, and not on the individual rates and times of infection. Second, the data clearly agree well with our analytic results for average outbreak size and epidemic size, confirming the correctness of our exact solution. The small disagreement between simulations and exact solution for $`s`$ close to the epidemic transition in the lower panel of the figure appears to be a finite size effect, due to the relatively small system sizes used in the simulations.

To emphasize the difference between our results and those for the equivalent fully mixed model, we compare the position of the epidemic threshold in the two cases. In the case $`\alpha =2`$, $`\kappa =10`$ (the middle curve in each frame of Fig. 1), our analytic solution predicts that the epidemic threshold occurs at $`T_c=0.329`$. The simulations agree well with this prediction, giving $`T_c=0.32(2)`$. By contrast, a fully mixed SIR model in which each infective individual transmits the disease to the same average number of others as in our network, gives a very different prediction of $`T_c=0.558`$.

## IV Correlated transmission probabilities

It is possible to imagine many cases in which the probabilities of transmission of a disease from an infective individual to those with whom he or she has connections are not iid random variables. In other words, the probabilities of transmission from a given individual to others could be drawn from different distributions for different individuals. This allows, for example, for cases in which the probabilities tend either all to be high or all to be low but are rarely a mixture of the two. In this section, we show how the model of Section III can be generalized to allow for this.

Suppose that the transmission rates $`r`$ for transmission from an infective individual $`i`$ to each of the $`k_i`$ others with whom they have connections are drawn from a distribution $`P_i(r)`$, which can vary from one individual to another in any way we like. Thus the a priori probability of transmission from $`i`$ to any one of his or her neighbors in the network is

$$T_i=1_0^{\mathrm{}}drd\tau P_i(r)P(\tau )\mathrm{e}^{r\tau }.$$
(35)
One could of course also allow the distribution from which the time $`\tau `$ is drawn to vary from one individual to another, although this doesn’t result in any functional change in the theory, so it would be rather pointless. In any case, the formalism developed here can handle this type of dependency perfectly well.

Following Eq. (13), we note that in the percolation representation of our model the distribution of the number of occupied edges leading from a particular vertex is now generated by the function

$`G_0(x;\left\{T_i\right\})`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=0}{\overset{N}{}}}{\displaystyle \underset{m=0}{\overset{k_i}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k_i}{m}}\right)T_i^m(1T_i)^{k_im}x^m`$ (36)
$`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=0}{\overset{N}{}}}(1+(x1)T_i)^{k_i}.`$
And similarly, the probability distribution of occupied edges leaving a vertex arrived at by following a randomly chosen edge is generated by

$$G_1(x;\left\{T_i\right\})=\frac{\underset{i}{}k_i(1+(x1)T_i)^{k_i1}}{_ik_i}.$$
(37)
Clearly these reduce to Eqs. (13) and (14) when $`T_i`$ is independent of $`i`$.

With these definitions of the basic generating functions, our derivations proceed as before. The complete distribution of the sizes of outbreaks of the disease, excluding epidemic outbreaks if there are any, is generated by

$$H_0(x;\left\{T_i\right\})=xG_0(H_1(x;\left\{T_i\right\});\left\{T_i\right\}),$$
(38)
where

$$H_1(x;\left\{T_i\right\})=xG_1(H_1(x;\left\{T_i\right\});\left\{T_i\right\}).$$
(39)
The average outbreak size when there is no epidemic is given by Eq. (22) as before, and the size of epidemics above the epidemic transition is given by Eqs. (25) and (26). The transition itself occurs when $`G_1^{}(1;\left\{T_i\right\})=1`$ and, substituting for $`G_1`$ from Eq. (37), we can also write this in the form

$$\underset{i=0}{\overset{N}{}}k_i[(k_i1)T_i1]=0.$$
(40)
In fact, it is straightforward to convince oneself that when the sum on the left-hand side of this equation is greater than zero epidemics occur, and when it is less than zero they do not.

For example, consider the special case in which the distribution of transmission rates $`P(r)`$ depends on the degree of the vertex representing the infective individual. One could imagine, for example, that individuals with a large number of connections to others tend to have lower transmission rates than those with only a small number. In this case $`T_i`$ is a function only of $`k_i`$ and hence we have

$`G_0(x;\left\{T_k\right\})`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=0}{\overset{N}{}}}(1+(x1)T_{k_i})^{k_i}`$ (41)
$`=`$ $`{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}p_k(1+(x1)T_k)^k,`$
and

$$G_1(x;\left\{T_k\right\})=\frac{\underset{k}{}kp_k(1+(x1)T_k)^{k1}}{_kkp_k},$$
(42)
where $`T_k`$ is the mean transmissibility for vertices of degree $`k`$.

One can also treat the case in which the transmissibility is a function of the number of connections which the individual being infected has. If the probability of transmission to an individual with degree $`k`$ is $`U_k`$, then we define

$`G_0(x;\left\{U_k\right\})`$ $`=`$ $`{\displaystyle \underset{k}{}}p_kx^k,`$ (43)
$`G_1(x;\left\{U_k\right\})`$ $`=`$ $`{\displaystyle \frac{\underset{k}{}kp_k[1+(x^{k1}1)U_k]}{_kkp_k}},`$ (44)
and then the calculation of cluster size distribution and so forth proceeds as before.

Further, one can solve the case in which probability of transmission of the disease depends on both the probabilities of giving it and catching it, which are arbitrary functions $`T_k`$ and $`U_k`$ of the numbers of connections of the infective and susceptible individuals. (This means that transmission from a vertex with degree $`j`$ to a vertex with degree $`k`$ occurs with a probability equal to the product $`T_jU_k`$.) The appropriate generating functions for this case are

$`G_0(x;\left\{T_k\right\},\left\{U_k\right\})`$ $`=`$ $`{\displaystyle \underset{k}{}}p_k(1+(x1)T_k)^k,`$ (45)
$`G_1(x;\left\{T_k\right\},\left\{U_k\right\})`$ $`=`$ (46)
$`{\displaystyle \frac{\underset{k}{}kp_k[1+((1+(x1)T_k)^{k1}1)U_k]}{_kkp_k}},`$
and indeed Eqs. (41) to (44) can be viewed as special cases of these equations when either $`T_k=1`$ or $`U_k=1`$ for all $`k`$. Note that $`G_0(x;\left\{U_k\right\})`$ and $`G_0(x;\left\{T_k\right\},\left\{U_k\right\})`$ are both independent of $`\left\{U_k\right\}`$, since the probability of a randomly-chosen infective individual having the disease is unity, regardless of the probability that they caught it in the first place.

As a concrete example of the developments of this section, consider the physically plausible case in which the transmissibility $`T`$ depends inversely on the degree of the infective individual: $`T_k=T_1/k`$. Then from Eq. (40) we find that there is epidemic behavior only if

$$T_1>\frac{z}{z1},$$
(47)
regardless of the degree distribution. Since $`T`$ lies strictly between zero and one however, this is impossible. In networks of this type, we therefore conclude that diseases cannot spread. Only if transmissibilities fall off slower than inversely with degree in at least some part of their range are epidemics possible. One plausible way in which this might happen is if $`T_k(T_0+k)^1`$. In this case it is straightforward to show that epidemics are possible for some degree distributions for some values of $`T_0`$.

Other extensions of the model are possible too. One area of current interest is models incorporating vaccination MN00a ; PV02 . Disease propagation on networks incorporating vaccinated individuals can be represented as a joint site/bond percolation process, which can also be solved exactly CNSW00 , both in the case of uniform independent vaccination probability (i.e., random vaccination of a population) and in the case of vaccination that is correlated with properties of individuals such as their degree (so that vaccination can be directed at the so-called core group of the disease-carrying network—those with the highest degrees).

## V Structured populations

The models we have studied so far have made use of simple unipartite graphs as the substrate for the spread of disease. These graphs may have any distribution we choose of the degrees of their vertices, but in all other respects are completely random. Many of the really interesting cases of disease spreading take place on networks that have more structure than this. Cases that have been studied previously include disease spreading among children who attend a common school and among patients in different wards of a hospital between whom pathogens are communicated by peripatetic caregivers Hyde01 . Here, we give just one example of disease spreading in a population with a very simple structure. The example we consider is the spread of a sexually transmitted disease. The important structural element of the population in this case is its division into men and women.

### V.1 Bipartite populations

Consider then a population of $`M`$ men and $`N`$ women, who have distributions $`p_j`$, $`q_k`$ of their numbers $`j`$ and $`k`$ of possibly disease-causing contacts with the opposite sex (connections in our nomenclature). In a recent study of 2810 respondents Liljeros et al. Liljeros01  recorded the numbers of sexual partners of men and women over the course of a year and found the distributions $`p_j`$, $`q_k`$ shown in Fig. 2. As the figure shows, the distributions appear to take a power-law form $`p_jj^{\alpha _m}`$, $`q_kk^{\alpha _f}`$, with exponents $`\alpha _m`$ and $`\alpha _f`$ that fall in the range $`3.1`$ to $`3.3`$ for both men and women <sup>1</sup><sup>1</sup>1One should observe that the network studied in Ref. Liljeros01,  is a cumulative network of actual sexual contacts—it represents the sum of all contacts over a specified period of time. Although this is similar to other networks of sexual contacts studied previously GAM89 ; Klovdahl94 , it is not the network required by our models, which is the instantaneous network of connections (not contacts—see Section II). While the network measured may be a reasonable proxy for the network we need, it is not known if this is the case.. (The exponent for women seems to be a little higher than that for men, but the difference is smaller than the statistical error on the measurement.)

We will assume that the disease of interest is transmitted primarily by contacts between men and women (true only for some diseases in some communities HY84 ), so that, to a good approximation, the network of contacts is bipartite, as shown in the inset of Fig. 2. That is, there are two types of vertices representing men and women, and edges representing connections run only between vertices of unlike kinds. With each edge we associate two transmission rates, one of which represents the probability of disease transmission from male to female, and the other from female to male. These rates are drawn from appropriate distributions as before, as are the times for which men and women remain infective. Also as before, however, it is only the average integrated probability of transmission in each direction that matters for our percolation model, so that we have two transmissibilities $`T_{mf}`$ and $`T_{fm}`$ for the two directions <sup>2</sup><sup>2</sup>2It is also worth noting that networks of sexual contacts observed in sociometric studies Klovdahl94  are often highly dendritic, with few short loops, indicating that the tree-like components of our percolating clusters may be, at least in this respect, quite a good approximation to the shape of real STD outbreaks..

We define two pairs of generating functions for the degree distributions of males and females:

$`f_0(x)`$ $`=`$ $`{\displaystyle \underset{j}{}}p_jx^j,f_1(x)={\displaystyle \frac{1}{\mu }}f_0^{}(x),`$ (48a)
$`g_0(x)`$ $`=`$ $`{\displaystyle \underset{k}{}}q_kx^k,g_1(x)={\displaystyle \frac{1}{\nu }}g_0^{}(x),`$ (48b)
where $`\mu `$ and $`\nu `$ are the averages of the two degree distributions, and are related by

$$\frac{\mu }{M}=\frac{\nu }{N},$$
(49)
since the total numbers of edges ending at male and female vertices are necessarily the same. Using these functions we further define, as before

$`f_0(x;T)`$ $`=`$ $`f_0(1+(x1)T),`$ (50a)
$`f_1(x;T)`$ $`=`$ $`f_1(1+(x1)T),`$ (50b)
$`g_0(x;T)`$ $`=`$ $`g_0(1+(x1)T),`$ (50c)
$`g_1(x;T)`$ $`=`$ $`g_1(1+(x1)T).`$ (50d)
Now consider an outbreak that starts at a single individual, who for the moment we take to be male. From that male the disease will spread to some number of females, and from them to some other number of males, so that after those two steps a number of new males will have contracted the disease, whose distribution is generated by

$$F_0(x;T_{mf},T_{fm})=f_0(g_1(x;T_{fm});T_{mf}).$$
(51)
For a disease arriving at a male vertex along a randomly chosen edge we similarly have

$$F_1(x;T_{mf},T_{fm})=f_1(g_1(x;T_{fm});T_{mf}).$$
(52)
And one can define the corresponding generating functions $`G_0`$ and $`G_1`$ for the vertices representing the females.

Using these generating functions, we can now calculate generating functions $`H_0`$ and $`H_1`$ for the sizes of outbreaks of the disease in terms either of number of women or of number of men affected. The calculation proceeds exactly as in the unipartite case, and the resulting equations for $`H_0`$ and $`H_1`$ are identical to Eqs. (17) and (18). We can also calculate the average outbreak size and the size of an epidemic outbreak, if one is possible, from Eqs. (22), (25), and (26). The average outbreak size for males, for example, is

$`s`$ $`=`$ $`1+{\displaystyle \frac{F_0^{}(1;T_{mf},T_{fm})}{1F_1^{}(1;T_{fm},T_{mf})}}`$ (53)
$`=`$ $`1+{\displaystyle \frac{T_{mf}T_{fm}f_0^{}(1)g_1^{}(1)}{1T_{mf}T_{fm}f_1^{}(1)g_1^{}(1)}}.`$
The epidemic transition takes place when $`F_1^{}(1;T_{mf},T_{fm})=1`$, or equivalently when

$$T_{mf}T_{fm}f_1^{}(1)g_1^{}(1)=1,$$
(54)
and hence the epidemic threshold takes the form of a hyperbola in $`T_{mf}`$$`T_{fm}`$ space:

$$T_{mf}T_{fm}=\frac{1}{f_1^{}(1)g_1^{}(1)}=\frac{\mu \nu }{_jj(j1)p_j_kk(k1)p_k}.$$
(55)
Note that this expression is symmetric in the variables describing the properties of males and females. Although we derived it by considering the generating function for males $`F_1`$, we get the same threshold if we consider $`G_1`$ instead. Eq. (53) is not symmetric in this way, so that the typical numbers of males and females affected by an outbreak may be different. On the other hand Eq. (55) involves the transmissibilities $`T_{mf}`$ and $`T_{fm}`$ only in the form of their product, and hence the quantities of interest are a function only of a single variable $`T_{mf}T_{fm}`$.

The generalizations of Section IV, where we considered transmission probabilities that vary from one vertex to another, are possible also for the bipartite graph considered here. The derivations are straightforward and we leave them as an exercise for the reader.

### V.2 Discussion

One important result that follows immediately from Eq. (55) is that if the degree distributions are truly power-law in form, then there exists an epidemic transition only for a small range of values of the exponent of the power law. Let us assume, as appears to be the case, that the exponents are roughly equal for men and women: $`\alpha _m=\alpha _f=\alpha `$. Then Eq. (55) tells us that the epidemic falls on the hyperbola $`T_{mf}T_{fm}=T_c^2`$, where

$$T_c=\frac{\zeta (\alpha 1)}{\zeta (\alpha 2)\zeta (\alpha 1)},$$
(56)
where $`\zeta (x)`$ is the Riemann $`\zeta `$-function. The behavior of $`T_c`$ as a function of $`\alpha `$ is depicted in Fig. 3. As the figure shows, if $`\alpha 3`$, $`T_c=0`$ and hence $`T_{mf}T_{fm}=0`$, which is only possible if at least one of the transmissibilities $`T_{mf}`$ and $`T_{fm}`$ is zero. As long as both are positive, we will always be in the epidemic regime, and this would clearly be bad news. No amount of precautionary measures to reduce the probability of transmission would ever eradicate the disease. (Lloyd and May LM01  have pointed out that a related result appears in the theory of fully mixed models, where a heterogeneous distribution of the infection parameter $`\beta `$ (see Eq. (1)) with a divergent coefficient of variation will result in the absence of an epidemic threshold. Pastor-Satorras and Vespignani PV01a  have made similar predictions using mean-field-like solutions for SIRS-type endemic disease models on networks with power-law degree distributions and a similar result has also been reported for percolation models by Cohen et al. CEBH00 .) Conversely, if $`\alpha >\alpha _c`$, where $`\alpha _c=3.4788\mathrm{}`$ is the solution of $`\zeta (\alpha 2)=2\zeta (\alpha 1)`$, we find that $`T_c=1`$ and hence $`T_{mf}T_{fm}=1`$, which is only possible if both $`T_{mf}`$ and $`T_{fm}`$ are 1. When either is less than 1 no epidemic will ever occur, which would be good news. Only in the small intermediate region $`3<\alpha <3.4788`$ does the model possess an epidemic transition. Interestingly, the real-world network measured by Liljeros et al. Liljeros01  appears to fall precisely in this region, with $`\alpha 3.2`$. If true, this would be both good and bad news. On the bad side, it means that epidemics can occur. But on the good side, it means that it is in theory possible to prevent an epidemic by reducing the probability of transmission, which is precisely what most health education campaigns attempt to do. The predicted critical value of the transmissibility is $`T_c=0.363\mathrm{}`$ for $`\alpha =3.2`$. Epidemic behavior would cease were it possible to arrange for the transmissibility to fall below this value.

Some caveats are in order here. The error bars on the values of the exponent $`\alpha `$ are quite large (about $`\pm 0.3`$ Liljeros01 ). Thus, assuming that the conclusion of a power-law degree distribution is correct in the first place, it is still possible that $`\alpha <3`$, putting us in the regime where there is always epidemic behavior regardless of the value of the transmissibility. (The error bars are also large enough to put us in the regime $`\alpha >\alpha _c`$ in which there are no epidemics. Empirical evidence suggests that the real world is not in this regime however, since epidemics plainly do occur.)

It is also quite possible that the distribution is not a perfect power law. Although the measured distributions do appear to have power-law tails, it seems likely that these tails are cut off at some point. If this is the case, then there will always be an epidemic transition at finite $`T`$, regardless of the value of $`\alpha `$. Furthermore, if it were possible to reduce the number of partners that the most active members of the network have, so that the cutoff moves lower, then the epidemic threshold rises, making it easier to eradicate the disease. Interestingly, the fraction of individuals in the network whose degree need change in order to make a significant difference is quite small. At $`\alpha =3`$, for instance, a change in the value $`\kappa `$ of the cutoff from $`\kappa =\mathrm{}`$ to $`\kappa =100`$ affects only 1.3% of the population, but increases the epidemic threshold from $`T_c=0`$ to $`T_c=0.52`$. In other words, targeting preventive efforts at changing the behavior of the most active members of the network may be a much better way of limiting the spread of disease than targeting everyone. (This suggestion is certainly not new, but our models provide a quantitative basis for assessing its efficacy.)

Another application of the techniques presented here is described in Ref. ANMS02, . In that paper we model in detail the spread of walking pneumonia (Mycoplasma pneumoniae) in a closed setting (a hospital) for which network data are available from observation of an actual outbreak. In this example, our exact solutions agree well both with simulations and with data from the outbreak studied. Furthermore, examination of the analytic solution allows us to make specific suggestions about possible new control strategies for M. pneumoniae infections in settings of this type.

## VI Conclusions

In this paper, we have shown that a large class of the so-called SIR models of epidemic disease can be solved exactly on networks of various kinds using a combination of mapping to percolation models and generating function methods. We have given solutions for simple unipartite graphs with arbitrary degree distributions and heterogeneous and possibly correlated infectiveness times and transmission probabilities. We have also given one example of a solution on a structured network—the spread of a sexually transmitted disease on a bipartite graph of men and women. Our methods provide analytic expressions for the sizes of both epidemic and non-epidemic outbreaks and for the position of the epidemic threshold, as well as network measures such as the mean degree of individuals affected in an epidemic.

Applications of the techniques described here are possible for networks specific to many settings, and hold promise for the better understanding of the role that network structure plays in the spread of disease.

###### Acknowledgements.

The author thanks Lauren Ancel, László Barabási, Duncan Callaway, Michelle Girvan, Catherine Macken, Jim Moody, Martina Morris, and Len Sander for useful comments. This work was supported in part by the National Science Foundation under grant number DMS–0109086.


format: text

text:
# RPA vs. exact shell-model correlation energies

## I introduction

One objective of nuclear theory, not reached yet, is an accurate method to compute the binding energy of nuclides globally–that is, from helium to the superheavies beyond uranium, and to the driplines relevant to nucleosynthesis. The most successful and widely referenced approaches are semiclassical, such as the finite-range droplet model MN95  and the extended Thomas-Fermi plus Strutinski-integral model AP95 . Exact, fully quantum mechanical ab initio microscopic calculations starting from experimental NN scattering have been successful benchmark ; Aeq8nuclei ; Aeq12nuclei  but are computationally too demanding for all but the lightest nuclei.

What about less ambitious but still microscopic approaches? The method of choice for a global description has been mean-field models such as Hartree-Fock (HF) with a phenomenological interaction (e.g. Skyrme) HFmass . For computational simplicity Hartree-Fock is easiest to apply to spherical nuclei, although advances in computers allow one to treat deformed HF as well; in that case the deformed HF state breaks rotational symmetry. Other symmetries that are broken by Hartree-Fock include translational invariance and (for Hartree-Fock-Bogoliubov) particle-number conservation. Mean-field theory thus leaves out both “kinematic” correlations, which restore broken symmetries and which must be corrected for approximately HFmass , and “dynamical” correlations, which arise from the residual two-body interaction.

The interacting shell model (SM), which diagonalizes an effective Hamiltonian in a basis of shell-model configurations, includes correlations, respects rotational invariance and other symmetries, and yields accurate binding energies, but at a computational price. Shell model calculations are currently limited to about $`10^{78}`$ basis states, which restricts one to light to medium nuclei ($`A<50`$) for full $`0\mathrm{}\omega `$ calculations. Furthermore, there is no global shell model interaction; instead one must derive or fit an effective interaction for each regime such as the $`sd`$-shell (with valence nucleons between magic numbers of 8 and 20) and the $`pf`$-shell (between magic numbers of 20 and 40).

The random phase approximation, or RPA, ring ; RoweRMP ; Marshalek  is an appealing middle ground because it adds correlations to the mean-field state and is numerically more tractable than full diagonalization in a SM basis; it also exactly separates out spurious modes induced by symmetry breaking. Because RPA implicitly treats 2-particle, 2-hole correlations in the ground state, one might imagine RPA to be an approximation to a more sophisticated microscopic model, a kind of cleverly truncated interacting shell model so to speak. Unfortunately, this image is erroneous because, unlike HF and the interacting SM, RPA is not a variational theory, and in fact, RPA can overestimate the ‘exact’ binding energy. (This of course supposes one has the ‘correct’ effective interaction in the first place, a big question we will not attempt to answer).

Surprisingly, tests of HF+RPA calculations of binding energies in the literature have been scant and all have been against simple models: the Lipkin-Meshkov-Glick (LMG) model LipkinNPA ; RowePR1968b , a schematic interaction in a small shell-model space RowePR1969  (here HF+RPA did poorly in describing the “exact” binding energy), a generalized 3-level LMG model bertschPRC2000  and a two-level pairing model bertschNPA2000 . Based particularly on the success of the latter two cases, Bertsch and Hagino advanced the idea of using the mean-field theory corrected through RPA for nuclear binding-energy systematics bertsch2000 .

This history motivates us to test the accuracy of RPA with respect to other microscopic calculations. We present a thorough study, comparing RPA to exact SM calculations in a large number of non-trivial systems.

The paper is organized as follows. Section II presents the model space and the interaction used, shortly reviews the approach and compares the RPA and exact correlation energies for several nuclei. We find that RPA reproduces well the ground state energy for a large number of nuclides. Because the HF+RPA occasionally misses the binding energy, we attempt in Section III to determine if the accuracy of the HF+RPA binding energy correlates with properties of the HF state, such as deformation geometry or the overlap of the HF state with the exact ground state. We find no such relationship. And because RPA is often said to “restore” broken symmetries, in section IV we compare the kinetic energy of the RPA zero modes (associated with symmetry restoration) with the “kinematic” correlation energy obtained by projecting the HF state onto good $`J`$. Finally, we summarize our conclusions in section V.

## II correlation energy in rpa

In the interacting shell model one restricts the model space to a limited set of valence fermion states outside an inert core, interacting via a residual two-body Hamiltonian. We consider nuclei throughout the $`sd`$-shell, with an 16O core and using the Wildenthal USD interaction wildenthal , and in the lower $`pf`$-shell, with a 40Ca core and using the modified KB3 interaction KB3 . For each nucleus investigated in this paper we are able to perform a full 0$`\mathrm{}\omega `$ diagonalization; in the same valence space we calculate the HF solution using separate Slater determinants for protons and neutrons. In this restricted model there is no mixing of radial degrees of freedom, the only degrees of freedom being angular. We emphasize that this is nonetheless a much more sophisticated testing ground for HF+RPA than the usual toy models.

Several textbooks discuss extensively the matrix formulation of RPA (see e.g., ring ), so we only summarize. One starts with a self-consistent mean-field solution, a Hartree-Fock Slater determinant, which can break symmetries of the Hamiltonian and which ignores correlations. The correlation energy is the additional binding energy that arises from particle-hole correlations. For example, if one performs a shell-model diagonalization in the full space, then the shell-model correlation energy, which we consider to be the exact correlation energy, is

$$E_{SM}^{corr}=E_{HF}E_{SM}.$$
(1)
To derive the correlation energy in RPA, one treats the energy surface around the Hartree-Fock minimum as a multi-dimensional harmonic oscillator, which allows one to compute excited states as 1-particle, 1-hole (1p-1h) configurations, as well as to compute corrections to the ground state energy. This last is approximately the correlation energy from an admixture of 2p-2h configurations. Because in the mapping to a multi-dimensional harmonic oscillator one has bosonized the energy surface, RPA does not rigorously enforce the Pauli principle, so RPA implicitly assumes the admixture of 2p-2h states into the ground state is small.

The resulting RPA matrix equation is well-known:

$$\left(\begin{array}{cc}A& B\\ B^{}& A^{}\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right)=\mathrm{\Omega }\left(\begin{array}{c}X\\ Y\end{array}\right),$$
(2)
with

$`A_{nj,mi}HF\left|[\widehat{c}_j^{}\widehat{c}_n,[\widehat{H},\widehat{c}_m^{}\widehat{c}_i]]\right|HF,`$
$`B_{nj,mi}HF\left|[[\widehat{H},\widehat{c}_n^{}\widehat{c}_j],\widehat{c}_m^{}\widehat{c}_i]\right|HF.`$ (3)
Most applications of RPA are to transition strengths and excitation energies, but one can also derive the RPA correlation energy:

$$E_{RPA}^{corr}=\frac{1}{2}\mathrm{Tr}(A)\frac{1}{2}\underset{i}{}\mathrm{\Omega }_i.$$
(4)
One can also compute the correlation energy in the response-function formulation of RPA ShimizuPRL2000 .

A particular feature of the RPA is the appearance of zero frequencies which correspond to conserved quantities, sometimes interpreted as restoration of broken symmetries RoweRMP ; Marshalek . In our restricted model space, the mean field can break the rotational symmetry; therefore, the number of zero RPA frequencies is related to the shape of the HF solution. Thus, when the HF solution is spherical, there are no zero modes; when axially symmetric, we get two zero eigenfrequencies, while for triaxial shapes we obtain three. The zero modes require careful treatment, and one can rewrite Eqn. (4):

$$E_{RPA}^{corr}=\underset{\nu (\mathrm{\Omega }_\nu 0)}{}\underset{mi}{}\mathrm{\Omega }_\nu |Y_{mi}^\nu |^2+\underset{\mu (\mathrm{\Omega }_\mu =0)}{}\frac{𝒫_\mu ^2}{2M_\mu },$$
(5)
where $`𝒫_\mu `$ are the canonical momentum operators corresponding to zero modes (for details, see ring ). Eqns. (4) and (5) are equivalent, but the latter explicitly segregates the contribution from zero modes, that is, from restored symmetries bertschPRC2000 ; bertschNPA2000 ; ring ; RoweRMP ; Marshalek . We compare the kinetic energy due to the zero modes against the correlation energy obtained through projection of the HF state onto good $`J`$ in section IV.

Ideally one would hope for a good agreement between $`E_{RPA}^{corr}`$ and $`E_{SM}^{corr}`$. Figures 1-5 compare the exact and RPA correlation energies for nuclides in the $`sd`$\- and $`pf`$-shells respectively. Note that it is not our intention at this moment to test against the experimental data. In the $`sd`$-shell all the nuclei are numerically tractable for a full diagonalization; for an easier comparison we grouped into even-even (Fig. 1), even-odd (Fig. 2), odd-odd (Fig. 3) and oxygen isotopes, that is, only neutrons active (Fig. 4). In the $`pf`$-shell the dimension of the SM basis rapidly increases with increasing the number of particles in the valence shell, limiting the number of exact cases we can hand; in Fig. 5(a) we present nuclides with both protons and neutrons in the active $`pf`$ space, while Fig. 5(b) contains calcium isotopes (neutrons only).

The general trend is a good agreement between $`E_{SM}^{corr}`$ and $`E_{RPA}^{corr}`$; in each group, except for oxygen, there is only one or two nuclides where RPA significantly overestimates or underestimates the correlation energy. In the $`sd`$-shell, if one excludes the oxygen results, the rms deviation for 41 nuclides is 870 keV; in the lower $`pf`$-shell, excluding calcium, the rms deviation is 480 keV for 11 nuclides. There was not a significant difference between even-even, odd-odd, and odd-even/even-odd nuclides. The single-species results were significantly worse, however: 1800 keV for 6 oxygen isotopes and 730 keV for 7 calcium isotopes. We search for potential explanations for this difference, as described in section III below, but found none.

Despite these limitations, based on the results presented in Figs. 1-5, one can assert that RPA seems to be a promising candidate for a global microscopic theory of the ground state energies, as advanced in bertsch2000 .

## III Analysis of RPA accuracy

Despite our overall good agreement, RPA sometimes fails to accurately reproduce the exact correlation energy. In this section we search for measures that correlate with the success or failure of RPA, which may in turn point towards possible solutions. For example, if RPA fails to describe pairing correlations, one may turn to Hartree-Fock-Bogoliubov + quasi-particle RPA. For another example, one might expect RPA, which implicitly assumes only small 2p-2h correlations, to be more accurate when the original HF state has a large overlap with the exact ground state wavefunction. If this were the case, one might correct through some self-consistent or renormalized RPA scheme.

To quantify our search use the ratio

$$R=\frac{E_{RPA}^{corr}}{E_{SM}^{corr}}.$$
(6)
$`R=1`$ when HF+RPA agrees with the exact shell-model binding energy.

We begin by considering how well RPA works with different geometries of the mean-field solution. In particular, how well does the RPA correlation energy account for (“restore”) broken symmetries? Table I lists the deformation geometry (spherical, prolate, oblate, triaxial) for various nuclides; we compute $`\beta `$, $`\gamma `$ by diagonalizing the mass quadrupole tensor of valence nucleons for the HF state deform\_method . $`R`$ is close to one for a wide range of deformation parameters. In fact, if anything RPA is less reliable for a spherical mean-field geometry, overestimating the correlation energy by as much as a factor of two. This result is surprising. In exact shell-model calculations, the excited states for these “spherical” nuclei are dominated by 1-particle 1-hole configurations, so one might expect RPA to be more successful than for “deformed” nuclei.

All of the spherical nuclei we found in the $`sd`$-shell were oxygen isotopes (valence neutrons only); this is not surprising as the proton-neutron interaction is well-known to induce deformation. To pursue this issue further, we took 28Si, which normally has an oblate mean-field solution, and lowered the $`0d_{5/2}`$ single-particle energy, until the mean-field solution became spherical (filled $`d_{5/2}`$ for both protons and neutrons). In Fig. 6 we plot $`R`$ as function of the $`d_{5/2}`$ single particle energy. Once again $`R`$ is closer to 1 for the deformed regime than for the spherical, or closed-shell, regime.

Now, as mentioned at the beginning of this section, one possible explanation for the poor results for spherical nuclides is inadequate treatment of pairing. Because the onset of pairing is non-perturbative one might not expect RPA to describe pairing correlations well. In order to test this idea, we computed the expectation value of the pairing Hamiltonian in both the exact shell model calculation and in HF+RPA JohnsonTBP . We found, surprisingly, better agreement for oxygen and calcium nuclides and considerably poorer results for nuclides with more deformation and better RPA correlations energies. Hence we cannot point conclusively to pairing as the culprit.

Another place RPA may stumble is in its implicit assumption that the admixture of 2p-2h correlations in the ground state is small. This shows up in two places, in the bosonization of the energy surface about the Hartree-Fock minimum, which does not enforce the Pauli principle, and, in the computation of the RPA A and B matrices, the replacement of the RPA wavefunction by the HF state in equation (3). To test this possibility, we calculate the overlap between the HF state and the exact wave function. This cannot be done by simply taking the projection

$$𝒪=HF|SM.$$
(7)
In our SM calculations, we took the basis states to have $`J_z=0`$ for $`A`$ even and $`J_z=1/2`$ for $`A`$ odd; the HF state is a superposition of states with good $`J_z`$ which are not restricted to the SM values. Therefore, any rotation of the Slater determinant can change the weight of $`J_z`$ components so that $`𝒪`$ is not invariant, although both $`E_{RPA}^{corr}`$ and $`E_{SM}^{corr}`$ are. We therefore define the overlap as

$$O=\mathrm{Max}\left(HF|SM\right),$$
(8)
where the maximum is taken when rotating about the $`y`$ axis. Note that when the exact ground state has $`J=0`$, $`𝒪`$ is independent of orientation.

Fig. 7 shows $`R`$ as function of the overlap $`O`$; here we have included only nuclides in the $`sd`$-shell, but the results for $`pf`$ are similar. $`O`$ varies from close to zero to close to one, but no compelling correlation with $`R`$ appears. That is to say, the accuracy of the RPA correlation energy does not seem predicated on how well the HF state approximates the exact ground state wavefunction.

A similar measure would be the amount of correlations in the RPA wavefunction, that is, $`N_{ph}^{RPA}=_\nu _i|Y_{mi}^\nu |^2`$. In principle, RPA is more accurate when $`N_{ph}^{RPA}`$ is small, but, as in Fig. 7, we find no correlation between $`N_{ph}^{RPA}`$ and $`R`$. Because of this, it is not not obvious that any renormalized, self-consistent, extended, etc., RPA scheme would in practice yield substantial improvements.

## IV Kinetic energy of spurious modes

One of the features that makes RPA appealing for correcting the mean-field results is its capacity to identify and separate out, as zero excitation energy, the modes associated with conserved quantities. Consequently, one can identify in Eq. (5) a term in the correlation energy corresponding to the kinetic energy of the spurious solution

$$KE_0=\underset{\mu (\mathrm{\Omega }_\mu =0)}{}\frac{𝒫_\mu ^2}{2M_\mu }.$$
(9)
On the other hand, the simplest way to correct the broken symmetries in the mean-field solution is to project it onto good quantum numbers. Table II presents the correlation energy obtained by projection of the HF Slater determinants onto good $`J`$, $`E_{proj}^{corr}`$, and the kinetic energy of the spurious modes given by Eq. (9); for comparison we have also included the exact correlation energy missing from the mean-field state. We conclude that $`KE_0`$ cannot be identified with the correlation energy obtained by direct restoration of symmetries obtained by projecting the HF state onto good quantum numbers. This may have relevance when if one goes to a multi-$`\mathrm{}\omega `$ space and wants to remove contamination by spurious center-of-mass motion.

## V summary and outlook

We have investigated the accuracy of HF+RPA for computing the nuclear binding energies, by comparing against the exact shell-model binding energies obtained through full $`0\mathrm{}\omega `$ diagonalization. In simpler models RowePR1968b ; bertschPRC2000 ; bertschNPA2000 , the RPA binding energy is generally, but not always RowePR1969 , satisfactory. We found this also to be true with our significantly more complicated shell-model Hamiltonians.

There were some non-trivial cases where RPA did not accurately estimate the exact binding energy; however, aside from the rather unhelpful observation that the inaccuracies tended to be worse for single-species cases, we could find no useful correlate with the accuracy of RPA. We did find, however, that RPA tended to yield very good results for deformed nuclides; despite this, however, the part of the RPA correlation energy identifiable with rotational kinetic energy had no relation with the energy gained from projection.

How could one make RPA more reliable? There are two obvious steps. The first is to attempt to build pairing directly into the mean-field state through Hartree-Fock-Bogolyubov and then computing the quasiparticle-RPA correlations energy. Because we found no relation between the accurate calculation of the expectation value of the pairing Hamiltonian and the accurate calculation of the total binding energy, however, we are far from confident that HFB+QRPA will yield significant gains.

The other possibility is to correct RPA for its violation of the Pauli principle as well as using the RPA ground state wavefunction function in computing A and B matrices in equation (3). The literature documents numerous proposals along these lines: self-consistent renormalization of the mean-field single-particle energies and the residual interaction RoweRMP ; scrpa ; treating the particle-hole amplitudes as variational parameters Dreizler ; rederivation of the RPA equations via the “number-operator method” to account for ground state correlations KleinNPA1991 ; and dressing the RPA phonons to account for ground state correlations DinhPRC2001 . None of these methods have become widely accepted, and they overlap each other to a large degree, although the literature does not clarify similarities and differences. As in the case of ungarnished RPA, these proposals for improving RPA have only been tested against toys such as the Lipkin model LipkinNPA . Perhaps by implementing these various proposals for full-scale shell model calculations one could compare their efficacy and practicality. This we leave to future work.

The U.S. Department of Energy supported this investigation through grant DE-FG02-96ER40985.