text:
# Epidemic dynamics in finite size scale-free networks

## Abstract

Many real networks present a bounded scale-free behavior with a connectivity cut-off due to physical constraints or a finite network size. We study epidemic dynamics in bounded scale-free networks with soft and hard connectivity cut-offs. The finite size effects introduced by the cut-off induce an epidemic threshold that approaches zero at increasing sizes. The induced epidemic threshold is very small even at a relatively small cut-off, showing that the neglection of connectivity fluctuations in bounded scale-free networks leads to a strong over-estimation of the epidemic threshold. We provide the expression for the infection prevalence and discuss its finite size corrections. The present work shows that the highly heterogeneous nature of scale-free networks does not allow the use of homogeneous approximations even for systems of a relatively small number of nodes.

In the last years it has been recognized that a large number of physical, biological, and social networks exhibits complex topological properties Amaral et al. (2000); Strogatz (2001). In particular many real world networks show the small-world phenomenon, related to a very small average path length between nodes Strogatz (2001); Watts and Strogatz (1998). More strikingly, in some cases this property is associated to a scale-free connectivity distribution, $`P(k)k^{2\gamma }`$, with $`0<\gamma 1`$, where $`k`$ is the number of links connected to a node Barabási and Albert (1999). This scale-free nature is associated to a large heterogeneity in the connectivity properties of the system. Since the second moment of the connectivity distribution $`k^2`$ is diverging when increasing the network size, the connectivity fluctuations in scale-free (SF) networks do not have an intrinsic bound and diverge in the infinite system size limit. Scale-free properties have been observed in many real systems such as the Internet Faloutsos et al. (1999); Caldarelli et al. (2000); Pastor-Satorras et al. (2001) and the World-Wide-Web Barabási and Albert (1999); Albert et al. (1999), food-webs, protein, and neural networks Albert and Barabási (2002). A very important example of scale-free networks is also found in the web of human sexual contacts Liljeros et al. (2001). This is a particularly relevant case since the unambiguous definition of contacts (links) is often missing in the analysis of social networks.

Since the Internet and the web of human sexual contacts appear to be scale-free, the study of epidemics and disease dynamics on SF networks is a relevant theoretical issue in the spreading of computer viruses and sexually transmittable diseases. In heterogeneous networks, it is well-known that the epidemic threshold decreases with the standard deviation of the connectivity distribution Hethcote and Yorke (1984); Anderson and May (1992), and this feature is amplified in SF networks, which have diverging connectivity fluctuations in the limit of infinite network size. Indeed, it was first noted in Ref. Pastor-Satorras and Vespignani (2001a) that, in infinite SF networks, epidemic processes do not possess an epidemic threshold below which diseases cannot produce a major epidemic outbreak or the inset of an endemic state. The absence of an intrinsic epidemic threshold has been found in both the susceptible-infected-susceptible (SIS) model Pastor-Satorras and Vespignani (2001a) and the susceptible-infected-removed (SIR) model May and Lloyd (2001); Moreno et al. (2001) in infinite SF networks. The immunization policies are as well very much affected by the SF nature of the connectivity distribution Pastor-Satorras and Vespignani (2002); Dezsö and Barabási (2001).

As customarily encountered in nonequilibrium statistical systems Marro and Dickman (1999), it has also been pointed out that in finite systems an epidemic threshold is induced by finite size effects May and Lloyd (2001). Real systems are actually made up by a finite number of individuals which is far from the thermodynamic limit. This finite population introduces a maximum connectivity $`k_c`$, depending on $`N`$, which has the effect of restoring a bound in the connectivity fluctuations, inducing in this way an effective nonzero threshold. More generally, we can consider a class of bounded scale-free (BSF) networks, in which the connectivity distribution has the form $`P(k)k^{2\gamma }f(k/k_c)`$, where the function $`f(x)`$ decreases very rapidly for $`x>1`$ Dorogovtsev and Mendes (2001). The cut-off $`k_c`$ can be due to the finite size of the network or to the presence of constraints limiting the addition of new links in an otherwise infinite network Amaral et al. (2000).

In this paper we present an analytical study of the SIS model in BSF networks with a generic connectivity exponent $`\gamma `$ ($`0<\gamma 1`$), focusing on the effects introduced by a finite cut-off $`k_c`$. We analyze the case of a hard cut-off, $`f(x)=\theta (1x)`$, where $`\theta (x)`$ is the Heaviside step function, as it happens in growing networks with a finite number of elements. We consider as well a soft exponential cut-off, $`f(x)=\mathrm{exp}(x)`$, as often found in systems where physical constraints are at play. We derive the behavior of the epidemic threshold as a function of $`k_c`$ and the network size $`N`$, and find that even for relatively small networks the induced epidemic threshold is much smaller than the epidemic threshold found in homogeneous systems. This confirms that the SF nature cannot be neglected in the practical estimates of epidemic and immunization thresholds in real networks. We also provide the explicit analytic form for the epidemic prevalence (density of infected individuals) in BSF networks. The results presented here can be readily extended to the SIR case.

In order to estimate the effect of $`k_c`$ in epidemics on BSF networks we will investigate the standard SIS model Diekmann and Heesterbeek (2000). This model relies on a coarse-grained description of individuals in the population. Namely, each node of the graph represents an individual and each link is a connection along which the infection can spread. Each susceptible (healthy) node is infected with rate $`\nu `$ if it is connected to one or more infected nodes. Infected nodes are cured and become again susceptible with rate $`\delta `$, defining an effective spreading rate $`\lambda =\nu /\delta `$ (without lack of generality, we set $`\delta =1`$). The SIS model does not take into account the possibility of individuals removal due to death or acquired immunization Diekmann and Heesterbeek (2000), and individuals run stochastically through the cycle susceptible $``$ infected $``$ susceptible. This model is generally used to study infections leading to endemic states with a stationary average density of infected individuals. In order to take into account the heterogeneity of SF networks, we have to relax the homogeneity assumption used in regular networks, and consider the relative density $`\rho _k(t)`$ of infected nodes with given connectivity $`k`$; i.e., the probability that a node with $`k`$ links is infected Pastor-Satorras and Vespignani (2001a). The dynamical mean-field equations can thus be written as

$$\frac{d\rho _k(t)}{dt}=\rho _k(t)+\lambda k\left[1\rho _k(t)\right]\mathrm{\Theta }(\rho (t)).$$
(1)
The first term in Eq. (1) considers infected nodes becoming healthy with unit rate. The second term represents the average density of newly generated infected nodes that is proportional to the infection spreading rate $`\lambda `$ and the probability that a node with $`k`$ links is healthy $`[1\rho _k(t)]`$ and gets the infection via a connected node. The rate of this last event is given by the probability $`\mathrm{\Theta }(\rho (t))`$ that any given link points to an infected node, which has the expression Pastor-Satorras and Vespignani (2001a)

$$\mathrm{\Theta }(\rho (t))=k^1\underset{k}{}kP(k)\rho _k(t).$$
(2)
By solving Eqs. (1) and (2) in the stationary state \[$`d\rho _k(t)/dt=0`$\] we obtain the self-consistency equation Pastor-Satorras and Vespignani (2001a)

$$\mathrm{\Theta }=k^1\underset{k}{}kP(k)\frac{\lambda k\mathrm{\Theta }}{1+\lambda k\mathrm{\Theta }},$$
(3)
where $`\mathrm{\Theta }`$ is now a function of $`\lambda `$ alone. The self-consistency Eq. (3) allows a solution with $`\mathrm{\Theta }0`$ and $`\rho _k0`$ only if the condition $`\lambda k^2/k1`$ is fulfilled Pastor-Satorras and Vespignani (2002), defining the epidemic threshold

$$\lambda _c=\frac{k}{k^2}.$$
(4)
In other words, if the value of $`\lambda `$ is above the threshold, $`\lambda \lambda _c`$, the infection spreads and becomes endemic. Below it, $`\lambda <\lambda _c`$, the infection dies out exponentially fast. This result implies that in infinite SF networks with connectivity exponent $`0<\gamma 1`$, for which $`k^2\mathrm{}`$, we have $`\lambda _c=0`$. This fact implies in turn that for any positive value of $`\lambda `$ the infection can pervade the system with a finite prevalence, in a sufficiently large network Pastor-Satorras and Vespignani (2001a). While this results is valid for infinite SF networks, $`k^2`$ assumes a finite value in BSF networks, defining an effective nonzero threshold due to finite size effects as usually encountered in nonequilibrium phase transitions Marro and Dickman (1999). This epidemic threshold, however, is not an intrinsic quantity as in homogeneous systems and it vanishes for a increasing network size or connectivity cut-off. In order to calculate the precise effects of a finite $`k_c`$ we consider two different cases of connectivity cut-off. At first instance we consider a soft exponential cut-off with characteristic connectivity $`k_c`$. This case corresponds to those real networks in which external factors set up an upper limit to the connectivity Amaral et al. (2000). The network can have an infinite number of elements but the power-law connectivity distribution decays exponentially for large values of $`k`$. In order to perform explicit calculations we use a continuous approximation that substitutes the connectivity by a real variable $`k`$ in the range $`[m,\mathrm{})`$, where $`m`$ is the minimum connectivity of the network. The connectivity probability distribution in this case is $`P(k)=Ak^{2\gamma }\mathrm{exp}(k/k_c)`$, where $`A`$ is a normalization factor. The effective nonzero epidemic threshold $`\lambda _c(k_c)`$ induced by the exponential cut-off is given by

$$\lambda _c(k_c)=\frac{_m^{\mathrm{}}k^{1\gamma }\mathrm{exp}(k/k_c)𝑑k}{_m^{\mathrm{}}k^\gamma \mathrm{exp}(k/k_c)𝑑k},$$
(5)
which, after integration, yields

$$\lambda _c(k_c)=k_c^1\frac{\mathrm{\Gamma }(\gamma ,m/k_c)}{\mathrm{\Gamma }(1\gamma ,m/k_c)},$$
(6)
where $`\mathrm{\Gamma }(x,y)`$ is the incomplete Gamma function Abramowitz and Stegun (1972). For large $`k_c`$ we can perform a Taylor expansion and retain only the leading term, obtaining for any $`0<\gamma <1`$

$$\lambda _c(k_c)\frac{1}{m\gamma \mathrm{\Gamma }(1\gamma )}(k_c/m)^{\gamma 1}.$$
(7)
The limit $`\gamma 1`$ in Eq. (6) corresponds to a logarithmic divergence, yielding at leading order $`\lambda _c(k_c)(m\mathrm{ln}(k_c/m))^1`$. In all cases we have that the epidemic threshold vanishes when increasing the characteristic cut-off. For large $`k_c`$, the average connectivity is virtually fixed and given by $`k=(\gamma +1)m/\gamma `$, for any $`\gamma >0`$. It is interesting thus to compare the intrinsic epidemic threshold obtained in homogeneous networks with negligible fluctuations and the nonzero effective threshold of BSF networks. The intrinsic epidemic threshold of homogeneous networks with constant node connectivity $`k`$ is given by $`\lambda _c^\mathrm{H}=k^1`$ Diekmann and Heesterbeek (2000); Pastor-Satorras and Vespignani (2001a). If we compare BSF and homogeneous networks with the same average connectivity $`k=(\gamma +1)m/\gamma `$ we obtain that the ratio between the epidemic thresholds is given by

$$\frac{\lambda _c(k_c)}{\lambda _c^\mathrm{H}}\frac{(\gamma +1)}{\gamma ^2\mathrm{\Gamma }(1\gamma )}(k_c/m)^{\gamma 1}.$$
(8)
This clearly shows that even in the case of a connectivity cut-off the effective epidemic threshold in BSF networks is much smaller than the intrinsic threshold obtained in regular networks. In Fig. 1 we plot the ratio obtained by using the full expression for $`\lambda _c(k_c)`$, Eq. (6). It is striking to observe that, even with relatively small cut-offs ($`k_c10^210^3`$), for $`\gamma 0.5`$ the effective epidemic threshold of BSF networks is smaller by a factor close to $`1/10`$ than the intrinsic threshold obtained on homogeneous networks.

As a second kind of finite size effect we consider the presence of a hard cut-off $`k_c`$. Since SF networks are often dynamically growing networks, this case represent a network which has grown up to a finite number of nodes $`N`$. The maximum connectivity $`k_c`$ of any node is related to the network age, measured as the number of nodes $`N`$, by the scaling relation Dorogovtsev and Mendes (2001)

$$k_cmN^{1/(1+\gamma )},$$
(9)
where $`m`$ is the minimum connectivity of the network.

In this case the network does not possess any node with connectivity $`k`$ larger than $`k_c`$, and we can think in terms of a hard cut-off. Using again the continuous $`k`$ approximation, the normalized connectivity distribution has now the form

$$P(k)=\frac{(1+\gamma )m^{1+\gamma }}{1(k_c/m)^{1\gamma }}k^{2\gamma }\theta (k_ck),$$
(10)
where $`\theta (x)`$ is the Heaviside step function. The finite size induced epidemic threshold $`\lambda _c(k_c)`$ is given by the expression

$$\lambda _c(k_c)=\frac{_m^{k_c}k^{1\gamma }𝑑k}{_m^{k_c}k^\gamma 𝑑k}.$$
(11)
Evaluating the above expression we obtain at leading order in $`k_c/m`$:

$$\lambda _c(k_c)\frac{1\gamma }{\gamma m}(k_c/m)^{\gamma 1}.$$
(12)
In this case the hard cut-off $`k_c`$ can be expressed as a function of the network size $`N`$ by using the scaling relation Eq. (9) and we can obtain the effective epidemic threshold as

$$\lambda _c(N)\frac{1\gamma }{\gamma m}N^{(\gamma 1)/(\gamma +1)}.$$
(13)
This expression is valid for any $`0<\gamma <1`$, while for $`\gamma =1`$ we obtain at the leading order the logarithmic behavior $`\lambda _c(N)2(m\mathrm{ln}(N))^1`$. Also in this case we have that the effective epidemic threshold is approaching zero for increasing network sizes, and it is worth comparing its magnitude with the corresponding intrinsic threshold in homogeneous networks with identical average connectivity. In Fig. 2 we report the ratio $`\lambda _c(N)/\lambda _c^\mathrm{H}`$ for different sizes of the SF network. It is striking to notice that for $`\gamma =0.5`$, small networks with $`N10^4`$ exhibit a finite size induced epidemic threshold that is close to be one order of magnitude smaller than the intrinsic epidemic threshold of a homogeneous network.

In order to find the prevalence behavior we have to solve Eq. (3) in the continuous approximation,

$$\mathrm{\Theta }=k^1_m^{\mathrm{}}kP(k)\frac{\lambda \mathrm{\Theta }k}{1+\lambda \mathrm{\Theta }k}𝑑k,$$
(14)
and use the value of $`\mathrm{\Theta }`$ to compute the density of infected sites $`\rho `$ as

$$\rho =\underset{k}{}\rho _kP(k)_m^{\mathrm{}}P(k)\frac{\lambda \mathrm{\Theta }k}{1+\lambda \mathrm{\Theta }k}𝑑k,$$
(15)
where $`P(k)`$ is given by Eq. (10). In the absence of any cut-off ($`k_c\mathrm{}`$) and in the thermodynamic limit ($`N\mathrm{}`$) the prevalence scales as $`\rho \lambda ^{1/(1\gamma )}`$ if $`0<\gamma <1`$, and as $`\rho \mathrm{exp}(1/m\lambda )`$ if $`\gamma =1`$ Pastor-Satorras and Vespignani (2001a). Accordingly with the absence of the epidemic threshold, the prevalence is null only if the spreading rate is $`\lambda =0`$. In the case of a hard cut-off we can integrate Eq. (14), neglecting terms of order $`(k_c/m)^\gamma `$ in the $`P(k)`$ distribution, to obtain:

$`\mathrm{\Theta }`$ $``$ $`\gamma m^\gamma \lambda \mathrm{\Theta }{\displaystyle _m^{k_c}}{\displaystyle \frac{k^\gamma }{1+\lambda \mathrm{\Theta }k}}𝑑k`$
$`=`$ $`F(1,\gamma ,1+\gamma ,[\lambda \mathrm{\Theta }m]^1)`$
$`(k_c/m)^\gamma F(1,\gamma ,1+\gamma ,[\lambda \mathrm{\Theta }k_c]^1),`$
where $`F`$ is the Gauss hyper-geometric function Abramowitz and Stegun (1972). For a fixed $`k_c`$ one can expand both hyper-geometric functions on the right hand side in the previous equation, keeping the most relevant terms in $`\mathrm{\Theta }`$ and considering afterwards the limit of large $`k_c`$. The final solution for $`\mathrm{\Theta }`$ is then given, at leading order in $`(k_c/m)`$, by

$$\mathrm{\Theta }\frac{1}{m\lambda ^2}\frac{2\gamma }{1\gamma }\left(\frac{k_c}{m}\right)^1\left[\lambda \frac{1\gamma }{\gamma m}\left(\frac{k_c}{m}\right)^{\gamma 1}\right].$$
(16)
By evaluating the integral in Eq. (15) and keeping the leading term in $`\mathrm{\Theta }`$ and $`k_c`$ we finally obtain the infection prevalence as

$$\rho \frac{(\gamma +1)(2\gamma )}{\lambda \gamma (1\gamma )}\left(\frac{k_c}{m}\right)^1\left[\lambda \frac{1\gamma }{\gamma m}\left(\frac{k_c}{m}\right)^{\gamma 1}\right].$$

Inserting the scaling relation Eq. (9) between the maximum connectivity $`k_c`$ and the network size $`N`$ we are led to the final expression

$$\rho N^{1/(\gamma +1)}(\lambda \lambda _c(N)).$$
(17)
That is, the finite size of the network induces a standard mean-field transition at the induced epidemic threshold $`\lambda _c(N)`$, given by Eq. (13). As can be seen from Eq. (17), however, the prevalence is depressed by a factor $`N^{1/(\gamma +1)}`$ from the corresponding value for a homogeneous network. The above calculations can be repeated along similar lines in the case of a soft exponential cut-off, obtaining similar results.

It is worth remarking that similar results hold as well for the SIR model. Despite this model confers permanent immunity and does not allow for a stationary state, the epidemic threshold over which an epidemic outbreak occurs has the same analytic form $`\lambda _c=k/k^2`$ Moreno et al. (2001). Thus, the present results for the effect of finite size and the induced epidemic threshold can be readily exported to the SIR case. The calculation of the epidemic prevalence is different due to the different evolution equations, but recovers the same onset of an induced mean-field transition at the effective threshold $`\lambda _c(N)`$. Interestingly, similar results have been obtained in the analysis of the resilience to damage of finite size scale-free networks Cohen et al. (2001); Callaway et al. (2000)

In conclusion, we have shown that the SF networks weakness to epidemic agents is also present in finite size networks. Using the homogeneity assumption in the case of SF networks will lead to a serious over-estimate of the epidemic threshold even for relatively small networks.

This work has been partially supported by the European Network Contract No. ERBFMRXCT980183. R.P.-S. acknowledges financial support from the Ministerio de Ciencia y Tecnología (Spain) and from the Abdus Salam International Centre for Theoretical Physics (ICTP), where part of this work was done. We thank F. Cecconi for helpful discussions.


format: text

text:
# Constraints on the Extremely High-Energy Cosmic Ray Accelerators from Classical Electrodynamics

## I Introduction

The existence of EHECRs with energies above $`10^{20}`$ eV is now a well-established fact (see, e.g. rev1  and rev2  for a review). Nearly perfect isotropy of the highest energy cosmic rays points to their likely extragalactic origin. On the other hand, the sources must be located within $`100`$ Mpc: the propagation distance for EHECRs of energy greater than $`5\times 10^{19}`$ eV is limited due to their interaction with the microwave background radiation Gr66 ; ZK , unless there are some yet unknown type(s) of particles that do not undergo such an interaction.

The presence of non-thermal particles is a usual feature of any plasma system. However, finding a way to accelerate these particles to $`10^{20}`$ eV is a serious theoretical challenge. Even if an electrodynamical system used to accelerate charged particles is sufficiently long-lived, the attainable energy is limited since the particles can only stay within the acceleration region until their gyroradius exceeds the size of the region. (This is known as Hillas criterion H84 .) Various types of model- and source-dependent radiative and collisional losses decrease the limiting particle energy furthermore. But even if one finds a way to circumvent the latter obstacle, there are radiative losses caused by the same electromagnetic field that is used to accelerate particles. For example, Greisen Gr65  considered the synchrotron losses and showed that for very small-sized accelerators, i.e., for those having very strong magnetic field, the synchrotron radiation drains energy so fast that even a particle accelerated at the largest possible rate for the magnetic-field dominated environments, $`qBc`$ ($`q`$ is the particle charge, $`B`$ the magnetic field strength, $`c`$ the velocity of light), is not able to achieve the energy given by the Hillas criterion. However, Greisen’s limit might seem to be not a general one, as in some objects accelerated particles move along the magnetic field lines (neutron stars are a perfect example of that), without suffering from the synchrotron losses.

In this paper, we relax the requirements to the acceleration sites as much as possible, assuming that inhomogeneity of electric and magnetic fields is small enough not to cause particle diffusion out of the acceleration site and taking into account only the radiation mechanism which gives the lowest energy loss rate. In fact, even if the synchrotron losses (or the electro-bremsstrahlung in the case of electric-field dominated environments) tend to zero, i.e., for particles moving along smooth field lines, there is still curvature radiation, which then limits the maximal particle energy. Furthermore, we generalize the analysis introducing acceleration efficiency parameter $`\eta `$, so that the acceleration rate is given by

$$\dot{\epsilon }=\eta qBc,$$
(1)
where $`\eta B=E_{\mathrm{eff}}`$ is the projection of an electric field on the particle trajectory, effectively averaged as the particle moves along this trajectory. Such a representation is useful when, e.g., diffusive shock acceleration scenario is brought into consideration.

We derive the ultimate lower limit for the total electromagnetic energy stored in an acceleration region and find the size and the field strength that are optimal with respect to minimization of this energy. Then we consider a possible gain from relativistic motion of an accelerator as a whole and discuss what kind of additional information one gets from analysis of hard $`\gamma `$-photons generated as a byproduct of cosmic ray acceleration. Finally, we make comparison of theoretical constraints with the parameters of known astrophysical objects.

## II Electrodynamical limitations

Let us consider a region of space of a size $`R`$ where particle acceleration occurs. A natural assumption is that there is nothing like linear accelerator geometry in this region, but instead the field lines of both electric and magnetic fields are curved with the curvature radius being of the order of $`R`$ or less.

If the magnetic field is much stronger than the electric field, $`BE`$, then a test particle is accelerated following a magnetic field line in its trajectory. If $`EB`$ but $`(\stackrel{}{E}\stackrel{}{B})EB`$, then particle’s trajectory depends on the radiative loss rate. High loss rate will force the particle to follow a field line (of an electric field in this case) since a transverse to the field line component of particle’s momentum will quickly decay due to an increase in radiative losses. In the opposite case, the field geometry will have little effect on particle’s trajectory which therefore will make an angle of $`45^{}`$ with field lines. The remaining case, where both Lorentz invariants $`(\stackrel{}{E}\stackrel{}{B})`$ and $`|E^2B^2|`$ are much smaller than $`E^2+B^2`$, can be reduced to one of the already mentioned situations by the appropriate Lorentz transformation to the reference frame in which the mean (after averaging in space and time) Poynting flux density becomes zero.

In what follows we assume that accelerator’s lifetime is at least as large as the acceleration timescale and use the most optimistic premise that the strongest of magnetic and electric fields is regular on scales $`R`$ and quasistatic. (Alternating field direction implies that particles escape more readily and that the time-averaged acceleration rate is smaller, while the radiative losses are roughly the same.) Then, the effect of small-scale magnetic field can be neglected if $`E>B`$, and in the case $`B>E`$ the factor $`\eta <1`$ takes into account all necessary averaging for a small-scale and/or variable electric field, so that $`E_{\mathrm{eff}}=\eta B`$.

The particle can be accelerated up to the terminal energy which is the smallest of either the work done by accelerating force in a time it takes for a particle to escape from the region, or the energy at which radiative losses balance the acceleration. The balance condition requires that the radiative loss rate $`\dot{\epsilon }_{\mathrm{rad}}`$ is equal to the energy gain rate, i.e.

$$\dot{\epsilon }_{\mathrm{rad}}=\frac{2}{3}\gamma ^4\frac{q^2}{R^2}c=\eta qBc$$
(2)
in the case of curvature-radiation dominated losses Ginz  (the particle is assumed to move along an arc of radius $`R`$), or

$$\dot{\epsilon }_{\mathrm{rad}}=\frac{2}{3}\gamma ^2\left(\frac{q^2}{mc^2}\right)^2c\left(B_{}^2+E_{}^2\right)=\eta qBc$$
(3)
in the case of synchrotron- or electro-bremsstrahlung dominated losses. Here $`\gamma `$ is the Lorentz factor of the particle, $`m`$ its mass, $`B_{}`$ and $`E_{}`$ are the r.m.s. field components perpendicular to the particle momentum.

The comparison of Eqs. (2) and (3) shows that the curvature radiation gives a more favorable (smaller) estimate for the radiative losses as long as the particle energy, $`\epsilon =\gamma mc^2`$, satisfies the condition

$$\epsilon <qR\sqrt{B_{}^2+E_{}^2}.$$
(4)
If the above condition were violated, it would mean one of the following. Either the particle is no longer bound in the acceleration region, i.e., its gyro-radius $`r_\mathrm{H}=\epsilon /(qB_{})`$ is larger than $`R`$ (when $`r_\mathrm{H}R`$, the very idea of curvature radiation becomes meaningless) and acceleration terminates upon the particle escape at $`\epsilon qRB_{}`$, or, if $`B_{}E_{}`$, the particle is unbound from the very beginning and its energy gain is limited by the difference of electric potentials across the acceleration region, $`\epsilon qRE_{\mathrm{eff}}=\eta qRB`$. This provides a generalized version of the Hillas criterion.

Hereafter we neglect numerical factors of $`1`$ assuming $`B=B_{}`$ and $`E=E_{}=E_{\mathrm{eff}}=\eta B`$. In general, one should bear in mind that for an inhomogeneous small-scale field the r.m.s. quantity $`E^2`$ must be replaced with $`E^2=E^2/\delta `$, where $`\delta <1`$ is a duty-factor. (The same is true for magnetic field, although the situations where $`E>B`$ and the magnetic field has the duty-factor $`\delta 1`$ seem quite unrealistic.) We assume below $`\delta =1`$, aiming to the most efficient acceleration. Thus, the condition given by Eq. (4) cannot be violated, as it satisfies the generalized Hillas criterion,

$$\epsilon _{\mathrm{max}}<qR\mathrm{max}\{B,E\}.$$
(5)
Since we are interested in the most efficient acceleration (i.e., the highest acceleration rate at smallest radiative losses), it is Eq. (2) that will be used below to put the lower limit on the radiative loss rate. Then, we conclude that the maximal energy of a particle is either determined by the balance condition given by Eq. (2),

$$\epsilon _{\mathrm{max}}^4=\frac{3\eta BR^2}{2q}\left(mc^2\right)^4,$$
(6)
or is limited by the Eq. (5), whichever gives the smaller value.

The overall energy of an electromagnetic field in a spherical region of radius $`R`$, capable of accelerating particles up to $`\epsilon _{\mathrm{max}}`$, is $`W_{\mathrm{em}}=R^3(B^2+E^2)/6`$. According to Eqs. (6) and (5), respectively, it is limited as follows:

$$W_{\mathrm{em}}>\frac{2}{27}\frac{q^2}{R}\left(\frac{\epsilon _{\mathrm{max}}}{mc^2}\right)^8\frac{1+\eta ^2}{\eta ^2}$$
(7a)
and
$$W_{\mathrm{em}}>\frac{R}{6}\left(\frac{\epsilon _{\mathrm{max}}}{q}\right)^2,$$
(7b)
Our interest is to lower the energy requirements taking $`R`$ as large as possible in Eq. (7a) and as small as possible in Eq. (7b). Doing this, we derive the optimal size

$$R^{(\mathrm{opt})}\frac{2}{3}\frac{\sqrt{1+\eta ^2}}{\eta }\frac{q^2\epsilon _{\mathrm{max}}^3}{\left(mc^2\right)^4},$$
(8)
such that both conditions become quantitatively the same. This finally gives the optimal (the minimum possible) estimate for the amount of electromagnetic energy stored in the acceleration region:

$$W_{\mathrm{em}}^{(\mathrm{opt})}\frac{1}{9}\frac{\sqrt{1+\eta ^2}}{\eta }\frac{\epsilon _{\mathrm{max}}^5}{\left(mc^2\right)^4}.$$
(9)
The corresponding optimal strength of the magnetic field is

$$B^{(\mathrm{opt})}\frac{3}{2}\frac{\eta }{1+\eta ^2}\frac{\left(mc^2\right)^4}{q^3\epsilon _{\mathrm{max}}^2},$$
(10)
and $`E^{(\mathrm{opt})}=\eta B^{(\mathrm{opt})}`$.

The energy requirements given by Eqs. (7) as functions of accelerator size are presented in Fig. 1. The required accelerator energy increases as the region size deviates from the optimal one. A smaller region would accelerate particles to the energy $`\epsilon _{\mathrm{max}}`$ only if the energy budget is increased by $`R^{(\mathrm{opt})}/R`$ times; in a larger region, the energy budget has to be increased by $`R/R^{(\mathrm{opt})}`$ times.

To have an idea of typical numbers, let us make the estimates for a proton of energy $`\epsilon _{\mathrm{max}}=10^{20}`$ eV. In this case, the total electromagnetic energy in the acceleration region should exceed $`W_{\mathrm{em}}^{(\mathrm{opt})}3\times 10^{51}`$ erg; the optimal size $`R^{(\mathrm{opt})}`$ is of the order of $`10^{17}`$ cm and the optimal magnetic field is of $`3`$ G. When making these estimates we assumed $`\eta =1`$ to lower the energy requirements. Taking $`\eta `$ still larger does not reduce them too much, by a factor of 2 at most.

All the results obtained are valid also for an accelerator that moves as a whole with the Lorentz factor $`\mathrm{\Gamma }1`$, if the quantities are measured in the comoving frame. However, it is more convenient to present the results in a form where $`W_{\mathrm{em}}^{(\mathrm{opt})}`$ and $`\epsilon _{\mathrm{max}}`$ are measured in the laboratory frame, while $`R^{(\mathrm{opt})}`$, $`B^{(\mathrm{opt})}`$ and $`E^{(\mathrm{opt})}`$ are measured in the comoving frame. Then

$$W_{\mathrm{em}}^{(\mathrm{opt})}\frac{1}{9\mathrm{\Gamma }^4}\frac{\sqrt{1+\eta ^2}}{\eta }\frac{\epsilon _{\mathrm{max}}^5}{\left(mc^2\right)^4},$$
(11a)
$$R^{(\mathrm{opt})}\frac{2}{3\mathrm{\Gamma }^3}\frac{\sqrt{1+\eta ^2}}{\eta }\frac{q^2\epsilon _{\mathrm{max}}^3}{\left(mc^2\right)^4},$$
(11b)
$$B^{(\mathrm{opt})}\frac{3\mathrm{\Gamma }^2}{2}\frac{\eta }{1+\eta ^2}\frac{\left(mc^2\right)^4}{q^3\epsilon _{\mathrm{max}}^2},$$
(11c)
where the primed quantities are measured in the comoving frame. A comparison of energy requirements given by Eqs. (7), generalized for the case $`\mathrm{\Gamma }1`$, for accelerators with different bulk Lorentz factors is presented in Fig. 2.

Apparently, the acceleration of EHECRs in ultrarelativistic bulk flows has an advantage of bringing the optimal source size to a more ”comfortable” (for short-lived phenomena like Gamma-Ray Bursts) range and reducing the energy requirements to the source. However, for wind-like relativistic flows the actual requirements are geometry-dependent. Indeed, for the causality reasons, the acceleration region does not occupy the whole sphere of radius $`R`$, but rather extends to a distance $`R^{}=R/\mathrm{\Gamma }`$ transverse to the radius and to a distance $`R^{}/\mathrm{\Gamma }`$ along it, so that the total energy within the volume of radius $`R`$ is of $`\mathrm{\Gamma }^4W_{\mathrm{em}}`$ for a wind with a broad beam pattern. Since the energy stored in the acceleration region, $`W_{\mathrm{em}}`$, is multiplied by $`\mathrm{\Gamma }^4`$ in this case, one can only gain from a more favorable ratio $`R^{}/R^{(\mathrm{opt})}`$ (cf. Eq. (11a)). A jet geometry is more favorable: with the jet opening angle $`1/\mathrm{\Gamma }`$ the energy requirements can be reduced by a factor as large as $`\mathrm{\Gamma }^2`$, with respect to the non-relativistic accelerator. To implement the advantage of relativistic bulk motion to the full extent, i.e., to get a gain by a factor of $`\mathrm{\Gamma }^4`$, one needs a transient jet source which lasts for a time $`R/(\mathrm{\Gamma }^2c)`$.

## III Constraints due to radiation from acceleration sites

Before moving to an overview of different source populations, let us discuss some implications regarding the acceleration timescale and the electromagnetic radiation accompanying the acceleration process.

When $`\eta 1`$, a particle is accelerated up to $`\epsilon _{\mathrm{max}}`$ in a time $`t_{\mathrm{acc}}R/c`$, which is less than or comparable to the dynamical timescale of an accelerator. However, at low efficiency $`\eta 1`$, the acceleration may take a longer time, i.e. $`\eta ^1(R/c)`$, if accelerator’s size is larger than $`R^{(\mathrm{opt})}`$, or $`\eta ^1(R/R^{(\mathrm{opt})})(R/c)`$ otherwise. Relativistic motion of an accelerator usually implies relativistic expansion. Consequently, the lifetime is of order $`R/c`$, so that small $`\eta `$ strongly disfavors such objects unless their size is $`\eta R^{(\mathrm{opt})}`$, that would allow a particle to attain the energy $`\epsilon _{\mathrm{max}}`$.

The acceleration of particles goes side by side with the emission of hard electromagnetic radiation. The lower limit on the energy lost to this radiation is given by the product of $`t_{\mathrm{acc}}`$ and the radiative loss rate (Eq. (2), where it is convenient to substitute $`\gamma ^3`$ from Eq. (8)). However, if $`t_{\mathrm{acc}}<R/c`$ one must integrate the loss rate over a time of order $`R/c`$, since a particle cannot leave the acceleration region in a shorter time. As a result of particle’s acceleration, the following fraction of its terminal energy $`\epsilon _{\mathrm{max}}`$ is lost for $`\gamma `$-radiation

$$\frac{E_{\mathrm{rad}}}{\epsilon _{\mathrm{max}}}\{\begin{array}{ccc}R^{(\mathrm{opt})}/R\hfill & \text{if}\hfill & R>R^{(\mathrm{opt})}\hfill \\ 1\hfill & \text{if}\hfill & R^{(\mathrm{opt})}>R>\eta R^{(\mathrm{opt})}\hfill \\ \eta R^{(\mathrm{opt})}/\sqrt{1+\eta ^2}R\hfill & \text{if}\hfill & R<\eta R^{(\mathrm{opt})}.\hfill \end{array}$$
(12)
Therefore, acceleration in small objects ($`R<R^{(\mathrm{opt})}`$) is always limited by radiative losses.

The spectral energy distribution of the accompanying radiation has a maximum (in the source frame) around

$$\epsilon _\gamma \mathrm{}\frac{\gamma ^3c}{R}\eta \frac{R^{(\mathrm{opt})}}{R}\frac{mc^2}{\alpha Z^2},$$
(13)
provided the accelerator size is less than $`R^{(\mathrm{opt})}`$; $`\alpha `$ is the fine-structure constant and $`Z`$ the ion charge in units of the elementary charge. The energy of curvature radiation photons given by Eq. (13) is related to the energy of synchrotron photons as $`\epsilon _\gamma =(r_\mathrm{H}/R)\epsilon _\gamma ^{\mathrm{sy}}`$.

The ratio of the luminosity in accompanying radiation to the luminosity in accelerated particles, $`L_\gamma /L_{\mathrm{acc}}`$, as well as the energy of individual photons, is reduced if the acceleration time is limited to $`R/c`$, like in the sources with relativistic bulk motion (particle diffusion may lead to analogous limit), but this is at expense of reducing $`\epsilon _{\mathrm{max}}`$. Since the radiative loss rate (Eq. (2)) and the energy of individual photons (Eq. (13)) scale as $`\gamma ^4`$ and $`\gamma ^3`$, while particle’s Lorentz factor increases linearly with time, the corresponding factors for $`L_\gamma /L_{\mathrm{acc}}`$, $`\epsilon _\gamma `$ and $`\epsilon _{\mathrm{max}}`$ are: $`(\eta R^{(\mathrm{opt})}/R)^4`$, $`(\eta R^{(\mathrm{opt})}/R)^3`$, and $`(\eta R^{(\mathrm{opt})}/R)`$ respectively, in the case where $`R^{(\mathrm{opt})}>R>\eta R^{(\mathrm{opt})}`$. Note that if the maximal energy were limited by the synchrotron rather than curvature losses, then the acceleration up to $`\epsilon _{\mathrm{max}}`$ would not be possible at all in the objects smaller than $`R^{(\mathrm{opt})}`$ if $`\eta 1`$, as noted in Gr65 .

The check whether the sky-averaged luminosity in accompanying electromagnetic emission exceeds the diffuse background is a powerful selection tool for possible acceleration sites. For example, the acceleration of protons with high efficiency, $`\eta 1`$, results in $`\gamma `$-radiation in energy range $`100`$ GeV (Eq. (13)). Such photons undergo reprocessing into $`10`$ MeV - $`100`$ GeV energy range through interaction with cosmic infrared and microwave backgrounds. The observed spectral energy distribution of diffuse $`\gamma `$-ray background, $`10^9`$ erg/(cm2 s sr)-1 gamma-back , is virtually constant in the $`10`$ MeV - $`100`$ GeV spectral domain, and the observed flux of $`10^{20}`$ eV cosmic rays is $`25\times 10^{12}`$ erg/(cm2 s sr)-1 rev2 , i.e., $`200500`$ times less. But the GZK effect limits the distance to which the sources of $`10^{20}`$ eV cosmic rays can be observed to $`30`$ Mpc, while the related $`\gamma `$-radiation is visible up to a redshift $`z1`$, i.e. the distance $`3`$ Gpc. If the accelerators are uniformly distributed in the Universe and no other contributions to the $`\gamma `$-ray background are taken into account, this implies $`L_\gamma /L_{\mathrm{acc}}5`$ and hence $`R/(\eta R^{(\mathrm{opt})})1/5`$. Since the spectral energy distribution of the $`\gamma `$-ray background is nearly constant in $`10`$ MeV - $`100`$ GeV range, the same limit applies for iron nuclei, although their individual $`\gamma `$-luminosity is concentrated in the lower-energy spectral domain $`10`$ GeV (see Eq. (13), which gives scaling $`\epsilon _\gamma A/Z^2`$).

The accompanying radiation may also help in identifying acceleration sites. The electrons are more suitable for this purpose as they can achieve radiative-loss limited regime of acceleration in an object with much smaller energy content. For example, there is an excess emission in $`100`$ MeV range from the Crab plerion, which was interpreted as the electron synchrotron radiation crab100 . Interestingly, the size of the plerion, $`10^{17}`$ cm, and its estimated magnetic field, $`10^4`$ G, are very close to the optimal values for electrons accelerated at the efficiency $`\eta =1`$. The electrons then attain the energy $`3\times 10^{15}`$ eV and the peak of their synchrotron spectrum is at $`100`$ MeV.

So far, we did not take into account the proton interaction with smallest-scale electromagnetic field, which results in photomeson reactions (at wavelengths $`<0.04`$ cm or photon energies $`>2\times 10^3`$ eV for $`10^{20}`$ eV proton). Such an interaction, with the cross-section $`\sigma _{p\gamma }2\times 10^{28}`$ cm2 and inelasticity $`f0.5`$, may in some cases lead to unacceptably high energy losses for the accelerated particles. For illustration, we make the estimate for $`10^{20}`$ eV protons and for the most favorable case $`\eta =1`$. In order to allow accelerated protons to escape without significant losses, the luminosity of the source at the maximum of integral photon spectrum must be

$$L<\pi \frac{Rc}{f\sigma _{p\gamma }}\overline{\epsilon }1.5\times 10^{44}\left(\frac{R}{10^{17}\text{cm}}\right)\left(\frac{\overline{\epsilon }}{1\text{eV}}\right)\text{erg/s},$$
(14)
where $`\overline{\epsilon }`$ is the location of the maximum.

## IV Overview of possible acceleration sites

Neutron stars. Because of their size, which is far too small compared to $`R^{(\mathrm{opt})}`$, the performance of neutron stars as cosmic ray accelerators is limited be the curvature losses (see, e.g., Ginz ). For $`10^{20}`$ eV EHECRs, the energy requirements to neutron star magnetospheres are unphysical: $`10^{62}`$ erg for protons and $`10^{51}`$ erg for iron nuclei. Acceleration of heavier nuclei, with mass $`A150`$ and charge $`Z50`$, to slightly smaller energy $`5\times 10^{19}`$ eV is marginally possible for rapidly rotating (i.e. $`\eta 1`$) neutron stars with magnetic field $`B10^{15}`$ G (so-called magnetars). However, the question whether such nuclei can be efficiently transported to neutron star surface and ejected from it remains open. In addition, the luminosity in accompanying $`10100`$ TeV curvature $`\gamma `$-radiation (which is likely to be reprocessed through electromagnetic cascade into $`1100`$ MeV range) would be 5 to 6 orders of magnitude higher than in the produced EHECRs. Particles accelerated near the light cylinder produce less accompanying radiation, but the energy requirements rise in this case. Therefore, magnetospheres of neutron stars have to be excluded from the list of EHECR sources.

A way out from the problem of catastrophic overproduction of the accompanying radiation is acceleration of particles within ultrarelativistic pulsar winds. However, the energy constraints, which now have to be reformulated in terms of required Poynting-flux luminosity $`L_{\mathrm{em}}\mathrm{\Gamma }^4W_{\mathrm{em}}c/R`$ for a wind with a broad beam pattern, remain extremely tough. The best case ($`R^{}=R^{(\mathrm{opt})}`$ and $`\eta =1`$) for $`10^{20}`$ eV protons corresponds to $`\mathrm{\Gamma }600`$ and $`L_{\mathrm{em}}10^{45}\mathrm{\Gamma }^2`$ erg/s. With iron nuclei the figures become $`\mathrm{\Gamma }60`$ and $`L_{\mathrm{em}}1.510^{42}\mathrm{\Gamma }^2`$ erg/s. In any case, the required luminosity is far beyond any of the observed values.

Black holes. Rapidly rotating black holes embedded in a magnetic field can generate electric field and accelerate particles (e.g., BHacc ; BHacc1 ). We make an (optimistic) estimate assuming that the magnetic field is roughly in equipartition with accreted plasma, an accretion rate corresponds to the Eddington regime, and acceleration occurs within $`3R_g`$, where $`R_g3\times 10^5`$ cm $`M/M_{}`$ is the Schwarzschild radius of the black hole of mass $`M`$. The energy $`W_{\mathrm{em}}\pi m_pc^2R_g^2/\sigma _T`$ is available in the acceleration region ($`m_p`$ is the proton mass, and $`\sigma _T6.6\times 10^{24}`$ cm2 the Thomson cross-section). This is to be compared with the required energy $`(R^{(\mathrm{opt})}/3R_g)W_{\mathrm{em}}^{(\mathrm{opt})}`$, which yields

$$\begin{array}{c}R_g>\left(\frac{2}{81\pi }\sigma _T\frac{1+\eta ^2}{\eta ^2}\frac{q^2\epsilon _{\mathrm{max}}^8}{\left(mc^2\right)^8m_pc^2}\right)^{1/3}\hfill \\ 3\times 10^{15}Z^{2/3}A^{8/3}\text{cm}.\hfill \end{array}$$
(15)
If the accelerated particles are protons, then only super-massive black holes with $`M>10^{10}M_{}`$ can meet the above requirement. On the other hand, a lighter hole of $`M>2\times 10^6M_{}`$ can accelerate iron nuclei, and, since the size of such a hole is close to the optimal one or exceeds it ($`R_gR^{(\mathrm{opt})}8\times 10^{12}`$ cm), this would not be at odds with the observed $`\gamma `$-ray background.

Active Galactic Nuclei (AGNs). This phenomenon is due to ejection of a relativistic jet by a massive black hole interacting with an accretion disk around it. Inner jets in AGNs have typically Lorentz factors between 5 and 10. The estimate of the magnetic field depend on the radiation model. For leptonic models which imply high-energy $`\gamma `$-ray production through inverse Compton scattering, the magnetic field is of $`0.1`$ G (e.g., AGNsmall ). However, in hadronic models, in particular in the extreme proton-synchrotron model, the field strength exceeds $`100`$ G AGNlarge . Furthermore, the size of acceleration region, $`10^{14}10^{16}`$ cm as inferred from variability timescale, is close to the optimal value for $`10^{20}`$ eV protons. This eliminates problems with $`\gamma `$-ray background and minimizes energy requirements that constrain the Poynting-flux luminosity of the source in the case of a relativistic outflow. This luminosity ($`\mathrm{\Gamma }W_{\mathrm{em}}c/R^{}`$) must be higher than $`10^{45}`$ erg/s (or a factor of $`(\theta \mathrm{\Gamma })^2`$ more if the opening angle of the jet, $`\theta `$, is larger than $`\mathrm{\Gamma }^1`$), which is reasonable for a typical AGN.

Large accelerators. There are three kinds of very large objects that were proposed to produce EHECRs: knots (or hot spots) in large-scale AGN jets, radio lobes (shocked regions), and galaxy clusters (e.g., Large ; H84 ; Bier ). Despite their sizes are far from the optimum: $`1`$ kpc, $`100`$ kpc, and $`1`$ Mpc respectively, the energy of a magnetic field is enough to bind $`10^{20}`$ eV protons. The limits on these objects are better formulated in terms of the requirements to the acceleration efficiency $`\eta `$, and they are likely to arise from particle diffusion. In all cases there is about 1000-fold excess of the energy stored in magnetic field with respect to that shown in Fig. 1. Therefore, if particles escape on the Bohm-diffusion timescale $`3R^2/(r_Hc)`$, a rather high acceleration efficiency $`\eta 0.03`$ is necessary. Even if the diffusion is strongly suppressed (e.g., special configuration of magnetic field lines), there are limits on $`\eta `$. It must be larger than $`(0.31)\times 10^2`$ in radio lobes and galaxy clusters, otherwise protons are accelerated so slowly that they start to lose energy interacting with the cosmic microwave background, which eventually stops acceleration. In hot spots, with their stronger magnetic field up to $`B10^3`$ G, $`\eta `$ can be as small as few $`\times 10^5`$, in which case their size becomes slightly smaller than the optimal one leading to a further increase in energy requirements.

Gamma-Ray Bursts (GRBs). These explosion-like events have been argued to be suitable places for acceleration of EHECRs (e.g., Wax ; Viet ), owing to their compactness (variability timescale less than a second), high apparent ($`4\pi `$) luminosity of the order of $`10^{52}`$ erg/s, and ultrarelativistic bulk flows with the Lorentz factors $`\mathrm{\Gamma }1001000`$, which are the key ingredients of any GRB model (see, e.g., GRBrev  for a review). Acceleration may take place either in internal shocks formed within the outflow (at a distance $`10^{10}10^{13}`$ cm from the central engine) by colliding shells of different Lorentz factors or in the external shock (at a distance $`10^{16}10^{17}`$ cm) at the interface between ejected material and interstellar gas. The distance from which the accelerated protons can escape is, however, limited at least by their interaction with the intense thermal photon field left after adiabatic expansion of relativistic fireball. The escape distance is (in the shock comoving frame)

$$R_{\mathrm{esc}}^{}0.08\frac{f\sigma _{p\gamma }L}{\pi \mathrm{\Gamma }^3cT}\frac{2\times 10^{17}\text{cm}}{\mathrm{\Gamma }^3}\frac{L}{10^{52}\text{erg/s}}\frac{3\text{MeV}}{T},$$
(16)
where $`L`$ is the GRB $`4\pi `$-luminosity and $`T`$ the temperature at the base of the fireball, $`T3`$ MeV. The Eqs. (11b) and (16) scale with the same power of $`\mathrm{\Gamma }`$, so that the ratio between $`R_{\mathrm{esc}}^{}`$ and $`R^{(\mathrm{opt})}`$ depends practically only on GRB luminosity, and the escape distance appears to be larger than the optimal size, unless a GRB has unusually low luminosity. Therefore, radiative losses are negligible and the maximal particle energy is $`\epsilon _{\mathrm{max}}^{}=\eta qB^{}c(R^{}/c)`$. With the magnetic field related to the GRB $`4\pi `$-luminosity via the parameter $`ϵ_m`$, which is the fraction of energy carried in the form of magnetic field, this gives

$$\epsilon _{\mathrm{max}}=\mathrm{\Gamma }\epsilon _{\mathrm{max}}^{}=\frac{\eta ϵ_m^{1/2}}{\mathrm{\Gamma }}q\sqrt{\frac{2L}{c}}\frac{2.5\times 10^{23}\text{eV}}{\mathrm{\Gamma }}.$$
(17)
Note that the maximal energy $`\epsilon _{\mathrm{max}}`$ decreases with increasing bulk Lorentz factor $`\mathrm{\Gamma }`$. This can be understood following the same reasoning as at the end of the Sect. II: changing the Lorentz factor allows one to gain only from a more favorable ratio $`R^{(\mathrm{opt})}/R^{}`$, but for GRBs the escape distance is larger than the optimal size, so that the increase in $`\mathrm{\Gamma }`$ makes things worse.

Since explanation of some GRB properties, such as rapid variability and the presence of emission in the $`10`$ GeV domain Hurley , requires $`\mathrm{\Gamma }>100`$, we conclude that GRBs might be capable of accelerating protons well above $`10^{20}`$ eV. However, to match the observed cosmic ray flux at least 50 per cent of the observed GRB energy has to be converted into EHECRs.

## V Conclusion

Radiation losses, either due to synchrotron or curvature radiation, play increasing role in the energy balance of accelerated particles with decreasing size (and accordingly increasing magnetic field) of the accelerator. In sufficiently compact objects the radiative losses limit maximal particle energy to the values that are smaller than those determined from the Hillas criterion. Thus, there are optimal size and magnetic field strength that minimize the requirements to the internal energy of the accelerator. For $`10^{20}`$ eV protons and acceleration efficiency $`\eta =1`$, the optimal size is $`R^{(\mathrm{opt})}10^{17}`$ cm and the optimal magnetic field is $`B^{(\mathrm{opt})}3`$ G. This corresponds to the total electromagnetic energy in the acceleration region $`W_{\mathrm{em}}^{(\mathrm{opt})}3\times 10^{51}`$ erg. Relativistic bulk motion reduces the energy requirements and the optimal comoving size, while increasing the optimal comoving magnetic field. The Lorentz factor $`\mathrm{\Gamma }=10`$ leads to $`R^{(\mathrm{opt})}10^{14}`$ cm and $`B^{(\mathrm{opt})}300`$ G, corresponding to the total electromagnetic energy $`W_{\mathrm{em}}^{(\mathrm{opt})}3\times 10^{47}`$ erg in the region of size $`R^{(\mathrm{opt})}`$.

Efficient acceleration in optimally-sized objects leads to significant radiative losses, which in this case occur in a regime transitional between the synchrotron and the curvature radiation. The emission is concentrated at photon energies of $`100300`$ GeV (protons) and $`50150`$ MeV (electrons) times the bulk Lorentz factor. This suggests the use of new generation of air-Cherenkov telescopes and GLAST, respectively, for a search for efficient accelerators.

The derived constraints seem to firmly exclude any known population of neutron stars as acceleration sites of EHECRs. Black holes also appear to be not able to accelerate protons to $`10^{20}`$ eV except for the supermassive ones with masses $`M>10^{10}M_{}`$, although acceleration of heavy nuclei by lighter holes is still a possibility. Other astrophysical objects considered above are able, in principle, to accelerate EHECRs, although the required efficiency $`\eta `$ is always comparable to unity. In the standard treatment of the problem, such an efficiency is difficult to achieve even in ultrarelativistic shocks rel-eta . In addition, energy losses due to photomeson reactions may be a serious problem for radiation-rich environments of supermassive black holes and for GRB sources.

The results obtained in this paper are the unavoidable limits, set by the laws of classical electrodynamics. Taking into account the effects of diffusion, the energy losses via photomeson process, the energy content of particles and waves in cosmic plasma, etc., should lead to still stronger limits. However, they cannot, of course, constrain the top-down models for the origin of EHECRs (e.g., top-down ), which make use of quantum effects and non-electromagnetic interactions.

In view of the development of advanced experiments which aim at detection of cosmic rays (CRs) above $`10^{21}`$ eV, it is worth noting what is the highest energy of EHECRs still consistent with the acceleration scenario of their origin. Plotting diagrams similar to Figs. 1 and 2, but for different values of $`\epsilon _{\mathrm{max}}`$, one can see that in GRBs and AGN inner jets protons can in principle be accelerated up to $`10^{21}`$ eV. At such energies, acceleration is likely to be radiative-loss limited in both cases. Above $`10^{21}`$ eV GRBs fail because of their insufficient duration and AGN inner jets – because of insufficient Poynting-flux luminosity. Hot spots, radio lobes and galaxy clusters can still work to $`(35)\times 10^{21}`$ eV under very speculative assumption that the magnetic field is ordered on all scales and the acceleration efficiency is $`\eta 1`$. In this case, acceleration is escape-limited. At the energies $`10^{22}`$ eV the CR primaries have to be heavy nuclei. In all the sources listed above the heavy nuclei are accelerated in the escape-limited regime, so that the attainable energy is roughly $`Z`$ times more than for protons. However, the nuclei of such energy are fragmented through interaction with the microwave background photons after traveling a distance of less than $`1`$ Mpc, that means they must be produced within the local group of galaxies, and GRBs is the only possibility to do this. Although the nuclei are easily fragmented in radiation-reach environments of GRBs, we have to conclude that formally the primaries with energy up to $`(23)\times 10^{22}`$ eV can be produced within the framework of acceleration scenario.


text:
# Evolving networks with distance preferences

## I Introduction

Graphs can be considered as substrata of dynamic networks, and so, several types of graph models have been proposed for capturing the properties of specific networks str ; bara ; doro . In particular, evolving networks can be modelled through growing graphs, i.e. graphs to which continuously new nodes (vertices) and new links (edges) are added. While regular graphs, i.e. ones where each node has the same connectivity pattern and where consequently interactions are local in nature and progress in a slow and orderly fashion from neighbor to neighbor, can exhibit subtle combinatorial patterns, for a realistic network model typically a certain amount of irregularity or randomness is needed. The prototypes here are the random graphs introduced by Erdös and Rényi where the connections between the nodes are chosen completely randomly bol . These exhibit quite interesting properties, but often real networks are not entirely random in this sense, but show some kind of regularity, not directly in their connectivity pattern, but with respect to some other variable or order parameter. Such a parameter can be a clustering coefficient, the average or maximal distance between nodes in the network (as measured by the minimal number of links separating them), the distribution of the number of links between the nodes, the correlation of such properties between neighboring nodes (i.e. those connected by a link of distance 1), or the first eigenvalue of the graph Laplacian which is relevant for synchronization properties throughout the network of dynamic activities at the individual nodes ws ; bara ; far ; go ; jo . Models have been proposed that capture some of these aspects. The small world networks introduced by Watts and Strogatz ws  are constructed from regular graphs by creating additional random links between nodes, with or without deleting some of the existing ones. Once a certain number of such new links has been introduced in proportion to the number of regular ones, distances in the graph get dramatically shortened, and, consequently, activity can spread quite rapidly from a localized source through the entire network. Another distinct feature of this model is that there is clustering which is absent in random models. Empirical evidence is available for the occurrence of clustering in real networks new1 . Another interesting model is the one of a scale free network as introduced by Barabási and Albert ba ; bara . This is a graph where new nodes are added and form a fixed number of links with the existing nodes not completely at random, but with a preference towards those nodes that already have more connections than other ones. More precisely, the probability with which existing nodes receives a link from a new node is proportional to the number of links it already possesses. The characteristic feature of the emerging graph here is that the number of nodes with a given number of links does not decrease exponentially as a function of the latter as for example in random graphs, but follows a power law – the reason why such a graph is called scale free. Such models can provide valuable insights into existing real networks, for example into patterns of social relations or spreading of diseases in the small world model, or the connection patterns of internet sites or flight connections between airports in the scale free model.

It is then a natural question whether there exists an encompassing scheme that on one hand can put these specific models into a more general perspective and that on the other hand can offer systematic tools for analyzing the dependencies among the various network features listed above. Ideally, these features should depend in an analyzable manner on certain parameters of the network construction, and so their interdependencies could then be studied in terms of relations between the parameters involved.

We attempt here to take a step in this direction by proposing a general scheme for constructing evolving networks. Our model is characterized by a distance preference function. This function specifies the probability in terms of the distance with which an existing node in the network receives a new link from a newly created node that already has formed one random link so as to attach it to the network and to define its distances to the other nodes. The number of links each node is allowed to make can be either fixed – as in our simulation results below – or also follow some random distribution. So, for example, we can stipulate that the shortest distances are always preferred. Thus, a node that is allowed to form a new link does so preferably to another node of distance 2, i.e. to a direct neighbor of a node that it is already attached to. This might constitute a useful model for the formation of social relationships (you want to become a friend of the friends of your present friends as the easiest or safest means of forming new relationships). Conversely, we might also stipulate that always the most distant nodes are the preferred recipients of new links. Obviously, one then expects that the resulting network has a quite short average distance between any two nodes, as in the small world and scale free models. In fact, however, our simulations demonstrate that directly selecting distances is not as efficient for reducing the average distance in the network as creating some highly connected nodes through which many shortest connections can go, as in the scale free model. More interestingly perhaps, one may even expect a certain tendency towards the scale free type when shortest distances are preferred. Namely, a node that is highly connected then has a greater chance of receiving a new link than a less well connected one, because the former has a greater chance of being a direct neighbor of another node that has received a previous link from a new node that is attaching itself to the network. Thus, we see the principle that the rich get richer that is characteristic for scale free networks also at work here, although in an indirect and somewhat mediated form. A conceptual advantage of this construction over the scale free one might be that here, for each link, we only need to evaluate local information, namely check those sites in its vicinity. More precisely, if we exclusively select sites of distance two as recipients of new links, then we only have to list all the neighbors of the present neighbors of the link forming node at each step. In contrast to this, for the scale free model, the complete connectivity pattern of any potential recipient anywhere in the network has to be evaluated. In general, in our scheme, whether we give preference to short distances or not, what is crucial for the decision about a new link is not an absolute property of the candidate as in the scale free model, but rather its relation, as expressed by the distance, to the link forming node. This may capture a property that is relevant in some applications.

On the other hand, the scheme where short distances are preferred should lead to more pronounced local clustering effects and larger average distances in the network than the scale free construction model. In this way, we can check that certain network properties are independent of or at least not strongly related to each other.

Of course, our scheme also includes the possibility that all distances are equally preferred. This should generate properties similar to a random network, although the construction is not entirely identical, because for a random graph, all nodes are considered equal, whereas here, only those of the same distance to the node forming links have equal recipient probabilities, because the distances need not be evenly distributed among the nodes.

We could also easily supplement our construction scheme by a rule for the deletion of links and/or nodes according to some criterion to be specified, as a means to stabilize the size of our network. This would allow a comparison of our model with other ones for evolving networks of given size range. Here, however, we do not pursue this aspect systematically.

## II Network construction

We start with a small network having $`m_0`$ nodes and then let it grow according to the following scheme. We fix a number $`m`$ as the number of connections each new node is allowed to establish to other nodes existing in the network; in principle, this number could also be randomly chosen from some distribution instead of being fixed, but, for simplicity, in our simulations, we only work with a fixed $`m=m_0`$, as this will probably not dramatically affect the resulting network properties. The crucial part of our scheme is the specification of a probability distribution $`p(d)`$ for the preferred distance to a node with which a new link is established. So, when a new node $`x_n`$ comes in, it is first allowed to make one connection with a randomly chosen node in the network, in order to attach it to the network. (We could also change this rule and let the first connection prefer well connected recipient nodes, as in the scale free model, but in the present paper, we do not perform numerical simulations for that rule.) This leaves us with $`m1`$ further links that it is allowed to establish. For the formation of any such link, we consider a node $`x`$ in the network and select it as the recipient of the new link with a probability given by $`p(d(x_n,x))`$. Of course, the formation of any new link changes the distances in the network and the creation of further links, until the allotted number $`m`$ of them has been formed from $`x_n`$, then is governed by the new distance pattern. Once $`x_n`$ is connected according to this scheme, we create a new node $`x_{n+1}`$ and repeat the procedure.

The distance preference function $`p(d)`$ encodes all the features of our construction. An important case is where this function is in fact deterministic, namely where only nodes of distance 2 from $`x_n`$ are allowed as link recipients, i.e. the ones that have the smallest possible distance from it (we are not allowing multiple links, and so no further link can be attached to a node at distance 1). Another deterministic choice of $`p(d)`$ would be to allow only recipients of maximal distance from $`x_n`$. This obviously makes the scheme computationally much more expensive than the exclusive selection of nodes at distance 2. More generally, we are interested in distance preference functions $`p(d)`$ that are decreasing functions of $`d`$, i.e. where short distances are preferred over large ones, but the latter can still be selected with a positive probability.

In our simulations as described in the Table 1, we consider the cases where the number of links that each new node is allowed to form is $`m=2,3,4,`$ and 5. We let the network grow until its size was 30,000 nodes when we evaluated the various parameters. We considered three different versions of the probability for the distances. In Model 1, we exclusively selected links to nodes of distance 2, i.e. we always formed triangles. In Model 2, we let the probability be proportional to the inverse distance. Thus, there was a (slight) preference for shorter distances over larger ones. In Model 3, in contrast to this, we let the preference function be proportional to the distance itself (scaled with the maximal distance in the network). Thus, there is a preference for larger over shorter distance. Our comparison models are the growing random graph model where all $`m`$ links are randomly connected (Model 4) and the scale free or real world model (Model 5).

In Table 1 we give the first eigenvalue $`\lambda _1`$, the clustering coefficient $`C`$, the mean path length $`L`$ and the second moment of degrees $`<k^2>`$, for different $`m`$ values, for Models 1 to 5. The discussion below will employ the simulation results for $`m=5`$; as one can see from the table, the results for $`m=3,4`$ are qualitatively similar but $`m=2`$ is slightly different. The table gives the averages over 10 simulations each; the standard deviations are quite small.

## III First eigenvalue

Spectral properties of small world, scale free and random graph models have been discussed in far ; go . The first (nonzero) eigenvalue of the graph Laplacian is the crucial parameter for the synchronization properties of activities at the network sites as systematically investigated in our previous work jo ; see also rang . We naturally assume here that the graph $`\mathrm{\Gamma }`$ under consideration is connected, as are the graphs resulting from our constructive scheme. Moreover they are symmetric because we consider undirected links. We label the nodes of $`\mathrm{\Gamma }`$ as $`x_1,x_2,\mathrm{},x_n`$, and we let $`k_i`$ denote the connectivity, i.e. the number of neighbors of the node $`x_i`$. The first eigenvalue is then given by

$$\lambda _1=\underset{u:\mathrm{\Gamma },{\scriptscriptstyle k_iu(x_i)}=0}{inf}\frac{_{x_ix_j}(u(x_i)u(x_j))^2}{k_iu(x_i)^2},$$
(1)
where $`x_ix_j`$ denotes that they are neighbors. We can now provide the following heuristic argument how the creation of a new link in the network affects $`\lambda _1`$ depending on the distance $`d(x,y)`$ between the two nodes $`x,y`$ before the link between them is formed. Namely, for any function $`u`$ as evaluated for the infimum in (1), the new link only creates an additional summand $`(u(x)u(y))^2`$ in the numerator while the denominator is left unchanged. As the difference in $`u`$ between neighbors is minimized for a first eigenfunction, the expected squared difference $`(u(x)u(y))^2`$ should be an increasing function of the distance between $`x`$ and $`y`$. Therefore, the value of a typical candidate function $`u`$ for the infimum in (1) should increase as a result of the new link in a manner that is positively correlated with the distance $`d(x,y)`$. Thus, if our scheme prefers larger distances the first eigenvalue should get larger than when we select short distances for new links. Of course, this fits well together with the fact that on one hand, a larger $`\lambda _1`$ facilitates synchronization across the network, and on the other hand, connecting nodes that had a large distance should have the effect of a more pronounced decrease of the average distance which in turn facilitates synchronization as well.

Our simulations (as described in the Table 1 ) yield that the first eigenvalue for Model 1 is .005 which is quite close to the value for a regular network. Thus, synchronization is quite difficult in such a network although the average or maximal distance in the network are quite low (as described below) and the degree distribution of the nodes is quite similar to the scale free case. In all the other models, $`\lambda _1`$ is substantially larger, namely around .39 for Models 2-4 and .41 for Model 5. It might be of some interest that it appears to be about the same or perhaps even slightly smaller in Model 2, where shorter distances are preferred, than in the random Model 4, which in turn has a smaller value than Model 3 with the preference for larger distances. Thus, the scale free model is the most easily synchronizable of the five, a not always desirable property.

## IV Clustering

If our distance preference is for the shortest possible distance, namely 2, then the emerging graph will contain many triangles, i.e. triples of nodes of mutual distance 1. As a consequence, we expect that the graph contains highly connected subclusters.

Also, since the creation of any new link increases the first eigenvalue, it has been suggested by Eckmann and Moses eck  to employ the number of triangles for defining some notion of curvature of a graph. This is based on an analogy with Riemannian geometry where the so-called Ricci curvature yields a lower bound for the eigenvalue of the Laplace-Beltrami operator (the Riemannian version of the Laplacian). In other words, the larger the curvature, the higher the expected value of the first eigenvalue. As our preceding heuristic analysis of the first eigenvalue of the graph Laplacian shows, however, there is a problem with the analogy between the number of triangles and the curvature. Namely, if we add a link to a given graph, then the expected increase in the eigenvalue is the higher, the larger the original distance between the two linked nodes was. In other words, when we select the new link so as to form a new triangle, the expected eigenvalue increase is smallest, or, when trying to pursue the analogy with Riemannian geometry, the additional curvature is least.

The clustering coefficient, $`C`$, of the graph is defined as follows new2 ,

$$C=\frac{3\times \text{(number of triangles on the graph)}}{\text{(number of connected triples of vertices)}},$$
(2)
where a “triangle” is a trio of vertices connected to each other and a “connected triple” is a vertex connected to an (unordered) pair of other vertices. For our choice $`m=5`$, for a regular network the value for $`C`$ is 2/3 (as the number of links of each node is constrained, not all the neighbors of a given node can be connected among each other, and so the value is smaller than 1 in any case). In our Model 1, the value 0.20 is quite high, as to be expected, whereas in all other Models, it is dramatically smaller. In fact, for Model 3 as well as for the random Model 4, it is even smaller than for the scale free Model 5. In particular, the difference between the Models 1 and 2 is striking here.

## V Distances

As already explained, the resulting average or maximal distance in our network should be smaller when large distances are preferred for the establishment of new links. However, this is not so easy to support through numerical simulations, as in any case, independently of the preference function adopted, our networks, like the small world and scale free ones, exhibit rather small maximal distances, say around the order of 4 or 5 for networks with ten or twenty thousand sites, and so the difference resulting from the preference function cannot be very pronounced.

There is one observation that can be made here, however. Namely, the direct preference for forming links to nodes at largest distance is not as efficient in reducing the average or maximal distance in the network as the more indirect scheme of preferential attachment to highly connected nodes employed in the Barabási-Albert model. This demonstrates the virtue of the latter model. In fact, the average distance $`L`$ between all possible pairs of nodes is smallest for that model, namely 4.06, around 4.6 for Models 2-4, and about 5.2 for Model 1. Not surprisingly, a preference for short connections leads to a larger average distance although the effect is by no means as pronounced as one might naively expect. It is surprising, however, that $`L`$ is slightly larger for Model 3 where large distances are preferred than for the random Model 4, and slightly smaller for Model 2 with its preference for shorter distances.

## VI Degree distribution

One of the distinguishing features of the scale free or real world model (Model 5) is that the distribution of the degrees of the nodes decays like a power law in contrast to the exponential of, for example, the random graph model. In Figure (1.a –1.e) we give the plots for degree distribution, $`P(k)`$ for models 1–5, respectively, with $`m=4`$. We find that in our Model 1, where exclusively short connections are selected once a node is anchored in the network, the degree distribution likewise follows a power laws, at least over most of its regime. (For $`m=3`$, we get a power law distribution only for some part of the distribution while the end decays exponentially.) Thus, our mechanism is capable of producing a network that exhibits a power law distribution of the degrees but that differs from the scale free model with respect to a number of distinctive other parameters, like first eigenvalue and synchronizability, clustering, average distance, etc. In particular, this feature is independent of those other features.

Models 2 and 3 show an exponential distribution as in the random model (Model 4). We also find that the distribution of the neighbor degrees (i.e. the sum of the degrees of all the neighbors of a given node, $`P(kk)`$) also partly follows a power law in our simulations for Models 1 and 5. In Fig. 2 we plot that for Model 1 with $`m=4`$.

## VII Correlations

We may ask whether our scheme leads to strong correlations between neighboring sites in the network, with regard to their connectivity. One possible source of such a correlation in connectivity could be a correlation in age. Namely, older nodes in the network have had more chances than younger ones of receiving a random connection from a new node, and so, the connectivity should be positively correlated with the age of a node. However, there is no direct reason why neighboring nodes should exhibit a pronounced age correlation.

Another line of reasoning can go as follows: If $`x_1`$ is a neighbor of a site $`x_2`$ of connectivity $`l`$, then if distance 2 is selected by our preference function, then $`x_2`$ has an $`l`$-fold chance of receiving the second connection that a new node $`x_n`$ is making, but the chances of $`x_1`$ to benefit from this and receive the third connection that $`x_n`$ is making is proportional to $`1/l`$ as it is facing the competition of the $`l1`$ other neighbors of $`x_2`$. Thus, the factors cancel, and here, we do not get an advantage for a node from being a neighbor of a well connected node. Of course, this heuristic argument does not take the triangle pattern in the network into account. We calculated the average of the square of the degrees of the nodes (second moment), $`<k^2>`$. The result is given in the last column of the table. The value of this parameter is around 250 for models 1 and 5 while for models 2, 3 and 4 it is almost half of that value.

## VIII Comparison with other recent network constructions

Dorogovtsev et al. dms ; doro  introduced a model which is similar to the special case of our Model 1, where each new node forms only two links and triangles are exclusively selected. They attach new nodes to the network with links to the two ends of some randomly chosen link already present in the network. This scheme depends on the distribution of links whereas the Model 1 depends on distribution of nodes, though in both cases triangles are formed.

Vazquez va  studied a network where the growth depends on the knowledge obtained by ’walking’ on it. It is a directed graph model unlike our model. New links are formed with a probability $`p`$ to a neighbour of a randomly linked node from the new node and this process is recursively continued. New nodes are added when there is no new link to form. Beyond a critical $`p`$-value it produces scale free network. Here when $`p=1`$, neighbours are preferred as in our Model 1 but the process continues recursively to produce a lot more links of longer distances.

Jin et al. jin  introduced a model with fixed number of vertices where the probability of formation of new links between two nodes depends preferentially on the number of mutual neighbours. There is a cutoff on the number of neigbours possible and a possibility for node removal. This model gives graphs with high clustering coefficient but there is no scale free degree distribution.

Holme and Kim kim  introduced a model that in some respects is similar to our Model 1. They let the first connection of a new node form according to preferential attachment as in the scale-free model and then introduce subsequent links that either form triangles or constitute once more preferential attachments, according to some random preference. The resulting network is again scale free. Their main result is that in a scale free network, the clustering coefficient can take different values (according to the strength of the triangle preference).

Klemm and Eguíluz kle  consider a growing network model based on the scale free paradigm, with the distinctive feature that older nodes become inactive at the same rate that new ones are introduced. This is interpreted as a finite memory effect, in the sense that older contributions tend to be forgotten when they are not frequently enough employed. This results in networks that are even more highly clustered than regular ones.

Davidsen et al. born  consider a network that rewires itself through triangle formation. Nodes together with all their links are randomly removed and replaced by new ones with one random link. The resulting network again is highly clustered, has small average distance, and can be tuned towards a scale free behavior.

## IX Conclusion and discussion

We have introduced a model for evolving networks where each new node, once it is (randomly) anchored to the network, forms further links according to some distance preference function, and we have compared simulation results for the evolved networks with those for two main types previously considered, namely the random graph model and the scale free or real world model of Barabási-Albert. We found that when always the shortest possible distances are selected for the recipients of new links, we get a highly clustered network which is difficult to synchronize, although it still has a relatively small average distance between nodes. It also exhibits a power law type behavior for the distribution of the degrees of the nodes comparable to the scale free model, although the underlying network forming mechanism is different, and, in particular, there is no explicit preference for highly connected nodes which is considered as the main reason for the power law behavior in the scale free model.

It has been shown that linear preferential attachment is a necessary condition for a growing power law network krap . To check this in our model, we calculated the attachment rate, $`\mathrm{\Pi }(k)`$, as a function of the degree $`k`$. To calculate this we used the method described in je . The attachment rate is numerically fitted with a power law in $`k`$ and we obtained the power equal to 1.0 for Model 1 (for $`m=5`$) and Model 5 and 0.0 for Models 2, 3 and 4. (In Model 1 for smaller values of $`m`$, this exponent is less than 1.) This indicates that there is preferential linear attachment in our Model 1 as in the case of BA model though we don’t explicity introduce that in our model. Surprisingly for Model 2, though it is similar to Model 1, the attachment rate is independent of the degree as indicated by the zero exponent of $`k`$. This explains why the degree distribution is similar to that of a random one. Even the small probability of attaching to second and higher order neighbours in Model 2 produces deviation from linear preferential attachment rate. The number of second and higher order neighbours are not linearly proportional to the number of first neighbours of a vertex in these models.

As the other network parameters are different from the scale free model, this shows that this feature is independent of clustering or synchronizability properties. For other distance preference functions, we found network parameters that were roughly comparable with the ones for a random graph network, and in fact regardless of whether our preference was proportional or inversely proportional the distance between the link forming node and the potential recipient.

###### Acknowledgements.

We thank the referee, Mark Newman, and several readers of our paper for useful comments.


text:
# Brownian dynamics approach to interacting magnetic moments

## Abstract

The question how to introduce thermal fluctuations in the equation of motion of a magnetic system is addressed. Using the approach of the fluctuation-dissipation theorem we calculate the properties of the noise for both, the fluctuating field and fluctuating torque (force) representation. In contrast to earlier calculations we consider the general case of a system of interacting magnetic moments without the assumption of axial symmetry. We show that the interactions do not result in any correlations of thermal fluctuations in the field representation and that the same widely used formula can be used in the most general case. We further prove that close to the equilibrium where the fluctuation-dissipation theorem is valid, both, field and torque (force) representations coincide, being different far away from it.

The problem of a correct introduction of temperature in the equation of motion of a magnetic system has gained much importance as a result of technological requirements of magnetic recording industry Weller ; Charap ; Bertram . This is associated with the need to perform calculations of magnetization dynamics at finite temperatures. Open problems include fast magnetization switching, thermal stability and magnetic viscosity, among others. The correct solution of the problem is still far from being understood. The main difference between the magnetic problem and the standard molecular dynamics approach is that the magnetic moment dynamics is governed by the Landau-Lifshitz equation which includes the precession of a magnetic moment around its internal field direction. It comprises coupled first-order equations for the magnetization components, and the requirement of conservation of the magnetization amplitude. As a consequence, no analogue of mass and kinetic energy exist in the system, thus making it impossible to introduce the temperature through this mechanism.

Consequently, the temperature is introduced through small deviations from the equilibrium configuration. Therefore, strictly speaking, this approach is only valid when these deviations are small and it cannot be used for fast magnetization switching.

Let us briefly summarize the original approach from W. Brown Brown ; Brown1 . The underlying equation of motion is the Landau-Lifshitz-Gilbert equation which can be written in the form

$$\frac{d\stackrel{}{\text{M}_i}}{d\tau }=\stackrel{}{\text{M}_i}\times \stackrel{}{\text{H}_i}\alpha \stackrel{}{\text{M}_i}\times [\stackrel{}{\text{M}_i}\times \stackrel{}{\text{H}_i}],$$
(1)
where

$$\tau =\frac{\gamma _0H_k}{M_s(1+\alpha ^2)}t,\stackrel{}{\text{H}}=\frac{\delta E^{}}{\delta \stackrel{}{\text{M}}}$$
(2)
$`\gamma _0`$ is the gyromagnetic ratio and $`\alpha `$ is the damping constant. The magnetic moment $`\stackrel{}{\text{M}}`$ is normalized to the saturation value $`M_s`$, and the internal field $`\stackrel{}{\text{H}}`$ is normalized to the anisotropy field $`H_k=2K/M_s`$. The energy $`E^{}=E/2KV`$, where $`K`$ is the anisotropy value and $`V`$ is the particle volume, contains all the necessary energy contributions: anisotropy, exchange, magnetostatic and Zeeman.

W. Brown proposed the inclusion of thermal fluctuations via a random field, added to the internal field, Eq. 2. For the calculation of the properties of the random field he outlined two methods: (i) based on the fluctuation-dissipation theorem (see also Berkov ) and (ii) by imposing the condition that the equilibrium solution of the correspondent Fokker-Planck equation is the Boltzmann distribution (see also Garcia-Palacios ). As a result of both the thermal field statistical properties are given by

$$\xi _i=0,\xi _i(0)\xi _j(\tau )=\frac{\alpha k_BT}{KV(1+\alpha ^2)}\delta _{ij}\delta (\tau ),$$
(3)
where $`i,j`$ denote Cartesian components $`x,y,z`$. Different approaches based, for example, on the Landau-Lifshitz rather than on the Landau-Lifshitz-Gilbert equation were also introduced Garcia-Palacios ; Moro .

However, the properties of the thermal noise, Eqs. 3, were derived only for one isolated particle. Moreover, Brown considered in his paper Brown  only the simplest axially symmetric case. Nevertheless, in the past the formulas above provided the basis for practically every numerical method Charap ; Lyberatos ; Hannay ; Nowak  for the computation of magnetization dynamics taking into account thermal fluctuations. But the investigated magnetic systems usually comprise interacting particles Braun ; Nakatani ; Boerner ; Zhang ; Hinzke ; Schrefl  due to magnetostatic and/or exchange couplings. For that case the thermal field may be expected to be influenced by correlations between different particles Hannay . Hence, it is necessary to generalize Brown’s result to the case of interacting magnetic moments. To the best of our knowledge, this has never been done before.

In what follows we start with the Brownian dynamics approach (see LandauLifshitz ) which was originally applied to magnetic systems by A. Lyberatos et al. Berkov ; Hannay . However, we consider the general case of an interacting system with a non-axially-symmetric potential. Following the standard approach, we introduce the temperature into the motion of the Brownian particles (i. e. the magnetic moments) as a result of the fluctuation-dissipation theorem. Consequently, this approach is only valid when small deviations from equilibrium are considered.

The general Langevin equation of motion is written in the form

$$\frac{dx_i}{dt}=\underset{j}{}\gamma _{ij}X_j+f_i,X_j=\frac{S}{x_j},$$
(4)
where the $`\gamma _{ij}`$ are the so-called kinetic coefficients, $`X_j`$ are variables which are thermodynamically conjugate to $`x_j`$, and $`S`$ is the entropy of the magnetic system. For a closed system in an external medium,

$$X_j=\frac{1}{k_BT}\frac{E}{x_j}.$$
(5)
In Eq. 4, $`f_i`$ is a random force representing thermal fluctuations in the system having the properties

$$f_i(t)=0\text{and}f_i(0)f_j(t)=\mu _{ij}\delta (t)$$
(6)
where

$$\mu _{ij}=\gamma _{ij}+\gamma _{ji}.$$
(7)
A linear equation of motion of the form

$$\frac{dx_i}{dt}=\underset{j}{}L_{ij}x_j,$$
(8)
with the associated energy

$$E=E_0+\frac{1}{2}\underset{i,j}{}A_{ij}x_ix_j,$$
(9)
can be rewritten as

$$\frac{dx_i}{dt}=\underset{j}{}L_{ij}x_j=\underset{j}{}\gamma _{ij}\frac{1}{k_BT}\underset{k}{}A_{kj}x_k,$$
(10)
so that the matrix $`L_{ik}`$ is related to the kinetic coefficients $`\gamma _{ik}`$ in the following way Berkov :

$$L_{ik}=\frac{1}{k_BT}\underset{j}{}\gamma _{ij}A_{kj}$$
(11)
In micromagnetics the motion of a magnetic moment $`𝐌`$ is governed by the deterministic LLG equation (Eq. 1). For the equilibrium state of the system Brown’s condition

$$\stackrel{}{\text{M}_i^0}\times \stackrel{}{\text{H}_i^0}=0$$
(12)
must be satisfied, implying that here $`\stackrel{}{\text{H}_i^0}`$ and $`\stackrel{}{\text{M}_i^0}`$ are parallel. Close to equilibrium, the LLG equation can be linearized using small deviations

$$\stackrel{}{\text{m}_i}=\stackrel{}{\text{M}_i}\stackrel{}{\text{M}_i^0},\stackrel{}{\text{h}_i}=\stackrel{}{\text{H}_i}\stackrel{}{\text{H}_i^0}$$
(13)
from their equilibrium values, yielding

$$\frac{dm_i}{dt}=\underset{j=1}{\overset{3N}{}}L_{ij}m_j.$$
(14)
Here, the indices $`i,j`$ count the particles sites $`1,\mathrm{},N`$ as well as their $`x,y,z`$ coordinates. The internal fields $`h_j`$ play the role of the variables which are thermodynamically conjugate to $`m_j`$,

$$X_j=\frac{1}{k_BT}\frac{E}{m_j}=\frac{2KV}{k_BT}h_j.$$
(15)
Thus, the LLG equation should be rewritten in the form

$$\frac{dm_i}{dt}=\frac{2KV}{k_BT}\underset{j=1}{\overset{3N}{}}\gamma _{ij}h_j$$
(16)
which is an easier way to calculate the kinetic coefficients than the use of Eq. 11. The representation of the LLG equation in the form of Eq. 4 means that in what follows the thermal fluctuations are introduced as a fluctuating torque (a generalized force rather than a field). Later we will show that in the linear approximation this is equivalent to the standard fluctuating field representation. Alternatively, Eq. 16 could be viewed as a polar representation of the magnetization vector $`m_i^1=\theta _i,m_i^2=\phi _i`$, in this case the conjugate variables are the polar projections of the internal fields $`(h_\theta ,h_\phi )`$ and the fluctuations $`f_i`$ will stand for the random field polar components. This latter approach was used originally by W. Brown Brown .

We continue by writing the energy of the system in the form

$$E^{}=\underset{i}{\overset{N}{}}\left(\stackrel{}{\text{M}_i}\stackrel{}{\text{H}_i}+\frac{\lambda }{2}\text{M}_i^2\right).$$
(17)
where $`\lambda `$ is the Lagrange multiplier. In the zero order approximation one obtains

$$\stackrel{}{\text{M}_i^0}=\frac{1}{\lambda }\stackrel{}{\text{H}_i^0}$$
(18)
which corresponds to Brown’s condition, Eq. 12. The linear approximation leads to the equilibrium condition

$$\stackrel{}{\text{M}_i^0}\stackrel{}{\text{h}_i}\stackrel{}{\text{H}_i^0}\stackrel{}{\text{m}_i}+\lambda \stackrel{}{\text{M}_i^0}\stackrel{}{\text{m}_i}=0,$$
(19)
which leaves us only the quadratic form for the energy expression near the equilibrium,

$$E^{}=E_0\underset{i}{\overset{N}{}}\left(\stackrel{}{\text{m}_i}\stackrel{}{\text{h}_i}\frac{\lambda }{2}\text{m}_i^2\right).$$
(20)
The general expression for magnetic energies is a quadratic form in terms of the magnetization (apart from the Zeeman term which is included in the equilibrium field $`\stackrel{}{\text{H}_i^0}`$ and condition 19). Therefore, it is reasonable to suppose that the total field can be expressed as

$$h_i^\alpha =\underset{j,\beta }{}B_{ij}^{\alpha \beta }m_j^\beta =h_{eff,i}^\alpha \lambda m_i^\alpha $$
(21)
where $`h_{eff,i}^\alpha `$ are the components of the effective field due to the different energy contributions and $`\lambda m_i^\alpha `$ is the field due to the kinematic interaction expressing the constraints (Lagrange multiplier). The value of the Lagrange multiplier is normally found from the equilibrium condition 19. However, its actual value is not necessary for calculations due to the fact that the LLG equation conserves the magnetization length. Latin indices represent the sites of the moments and the Greek ones the magnetization components $`x,y,z`$. In this case the final expression for the energy (Eq. 20) takes the form

$$E^{}=E_0\underset{i,j,\alpha ,\beta }{}(B_{ij}^{\alpha \beta }\delta _{ij}\delta _{\alpha \beta })m_i^\alpha m_j^\beta .$$
(22)
The expressions for the kinetic coefficients could be obtained by using directly the expression 11. In this approach it seems that the final result could also include correlations between different particles Hannay . But this is not the case: the kinetic coefficients can be obtained much easier representing the linearized LLG equation in the form of Eq. 16 yielding

$`\gamma _{ij}^{xx}`$ $`=`$ $`{\displaystyle \frac{\alpha k_BT}{2KV}}\left[(M_i^{0,y})^2+(M_i^{0,z})^2\right]\delta _{ij}`$ (23)
$`\gamma _{ij}^{xy}`$ $`=`$ $`{\displaystyle \frac{k_BT}{2KV}}\left[M_i^{0,z}+\alpha M_i^{0,x}M_i^{0,y}\right]\delta _{ij}`$ (24)
$`\gamma _{ij}^{yx}`$ $`=`$ $`{\displaystyle \frac{k_BT}{2KV}}\left[M_i^{0,z}+\alpha M_i^{0,x}M_i^{0,y}\right]\delta _{ij}`$ (25)
$`\gamma _{ij}^{yy}`$ $`=`$ $`{\displaystyle \frac{\alpha k_BT}{2KV}}\left[(M_i^{0,x})^2+(M_i^{0,z})^2\right]\delta _{ij}.`$ (26)
The other coefficients can be obtained by symmetry. Note, that there are no correlations between different particles. Also, the kinetic coefficients have obviously reversible parts (coming from rotation) and irreversible parts (from damping). The reversible antisymmetric parts do not contribute to the thermal fluctuations after adding the kinetic coefficients to calculate the matrix $`\mu `$ from Eq. 7, yielding

$`\mu _{ij}^{xx}`$ $`=`$ $`{\displaystyle \frac{\alpha k_BT}{KV}}\left[(M_i^{0,y})^2+(M_i^{0,z})^2\right]\delta _{ij}`$ (27)
$`\mu _{ij}^{xy}`$ $`=`$ $`{\displaystyle \frac{\alpha k_BT}{KV}}M_i^{0,x}M_i^{0,y}\delta _{ij}.`$ (28)
Once again, the others can be obtained by symmetry. Note that in a general system of coordinates there are correlations between different magnetization components but no correlations between different particles. However, if we set the local coordinate system such that the $`z`$ axis coincides with the equilibrium magnetization direction, $`M_i^{0,x}=0,M_i^{0,y}=0,M_i^{0,z}=1`$, these correlations disappear and we have the same thermal fluctuations in $`x`$ and $`y`$ directions but no fluctuations in $`z`$ direction,

$$\mu _{ij}^{xx}=\mu _{ij}^{yy}=\frac{\alpha k_BT}{KV}\delta _{ij}\text{and}\mu _{ij}^{zz}=0.$$
(29)
Thus, the torque fluctuations produce effectively correlations and different values of thermal fluctuations in all other systems of coordinates different from the global one, where one of the axes is parallel to the equilibrium magnetization direction and where the equation of motion for this component disappears.

It is customary to introduce thermal fluctuations in the field components (see Lyberatos  and originally W. F. Brown Brown ) instead of the torque fluctuations as derived above. This has its origin in the representation of the LLG equation in a spherical system of coordinates in form of Eq. 4. However, in both of these papers above only the axially symmetric case without interactions was considered. The big difference between these two approaches is the *multiplicative* character of the field noise versus the *additive* noise of the torque. This turns out to be important for larger magnetization deviations. But first we will show that in the global coordinate system both approaches, torque and field, give the same result, as long as the magnetization deviations from the equilibrium are small.

Let us use a decomposition of the field components according to $`H^iH^i+\xi ^i`$, where $`\xi ^i`$ are the components of the fluctuation part of the field. When this is done we obtain the following expansion of the equations of motion,

$`{\displaystyle \frac{dM^i}{d\tau }}`$ $`=`$ $`\epsilon ^{ijk}M^jH^k\alpha H^m\left[M^mM^i\delta ^{mi}\right]`$ (30)
$`\epsilon ^{ijk}M^j\xi ^k\alpha \xi ^m\left[M^mM^i\delta ^{mi}\right]`$
$`=`$ $`A^i(M^n,H^l)+B^{ij}(M^n)\xi ^j.`$
Furthermore, in the global system of coordinates we linearize the magnetization by the decomposition $`M^iM_0^i+m^i`$, where $`m^i`$ are *s*mall fluctuations around the equilibrium values $`M_0^i`$, and apply the constraint condition, $`|\stackrel{}{\text{M}}|=1`$. For simplicity below we drop in the formulas the particle index $`i`$. The components in the specified coordinate system are then

$`{\displaystyle \frac{dm^x}{d\tau }}`$ $`=`$ $`A^x(\stackrel{}{m})\left(m^y+\alpha m^x\right)\xi ^z+f^x,`$ (31)
$`{\displaystyle \frac{dm^y}{d\tau }}`$ $`=`$ $`A^y(\stackrel{}{m})+\left(m^x\alpha m^y\right)\xi ^z+f^y,`$ (32)
$`{\displaystyle \frac{dm^z}{d\tau }}`$ $`=`$ $`\left(m^x+m^y\right)\left(\xi ^x\alpha \xi ^y\right),`$ (33)
where $`\stackrel{}{A}(\stackrel{}{m})`$ stands for the linearized deterministic part of the LLG equation and

$$f^x=\xi ^y+\alpha \xi ^x,f^y=\xi ^y+\alpha \xi ^x.$$
(34)
The constraint condition implies that in a first order approximation it is $`m^z(\tau )=0,\tau `$. This is compatible with Eq. 33 only if the field fluctuations $`\xi ^i`$ can be considered to be small quantities, in which case products of the $`\xi ^i`$ with the $`m^i`$ can be ignored in Eq.(31). These equations suggest that the field fluctuations contribute *additively*. From Eqs. 34 one can also obtain Brown’s formulas for the field fluctuations (Eq. 3).

It is important to note that the last equation is satisfied for any fluctuation field value due to the character of the LLG equation. Thus the $`f^z`$ value (or $`\xi ^z`$) is in this case undefined. In any case the component $`\xi ^z`$ is not efficient since it acts parallel to the magnetization direction. The assumption made in the paper of A. Lyberatos and R. Chantrell Lyberatos  is that the field components are isotropic and that

$$\xi ^x=\xi ^y=\xi ^z.$$
(35)
This assumption in the global system of coordinates (where the fluctuation-dissipation theorem is applied) leads to the remarkable symmetry (35) of the field components in all the systems of coordinates and to the absence of correlations. Furthermore, it is assumed that this property is valid through the magnetization reversal.

For the torque fluctuations the reasonable hypothesis to mimic the field ones would be the assumption that there are never torque (force) fluctuations along the magnetization direction. In this case the correlations between different noise components would appear in all other systems of coordinates different from the global one. While equivalent near the equilibrium, these two approaches will be different far from it. At this point, we would like to restate that the whole theory is valid for small fluctuations around the equilibrium where both approaches coincide.

In conclusion, the application of the Brownian dynamics approach to the motion of a magnetic system shows that interactions do not introduce correlations into thermal fluctuations introduced as both, either a fluctuating torque or a fluctuating field. Correlations may appear between different magnetization components as a result of the conservation of the value of the magnetic moment. The reasonable hypothesis that all the fluctuating field components are equivalent leads to Brown’s well-known formulas for the fluctuating fields values without correlations. This validates all previously done micromagnetic calculations where this kind of assumption was made.

###### Acknowledgements.

OC acknowledges the hospitality and support from Durham University, UK, Duisburg University, Germany, and Seagate Research Center, Pittsburg, USA, where a part of this work was done. RS-R and UN thank Durham University, UK, for hospitality and support. MAW thanks ICMM, Madrid, Spain for hospitality and support, and acknowledges the EPSRC who supported this project under grant No. R040 318. The authors acknowledge useful discussions with A. Lyberatos (Seagate Research Center, Pittsburgh, USA).