format: text

text:
# Statistical Mechanics of Canonical-Dissipative Systems and Applications to Swarm Dynamics

## I Introduction

The collective motion of biological entities, like schools of fish, flocks of birds, herds of hoof animals, or swarms of insects, recently also attracted the interest of physicists. Here, the question how a long-range order between the moving entities can be established is of particular interest. Consequently, some of the more biologically centered questions of swarming behavior, namely about the resons of swarming or the group size dependence, have been dropped so far in physical swarm models. The main focus was rather on the emergence of coherent motion in a “swarm” of locally or globally coupled particles.

For the coupling different assumptions have been proposed (cf. also Sect. IV), such as the coupling of the particles’ individual orientation (i.e. direction of motion) to the mean orientation of the swarm (czirok-et-96, ; czirok-vicsek-00, ), or the coupling of the particles’ individual position to the mean position (center of mass) of the swarm (mikhailov-zanette-99, ). On the other hand, also some more local couplings have been considered, such as the coupling of the particle’s individual velocity to a local average velocity (toner-tu-95, ; vicsek-et-95, ; czirok-et-99, ; czirok-vicsek-00, ). A different class of models further assumes a local coupling of the particles via a self-consistent field that has been generated by them (lsg-mieth-rose-malch-95, ; lsg-fs-mieth-97, ). This models the case of chemical communication between the particles as widely found in biology. For example, the streaming behaviour of myxobacteria (stevens-95, ; ben-jacob-et-95-b, ; hoefer-et-95, ), or the directed motion of ants (calenbuhr-deneubourg-90, ; edelstein-keshet-94, ; fs-lao-family-97, ) or the coherent movement of cells (schienbein-gruler-95, ) have been described and modeled by means this kind of local coupling.

But long-range or short-range coupling of the particles is only one of the prerequisites that account for swarming. Another one is the *active motion* of the particles. Of course, particles can also move passively, driven by thermal noise, by convection, currents or by external fields. This kind of driving force however does not allow the particle to change its direction of motion, or velocity etc. *itself*. Recent models of *self-driven* particles which are used to simulate swarming behavior vicsek-et-95 ; albano-96 ; helbing-vicsek-99  usually just postulate that the entities move with a certain non-zero velocity, without considering the *energetic* implications of active motion. In order to do so, we need to consider that the many-particle system is basically an *open* system which is driven into non-equilibrium.

To this end, our approach to swarming is based on the theory of *canonical-dissipative systems*. This theory – which is not so well-known even among experts – results from an extension of the statistical physics of Hamiltonian systems to a special type of dissipative systems, the so-called canonical-dissipative systems graham-73 ; haken-73 ; hongler-ryter-78 ; eb-engelh-80 ; eb-81 ; feistel-eb-89 ; eb-00-b . The term *dissipative* means here that the system is non-conservative and the term *canonical* means, that the dissipative as well as the conservative parts of the dynamics are both determined by a Hamilton function $`H`$ (or a larger set of invariants of motion, see Sect. II for details).

This special assumption allows in many cases exact solutions for the distribution functions of many-particle systems, even in far-from-equilibrium situations. This was known already to pioneers as Poincare, Andronov and Bautin who gave exact solutions for a special class of nonlinear oscillators (andronov-witt-chaikin-65, ). In recent work (makarov-eb-velarde-00, ) the properties of canonical-dissipative systems were used to find exact solutions for *dissipative solitons* in Toda lattices which are pumped with free energy from external sources. The starting point was the well-known Toda theory of soliton solutions in one-dimensional lattices with a special nonlinear potential (exponential for compression and linear for expansion) of the springs toda-81 ; toda-83 , which was then extended towards a non-conservative system graham-73 ; eb-81 ; feistel-eb-89 . In another work eb-00  an application to systems of Fermi and Bose particles was given.

In this paper, we will apply the theory of canonical-dissipative systems to the dynamics of swarms. We start with an outline of the general theory in Sect. II, describing first the deterministic and then the stochastic approach. As a first step towards the dynamics of swarms, in Sect. III we investigate the *energetic conditions* of swarming, i.e. we discuss the conditions for active motion and their impact on the distribution function and the mean squared displacement of the particles. In Sect. IV we introduce *global interactions* between the particles which may account for swarming, i.e. for the maintenance of *coherent* motion of the ensemble. The global interactions are choosen with respect to the general theory of canonical-dissipative systems. By means of computer simulations we show how different kinds of coupling affect the overall behavior of the swarm. In Sect. V we conclude with some ideas how to generalize the approach presented in this paper.

## II General Theory of Canonical-Dissipative Systems

### II.1 Dynamics of Canonical-Dissipative Systems

Let us consider a mechanical many-particle sytem with $`f`$ degrees of freedom $`i=1,\mathrm{},f`$ and with the Hamiltonian $`H(q_1\mathrm{}q_f,p_1\mathrm{}p_f)`$ The corresponding equations of motion are

$$\frac{dq_i}{dt}=\frac{H}{p_i};\frac{dp_i}{dt}=\frac{H}{q_i}$$
(1)
Each solution of the system of equations (1)

$$p_i=p_i(t);q_i=q_i(t)$$
(2)
defines a trajectory on the plane $`H=E=`$const. This trajectory is determined by the initial conditions, and also the energy $`E=H(t=0)`$ is fixed due to the initial conditions. We construct now a *canonical-dissipative system* with the same Hamiltonian, where $`g(H)`$ denotes the *dissipation function*:

$$\frac{dp_i}{dt}=\frac{H}{q_i}g(H)\frac{H}{p_i}$$
(3)
In order to elucidate this kind of canonical-dissipative dynamics, we will consider different examples for the dissipation function in the following. In general we will only assume that $`g(H)`$ is a nondecreasing function of the Hamiltonian. The canonical-dissipative system, eq. (3), graham-73 ; eb-81 ; feistel-eb-89  does not conserve the energy because of the following relation:

$$\frac{dH}{dt}=g(H)\underset{i}{}\left(\frac{H}{p_i}\right)^2$$
(4)
Whether the total energy increases or decreases, consequently depends on the form of the dissipation function $`g(H)`$. In the simplest case, we may consider a constant friction:

$$g(H)=\gamma _0>0$$
(5)
As far as $`g(H)`$ is always positive, the energy always decays.

In a more interesting case the dissipative function $`g(H)`$ has a root for a given energy $`E_1`$: $`g(E_1)=0`$. Let us further assume that at least in certain neighbourhood of $`E_1`$ the function $`g(H)`$ is increasing. With these assumptions the states with $`H<E_1`$ are pumped with energy due to the negative dissipation, while energy is extracted from the states with $`H>E_1`$. That means, any given initial state with $`H(t=0)<E_1`$ will increase its energy until it reaches the shell $`H(t)=E_1`$, while any given initial state with $`H(t=0)>E_1`$ will decrease its energy until the shell $`H(t)=E_1`$ is reached, too. Therefore the solution of eq. (4) converges to the energy surface $`H=E_1`$. On this surface the solution of eq. (3) agrees with one of the possible solutions of the original Hamiltonian eq. (1) for $`H=E_1`$.

The simplest ansatz for $`g(H)`$ with a root $`g(H)=0`$ is a *linear* dissipation function:

$$g(H)=C(HE_1)$$
(6)
With respect to eq. (4) and the discussion above, it is obvious for this case that the process comes to rest when the shell $`H(t)=E_1`$ is reached. The relaxation time is proportional to the constant $`C^1`$.

We note that the linear dissipation function eq. (6) has found applications in Toda chains makarov-eb-velarde-00 . Here, the fact that on the shell $`H=E_1`$ the trajectory should obey the original *conservative* canonical dynamics, has been used to derive exact solutions for canonical-*dissipative* Toda systems.

A rather general *nonlinear* and non-decreasing dissipation function, which has been proposed in eb-00 , reads:

$$g(H)=\gamma _0\gamma _1\frac{(1+A)}{1+A\mathrm{exp}(\beta H)}$$
(7)
Here $`\gamma _0>0`$ represents the normal positive friction, $`\gamma _1>0`$ represents a kind of negative friction. $`A`$ is a dimensionless constant and $`\beta `$ is a parameter with the meaning of a reciprocal temperature. For $`\gamma _1\gamma _0`$ the friction is always positive, i.e. energy is extracted. For the opposite case, $`\gamma _1>\gamma _0`$, we have negative friction and the system is pumped with energy at least in some parts of the phase space. This allows to drive the system into situations far from equilibrium, therefore in the following we may assume $`\gamma _1>\gamma _0>0`$.

In the limit $`\beta 0`$ and $`A0`$ the dissipative function eq. (7) reduces to the linear case, eq. (6), discussed above. On the other hand, for small values of $`\beta `$ and finite $`A`$ we get

$$g(H)=\gamma _0\gamma _1\frac{(1+A)}{1+A+A\beta H}$$
(8)
which yields $`g(H)(\gamma _0\gamma _1)<0`$ for $`H0`$, and $`g(H)\gamma _0`$, if $`H\mathrm{}`$. This way, the existence of a root $`g(E_1)=0`$ is always guaranteed.

Before investigating a special case for $`g(H)`$ in Sect. III, let us introduce a generalization of the formalism. Instead of a driving function $`g(H)`$ which depends only on the Hamiltonian we may include the dependence on a larger set of invariants of motion, $`I_0,I_1,I_2,\mathrm{},I_s`$, for example:

* Hamilton function of the many-particle system: $`I_0=H`$
* total momentum of the many-particle system: $`I_1=P`$
* total angular momentum of the many-particle system: $`I_2=L`$

The dependence on this larger sets of invariants may be considered by defining a *dissipative potential* $`G(I_0,I_1,I_2,\mathrm{})`$. The canonical-dissipation equation of motion, eq. (3) is then generalized towards:

$$\frac{dp_i}{dt}=\frac{H}{q_i}\frac{G(I_0,I_1,I_2,\mathrm{})}{p_i}$$
(9)
Using this generalized canonical-dissipative formalism, by an appropriate choice of the dissipative potential $`G`$ the system may be driven to particular subspaces of the energy surface, e.g. the total momentum or the angular momentum may be prescribed. Different examples for this will be also discussed in Sect. IV.

### II.2 Stochastic Theory of Canonical-Dissipative Systems

We will now investigate an approach to the stationary probabilities which is based on Langevin and Fokker-Planck equations. The Langevin equations are obtained by adding a white noise term $`\xi _i(t)`$ to the deterministic eq. (3):

$$\frac{dp_i}{dt}=\frac{H}{q_i}g(H)\frac{H}{p_i}+\left(2D(H)\right)^{1/2}\xi _i(t)$$
(10)
The essential assumption is that both the strength of the noise, expressed in terms of $`D(H)`$ and the dissipation, expressed in terms of the friction function $`g(H)`$ depend only on the Hamiltonian $`H`$.

With respect to eq. (10), the corresponding Fokker-Planck equation for the probability distribution $`\rho (q_1\mathrm{}q_f,p_1\mathrm{}p_f)`$ of the many-particle system reads:

$`{\displaystyle \frac{\rho }{t}}+{\displaystyle p_i\frac{\rho }{q_i}}{\displaystyle \frac{H}{p_i}\frac{\rho }{p_i}}=`$
$`{\displaystyle \frac{}{p_i}\left[g(H)\rho +D(H)\frac{\rho }{p_i}\right]}`$ (11)
An exact stationary solution of eq. (II.2) reads:

$$\rho ^0(q_1\mathrm{}q_f,p_1\mathrm{}p_f)=Q^1\mathrm{exp}\underset{0}{\overset{H}{}}𝑑H^{}\frac{g(H^{})}{D(H^{})}$$
(12)
The derivative of $`\rho ^0`$ vanishes if $`g(H)=0`$, which means that the probability distribution is maximal at the surface $`H=E_1`$.

For the special case of a linear dissipation function, eq. (6), we find the stationary solution:

$$\rho ^0(q_1\mathrm{}q_f,p_1\mathrm{}p_f)=Q^1\mathrm{exp}\left(\frac{cH(2E_1H)}{2D}\right)$$
(13)
In the limit of very strong pumping, $`HE_1`$, this probability distribution reduces to a kind of microcanonical ensemble corresponding to the energy $`E_1`$.

The existence of exact solutions for the probability distribution, eq. (12), allows to derive several thermodynamic functions for the many-particle system, such as the *mean energy*:

$$U=Q^1𝑑HJ(H)H\mathrm{exp}\underset{0}{\overset{H}{}}𝑑H^{}\frac{g(H^{})}{D(H^{})}$$
(14)
where the Jacobian $`J(H)`$ is defined by:

$$dq_1\mathrm{}dq_fdp_1\mathrm{}dp_f=J(H)dH$$
(15)
The *entropy* follows from the Gibbs formula which yields here:

$`S=+k_B\mathrm{ln}Qk_B{\displaystyle }dHJ(H)\times `$
$`\times \left[Q^1\mathrm{exp}{\displaystyle \underset{0}{\overset{H}{}}}{\displaystyle \frac{g(H)}{D(H)}}\right]{\displaystyle \underset{0}{\overset{H}{}}}𝑑H^{}{\displaystyle \frac{g(H^{})}{D(H^{})}}`$ (16)
Further, the system has a Lyapunov functional $`K`$ which is provided by the Kullback entropy:

$$K[\rho ,\rho ^0]=𝑑q_1\mathrm{}𝑑q_f𝑑p_1\mathrm{}𝑑p_f\rho \mathrm{ln}\left(\frac{\rho }{\rho ^0}\right)$$
(17)
This gives explicitely:

$$K[\rho ,\rho _0]=\mathrm{ln}Q\frac{S}{k_B}+\underset{0}{\overset{H}{}}\frac{g(H)}{D(H)}$$
(18)
This functional is always nonincreasing.

If we want to generalize the description by means of the dissipative potential $`G(I_0,I_1,I_2,\mathrm{})`$ as used in eq. (9), the situation is more difficult. But at least if the noise strength $`D(H)`$ is a constant ($`D=`$const.), the corresponding Fokker-Planck equation

$`{\displaystyle \frac{\rho }{t}}+{\displaystyle p_i\frac{\rho }{q_i}}{\displaystyle \frac{H}{p_i}\frac{\rho }{p_i}}=`$
$`{\displaystyle \frac{}{p_i}\left[\frac{G(I_0,I_1,I_2,\mathrm{})}{p_i}\rho +D\frac{\rho }{p_i}\right]}`$ (19)
is still solvable. Due to the invariant character of the $`I_k`$, the l.h.s. of eq. (II.2) disappears for all functions of $`I_k`$. Therefore we have to search only for a function $`G(I_0,I_1,I_2,\mathrm{})`$, for which also the collision term disappears. This way we find the stationary solution of eq. (II.2):

$$\rho ^0(q_1\mathrm{}q_f,p_1\mathrm{}p_f)=Q^1\mathrm{exp}\left(\frac{G(I_0,I_1,I_2,\mathrm{})}{D}\right)$$
(20)
The derivative of $`\rho _0`$ vanishes if $`G(I_0,I_1,I_2,..)=\mathrm{min}.`$ which means that the probability is maximal on the attractor of the dissipative motion.

## III Energetic Conditions of Swarming

In this section, we want to apply the general description outlined above to swarm dynamics. Basically, at a certain level of abstraction a swarm can viewed as a many-particle system with some additonal coupling which would account for the typical correlated motion of the entities. In addition, also some energetic conditions must be satisfied in order to keep the swarm moving. Thus, *active motion* of the particles with a non-zero velocity is another basic ingredient for swarming. Since the conditions for active motion are dropped in many of the currently discussed swarm models albano-96 ; czirok-et-99 ; vicsek-et-95 , we first want to investigate the energetic conditions of swarming, before turning to the second ingredient, coupling of individual motion.

### III.1 Conditions for Active Motion

If we omit interactions between the particles, the Hamiltonian of the many-particle system is of the simple form:

$$H=\underset{i=1}{\overset{N}{}}H_i=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2m}$$
(21)
This allows to reduce the description level from the $`N`$-particle distribution function to the one-particle distribution function:

$$\rho (q_1\mathrm{}q_f,p_1\mathrm{}p_f)=\underset{i=1}{\overset{N}{}}\rho (𝒓_𝒊,𝒑_𝒊)$$
(22)
where the variable $`𝒓_i`$ is now used for the space coordinate of the particle $`i`$, and $`𝒑_i`$ stands for the momentum or, if $`m=1`$ is used in the following, for the *velocity* of the particle, respectively. The center of mass of the swarm is defined as

$$𝑹=\frac{1}{N}𝒓_i$$
(23)
whereas the mean momentum $`𝑷`$ and the mean angular momentum $`𝑳`$ are defined as:

$`𝑷`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle 𝒑_i}`$ (24)
$`𝑳`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle 𝑳_i}={\displaystyle \frac{1}{N}}{\displaystyle 𝒓_i\times 𝒑_i}`$ (25)
With $`g(H)=g(p_i^2)`$ and $`D=\mathrm{const}.`$, the Langevin eq. (10) for each particle $`i`$ reads in the absence of an external potential:

$$\dot{𝒓_i}=𝒑_i;\dot{𝒑}_i=g(p_i^2)𝒑_i+(2D)^{1/2}𝝃_i(t)$$
(26)
where the strength $`D`$ of the stochastic force results from the Einstein relation:

$$D=\gamma _0k_BT$$
(27)
$`T`$ is the temperature and $`k_B`$ is the Boltzmann constant.

For the dissipation function $`g(p^2)`$, we use again the general nonlinear ansatz of eq. (7) in the limit of small values of $`\beta `$ and finite $`A`$, eq. (8). By means of the substitutions

$$(1+A)=c;\gamma _1(1+A)=d_2s_0;\beta A=2d_2$$
(28)
and with $`H`$ of eq. (21), eq. (8) can be written in the form:

$$g(p^2)=\gamma _0\frac{s_0d_2}{c+d_2p^2}$$
(29)
We note that all noninteracting systems with $`g=g(p^2)`$ are of canonical-dissipative type.

Eq. (29) agrees with the velocity dependent non-linear friction function previously used in a model of active Brownian particles fs-eb-tilch-98-let ; eb-fs-tilch-99 . These are driven Brownian particles, which move due to the influence of a stochastic force, but additionally are pumped with energy due to a velocity-dependent dissipation function $`g(p^2)`$, eq. (29) which is plotted in Fig. 1.

The second term of the r.h.s. in eq. (29) which results from the pumping of energy, has been physically substantiated in our earlier work fs-eb-tilch-98-let ; eb-fs-tilch-99 . We have assumed that the Brownian particles are able to take up energy from the environment at a constant rate $`s_0`$, which can be stored in an *internal depot*, $`e`$. The internal energy can be converted into kinetic energy with a velocity dependent rate $`d(𝒑)=d_2p^2`$, which results in an additional acceleration of the Brownian particle in the direction of movement. The value of the internal energy depot may be further decreased due to internal dissipation processes, described by the constant $`c`$. Thus the resulting balance equation for the energy depot reads:

$$\frac{de}{dt}=s_0ce(t)d_2p^2e(t)$$
(30)
If we assume that the internal energy depot relaxes fast compared to the motion of the particle, we find the quasistationary value:

$$e_0=\frac{s_0}{c+d_2p^2}$$
(31)
which eventually leads to the second term of the r.h.s. in eq. (29).

Dependent on the parameters $`\gamma _0`$, $`d_2`$, $`s_0`$, $`c`$, the dissipation function, eq. (29), may have a zero, where the friction is just compensated by the energy supply. It reads in the considered case:

$$p_0^2=v_0^2=\frac{s_0}{\gamma _0}\frac{c}{d_2}$$
(32)
We see that for $`p<p_0`$, i.e. in the range of small momentums pumping due to negative friction occurs, as an additional source of energy for the Brownian particle. Hence, slow particles are accelerated, while the motion of fast particles is damped.

For $`s_0d_2<\gamma _0c`$, we find no real-valued root of eq. (32). This is the case of subcritical pumping, where the particle will move more or less like a simple Brownian particle. However, provided the existence of a non-zero momentum $`p_0`$, i.e. for a supercritical pumping, the particle will be able to move in a “high velocity” or *active* mode tilch-fs-eb-99 ; fs-tilch-eb-00 , which displays several non-trivial features of motion, as will be shown by means of computer simulations in the next section.

### III.2 Distribution Function and Dissipative Potential

Due to the pumping mechanism discussed above, the conservation of energy clearly does not hold for the particle, i.e. we now have a non-equilibrium, canonical-dissipative system as discussed in Sect. II. This results in deviations from the known Maxwellian velocity distribution of an equilibrium canonical system.

As pointed out in Sect. II.2, the probability density for the velocity $`\rho (𝒑,t)`$ obeys the Fokker-Planck eq. (II.2), which reads for the special case of the dissipation function, eq. (29), and in the absence of an external potential:

$$\frac{\rho (𝒑,t)}{t}=\frac{}{𝒑}\left[\left(\gamma _0\frac{d_2s_0}{c+d_2p^2}\right)𝒑\rho (𝒑,t)+D\frac{\rho (𝒑,t)}{𝒑}\right]$$
(33)
We mention that Fokker-Planck equations with nonlinear friction functions are discussed in detail in (klimontovich-95, ).

The stationary solution of eq. (33) is given by eq. (20), which reads in the considered case explicitely:

$`\rho ^0(𝒑)`$ $`=`$ $`C_0\mathrm{exp}\left({\displaystyle \frac{G_0(p^2)}{D}}\right)`$ (34)
$`=`$ $`C_0\left(1+{\displaystyle \frac{d_2p^2}{c}}\right)^{\frac{s_0}{2D}}\mathrm{exp}\left({\displaystyle \frac{\gamma _0}{2D}}p^2\right)`$ (35)
where $`C_0`$ results from the normalization condition. $`G_0(p^2)`$ is the special form of the dissipation potential $`G(I_0,I_1,I_2,\mathrm{})`$ considering only $`I_0=H`$ as the invariant of motion, and further $`H`$ as given by eq. (21). It reads explicitely:

$$G(I_0)=G_0(p^2)=\gamma _0\frac{p^2}{2}\frac{s_0}{2}\mathrm{ln}\left(1+\frac{d_2p^2}{c}\right)$$
(36)
Compared to the Maxwellian velocity distribution of “simple” Brownian particles, a new prefactor appears now in eq. (35) which results from the additional pumping of energy. For a subcritical pumping, $`s_0d_2<c\gamma _0`$, where we do not find a real-valued root of the dissipation function, eq. (29), only an *unimodal velocity distribution* results, centered around the maximum $`𝒑_0=0`$. However, for supercritical pumping, $`s_0d_2>c\gamma _0`$, if the root of $`\gamma (p^2)`$ is real, we find a *crater-like velocity distribution*, which indicates strong deviations from the Maxwell distribution erdmann-et-00 .

This is also shown in Fig. 2 which presents computer simulations of the velocity distribution of 10.000 particles after a sufficiently long time (only the x-dimension of the 2-d simulation is shown). For the supercritical case, two distinct peaks of the velocity distribution are found at $`p_x=\{0.63,+0.63\}`$. The values of these maxima agree with the deterministic result for the stationary velocity, eq. (32).

We note that non-Maxwellian velocity distributions for active motion have been also observed experimentally in cells, such as granulocytes franke-gruler-90 ; schienb-gruler-93 .

### III.3 Mean Squared Displacement and Stationary Values

As Fig. 2 shows, the momentum distribution is centered around $`𝒑=0`$ both for subcritical and supercritical pumping. If we consider a nearly spherical swarm of particles in the two-dimensional space $`𝒓=\{x,y\}`$ as in the computer simulations in this section, its center of mass, eq. (23) and mean momentum, eq. (24) will come to rest. Thus they are not affected by the pumping, but other quantities are, such as the mean squared displacement:

$$\mathrm{\Delta }R^2(t)=\left(\frac{1}{N}\left[𝒓_i(t)𝒓_i(0)\right]\right)^2$$
(37)
In the limit of pure Brownian motion, it is known that the mean squared displacement increases in time as:

$$\mathrm{\Delta }R^2(t)=4D_rt$$
(38)
where $`D_r=k_BT/\gamma _0=D/\gamma _0^2`$ is the spatial diffusion coefficient. Thus, eq. (38) will be the lower limit for subcritical pumping of the particles. Contrary, in the case of supercritical pumping it has been shown mikhailov-meinkoehn-97 ; erdmann-et-00  that the mean squared displacement will grow in time approximately as

$$\mathrm{\Delta }R^2(t)=\frac{2v_0^4}{D}t$$
(39)
where $`v_0`$ is given by eq. (32). Consequently, the diffusion coefficient $`D_r`$ in eq. (38) has for the case of supercritical pumping to be replaced by an *effective* spatial diffusion coefficient:

$$D_r^{\mathrm{eff}}=\frac{v_0^4}{2D}=\frac{1}{2D}\left(\frac{s_0}{\gamma _0}\frac{c}{d_2}\right)^2.$$
(40)
This result holds for noninteracting particles in the limit of relatively weak noise intensity $`D`$ and/or strong pumping and will therefore give an upper limit for $`\mathrm{\Delta }R^2(t)`$. We note the high sensitivity with respect to noise expressed in the scaling with $`(1/D)`$.

Fig. 3 shows the mean squared displacement of a swarm of 2.000 particles both for the case of sub- and supercritical pumping together with the theoretical results of eqs. (38), (39). We see that for long times the computer simulations for supercritical pumping agree very well with eq. (39).

Another quantity affected by the sub/supercritical pumping is the stationary velocity $`p_0^2=v_0^2`$, eq. (32). In 2d, these stationary velocities define a cylinder, $`p_x^2+p_y^2=p_0^2`$, in the four-dimensional state space $`\{x,y,p_x,p_y\}`$ which attracts all deterministic trajectories of the dynamic system erdmann-et-00 . Fig. 4 shows the results of computer simulations for $`p^2(t)`$ for the case of supercritical pumping. The convergence toward the theoretical result, eq. (32) can be clearly observed.

## IV Globally Coupled Swarms

So far we have neglected any coupling within the many-particle ensemble. This leads to the effect that the swarm eventually disperses in the course of time, whereas a “real” swarm would maintain its coherent motion. A common way to introduce correlations between the moving particles in physical swarm models is the coupling to a mean value. For example, in (czirok-et-96, ; czirok-vicsek-00, ) the coupling of the particles’ individual *orientation* (i.e. direction of motion) to the mean orientation of the swarm is discussed. Other versions assume the coupling of the particles’ velocity to a *local average velocity*, which is calculated over a space interval around the particle czirok-et-99 ; czirok-vicsek-00 .

### IV.1 Coupling to the Center of Mass

In this paper, we are mainly interested in *global couplings* of the swarm, which fit it into the theory of canonical-dissipative systems outlined in Sect. II. As the most simple case we may first discuss the global coupling of the swarm to the *center of mass*, eq. (23). That means the particle’s position $`𝒓_i`$ is related to the mean position of the swarm $`𝑹`$ via a potential $`U(𝒓_i,𝑹)`$. For simplicity, we may assume a parabolic potential:

$$U(𝒓_i,𝑹)=\frac{a}{2}\left(𝒓𝑹\right)^2$$
(41)
The harmonic potential generates a force directed to the center of mass which can be used to control the dispersion of the swarm. It reads in the considered case:

$$\mathbf{}U(𝒓)=a\left(𝒓_i𝑹\right)=\frac{a}{N}\underset{j=1}{\overset{N}{}}\left(𝒓_i𝒓_j\right)$$
(42)
With eq. (42), the corresponding Langevin eq. (26) of the many-particle system reads explicitely:

$$\dot{𝒑}_i=g(p_i^2)𝒑_i\frac{a}{N}\underset{j=1}{\overset{N}{}}\left(𝒓_i𝒓_j\right)+(2D)^{1/2}𝝃_i(t)$$
(43)
Hence, in addition to the dissipation function there is now an attractive force between each two particles $`i`$ and $`j`$ which depends linearly on the distance between them. With respect to the harmonic interaction potential eq. (41), we call such a swarm a *harmonic* swarm eb-fs-01-pas .

Strictly speaking, the dynamical system of eq. (43) is not a canonical-dissipative one, but as shown in (eb-fs-tilch-99, ) it may be reduced to this type by some approximations, which will be also discussed below. We note that this kind of swarm model has been previously investigated in mikhailov-zanette-99  for the one-dimensional case, however with a different dissipation function $`g(p^2)`$, for which we use eq. (29) again. Obviously, as shown in Sect. III.2 swarming will occur only for supercritical conditions.

With the assumed coupling to the center of mass, $`𝑹`$, the motion of the swarm can be considered as a superposition of two motions: (i) the motion of the center of mass itself, and (ii) the motion of the particles relative to the center of mass. Taking into account that the noise acting on the different particles is not correlated, the center of mass for the assumed coupling obeys a force-free motion,

$$\dot{𝑹}=𝑷;\dot{𝑷}=\frac{1}{N}\underset{i=1}{\overset{N}{}}g(p_i^2)𝒑_i$$
(44)
Because of the nonlinearities in the dissipation function $`g(p^2)`$ both motions (i) and (ii) cannot be simply separated. The term $`g(p^2)`$ vanishes only for two cases: the trivial one which is free motion without dissipation/pumping, or the case of supercritical pumping where $`p_i^2=p_0^2`$, eq. (32) for each particle. Then, the mean momentum becomes an invariant of motion, $`𝑷(t)=𝑷_0=\mathrm{const}.`$ and the center of mass moves according to $`𝑹(t)=𝑹(0)+𝑷_0(t)`$. This behavior may also critically depend on the initial conditions of the particles, $`𝒑_i(0)`$, and shall be investigated in more detail now.

In mikhailov-zanette-99  an approximation for the mean velocity $`𝑷(t)`$ of the swarm in one dimension is discussed which shows the existence of two different asymptotic solutions dependent on the noise intensity $`D`$ and the initial momentum $`𝒑_i(0)`$ of the particles. Below a critical noise intensity $`D_c`$, the initial condition $`p_i^2(0)>p_0^2`$ leads to a swarm the center of which travels with a constant non-trivial mean velocity, while for the initial condition $`p_i^2(0)<p_0^2`$ the center of the swarm is at rest.

We can confirm these findings by means of two-dimensional computer simulations presented in Fig. 5 and Fig. 6, which show the mean squared displacement, the average squared velocity of the swarm and the squared mean velocity of the center of mass for the two different initial conditions.

For (a) $`p_i^2(0)>p_0^2`$ we find a continuous increase of $`\mathrm{\Delta }R^2(t)`$ (Fig. 5 a), while the velocity of the center of mass reaches a constant value: $`P^2(t)=\left(N^1_i𝒑_i(t)\right)^2p_0^2`$, known as the stationary velocity of the force-free case, eq. (32). The average squared velocity reaches a constant non-trivial value, too, which depends on the noise intensity and the initial conditions, $`p_i^2(0)>p_0^2`$, i.e. on the energy initially supplied (cf Fig. 6 top).

For (b) $`p_i^2(0)<p_0^2`$ we find that the mean squared displacement after a transient time reaches a constant value, i.e. the center of mass comes to rest (Fig. 5 b), which corresponds to $`P^2(t)0`$ in Fig. 6 (bottom). In this case, however, the averaged squared velocity of the swarm reaches the known stationary velocity, $`p^2(t)=1/Np_i^2(t)p_0^2`$. Consequently, in this case the energy provided by the pumping goes into the motion *relative to the center of mass* while the motion of the center of mass is damped out (cf. also (eb-fs-01-pas, )). Thus, in the following we want to investigate the relative motion of the particles in more detail.

Using relative coordinates, $`\{x_i,y_i\}𝒓_i𝑹`$, the dynamics of each particle in the two-dimensional space is described by four coupled first-order differential equations:

$`\dot{x}_i`$ $`=`$ $`p_{xi}P_x`$
$`\dot{p}_{xi}\dot{P}_x`$ $`=`$ $`g\left(p_i^2\right)p_{xi}ax_i+(2D)^{1/2}\xi _i(t)`$
$`\dot{y}_i`$ $`=`$ $`p_{yi}P_y`$
$`\dot{p}_{yi}\dot{P}_y`$ $`=`$ $`g\left(p_i^2\right)p_{yi}ay_i+(2D)^{1/2}\xi _i(t)`$ (45)
For $`𝑷=0`$, i.e. for the initial conditions $`p_i^2<p_0^2`$ and sufficiently long times, this dynamics is equivalent to the motion of free (or uncoupled) particles in a parabolic potential $`U(x,y)=a(x^2+y^2)/2`$ with the origin $`\{0,0\}`$. Thus, within this approximation the system becomes a canonical-dissipative system again, in the strict sense used in Sect. II.

Fig. 7 presents computer simulations of eq. (45) for the relative motion of the particle swarm in the parabolic potential. <sup>1</sup><sup>1</sup>1The reader is invited to view a *movie* of the respective computer simulations (t=0–130), which can be found at http://ais.gmd.de/frank/swarm1.html  (Note that in this case all particles started from the same position slightly *outside* the origin of the parabolic potential. This has been chosen in order to make the evolution of the different branches more visible.) As the snapshots of the spatial dispersion of the swarm show, we find after an inital stage the occurence of two branches of the swarm which results from a *spontaneous symmetry break* (cf. Fig. 7 top). These two branches will after a sufficiently long time move on two limit cycles (as already indicated in Fig. 7 bottom). One of these limit cycles refers to the left-handed, the other one to the right-handed direction of motion in the 2d-space.

This finding also agrees with the the theoretical investigations of the deterministic case eb-fs-tilch-99  which showed the existence of a limit cycle with the amplitude

$$r_0=\left|p_0\right|a^{1/2}$$
(46)
provided the relation $`s_0d_2>\gamma _0c`$ is fulfilled. In the small noise limit the radius of the limit cycles shown in Fig. 7 agrees with the value of $`r_0`$. Further, Fig. 6 has shown that the averaged squared velocity $`p^2(t)`$ of the swarm indeed approaches the theoretical value of eq. (32).

The existence of two opposite rotational directions of the swarm can be also clearly seen from the distribution of the angular momenta $`L_i`$ of the particles. Fig. 8 shows the existence of a bimodal distribution for $`\rho (L)`$. (The observant reader may notice that each of these peaks actually consists of two subpeaks resulting from the initial conditions, which are still not forgotten at $`t=99`$). Each of the main peaks is centered around the theoretical value

$$\left|𝑳\right|=L_0=r_0p_0$$
(47)
where $`r_0`$ is given by eq. (46) and $`p_0`$ is given by eq. (32).

The emergence of the two limit cycles means that the dispersion of the swarm is kept within certain spatial boundaries. This occurs after a transient time used to establish the correlation between the individual particles. In the same manner as the motion of the particles becomes correlated, the motion of the center of mass is slowed down until it comes to rest, as already shown in Fig. 6.

This however is not the case if the initial conditions $`p_i^2(0)>p_0^2`$ are chosen. Then, the energy provided by the pumping does not go completely into the relative motion of the particles and the establishment of the limit cycles as discussed above. Instead, the center of mass keeps moving as shown in Fig. 5, while the swarm itself does not establish an internal order. Fig. 9 displays a snapshot of the relative positions of the particles in this case (note the different scales of the axes compared to Fig. 7).

### IV.2 Coupling via Mean Momentum and Mean Angular Momentum

In the following we want to discuss two other ways of *global* coupling of the swarm which fit into the general framework of canonical-dissipative systems outlined in Sect. II. There, we have introduced a dissipative potential $`G(I_0,𝑰_1,𝑰_2,\mathrm{})`$ which depends on the different invariants of motion, $`𝑰_i`$. So far, we have only considered $`I_0=H`$, eq. (21) in the swarm model. If we additionally include the mean momentum $`𝑰_1=𝑷`$, eq. (24) as the first invariant of motion, the dissipative potential may read:

$`G(I_0,𝑰_1)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}G_0(p_i^2)+G_1(𝑷)`$ (48)
$`G_1(𝑷)`$ $`=`$ $`{\displaystyle \frac{C_P}{2}}\left({\displaystyle \underset{i=1}{\overset{N}{}}}𝒑_iN𝑷_\mathrm{𝟏}\right)^2`$ (49)
Here, $`G_0(p^2)`$ is given by eq. (36). The stationary solution of the probability distribution, $`\rho ^0(𝒑)`$ is again given by eq. (20). In the absence of an external potential $`U(𝒓)`$, the Langevin equation that corresponds to the dissipation potential of eq. (48) reads now:

$`\dot{𝒓_i}`$ $`=`$ $`𝒑_i`$
$`\dot{𝒑}_i`$ $`=`$ $`g(p_i^2)𝒑_iC_P\left({\displaystyle 𝒑_i}N𝑷_\mathrm{𝟏}\right)`$ (50)
$`+(2D)^{1/2}𝝃_i(t)`$
The term $`G_1(𝑷)`$ is choosen that way that it may drive the system towards the prescribed momentum $`𝑷_1`$, where the relaxation time is proportional to $`C_P^1`$. If we would have a vanishing dissipation function, i.e. $`g=0`$ for $`p_i^2=p_0^2`$, it follows from eq. (IV.2) for the mean momentum:

$$𝑷(t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}𝒑_i(t)=𝑷_1+\left[𝑷(0)𝑷_1\right]e^{C_pt}$$
(51)
The existence of two terms $`G_0`$ and $`G_1`$ however could lead to competing influences of the resulting forces, and a more complex dynamics of the swarm results. As before, this may also depend on the initial conditions, i.e. $`P_1^2p_0^2`$ or $`P_1^2p_0^2`$.

Fig. 10 shows the squared velocity of the center of mass, $`P^2(t)`$ and the averaged squared velocity of the swarm $`p^2(t)`$ for $`P_1^2p_0^2`$, We find an intermediate stage, where both velocities are equal, before the global coupling drives the mean momentum $`𝑷`$ towards the prescribed value $`𝑷_1`$, i.e. $`P^2(t)(P_{1x}^2+P_{1y}^2)`$. On the other hand, $`p^2(t)p_0^2`$, as we have found before for the force-free case and for the linearly coupled case for similar initial conditions. The noticeable decrease of $`P^2`$ after the initial time lag can be best understood by looking at the spatial snapshots of the swarm provided in Fig. 11. For $`t=10`$, we find a rather compact swarm where all particles move into the same (prescribed) direction. For $`t=50`$, the correlations between the particles have already become effective, which means the swarm begins to establish a circular front, which however does not become a full circle <sup>2</sup><sup>2</sup>2The movie that shows the respective computer simulations for t=0–100 can be found at http://ais.gmd.de/frank/swarm2.html. Eventually, we find again that the energy provided by the pumping goes into the motion of the particles relative to the center of mass, while the motion of the center of mass itself is driven by the prescribed momentum.

For the initial condition $`P_1^2p_0^2`$ the situation is different again, as Fig. 12 shows. Apparently, both curves are the same for a rather small noise intensity, i.e. $`P^2(t)=p^2(t)`$ are both equal, but different from $`p_0^2`$, eq. (32) and the prescribed momentum $`P_1^2`$. This can be only realized if all particles move in parallel into the same direction. Thus, a snapshot of the swarm would much look like the top part of Fig. 11.

Eventually, we may also use the second invariant of motion, $`𝑰_2=𝑳`$, eq. (25) for a global coupling of the swarm. In this case, the dissipative potential may be defined as follows:

$`G(I_0,𝑰_2)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}G_0(p_i^2)+G_2(𝑳)`$ (52)
$`G_2(𝑳)`$ $`=`$ $`{\displaystyle \frac{C_L}{2}}\left({\displaystyle 𝒓_i\times 𝒑_i}N𝑳_\mathrm{𝟏}\right)^2`$ (53)
$`G_0(p^2)`$ is again given by eq. (36). The term $`G_2(𝑳)`$ shall drive the system to a prescribed angular momentum $`𝑳_1`$ with a relaxation time proportional to $`C_L^1`$.

$`𝑳_1`$ can be used to break the symmetry of the swarm towards a prescribed rotational direction. In Sect. IV.1 we have observed the spontaneous occurence of lefthand and righthand rotations of a swarm of linearly coupled particles. Without an additional coupling, both rotational directions are equally probable in the stationary limit. Considering both the parabolic potential $`U(𝒓,𝑹)`$, eq. (41) and the dissipative potential eq. (52), the corresponding Langevin equation may read now:

$`\dot{𝒓_i}`$ $`=`$ $`𝒑_i`$
$`\dot{𝒑}_i`$ $`=`$ $`g(p_i^2)𝒑_i{\displaystyle \frac{a}{N}}{\displaystyle \underset{j=1}{\overset{N}{}}}\left(𝒓_i𝒓_j\right)`$ (54)
$`+C_L𝒓_i\times \left({\displaystyle 𝒓_i\times 𝒑_i}N𝑳_\mathrm{𝟏}\right)`$
$`+(2D)^{1/2}𝝃_i(t)`$
The computer simulations shown in Fig. 13 clearly display an unimodal distribution of the angular momenta $`L_i`$ of the particles, which can be compared to Fig. 8 without coupling to the angular momentum. Consequently, we find in the long time limit only one limit cycle corresponding to the movement of the swarm into the same rotational direction. The radius $`r_0`$ of the limit cycle is again given by eq. (46).

We would like to add that also in this case the dynamics depends on the initial condition, $`𝑳_1`$. For simplicity, we have assumed here $`\left|𝑳_1\right|=L_0=r_0p_0`$, eq. (47) which is also reached by the mean angular momentum $`𝑳`$ in the course of time (cf. Fig. 13). For initial conditions $`\left|𝑳_\mathrm{𝟏}\right|L_0`$, there is of course no need for the rotation of *all* particles into the same direction. Hence we observe both left- and righthanded rotations of the particles with different shares, so that the mean angular momentum is still $`𝑳𝑳_1`$. This results in a broader distribution of the angular momenta of the particles instead of the clear unimodal distribution shown in Fig. 13. For initial conditions $`\left|𝑳_\mathrm{𝟏}\right|L_0`$ on the other hand, the stable rotation of the swarm breaks down after some time, since the driving force $`𝑳𝑳_1`$ tends to destabilize the attractor $`𝑳𝑳_0`$. This effect will be investigated in a forthcoming paper, together with some combined effects of the different global couplings.

## V Discussion

Eventually, we can also combine the different global couplings discussed above by defining the dissipation potential as:

$`G(I_0,𝑰_1,𝑰_2)`$ $`=`$ $`G(p^2,𝑷,𝑳)`$ (55)
$`=`$ $`G_0(p^2)+G_1(𝑷)+G_2(𝑳)`$
$`G_0(p^2)`$ is given by eq. (36), $`G_1(𝑷)`$ by eq. (49) and $`G_2(𝑳)`$ by eq. (53). Considering further an additional – external or interaction – potential, the corresponding Langevin equation can be written in the more general form:

$`\dot{𝒓_i}`$ $`=`$ $`𝒑_i`$
$`\dot{𝒑}_i`$ $`=`$ $`g(p_i^2)𝒑_i{\displaystyle \frac{}{𝒓_i}}\left[\alpha _0U(𝒓,𝑹)\right]`$ (56)
$`{\displaystyle \frac{}{𝒑_i}}\left[\alpha _1\left(𝑷𝑷_\mathrm{𝟏}\right)^2+\alpha _2\left(𝑳𝑳_\mathrm{𝟏}\right)^2\right]`$
$`+(2D)^{1/2}𝝃_i(t)`$
The mean momentum $`𝑷`$ and mean angular momentum $`𝑳`$ are given by eqs. (24), (25); whereas the constant vectors $`𝑷_1`$ and $`𝑳_1`$ are used to break the spatial or rotational symmetry of the motion toward a preferred direction. The different constants $`\alpha _i`$ may decide whether the respective influence of the conservative or dissipative potential is effective or not, they further determine the time scale when the global coupling becomes effective. The term $`g(p^2)`$, eq. (29) on the other hand considers the energetic conditions for the active motion of the swarm, i.e. it determines whether the particle of the swarm are able to “take off” at all.

The combination of the different types of coupling may lead to a rather complex swarm dynamics, as already indicated in the examples discussed in this paper. In particular, we note that the different terms may have competing influences on the swarm, which then would lead to a “frustrated” dynamics with many possible attractors.

In this paper, we have basically restricted the investigation of the swarm dynamics to global couplings which fit into (or can be reduced within some approximations to) the general outline of canonical-dissipative systems. Finally, we want to add some comments on that. On one hand, it is possible to extend this kind of approach also to other invariants of motion, this way e.g. covering previous investigations of swarms coupled via the mean orientation of the particles (czirok-et-96, ; czirok-vicsek-00, ). On the other hand, we want to underline that canonical-dissipative systems are a theoretical class of models, where *both* conservative and dissipative elements of the dynamics are determined by invariants of the mechanical motion. Thus, from this perspective, a more realistic swarm dynamics may be also based on less restrictive assumptions.

The advantage to use canonical-dissipative systems as a framework for swarm dynamics is given by the fact that in many cases the rather complex dynamics can be mapped to an analytically tractable model. With the Hamiltonian theory of many-particle systems as a starting point, we are able to extend known solutions for conservative systems to non-conservative systems. This allows us to to construct a canonical-dissipative system, the solutions of which converges to the solution of the conservative system with given energy. That means, for given initial conditions, we could predict the asymptotic solution of the canonical-dissipative dynamics by means of the solutions of the Hamiltonian equations on the respective energy surface.

In addition to these theoretical considerations which have a value of their own, we want to note that the framework of canonical-dissipative systems still covers important features of real (biological) systems, such as energy take-up and dissipation. The general description outlined in this paper allows us to gradually add more and more complexity to the swarm model, this way bridging the gap between a known (physical) dynamics to a more complex (biological) dynamics (eb-fs-01-pas, ; eb-fs-01, ). Some hints for this shall be given at the end.

On the level of the “individual” particles, the key dynamics of the model is given by a modified Langevin equation that, in addition to stochastic influences, also considers other forces on the particle, resulting e.g. from external potentials, interactions, use of stored energy (eb-fs-tilch-99, ), or influences of a self-consistent field (lsg-fs-mieth-97, ; lsg-mieth-rose-malch-95, ) that already exceed the framework of canonical-dissipative systems.

The consideration of an *external* potential also allows to model the *spatial environment* of the swarm, for instance to consider *obstacles* (fs-eb-tilch-98-let, ). Additionally, we can also consider that the pumping of energy for the particles is restricted to certain spatial domains which model *food sources*. In this case the dissipation function $`g(H)`$ also becomes a spatial function. Such an extension has been already discussed in (steuern-et-94, ; fs-eb-tilch-98-let, ) and can be also implemented in the swarm model discussed here. Eventually, we note that the genuin particle-based approach to collective phenomena used in this paper is not restricted to biological systems, but also applicable for describing and simulating complex interactive systems in a wide range of applications, even in economics and social systems (fs-book-01, ).

###### Acknowledgements.

The authors thank J. Dunkel and U. Erdmann for discussions.


text:
# Resonance accretion on the magnetized rotating neutron star

## 1 Introduction

Rossi X-Ray Timing Explorer (RXTE) in 1996 discovered the quasi-periodic oscillations (QPO) of X-ray emissions in low-mass X-ray binaries (LMXB) in the band of $`1`$ kHz frequencies (van der Klis et al. 1996; van der Klis 2000). They are narrow, the ratio of the line’s width to the central frequency can be as low as $`10^2`$, and in most cases there were observed two lines. The existence of such lines means that the accretion disk approaches very close to the surface of the rotating neutron star. Indeed, the Keplerian frequency

$$\nu _\mathrm{K}=\frac{1}{2\pi }\left(\frac{MG}{r^3}\right)^{1/2},$$
(1)
where $`M`$ is the neutron star’s mass and $`r`$ is the distance from the star’s center, is of the order of 1 kHz at the distance of $`1.710^6`$ cm for the typical neutron star mass $`M1.4M_{}`$. But estimated neutron star radius $`R`$ is of about $`10^6`$ cm, and we see that the inner edge of the accretion disk can approach almost to the star surface. If we estimate the kinetic energy of a proton near the star’s radius ($`rR`$) it turns out of $`310^5`$ erg or 20 MeV. For the characteristic proton’s density in the accretion flow near the star surface, corresponding to the accretion rate $`\dot{M}10^{17}`$ g/sec ($`L_xL_{\mathrm{Edd}}`$), $`n_\mathrm{p}810^{17}`$ cm-3, we obtain the energy density of the accreting plasma of

$$2.510^{13}\mathrm{erg}/\mathrm{cm}^3.$$

It seams that for the existence of the disk in the region near the star’s surface, the star’s magnetic field should not exceed the value of

$$\frac{B^2}{8\pi }<,\mathrm{or}B<2.510^7\mathrm{G}.$$

But this low value of magnetic field on the neutron star surface contradicts with the magnitude of the magnetic field observed for the milliseconds pulsars $`B10^{910}`$ G, which are considered originate from the LMXB. But these estimations are incorrect. The plasma destroys the magnetic field structure when its thermal energy exceeds the magnetic field energy, but not when its kinetic flow energy exceeds the magnetic one. Moving across the magnetic field lines plasma polarizes so that the electric field across the magnetic one appears. In the case of a thin disk the radial electric field

$$E_r=\frac{1}{c}[v_\mathrm{K}B]$$
(2)
is enough to carry the flow across the magnetic field with the Keplerian velocity $`v_\mathrm{K}=2\pi \nu _\mathrm{K}r`$. For the dipole magnetic field ($`Br^3`$) the disk of the temperature $`T`$ destroys the star’s magnetic field at the distances

$$r>r_A=R\left[\frac{B^2R^2}{4\dot{M}(T/m_\mathrm{p})^{1/2}}\right]^{1/4}$$
(3)
If we substitute to (3), the characteristic energy of X-ray emission ($`1`$ keV) as the maximum disk temperature we obtain

$$r>20R\left(\frac{B}{10^9\mathrm{G}}\right)^{1/2}.$$

Thus, we see that in the region where QPO of 1 kHz frequency can originate the structure of the magnetic field of the star is not disturbed by the accreting plasma. We can consider that the structure of the magnetic field near the star surface($`r10R`$) is given by the interior of the neutron star and has the dipole configuration. Besides, the accretion changes the initial neutron star magnetic field $`B10^{12}`$ G up to $`B10^9`$ G. And it is naturally to consider the direction of the dipole magnetic momentum $`\mu `$ coinciding with the axis of star’s rotation $`\mathrm{\Omega }`$.

We assume that the accreting disk approaches close to the neutron star surface. Plasma flows from the disk to the neutron star along magnetic field lines and form so-called “magnetic field plasma tubes”, which are shaped as magnetic field lines, are filled with plasma and have a small thickness.

The narrow structure of the QPO of 1 kHz frequency suggests the resonance nature of its origin. We have the idea that it is a resonance between Keplerian rotation of the disk and an eigen modes of the oscillation of magnetic field tubes, frozen to the conducting star surface.

Such model is more close to the sonic point beat frequency model (Miller et al. 1998), where the resonance takes place not with the Keplerian rotation of the disk, but with the radial motion of accreting flow coinciding with the sonic velocity. The other model based on the eigen modes of the disk in the rotating frame was proposed by Titarchuk & Osherovich (1999). There is more exotic mechanism proposed by Klein et al. (1998) where QPOs are the oscillations of photon’s bubbles.

## 2 Magnetic tubes oscillations

In this section we discuss the oscillations of the magnetic field lines with plasma frozen to them. Discovering this model it is convenient to use the dipole coordinates ($`l,a`$) introduced by the definition:

$`B`$ $`=`$ $`l;`$
$`B`$ $`=`$ $`[a,\varphi ].`$ (4)
$`l`$ is the coordinate along the magnetic field line

$$l=\frac{(\mu r)}{r^3}=\mu \frac{z}{(\varrho ^2+z^2)^{3/2}},$$
(5)
$`a`$ numerates the magnetic field surfaces

$$a=\mu \frac{\varrho ^2}{(\varrho ^2+z^2)^{3/2}}.$$
(6)
$`a=\mathrm{const}`$ is the equation for the magnetic field lines

$$z^2=(\frac{\mu }{a})^{2/3}\varrho ^{4/3}\varrho ^2,\text{or}r=\frac{\mu }{a}\mathrm{sin}^2\theta .$$
(7)
Here $`(\varrho ,z)`$ are the cylindrical coordinates, $`(r,\theta )`$ are the spherical ones. The value of $`a`$ is also connected with the maximum deflection $`\varrho _{max}`$ of the dipole magnetic field line from the star’s center $`a=\mu /\varrho _{max}`$. For different magnetic surfaces the coordinate $`a`$ changes from 0 to $`\mu /R`$. The coordinate $`l`$ changes for a given $`a`$, from $`(\mu ^2\mu Ra)^{1/2}/R^2`$ to $`(\mu ^2\mu Ra)^{1/2}/R^2`$. The point $`l=0`$ corresponds to the top of the magnetic field line. Three coordinates $`(l,a,\varphi )`$ ($`\varphi `$ is the azimuthal angle) are the curvilinear orthogonal coordinate system.

Because in the near ($`r<10R`$) region $`B^2/8\pi `$ the plasma can be considered as cold, the plasma density is proportional to the magnetic field strength, $`n_\mathrm{p}=\lambda B`$. Here $`\lambda `$ is the constant not depending on the coordinates. The equations of the ideal MHD in rotating with angular velocity $`\mathrm{\Omega }`$ frame ($`\mathrm{\Omega }`$ is the star rotation frequency) can be presented for the electric current $`j`$ as follows

$$j=\frac{v_A^2}{\omega ^2B}\mathrm{rot}\mathrm{rot}\left\{\frac{1}{B}\left[B[jB]\right]\right\},$$
(8)
where $`v_A`$ is the Alfven velocity

$$v_A^2=B^2/4\pi n_\mathrm{p}m_\mathrm{p}=B/4\pi \lambda m_\mathrm{p}.$$

We consider here that all the perturbed values $`j`$, $`\delta B`$, $`E`$ depend on the time $`t`$ and the azimuthal angle $`\varphi `$ as $`\mathrm{exp}\left\{i\omega t+im\varphi \right\}.`$

The electric current can be presented through its covariant components $`j_1,j_2,j_3`$

$$j=j_1l+j_2a+j_3\varphi ,$$

for which from (8) we obtain

$$\begin{array}{c}\hfill j_2\left(1\frac{v_A^2}{\omega ^2\varrho ^2}m^2\right)+\frac{v_A^2}{\omega ^2\varrho ^2B}\frac{}{l}\left[\varrho ^2B^2\frac{}{l}\left(j_2B\right)\right]=\\ \hfill im\frac{v_A^2}{\omega ^2\varrho ^2B}\frac{}{a}\left(j_3B\right);\end{array}$$
(9)
$$\begin{array}{c}\hfill j_3+\frac{v_A^2\varrho ^2B}{\omega ^2}\left[\frac{^2}{a^2}\left(j_3B\right)+\frac{}{l}\frac{1}{\varrho ^2}\frac{}{l}\left(j_3B\right)\right]=\\ \hfill im\frac{v_A^2\varrho ^2B}{\omega ^2}\frac{}{a}\left(j_2B\right);\end{array}$$
(10)
$$j_1=\frac{v_A^2}{\omega ^2B}\left[\frac{}{a}\varrho ^2B^2\frac{}{l}\left(j_2B\right)+\frac{im}{\varrho ^2}\frac{}{l}\left(j_3B\right)\right].$$
(11)
The boundary conditions for (9, 10, 11) are on the star surface $`r=R`$ . Due to the ideal conductivity of neutron star matter the tangential components of the electric field have to be zero. The electric field of the disturbance is proportional to the electric current

$$E=\frac{iB}{\omega c^2\lambda }\left(j_2a+j_3\varphi \right).$$

It gives the boundary conditions for $`j_2`$ and $`j_3`$

$$\begin{array}{c}j_2=j_3=0\text{ at }r=R,\text{ or }\hfill \\ \hfill l=\pm l_m=\pm \left(\mu ^2\mu Ra\right)^{1/2}/R^2.\end{array}$$
(12)
Equation (9) describes the Alfven type oscillations because it contains only the second derivatives over longitudinal coordinate $`l`$ , but equation (10) is that for the magnetosonic waves (also the derivatives over $`a`$ coordinate). In the case of nonaxial symmetric oscillations ($`m0`$) these modes are coupled. We consider the case when the accretion disk excites the MHD waves in the magnetosphere with the dipole magnetic field. Rotating with the Keplerian frequency $`\omega _\mathrm{K}`$ different from the star rotation frequency $`\mathrm{\Omega }`$, the conducting disk deflects the magnetic field lines in equatorial plane to $`\varphi `$-direction. Because this deflection is symmetric over $`\varphi `$ we consider here only the case $`m=0`$. As Alfven and magnetosonic oscillations are independent and Alfven wave has only the $`\varphi `$ component of the disturbed magnetic field

$$j_2=\frac{c}{4\pi \varrho ^2}\frac{\delta B_3}{l},$$

$$\delta B=\delta B_1l+\delta B_2a+\delta B_3\varphi ,\delta B_3=\varrho \delta B_\varphi ,$$

while the magnetosonic has two, $`\delta B_1`$ and $`\delta B_2`$

$$j_3=\frac{cB^2\varrho ^2}{4\pi }\left(\frac{\delta B_2}{l}\frac{\delta B_1}{a}\right).$$
(13)
For this reason we choose only the Alfven wave ($`\delta B_1=\delta B_2=0,\delta B_\varphi 0`$)

$$j_2+\frac{v_A^2}{\omega ^2\varrho ^2B}\frac{}{l}\left[\varrho ^2B^2\frac{}{l}\left(j_2B\right)\right]=0.$$
(14)
The boundary conditions

$$j_2\left(a,l=\pm l_m\right)=0$$

imply the quantization of eigen frequencies of the Alfven modes

$$\omega =\omega _j(a),j=1,2,\mathrm{}$$

After parameterization the eigen frequencies can be presented in the form

$$\omega _j=\omega _0\omega _j^{}(f),$$
(15)
where

$$\omega _0=\left(\frac{\mu }{4\pi R^5\lambda m_\mathrm{p}}\right)^{1/2},$$

is the characteristic Alfven frequency, and $`\omega _j^{^{}}`$ is the dimensionless amplitude depending on the dimensionless parameter $`f=\varrho _{\mathrm{max}}/R`$, which means the height of the dipole magnetic field line in the star radius units. The calculation for the general mode $`j=1`$ gives the dependence of on $`f`$ as pictured on Fig. 2. At large $`f`$ the eigen frequency falls down like $`\omega _\mathrm{A}^{^{}}f^{5/2}`$. The reason of this behavior is that the main contribution to the frequency gives the wave structure in the top of magnetic field tube. Indeed, multiplying the (14) on to $`\left(B^2\varrho ^2\omega ^2/v_A^2\right)j_2^{}`$ and integrating it over $`l`$ from$`l_m`$ to $`+l_m`$, we get

$$\omega ^2=\frac{_{l_m}^{l_m}\varrho ^2\left|B\frac{}{l}\left(Bj_2\right)\right|^2dl}{_{l_m}^{l_m}\frac{\varrho ^2}{v_A^2}\left|\left(Bj_2\right)\right|^2dl}.$$
(16)
From (16) it is seen that $`\omega \left(v_\mathrm{A}\right)_{\mathrm{min}}/L`$, where $`L`$ is the length of the magnetic field line. Because of $`v_A=\left(B/4\pi \lambda m_\mathrm{p}\right)^{1/2}f^{3/2}`$ and $`Lf`$ we obtain the calculated dependence

$$\omega f^{5/2}(f1).$$

As for the values of $`f`$ near the star surface $`|f1|1`$, then $`L(f1)^{1/2}`$ and $`\omega (f1)^{1/2}(f1)`$ as imagined on the Fig. 2. For the same reason the more high modes have the frequencies very close to the harmonics of the general mode

$$\omega _jj\omega _1,\omega _\mathrm{A}\omega _1.$$
(17)
For the plasma density of $`n_\mathrm{p}810^{17}`$ cm-3 near the star surface and magnetic momentum $`\mu 10^{27}`$ G$``$cm3 ($`B10^9`$ G) the frequency $`\omega _\mathrm{A}`$ is of the order of $`2\pi `$ kHz at the distance $`r5.6R`$. It means that QPO of 1 kHz can appear as the result of the resonance of the disk rotation and Alfven eigen modes excited by the disk

$$i(\omega _\mathrm{K}\mathrm{\Omega })=\omega _jj\omega _\mathrm{A}.$$
(18)
Here $`i`$ is the integer number. The resonance condition (18) defines the magnetic field tubes, where the resonance occur, $`f_{i,j}`$ :

$$\omega _\mathrm{K}(f_{i,j})\mathrm{\Omega }=\frac{j}{i}\omega _\mathrm{A}(f_{i,j}).$$
(19)
The solution of (18) is pictured on the Fig. 3. In general case we have two solutions $`f_1`$ and $`f_2`$ ($`f_2>f_1`$ ). Their positions are not constants in time because of the changing of the accretion rate and as a consequence of changing of the plasma density $`n_\mathrm{p}`$ . The frequency of the Alfven eigen modes $`\omega _jn_\mathrm{p}^{1/2}`$. It means that higher accretion rate (or $`n_\mathrm{p}`$) larger $`f_1`$ and smaller $`f_2`$. Accordingly $`\omega _\mathrm{K}(f_1)`$ becomes larger, but $`\omega _\mathrm{K}(f_2)`$ becomes smaller. The observed changes of QPO frequency correspond to the behavior of the root $`f_1`$ . Apart, it is observed the clear proportionality between intensity of X-ray emission and the frequency of QPO of 1 kHz band. We consider that the resonance at $`f=f_1`$ is more strong than at $`f=f_2`$ because of

* it is more close to the star surface where the dipole magnetic field is large; and
* usually $`f_2`$ is large enough, so that $`\omega _\mathrm{K}(f_2)`$ falls down from the 1 kHz frequencies.

Further we will discuss only the high frequency solution $`f=f_1`$ . How the QPO frequencies are connected with discussed resonance? At first the excited Alfven oscillations will destroy the accretion disk. The magnetic field gives the drift velocity $`\delta v_z`$ for the disk plasma particles

$$\delta v_z=c\frac{\left[E_r\delta B_\varphi \right]}{B^2};\delta v_z=v_\mathrm{K}\frac{\delta B_\varphi }{B},$$
(20)
(see $`E_r`$ from (2)). The disk matter will be pushed from equatorial plane in the direction along the dipole magnetic line to the star surface. And the accretion flow will be modulated by the frequency $`\omega _j(f_1)+\mathrm{\Omega }`$. The point $`f=f_1`$ will be factually the inner edge of the accretion disk. In observed X-ray emission we can see several characteristic frequencies:

1. $`\omega _\mathrm{K}(f_1)`$, it is the modulation of the X-ray flux by the inhomogeneity of the inner edge of the disk;
2. $`\omega _\mathrm{K}(f_1)\mathrm{\Omega }`$, it is the beating frequency corresponding to the scattering of X-rays from the sources on the star surface, rotating with $`\mathrm{\Omega }`$ , by the inner edge of the disk;
3. $`\omega _j(f_1)+\mathrm{\Omega }=i\omega _\mathrm{K}(f_1)+(1i)\mathrm{\Omega }`$, it is the frequency of the accretion flow modulated by Alfven oscillations;
4. different combinations of frequencies 1, 2, 3.

## 3 Sonic oscillation

It is naturally to consider that the plasma flow along the tube $`f=f_1`$ is not stationary. The sound mode can be excited.

Let us consider plasma as isothermic with temperature $`T`$. For disturbance of concentration $`\delta n`$ along the magnetic tube we receive from hydrodynamics

$$\delta n+\frac{T}{m_\mathrm{p}\omega ^2}\mathrm{div}\mathrm{grad}\delta n=0.$$
(21)
In dipole coordinates it looks like

$$\delta n+\frac{T}{m_\mathrm{p}\omega ^2}B^2\frac{^2}{l^2}\delta n=0.$$
(22)
The boundary condition is on the star surface. The plasma velocity component which is normal to the surface must be equal to zero.

$$v=i\frac{T}{m_\mathrm{p}n_\mathrm{p}\omega }\delta n|_{l=\pm l_m}=0.$$
(23)
The first harmonics corresponding to the solution of this equation with boundary condition can be presented as

$$\omega (f)=\omega _\mathrm{s}\omega _\mathrm{s}^{}(f),$$
(24)
where $`\omega _\mathrm{s}^2=T/m_\mathrm{p}R^2`$ and $`\omega _\mathrm{s}^{}(f)`$ is dimensionless function on $`f`$, which is presented on the Fig. 4, $`\omega _\mathrm{s}^{}(f)`$ is of the order of $`1`$ when $`f2`$. In the region of $`f1`$ the function $`\omega _\mathrm{s}^{}(f)`$ can be approximated by the expression $`2/(\pi \sqrt{f1})`$. So, $`\omega _0`$ is the charachteristic frequency for sonic oscillation. If we take $`T=1`$ keV and $`R=10`$ km then $`\omega _0=30`$ s-1 or $`\nu _0=4.9`$ Hz. Thus, sonic oscillation can be found on $`10`$ Hz band that is observed simaltaneously with the QPO of $`1`$ kHz band.

## 4 Lense-Thirring precession

While the accreting disk is spinning at the distance $`r`$ about the neutron star with Keplerian frequency, it precesses about its equilibrium position with Lense-Thirring frequency due to Kerr metrics (Thorne et al. 1986; Morsink & Stella 1999)

$$\omega _{\mathrm{LT}}=a\frac{r_\mathrm{g}}{r^3+\left(a/c\right)^2(r_\mathrm{g}+r)},$$
(25)
where $`a=J/M`$ ($`J=I\mathrm{\Omega }`$ is the star’s angular momentum and $`I=2MR^2/7`$ is the momentum if inertia of the star), $`r_\mathrm{g}=2GM/c^2`$ is the gravitational radius of the neutron star of mass $`M`$ and radius $`R`$.

In the system of units where $`M`$ is measured in masses of sun and $`R`$ in $`10`$ km the Lense-Thirring frequency can be presented in the form of

$$\omega _{\mathrm{LT}}=8.4410^2\frac{M\mathrm{\Omega }}{R}f^3.$$
(26)
For $`M/R1`$, $`f1`$ and $`\mathrm{\Omega }2\pi 300`$ s-1 $`\omega `$ is of the order of $`2\pi 25`$ s-1. The second harmonics of Lense-Thirring precession is usually observed (Stella & Vietri 1999). In our model it can be illustrated by the fact that some mass is thrown out the disk each time it is in extreme position, i.e. with frequency $`2\omega _{\mathrm{LT}}`$ (see Fig. 5).

## 5 Scorpius X-1

We have applied our speculations for Sco X-1 (van der Klis et al. 1996). There were observed 4 QPOs: $`6`$$`20`$ Hz normal/flaring quasi-periodic oscillations; $`45`$ Hz QPOs and two kHz QPOs. According to our model we ascribed the $`6`$$`20`$ Hz QPOs to the sound oscillations, $`45`$ Hz QPOs to Lense-Thirring oscillations, and two kHz QPOs to Alfven and Keplerian oscillations correspondingly.

Below we use the system of units where temperature is measured in eV, masses of objects – in masses of sun and distances – in $`10`$ km. Thus, we can write

for sonic: $`\omega _\mathrm{s}(f)`$ $`=0.98\sqrt{{\displaystyle \frac{T}{R^2}}}\omega _\mathrm{s}^{}(f),`$ (27)
for Keplerian: $`\omega _\mathrm{K}(f)`$ $`=11.510^3\sqrt{{\displaystyle \frac{M}{R^3}}}f^{3/2},`$ (28)
where $`f=r/R`$ – a parameter that characterizes the distance between the disk and neutron star surface.

For every pair of data $`(\omega _\mathrm{K},\omega _\mathrm{s})`$ we can build a curve in $`({\displaystyle \frac{M}{R^3}};{\displaystyle \frac{T}{R^2}})`$ space. This results in the graph on Fig. 6.

The axis of abscissae is the space for possible meanings of the ratio of $`M/R^3`$. For each of possible meanings of this ratio we can find the temperature $`T`$ (or, to be more accurate, the ratio $`T/R^2`$) for each experiment on determining the frequencies $`\omega _\mathrm{K}`$ and $`\omega _\mathrm{s}`$. So we can evaluate the range of the change of temperature $`T`$. One can find it remarkable that all the curves cross each other in almost one point! So we can suppose that the temperature of the disk is approximately constant. From the graph we determine $`M/R^3=0.374\pm 0.004`$ and $`T/R^2=16.8\pm 0.8`$.

Méndez & van der Klis (1999) have observed that the separation between the two high-frequency QPO peaks changes from one observation to another. In our work we considered the equation

$$i(\omega _\mathrm{K}\mathrm{\Omega })=\omega _jj\omega _\mathrm{A}.$$
(29)
In general case the peak separation will be

$$\omega _\mathrm{K}\omega _\mathrm{A}=\frac{j}{i}\omega _\mathrm{A}+\mathrm{\Omega }\omega _\mathrm{A}=\left(\frac{j}{i}1\right)\omega _\mathrm{A}+\mathrm{\Omega }.$$
(30)
It is easy to see, that when $`ij`$ the peak separation will not remain constant. For example, if $`i=2`$ and $`j=1`$ then $`\omega _\mathrm{K}\omega _\mathrm{A}=\omega _\mathrm{A}/2+\mathrm{\Omega }`$. If $`i=1`$ and $`j=2`$, then $`\omega _\mathrm{K}\omega _\mathrm{A}=\omega _\mathrm{A}+\mathrm{\Omega }`$. Here we must mention that Alfven eigen frequency depends on the parameter $`\lambda `$ that changes when $`f`$ will change and, therefore, when $`f`$ changes, numbers $`i`$ and $`j`$ can change. As a result, the dependence of the peak separation on $`\omega _\mathrm{K}`$ will be a piece linear function. The dependencies of separation frequency on the lower one for different objects is shown on the Fig. 7. For Sco X-1 the dependence can be well approximated by the linear function (30) with $`j/i=3/4`$ and $`\mathrm{\Omega }=2\pi 450`$ s-1.

So, comparing the $`45`$ Hz QPO with Lense-Thirring precession we determine the ratio $`M/R0.60.`$ Finally $`M=0.75M_{}`$, $`R=12.6`$ km, $`T=27`$ eV. The point $`(M,R)`$ on the gravitational mass – radius diagram is shown on the Fig. 8.

For this $`M,R,T`$ we can find that $`\omega _{0\mathrm{s}}=4.16`$ s-1 or $`\nu _{0\mathrm{s}}=0.66`$ Hz. This low value can be explained by the fact that disk approaches too close to the neutron star (down to $`f=1.01`$) and in this region of $`f`$ function $`\omega _\mathrm{s}^{}(f)`$ takes values much greater than $`1`$ (e.g. $`\omega _\mathrm{s}^{}(1.01)=15.7`$).

At first sight, $`T`$ must be equal to the radiation temperature of the neutron star, which can be determined from the X-ray color-color diagram. We put it on one graph from Scorpius X-1 and from black body (see Fig. 9). The normal branch of that from Scorpius almost coincide with that of black body. And the temperature of the black body varies from $`1.7`$ keV up to $`1.95`$ keV. But if we estimate the free path of the photons from the neutron star $`l1/\sigma n_\mathrm{p}20`$ km we find that it is greater then $`R`$ and so, radiation does not heat the inner part of the disk.

###### Acknowledgements.

This work is supported by Russian Fundamental Research Foundation Grant No 99-02-17184.


text:
# Search in Power-Law Networks

## 1 Introduction

A number of large distributed systems, ranging from social  to communication  to biological networks  display a power-law distribution in their node degree. This distribution reflects the existence of a few nodes with very high degree and many with low degree, a feature not found in standard random graphs . An illustration of the power-law nature of such networks is given by the AT&T call graph. A call graph is a graph representation of telephone traffic on a given day in which nodes represent people and links the phone calls among them.

As is shown in Figure 1, the out-link degree distribution for a massive graph of telephone calls between individuals is power-law, with an exponent of approximately 2.1. The same distribution is obtained for the case of in-links. This power-law in the link distribution reflects the presence of central individuals who interact with many others on a daily basis and play a key role in relaying information.

While recent work has concentrated on the properties of these power-law networks and how they are dynamically generated , there remains the interesting problem of finding efficient algorithms for searching within these particular kinds of graphs. Recently, Kleinberg  studied search algorithms in a graph where nodes are placed on a 2-D lattice and each node has a fixed number of links whose placement is correlated with lattice distance to the other nodes. Under a specific form of the correlation, an algorithm with knowledge of the target’s location can find the target in polylogarithmic time.

In the most general distributed search context however, one may have very little information about the location of the target. Increasingly a number of pervasive electronic networks, both wired and wireless, make geographic location less relevant. A particularly interesting example is provided by the recent emergence of peer-to-peer networks, which have gained enormous popularity with users wanting to share their computer files. In such networks, the name of the target file may be known, but due to the network’s ad hoc nature, until a real-time search is performed the node holding the file is not known. In contrast to the scenario considered by Kleinberg, there is no global information about the position of the target, and hence it is not possible to determine whether a step is a move towards or away from the target. One simple way to locate files, implemented by Napster, is to use a central server which contains an index of all the files every node is sharing as they join the network. This is the equivalent of having a giant white pages for the entire United States. Such directories now exist online, and have in a sense reduced the need to find people by passing messages. But for various reasons, including privacy and copyright issues, in a peer-to-peer network it is not always desirable to have a central server.

File-sharing systems which do not have a central server include Gnutella and Freenet. Files are found by forwarding queries to one’s neighbors until the target is found. Recent measurements of Gnutella networks  and simulated Freenet networks  show that they have power-law degree distributions. In this paper, we propose a number of message-passing algorithms that can be efficiently used to search through power law networks such as Gnutella. Like the networks that they are designed for, these algorithms are completely decentralized and exploit the power law link distribution in the node degree. The algorithms use local information such as the identities and connectedness of their neighbors, and their neighbors’ neighbors, but not the target’s global position. We demonstrate that our search algorithms work well on real Gnutella networks, scale sub-linearly with the number of nodes, and may help reduce the network search traffic that tends to cripple such networks.

The paper is organized as follows. In Section 2, we present analytical results on message-passing in power-law graphs, followed by simulation results in Section 3. Section 4 compares the results with Poisson random graphs. In Section 5 we consider the application of our algorithms to Gnutella, and Section 6 concludes.

## 2 Search in Power-Law Random Graphs

In this section we use the generating function formalism introduced by Newman  for graphs with arbitrary degree distributions to analytically characterize serach-cost scaling in power-law graphs.

### 2.1 Random Walk Search

Let $`G_0(x)`$ be the generating function for the distribution of the vertex degrees $`k`$. Then

$$G_0(x)=\underset{0}{\overset{\mathrm{}}{}}p_kx^k$$
(1)
where $`p_k`$ is the probability that a randomly chosen vertex on the graph has degree $`k`$.

For a graph with a power-law distribution with exponent $`\tau `$, minimum degree $`k=1`$ and an abrupt cutoff at $`m=k_{max}`$, the generating function is given by

$$G_0(x)=c\underset{1}{\overset{m}{}}k^\tau x^k$$
(2)
with $`c`$ a normalization constant which depends on $`m`$ and $`\tau `$ to satisfy the normalization requirement

$$G_0(1)=c\underset{1}{\overset{m}{}}k^\tau =1$$
(3)
The average degree of a randomly chosen vertex is given by

$$<k>=\underset{1}{\overset{m}{}}kp_k=G_0^{}(1)$$
(4)
Note that the average degree of a vertex chosen at random and one arrived at by following a random edge are different. A random edge arrives at a vertex with probability proportional to the degree of the vertex, i.e. $`p^{}(k)kp_k`$. The correctly normalized distribution is given by

$$\frac{_kkp_kx^k}{_kkp_k}=x\frac{G_0^{}(x)}{G_0^{}(1)}$$
(5)
If we want to consider the number of outgoing edges from the vertex we arrived at, but not include the edge we just came on, we need to divide by one power of $`x`$. Hence the number of new neighbors encountered on each step of a random walk is given by the generating function

$$G_1(x)=\frac{G_0^{}(x)}{G_0^{}(1)}$$
(6)
where $`G_0^{}(1)`$ is the average degree of a randomly chosen vertex as mentioned previously.

In real social networks it is reasonable that one one would have at least some knowledge of one’s friends’ friends. Hence, we now compute the distribution of second degree neighbors. The probability that any of the 2nd degree neighbors connect to any of the first degree neighbors or to one another goes as $`N^1`$ and can be ignored in the limit of large $`N`$. Therefore, the distribution of the second neighbors of the original randomly chosen vertex is determined by

$$\underset{k}{}p_k[G_1(x)]^k=G_0(G_1(x))$$
(7)
It follows that the average number of second degree neighbors is given by

$$z_{2A}=[\frac{}{x}G_0(G_1(x))]_{x=1}=G_0^{}(1)G_1^{}(1)$$
(8)
Similarly, if the original vertex was not chosen at random, but arrived at by following a random edge, then the number of second degree neighbors would be given by

$$z_{2B}=[\frac{}{x}G_1(G_1(x))]_{x=1}=[G_1^{}(1)]^2$$
(9)
In both Equation 8 and Equation 9 the fact that $`G_1(1)=1`$ was used.

Both these expressions depend on the values $`G_0^{}(1)`$ and $`G_1^{}(1)`$ so let’s calculate those for given $`\tau `$ and $`m`$. For simplicity and relevance to most real-world networks of interest we assume $`2<\tau <3`$.

$$G_0^{}(1)=\underset{1}{\overset{m}{}}ck^{1\tau }_1^mx^{\tau 1}𝑑x=\frac{1}{\tau 2}(1m^{2\tau })$$
(10)
$`G_1^{}(1)`$ $`=`$ $`{\displaystyle \frac{1}{G_0^{}(1)}}{\displaystyle \frac{}{x}}{\displaystyle \underset{1}{\overset{m}{}}}ck^{1\tau }x^{k1}`$ (11)
$`=`$ $`{\displaystyle \frac{1}{G_0^{}(1)}}{\displaystyle \underset{2}{\overset{m}{}}}ck^{1\tau }(k1)x^{k2}`$ (13)
$``$ $`{\displaystyle \frac{1}{G_0^{}(1)}}({\displaystyle \frac{m^{3\tau }(\tau 2)2^{2\tau }(\tau 1)+m^{2\tau }(3\tau )}{(\tau 2)(3\tau )}}`$ (14)
for large cutoff values $`m`$. Now we impose the cutoff of Aiello et al.  at $`mN^{1/\tau }`$. Since $`m`$ scales with the size of the graph $`N`$ and for $`2<\tau <3`$ the exponent $`2\tau `$ is negative, we can neglect terms constant in $`m`$. This leaves

$$G_1^{}(1)=\frac{1}{G_0^{}(1)}\frac{m^{3\tau }}{(3\tau )}$$
(15)
Substituting into Equation 8 (the starting node is chosen at random) we obtain

$$z_{2A}=G_0^{}(1)G_1^{}(1)m^{3\tau }$$
(16)
We can also derive $`z_{2B}`$, the number of 2nd degree neighbors encountered as one is doing a random walk on the graph.

$$z_{2B}=[G_1^{}(1)]^2=[\frac{\tau 2}{1m^{2\tau }}\frac{m^{3\tau }}{3\tau }]^2$$
(17)
Letting $`mN^{1/\tau }`$ as above, we obtain

$$z_{2B}N^{2(\frac{3}{\tau }1)}$$
(18)
and the scaling of the number of steps required becomes

$$sN^{3(12/\tau )}$$
(19)
In the limit $`\tau 2`$, equation 17 becomes

$$z_{2B}\frac{N}{ln^2(N)}$$
(20)
and the scaling of the number of steps required is

$$sln^2(N)$$
(21)
### 2.2 Search utilizing high degree nodes

Random walks in power-law networks naturally gravitate towards the high degree nodes, but an even better scaling is achieved by intentionally choosing high degree nodes. For $`\tau `$ sufficiently close to $`2`$ one can walk down the degree sequence, visiting the node with the highest degree, followed by a node of the next highest degree, etc. Let $`ma`$ be the degree of the last node we need to visit in order to scan a certain fraction of the graph. We make the self-consistent assumption that $`a<<`$m, i.e. the degree of the node hasn’t dropped too much by the time we have scanned a fraction of the graph. Then the number of first degree neighbors scanned is given by

$$z_{1D}=_{ma}^mNk^{1\tau }𝑑kNam^{1\tau }$$
(22)
The number of nodes having degree between $`ma`$ and $`m`$, or equivalently, the number of steps taken is given by $`_{ma}^mk^\tau a`$. The number of second degree neighbors when one follows the degree sequence is given by:

$$z_{1D}G_1^{}(1)Nam^{2(2\tau )}$$
(23)
which gives the number of steps required as

$$sm^{2(\tau 2)}N^{2\frac{4}{\tau }}$$
(24)
We now consider when and why it is possible to go down the degree sequence. We start with the fact that the original degree distribution is a power-law:

$$p(x)=(\underset{1}{\overset{m}{}}x^\tau )^1x^\tau $$
(25)
where $`m=N^{1/\tau }`$ is the maximum degree. A node chosen by following a random link in the graph will have its remaining outgoing edges distributed according to

$$p^{}(x)=(\underset{0}{\overset{m1}{}}(x+1)^{(1\tau )})^1(x+1)^{(1\tau )}$$
(26)
At each step one can choose the highest degree node among the $`n`$ neighbors. The expected number of the outgoing edges of that node can be computed as follows. In general, the cumulative distribution (CDF) $`P_{max}(x,n)`$ of the maximum of n random variables can be expressed in terms of the CDF $`P(x)=_0^xp(x)\mathrm{d}x`$ of those random variables: $`P_{max}(x,n)=P(x)^n`$. This yields

$$p_{max}^{}(x,n)=n(1+x)^{1\tau }(\tau 2)(1(x+1)^{2\tau })^{n1}(1N^{2/\tau 1})^n$$
(27)
for the distribution of the number of links the richest neighbor among n neighbors has.

Finally, the expected degree of the richest node among n is given by

$$E[x_{max}(n)]=\underset{0}{\overset{m1}{}}xp_{max}^{}(x,n)$$
(28)
We numerically integrated the above equation to derive the ratio between the degree of a node and the expected degree of its richest neighbor. The ratio is plotted in Figure 2. For a range of exponents and node degrees, the expected degree of the richest neighbor is higher than the node itself. However, eventually (the precise point depends strongly on the power-law exponent), the probability of finding an even higher degree node in a neighborhood of a very high degree node starts falling.

What this means is that one can approximately follow the degree sequence across the entire graph for a sufficiently small graph or one with a power-law exponent close to 2 ($`2.0<\tau <2.3`$). At each step one chooses a node with degree higher than the current node, quickly finding the highest degree node. Once the highest degree node has been visited, it will be avoided, and a node of approximately second highest degree will be chosen. Effectively, after a short initial climb, one goes down the degree sequence. This is the most efficient way to do this kind of sequential search, visiting highest degree nodes in sequence.

## 3 Simulation

We used simulations of a random social network with a power-law link distribution of $`\tau =2.1`$ to validate our analytical results. As in the analysis above, a simple cutoff at $`mN^{1/\tau }`$ was imposed. The expected number of nodes among N having exactly the cutoff degree is 1. No nodes of degree higher than the cutoff are added to the graph. In real world graphs one of course does observe nodes of degree higher than this imposed cutoff, so that our simulations become a worse case scenario. Once the graph is generated, the largest connected component (LCC) is extracted, that is the largest subset of nodes such that any node can be reached from any other node. For $`2<\tau <3.48`$ a giant connected component exists , and all our measurements are performed on the LCC. We observe that the LCC contains the majority of the nodes of the original graph and most of the links as well. The link distribution of the LCC is nearly identical to that of the original graph with a slightly smaller number of 1 and 2 degree nodes.

Next we apply our message passing algorithm to the network. Two nodes, the source and the target, are selected at random. At each time step the node which has the message passes it on to one of its neighbors. The process ends when the message is passed on to a neighbor of the target. Since each node knows the identity of all of its neighbors, it can pass the message directly to the target if the target happens to be one of it’s neighbors. The process is analogous to performing a random walk on a graph, where each node is ’visited’ as it receives the message.

There are several variants of the algorithm, depending on the strategy and the amount of local information available.

1. The node can pass the message on to one of its neighbors at random, or it can avoid passing it on to a node which has already seen the message.
2. If the node knows how many neighbors it has, it can choose to pass the message onto the neighbor with the most neighbors.
3. The node may know only its neighbors, or it may know who its neighbors’ neighbors are. In the latter case it would pass the message onto a neighbor of the target.

In order to avoid passing the message to a node that has already seen the message, the message itself must be signed by the nodes as they receive the message. Further, if a node has passed the message, and finds that all of its neighbors are already on the list, it puts a special mark next to its name, which means that it is unable to pass the message onto any new node. This is equivalent to marking nodes as follows:

Node hasn’t been visited.

Node has been visited, but all its neighbors have not.

Node and all its neighbors have been visited already.

Here we compare two strategies. The first performs a random walk, where only retracing the last step is disallowed. In the message passing scenario, this means that if Bob just received a message from Jane, he wouldn’t return the message to Jane if he could pass it to someone else. The second strategy is a self avoiding walk which prefers high degree nodes to low degree ones. In each case both the first and second degree neighbors are scanned at each step.

Figure 4 shows the scaling of the average search time with the size of the graph for the two strategies. The scaling (exponent 0.79 for the random walk and 0.70 for the high degree strategy) is not as favorable as in the analytic results derived above ($`0.14`$ for the random walk and 0.1 for the high degree strategy when $`\tau =2.1`$) .

Consider, on the other hand, the number of steps it takes to cover half the graph. For this measure we observe a scaling which is much closer to the ideal. As shown in Figure 4, the cover time scales as $`N^{0.37}`$ for the random walk strategy vs. $`N^{0.15}`$ from Equation 19. Similarly, the high degree strategy cover time scales as $`N^{0.24}`$ vs. $`N^{0.1}`$ in Equation 24.

The difference in the value of the scaling exponents of the cover time and average search time implies that a majority of nodes can be found very efficiently, but others demand high search costs. Figures 6 and 6 show just that. A large portion of the 10,000 node graph is covered within the first few steps, but some nodes take as many steps or more to find as there are nodes total. Take for example the high degree seeking strategy. About 50% of the nodes are scanned within the first 10 steps (meaning that it would take about $`10+2=12`$ hops to reach 50% of the graph). However, the skewness of the search time distribution bring the average number of steps needed to 217.

Figure 8 illuminates why some nodes take such a long time to find. It shows the color of nodes visited in a random walk on a $`N=1000`$ node power-law graph with exponent 2.1 and an abrupt cutoff at $`N^{1/2.1}`$. The number of nodes of each color encountered in 50 step segments is recorded in the bar for that time period.

The brief initial period where many white nodes are visited is followed by a long period dominated by gray and black nodes. The random walk starts retracing its steps, only occasionally encountering a new node it hasn’t seen. Close to the end of the walk, even the gray nodes have been exhausted as the random walk repeatedly visits back nodes before it finds the final few nodes.

A self-avoiding high-degree seeking random walk is an improvement over the random walk described above, but still cannot avoid retracing its steps. Figure 8 shows the corresponding bar chart in 50 step increments. Comparing figures 8 and 8 we observe that the self-avoiding strategy is somewhat effective, with the total number of steps needed to cover the graph about 13 times smaller, and the fraction of visits to gray and black nodes significantly reduced. Figure 9 shows the degree in addition to the color of the node visited at every step. It demonstrates that one is able to follow the degree sequence initially, going from high degree white nodes to low degree white nodes.

Although the strategy of revisiting nodes modifies the scaling behavior, it is not the process of revisiting itself which necessarily changes the scaling. If nodes were uniformly linked, at every step the number of new nodes seen would be proportional to the number of unexplored nodes in the graph. The factor by which the search is slowed down through revisits would be independent of the size of the graph. Hence, revisiting alone does not account for the difference in scaling.

The reason why the simulated scaling exponents for these search algorithms do not follow the ideal is the same reason why power-law graphs are so well suited to search: the link distribution is extremely uneven. A large number of links point to only a small subset of high degree nodes. When a new node is visited, its links do not let us uniformly sample the graph, they preferentially lead to high degree nodes, which have likely been seen or visited in a previous step. This would not be true of a Poisson graph, where all the links are randomly distributed and hence all nodes have approximately the same degree. We will explore and contrast the search algoritm on a Poisson graph in the following section.

## 4 Comparison with Poisson Distributed Graphs

In a Poisson random graph with $`N`$ nodes and $`z`$ edges, the probability $`p=z/N`$ of an edge between any two nodes is the same for all nodes. The generating function $`G_0(x)`$ is given by :

$$G_0(x)=e^{z(x1)}$$
(29)
In this special case $`G_0(x)=G_1(x)`$, so that the distribution of outgoing edges of a node is the same whether one arrives at the vertex by following a link or picks the node at random. This makes analysis search in a Poisson random graph particularly simple. The expected number of new links encountered at each step is a constant $`p`$. So that the number of steps needed to cover a fraction $`c`$ of the graph is $`s=cN/p`$. If $`p`$ remains constant as the size of the graph increases, the cover time scales linearly with the size of the graph. This has been verified via simulation as shown in Figure 10.

In our simulations the probability $`p`$ grows slowly towards its asymptotic value as the size of the graph is increased because of the particular choice of cutoff at $`mN^{(1/\tau )}`$ for the power-law link distribution. We generated Poisson graphs with the same number of nodes and links for comparison. Within this range of graph sizes, growth in the average number of links per node appears as $`N^{0.6}`$, making the average number of 2nd degree neighbors scale as $`N^{0.15}`$. This means that the scaling of the cover time scales as $`N^{0.85}`$, as shown in Figure 10. Note how well the simulation results match the analytical expression. This is because nodes can be approximately sampled in an even fashion by following links.

The reason why the cover time for the Poisson graph matches the analytical prediction and the power-law graph does not is illustrated in Figure 11. If links were approximately evenly distributed among the nodes, then if at one point in the search 50% of the graph has already been visited, one would expect to revisit previously seen nodes about 50% of the time. This is indeed the case for the Poisson graph. However, for the power-law graph, when 50% of the graph has been visited, nodes are revisited about 80% of the time, which implies that the same high degree nodes are being revisted before new low-degree ones. It is this bias which accounts for the discrepancy between the analytic scaling and the simulated results in the power-law case.

However, even the simulated $`N^{0.35}`$ scaling for a random, minimally self-avoiding strategy on the power-law graph out-performs the ideal $`N^{0.85}`$ scaling for the Poisson graph. It’s also important to note that the the high degree node seeking strategy has a much greater success in the power-law graph because it relies heavily on the fact that the number of links per node varies considerably from node to node. To illustrate this point, we executed the high degree seeking strategy on two graphs, Poisson and power-law, with the same number of nodes, and the same exponent $`\tau =2`$. In the Poisson graph, the variance in the number of links was much smaller, making the high degree node seeking strategy comparatively ineffective as shown in Figure 12.

In the power-law graph we can start from a randomly chosen node. In this case the starting node has only one link, but two steps later we find ourselves at a node with the highest degree. From there, one approximately follows the degree sequence, that is, the node richest in links, followed by the second richest node, etc. The strategy has allowed us to scan the maximum number of nodes in the minimum number of steps. In comparison, the maximum degree node of the exponential graph is 11, and it is reached only on the 81st step. Even though the two graphs have a comparable number of nodes and edges, the exponential graph does not lend itself to quick search.

## 5 Gnutella

Gnutella is a peer-to-peer filesharing system which treats all client nodes as functionally equivalent and lacks a central server which can store file location information. This is advantageous because it presents no central point of failure. The obvious disadvantage is that the location of files is unknown. When a user wants to download a file, she sends a query to all the nodes within a neighborhood of size $`ttl`$, the time to live assigned to the query. Every node passes on the query to all of its neighbors and decrements the $`ttl`$ by one. In this way, all nodes within a given radius of the requesting node will be queried for the file, and those who have matching files will send back positive answers.

This broadcast method will find the target file quickly, given that it is located within a radius of $`ttl`$. However, broadcasting is extremely costly in terms of bandwidth. Every node must process queries of all the nodes within a given $`ttl`$ radius. In essence, if one wants to query a constant fraction of the network, say 50%, as the network grows, each node and network edge will be handling query traffic which is proportional to the total number of nodes in the network.

Such a search strategy does not scale well. As query traffic increases linearly with the size of Gnutella graph, nodes become overloaded as was shown in a recent study by Clip2 . 56k modems are unable to handle more than 20 queries a second, a threshold easily exceeded by a network of about 1,000 nodes. With the 56k nodes failing, the network becomes fragmented, allowing users to query only small section of the network.

The search algorithms described in the previous sections may help ameliorate this problem. Instead of broadcasting a query to a large fraction of the network, a query is only passed onto one node at each step. The search algorithms are likely to be effective because the Gnutella network has a power-law connectivity distribution as shown in Figure 13.

Typically, a Gnutella client wishing to join the network must find the IP address of an initial node to connect to. Currently, ad hoc lists of ”good” Gnutella clients exist . It is reasonable to suppose that this ad hoc method of growth would bias new nodes to connect preferentially to nodes which are already fairly well-connected, since these nodes are more likely to be ”well-known”. Based on models of graph growth  where the ”rich get richer”, the power-law connectivity of ad hoc peer-to-peer networks may be a fairly general topological feature.

By passing the query to every single node in the network, the Gnutella algorithm fails to take advantage of its connectivity distribution. To implement our algorithm the Gnutella clients must be modified to keep lists of the files stored by their first and second degree neighbors have<sup>1</sup><sup>1</sup>1This idea has already been implemented by Clip2 in a limited way. 56k modem nodes attach to a high bandwidth Reflector  node which stores the filenames of the 56k nodes and handles queries on their behalf.. This information must be passed at least once when a new node joins the network, and it may be necessary to periodically update the information depending on the typical lifetime of nodes in the network. Instead of passing the query to every node, queries are only passed along to the highest degree nodes. The IP numbers of the nodes already queried are appended to the query, and they are avoided.

The modified algorithm places an additional cost on every node, that of keeping track of the filenames of its neighbors’ files. Since network connections saturated by query traffic are a major weakness in Gnutella, and since computational and storage resources are likely to remain much less expensive than bandwidth, such a tradeoff is readily made. However, now instead of every node having to handle every query, queries are routed only through high connectivity nodes. Since nodes can select the number of connections that they allow, high degree nodes are presumably high bandwidth nodes that can handle the query traffic. The network has in effect created local directories valid within a two link radius. It is resilient to attack because of the lack of a central server. As for power-law networks in general , the network is more resilient than random graphs to random node failure, but less resilient to attacks on the high degree nodes.

Figure 14 shows the success of the high degree seeking algorithm on the Gnutella network. We simulated the search algorithm on a crawl by Clip2 of the actual Gnutella network of approximately 700 nodes. Assuming that every file is stored on only one node, 50% of the files can be found in 8 steps or less. Furthermore, if the file one is seeking is present on multiple nodes, the search will be even faster.

To summarize, the power-law nature of the Gnutella graph means that these search algorithms can be effective. As the number of nodes increases, the (already small) number of nodes that will need to be queried increases sub-linearly. As long as the high degree nodes are able to carry the traffic, the Gnutella network’s performance and scalability may improve by using these search strategies.

## 6 Conclusion

In this paper we have shown that local search strategies in power-law graphs have search costs which scale sub-linearly with the size of the graph, a fact that makes them very appealing when dealing with large networks. The most favorable scaling was obtained by using strategies which preferentially utilize the high connectivity nodes in these power-law networks. We also established the utility of these strategies for searching on the Gnutella peer-to-peer network.

It may not be coincidental that several large networks are structured in a way that naturally facilitates search. Rather, we find it likely that these networks could have evolved to facilitate search and information distribution. Networks where locating and distributing information, without perfect global information, plays a vital role tend to be power-law with exponents favorable to local search.

For example, large social networks, such as the AT&T call graph and the collaboration graph of film actors, have exponents in the range ($`\tau =2.12.3`$) which according to our analysis makes them especially suitable for searching using our simple, local algorithms. Being able to reach remote nodes by following intermediate links allows communication systems and people to get to the resources they need and distribute information within these informal networks. At the social level, our analysis supports the hypothesis that highly connected individuals do a great deal to improve the effectiveness of social networks in terms of access to relevant resources .

Furthermore, it has been shown that the Internet backbone has a power-law distribution with exponent values between 2.15 and 2.2 , and web page hyperlinks have an exponent of 2.1 . While in the Internet there are other strategies for finding nodes, such as routing tables and search engines, one observes that our proposed strategy is partially used in these systems as well. Packets are routed through highly connected nodes, and users searching for information on the Web turn to highly connected nodes, such as directories and search engines, which can bring them to their desired destinations.

On the other hand, a system such as the power grid of the western United States, which does not serve as a message passing network, has an exponent $`\tau 4`$ . It would be fairly difficult to pass messages in such a network without knowing the target’s location.

Acknowledgement

We would like to thank Clip2 for the use of their Gnutella crawl data.