text:
# Qualms regarding “Dynamical Foundations of Nonextensive Statistical Mechanics ” by C. Beck (cond-mat/0105374)

## Abstract

The derivation of Student’s pdf from superstatistics is a mere coincidence due to the choice of the $`\chi ^2`$ distribution for the inverse temperature which is actually the Maxwell distribution for the speed. The difference between the estimator and the variance introduces a fluctuating temperature that is generally different than the temperature of the variance. Fluctuations in the energy of a composite system are transferred to fluctuations in the temperature via Lévy’s transform. The randomization of the operational parameter leading to a subordinated process requires the same number of degrees of freedom in the directing as in the original process. In place of Student’s pdf we find the subordinated process to Maxwell’s speed distribution has a Pareto distribution which is interpreted as the density of the length of the random velocity vector.

There is always a certain delight to see old ideas from a new perspective, but there is no delight to see professed ‘new’ ideas that are not new. The notion of ‘superstatistics’ grew out of an idea that if one began with a stationary Maxwell distribution and introduced a *randomized operational temperature*, then a new process could be derived which has something to do with Tsallis’ distribution , which he obtained as a stationary condition from maximizing a nonadditive entropy of degree-$`q`$ with respect to the energy constraint . The generalization to ‘superstatistics’ then consisted in replacing the exponential in the Maxwell distribution by $`e^{\beta H}`$, where $`H`$ is any generic hamiltonian, and using an arbitrary probability distribution for the conjugate $`\beta `$ one could obtain a new probability distribution for $`H`$ by integrating the product of the two probability distributions over all values of $`\beta `$. We have questioned the physical justification of this procedure elsewhere .

The idea of randomizing temperature is not new, as Beck would have us believe; it can be found in ref. . It is based on subordination, which can be found in Feller’s book . If $`X(t)`$ represents a brownian motion then by *randomizing the operational time* $`t`$, a variety of new processes can be derived. That is, if the displacements of brownian motion are governed by the probability density function (pdf),

$$f_{t_0}(x)=\frac{1}{\sqrt{2\pi t_0}}e^{x^2/2t_0},$$
(1)
where $`x`$ is the random variable and $`t_0`$ the parameter, then exchanging their roles leads to a Lévy pdf,

$$f_{x_0}(t)=\frac{x_0}{\sqrt{2\pi t^3}}e^{x_0^2/2t},$$
(2)
for the randomized time, $`t`$, where $`x_0`$ is the parameter. The two pdf are related by the transform $`x^2/t_0=x_0^2/t`$. The Cauchy process is then said to be subordinated to brownian motion

$$f_{x_0}(x)=_{\mathrm{}}^+\mathrm{}f_t(x)f_{x_0}(t)𝑑t=\frac{x_0}{2\pi }_{\mathrm{}}^+\mathrm{}\frac{1}{t^2}e^{(x^2+x_0^2)/2t}𝑑t=\frac{1}{\pi }\frac{x_0}{x_0^2+x^2},$$

where the two processes are evaluated at the same time.

The Maxwell speed, $`U(T)`$, in a single dimension, and the temperature, $`T`$, stand in the same relation as the displacement of a brownian motion process, $`X(t)`$, and the time, $`t`$ . The Cauchy distribution for the speed is therefore subordinated to the Maxwell distribution in one dimension .

Therefore, Beck  does not have to suppose that the fluctuations in inverse temperature, $`\beta `$ follow a $`\chi ^2`$ distribution given by his equation (2),

$$f_{\beta _0}(\beta )=\frac{1}{\mathrm{\Gamma }(n/2)}\left(\frac{n}{2\beta _0}\right)^{n/2}\beta ^{n/21}e^{n\beta /2\beta _0},$$
(3)
which we would prefer to write as

$$f_{u_0}(\beta )=\frac{1}{\mathrm{\Gamma }(n/2)}\left(\frac{u_0^2}{2}\right)^{n/2}\beta ^{n/21}e^{{\scriptscriptstyle \frac{1}{2}}\beta u_0^2},$$
(4)
on the strength of equipartition, $`u_0^2\beta _0=n`$. By what we shall show is an incorrect association of Tsallis’ distribution with the subordinated process, Beck finds $`n=2/(q1)`$, where $`q`$ is the degree of the nonadditive entropy, or, what is commonly known as the Tsallis exponent. A single mode of (3) exists at $`\beta =(12/n)\beta _0`$. The single mode therefore vanishes for all $`q2`$. We should expect a qualitative change in the behavior of the process for $`q<2`$ and $`q2`$. That no behavioral change has been predicted, already casts doubts on the association of half the number of the degrees of freedom with the inverse of $`(q1)`$.

It follows from the Lévy transform

$$\beta u_0^2=\beta _0u^2$$
(5)
that the $`n`$-dimensional Maxwell speed distribution

$$f_{\beta _0}(u)=\frac{2}{\mathrm{\Gamma }(n/2)}\left(\frac{\beta _0}{2}\right)^{n/2}u^{n1}e^{{\scriptscriptstyle \frac{1}{2}}\beta _0u^2},$$
(6)
corresponds to the $`\chi ^2`$ distribution for inverse temperature, (4), and not Beck’s one dimensional Maxwell pdf (6)\[eqn. (14) below\]. What has happened is that, in a space of $`n`$ dimensions, we began with a Markov process, $`U(\beta )`$, whose stationary transition probability is given by (6). A whole host of new processes can be derived by randomizing the temperature, or its inverse, $`\beta >0`$, corresponding to a new random variable $`𝔅(\beta )`$ whose pdf is (4). The process $`𝔅(\beta )`$ is called the ‘directing’ process \[6, p. 347\]\[5, p. 56\]. The pdf of the new process $`U(𝔅(\beta ))`$ will be given by

$$f_{u_0}(u)=_0^{\mathrm{}}f_\beta (u)f_{u_0}(\beta )𝑑\beta =\frac{2}{B(n/2,n/2)}\frac{(u/u_0)^{n1}}{(1+(u/u_0)^2)^n}\frac{1}{u_0},$$
(7)
where we have set $`\beta =\beta _0`$ so as to ensure that the two processes are in thermal equilibium, and $`B(,)`$ is the beta function. Contrary to what has been claimed , (7) is *not* Tsallis’ canonical probability distribution \[2, eqn (9)\]

$$f_{\beta _0}(u)=\frac{\mathrm{\Gamma }\left(\frac{n+1}{2}\right)}{\sqrt{n/2}\mathrm{\Gamma }(n/2)}\sqrt{\frac{\beta _0}{2\pi }}\left(1+\frac{\beta _0u^2}{n}\right)^{(n+1)/2}.$$
(8)
From symmetry considerations, all odd order moments of (8) are zero. In particular, the second moment is

$$\overline{u^2}=\frac{n\beta _0^1}{n2},$$
(9)
which does not reflect equipartition that follows from (6) at all \[*cf*., eqn (18) below \].

For a generic value of $`n`$, (7) is a special case of the inverse beta pdf

$$f(x)=\frac{1}{B(n/2,n/2)}\frac{x^{n/21}}{(1+x)^n},$$
(10)
where $`x=(u/u_0)^{{\scriptscriptstyle \frac{1}{2}}}`$. In economics, (10) is known as the Pareto pdf, since “it was thought (rather naively from a modern statistical standpoint) that income distributions should have a tail with a density $`Ax^\alpha `$ as $`x\mathrm{}`$, and (10) fulfills this requirement”\[6, p. 50\].

In $`n=1`$ dimension, (7) is the half-Cauchy pdf for a positive variate. *Just as the Cauchy process for the displacement is subordinated to brownian motion in one dimension when time is randomized, the Cauchy process for the kinetic energy is subordinated to the Maxwellian in a single dimension when the temperature, or its inverse, is randomized* .

In $`n=3`$ dimensions, (7) is \[5, p. 95\]

$$f_{n=3}(v)=\frac{8}{\pi }\frac{v^2}{(1+v^2)^3}.$$
(11)
This is the density of the length of a random velocity vector $`v=\sqrt{v_x^2+v_y^2+v_z^2}`$ in $`\mathrm{}^3`$ dimensions. It is related to the density

$$f_{n=1}(v)=\frac{1}{B(\frac{1}{2};m\frac{1}{2})}\frac{1}{(1+v^2)^m},$$
(12)
in a fixed direction, by \[6, p. 32\]

$$f_{n=3}(v)=vf_{n=1}^{}(v),$$
(13)
where the prime denotes differentiation with respect to the argument. The pdf (12) has $`2m1`$ degrees of freedom. For $`m=1`$, (12) reduces to the Cauchy pdf, while for $`m=2`$, implying that the distribution has $`3`$ degrees of freedom, we obtain (11) from (13). Hence, (11) is the density of the length $`v`$ of a random velocity vector in $`n=3`$ dimensions \[5, p. 95\].

Although (3) is valid for any $`n>0`$, Beck’s ‘conditional’ pdf for the speed, (6), is our (6) for $`n=1`$. For any $`n1`$ there is an incompatibility in the dimensions of the two distributions. This is responsible for the seemingly close appearance of the pdf of the subordinated process with the Tsallis distribution. In other words, if Beck’s (6) is meant to be the speed distribution, it does not reflect the dimensionality of the $`\chi ^2`$ pdf for $`\beta `$, (3). Rather, Beck finds that the subordinated process is governed by Student’s $`t`$ distribution, (8). The first factor tends to unity as $`n\mathrm{}`$, and for every fixed $`u`$

$$\frac{n+1}{2}\mathrm{log}\left(1+\frac{\beta _0u^2}{n}\right)\frac{1}{2}\beta _0u^2$$

so that in the limit as $`n\mathrm{}`$ we have

$$f_{\beta _0}(u)=\sqrt{\frac{\beta _0}{2\pi }}e^{{\scriptscriptstyle \frac{1}{2}}\beta _0u^2}.$$
(14)
This is the $`n=1`$ Maxwellian of (6) which is the conditional pdf (6) that Beck started out with!

What has happened is that the derivation of the Student pdf (8) *formally* resembles subordination, but has really nothing to do with it. Suppose that $`n+1`$ random variables $`u`$ and $`u_1,u_2,\mathrm{},u_n`$ are independent and identically distributed according to the normal distribution with zero mean and standard deviation $`1/\sqrt{\beta _0}`$. The distribution of $`u`$ is given by (14), while the distribution of the square root of the average of the sum of squares, $`v=\sqrt{\frac{1}{n}_{i=1}^nu_i^2}`$ is

$$f_{\beta _0}(v)=\frac{2}{\mathrm{\Gamma }(n/2)}\left(\frac{n\beta _0}{2}\right)^{n/2}v^{n1}e^{{\scriptscriptstyle \frac{1}{2}}n\beta _0v}.$$
(15)
On account of the independence of the random variables $`u`$ and $`v`$, their joint pdf will be given by

$$f_{\beta _0}(u)f_{\beta _0}(v)=\frac{2}{\mathrm{\Gamma }(n/2)}\sqrt{\frac{2}{\pi }}\left(\frac{n}{2}\right)^{n/2}\beta _0^{(n+1)/2}v^{n1}e^{{\scriptscriptstyle \frac{1}{2}}\beta _0(u^2+nv^2)}.$$
(16)
The probability that $`u/vt`$ is the integral of (16) over the region $`v>0`$ and $`utv`$. Introducing the transformation $`u=xy`$ and $`v=y`$, whose Jacobian is $`y`$, gives the distribution

$$F(t)=\frac{2}{\mathrm{\Gamma }(n/2)}\sqrt{\frac{2}{\pi }}\left(\frac{n}{2}\right)^{n/2}\beta _0^{(n+1)/2}_{\mathrm{}}^t𝑑x_0^{\mathrm{}}𝑑yy^ne^{{\scriptscriptstyle \frac{1}{2}}\beta _0n(1+x^2/n)y^2}.$$

Its derivative yields Student’s pdf

$$f(t)=\frac{1}{\sqrt{n\pi }}\frac{\mathrm{\Gamma }\left(\frac{n+1}{2}\right)}{\mathrm{\Gamma }(n/2)}\left(1+\frac{t^2}{n}\right)^{(n+1)/2}.$$
(17)
It is just a mere coincidence that the $`\chi ^2`$ pdf which Beck chose to represent fluctuations in the inverse temperature, (3), has the same form as Maxwell’s density for the speed $`\sqrt{\frac{1}{n}_{i=1}^nu_i^2}`$. Quite surprisingly, the Student distribution (17) is independent of the variance, $`1/\beta _0`$. Although this is important insofar as testing hypotheses of the mean of a population, because it is independendent of the variance, is concerned, it nevertheless implies that information has been lost \[*cf*., (21) and (23) below\].

The reason why Beck has $`\beta _0`$ as a parameter in his pdf (8) is that he obtained

$`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{1}{\mathrm{\Gamma }(n/2)}}\left({\displaystyle \frac{n}{2\beta _0}}\right)^{n/2}{\displaystyle _0^{\mathrm{}}}\beta ^{(n1)/2}e^{{\scriptscriptstyle \frac{1}{2}}\beta (u^2+n/\beta _0)}𝑑\beta `$
$`=`$ $`{\displaystyle \frac{1}{B(\frac{1}{2},n/2)}}\sqrt{{\displaystyle \frac{\beta _0}{n}}}\left(1+{\displaystyle \frac{\beta _0u^2}{n}}\right)^{(n+1)/2},`$
using (3) instead of using (4). Had he done so, he would have obtained

$$f_{u_0}(u)=\frac{1}{B(\frac{1}{2},n/2)}u_0^1\left(1+\frac{u^2}{u_0^2}\right)^{(n+1)/2}.$$

The two are the same if equipartition holds,

$$u_0=\sqrt{\underset{i=1}{\overset{n}{}}u_i^2}=\sqrt{\frac{n}{\beta _0}}.$$
(18)
In general, the estimator of the variance, $`\frac{1}{n}_{i=1}^nu_i^2`$, will be different from the variance $`\beta _0^1`$. The estimator is $`\beta ^1`$, and this temperature will usually be different from the temperature of the reservoir, $`\beta _0^1`$. In other words, the statistic for the Student distribution, $`t=u/\sqrt{_{i=1}^nu_i^2}`$, is not distributed as a normal random variable, as it would be had the statistic been given by $`t=u/\sqrt{n/\beta _0}`$. The actual variance can be completely unknown, and it suffices the sample variance alone in order to make statistical predictions. Nevertheless, we expect, that as the number of degrees of freedom increases the distribution of $`t`$ will be very similar to that of a standard normal variable (*i.e.*, it is a consistent estimator). Accordingly, we can identify the $`X_i^2`$ in Beck’s eqn (3), which are distributed as $`\chi ^2`$, with the kinetic energies of the individual particles, and the left-hand side should be $`n/\beta `$, and not $`\beta `$ itself.

The idea of ‘superstatistics’ is that the exponential pdf

$$f_{\beta _0}(E)=\frac{e^{\beta _0E}\rho (E)}{Z(\beta _0)}$$
(19)
is actually a conditional pdf, and specifying a *normalized* pdf for what was previously a mere parameter, say, by the $`\chi ^2`$ pdf

$$f_{E_0}(\beta )=\frac{1}{\mathrm{\Gamma }(n/2)}E_0^{n/2}\beta ^{n/21}e^{\beta E_0}.$$
(20)
leads to a new distribution of $`E`$ which is precisely the Tsallis pdf. However, (20) is no arbitrarily chosen pdf for the inverse temperature, but, rather the directing process which results when the parameter of the original distribution is randomized.

Consider the usual case of a power law for the structure function ,

$$\rho (E)=\frac{E^{m1}}{\mathrm{\Gamma }(m)},$$

which leads to equipartition because the partition function is $`Z(\beta )=\beta ^m`$. Then the process which is subordinated to the exponential pdf (19) has a pdf

$$f_{E_0}(E)=_0^{\mathrm{}}f_\beta (E)f_{E_0}(\beta )𝑑\beta =\frac{1}{B(m,n/2)}\frac{(E/E_0)^{m1}}{(1+E/E_0)^{m+n/2}}\frac{1}{E_0}$$
(21)
which is still the inverted beta pdf \[*cf*., eqn (7)\]. This is because the Lévy transform,

$$\beta _0E=\beta E_0,$$
(22)
in (19) produces a randomized inverse temperature which is distributed according to (20) with $`n=2m`$ degrees of freedom. Consequently, the inverse beta pdf (21) is the Fisher-Snedecor pdf with equal degrees of freedom, $`m=n/2`$, or what is commonly known as the Pareto distribution. The transformation to a new process by randomizing the operational temperature cannot change the number of degrees of freedom of the system.

‘F-superstatistics’ is not a directing process, but, rather, a subordinated process that is obtained from (4) and (6) by averaging the product of these densities over a common value of the kinetic energy, *viz*.

$$f_{\beta _0}(\beta )=_0^{\mathrm{}}f_{\beta _0}(u)f_u(\beta )𝑑u=\frac{1}{B(n/2,n/2)}\frac{(\beta /\beta _0)^{n/21}}{(1+\beta /\beta _0)^n}\frac{1}{\beta _0}.$$
(23)
On account of the necessity that both the original and directing processes have the same number of degrees of freedom, (23) is a symmetrical inverse beta distribution, which necessarily has *the same number of degrees of freedom*. Comparing (23) with (21), it is easily seen that they are interchangeable on the strength of the Lévy transformation (22). Reference to a Tsallis distribution in $`\beta `$ space is completely inappropriate.

The symmetry of the inverse beta pdf (23) means that we are dealing with a composite system comprised of two subsystems with the same number of degrees of freedom. The total energy $`E_t=E_0+E`$, formed from two subsystems with energies $`E`$ and $`E_0`$, is fixed. The inverse beta pdf (21) is thus transformed into the beta pdf, which can be written as the composition law \[5, p. 80\]

$$\rho (E_t)=_0^{E_t}\rho (E_tE)\rho (E)𝑑E=\frac{E_t^{2m1}}{\mathrm{\Gamma }(2m)}$$

for the structure function.

Thus, we conclude that the statement “Of course, other distribution functions $`f(\beta )`$ can also be considered which may lead to other generalized statistics”  is devoid of meaning when the original process is distributed according to (19). More importantly, as we have concluded elsewhere, “subordination can be considered as the probabilistic origin of power laws in physics ”\[5, p. 58\].


format: text

text:
# References

The Calculation of Critical Parameters in SU(2) Gauge Theory with Kouvel-Fisher Method

O.A. Mogilevsky

Bogolyubov Institute for Theoretical Physics, 14-B Metrologichna Str., Kiev 03143, Ukraine

## Abstract

We calculate the critical coupling $`4/g_c^2`$ and critical exponent $`\beta `$ for the order parameter in SU(2) lattice gauge theory by applying of the finite size scaling technique and the method proposed by Kouvel and Fisher for analysis of experimental data. In contrast to the standard finite size scaling approach, this method allows to determine simultaneously both $`g_c^2/4`$ and $`\beta `$ as two parameters of the linear fit to the Monte-Carlo data.

Lattice calculations by Monte-Carlo (MC) method have so far been the only method for studying the temperature phase transition in gauge theories from the first principles. By their nature, MC calculations are carried out on the finite lattices. This is a most serious drawback of a computer simulations, because no finite system can exhibit a true phase transition. Nevertheless, finite systems remind of phase transitions and systematic studies of these pseudo-transitions as function of system size may reveal information about the phase transition in the thermodynamic limit. One way to do this particularly in order to evaluate the critical exponents at the theory is to use finite size scaling (FSS).

FSS was proposed in statistical physics  and applied to high precision MC data of Ising  and other models. The validity of this method for gauge theories has been demonstrated through the investigation of critical properties of SU(2) lattice gauge theories at finite temperature \[5–9\].

In this letter we propose the method of analysis of MC data which is the combination of FSS technique and so-called Kouvel-Fisher method for analysis of experimental data in solid state physics . We consider MC data for SU(2) order parameter on the lattices $`8^3\times 4`$, $`12^3\times 4`$ and $`18^3\times 4`$. In contrast to the standard FSS approach, our method allows to determine the critical coupling $`4/g_\mathrm{c}^2`$ and critical exponent $`\beta `$ as two parameters of the linear fit.

The partition function for SU(2) gauge theory on $`N_\sigma ^3\times N_\tau `$ lattice is

$$Z=\underset{(\mu ,\nu )}{}dU_{\mu ,\nu }e^{S(U)},$$
(1)
where $`U_{\mu ,\nu }`$ are the SU(2) matrices on the link $`(\mu ,\nu )`$, $`S`$ is the standard Wilson action

$$S(U)=\frac{4}{g^2}\underset{p}{}\left(1\frac{1}{2}TrUp\right).$$
(2)
$`Up`$ is a product of link matrices around a plaquette. The temperature $`T`$ (volume $`V`$) is defined by size of the lattice in the time-like (space-like) direction ($`a`$ is a lattice spacing)

$$T=1/\left(N_\tau a\right),V=\left(N_\sigma a\right)^3.$$
(3)
The order parameter of the deconfinement phase transition on an infinite volume lattice is the expectation value of lattice average of the Polyakov loop

$$L=\frac{1}{N_\sigma ^3}\underset{𝐧}{}L_𝐧,$$
(4)
where

$$L_𝐧=\frac{1}{2}Tr\underset{\tau =1}{\overset{N_\tau }{}}U_{\tau ,𝐧;0}$$
(5)
and $`U_{\tau ,𝐧;0}`$ is the SU(2) link matrix at the point $`(𝐧,\tau )`$ in the time-like direction. On the finite lattice one considers as the order parameter $`|L|`$, since in this case $`L`$ is always zero due to system flips between the two ordered states.

According to FSS theory the scaling function for the order parameter on the lattice $`N_\sigma ^3\times N_\tau `$, $`N_\tau `$ fix is given by

$$N_\sigma ^{\frac{\beta }{\nu }}|L|=Q_L(x,g_1N_\sigma ^{y_1}),$$
(6)
where $`\beta `$, $`\nu `$ are the critical exponents for order parameter and correlation length correspondingly, $`g_1`$ is the irrelevant scaling field, $`y_11`$. The scaling variable $`x`$ is defined as

$$x=tN_\sigma ^{\frac{1}{\nu }},t=\frac{g_\mathrm{c}^2/4g^2/4}{g_\mathrm{c}^2/4}.$$
(7)
Our definition of reduce temperature $`t`$ differs from one used in \[7–10\]

$$t=\frac{4/g^24/g_\mathrm{c}^2}{4/g_\mathrm{c}^2}$$
(8)
and seems to be more natural because the coupling constant $`g^2/4`$ plays just the role of a ”temperature” of the effective theory which can be obtained from the partition function (1) by integration on all link variables except those for Polyakov loops. It is clear that this difference disappears in the limit $`g^2/4g_\mathrm{c}^2/4`$ .

At large $`N_\sigma `$ the scaling function $`Q_L`$ must behave as

$$Q_L(x)^\beta .$$
(9)
This behaviour may be observed only if $`x`$ becomes ”sufficiently large”. Deviations from (9) are in close vicinity of the critical point due to finite size rounding and further away from $`x=0`$ due to correction to scaling. Expanding $`Q_L`$ in (6) around $`x=0`$ at large $`N_\sigma `$ results in

$$N_\sigma ^{\frac{\beta }{\nu }}|L|=a+\frac{b}{N_\sigma ^{y_1}}+\left(c+\frac{d}{N_\sigma ^{y_1}}\right)x,$$
(10)
where $`a`$, $`b`$, $`c`$, $`d`$ are unknown constants. Dropping the term $`O(N_\sigma ^{y_1})`$ one can obtain

$$N_\sigma ^{\frac{\beta }{\nu }}|L|=a+cx,$$
(11)
what is exactly the approach of Refs. \[7–10\]. On the other hand at $`x`$ ”sufficiently large” the $`x`$-dependence is drastically changed. In Ref.  the scaling function $`Q_L`$ has been calculated in the form

$$Q_L=Ax^\beta \left[1+\frac{B}{N_\sigma ^{y_1}}x^{y_1\nu }\right]$$
(12)
by the data collapsing method. It has been shown that the correction to scaling term in (12) is negligibly small as compared to the leading term on the interval $`0,5<x<2`$. In both cases the SU(2) Monte-Carlo data for the order parameter have been used which were computed on $`N_\sigma ^3\times N_\tau `$ lattices with $`N_\sigma `$=8,12,18,26 and $`N_\tau =4`$ and have already reported on in Refs. \[5–10\]. The numerical value of the critical exponent $`\beta `$ has been found which is compatible with the corresponding value $`\beta =0,327`$ of the 3-$`d`$ Ising model .

These two procedures have two disadvantages: i) $`3d`$ Ising model value of $`\nu =0,631`$ has to be used; ii) the very accurate calculation of critical coupling $`g_\mathrm{c}^2/4`$ is necessary. The latter is determining by the fourth cumulant of $`|L|`$ as in Refs.  or by the $`\chi ^2`$-method as proposed in Ref. . Even small shift leads to the additional uncertainty in the determination of critical exponent $`\beta `$, because it depends very sensitively on the value inserted for $`g_\mathrm{c}^2/4`$.

We propose below the method for the determination of the critical exponent $`\beta `$ which does not require to know of $`g_\mathrm{c}^2/4`$ from the very beginning. No information about the critical exponent $`\nu `$ is also required. Following the Kouvel-Fisher idea which was used for analysis of the experimental data in solid state physics  we consider the function

$$K(x)=\frac{Q_L}{Q_L/x}.$$
(13)
According to the results of Ref.  in the region $`0,5<x<2`$ it is easy to find from (12)

$$K(x)=\frac{x}{\beta }\left[1+\frac{\mathrm{const}}{N_\sigma ^{y_1}}x^{y_1\nu }\right].$$
(14)
Dropping the correction to scaling term we have in this region with the same accuracy as in (11)

$$K(x)=\frac{x}{\beta }.$$
(15)
Taking into account the identity

$$\frac{Q_L}{x}=\frac{g_\mathrm{c}^2}{4}N_\sigma ^{\frac{(1\beta )}{\nu }}\frac{|L|}{(g^2/4)}$$
(16)
we have immediately

$$K(x)=\frac{x}{\beta }=\frac{N_\sigma ^{\frac{1}{\nu }}}{g_\mathrm{c}^2/4}\frac{|L|}{|L|/(g^2/4)}.$$
(17)
Finally, we obtain the simple but informative expression

$$\frac{|L|}{|L|/(g^2/4)}=\frac{1}{\beta }\left(\frac{g^2}{4}\frac{g_\mathrm{c}^2}{4}\right).$$
(18)
Eq. (18) means that for the lattices with $`N_\tau `$ fix and various $`N_\sigma `$ there are the intervals of $`g^2/4`$ where MC data for different $`N_\sigma `$ must arrange around the same straight line. The inverse slop of this line represents the critical exponent $`\beta `$ and the intersection with $`g^2/4`$ axis yields the critical coupling.

In order to verify this approach we have used the high precision MC data produced by Bielefeld group on the lattices $`8^3\times 4`$, $`12^3\times 4`$ and $`18^3\times 4`$ . The calculation of derivatives with respect to $`g^2/4`$ has been made by using the parabolic fit for $`|L|`$ (every parabolae includes five MC points). This fit does not smooth the statistical fluctuations, but connects the data continuously and with continuous first and second derivatives. The least minimal $`\chi ^2`$-fit of Eq. (18) ($`\chi ^2/N=0,62`$) includes 7 MC points in the region $`0.396<g^2/4<0.426`$ ($`2.345<4/g^2<2.527`$) for $`N_\sigma =8`$, 24 MC points in the region $`0.414<g^2/4<0.430`$ ($`2.325<4/g^2<2.415`$) for $`N_\sigma =12`$ and 16 MC points in the region $`0.427<g^2/4<0.433`$ ($`2.310<4/g^2<2.340`$) for $`N_\sigma =18`$. The result is

$$\beta =0.326\pm 0.003;4/g_\mathrm{c}^2=2.2988\pm 0.0007.$$
(19)
We see that the critical exponent $`\beta `$ is in excellent agreement with $`3d`$ Ising model value and the inverse critical coupling $`4/g_\mathrm{c}^2`$ coincides with that value which has been calculated in Refs. \[7–10\] by other method.

In conclusion the following comments are in order. We have shown that the Kouvel-Fisher modification of FSS gives the extremely useful method for the calculation of critical exponent $`\beta `$ in SU(2) lattice gauge theory. This method allows to investigate the finite lattice system at small but finite $`t=(g^2/4g_\mathrm{c}^2/4)/(g_\mathrm{c}^2/4)`$ and large $`N_\sigma `$ with infinite volume formulae (18). The critical coupling $`g_\mathrm{c}^2/4`$ and critical exponent $`\beta `$ are simultaneously determined as two parameters of the linear fit. It is clear that this approach may be applied to the calculation of all critical exponents in SU(2) gauge theory and to other models which undergo the second-order phase transition.


text:
# Observation of Reduced Three-Body Recombination in a Fermionized 1D Bose Gas

Permanent address: \] Laboratoire de Physique des Lasers (UMR 7538), Université Paris 13. Permanent address: \] Department of Physics, University of Maryland, College Park, MD. To whom correspondence should be addressed: \]trey@nist.gov

## Abstract

We investigate correlation properties of a one-dimensional interacting Bose gas by loading a magnetically trapped 87Rb Bose-Einstein condensate into a deep two-dimensional optical lattice. We measure the three-body recombination rate for both the BEC in the magnetic trap and the BEC loaded into the optical lattice. The recombination rate coefficient is a factor of seven smaller in the lattice, which we interpret as a reduction in the local three-body correlation function in the 1D case. This is a signature of correlation intermediate between that of the uncorrelated phase coherent 1D mean-field regime and the strongly correlated Tonks-Girardeau regime.

The majority of experiments with quantum degenerate gases have been performed in the weakly interacting limit, on Bose-Einstein condensates (BECs) characterized by long-range phase coherence and well described by the mean-field Gross-Pitaevskii (GP) equation stringari . While the success of the GP equation in accounting for many experimental results has been spectacular, it has also led to the search for physics beyond mean-field theory. As in condensed matter physics, there is now great interest in highly correlated systems, where mean-field approaches are inapplicable. Progress toward such correlated systems includes the recent observation of the Mott-insulator transition in BECs loaded into optical lattices Greiner02  and the use of Feshbach resonances to increase interactions between atoms bosanova . Here we present evidence of strong correlations in a 1D degenerate Bose gas as reflected in a reduction of three-body recombination.

The role of fluctuations and correlations in Bose gases increases with reduced dimension. In homogenous systems, BEC is only possible in 3D. In 2D, a Kosterlitz-Thouless transition occurs, and in 1D there is no finite temperature transition phasetransition . By contrast, BEC is possible in 1, 2 and 3D for trapped systems Bagnato . Trapped 1D systems with $`\delta `$-function repulsive interactions are particularly interesting, in that for high density the ground state is a condensate, while in the low density limit the ground state is a highly correlated state known as a Tonks gas LL . This ground state is an example of “fermionization,” where the repulsive interactions mimic the Pauli exclusion principle. Indeed, the low energy excitation spectrum and the correlation functions are identical to those of non-interacting fermions, and the many-body wave function of the Bose gas is equal to the absolute value of the fermionic wave function girardeau .

For such a 1D Bose gas, the degree of correlation depends on the ratio between two energies: the repulsive energy of uncorrelated atoms at a given density, $`E_{\text{unc}}=gn_{\text{1D}}`$, and the quantum kinetic energy needed to correlate particles by localizing them with respect to each other on the order of the mean inter-particle distance $`d`$, $`E_{\text{cor}}=\mathrm{}^2/2md^2`$. Here, $`g`$ is the strength of the $`\delta `$-function interaction, $`m`$ is the atomic mass, and $`n_{\text{1D}}=1/d`$ is the 1D density. A single parameter $`\gamma E_{\text{unc}}/4E_{\text{cor}}=mg/2n_{1\text{D}}\mathrm{}^2`$ entirely characterizes a homogeneous 1D gas with repulsive short range interactions. Ref. LL  provides the exact eigenstate solutions for all $`\gamma `$. Recently, these results have been extendend to harmonically trapped gases, addressing, e.g., the excitation spectrum menotti02 , the shape of the trapped gas dunjko01 , and correlation properties gangardt03 .

The many-body ground state has two limiting forms. In the Tonks-Girardeau (TG) regime, where $`\gamma 1`$, the ground state becomes correlated in order to minimize the interaction energy and the bosons become impenetrable, behaving like fermions as described in girardeau . The second and higher-order local correlation functions $`g_i`$ vanish kheruntsyan , meaning that no more than one particle can be found at a given position. On the other hand, in the mean-field (MF) regime when $`\gamma 1`$, the GP equation describes the system well. In this regime the healing length, $`l_h=\mathrm{}/\sqrt{mgn_{\text{1D}}}`$, is much larger than the mean inter-particle distance. Note the counter-intuitive result that the system reaches the correlated regime for low 1D densities, contrary to the 3D case where $`n_{3\text{D}}a_s^3=(\gamma /2\pi )^31`$ corresponds to the correlated regime. (Here $`a_s`$ is the zero-energy 3D scattering length vankempen02 , and $`\gamma `$ is the appropriate energy ratio in 3D.)

To probe correlations we measure three-body recombination rates (proportional to the local third-order correlation function $`g_3`$) of 1D gases produced in a 2D optical lattice. This technique was used in Ref. burt97  to demonstrate that there is a reduction of $`g_3`$ in a 3D BEC by a factor of six compared to a thermal gas. We observe a further reduction of three-body recombination in a 1D gas compared to the 3D BEC situation. Even though $`\gamma 0.5`$ for our system, which is far from the TG regime, this is a signature that the correlations are significant due to the “fermionization” of the particles.

We realize a 1D gas by confining a 3D gas sufficiently tightly in two directions that the radial confinement energy $`\mathrm{}\omega _{}`$ is much larger than all other relevant energies in the system: $`E_{\text{unc}}`$, $`E_{\text{cor}}`$, $`k_BT`$, and the axial trapping energy $`\mathrm{}\omega _z`$. Since $`a_s`$ is much smaller than $`a_{}=\sqrt{\mathrm{}/m\omega _{}}`$ ($`a_s/a_{}0.1`$ in our system), the atom-atom interaction strength is largely determined by $`a_s`$, with only a small correction due to confinement olshanii98 ; eite : $`a_{\text{eff}}=a_s/(11.46a_s/\sqrt{2}a_{})`$. There is no excitation in the radial direction and by integrating over the radial coordinates, one can show olshanii98  that the system is formally equivalent to a true 1D gas with interaction strength $`g=4\mathrm{}^2a_{\text{eff}}/ma_{}^2`$, so that $`\gamma =2a_{\text{eff}}/(n_{\text{1D}}a_{}^2)`$. The 3D density is related to the effective 1D density by $`n_{1\text{D}}=1/d=\pi a_{}^2n_{3\text{D}}`$.

Our approach is to load a BEC into the ground state of a deep 2D optical lattice so that the BEC is divided into an array of independent 1D quantum gases, each tubular lattice site acting as a highly anisotropic trap. Our experimental apparatus has been described elsewhere peil03 . We achieve BEC with up to $`N_0=5\times 10^5`$ atoms in the $`(F,m_F)=(1,1)`$ hyperfine state of 87Rb (for which $`a_s=5.313`$ nm vankempen02 ). A Ioffe-Pritchard trap confines the atoms with initial “tight” trap frequencies of $`\nu _x=\nu _z=210`$ Hz, and $`\nu _y=8.2`$ Hz, giving a peak atomic density of up to $`3\times 10^{14}`$ cm-3. Before applying the optical lattice, we adiabatically lower $`\nu _x`$ and $`\nu _z`$ to a “weak” trap frequency of 28 Hz, resulting in peak densities of $`5\times 10^{13}`$ cm-3.

We create a 2D optical lattice from two independent, retro-reflected 1D lattices which lie in the $`xy`$-plane and intersect at an angle of 80 degrees. The independent 1D lattices are detuned from each another by 5 MHz. All beams derive from a Ti:sapphire laser operating at $`\lambda =810.08`$ nm (detuned below both $`5S5P`$ transitions at 795 nm and 780 nm), and the polarizations of the lattice beams are in the $`xy`$-plane. Each 1D lattice is measured depth  to be 29(1) $`E_R`$ deep (where $`E_R=h^2/2m\lambda ^2`$stddev . At each lattice site the ground state of the radial motion is well approximated by a gaussian wavefunction with $`a_{}=58.5(5)\mathrm{nm}`$ corresponding to an effective $`\omega _{}/2\pi =`$ 33.8(6) kHz. By observing dipole oscillations following a sudden, brief displacement of the trap center, we measure the axial frequency along the tubes to be $`\omega _z/2\pi =`$ 55.9(6) Hz. This frequency results from the combined effect of the magnetic trap and the dipole potential of the lattice beams along the tubes. To load the atoms into the lattice, the laser light is increased over 200 ms with an approximately half gaussian shape (rms width 70 ms), which is adiabatic with respect to all vibrational excitations. We estimate that the interaction-free tunnelling time from one lattice site to the next for a 29 $`E_R`$ lattice is $``$150 ms. Although this is shorter than the time of the experiment (up to 12 s), it corresponds to an energy much smaller than the interaction energy in the tubes, and should not modify the local 1D correlation properties pedri03 .

To measure the reduction of $`g_3`$ due to correlations, we observe the corresponding reduction in the three-body recombination rate coefficient. The local three-body recombination rate (in either 1D or 3D) is proportional to the cube of the local density. For 87Rb, it is known that two-body losses burt97 ; soding99 , including photoassociation at 810 nm carl , are very small K2note . Our model, therefore, includes only one-body and three-body processes so that the total number of atoms $`N`$ decays according to

$$\frac{dN}{dt}=K_1NK_3^{\text{1D}}n_{\text{3D}}^3𝑑V.$$
(1)
We account for atomic redistribution during decay through the evolution of the density profile. Determination of the three-body recombination rate coefficient $`K_3^{\text{1D}}`$ requires an accurate estimate of the density, which we ascertain from a measurement of the number of trapped atoms as a function of time, along with a determination of the size and shape of the atom cloud.

Figure 1 shows the total number of atoms as a function of time $`t`$ in the lattice, obtained by absorption imaging 34 ms after release from the lattice and magnetic trap. We calibrate our absorption measurements by comparing the observed expansion of a released condensate to the known number-dependent expression for the expansion castin96 . (The inferred absorption cross section agrees ($``$10%) with one calculated from the steady state Zeeman sub-level distribution resulting from optical pumping.) The number of atoms in the BEC fluctuates by less than 20% from shot to shot. We automate the experiment to produce a BEC every minute and each data point is typically an average of five measurements. In order to minimize systematic effects due to long term thermal drift of the trap coils, we run current in the magnetic trap after the imaging such that the total time that the magnetic trap is on is the same for each measurement. While the atoms are in the lattice, we apply a radio frequency shield soding99  tuned 500 kHz above the minimum of the trap to reduce heating without significantly increasing the loss of atoms trapped in the lattice.

We measure the size of the lattice-trapped cloud in the $`xy`$-plane by phase contrast imaging. The initial column density distribution of the cloud is well described by an integrated Thomas-Fermi (TF) profile, with radii of $`R_x=13.1(5)`$ $`\mu `$m and $`R_y=22.5(10)`$ $`\mu `$m so that the observed number of atoms per tube at $`(x,y)`$ is well described by $`𝒩_{\text{tube}}=𝒩_{\text{max}}(1(x/R_x)^2(y/R_y)^2)^{3/2}`$, where $`𝒩_{\text{max}}=5N_0\lambda ^2/8\pi R_xR_y`$ is the number of atoms in the central tube and $`\lambda /2`$ is the spacing of the tubes. Based on the initial total number and the measured sizes of the cloud, we determine $`𝒩_{\text{max}}=230(40)`$. In the $`xz`$-plane, we measure the size of the cloud by absorption imaging. The initial $`xz`$ density distribution is also described by a TF profile, of radii $`R_x=15(2)`$ $`\mu `$m (in agreement with our phase contrast $`xy`$-measurement) and $`R_z=17(2)`$ $`\mu `$m. For our parameters, the atom distribution along the tubes (along z) is not expected to deviate significantly from a TF profile dunjko01 ; indeed $`R_z`$ agrees with the 1D TF value calculated based on $`𝒩_{\text{max}}`$. We note that the peak density is $`1\times 10^{15}`$ cm-3, which would lead to rapid three-body loss in a 3D system.

We observe that the cloud slowly expands in the $`z`$-direction over the course of the measurement: the cloud expands by 5 $`\mu `$m in 2 s while the $`R_x`$ and $`R_y`$ radii remain constant. This expansion is consistent with a 1 kHz/s rate of energy increase, which is less than the expected rate due to spontaneous emission with full equilibration between the radial and longitudinal degrees of freedom, but more than the expected rate with no equilibration. Regardless of the mechanism, this expansion reduces the density only modestly during the first two seconds, when most of the three-body decay occurs. We take this into account empirically in modelling the decay.

To model the decay using Eq. 1, we assume an overall 3D TF density profile with gaussian radial distributions within each tube. In addition, for simplicity of modelling we assume $`K_3^{\text{1D}}`$ to be a constant (see below). With these approximations, Eq. 1 becomes:

$$\frac{dN}{dt}=K_1N\alpha (t)K_3^{\text{1D}}N^3,$$
(2)
where $`\alpha (t)=(25/896\pi ^4)(\lambda ^2/a_{}^2R_xR_yR_z(t))^2`$. The radii $`R_x`$ and $`R_y`$ are kept constant at their measured values, and $`R_z(t)`$ grows linearly in time at the measured rate of 2.5 $`\mu `$m/s. This differential equation has an analytic solution which gives the total number as a function of time, to which we fit the data of Fig. 1.

With this analysis, we determine $`K_3^{\text{1D}}=1.2(7)\times 10^{30}`$ cm6 s-1 and $`K_1=0.16(2)`$ s-1. This result is relatively insensitive to the specific model used for atomic spatial redistribution during decay, and variations among realistic models fall within the quoted uncertainties. We attribute $`K_1`$ mainly to optical pumping to states other than the original $`(1,1)`$ state, which are untrapped in the combined optical, magnetic and gravitational potential. For the depth and detuning of our lattice, the majority of photon scattering events return the atoms to the original $`(1,1)`$ Zeeman sub-level and do not contribute to $`K_1`$. The calculated loss rate for a 29 $`E_R`$ lattice is 0.17 s-1, in good agreement with our measured $`K_1`$.

To determine the reduction of three-body recombination in 1D, we must compare $`K_3^{\text{1D}}`$ to $`K_3^{\text{3D}}`$. A comparison in the same apparatus reduces systematic uncertainty due to our 15% number uncertainty. (The effect of this systematic uncertainty is not eliminated entirely because the power law dependence of $`dN/dt`$ on $`N`$ is different in 1D and 3D.) We therefore repeat our experiments in a tight magnetic trap in the absence of a lattice, similar to Refs. burt97 ; soding99 . (But see arescale .) For the (1,-1) state, we measure F2note  $`K_3^{\text{3D}}=8.3(20)\times 10^{30}`$ cm6 s-1 (see Fig. 2), which is in agreement with the value of $`5.8(1.9)\times 10^{30}`$ cm6 s-1 measured in burt97 .

Comparing our measurements, we find that the ratio of the three-body decay coefficients in 1D and 3D is 0.14(9). This represents a factor of seven reduction in $`g_3^{\text{1D }}`$ over $`g_3^{\text{3D }}`$, a clear signature of correlations.

For comparison of the observed reduction in $`g_3`$ with theory, we calculate $`\gamma `$ at the center of each tube. From the experimentally determined density distribution we find at $`t=0`$ that $`\gamma >0.34`$, with 80% of the atoms having $`0.34<\gamma <0.65`$, and median value $`\gamma _\text{m}=0.45`$. We do not expect correlations to vary significantly over this range of $`\gamma `$ kheruntsyan , so the assumption of using a single average $`K_3^{\text{1D}}`$ in the model should be reasonable.

In Fig. 3 we compare the measured reduction in the three-body loss rate coefficient with theoretical estimates for $`g_3`$ in an interacting 1D gas at $`T=0`$. There is currently no theory for $`g_3`$ in the regime intermediate to the MF and TG limits ($`\gamma 1`$). We have plotted an approximate $`g_3`$ in the intermediate regime, using the expression $`g_3(g_2)^3`$ shlyapnikov , and the value of $`g_2`$ for a homogenous system from Ref. kheruntsyan . The value of $`g_2`$ is expected to be insensitive to $`T`$ for $`TT_d𝒩_{\text{max}}\mathrm{}\omega _z/k_B`$ in this range of $`\gamma `$. From the measured values of $`𝒩`$ and $`\omega _z`$ we estimate the distribution of degeneracy temperatures, finding at $`t=0`$ a peak value of $``$13 kHz and a median value of $``$9 kHz. While it is difficult to measure the temperature in our system, the measured size at $`t=0`$ is consistent with zero-temperature TF theory, and is certainly much less than $`m\omega _z^2R_z^2/2`$, which during the first several seconds does not exceed $`6`$ kHz.

The reduction in $`K_3^{\text{1D}}`$ relative to $`K_3^{\text{3D}}`$ indicates that we are beyond the mean-field regime, and is a signature of the “fermionization” of the bosonic particles in 1D. The fundamental effects of low dimensionality on the correlation properties of a quantum Bose gas are also of practical interest, given the interest in the physics of “atom lasers” loaded in waveguides. In addition, these experiments indicate that high density, strongly correlated 1D systems can be realized without fast decay due to three-body recombination.

We acknowledge helpful conversations with M. Olshanii, G. Shlyapnikov and K. Kheruntsyan, and partial support from ARDA, ONR and NASA.


text:
# Dynamics of Rumor Spreading in Complex Networks

## I Introduction

During the last years, many systems have been analyzed from the perspective of graph theory doro ; bara02 . It turns out that seemingly diverse systems such as the Internet, the World Wide Web (WWW), metabolic and protein interaction networks and food webs, to mention a few examples, share many topological properties strogatz . Among these properties, the fact that one can go from one node (or element) of the network to another node passing by just a few others is perhaps the most popular property, known as “six degrees of separation” or small-world (SW) property strogatz ; ws98 . The SW feature has been shown to improve the performance of many dynamical processes as compared to regular lattices; a direct consequence of the existence of key shortcuts that speed up the communication between otherwise distant nodes and of the shorter path length among any two nodes on the net doro ; bara02 ; strogatz .

However, it has also been recognized that there are at least two types of networks fulfilling the SW property but radically different as soon as dynamical processes are ran on top of them. The first type can be called “exponential networks” since the probability of finding a node with connectivity (or degree) $`k`$ different from the average connectivity $`k`$ decays exponentially fast for large $`k`$ ama00 . The second kind of networks comprises those referred to as “scale-free” (SF) networks bar99 . For these networks, the probability that a given node is connected to $`k`$ other nodes follows a power-law of the form $`P(k)k^\gamma `$, with the remarkable feature that $`\gamma 3`$ for most real-world networks doro ; bara02 .

The heterogeneity of the connectivity distribution in scale-free networks greatly impacts the dynamics of processes that they support. One of the most remarkable examples is that an epidemic disease will pervade in an infinite-size SF network regardless of its spreading rate pv01a ; moreno02 ; virusreview ; av03 ; n02b . The change in the behavior of the processes is so radical in this case that it has been claimed that the standard epidemiological framework should be carefully revisited. This might be bad news for epidemiologists, and those fighting natural and computer viruses. On the other hand, in a number of important technological and commercial applications, it is desirable to spread the “epidemic” as fast and as efficient as possible, not to prevent it from spreading. Important examples of such applications are epidemic (or rumor-based) protocols for data dissemination and resource discovery on the Internet vogels ; ep3 ; p2p ; ep1 , and marketing campaigns using rumor-like strategies (viral marketing).

The above applications, and their dynamics, have passed almost unnoticed zanette ; liu  to the physics community working on complex networks despite the fact that they have been extensively studied by computer scientists and sociologists ep1 ; dkbook . The problem here consists of designing an epidemic (or rumor-mongering) algorithm in such a way that the dissemination of data or information from any node of a network reaches the largest possible number of remaining nodes. Note that in this case, in contrast to epidemic modeling, one is free to design the rules of epidemic infection in order to reach the desired result, instead of having to model an existing process. Furthermore, in a number of applications, such as peer-to-peer file sharing systems p2p  built on top of the Internet and grid computing grid , the connectivity distribution of the nodes can also be changed in order to maximize the performance of such protocols.

In this paper we study in detail the dynamics of a generic rumor model dk64  on complex scale-free topologies through analytic and numerical studies, and investigate the impact of the interaction rules on the efficiency and reliability of the rumor process. We first solve the model analytically for the case of exponential networks in the infinite time limit and then introduce a stochastic approach to deal with the numerical solution of the mean-field rate equations characterizing the system’s dynamics. The method pre98 ; jgr ; mgp03  is used to obtain accurate results for several quantities when the topology of random SF networks is taken into account, without using large and expensive Monte Carlo (MC) simulations. The rest of the paper is organized as follows. Section II is devoted to introducing the rumor model and to derive the mean-field rate equations used throughout the paper. In Section III we deal with the stochastic approach, and compare its performance with analytical and MC calculations in homogeneous systems. We extend the method to the case of power-law distributed networks and present the results obtained for this kind of networks in Section IV and V. Finally, the paper is rounded off in the last Section, where conclusions are given.

## II Rumor Model in Homogeneous Networks

The rumor model is defined as follows. Each of the $`N`$ elements of the network can be in three different states. Following the original terminology and the epidemiological literature dkbook , these three classes correspond to ignorant, spreader and stifler nodes. Ignorants are those individuals who have not heard the rumor and hence they are susceptible to be informed. The second class comprises active individuals that are spreading the rumor. Finally, stiflers are those who know the rumor but that are no longer spreading it. The spreading process evolves by directed contacts of the spreaders with others in the population. When a spreader meets an ignorant the last one turns into a new spreader with probability $`\lambda `$. The decay of the spreading process may be due to a mechanism of “forgetting” or because spreaders learn that the rumor has lost its “news value”. We assume this latter hypothesis as the most plausible so that the contacting spreaders become stiflers with probability $`\alpha `$ if they encounter another spreader or a stifler. Note that as we are designing our rumor strategy in such a way that the fraction of the population which ultimately learns the rumor be the maximum possible, we have assumed that contacts of the type spreader-spreader are directed, that is, only the contacting individual loses the interest in propagating the rumor further. Therefore, there is no double transition to the stifler class.

In a homogeneous system, the original rumor model due to Daley and Kendall dk64  can be described in terms of the densities of ignorants, spreaders, and stiflers, $`i(t)`$, $`s(t)`$, and $`r(t)`$, respectively, as a function of time. Besides, we have the normalization condition,

$$i(t)+s(t)+r(t)=1.$$
(1)
In order to obtain an analytical insight and a way to later test our numerical approach, we first study the rumor model on top of exponentially distributed networks. These include models of random graphs as well as the Watts and Strogatz (WS) small-world model ws98 ; strogatz . This model produces a network made up of $`N`$ nodes with at least $`m`$ links to other nodes. The resulting connectivity distribution in the random graph limit of the model ws98  takes the form

$$P(k)=\frac{m^{km}}{(km)!}e^m,$$

which gives an average connectivity $`k=2m`$. Hence, the probability that a node has a degree $`kk`$ decays exponentially fast and the network can be regarded as homogeneous.

The mean-field rate equations for the evolution of the three densities satisfy the following set of coupled differential equations:

$`{\displaystyle \frac{di(t)}{dt}}`$ $`=`$ $`\lambda ki(t)s(t),`$ (2)
$`{\displaystyle \frac{ds(t)}{dt}}`$ $`=`$ $`\lambda ki(t)s(t)\alpha ks(t)[s(t)+r(t)],`$ (3)
$`{\displaystyle \frac{dr(t)}{dt}}`$ $`=`$ $`\alpha ks(t)[s(t)+r(t)],`$ (4)
with the initial conditions $`i(0)=(N1)/N`$, $`s(0)=1/N`$ and $`r(0)=0`$. The above equations state that the density of spreaders increases at a rate proportional to the spreading rate $`\lambda `$, the average number of contacts of each individual $`k`$ and to the densities of ignorant and spreader individuals, $`i(t)`$ and $`s(t)`$, respectively. On the other hand, the annihilation mechanism considers that spreaders decay into the stifler class at a rate $`\alpha k`$ times the density of spreaders and of non-ignorant individuals $`1i(t)=s(t)+r(t)`$.

The system of differential equations (2-4) can be analytically solved in the infinite time limit when $`s(\mathrm{})=0`$. Using equation (1), we have that $`_0^{\mathrm{}}s(t)𝑑t=r_{\mathrm{}}=lim_t\mathrm{}r(t)`$. Introducing the new variable $`\beta =1+\lambda /\alpha `$ we obtain the transcendental equation,

$$r_{\mathrm{}}=1e^{\beta r_{\mathrm{}}}.$$
(5)
Equation (5) always admits the trivial solution $`r_{\mathrm{}}=0`$, but at the same time it also has another physically relevant solution for all values of the parameters $`\lambda `$ and $`\alpha `$. This can be easily appreciated since the condition,

$$\frac{d}{dr_{\mathrm{}}}\left(1e^{\beta r_{\mathrm{}}}\right)|_{r_{\mathrm{}}=0}>1,$$
(6)
reduces to $`\lambda /\alpha >0`$. That is, there is no “rumor threshold” contrary to the case of epidemic spreading moreno02 . This strikingly different behavior does not come from any difference in the growth mechanism of $`s(t)`$ $``$the two are actually the same$``$, but from the disparate rules for the decay of the spreading process.

On the other hand, this result also points out that a mathematical model for the spreading of rumors can be constructed in many different ways. The results of this paper, however, indicate that the presence of spreader annihilation terms due to spreader-spreader and spreader-stifler interactions is very relevant for practical implementations mnv03 ; nm03 . We shall come back to this point later on.

## III Stochastic Numerical Approach

Recently mgp03 , we have introduced a numerical technique pre98  to deal with the mean-field rate equations appearing in epidemic-like models. It solves the differential equations by calculating the passage probabilities for the different transitions. The main advantage of this method, as compared to MC simulations, is its modest memory and CPU time requirements for large system sizes. Besides, we do not have to generate any network. Instead, we produce a sequence of integers distributed according to the desired connectivity distribution $`P(k)`$. The numerical procedure here proceeds as follows. At each time step until the end of the rumor spreading process, the following steps are performed:

1. Identify from the mean-field rate equations the transition probabilities per time unit from one state into the following one, that is, from the $`i`$ class to the $`s`$ class, $`W_{is}`$, and finally to the $`r`$ class, $`W_{sr}`$.
2. Calculate the mean time interval, $`\tau `$, for one transition to occur. This is determined as the inverse of the sum of all the transition probabilities; $`\tau =1/(W_{is}+W_{sr})`$.
3. Stochastically decide what transition will actually take place. This is done by deciding that the probabilities for both transitions are given by $`\mathrm{\Pi }_{is}=W_{is}\tau `$ and $`\mathrm{\Pi }_{sr}=W_{sr}\tau `$, respectively, materializing the choice by generating a random number between 0 and 1.

The numerical algorithm described above does not depend on the topological features of the network on top of which the rumor dynamics is taking place. Indeed, all the topological information, including correlations, enters in the computation of the transition probabilities. We should note here that the present results are obtained for uncorrelated networks. The method could also be applied to correlated networks without explicit generation of them. In that case, one should work with the two point correlation function $`P(k,k^{})`$ mgp03  instead of using $`P(k)`$. On the other hand, a correlated network could be built up as in vw .

In order to gain confidence with the method and to show its soundness, we show in Table 1 the values of $`r_{\mathrm{}}`$ obtained from Eq. (5), MC simulations and the stochastic approach for homogeneous networks. In this case, the transition probabilities are the same for all the elements within a given class ($`i`$, $`s`$ or $`r`$) irrespective of their actual connectivities. From equations (2-4) we get

$`W_{is}(t)`$ $`=`$ $`N\lambda ki(t)s(t),`$ (7)
$`W_{sr}(t)`$ $`=`$ $`N\alpha ks(t)[s(t)+r(t)],`$ (8)
for the transitions from the ignorant to the spreader class and from the spreader to the stifler class, respectively.

It can be seen from Table 1 that the difference between the SNA result and the MC simulations is less that $`1.4\%`$, indicating the reliability of the SNA approach. The remaining small differences between the SNA and the MC results is mainly due to the fact that the homogeneous SNA model does not take into account the exponentially decaying fluctations in the connectivity of WS networks. On the other hand, MC simulations of the rumor dynamics for a network made up of $`N=10^4`$ nodes, averaged over at least 10 different network realizations and 1000 iterations, took several hours. Eventually, this method takes up to a few days when increasing the system size and decreasing the value of $`\alpha `$. On the contrary, the stochastic approach is very fast. Indeed, for the same parameter values, the numerical simulation takes around 5 minutes CPU time in a 2.0Ghz-P4 PC. Therefore, having such a method will allow us to scrutinize very efficiently and accurately the whole phase diagram and time profiles of the process under study. In what follows, we analyze in detail the dynamics of the rumor spreading process by numerically solving the mean-field rate equations for SF networks.

## IV power-law distributed networks

The heterogeneity of the connectivity distribution inherent to SF networks significantly affects the dynamical evolution of processes that take place on top of these networks pv01a ; moreno02 ; virusreview ; av03 ; n02b ; havlin01 ; newman00 ; bar00 ; mv02 . We have learned in recent years that the fluctuations of the connectivity distribution, $`k^2`$, can not be neglected even for finite size systems virusreview . Thus, the system of differential equations (2-4) should be modified accordingly. In particular, we should take into account that nodes could not only be in three different states, but also they belong to different connectivity classes $`k`$. Let us denote by $`i_k(t)`$, $`s_k(t)`$, and $`r_k(t)`$ the densities of ignorants, spreaders and stiflers with connectivity $`k`$, respectively. In addition, we have that $`i_k(t)+s_k(t)+r_k(t)=1`$. The mean-field rate equations now read as,

$`{\displaystyle \frac{di_k(t)}{dt}}`$ $`=`$ $`\lambda ki_k(t){\displaystyle \underset{k^{}}{}}{\displaystyle \frac{k^{}P(k^{})s_k^{}(t)}{k}},`$ (9)
$`{\displaystyle \frac{ds_k(t)}{dt}}`$ $`=`$ $`\lambda ki_k(t){\displaystyle \underset{k^{}}{}}{\displaystyle \frac{k^{}P(k^{})s_k^{}(t)}{k}}`$ (10)
$`\alpha ks_k(t){\displaystyle \underset{k^{}}{}}{\displaystyle \frac{k^{}P(k^{})[s_k^{}(t)+r_k^{}(t)]}{k}},`$
$`{\displaystyle \frac{dr_k(t)}{dt}}`$ $`=`$ $`\alpha ks_k(t){\displaystyle \underset{k^{}}{}}{\displaystyle \frac{k^{}P(k^{})[s_k^{}(t)+r_k^{}(t)]}{k}},`$ (11)
where $`P(k)`$ is the connectivity distribution of the nodes and $`_k^{}k^{}P(k^{})s_k^{}(t)/k`$ is the probability that any given node points to a spreader. We start from a randomly selected spreader and all the remaining nodes in the ignorant class. The summation in Eq. (10) stands for the probability that a node points to a spreader or a stifler. Note that as before, we do not allow for double transitions from the spreader to the stifler class. Next, we compute the respective transition probabilities. In this case, we should also consider that transitions from one state into another also take place within connectivity classes. Thus, the transition probabilities depend on $`k`$ as well. From Eq. (10) we obtain,

$`W_{is}(t,k)`$ $`=`$ $`\lambda kNP(k)i_k(t){\displaystyle \underset{k^{}}{}}{\displaystyle \frac{k^{}P(k^{})s_k^{}(t)}{k}},`$ (12)
$`W_{sr}(t,k)`$ $`=`$ $`\alpha kNP(k)s_k(t){\displaystyle \underset{k^{}}{}}{\displaystyle \frac{k^{}P(k^{})[s_k^{}(t)+r_k^{}(t)]}{k}},`$
where all the topological information is contained. Finally for the mean time interval after $`i1`$ transitions, $`\tau `$, we find at each time step

$$\tau =\frac{1}{W_{is}(t)+W_{sr}(t)},$$
(14)
with $`W_{is}(t)=_kW_{is}(t,k)`$, $`W_{sr}(t)=_kW_{sr}(t,k)`$ and $`t=_j^{i1}\tau _j`$, where the $`\tau _j`$s are the mean times of the $`i1`$ previous transitions. At this point, the identification of what transition takes place and which connectivity class is affected proceeds as defined in step $`3`$ of the previous section.

## V results and discussion

The stochastic method described above can be used to explore several quantities characterizing the dynamics of the rumor spreading process. Throughout the rest of the paper we set $`\lambda =1`$ without loss of generality and vary the value of $`\alpha `$. We first generated a sequence of integers distributed according to $`P(k)k^\gamma `$ with $`\gamma =3`$ and $`k=6`$. As initial condition we use $`r_k(t=0)=0`$, and

$$s_k(t)=\{\begin{array}{cc}\frac{1}{NP(k)}\hfill & k=k_i\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(15)
where $`k_i`$ is the connectivity of the randomly chosen initial spreader. The results are then averaged over at least $`1000`$ different choices of $`k_i`$.

One of the most important practical aspects of any rumor mongering process is whether or not it reaches a high number of individuals. This magnitude is simply given by the final density of stiflers and is called reliability of the rumor process. However, it is also of great importance for potential applications that higher levels of reliability are reached as fast as possible, which constitute a practical measure of the cost associated to such levels of stiflers. For example, in technological applications, where one may consider several strategies mnv03 ; nm03 , it is possible to define a key global quantity, the efficiency of the process, which is the ratio between the reliability and the traffic imposed to the network. For these applications it is not only important to have high levels of reliability but also to achieve these with the lowest possible load resulting from the epidemic protocol’s message passing traffic. This is important in order to avoid network congestion and also to reduce the amount of processing power used by nodes participating in the rumor process.

In order to analyze, from a global perspective, this trade-off between reliability and cost, we use time as a practical measure of efficiency. We call a rumor process less efficient than another if it needs more time to reach the same level of reliability. Figure 1 shows the time evolution of the density of stiflers for several values of the parameter $`\alpha `$. It turns out, as expected, that the number of individuals who finally learned the rumor increases as the probability of becoming stifler decreases. On the other hand, the time it takes for $`R(t)`$ to reach its asymptotic value slightly increases with $`\alpha ^1`$, but clear differences do not arise for the two extreme values of $`\alpha `$. In fact, for a given time after the beginning of the rumor propagation, the density of stiflers scales with the inverse of $`\alpha `$. This behavior is further corroborated in the inset, where the growth of the density of spreaders as time goes on is shown for the same values of the parameter $`\alpha `$. While the peaks of the curves get larger and larger, the times at which the maxima are reached are of the same order of magnitude and thus the meantimes of the spreading processes do not differ significantly.

Figure 2 shows another aspect worth taking into account when dealing with rumor algorithms. For a given level of reliability, it is also of interest to know the distribution of ignorants (or stiflers) by classes $`k`$. The figure shows a coarse-grained picture of Fig. 1, where the density of ignorants $`i_k`$ according to the connectivity of the individuals has been represented for different values of $`\alpha `$. The results indicate that the probability of having an ignorant with a connectivity $`k`$, at the end of the rumor propagation, decays exponentially fast with a sharp cut-off $`k_c`$ for large connectivity values, which depends on $`\alpha `$. In fact, $`k_c`$ is always well below the natural cut-off of the connectivity distribution ($`10^2`$) even for small values of $`\alpha `$. This implies that hubs effectively learn the rumor.

We can further scrutinize the dynamics of the rumor spreading process by looking at the final density of stiflers when the initial spreader has a given connectivity $`k_i`$. Figure 3 represents the reliability as a function of time (in units of $`\alpha ^1`$) when the rumor starts propagating from a node of connectivity $`k_i=k_{min}=3`$, $`k_i=k=6`$, $`k_i=20`$ and $`k_i=k_{max}280`$ for two different values of $`\alpha `$: $`0.1`$ (main figure) and $`1.0`$ (inset). Interestingly, the final value of $`R(t)`$ does not depend on the initial seed, but reaches the same level irrespective of the connectivity of the very first spreader $`k_i`$. This is a genuine behavior of the rumor dynamics and is the opposite to what has been observed in other epidemic models like the SIR model moreno02 , where the final number of recovered individuals strongly depends on the connectivity of the initially infected individuals. However, a closer look to the spreading dynamics tell us that not all is the same for different initial spreaders.

The figure also indicates that as the connectivity of the seed is increased, the time it takes for the rumor to reach the asymptotic value decreases, so that for a fixed time length the number of individuals in the stifler class is higher when $`k_i`$ gets larger. This feature suggests an interesting alternative for practical applications: start propagating the rumor from the most connected nodes. Even in case that no direct link exits between a node that is willing to spread an update and a hub, a dynamical (or temporal) shortcut to a well connected node could be created in order to speed up the process. With this procedure, the density of stiflers at the intermediate stages of the spreading process could be as much different as $`30\%`$ for moderate values of $`\alpha `$. This translates in less costs, because one can always implement an algorithm that kill off the actual spreading when a given level of reliability is reached. Note, however, that this behavior slightly depends on $`\alpha `$, being the differences always appreciable, but more important as $`\alpha `$ increases.

Finally, we have exploited the fastness of the stochastic approach used here to explore the consequences of implementing three different annihilation rules for the rumor spreading decay. In particular, we consider that the spreading process dies out proportionally only to the number of spreaders ($`ss`$ interactions) or to the number of stiflers ($`sr`$ interactions). This modifies the terms entering in the sum of Eqs. (10-11) so that now the transition probabilities from the $`s`$ into the $`r`$ class read

$`W_{sr}^{ss}(t,k)`$ $`=`$ $`\alpha kNP(k)s_k(t){\displaystyle \underset{k^{}}{}}{\displaystyle \frac{k^{}P(k^{})s_k^{}(t)}{k}},`$ (16)
$`W_{sr}^{sr}(t,k)`$ $`=`$ $`\alpha kNP(k)s_k(t){\displaystyle \underset{k^{}}{}}{\displaystyle \frac{k^{}P(k^{})r_k^{}(t)}{k}},`$ (17)
respectively. Table 2 summarizes the reliability of the process as a function of $`\alpha `$ for the three mechanisms considered note1 . The results indicate that in all variants, the final density of stifler individuals is higher than for the “classical” setting. However, in order to evaluate the efficiency of the process from a global perspective, we must look at the time evolution of the densities as we did before.

In Figs. 4 and 5 we have represented the time (in units of $`\alpha ^1`$) profiles of $`R(t)`$ and $`S(t)`$ for each decay term and several values of $`\alpha `$. From the figures, it is clear that while the final density of stiflers increases when modifying the original decay rules, the time needed to reach such high levels of reliability also increases. This is due to the fact that the tails of the densities of spreaders decay more slowly than before. In particular, it is noticeable that when only spreader-spreader interactions are taken into account in the decay mechanism, the lifetime of the propagation process is more than two times longer than for the other two settings. This means that this implementation is not very suitable for practical applications as the costs associated to the process rise as well. On the other hand, the performance of the spreader-stifler setting seems to depend on the value of $`\alpha `$ in such a way that it is more efficient at both the reliability level and time consumption for a large $`\alpha `$, but not in the middle region of the parameter space. In summary, the present results support that the original model works quite well under any condition, while other variants can be considered depending on the value of $`\alpha `$ used and the type of applications they are designed for.

## VI conclusions

In this paper, we have analyzed the spreading dynamics of rumor models in complex heterogeneous networks. We have first introduced a useful stochastic method that allows us to obtain meaningful time profiles for the quantities characterizing the propagation process. The method is based on the numerical solution of the mean-field rate equations describing the model, and contrary to Monte Carlo simulations, there is no need of generating explicitly the network. This allows to save memory and a fast exploration of the whole evolution diagram of the process.

The kind of processes studied here are of great practical importance since epidemic data dissemination might become the standard practice in multiple technological applications. The results show that there is a fragile balance between different levels of reliability and the costs (in terms of time) associated to them. In this sense, our study may open new paths in the use of rumor-mongering process for replicated database maintenance, reliable group communication and peer to peer networks dave ; deering ; vogels ; ep3 ; p2p ; ep1 . Besides, as shown here, the behavior and features of the different algorithms one may implement are not trivial and depend on the type of mechanisms used for both the creation and the annihilation terms. It is worth noting here that we have studied the simplest possible set of rumor algorithms, but other ingredients such as memory must be incorporated in more elaborated models mnv03 ; nm03 .

Of further interest would be a more careful exploration of the possibility of using dynamical shortcuts for a more efficient spreading of the updates. Our results suggest that it would be more economic to start from hubs and then kill off the updating process when a given level of reliability is reached than starting at random and letting the process dies out by itself. Preliminary studies of more elaborated models aimed at implementing a practical protocol confirm our results nm03 . This feature is specially relevant for the understanding and modeling of social phenomena such as the spreading of new ideas or the design of efficient marketing campaigns.

###### Acknowledgements.

We would like to thank A. Vespignani for many useful discussions and comments. Y. M. thanks M. Vázquez-Prada for helpful comments and the hospitality of BT Exact, U.K. where parts of this work were carried out. Y. M. acknowledges financial support from the Secretaría de Estado de Educación y Universidades (Spain, SB2000-0357). This work has been partially supported by the Spanish DGICYT project BFM2002-01798.


text:
# Adaptation and enslavement in endosymbiont-host associations

## I  Introduction

Cooperation between endosymbionts and their hosts is most striking in the common case where the symbiont can survive in a free-living form and each host is infected anew in its infancy Maynard-Smith and Szathmary , Frank , Herre et al. . If the host is unable to respond to individual symbionts who do not contribute to the association, then free-riders will prosper and the symbiont will ultimately become parasitic on its host Bremermann and Pickering , Knolle . On the other hand, if the host is able to punish defecting symbionts, the expectation is that a “you scratch my back, and I’ll scratch yours” relationship will develop Trivers , Axelrod . It is supposed that the reciprocation of both cooperation and defection maintains a mutualistic relationship, from which both parties benefit. This appears to be supported by game-theoretic analyses of simple models of cooperation, with tit-for-tat strategies being successful in Prisoner’s Dilemma scenarios Axelrod and Hamilton , Nowak and Sigmund , Dugatkin . Although tit-for-tat emerges as the winning strategy under a wide range of conditions, most studies of cooperation are restricted to intra-specific competition. It is not clear that this can be generalized to the host-endosymbiont system, where the Prisoner’s Dilemma occurs between two species.

In natural systems, it is found that many horizontally transmitted non-obligate endosymbionts apparently generate large benefits for their hosts. Examples include nitrogen fixation by rhizobial bacteria in legumes Sprent et al. ; enhanced uptake of nutrients due to mycorrhizal fungi Smith and Read ; carbon uptake by green algae in many aquatic invertebrates, such as the zooxanthellae of corals Veron ; and bioluminescence provided by bacteria in squid and fish McFall-Ngai . Surprisingly, in very few cases have endosymbionts been shown to benefit significantly from their interactions with host organisms Douglas and Smith . For example, *Rhizobia* reproduce happily enough when free-living in the soil, but most of these bacteria hardly reproduce at all once inside the root nodules of legumes Sprent et al. , Denison . For the putative benefits of symbiotic life as zooxanthellae, dinoflagellates give up their cell wall and their flagella, sacrifice most of their photosynthetic products, and reduce their reproductive rate Saffo . The non-obligate forms of endomycorrhizal fungi can be readily grown alone, and there is no evidence for a fitness-related benefit from symbiosis for them Smith and Read . In particular, there is no evidence that carbon flows from the plant to the fungus in return for the nutrients these endosymbionts provide. Thus notions of host and parasite appear to be turned on their head, and this asymmetry is not explained by the standard tit-for-tat solution to the problem of cooperation.

Most previous models of cooperation have only represented one population. In this paper we present a simple model of a two-species association based on reciprocal altruism, and explore whether it can explain both the maintenance of symbiotic relationships and the apparent enslavement of the symbionts by their hosts. In keeping with much previous work on the evolution of cooperation, we formalize the Prisoner’s Dilemma in terms of the costs, $`c`$, and benefits, $`b`$, of cooperative acts. The species have variable degrees of cooperation, depending on the value of these parameters Frean , Wahl and Nowak \[1999a, b\], Killingback and Doebli . The degree of cooperation adopted by each species is subject to evolutionary change. The challenge is then to demonstrate that a mutualistic relationship is stable in the face of mutations in the parameters governing cooperation.

A core asymmetry between the host and the endosymbionts is that the symbionts have a shorter generation time. This rapid turnover allows their population to respond quickly to changes in the hosts’ behavior. A consequence of this differential adaptive rate is that the tit-for-tat strategy is no longer optimal. Rather the endosymbiont becomes fully cooperative and the host becomes miserly, adopting a level of cooperation just high enough to maintain the symbiont cooperation. The differential rates of adaptation appear to lead directly to the differential pay-offs inferred from studies of host-endosymbiont systems.

## II  Misers and Slaves

In this section we consider just a pair of agents that interact, to illustrate the idea that a slower rate of adaptation can be advantageous. This appears closely related to the Red King effect Bergstrom and Lachmann , although that model specifically avoids treating the Prisoner’s Dilemma scenario, in which there is a constant temptation to defect on cooperative agreements. In later sections we apply this notion to a model of endosymbiosis.

Consider two agents $`S`$ and $`H`$, later we will interpret these as symbiont and host. Suppose that $`S`$ adopts some behaviour to a degree $`s`$ at a cost to itself of $`c_s`$ per unit, and a benefit to $`H`$ of $`b_{sh}`$ per unit. For example, this could be nitrogen fixation by *Rhizobia* which makes nitrogen available to legumes in a useful form, but which entails a metabolic cost to the *Rhizobia*. Suppose that $`H`$ can also adopt some behaviour $`h`$ between zero and one. This benefits $`S`$ by $`b_{hs}`$ at a cost of $`c_h`$ per unit. The net payoff to $`S`$ from the association is then

$$w_s=b_{hs}hc_ss$$
(1)
while that to $`H`$ is

$$w_h=b_{sh}sc_hh$$
(2)
For the purposes of our model we take $`s`$ and $`h`$ to be bounded below by zero and above by one. That is, we assume there is some natural limit to the amount of help that one organism can provide another. If benefits exceed costs for both parties, the maximum total payoff, $`w_s+w_h`$, is achieved when both $`s`$ and $`h`$ are at their maximum values. However any non-zero $`s`$ or $`h`$ only benefits the other party, and so is an altruistic act of cooperation. By reducing the level of cooperation, an agent reduces the costs, and so there is an incentive for each agent to defect. This selfish behavior results in the payoff to each agent decreasing and we have what is known in game theory terms as a Prisoner’s Dilemma Axelrod .

To provide an incentive for $`S`$ to cooperate, $`H`$’s behaviour needs to depend on $`s`$ in some way, and the simplest non-trivial contingency is linear, $`h=H_0+(H_1H_0)s`$. Here $`H_0`$ is the response to complete defection ($`s=0`$) and $`H_1`$ is the response to complete cooperation ($`s=1`$). Parameters $`H_0`$ and $`H_1`$ lie between zero and one, ensuring that $`h`$ also lies in this range. The difference $`H_1H_0`$ can be thought of as the ‘responsiveness’ of $`H`$. Completely unresponsive agents have $`H_0=H_1`$ , extreme cases being the naive cooperator Always Cooperate (1,1) and the stalwart defector Always Defect (0,0). For unresponsive agents, $`h`$ is independent of $`s`$ and the gradient $`w_s/s`$ is $`c_s`$ , so $`S`$ is motivated to decrease its level of cooperation. In contrast, the Tit-for-Tat strategy ($`H_0`$, $`H_1`$) = (0,1) is as responsive as possible (see fig. 1). In this case we have $`h=s`$ and $`w_s/s=b_{hs}c_s>0`$, so $`S`$ should cooperate more to increase its payoff. Notice that this argument does not depend on specifying the details of $`S`$’s strategy (meaning the way it arrives at the value $`s`$, given $`h`$) - it merely indicates whether $`S`$ should increase or decrease the level of cooperation it adopts with $`H`$ in pursuit of greater payoffs.

More generally, substituting the expression for $`h`$ into the payoff for $`S`$ and differentiating with respect to $`s`$, shows that $`S`$ is motivated to cooperate as much as it can, provided that the following ‘incentive condition’ is met,

$$H_1H_0>\frac{c_s}{b_{hs}}.$$
(3)
This defines two regions on a plot of $`H_1`$ versus $`H_0`$ (see fig. 2). If $`H`$’s responsiveness is too low, $`S`$ would benefit from lowering its level of cooperation, while if $`H`$’s responsiveness is above the threshold $`c_s/b_{hs}`$, $`S`$ should become as cooperative as it can. Now consider a scenario in which $`S`$ adapts its strategy over time in pursuit of higher payoffs. If the incentive condition (eq. 3) is not met, $`s`$ will decrease to zero. In particular, any flat response ($`H_0=H_1`$) by $`H`$ will lead to pure defection by $`S`$, and consequently cannot account for the persistence of a mutualistic association between them. Conversely, if the incentive condition is met, $`s`$ will rise to one.

Now consider the effect of this on $`H`$. The payoff to an $`H`$ agent whose responsiveness is above the threshold will be $`b_{sh}c_hH_1`$ (since $`s=1`$), compared to $`c_hH_0`$ for an $`H`$ below the threshold (since then $`s=0`$). Since benefits outweigh costs $`b_{sh}c_hH_1>0`$ and so all $`H`$ strategies above the threshold out-compete those below it. The best possible strategy for $`H`$ will be the one with lowest $`H_1`$ that still obeys the incentive condition. This strategy has $`H_0=0`$ and $`H_1`$ just above $`c_s/b_{hs}`$ (it needs to be just above to ensure a strictly positive gradient for $`s`$). Equivalently, the optimal strategy for $`H`$ is to cooperate with $`S`$ at a level

$$h_{\mathrm{miser}}=\left(\frac{c_s}{b_{hs}}s\right)^+$$
(4)
where the superscript ($`+`$) means ‘just above’. The conclusion is that the best strategy to adopt when faced with an optimally adaptive co-player is that of a ‘miser’. The miser reacts with just enough cooperation to encourage the other agent to cooperate more. With $`H_1`$ at $`c_s/b_{hs}`$ the miser collects a payoff of $`b_{sh}c_sc_h/b_{hs}`$, which is not much less than the maximum payoff of $`b_{sh}`$ obtainable by pure defection against the sucker strategy Always Cooperate. By contrast, the highly adaptive agent is induced to cooperate at the maximal level and for this receives a payoff only marginally above zero. If $`H`$ is able to adapt but does so at a significantly lower rate than $`S`$, we might expect a miser-slave relationship to arise and be stable. This phenomenon has been noted in simulations of the prisoner’s dilemma, although in that case the model involved was not amenable to analysis Doebli and Knowlton .

It might be argued at this point that the preceding model is biased, in that the $`H`$ strategies are responsive whereas $`S`$ just assumes a fixed level of cooperation. The $`S`$’s are not given access to the same repertoire of strategies as the $`H`$ strategies, and so perhaps their enslavement is due merely to this built-in deficiency. To show that it occurs in the same way when *both* players may change their strategies, introduce the parameters $`S_0`$ and $`S_1`$. The level of cooperation is $`s=S_0+(S_1S_0)h`$ and we have a simple linear dynamical system in which $`h`$ depends on $`s`$, which depends on $`h`$ itself. The dynamics have a single attractor that is approached rapidly, regardless of whether $`h`$ and $`s`$ are updated synchronously or asynchronously Frean . This point is

$`h^{}=`$ $`\alpha [H_0+S_0(H_1H_0)],`$ (5)
$`s^{}=`$ $`\alpha [S_0+H_0(S_1S_0)],`$ (6)
where $`\alpha =[1(H_1H_0)(S_1S_0)]^1`$. Geometrically, the point ($`h^{}`$, $`s^{}`$) is simply the intersection of the two strategies plotted on the unit square. On fig. 1, an $`S`$ strategy maps $`h`$ back to $`s`$ and must be a line from the top of the figure to its base, crossing the $`H`$ strategy at one point. Successive updates move the joint state towards this intersection point and convergence to equilibrium is rapid (e.g. of the order of 10 to 20 responses). The final levels determine the payoff that agents receive, namely

$`w_h^{}`$ $`=`$ $`b_{hs}s^{}c_hh^{},`$ (7)
$`w_s^{}`$ $`=`$ $`b_{sh}h^{}c_ss^{},`$ (8)
and similarly for $`S`$. If the player $`S`$ evolves its strategy sufficiently rapidly, then it is straightforward to show that the incentive condition is unchanged by this extension.

## III Hosts and Symbionts

So far we have discussed just two agents and assumed that both adapt their behaviour, with one of them ($`S`$) adapting much more rapidly. It is interesting to speculate whether this phenomenon might occur between two co-evolving populations, and specifically between hosts and their endosymbionts. In an endosymbiotic relationship, the lifespan of the host is longer than that of the symbiont. Each of the partners generates payoffs via interactions with the other type, yet uses these to competitively replace only their own type. A single host typically contains many endosymbionts, and these replicate within the host for many generations. In the remainder of this paper we use computer simulations to explore the idea that the preceding analysis of agents $`H`$ and $`S`$ applies similarly to hosts and endosymbionts. The intent is not to present a detailed model of a specific host-endosymbiont system, but to examine in general terms whether this mechanism can be applicable.

Our model (shown schematically in fig. 3) is intended to capture the simplest aspects of endosymbiotic associations. We assume purely horizontal transmission of symbionts, since this is the situation that is most puzzling. Each host organism acquires endosymbionts from the environment at some point in its early life, and a symbiotic association is formed anew with each generation. The host harbors its symbionts for some period, during which the symbionts reproduce with mutations. We assume that symbionts compete for limited resources within each host, leading to traits associated with greater symbiont payoffs becoming more prevalent over time. Hosts are also in a competitive environment, albeit one that occurs over a longer time-scale. For simplicity we assume that variations in fitness are due solely to the host-endosymbiont association.

Specifically, we model a population of 100 individuals of a host species denoted $`H`$, each with parameters ($`H_0`$, $`H_1`$). Each host has an exclusive association with its own population of 25 symbionts, denoted $`S`$, and each such symbiont has a strategy ($`S_0`$, $`S_1`$). In our simulations the population sizes are fixed: introducing fluctuating populations would require additional assumptions (about carrying capacities for instance).

All the evolvable parameters ($`H_0`$, $`H_1`$ for hosts and $`S_0`$, $`S_1`$ for symbionts) were started with random values between zero and one. Each host applied its strategy of contingent cooperation in interacting with each of its symbionts on an individual basis, resulting in payoffs as described above. For simplicity, the simulations shown here used the same costs and benefits for both types of creature ($`c_s=c_h`$ and $`b_{hs}=b_{sh}`$). A very simple model of competitive replacement was used to simulate the effect of natural selection, as follows. For symbionts, two individuals were chosen at random from a single host. The payoff to each individual was calculated using eqs. (8), but with a small amount of random noise (uniformly distributed in the range $`\pm 0.01`$) being added. The symbiont with lower fitness (payoff) was then deleted and replaced by the other. Rather than perfect replication of the better strategy, mutations were introduced by adding random noise uniform in the range $`\pm 0.01`$ to each parameter as it was copied, while constraining them to remain within zero and one. This allowed new strategies to arise that were not present in the original population. Exactly the same procedure was applied to competition between hosts. We refer to $`N`$ replacements, where $`N`$ is the size of the relevant population, as one generation. Endosymbionts underwent 25 such generations per host replacement, and the simulation was run for 100 host generations. The fitness of a given host was taken to be the sum of the payoffs from its interactions with all its endosymbionts at the end of 100 endosymbiont generations. Each new host began life with randomly generated new endosymbionts. This reflects the fact that the most difficult case to understand is horizontal transmission, where symbionts exist in a free-living form and infect hosts anew with each generation. Horizontal transmission also appears to be the most common mode of maintaining symbiont-host associations Maynard-Smith and Szathmary .

Within each host, the symbionts evolve towards either full cooperation or defection, depending on whether the host strategy is in region I or II of fig. 2. The examples in fig. 4 show the evolution of the symbionts within two hosts, one slightly more cooperative than the miser strategy, and one slightly less cooperative. There is initially a rapid change in the symbiont strategies, towards $`S_0=1`$ and $`S_0=0`$, respectively. The change in $`S_1`$ is slower, because in the example given $`H_0`$ is small. As expected, the payoff to the symbionts within a host increases with time.

At the end of the full simulation the miser strategy dominated the host population and most endosymbionts adopted strategies very close to unconditional cooperation (see fig. 5). After an initial transient, lasting for around 20 host generations, the average payoff to the endosymbionts remained at less than one tenth of the averaged payoff to the hosts. Note that the payoff to the symbionts decreases as the hosts evolve. Occasionally a mutation would take a host across the boundary into the region which is too non-cooperative, this host’s symbionts would then stop cooperating, the payoff to the host would fall, and it would be eliminated from the population through the selection process. The asymmetry in the rewards of the association remains even at low cost-to-benefit ratios (fig. 6). As $`c`$ approaches zero there is little cost to cooperation and hence little incentive to defect, and yet the degree of miserliness becomes increasingly severe. Similar results are achieved if the symbionts do not respond to the host (i.e. if the constraint $`S_0=S_1`$ is applied). In order to maintain the cooperation of the symbionts the host strategy must be responsive. In a simulation which imposed both $`S_0=S_1`$ and $`H_0=H_1`$ the hosts and the symbionts both rapidly adopted strategies close to All Defect, with little payoff to either party. Similarly, if the host responds not to the individual symbionts but rather to an averaged level of cooperation, then mutual defection rapidly ensues.

## IV Three Moderating Influences

The argument given here makes a number of simplifying assumptions, some of which accentuate the degree of miserliness. In real systems we might expect somewhat less extreme outcomes for the following reasons.

Firstly, the model assumes that once symbionts enter a host they are trapped there, at least until the host dies. Clearly things will be different if they’re free to go. One can draw the analogy with sets of companies and their employees: for the worker there is a tradeoff between staying (with poor but secure working conditions) and going in search of greater prosperity with another employer, entailing some risk associated with the transfer (i.e. the costs and risks of transit itself, and the possibility that the new host/employer is even worse). On the other hand, this implies some selection pressure in favour of barriers or other ‘exit costs’ that hosts might be expected to impose on their workforces. For example, the ‘walling in’ of *Rhizobia* in the root nodules of legumes may merely enhance the uptake of fixed nitrogen by the plant, but it also could be argued that it plays a role in preventing the bacteria from leaving.

Secondly, all hosts in the simulation contain the same number of symbionts. One might instead expect the carrying capacity within more miserly hosts to be lower, indirectly decreasing their fitness. In general, hosts should evolve to be as miserly as possible while still retaining sufficient numbers of sufficiently compliant endosymbionts.

Finally, the model of selection employed here is a brutal one, and this makes the division between regions I and II (fig. 2) particularly stark. As a host approaches the miser strategy the gradient in fitness being followed by its endosymbiont population becomes flatter and flatter, and at the miser strategy itself the gradient is zero and so they spread out everywhere. If the endosymbionts aren’t responding strongly to the host strategy, the host’s gradient disappears as well. Again we have a trade-off: a given host wants to be as miserly as possible while ensuring that differential incentives continue to encourage cooperation strongly enough.

Against these points, however, the inheritance regime used in our simulations makes it difficult for hosts to evolve cooperation in their symbionts. Here each host begins life with an assortment of endosymbionts with randomly chosen strategies. Increased levels of symbiont cooperation are never passed on - they must be evolved anew with each generation. More realistic models might allow the new generation of hosts to begin with endosymbionts chosen from the generation before. This could take the form of vertical inheritance (in effect, perpetuating the ‘walling in’ across generations), or else symbionts might be expelled to the external environment and reacquired by new hosts. In either case, it will be easier to evolve full cooperation from symbionts, enhancing the tendency towards miserliness in hosts.

## V Discussion

A theme of recent work on symbiosis has been the surprising richness and complexity of the entanglement between host and symbiont, in which the partners undergo substantial metabolic and morphological changes mediated by a complex series of mutual signals, even modulating one another’s gene expression Baumann . It is clear that a model of the simplicity described here ignores many potentially relevant biological details. Its role is to show that, in principle, a low pay-off by the host may be sufficient to coerce a population of endosymbionts into nearly complete cooperation. An assumption of both the theoretical analysis and the model is that the host is able to respond to individual symbionts, preferentially rewarding those who cooperate. Perhaps the complexity of the real-life symbiotic dialogue concerns the host’s efforts to detect and respond at a local scale to the efforts of its symbiont passengers. The punishment of defecting symbionts may not require the host to differentiate individual symbionts, however. Suppose, for example, that a host species is well served by symbionts that convert some compound $`X`$ into some other compound $`Y`$. Even if the host is unable to reward the production of $`Y`$ at a local scale, it can punish endosymbionts that don’t carry out the reaction by flooding the area with something that’s toxic to them in the presence of $`X`$, or toxic in the absence of $`Y`$. A mechanism like this potentially allows a single host to interact with an essentially arbitrary number of symbionts at the same time.

The theory proposed here accounts for both the persistence of host-endosymbiont mutualisms and the observed asymmetry in their benefits, as well as suggesting a role for the complex host structures now being found in many associations. The very different dynamics found in a two-species interaction, when compared with the usual iterated Prisoner’s Dilemma, suggests that there is still much to be discovered about cooperative associations.

###### Acknowledgements.

The authors acknowledge the financial support of the Marsden Fund, which is administered by the Royal Society of New Zealand.