text:
# Radiation drag effects on magnetically-dominated outflows around compact objects

## 1 Introduction

Magnetohydrodynamic (MHD) models are now developed intensively in theories of the magnetospheres of rotating supermassive black holes ($`M10^8\text{}10^9M_{}`$, $`B_010^4`$ G), which are believed to reside in central engines of Active Galactic Nuclei (AGNs) and quasars (Begelman, Blandford & Rees 1984; Thorne, Price & Macdonald 1986). In particular, it is the MHD model that is the most promising in the problem of the origin and stability of jets. Indeed, the MHD approach explains both the energy source (the rotational energy of the compact object) and the mechanism of the energy and angular momentum loss (for an overview, see e.g. Blandford 2002). Observational evidence in favor of MHD models was recently found in the possible presence of toroidal magnetic fields in jets (Gabuzda et al 1992; Gabuzda et al 1999). Magnetically-dominated outflows are also believed to be responsible for the energy transport in cosmological gamma-ray bursts (Mészáros & Rees 1997; Lee, Wijers & Brown 1998; van Putten & Levinson 2003), when energy is released in the merging of black holes or neutron stars ($`MM_{}`$, $`B_010^{15}`$ G).

It has been suggested that the density of photons in the vicinity of the central engine is so high that they may drastically change the characteristics of the ideal MHD outflow. For example, they may result in extensive $`e^+e^{}`$ pair creation (Svensson 1984), acceleration of low-energy pairs by the radiation drag force (Phinney 1982; Turolla, Nobili & Calvani 1986; Beloborodov 1999) and deceleration of high-energy particles (Melia & Königl 1989; Sikora et al 1996). In other words, a self-consistent consideration should take the drag force into account.

So far the two processes – the ideal MHD acceleration and the action of external photons – have been considered independently. The first step to combine them was made by Li, Begelman & Chiueh (1992). In particular, it was shown how the equations can be integrated in a conical geometry (which is impossible in the general case). On the other hand, the consideration was performed within the approximation of a fixed poloidal magnetic field. Under this assumption the fast magnetosonic surface of a cold flow is shifted to infinity (Michel 1969; Kennel, Fujimura & Okamoto 1976; Lery et al 1998). As a result, it was impossible to analyze the effects of radiation drag on the position of the fast magnetosonic surface and the properties of the outflow outside this surface.

The main goal of this paper is to determine more carefully the radiation drag effects on a magnetically-dominated outflow. To describe analytically the effects of radiation drag, including simultaneously the disturbance of the magnetic surfaces we consider here a quasi-monopole outflow. For AGNs, such geometry in the immediate vicinity of the central engine was recently confirmed by direct observations (Junior, Biretta & Livio 1999). In other words, in the zeroth approximation (i.e., without drag) we use the analytical solution for a magnetically-dominated MHD outflow (Beskin, Kuznetsova & Rafikov 1998, hereafter Paper I), in which the fast magnetosonic surface is located at a finite distance from the origin.

For simplicity we consider the following model of the radiation field in the vicinity of the central engine. First, we notice that for ultra-relativistic particles the energy of a photon propagating nearly along the particle trajectory remains almost the same after a collision. This means that the drag force from these photons is small. Thus, only the isotropic component of the photon field contributes substantially to the drag force. Hence, in our geometry with a strong central source of photons and a monopole outflow of particles, only a small fraction of photons (the isotropic component of the photon field) interacts efficiently with the particles, producing inverse Compton photons with energies $`_{\mathrm{IC}}\gamma ^2_{\mathrm{ph}}`$.

The isotropic component can be produced, firstly, by the outer part of the accretion disk and, secondly, by external sources. It can be modeled as

$$U=U_{\mathrm{iso}}=U_\mathrm{A}\left(\frac{r}{R_\mathrm{L}}\right)^n+U_{\mathrm{ext}},$$
(1)
where $`R_\mathrm{L}=c/\mathrm{\Omega }`$ is the radius of the light cylinder, $`U_\mathrm{A}=U(R_\mathrm{L})=L_{\mathrm{tot}}/(4\pi R_\mathrm{L}^2c)`$, and $`n3`$ (for more details see, e.g., Sikora et al 1996). Here the first term describes the radiation from the outer parts of the disk, $`r_{\mathrm{rad}}>r`$, while the second one corresponds to the homogeneous external radiation. For AGNs this can be due to clouds located at a distance $`r_{\mathrm{cloud}}1`$pc from the central engine and reradiating $`kL_{\mathrm{tot}}`$ of the total luminosity ($`k10\%`$). In this case

$$U_{\mathrm{ext}}=k\frac{L_{\mathrm{tot}}}{4\pi r_{\mathrm{cloud}}^2c}.$$
(2)
However, this model only makes physical sense at distances less than $`r_{\mathrm{cloud}}`$, and the term vanishes at larger distances.

Finally, as some arguments exist both in favor of (Reynolds et al 1996; Hirotani et al 1999) and against (Sikora & Madejski 2000) the leading role of $`e^+e^{}`$ plasma in relativistic jets, in what follows we consider both electron-positron and electron-proton outflows.

In Section 2 we formulate the basic equations describing a quasi-monopole outflow of relativistic plasma in two-fluid and one-fluid approximations. Then in Section 3 we analyze the main properties of an electron-positron outflow. A similar analysis for electron-proton plasma is produced in Section 4. Finally in Section 5 we consider the effects of radiation drag for real astrophysical objects.

## 2 Basic Equations

### 2.1 The Two-Fluid Description

We consider an axisymmetric, stationary outflow of two-component cold plasma from the magnetosphere of a rotating body with a split monopole poloidal magnetic field. This geometry can be realized in the presence of an accretion disk separating the ingoing and outgoing magnetic fluxes (Blandford & Znajek 1977). In the hydrodynamic approximation, the structure of the flow is described by Maxwell’s equations and the separate equations of motion for positively and negatively charged particles:

$`𝑬=4\pi \rho _\mathrm{e},\times 𝑬=0,`$
$`𝑩=0,\times 𝑩={\displaystyle \frac{4\pi }{c}}𝒋,`$ (3)
$`(𝒗^\pm )𝒑^\pm =\pm e\left(𝑬+{\displaystyle \frac{𝒗^\pm }{c}}\times 𝑩\right)+𝑭_{\mathrm{drag}}^\pm .`$
Here $`𝑬`$ and $`𝑩`$ are the electric and magnetic fields, $`\rho _\mathrm{e}=e(n^+n^{})`$ and $`𝒋=e(n^+𝒗^+n^{}𝒗^{})`$ are the charge and current densities, and $`𝒗^\pm `$ and $`𝒑^\pm `$ are the velocities and momenta of the charged particles. Finally, $`𝑭_{\mathrm{drag}}^\pm `$ are the radiation drag forces which, for an isotropic photon field, have the form

$$𝑭_{\mathrm{drag}}^\pm =\frac{4}{3}\frac{𝒗^\pm }{|v^\pm |}\left(\frac{m_\mathrm{e}}{m_\pm }\right)^2\sigma _\mathrm{T}U_{\mathrm{iso}}(\gamma ^\pm )^2.$$
(4)
Here $`\gamma ^\pm `$ are the Lorentz factors, $`m_+`$ is the mass of the positively charged particles (positrons or protons) and $`m_{}=m_\mathrm{e}`$ is the mass of the negatively charged particles. $`\sigma _\mathrm{T}=(8\pi /3)(e^2/m_\mathrm{e}c^2)^2`$ is the Thompson cross section. To close system (3) the continuity equations $`(n^\pm 𝒗^\pm )=0`$ should be added. It is enough to add the equation for one component, e.g.

$$(n^+𝒗^+)=0,$$
(5)
because the continuity equation for the second component then follows from (5) and the current continuity equation $`𝒋=0`$ (which, in turn, results from Maxwell’s equation $`\times 𝑩=(4\pi /c)𝒋`$).

In the limit of infinite particle energy,

$$\gamma =\mathrm{},v_r^{(0)}=c,v_\phi ^{(0)}=0,v_\theta ^{(0)}=0,$$
(6)
so that

$$\rho _\mathrm{e}=\rho _s\frac{R^2}{r^2}\mathrm{cos}\theta ,j_r=\rho _sc\frac{R^2}{r^2}\mathrm{cos}\theta ,j_\theta =j_\phi =0,$$
(7)
the monopole poloidal magnetic field

$$B_r^{(0)}=B_0\frac{R^2}{r^2},B_\theta ^{(0)}=0,$$
(8)
is an exact solution of Maxwell’s equations. In this case,

$`B_\phi ^{(0)}`$ $`=`$ $`E_\theta ^{(0)}=B_0{\displaystyle \frac{\mathrm{\Omega }R}{c}}{\displaystyle \frac{R}{r}}\mathrm{sin}\theta ,`$ (9)
$`E_r^{(0)}`$ $`=`$ $`E_\phi ^{(0)}=0,`$ (10)
which is just the well-known Michel (1973) solution. Here $`B_0`$ and $`\rho _s`$ are the magnetic field and charge density on the surface $`r=RR_\mathrm{L}`$, and the angular velocity is $`\mathrm{\Omega }=2\pi c|\rho _s|/B_0`$. The limit $`\gamma \mathrm{}`$ corresponds to zero particle mass in the force-free approximation.

It is also convenient to introduce the electric field potential $`\mathrm{\Phi }(r,\theta )`$, so that $`𝑬=\mathrm{\Phi }`$ and

$$\mathrm{\Phi }^{(0)}=\frac{\mathrm{\Omega }R^2B_0}{c}\mathrm{cos}\theta ,$$
(11)
and the flux function $`\mathrm{\Psi }(r,\theta )`$, so that the poloidal magnetic field

$$𝑩_\mathrm{p}=\frac{\mathrm{\Psi }\times 𝒆_\phi }{2\pi r\mathrm{sin}\theta },$$
(12)
and $`\mathrm{\Psi }^{(0)}=2\pi B_0R^2(1\mathrm{cos}\theta )`$.

Then the dimensionless corrections $`\eta ^\pm (r,\theta )`$, $`𝝃^\pm (r,\theta )`$, $`\delta (r,\theta )`$, $`\epsilon f(r,\theta )`$ and $`\zeta (r,\theta )`$ for the case $`vc`$ can be introduced in the following form:

$`n^+`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }B_0}{2\pi ce}}{\displaystyle \frac{R^2}{r^2}}\left[\lambda {\displaystyle \frac{1}{2}}\mathrm{cos}\theta +\eta ^+(r,\theta )\right],`$ (13)
$`n^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }B_0}{2\pi ce}}{\displaystyle \frac{R^2}{r^2}}\left[\lambda +{\displaystyle \frac{1}{2}}\mathrm{cos}\theta +\eta ^{}(r,\theta )\right],`$ (14)
$`v_r^\pm `$ $`=`$ $`c\left[1\xi _r^\pm (r,\theta )\right],v_\theta ^\pm =c\xi _\theta ^\pm (r,\theta ),v_\phi ^\pm =c\xi _\phi ^\pm (r,\theta ),`$ (15)
$`\mathrm{\Phi }(r,\theta )`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }R^2B_0}{c}}\left[\mathrm{cos}\theta +\delta (r,\theta )\right],`$ (16)
$`\mathrm{\Psi }(r,\theta )`$ $`=`$ $`2\pi B_0R^2\left[1\mathrm{cos}\theta +\epsilon f(r,\theta )\right],\text{and thus}`$ (17)
$`B_r`$ $`=`$ $`B_0{\displaystyle \frac{R^2}{r^2}}\left(1+{\displaystyle \frac{\epsilon }{\mathrm{sin}\theta }}{\displaystyle \frac{f}{\theta }}\right),`$ (18)
$`B_\theta `$ $`=`$ $`\epsilon {\displaystyle \frac{B_0R^2}{r\mathrm{sin}\theta }}{\displaystyle \frac{f}{r}},`$ (19)
$`B_\phi `$ $`=`$ $`B_0{\displaystyle \frac{R\mathrm{\Omega }}{c}}{\displaystyle \frac{R}{r}}\left[\mathrm{sin}\theta \zeta (r,\theta )\right],`$ (20)
$`E_r`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }B_0R^2}{c}}{\displaystyle \frac{\delta }{r}},`$ (21)
$`E_\theta `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }R^2B_0}{cr}}\left(\mathrm{sin}\theta {\displaystyle \frac{\delta }{\theta }}\right).`$ (22)
Here $`\lambda 1`$ is the multiplication parameter ($`\lambda =en_s/|\rho _s|`$, where $`n_s`$ is the number density of particles on the surface $`r=R`$). For $`\lambda <1`$ the approach under consideration is not valid. In what follows we consider for simplicity the case $`\lambda =`$ const. Such a choice corresponds to a constant particle-to-magnetic flux ratio $`\kappa =`$ const.

Switching to dimensionless variables, we use below the dimensionless radius $`x=r/R_\mathrm{L}=r\mathrm{\Omega }/c`$ and the dimensionless drag force:

$$F_\mathrm{d}^\pm =\frac{4}{3}\frac{\sigma _\mathrm{T}U_{\mathrm{iso}}}{\mathrm{\Omega }m_\mathrm{e}c}(\gamma ^\pm )^2.$$
(23)
Now, substituting (13)–(22) into (3) we obtain, to first order in all the correcting functions, the following system of equations:

$`{\displaystyle \frac{1}{\mathrm{sin}\theta }}{\displaystyle \frac{}{\theta }}(\zeta \mathrm{sin}\theta )=2(\eta ^+\eta ^{})2\left[\left(\lambda {\displaystyle \frac{1}{2}}\mathrm{cos}\theta \right)\xi _r^+\left(\lambda +{\displaystyle \frac{1}{2}}\mathrm{cos}\theta \right)\xi _r^{}\right],`$ (24)
$`2(\eta ^+\eta ^{})+{\displaystyle \frac{}{x}}\left(x^2{\displaystyle \frac{\delta }{x}}\right)+{\displaystyle \frac{1}{\mathrm{sin}\theta }}{\displaystyle \frac{}{\theta }}\left(\mathrm{sin}\theta {\displaystyle \frac{\delta }{\theta }}\right)=0,`$ (25)
$`{\displaystyle \frac{\zeta }{x}}={\displaystyle \frac{2}{x}}\left[\left(\lambda {\displaystyle \frac{1}{2}}\mathrm{cos}\theta \right)\xi _\theta ^+\left(\lambda +{\displaystyle \frac{1}{2}}\mathrm{cos}\theta \right)\xi _\theta ^{}\right],`$ (26)
$`{\displaystyle \frac{\epsilon }{\mathrm{sin}\theta }}{\displaystyle \frac{^2f}{x^2}}{\displaystyle \frac{\epsilon }{x^2}}{\displaystyle \frac{}{\theta }}\left({\displaystyle \frac{1}{\mathrm{sin}\theta }}{\displaystyle \frac{f}{\theta }}\right)={\displaystyle \frac{2}{x}}\left[\left(\lambda {\displaystyle \frac{1}{2}}\mathrm{cos}\theta \right)\xi _\phi ^+\left(\lambda +{\displaystyle \frac{1}{2}}\mathrm{cos}\theta \right)\xi _\phi ^{}\right],`$ (27)
$`{\displaystyle \frac{}{x}}\left(\xi _\theta ^+\gamma ^+\right)+{\displaystyle \frac{\xi _\theta ^+\gamma ^+}{x}}=\xi _\theta ^+F_\mathrm{d}^+\left({\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}\right)^3+4\lambda \sigma \left({\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}\right)\left({\displaystyle \frac{1}{x}}{\displaystyle \frac{\delta }{\theta }}+{\displaystyle \frac{\zeta }{x}}{\displaystyle \frac{\mathrm{sin}\theta }{x}}\xi _r^++{\displaystyle \frac{1}{x^2}}\xi _\phi ^+\right),`$ (28)
$`{\displaystyle \frac{}{x}}\left(\xi _\theta ^{}\gamma ^{}\right)+{\displaystyle \frac{\xi _\theta ^{}\gamma ^{}}{x}}=\xi _\theta ^{}F_\mathrm{d}^{}4\lambda \sigma \left({\displaystyle \frac{1}{x}}{\displaystyle \frac{\delta }{\theta }}+{\displaystyle \frac{\zeta }{x}}{\displaystyle \frac{\mathrm{sin}\theta }{x}}\xi _r^{}+{\displaystyle \frac{1}{x^2}}\xi _\phi ^{}\right),`$ (29)
$`{\displaystyle \frac{}{x}}\left(\gamma ^+\right)=F_\mathrm{d}^+\left({\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}\right)^3+4\lambda \sigma \left({\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}\right)\left({\displaystyle \frac{\delta }{x}}{\displaystyle \frac{\mathrm{sin}\theta }{x}}\xi _\theta ^+\right),`$ (30)
$`{\displaystyle \frac{}{x}}\left(\gamma ^{}\right)=F_\mathrm{d}^{}4\lambda \sigma \left({\displaystyle \frac{\delta }{x}}{\displaystyle \frac{\mathrm{sin}\theta }{x}}\xi _\theta ^{}\right),`$ (31)
$`{\displaystyle \frac{}{x}}\left(\xi _\phi ^+\gamma ^+\right)+{\displaystyle \frac{\xi _\phi ^+\gamma ^+}{x}}=\xi _\phi ^+F_\mathrm{d}^+\left({\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}\right)^3+4\lambda \sigma \left({\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}\right)\left(\epsilon {\displaystyle \frac{1}{x\mathrm{sin}\theta }}{\displaystyle \frac{f}{x}}{\displaystyle \frac{1}{x^2}}\xi _\theta ^+\right),`$ (32)
$`{\displaystyle \frac{}{x}}\left(\xi _\phi ^{}\gamma ^{}\right)+{\displaystyle \frac{\xi _\phi ^{}\gamma ^{}}{x}}=\xi _\phi ^{}F_\mathrm{d}^{}4\lambda \sigma \left(\epsilon {\displaystyle \frac{1}{x\mathrm{sin}\theta }}{\displaystyle \frac{f}{x}}{\displaystyle \frac{1}{x^2}}\xi _\theta ^{}\right).`$ (33)
Here

$$\sigma =\frac{\mathrm{\Omega }eB_0R^2}{4\lambda m_\mathrm{e}c^3}1$$
(34)
is the Michel (1969) magnetization parameter describing the particle-to-electromagnetic energy flux ratio $`W_{\mathrm{part}}/W_{\mathrm{em}}=(m_\mathrm{p}/m_\mathrm{e})\gamma /\sigma `$. Hence, for a Poynting-dominated flow we have $`\gamma \sigma `$. As we see, the disturbances of the particle density $`\eta ^+`$ and $`\eta ^{}`$ enter equations (24)–(33) only in the combination $`\eta ^+\eta ^{}`$. Therefore, the system (24)–(33) is closed, but the equation of mass continuity (5) is necessary to determine $`\eta ^+`$ and $`\eta ^{}`$ separately. The system (24)–(33) differs from the one considered by Beskin & Rafikov (2000, hereafter Paper II) only by the additional drag terms in the r.h.s. of (28)–(33).

### 2.2 The One-Fluid Limit

Equations (24)–(33) describe the flow in the two-fluid approximation. We now reduce the complete system of equations (24)–(33) to consider the one-fluid approximation. In the case of an electron-proton outflow there is a small parameter $`m_\mathrm{e}/m_\mathrm{p}10^3`$ which allows us to neglect the electron mass and thus to proceed in a standard way to the one-fluid approximation. But as demonstrated in Paper II, this can be done in the electron-positron case for a magnetically-dominated ($`\sigma 1`$) dense ($`\lambda 1`$) plasma as well. Because the nonhydrodynamic components of the velocity are small in this case (cf. Melatos & Melrose 1996),

$$\frac{\mathrm{\Delta }\xi _r^\pm }{\xi _r}\lambda ^1\sigma ^{2/3},\frac{\mathrm{\Delta }\xi _\theta ^\pm }{\xi _\theta }\lambda ^1,\frac{\mathrm{\Delta }\xi _\phi ^\pm }{\xi _\phi }\lambda ^1\sigma ^{2/3},$$
(35)
we can set $`\xi _i^+=\xi _i^{}=\xi _i`$ ($`i=r,\theta ,\phi `$), where $`\xi _i`$ is the hydrodynamic velocity. In this limit we also have

$$\frac{\delta \epsilon f}{\epsilon f}\lambda ^2\sigma ^{2/3}.$$
(36)
As a result, for a magnetically-dominated outflow with $`\sigma 1`$ and $`\lambda 1`$ in the one-fluid approximation,

$$\delta =\epsilon f.$$
(37)
Finally, equations (28) and (29) together with (37) in the limit $`\lambda \sigma 1`$ under consideration give another useful one-fluid relation:

$$\epsilon \frac{1}{x}\frac{f}{\theta }+\frac{\zeta }{x}\frac{\mathrm{sin}\theta }{x}\xi _r+\frac{1}{x^2}\xi _\phi =0,$$
(38)
this equation being the same as in the drag-free case. On the other hand, as demonstrated in Paper II,

$$\frac{\xi _\theta }{\xi _\phi }\sigma ^{1/3}.$$
(39)
Hence, in the one-fluid approximation one can set $`\xi _\theta =0`$ so that

$$\gamma ^2=\frac{1}{2\xi _r\xi _\phi ^2}.$$
(40)
It is necessary to stress that in some cases the monopole geometry allows one to consider separately the set of equations describing the particle energy and the set of equations resulting in the Grad-Shafranov (GS) equation, which determines the disturbance of magnetic surfaces. Thus we can consider particle energy without formulating here the general form of the GS equation. Some asymptotics of the GS equation are discussed below.

## 3 The Electron-Positron Outflow

### 3.1 Integrals of motion

In this section we consider the properties of the electron-positron outflow when $`m_\mathrm{p}=m_\mathrm{e}=m`$. In our simple geometry with a (split) monopole poloidal magnetic field in the zeroth approximation the particle motion can be considered as radial. It allows us to integrate some equations along the $`r`$ coordinate axis. This feature was first demonstrated by Li et al (1992).

Indeed, combining (30) and (31) with (26) one can obtain in the one-fluid approximation

$$\zeta =\frac{l(\theta )}{\mathrm{sin}\theta }+\frac{2\epsilon }{\mathrm{tan}\theta }f\frac{1}{\sigma \mathrm{sin}\theta }(\gamma \gamma _{\mathrm{in}})\frac{1}{\sigma \mathrm{sin}\theta }l_\mathrm{A}_{x_0}^xu(x^{})\gamma ^2(x^{})dx^{}.$$
(41)
Here $`\gamma _{\mathrm{in}}`$ is the Lorentz factor near the origin, $`x_0=\mathrm{\Omega }R/c`$, and the integration constant $`l(\theta )`$ describes the disturbance of the electric current $`I(R,\theta )=I_A\left[\mathrm{sin}^2\theta +l(\theta )\right]`$ on the surface $`r=R`$. The integration constant must be determined from the critical conditions on singular surfaces. The compactness parameter $`l_\mathrm{A}`$ is defined by the relation

$$l_\mathrm{A}=\frac{4}{3}\frac{\sigma _\mathrm{T}U_\mathrm{A}}{m_\mathrm{e}c\mathrm{\Omega }},$$
(42)
where $`U_\mathrm{A}=U(R_\mathrm{L})`$ is again the first term in (1) on the light cylinder $`r=R_\mathrm{L}`$. Finally, we also denote

$$l_{\mathrm{ext}}=\frac{4}{3}\frac{\sigma _\mathrm{T}U_{\mathrm{ext}}}{mc\mathrm{\Omega }},$$
(43)
so that (1) for $`x>1`$ can be rewritten as

$$u(x)=x^n+l_{\mathrm{ext}}/l_\mathrm{A}.$$
(44)
For $`x<1`$ it is natural to set $`u(x)=1`$.

Expression (41) is just another way of writing the diminishing of the Bernoulli integral $`E_\mathrm{B}`$ along a magnetic field line as a result of the drag force. Within the Grad-Shafranov approach the full energy loss is determined by $`W=E_\mathrm{B}d\mathrm{\Psi }`$. Thus relation (41) can be rewritten as

$$E_\mathrm{B}(r)=E_\mathrm{B}(R)\frac{\lambda \mathrm{\Omega }}{2\pi e}_R^rF_{\mathrm{drag}}dr,$$
(45)
where the energy flux per unit magnetic flux $`E_\mathrm{B}`$ has the standard form $`E_\mathrm{B}=\mathrm{\Omega }I/2\pi c+\gamma mc^2\kappa `$. Here $`I`$ is the total electric current inside the magnetic tube and $`\kappa `$ is the particle-to-magnetic flux ratio (see e.g. Beskin 1997 for details).

In what follows we consider the case $`\gamma _{\mathrm{in}}1`$, i.e.

$$\gamma _{\mathrm{in}}^3\sigma ,$$
(46)
when additional acceleration of particles inside the fast magnetosonic surface takes place (see e.g. Paper I). In the case $`\gamma _{\mathrm{in}}^3\sigma `$ corresponding to ordinary pulsars, the particle energy remains constant ($`\gamma =\gamma _{\mathrm{in}}`$) on any way up to the fast magnetosonic surface (Bogovalov 1997).

The other two integrals of motion, namely the conservation of angular momentum separately for electrons and positrons, can be obtained from equations (30)–(33):

$`\gamma ^+(1x\mathrm{sin}\theta \xi _\phi ^+)=\gamma _{\mathrm{in}}^+l_\mathrm{A}{\displaystyle _{x_0}^x}u(x^{})(1x^{}\mathrm{sin}\theta \xi _\phi ^+)(\gamma ^+)^2dx^{}4\lambda \sigma (\delta \epsilon f),`$ (47)
$`\gamma ^{}(1x\mathrm{sin}\theta \xi _\phi ^{})=\gamma _{\mathrm{in}}^{}l_\mathrm{A}{\displaystyle _{x_0}^x}u(x^{})(1x^{}\mathrm{sin}\theta \xi _\phi ^{})(\gamma ^{})^2dx^{}+4\lambda \sigma (\delta \epsilon f).`$ (48)
In the one-fluid approximation this gives

$$\gamma (1x\mathrm{sin}\theta \xi _\phi )=\gamma _{\mathrm{in}}l_\mathrm{A}_{x_0}^xu(x^{})(1x^{}\mathrm{sin}\theta \xi _\phi )\gamma ^2(x^{})dx^{}.$$
(49)
At large distances, where $`\gamma \gamma _{\mathrm{in}}`$, we have

$$\xi _\phi \frac{1}{x\mathrm{sin}\theta }.$$
(50)
It can be seen from (49) that the presence of the drag force (the term proportional to $`l_\mathrm{A}`$) makes this approximation more accurate.

### 3.2 Fast Magnetosonic Surface

Substituting $`\zeta `$ from (41), $`\xi _r`$ from (40), and $`\xi _\phi `$ from (49) into (38) we get the following equation to determine the position of the fast magnetosonic surface $`r=r_\mathrm{F}`$:

$`\epsilon {\displaystyle \frac{f}{\theta }}+{\displaystyle \frac{2\epsilon }{\mathrm{tan}\theta }}f+{\displaystyle \frac{l(\theta )}{\mathrm{sin}\theta }}{\displaystyle \frac{1}{\sigma \mathrm{sin}\theta }}\left(\gamma \gamma _{\mathrm{in}}+J\right)\mathrm{sin}\theta \left[{\displaystyle \frac{1}{2\gamma ^2}}+{\displaystyle \frac{\left(\gamma \gamma _{\mathrm{in}}+J_1\right)^2}{2x^2\gamma ^2\mathrm{sin}^2\theta }}\right]+{\displaystyle \frac{\gamma \gamma _{\mathrm{in}}+J_1}{x^2\gamma \mathrm{sin}\theta }}=0.`$ (51)
Here

$`J`$ $`=`$ $`l_\mathrm{A}{\displaystyle _{x_0}^x}u(x^{})\gamma ^2(x^{})dx^{},`$ (52)
$`J_1`$ $`=`$ $`l_\mathrm{A}{\displaystyle _{x_0}^x}u(x^{})\gamma ^2(x^{})(1x^{}\mathrm{sin}\theta \xi _\phi )dx^{}.`$ (53)
But, in accordance with (49) and thereafter, the terms $`J_1`$ and $`\gamma _{\mathrm{in}}`$ can be neglected far from the origin of the flow $`rR_\mathrm{L}`$. Also, since $`l(\theta )\sigma ^{4/3}`$ (see Paper I), the term with $`l(\theta )`$ can be omitted as well. Hence, one can write down the following algebraic equation for the Lorentz factor $`\gamma `$:

$$\gamma ^3\sigma \left(P+\frac{1}{2x^2}\right)\gamma ^2+\frac{1}{2}\sigma \mathrm{sin}^2\theta =0,$$
(54)
where

$$P=\frac{J}{\sigma }+2\epsilon f\mathrm{cos}\theta \epsilon \mathrm{sin}\theta \frac{f}{\theta }.$$
(55)
This differs from the drag-free case by the additional term $`J/\sigma `$.

Equation (54) allows us to determine the position of the fast magnetosonic surface $`r=r_\mathrm{F}`$ and the energy of particles on this surface $`\gamma _\mathrm{F}=\gamma (r_\mathrm{F})`$. Indeed, as the fast magnetosonic surface is the $`X`$–point, we find the exact solution for coinciding roots

$$\gamma _\mathrm{F}=\sigma ^{1/3}\mathrm{sin}^{2/3}\theta ,$$
(56)
which does not depend at all on $`P`$, or thus on the drag force. On the other hand, both the numerator and the denominator of the derivative $`\mathrm{d}\gamma /\mathrm{d}x`$ must be equal to zero on the surface $`r=r_\mathrm{F}`$:

$$x\frac{\mathrm{d}\gamma }{\mathrm{d}x}=\frac{\gamma \sigma \left(x\mathrm{d}P/\mathrm{d}xx^2\right)}{3\gamma \sigma \left(2P+x^2\right)}.$$
(57)
It gives

$$(P+x^2)_\mathrm{F}\sigma ^{2/3},|P|_\mathrm{F}|x^2|_\mathrm{F}.$$
(58)
Now using (55), one can see that the conditions for weak and strong drag are, respectively, $`l_\mathrm{A}l_{\mathrm{cr}}`$ and $`l_\mathrm{A}l_{\mathrm{cr}}`$, where

$$l_{\mathrm{cr}}=\sigma ^{1/3}.$$
(59)
In other words, for $`l_\mathrm{A}l_{\mathrm{cr}}`$ the flow remains the same as in the drag-free case, and

$`x_\mathrm{F}`$ $``$ $`\sigma ^{1/3}\mathrm{sin}^{1/3}\theta ,`$ (60)
$`(\epsilon f)_\mathrm{F}`$ $``$ $`\sigma ^{2/3}.`$ (61)
On the other hand, for a high enough photon density $`l_\mathrm{A}l_{\mathrm{cr}}`$ one can obtain, for $`n=3`$,

$`x_\mathrm{F}`$ $``$ $`\left({\displaystyle \frac{\sigma }{l_\mathrm{A}}}\right)^{1/2}=\sigma ^{1/3}\left({\displaystyle \frac{l_{\mathrm{cr}}}{l_\mathrm{A}}}\right)^{1/2},`$ (62)
$`(\epsilon f)_\mathrm{F}`$ $``$ $`{\displaystyle \frac{l_\mathrm{A}}{\sigma }}.`$ (63)
Thus, as we can see, the energy $`\gamma _\mathrm{F}mc^2`$ do not depend on the drag. The disturbance of the magnetic surfaces increases with increasing $`l_\mathrm{A}`$, but remains small for $`l_\mathrm{A}\sigma `$.

For a very dense photon field, according to (41) and (63), the disturbance of the magnetic surfaces becomes of order unity for $`l_\mathrm{A}l_{\mathrm{max}}`$ where

$$l_{\mathrm{max}}=\sigma ,$$
(64)
i.e., when the disturbance of the Bernoulli integral $`\mathrm{\Delta }E_\mathrm{B}=E_\mathrm{B}(r)E_\mathrm{B}(R_\mathrm{L})`$ (45) is on the order of the total energy flux $`E_\mathrm{B}`$, or in other words when the drag force substantially diminishes the total energy flux of the flow. But this means that the disturbance of the magnetic surfaces becomes large when the energies of particles in the vicinity of the fast magnetosonic surface become nonrelativistic. This feature is well known for a drag-free flow both within numerical (Sakurai 1985; Bogovalov 1997) and analytical (Bogovalov 1992; Tomimatsu 1994; Paper I) considerations. Of course, our analysis within the small disturbance approach can only demonstrate the tendency. For this reason, within our approach $`\sigma =`$ const. In reality, since $`\sigma =E_\mathrm{B}/(mc^2\kappa )`$, a decrease of the energy flux results in a decrease of $`\sigma `$ as well.

Finally, it is necessary to stress that the results depend significantly on the assumed model of the isotropic photon density. In particular, for the general power-law dependence of the photon density for $`n<3`$ in (44), one can obtain for the position of the fast magnetosonic point and the disturbance of the magnetic surfaces

$`x_\mathrm{F}`$ $``$ $`\sigma ^{1/3}\left({\displaystyle \frac{l_{\mathrm{cr}}}{l_\mathrm{A}}}\right)^{1/(5n)},`$ (65)
$`(\epsilon f)_\mathrm{F}`$ $``$ $`\left({\displaystyle \frac{l_\mathrm{A}}{\sigma }}\right)^{2/(5n)}`$ (66)
instead of (62) and (63). This takes place for $`l_\mathrm{A}>l_{\mathrm{cr}}`$, where $`l_{\mathrm{cr}}`$ is the photon density capable of disturbing the magnetically-dominated outflow:

$$l_{\mathrm{cr}}=\sigma ^{(n2)/3}.$$
(67)
A large disturbance $`(\epsilon f)_\mathrm{F}1`$ can only be realized for a very high photon density $`l_\mathrm{A}>l_{\mathrm{max}}`$, where again

$$l_{\mathrm{max}}=\sigma .$$
(68)
For $`n3`$ (when the value J defined by (52) is determined by the lower limit of integration and does not depend on $`n`$) we have $`(\epsilon f)_\mathrm{F}l_\mathrm{A}/\sigma `$, $`x_\mathrm{F}(\sigma /l_\mathrm{A})^{1/2}`$, $`l_{\mathrm{cr}}=\sigma ^{1/3}`$, and $`l_{\mathrm{max}}=\sigma `$, i.e. same as for $`n=3`$. According to (41), $`(\epsilon f)_\mathrm{F}1`$ occurs when $`\mathrm{\Delta }E_\mathrm{B}/E_\mathrm{B}1`$ for all values of $`n`$.

### 3.3 The Structure of the Flow at Small and Large Distances

#### 3.3.1 Inner region

We first consider the structure of the flow well within the fast magnetosonic surface $`rr_\mathrm{F}`$. As one can see from (54), for $`rr_\mathrm{F}`$

$$\gamma \frac{\mathrm{sin}\theta }{\left(2P+1/x^2\right)^{1/2}}.$$
(69)
In particular, for $`l_\mathrm{A}l_{\mathrm{cr}}`$ one can neglect the term $`P`$ in the denominator. On the other hand, for $`l_\mathrm{A}l_{\mathrm{cr}}`$ in the immediate vicinity of the fast magnetosonic surface the negative term $`J/\sigma `$ in $`P`$ should to be included into consideration. As a result, the particle energy increases abruptly up to the value $`\sigma ^{1/3}`$. Nevertheless, for $`rr_\mathrm{F}`$ the value $`P`$ can be omitted for arbitrary compactness parameter $`l_\mathrm{A}`$. It also means that for $`rr_\mathrm{F}`$ it is possible to neglect the first two terms in (38). Thus, in the internal region $`rr_\mathrm{F}`$ we have

$$\xi _r=\frac{\xi _\phi }{x\mathrm{sin}\theta }.$$
(70)
Now using relations (50) and (70), one can obtain

$`\gamma ^2`$ $`=`$ $`\gamma _{\mathrm{in}}^2+x^2\mathrm{sin}^2\theta x^2\mathrm{sin}^2\theta ,`$ (71)
$`\xi _\phi `$ $`=`$ $`{\displaystyle \frac{\sqrt{\gamma _{\mathrm{in}}^2+x^2\mathrm{sin}^2\theta }\gamma _{\mathrm{in}}}{x\mathrm{sin}\theta \sqrt{\gamma _{\mathrm{in}}^2+x^2\mathrm{sin}^2\theta }}}{\displaystyle \frac{1}{x\mathrm{sin}\theta }},`$ (72)
$`\xi _r`$ $`=`$ $`{\displaystyle \frac{\sqrt{\gamma _{\mathrm{in}}^2+x^2\mathrm{sin}^2\theta }\gamma _{\mathrm{in}}}{x^2\mathrm{sin}^2\theta \sqrt{\gamma _{\mathrm{in}}^2+x^2\mathrm{sin}^2\theta }}}{\displaystyle \frac{1}{x^2\mathrm{sin}^2\theta }},`$ (73)
in full agreement with the ideal MHD approximation. Hence, one can conclude that in the internal region $`r<r_\mathrm{F}`$ the radiation drag does not affect the particle motion. Here the universal value of the Lorentz factor (71) corresponds to the drift velocity, and it depends on no external disturbances. Indeed, using the frozen-in condition $`𝑬+𝒗\times 𝑩/c=0`$, we obtain for the drift velocity

$$U_{\mathrm{dr}}^2=c^2\frac{𝑬^2}{𝑩^2}=c^2\left(\frac{B_\phi ^2}{E^2}+\frac{B_r^2}{E^2}\right)^1.$$
(74)
In our case, however, according to (16)–(22), we have

$$B_\phi ^2E^2,B_r^2E^2/(x^2\mathrm{sin}^2\theta ).$$
(75)
These relations immediately lead to the previously-stated asymptotic behavior (71). In particular, this means that disturbance of the magnetic surfaces plays no role in the determination of the particle energy. For this reason, it is possible not to consider the radiation drag corrections to the field structure for $`rr_\mathrm{F}`$.

#### 3.3.2 Outer region

In the other limit, well outside the fast magnetosonic surface ($`rr_\mathrm{F}`$) equations (41), (38), (24)–(25), and (40) can be rewritten in the form

$`\zeta =\epsilon {\displaystyle \frac{2}{\mathrm{tan}\theta }}f{\displaystyle \frac{1}{\sigma \mathrm{sin}\theta }}\gamma {\displaystyle \frac{1}{\sigma \mathrm{sin}\theta }}l_\mathrm{A}{\displaystyle _{x_0}^x}u(x^{})\gamma ^2(x^{})dx^{},`$ (76)
$`\zeta =\epsilon {\displaystyle \frac{f}{\theta }}+\mathrm{sin}\theta \xi _r,`$ (77)
$`\epsilon {\displaystyle \frac{}{x}}\left(x^2{\displaystyle \frac{f}{x}}\right){\displaystyle \frac{1}{\mathrm{sin}^3\theta }}{\displaystyle \frac{}{\theta }}\left(\mathrm{sin}^4\theta \xi _r\right)=0,`$ (78)
$`\gamma ^2={\displaystyle \frac{1}{2\xi _r}}.`$ (79)
Here we neglect the small values $`l(\theta )`$, $`\gamma _{\mathrm{in}}/\gamma `$ and $`\xi _\phi `$, and set $`\delta =\epsilon f`$. As we can see, the only correction compared to the drag-free case is the last additional term in (76). Without the drag force the system (76)–(79) results in a very slow increase of the particle energy ($`\gamma \mathrm{ln}^{1/3}r`$) and actually in the absence of collimation, $`\epsilon f\sigma ^{2/3}\mathrm{ln}^{1/3}r`$ (Tomimatsu 1994; Paper I). Unfortunately, because of the nonlinearity of the system (76)–(79), in the general case it is impossible to reduce it to the GS equation for the magnetic disturbance $`\epsilon f`$ only. For this reason, in what follows we present only an asymptotic representation of the GS equation.

a) Low photon density $`l_\mathrm{A}l_{\mathrm{cr}}`$.

In this case the action of drag results in a small correction to the particle energy. Only at very large distances, where the first term in (44) can be neglected, does the drag significantly change the energy of particles. Indeed, for the decreasing component of the isotropic photon density ($`u(x)x^n`$ in (44)) and for an almost constant particle energy the drag term in (76) does not increase with the distance $`r`$. Here the increase of the drag term is due to the homogeneous part of the photon density ($`l_{\mathrm{ext}}/l_\mathrm{A}`$ in (44)). Clearly, this can be realized at large enough distances ($`r>R_\mathrm{L}(l_\mathrm{A}/l_{\mathrm{ext}})\sigma ^{(3n)/3}`$), where the contribution from external photons is the leading one. Thus, we see that the action of the drag force can be significant for a high enough density of external isotropic photons in the vicinity of the compact object.

Now we consider this asymptotic region in more detail. Neglecting the term containing $`\gamma `$ in (76) and the first term in (44), one can obtain for the generalized GS equation

$$\epsilon \frac{}{x}\left(x^2\frac{f}{x}\right)\frac{1}{\mathrm{sin}^3\theta }\frac{}{\theta }\left[2\epsilon \mathrm{sin}^2\theta \mathrm{cos}\theta f\epsilon \mathrm{sin}^3\theta \frac{f}{\theta }\frac{\mathrm{sin}^2\theta }{\sigma }l_{\mathrm{ext}}_{x_0}^x\gamma ^2(x^{})dx^{}\right]=0.$$
(80)
One can seek the solution of this equation in the form

$$\epsilon f(x,\theta )x^{\alpha _\delta },\gamma (x,\theta )x^{\alpha _\gamma }.$$
(81)
Substituting these expressions into (80) we obtain

$`\epsilon f(x,\theta )`$ $`=`$ $`k_1(\theta )l_{\mathrm{ext}}^{1/2}\sigma ^{1/2}x^{1/2},`$ (82)
$`\gamma (x,\theta )`$ $`=`$ $`k_2(\theta )l_{\mathrm{ext}}^{1/4}\sigma ^{1/4}x^{1/4},`$ (83)
where $`k_1(\theta )k_2(\theta )1`$ describe the $`\theta `$ dependence. This takes place for $`r>r_\mathrm{h}`$, where

$$r_\mathrm{h}=R_\mathrm{L}l_{\mathrm{ext}}^1\sigma ^{1/3}.$$
(84)
For $`r_\mathrm{F}<r<r_\mathrm{h}`$ the action of drag is negligible.

As we see, for $`r>r_\mathrm{h}`$ the drag force results in a decrease of the particle energy and an additional collimation of the magnetic surfaces outside the fast magnetosonic surface. These asymptotic solutions, however, do not seem to be realized in reality. Indeed, as one can check, $`r_\mathrm{h}(10^610^8)R_\mathrm{L}>r_{\mathrm{cloud}}`$. But for $`r>r_{\mathrm{cloud}}`$ the model that we adopted for the homogeneous component of the photon density is no longer valid and thus the asymptotic solutions (82)–(83) do not apply. This conclusion is made under the assumption that the homogeneous component is produced by the reradiating clouds surrounding the central engine. One could imagine producing a homogeneous component of the photon density in the host galaxy. However, in this case the photon density $`U_{\mathrm{ext}}`$ is too low and thus the collimation distance $`r_\mathrm{h}`$ is too large for our approach to be valid in the region $`rr_\mathrm{h}`$.

The characteristic radial dependence of the particle energy in the presence of the drag is demonstrated in Fig. $`1a`$. One can conclude that for a low photon density the action of the drag force is very weak unless photon field is present out to distances $`rr_\mathrm{h}`$. In the latter case the drag force efficiently reduces the particle energy beyond $`rr_\mathrm{h}`$.

It is necessary to stress here an important property of the fast magnetosonic surface for $`l_\mathrm{A}<l_{\mathrm{cr}}`$. Introducing a small disturbance of the particle energy resulting from the drag

$$\gamma =\gamma _0+\mathrm{\Delta }\gamma ,$$
(85)
where $`\gamma _0`$ is the Lorentz factor of particles without the drag force, one can obtain from (76)–(79)

$$\mathrm{sin}^2\theta \frac{\mathrm{\Delta }\gamma }{\gamma _0^3}+\epsilon \frac{\mathrm{\Delta }f}{\theta }2\epsilon \mathrm{cos}\theta \mathrm{\Delta }f+\frac{1}{\sigma }\mathrm{\Delta }\gamma =\frac{J}{\sigma }.$$
(86)
As equation (78) results in $`\epsilon \mathrm{\Delta }f\mathrm{\Delta }\gamma /\gamma _0^3`$, we see that outside the fast magnetosonic surface (where $`\gamma _0>\sigma ^{1/3}`$)

$$mc^2\mathrm{\Delta }\gamma F_{\mathrm{drag}}dr.$$
(87)
On the other hand, within the fast magnetosonic surface we have

$$mc^2|\mathrm{\Delta }\gamma |F_{\mathrm{drag}}dr,$$
(88)
in full agreement with (71). This means that inside the fast magnetosonic surface the decrease of the total energy flux $`E_\mathrm{B}`$ results from the decrease of the electromagnetic energy flux rather than the particle flux. Outside the fast surface the photons act upon the particle directly. In other words, for $`l_\mathrm{A}<l_{\mathrm{cr}}`$ the fast magnetosonic surface separates the regions in space where particles are strongly or weakly frozen into the electromagnetic field.

b) High photon density $`l_\mathrm{A}>l_{\mathrm{cr}}`$.

For a high photon density $`l_\mathrm{A}>l_{\mathrm{cr}}`$ the drag force significantly changes the energy of the outgoing particles. As is shown in Fig. $`1b`$, the increase of the particle energy continues outside the fast magnetosonic surface up to the maximum value

$$\gamma _{\mathrm{max}}\left(\frac{l_\mathrm{A}}{l_{\mathrm{cr}}}\right)^{4(n2)/(5nn^2)}\sigma ^{1/3}.$$
(89)
For $`n<2`$ there is no additional acceleration. Only at larger distances ($`r>r_{\mathrm{cr}}`$) does the radial dependence of the energy of the outgoing plasma become similar to the previous case.

Indeed, outside the fast magnetosonic surface one can neglect the term $`\mathrm{sin}\theta \xi _r`$ in (77). As a result, the system (76)–(79) can be rewritten in the form

$`\epsilon {\displaystyle \frac{f}{\theta }}\epsilon {\displaystyle \frac{2}{\mathrm{tan}\theta }}f={\displaystyle \frac{1}{\sigma \mathrm{sin}\theta }}\gamma {\displaystyle \frac{1}{\sigma \mathrm{sin}\theta }}l_\mathrm{A}{\displaystyle _{x_0}^x}u(x^{})\gamma ^2(x^{})dx^{},`$ (90)
$`\epsilon {\displaystyle \frac{}{x}}\left(x^2{\displaystyle \frac{f}{x}}\right){\displaystyle \frac{1}{2\mathrm{sin}^3\theta }}{\displaystyle \frac{}{\theta }}\left({\displaystyle \frac{\mathrm{sin}^4\theta }{\gamma ^2}}\right)=0,`$ (91)
which results in

$$\frac{}{x}\left(x^2\frac{\gamma }{x}\right)+l_\mathrm{A}\frac{}{x}\left(x^{2n}\gamma ^2\right)=\frac{\sigma }{2\mathrm{sin}^3\theta }\left[2\mathrm{cos}\theta \frac{}{\theta }\left(\frac{\mathrm{sin}^4\theta }{\gamma ^2}\right)\mathrm{sin}\theta \frac{^2}{\theta ^2}\left(\frac{\mathrm{sin}^4\theta }{\gamma ^2}\right)\right].$$
(92)
For $`l_\mathrm{A}>l_{\mathrm{cr}}`$ the Lorentz factor $`\gamma `$ increases with $`x`$, and the r.h.s. in (92) can be omitted. This gives

$$x^2\frac{\gamma }{x}+l_\mathrm{A}x^{2n}\gamma ^2=C(\theta ).$$
(93)
Here the integration constant $`C(\theta )l_\mathrm{A}x_\mathrm{F}^{2n}\gamma _\mathrm{F}^2`$ for $`l_\mathrm{A}l_{\mathrm{cr}}`$ actually does not depend on the boundary conditions $`\mathrm{d}\gamma /\mathrm{d}x`$ at $`x=x_\mathrm{F}`$. In what follows we do not consider the $`\theta `$-dependence.

Now analyzing equation (93), one can find that for $`x_\mathrm{F}xx_{\mathrm{cr}}`$, where

$$x_{\mathrm{cr}}\left(\frac{l_\mathrm{A}}{l_{\mathrm{cr}}}\right)^{8/(5nn^2)}x_\mathrm{F},$$
(94)
the particle energy increases as

$$\gamma (x)\sigma ^{1/3}\left(\frac{x}{x_\mathrm{F}}\right)^{(n2)/2}$$
(95)
for arbitrary $`n`$. As to the disturbance of magnetic surfaces $`\epsilon f`$, it remains approximately constant: $`\epsilon f(l_\mathrm{A}/\sigma )^{2/(5n)}`$.

On the other hand, in the saturation region $`xx_{\mathrm{cr}}`$, where one can now neglect the term $`l_\mathrm{A}x^{2n}\gamma ^2`$ in (93), we have

$$\gamma (x)=C\left(\frac{1}{x_\mathrm{b}}\frac{1}{x}\right),$$
(96)
where $`x_\mathrm{b}x_{\mathrm{cr}}`$. One can easily check that the maximum energy $`\gamma _{\mathrm{max}}C/x_{\mathrm{cr}}\sigma ^{1/3}(x_{\mathrm{cr}}/x_\mathrm{F})^{(n2)/2}`$ corresponds to (89). As demonstrated in Fig. 2., the analytical estimate (89) is in good agreement with the numerical integration of equation (93).

#### 3.3.3 Physical interpretation of particle acceleration

At first glance, it is quite unexpected that a drag force can result in an acceleration of particles (so even the word ‘drag’ itself is not appropriate any more). As was demonstrated, additional acceleration can be realized for a highly magnetized outflow when the photon density decreases rapidly with radius ($`n>2`$). This acceleration occurs as a result of the action of the drag force onto magnetic surfaces. Thus, acceleration is only obtained when the disturbance of the magnetic surfaces is taken into account self-consistently.

To understand the nature of the additional particle acceleration in the supersonic region $`r>r_\mathrm{F}`$, it is necessary to return to the system (90)–(91). As one can see, equation (91) actually plays the role of the GS equation, describing the force balance in the direction perpendicular to the magnetic surfaces. It contains no drag term, because the drag force acts along the magnetic surfaces. Indeed, as was shown before (Bogovalov 1998; Okamoto 1999; Beskin & Okamoto 2000), in the asymptotic region $`rr_\mathrm{F}`$ the transfield GS equation can be written down as $`_\mathrm{c}=_{\mathrm{em}}`$, i.e., as a competition between the centrifugal volume force

$$_\mathrm{c}=\frac{nmc^2\gamma +S/c}{R_\mathrm{c}}$$
(97)
and the electromagnetic volume force

$$_{\mathrm{em}}=\rho _\mathrm{e}E_\theta +\frac{1}{c}j_{}B_\phi \frac{1}{8\pi r\mathrm{sin}^2\theta }\frac{}{\theta }\left[(B_\phi ^2E_\theta ^2)\mathrm{sin}^2\theta \right].$$
(98)
Here $`S(c/4\pi )E_\theta B_\phi `$ is the Poynting flux. Now using the expression for the curvature radius $`R_\mathrm{c}`$ (Begelman & Li 1994; Beskin & Okamoto 2000),

$$\frac{1}{R_\mathrm{c}}=\frac{\epsilon }{r\mathrm{sin}\theta }\frac{}{r}\left(r^2\frac{f}{r}\right),$$
(99)
and the condition

$$B_\phi ^2E_\theta ^2\frac{1}{\gamma ^2}B_\phi ^2$$
(100)
resulting from the relativistic Bernoulli equation (see e.g. Bogovalov 1998), for a magnetically-dominated flow $`Snmc^3\gamma `$ we recover relation (91).

On the other hand, the Bernoulli equation (90) describes the change of the total energy flux (45) due to the drag force. As was shown above, for $`l_\mathrm{A}>l_{\mathrm{cr}}`$ in the vicinity of the fast magnetosonic surface the leading terms in the energy equation (90) are those containing $`\epsilon f`$ and $`l_\mathrm{A}`$. This means that here the drag force again acts mainly on the electromagnetic (Poynting) flux $`S`$. The drag diminishes the $`\theta `$-component of the electric field $`E_\theta `$, i.e. it disturbs the equipotential surfaces $`\delta (r,\theta )=`$ const. But since in the one-fluid approximation the magnetic surfaces $`\epsilon f(r,\theta )=`$ const must follow the equipotential surfaces, the decrease of the electric field $`E_\theta `$ results in a change of the curvature of magnetic field lines as well. As the condition $`_\mathrm{c}=_{\mathrm{em}}`$ can now be rewritten in the form

$$\gamma ^2\frac{4\pi j_{}}{cE_\theta }R_\mathrm{c}\frac{R_\mathrm{c}}{r},$$
(101)
we see that the increase of the curvature radius $`R_\mathrm{c}`$ faster than $`r`$ results in an increase of the particle energy. In our case such an acceleration is due to the rapid decrease ($`n>2`$) of the isotropic photon density with the distance $`r`$.

Thus, we see that for high enough photon density there is an additional acceleration of outgoing plasma outside the fast magnetosonic surface. On the other hand, as for $`l_\mathrm{A}>l_{\mathrm{max}}`$ (64) the drag force significantly diminishes the total energy flux $`E`$, one can conclude that this acceleration may only take place up to the energy

$$_{\mathrm{sup}}\sigma ^{1/3}\left(\frac{l_{\mathrm{max}}}{l_{\mathrm{cr}}}\right)^{4(n2)/n(5n)}m_\mathrm{e}c^2$$
(102)
for $`2<n<3`$. On the other hand, $`_{\mathrm{sup}}\sigma ^{7/9}m_\mathrm{e}c^2`$ for $`n3`$. This energy is always much lower than $`\sigma m_\mathrm{e}c^2`$, corresponding to the total conversion of electromagnetic energy intoparticle energy.

## 4 The Electron-Proton Outflow

In this section we briefly consider the results of the analysis of the electron-proton outflow. As the procedure is quite similar to the electron-positron case, we present here the principal relations only.

The main difference from the electron-positron case that occurs due to the large mass ratio is that the drag force acts on the electron component only, but the mass flow is determined entirely by protons. As a result, the energy conservation law has the form

$`\zeta ={\displaystyle \frac{2\epsilon }{\mathrm{tan}\theta }}f{\displaystyle \frac{(\lambda +\mathrm{cos}\theta /2)\gamma ^{}+(m_\mathrm{p}/m_\mathrm{e})(\lambda \mathrm{cos}\theta /2)\gamma ^+}{2\lambda \sigma \mathrm{sin}\theta }}`$
$`{\displaystyle \frac{l_\mathrm{A}}{2\lambda \sigma \mathrm{sin}\theta }}{\displaystyle _{x_0}^x}\left[\left(\lambda +{\displaystyle \frac{1}{2}}\mathrm{cos}\theta \right)(\gamma ^{})^2+\left({\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}\right)^2\left(\lambda +{\displaystyle \frac{1}{2}}\mathrm{cos}\theta \right)(\gamma ^+)^2\right]u(x^{})\gamma ^2(x^{})dx^{}`$
$`{\displaystyle \frac{2\epsilon }{\mathrm{tan}\theta }}f\left({\displaystyle \frac{m_\mathrm{p}}{m_\mathrm{e}}}\right){\displaystyle \frac{\gamma }{2\sigma \mathrm{sin}\theta }}{\displaystyle \frac{1}{2\sigma \mathrm{sin}\theta }}l_\mathrm{A}{\displaystyle _{x_0}^x}u(x^{})\gamma ^2(x^{})dx^{},`$ (103)
where now $`\gamma =\gamma ^+`$. The conservation of angular momentum (47)–(48) looks like

$`\gamma ^+(1x\mathrm{sin}\theta \xi _\phi ^+)`$ $`=`$ $`\gamma _{\mathrm{in}}^+l_\mathrm{A}\left({\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}\right)^3{\displaystyle _{x_0}^x}u(x^{})(1x^{}\mathrm{sin}\theta \xi _\phi ^+)(\gamma ^+)^2dx^{}4\lambda \sigma {\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}(\delta \epsilon f),`$ (104)
$`\gamma ^{}(1x\mathrm{sin}\theta \xi _\phi ^{})`$ $`=`$ $`\gamma _{\mathrm{in}}^{}l_\mathrm{A}{\displaystyle _{x_0}^x}u(x^{})(1x^{}\mathrm{sin}\theta \xi _\phi ^{})(\gamma ^{})^2dx^{}+4\lambda \sigma (\delta \epsilon f),`$ (105)
resulting in $`\xi _\phi =1/x\mathrm{sin}\theta `$. Since for $`rr_\mathrm{F}`$ we again have $`\xi _r=\xi _\phi /(x\mathrm{sin}\theta )`$, we return to the universal dependence

$$\gamma x\mathrm{sin}\theta .$$
(106)
Equation (54) determining the energy on the fast magnetosonic surface now has the form

$`\gamma ^32\sigma {\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}\left[2\epsilon f\mathrm{cos}\theta \epsilon \mathrm{sin}\theta {\displaystyle \frac{f}{\theta }}{\displaystyle \frac{1}{2}}{\displaystyle \frac{l_\mathrm{A}}{\sigma }}{\displaystyle _{x_0}^x}u(x^{})\gamma ^2(x^{})dx^{}+{\displaystyle \frac{1}{2x^2}}\right]\gamma ^2+{\displaystyle \frac{m_\mathrm{e}}{m_\mathrm{p}}}\sigma \mathrm{sin}^2\theta =0.`$ (107)
Hence, for arbitrary $`l_\mathrm{A}`$ the Lorentz factor of particles can be presented in the form

$$\gamma _\mathrm{F}=\left(\frac{2m_\mathrm{e}}{m_\mathrm{p}}\right)^{1/3}\sigma ^{1/3}\mathrm{sin}^{2/3}\theta .$$
(108)
Now, for $`l_\mathrm{A}<l_{\mathrm{cr}}^{(\mathrm{p})}`$ the position of the fast magnetosonic surface and the disturbance of magnetic surfaces on the fast magnetosonic surface are

$`x_\mathrm{F}`$ $``$ $`\left({\displaystyle \frac{2m_\mathrm{e}}{m_\mathrm{p}}}\right)^{1/3}\sigma ^{1/3}\mathrm{sin}^{1/3}\theta ,`$ (109)
$`(\epsilon f)_\mathrm{F}`$ $``$ $`\left({\displaystyle \frac{2m_\mathrm{e}}{m_\mathrm{p}}}\right)^{2/3}\sigma ^{2/3}.`$ (110)
Here for $`n<3`$ we have

$$l_{\mathrm{cr}}^{(\mathrm{p})}=\left(\frac{m_\mathrm{p}}{m_\mathrm{e}}\right)^{(5n)/3}\sigma ^{(n2)/3},$$
(111)
and for $`n3`$ we have $`l_{\mathrm{cr}}^{(\mathrm{p})}=(m_\mathrm{p}/m_\mathrm{e})^{2/3}\sigma ^{1/3}`$. On the other hand, for $`l_\mathrm{A}>l_{\mathrm{cr}}^{(\mathrm{p})}`$ we return to the relations $`x_\mathrm{F}=(\sigma /l_\mathrm{A})^{1/2}`$, $`(\epsilon f)_\mathrm{F}=l_\mathrm{A}/\sigma `$ for $`n3`$, but for $`n<3`$ now

$`x_\mathrm{F}`$ $``$ $`\left({\displaystyle \frac{\sigma }{l_\mathrm{A}}}\right)^{1/(5n)},`$ (112)
$`(\epsilon f)_\mathrm{F}`$ $``$ $`\left({\displaystyle \frac{l_\mathrm{A}}{\sigma }}\right)^{2/(5n)}.`$ (113)
As we see, for the disturbance of magnetic surfaces and the position of the fast magnetosonic surface remain the same as for the electron-positron outflow. In particular, the terminating compactness parameter $`l_{\mathrm{max}}=\sigma `$ is again determined by electrons.

Thus, one can conclude that inside the fast magnetosonic surface $`r<r_\mathrm{F}`$ the Lorentz factors of all particles are given by the universal relation (106). Hence, our results demonstrate that the drag force does not affect the structure of the flow inside the fast magnetosonic surface for both the electron-positron and the electron-proton cases. But the position of the fast magnetosonic surface and the critical value of the photon density depend on the mass of the outgoing particles, and for the electron-proton case they are determined by protons.

At large distances $`r>r_\mathrm{F}`$ for $`l_\mathrm{A}<l_{\mathrm{cr}}^{(\mathrm{p})}`$ the particle energy actually remains the same as on the fast magnetosonic surface ($`\gamma (r>r_\mathrm{F})\gamma _\mathrm{F}`$), while for $`l_\mathrm{A}>l_{\mathrm{cr}}^{(\mathrm{p})}`$ the Lorentz factor increases up to the maximum value

$$\gamma _{\mathrm{max}}\gamma _\mathrm{F}\left(\frac{l_\mathrm{A}}{l_{\mathrm{cr}}^{(\mathrm{p})}}\right)^{4(n2)/(5nn^2)}$$
(114)
for $`n<3`$. This takes place up to the distance $`rr_{\mathrm{cr}}`$, where

$$r_{\mathrm{cr}}=r_\mathrm{F}\left(\frac{l_\mathrm{A}}{l_{\mathrm{cr}}^{(\mathrm{p})}}\right)^{8/(5nn^2)}.$$
(115)
Clearly, this additional acceleration can be realized only for $`\sigma >m_\mathrm{p}/m_\mathrm{e}`$. Then, we have for $`n=3`$

$$_{\mathrm{sup}}\left(\frac{m_\mathrm{p}}{m_\mathrm{e}}\right)^{2/9}\sigma ^{7/9}m_\mathrm{e}c^2.$$
(116)
## 5 Discussion and Astrophysical Applications

We have demonstrated how for a simple geometry it is possible to determine a small radiation-drag-force correction to the one-fluid ideal MHD outflow. The disturbance of magnetic surfaces was self-consistently taken into consideration. As a result, it is possible to characterize the general influence of the drag action on the magnetic field structure for an ideal magnetically-dominated quasi-monopole cold outflow and to determine under what circumstances radiation drag is important.

As demonstrated above, the characteristics of the flow are determined by two main parameters, namely the compactness parameter $`l_\mathrm{A}`$ (42) (which is proportional to the photon density) and the magnetization parameter $`\sigma `$ (34). If the photon density is low, so that the compactness parameter is small $`l_\mathrm{A}l_{\mathrm{cr}}(\sigma )`$, the action of the drag force is negligible, while for a high photon density $`l_\mathrm{A}l_{\mathrm{cr}}(\sigma )`$, particles are additionally accelerated outside the fast magnetosonic surface.

In particular, for $`l_{\mathrm{cr}}l_\mathrm{A}l_{\mathrm{max}}`$ the increase of the drag force results in an increase of the outgoing plasma energy $`_{\mathrm{max}}\gamma _{\mathrm{max}}m_{\mathrm{e},\mathrm{p}}c^2`$, but the disturbance of magnetic surfaces is small ($`\epsilon f1`$). For $`l_\mathrm{A}l_{\mathrm{max}}`$ an increase of the photon density results in the increase of collimation up to values $`\epsilon f1`$, but the particle energy remains near the saturation value $`_{\mathrm{sup}}`$. Finally, for a very high photon density $`l_\mathrm{A}l_{\mathrm{max}}`$ an effective collimation of magnetic surfaces becomes possible, but in this case the drag force substantially diminishes the flux of electromagnetic energy inside the fast magnetosonic surface. As a result, for $`l_\mathrm{A}l_{\mathrm{max}}`$ almost all the energy of the electromagnetic field is lost via the inverse Compton interaction of particles with external photons. For this reason, the very existence of a magnetically-dominated flow becomes impossible. The dependence of the maximum particle energy $`_{\mathrm{max}}=\gamma _{\mathrm{max}}m_{\mathrm{e},\mathrm{p}}c^2`$ on the compactness parameter $`l_\mathrm{A}`$ is shown in Fig. 3.

We now consider several astrophysical applications.

### 5.1 Active Galactic Nuclei

For AGNs (the central engine is assumed to be a rotating black hole with mass $`M10^9M_{}`$, $`R10^{14}`$ cm, the total luminosity $`L10^{45}`$ erg s<sup>-1</sup>, $`B_010^4`$ G) the compactness parameter $`l_\mathrm{A}`$ (42) can be evaluated as

$$l_\mathrm{A}30M_9^1\left(\frac{\mathrm{\Omega }R}{c}\right)L_{45}.$$
(117)
In the Michel magnetization parameter $`\sigma `$ (34)

$$\sigma 10^{14}\lambda ^1M_9B_4\left(\frac{\mathrm{\Omega }R}{c}\right),$$
(118)
the main uncertainty comes from the multiplication parameter $`\lambda `$, i.e., in the particle number density $`n`$. Indeed, for an electron-positron outflow this value depends on the efficiency of pair creation in the magnetosphere of a black hole, which is still undetermined. In particular, this process depends on the density and energies of the photons in the immediate vicinity of the black hole. As a result, if the hard-photon density is not high, then the multiplication parameter is small ($`\lambda 10100`$; Beskin, Istomin & Pariev 1992; Hirotani & Okamoto 1998). In this case for $`(\mathrm{\Omega }R/c)0.1`$$`0.01`$ we have $`\sigma 10^910^{12}`$, so that $`l_{\mathrm{cr}}10^3`$$`10^4`$. On the other hand, if the density of photons with energies $`_\gamma >1`$MeV is high enough, direct particle creation $`\gamma +\gamma e^++e^{}`$ results in an increase of the particle density (Svensson 1984). This gives $`\sigma 1010^3`$, and hence $`l_{\mathrm{cr}}10`$ for an electron-positron outflow.

From a theoretical point of view, the most interesting result here is the possibility of an additional acceleration of particles outside the fast magnetosonic surface. Indeed, for a high enough photon density ($`l_\mathrm{A}10100`$, i.e., for $`L10^{46}10^{48}`$ erg s<sup>-1</sup>) and a small magnetization parameter $`\sigma 10100`$, the compactness parameter $`l_\mathrm{A}`$ can exceed the critical value $`l_{\mathrm{cr}}`$ for an electron-positron outflow. In this case, according to Fig. $`1b`$, our analysis suggests that the kinetic luminosity of the relativistic jet should be proportional to $`l_\mathrm{A}^{2/3}L_{\mathrm{tot}}`$, where $`L_{\mathrm{tot}}`$ in the total luminosity of the central engine. Kinetic luminosity is not easily determined from observations. However, observational evidence suggests that the radio luminosity of the jets is positively correlated with the luminosity of the central engine and the scatter of this correlation decreases towards larger luminosities (Baum, Zirbel & O’Dea 1995).

For an electron-proton outflow the magnetization parameter (34) can be rewritten in the form (Camenzind 1990)

$$\sigma =\frac{m_\mathrm{p}}{m_\mathrm{e}}\left(\frac{\mathrm{\Omega }R}{c}\right)^2\frac{B_0^2R^2}{c\dot{M}}3\times 10^4\left(\frac{\mathrm{\Omega }R}{c}\right)^2B_4^2M_9^2\left(\frac{\dot{M}}{0.1M_{}/\mathrm{yr}}\right)^1.$$
(119)
Here $`\dot{M}=4\pi nm_\mathrm{p}R^2c`$ is the mass ejection rate. Hence, for a high ejection rate ($`\dot{M}>0.1M_{}`$ yr <sup>-1</sup>) the magnetization parameter $`\sigma <m_\mathrm{p}/m_\mathrm{e}`$. In this case there is no acceleration of plasma. On the other hand, for low ejection rate $`\dot{M}<0.1M_{}/\mathrm{yr}`$ the magnetization parameter becomes too large for the drag force to be efficient.

Thus, the drag force can substantially disturb the MHD parameters of a Poynting-dominated outflow only for a very high luminosity of the central engine ($`L_{\mathrm{tot}}10^{45}`$ erg s <sup>-1</sup>) and only for an electron-positron outflow. In all other cases the action of the drag force remains negligible. In particular, the additional acceleration of particles outside the fast magnetosonic surface is not efficient.

### 5.2 Cosmological Gamma-Ray Bursts

For cosmological gamma-ray bursts (the central engine is represented by the merger of very rapidly orbiting neutron stars or black holes with $`MM_{}`$, $`R10^6`$ cm, total luminosity $`L10^{52}`$ erg s<sup>-1</sup>, $`B_010^{15}`$ G; see, e.g., Lee et al 2000 for details) the compactness parameter $`l_\mathrm{A}`$ is extremely large:

$$l_\mathrm{A}10^{17}\left(\frac{\mathrm{\Omega }R}{c}\right)L_{52}.$$
(120)
On the other hand, even for a superstrong magnetic field of $`B_010^{15}`$ G (which is necessary to explain the total energy release) the magnetization parameter $`\sigma `$ is small ($`\sigma <110`$), because within this model the magnetic field itself is secondary and its energy density cannot exceed the plasma energy density. Thus, one can conclude that for these characteristics of cosmological gamma-ray bursts the density of photons is very high so that $`l_\mathrm{A}l_{\mathrm{max}}`$ and the drag force can make it difficult to form a Poynting-dominated outflow. A self-consistent analysis should include into consideration other physical processes such as high optical thickness resulting in the diminishing of the photon density, radiation and particle pressure, etc. Nevertheless, in our opinion, our conclusion may substantially restrict some recent models of cosmological gamma-ray bursts.

### 5.3 Radio Pulsars

For radio pulsars the central engine is a rotating neutron star with $`MM_{}`$, $`R10^6`$ cm, total luminosity of the surface $`L_\mathrm{X}10^{33}`$$`10^{37}`$ erg s<sup>-1</sup>, and $`B_010^{12}`$ G. In this case the magnetization parameter $`\sigma 10^4`$$`10^6`$, corresponding to relativistic electron-positron plasma, is known with rather high accuracy (see, e.g., Bogovalov 1997). This gives $`l_{\mathrm{cr}}10^2`$$`10^3`$, and the compactness parameter

$$l_\mathrm{A}\left(\frac{\mathrm{\Omega }R}{c}\right)L_{35}$$
(121)
remains small ($`<1`$) even for the most energetic ($`L_\mathrm{X}10^{37}`$ erg s <sup>-1</sup> ) fast ($`\mathrm{\Omega }R/c10^2`$) pulsars like Crab and Vela. Thus, one can conclude that the drag force does not substantially disturb the magnetically-dominated outflow from radio pulsars.

Thus, the drag force does not affect the wind characteristics (particle energy, magnetic field structure, etc.) of pulsars. However, interaction of outflowing relativistic particles with thermal photons can be important in other ways. In the wind region ($`rR_\mathrm{L}`$) even a weak interaction with photons can result in a detectable flux of inverse Compton gamma-ray photons (Bogovalov & Aharonian 2000). On the other hand, near the surface of the star ($`rR_\mathrm{L}`$), inverse Compton photons are important in the pair creation process (Kardashev, Mitrofanov & Novikov 1984; Zhang & Harding 2000).

## Acknowledgments

The authors are grateful to K.A.Postnov for fruitful discussions. We thank the referee for a detailed report and useful comments. VSB thanks Observatory Paris-Meudon for hospitality. This work was partially supported by Russian Foundation for Fundamental Research (Grant 02-02-16762).


format: text

text:
# Cooperative behavior in a spatial model of “commons”

## 1 Introduction

Selfishness is regarded as the rational behavior of agents in various social and economic studies. However, there are cases where altruistic or cooperative tendency can survive. The game of “prisoners dilemma” is one such case where cooperative behavior is necessary to maximize ones payoff . Axelrod and Hamilton have shown that selfish agents with finite memory will eventually learn to cooperate if the game is played many times . In a later work, May and Nowak showed that in a spatially extended population, altruistic or cooperative tendency can survive by forming spatial territories. A single altruist is not strong but it becomes strong in a group because of its cooperative tendency. Thus once a domain of altruists is formed in space it can survive against the selfish members of the species. In some biological systems altruistic behavior is observed instead of purely selfish behavior .

The problem of “tragedy of commons” is the other such case where the selfish tendency to maximize one’s gain leads to a degeneration of the common resource, and thus reduces everyones gain . We study a spatial version of this problem in which we assume a resource spread over all the sites of the lattice. If an uncontrolled exploitation by a few agents leads to a local degeneration of the resource; its negative effect is felt only in an immediate neighborhood. Similarly, if the resource is maintained in a good condition, its positive effect is also felt locally. Thus, the resource is shared only locally. This is the situation observed in case of most natural resources.

We model this situation with evolution rules similar to those of May and Nowak. In our model, instead of only two types of agents, there are many types of agents whose cooperative tendencies are labeled by a real variable between $`0`$ and $`1`$. The agents of type $`1`$ are totally cooperative and the agents of type $`0`$ are totally selfish, and the agents of type between $`0`$ and $`1`$, represent a mixed strategy. In the steady state, the density profile of agents develops two peaks, one near $`0`$ representing the defector band and the other near $`1`$, representing the cooperator band. Agents from these two bands form a complex spatial pattern (see Fig. 1). Members of the defector band form thin lines whereas members of the cooperator band form compact domains, which we will also call colonies. It should be emphasized that the colonies are formed only by cooperators, demonstrating the advantage of cooperative over selfish behavior. This pattern is similar to the one found by May and Nowak in certain range of the parameter used in their model.

The complex pattern formed in this model has some interesting statistical properties. If we allow mutation, i.e., a change of strategy, of agents at a moderate rate, then we find that the model self-organizes to a critical steady state. The power spectrum of fluctuation of the activity has a $`1/f^\alpha `$ form with $`\alpha =1.30`$. The critical behavior is lost for very slow mutation rate. Thus this case of self-organized criticality is different from that of the sandpile model which is critical at very slow drive rate .

The moderate mutation rate helps evolution of cooperation in the initial stage. But, the large critical fluctuations restrict the cooperative tendency to a maximum of $`8/13`$. At very slow mutation rate an interesting process of evolution starts. A more stable steady state helps the growth of cooperative tendency higher than $`8/13`$.

## 2 “Tragedy of commons”

Consider a resource shared by a group of $`N`$ agents. When the exploitation rate $`X`$ is very slow compared to the regeneration rate of the resource, the yield $`P(X)`$ is proportional to $`X`$. However, when $`X`$ becomes approximately equal to the regeneration rate the yield $`P(X)`$ reaches a saturation point. If the resource is exploited at an even faster rate, it degenerates, and $`P(X)`$ decreases. Thus the yield $`P(X)`$, at the exploitation rate $`X`$ has the form shown in Fig. 2. $`P(X)`$ is maximum at $`X_c`$, and goes to zero at $`X_m`$ (the resource is fully destroyed at $`X_m`$). For definiteness, we assume the form $`P(X)=X(2\frac{X}{N})`$. So $`X_c=N`$, and $`X_m=2N`$. The income of the $`i`$th agent $`P_i(x_i)`$ is proportional to his exploitation rate $`x_i`$.

$$P_i(x_i)=\frac{x_i}{X}P(X)=x_i(2\frac{X}{N}),$$
(1)
where $`x_i=X`$.

Consider the case when all the agents operate at a rate $`x=1`$, the income of each of them is $`1`$, and net income of the group is maximum. Now, let the $`jth`$ agent increases his exploitation rate from $`1`$ to $`2`$. Then $`X=N+1`$, $`P_i=11/N`$, for $`ij`$, and $`P_j=2(11/N)`$. Thus the $`j`$th agent almost doubles his income, while the other members incur a small loss of order $`1/N`$. The net income of the group also decreases by a fraction of order $`1/N^2`$. Therefore, if the agents try to maximize their personal income, the system is pushed beyond the optimum exploitation rate $`X_c`$, towards $`X_m`$, where the income of all agents goes to zero, and the “tragedy of commons” occurs .

One observes that in this case, there is a conflict between individual interest and global interest ( and eventually the individual interest has conflict with itself). Thus if we try to maximize the income of each agent by the set of equations of form $`\frac{dx_i}{dt}=\mathrm{\Gamma }\frac{dP_i(x_i)}{dx_i}+\eta `$, we find that this set of equations, instead of maximizing the $`P_i`$’s, lead to a state of lower $`P_i`$’s. Hence, the separate optimization of the individual incomes, does not work.

In the physical systems this kind of a situation can be seen in case of granular media. If sand is allowed to flow in presence of the gravitational field, the individual tendency of grains to minimize their potential energy may lead to a jammed configuration. Thus individual tendency has conflict with itself.

Here, we try to understand that how cooperation can evolve so that the “tragedy of commons” can be averted.

## 3 Model of locally shared resource (LSR)

We study a model of locally shared resource (LSR). This is defined on a square lattice, with a real variable $`s_i[0,1]`$, called the cooperative level of the occupant of the $`i`$th site. This represents the cooperative tendency of an agent. The exploitation rate depends on level $`s_i`$ through, $`x_i=2s_i`$ ( clearly $`x_i[1,2]`$). The higher $`s_i`$ implies higher cooperative tendency and lower $`x_i`$. The agents of $`s<\frac{1}{2}`$, will be called “defectors” or “selfish”, and those of $`s\frac{1}{2}`$ will be called “cooperators” or “altruists”. The resource near the site $`i`$ is shared in a neighborhood $`D(i)`$, which consists of the occupants of the $`i`$th site and its nearest and next nearest neighbors. Thus $`N=9`$ in two dimensions. The payoff of the occupant of the $`i`$th site can be expressed in terms of $`s_i`$ using (1) and the relation $`x_i=2s_i`$.

$$P_i=\frac{1}{N}(2s_i)\underset{jD(i)}{}s_j.$$
(2)
The evolution rule consists of two parts:

1. Mutual struggle: After one time step, the site $`i`$ will be occupied by the neighbor of highest payoff, i.e., $`s_is_h`$, where $`hD(i)`$ has highest income $`P_h`$.
2. Mutation: The levels of a small fraction $`p`$ of agents are changed to randomly chosen levels.

We study this model with both periodic and fixed boundary conditions. In case of fixed boundary condition, the set $`D(i)`$ has only four elements for the sites on the corners and it has six elements for the sites on the edges. In case of periodic boundary condition, $`D(i)`$ contains nine elements for all $`i`$.

The first rule of evolution in our model is same as that in the model of May and Nowak. However in their model, there are only two types of agents: One who always cooperates ($`s=1`$) and the other who always defects ($`s=0`$). The game of prisoners’ dilemma is played between two adjacent neighbors. If both cooperate then the payoff of each of them is 1. If one defects while the other cooperates then the payoff of the defector is $`b>1`$, while the payoff of the cooperator is 0. If both defect then the payoff of each of them is 0. If agent of the $`i`$th site plays this game with all his eight neighbors, then his payoff is

$$P_i=[b(b1)s_i]\underset{jD(i)i}{}s_j.$$
(3)
In the LSR model, there are infinitely many types of agents labeled by a continuous variable $`s`$. If we consider a population of only two levels $`s^1`$, and $`s^2`$, then this model is similar to the model of May and Nowak. Note that the two models are not identical because the summation in (2) for the LSR model includes $`i`$. The effective $`b`$ can be defined as the ratio of payoffs of $`s^1`$ and $`s^2`$ in identical neighboring configurations. Ignoring $`O(1/N)`$ correction, we can write $`b(s^1,s^2)(2s^1)/(2s^2)`$. However, it is important to point out that $`b`$ alone does not determine the behavior of the model.

## 4 Density profile in LSR model

We simulate this model for different mutation rates on lattices of various sizes. The statistical behavior of the model does not depend on the boundary condition. However, it depends on the mutation rate. For moderate mutation, i.e. $`p>0.0005`$ approximately, the model shows critical behavior. In this case, the average cooperative level $`s0.50`$. In case of very slow mutation rate, i.e. $`p<0.0005`$ approximately, the model slowly develops higher cooperative regions, and the average cooperative level rises to $`s0.53`$.

In both cases the highest cooperative members (i.e. $`s=1`$) do not survive, but a little practical agents from a band of lower cooperative levels, i.e. $`\frac{8}{13}<s<\frac{5}{7}`$, dominate in the steady state (see the second peak in the density plot, Fig. 3, and second step in the cumulative plot, Fig. 4.). We call it the cooperator band (C-band). This band constitutes around $`71\%`$ of the population. In case of a moderate mutation rate this band shrinks to just one level $`s=\frac{8}{13}`$. In case of a slow mutation rate the higher levels also form colonies. As a result the population of cooperators is evenly distributed in range $`\frac{8}{13}s\frac{5}{7}`$.

Around $`25\%`$ of the population come from the defector band (D-band), corresponding to levels, $`0.1<s<0.3`$ approximately (see the first peak in density plot, Fig. 3, and the first step in the cumulative plot, Fig. 4.). The level lying between these two bands, i.e. $`0.3<s<0.6`$ approximately, are almost absent from the population. Thus we find that in steady state the density profile of cooperative levels develops two peaks.

The C-band is very sharp, while the D-band is wide and contains more levels. This is because the spatially close levels have very high mutual competition. As a result the levels of C-band eliminate each other and in the end only a few levels survive. The levels of D-band have low density. They exist as scattered separated lines. A defector of one level hardly encounters a defector of another level. Therefore, they are not able to eliminate each other.

A good deal of understanding about the density profile of agents of various levels can be obtained by considering the stability of some of the special cases.

Edge stability: Consider the domains of two types of agents $`s`$ and $`s^{}`$, with $`s>s^{}`$, separated along a straight line (see Fig. 5). The agents of level $`s`$ is represented by white color, and the agents of level $`s^{}`$ is represented by gray color. The agent $`1`$ is among the strongest white agent which can have influence on the domain boundary. The agent $`2`$ is the strongest white agent near the domain boundary. And, the agent $`3`$ is the strongest gray agent near the domain boundary. The payoffs of these three agents are, $`P_1=(2s)s`$, $`P_2=\frac{1}{9}(2s)(6s+3s^{})`$, and $`P_3=\frac{1}{9}(2s^{})(3s+6s^{})`$. Now there are three possibilities:

* $`P_2>P_3`$, i.e. $`s^{}<1s`$. In this case the white agent $`2`$ can take over the gray agent $`3`$, and the edge will move towards the domain of the level $`s^{}`$. Hence, the level $`s`$ will be called edge-strong against the level $`s^{}`$. The levels $`s<\frac{1}{2}`$ are edge-strong against all the lower levels.
* $`P_1<P_3`$, i.e. $`2\frac{3}{2}s<s^{}<s`$. In this case the gray agent $`3`$ can take over his neighboring white agents including agent $`2`$, and the edge will move towards the domain of the level $`s`$, so the level $`s`$ is edge-weak against the level $`s^{}`$. The levels $`s<\frac{4}{5}`$ are never edge-weak against any of the lower levels.
* $`P_2<P_3`$, but $`P_1>P_3`$, i.e. $`1s<s^{}<2\frac{3}{2}s`$. In this case the edge will be stable, so the level $`s`$ will be called edge-stable against the level $`s^{}`$.

Corner stability: Consider now the case in which a domain of level $`s^{}`$ forms a corner in a domain of level $`s`$, and $`s>s^{}`$ (see Fig. 5). Here again we select three agents using the same criterion as in the previous case. The payoffs of these three agents are, $`P_1=(2s)s`$, $`P_2=\frac{1}{9}(2s)(8s+s^{})`$, and $`P_3=\frac{1}{9}(2s^{})(5s+4s^{})`$. Again there are three possibilities:

* $`P_2>P_3`$, i.e. $`s^{}<\frac{1}{2}(34s)`$. In this case the corner of $`s^{}`$ will be decimated. Hence, the level $`s`$ will be called corner-strong against the level $`s^{}`$. The levels $`s<\frac{1}{2}`$ are corner-strong against all the lower levels, and the levels $`s>\frac{3}{4}`$ are never corner-strong against any of the lower levels .
* $`P_1<P_3`$, i.e. $`2\frac{9}{4}s<s^{}<s`$. The domain of $`s^{}`$ will grow near the corner, so the level $`s`$ is corner-weak against the level $`s^{}`$. The levels $`s<\frac{8}{13}`$ are never corner-weak against any of the lower level.
* $`P_2<P_3`$, but $`P_1>P_3`$, i.e. $`\frac{1}{2}(34s)<s^{}<2\frac{9}{4}s`$. In this case corner will be stable, so the level $`s`$ will be called corner-stable against the level $`s^{}`$.

Step stability: Consider now the case in which a domain of level $`s^{}`$ forms a step in a domain of level $`s`$, and $`s>s^{}`$. The payoffs of the three agents marked in Fig. 5 are, $`P_1=(2s)s`$, $`P_2=\frac{1}{9}(2s)(8s+s^{})`$, and $`P_3=\frac{1}{9}(2s^{})(4s+5s^{})`$. Again there are three possibilities:

* $`P_2>P_3`$, i.e. $`s^{}<\frac{8}{5}(1s)`$. In this case the agent $`3`$ will be taken over by the agent $`2`$. Thus the domain of level $`s`$ will grow along the step. The level $`s`$ will be called step-strong against the level $`s^{}`$. The levels $`s<\frac{8}{13}`$ are step-strong against all the lower levels.
* $`P_1<P_3`$, i.e. $`2\frac{9}{5}s<s^{}<s`$. The domain of $`s^{}`$ will grow, so the level $`s`$ is step-weak against the level $`s^{}`$. The levels $`s<\frac{5}{7}`$ are never step-weak against any of the lower level.
* $`P_2<P_3`$, but $`P_1>P_3`$, i.e. $`\frac{1}{2}(34s)<s^{}<2\frac{9}{4}s`$. In this case step will be stable, so the level $`s`$ will be called step-stable against the level $`s^{}`$.

From this simplistic analysis of corner edge and step stability, one can understand the density profile shown in Fig. 4.

The levels $`s<\frac{1}{2}`$ do not have any threat from the lower levels, but they are edge-weak and corner-weak against the higher levels. Therfore they cannot form colonies and survive just as thin lines.

The levels $`\frac{1}{2}<s<\frac{8}{13}`$ are not edge-weak or corner-weak against any level. Therefore they form colonies in the initial stage. However, because of their step-weakness, the colonies are slowly decimated by the higher levels.

The level $`s=\frac{8}{13}`$ is not weak against any level. Therefore, in case of moderate mutation rate there is a huge jump in cumulative plot at this level (see Fig. 4). It constitutes $`71\%`$ of the whole population.

The levels $`\frac{8}{13}<s<\frac{5}{7}`$ are not edge-weak or step-weak against any level. But they are corner-weak against some of the lower levels. In case of moderate mutation rate these levels are decimated by the level $`s=\frac{8}{13}`$ (so there is a huge peak at $`s=\frac{8}{13}`$). However, in case of very slow mutation rate these levels can form some colonies by forming a protective coat of defectors near the corners. Therefore, in this case the population of cooperators is evenly distributed in the whole range (see Fig. 4, the second step is broadened).

The levels $`s>\frac{5}{7}`$ are very weak against the lower levels. They are not even stable as thin lines. Therefore they are completely wiped out from the population. Their high cooperative tendency works against their survival.

## 5 Moderate mutation rate

The behavior of this model depends on the mutation rates. For moderate mutation rate, i.e., $`p>0.0005`$ approximately, the model shows critical fluctuation in time. We study the power spectrum of the number of active sites (sites changing their levels) with both fixed and periodic boundary conditions. We find that the power spectrum does not depend on the boundary condition and it shows $`1/f^\alpha `$ behavior, with $`\alpha 1.30`$ (see Fig. 6, and 7).

The mutation of the occupant of a site creates an avalanche. In case of the moderate mutation rate, the system does not get enough time to relax, and these avalanches overlap with each other leading to a bigger avalanche. The overlap of these avalanches cannot be regarded as uncorrelated superposition of smaller avalanches because of the temporal correlation. These big avalanches are responsible for the nontrivial behavior of the power spectrum. In case of the slow mutation rate, i.e., $`p<0.0005`$, the system gets enough time to relax, therefore the avalanches do not overlap. In this case, the power spectrum has $`1/f^2`$ behavior. We note that self-organized criticality observed in this model, in case of moderate mutation rate, is different from that of the sandpile model, in which the avalanches do not overlap because of the slow drive rate .

The “$`1/f`$” form of the power spectrum alone does not imply that the model has long range temporal correlation. This kind of behavior is seen in models with no translational invariance . In this case the autocorrelation may have exponential form with different relaxation time at different point in space. So the power spectrum at each point in space has the form $`1/f^2`$, i.e., the power spectrum of a random walk. Averaging over space, one obtains a power spectrum of form $`1/f^\alpha `$, with a nontrivial value of $`\alpha `$, which is less than $`2`$. By this mechanism even a noncritical models can show a $`1/f`$ noise . Thus $`0<\alpha <2`$, does not imply long rang autocorrelation.

In our model there is a spatial pattern (Fig. 1) which could lead to spatially varying relaxation time. However, this pattern slowly changes in time because of the random mutations. To check for the spatially varying relaxation time, we study the waiting time (i.e., the time taken between two consecutive changes of level), with both types of boundary conditions. For fixed boundary condition we select three different sites on lattice. The first site is near the center, the second is near a boundary but away from a corner, and the third is near the a corner. The frequency distribution of the waiting time in the first two cases show a fairly good power law behavior in two decades (see Fig. 8). It can be fitted to the form: $`n(t)t^\beta `$, with $`\beta =1.25`$. At the site near a corner, the waiting time is large and does not fit well with a power law. However, the contribution to the power spectrum from such a sites is negligible. For periodic boundary condition, the waiting time distribution is similar to the one obtained for the site in the bulk with fixed boundary condition (Fig. 9). Thus in both cases, the relaxation time does not significantly depend on the position. Thus the “$`1/f`$” behavior of the power spectrum is due to the temporal correlation of the steady state dynamics.

We have found similar behavior in the power spectrum of the average cooperative level $`S`$. Thus the exponent $`\alpha `$ of the power spectrum appears to be universally in this model.

Our model has many periodically changing subconfigurations in the steady state. For this reason, analysis in terms of avalanches produced by small local perturbations is forbiddingly difficult. Therefore we look for a different method to investigate the dynamical behavior. The analysis of the waiting time gives an alternative way to explore the critical dynamics of the model.

## 6 Slow mutation rate

In case of a slow mutation rate, i.e., $`p<0.0005`$ approximately, the model slowly evolves colonies of higher cooperation level, i.e., $`s>8/13`$ (see Fig. 3 and 4. The defectors play an important role in protecting these colonies against the raid of the level $`s=8/13`$. The concave part (corner) in the boundary of the colony of higher cooperators is weak against the immediate lower levels. However, one defector at the corner makes the colony of higher cooperators stable.

The higher cooperative levels are not able to grow in the initial stage. They are mainly threatened by the immediate lower levels, rather than the defectors. Slowly a small domain of higher level is formed by mutation near the boundary of the colony of lower cooperators. This process is normally very slow and may take hundreds of time steps. The defectors in the neighborhood of these cooperators become strong and may advance toward the colony of $`s=8/13`$, whenever there is a fluctuation at the boundary because of mutation. Behind the line of these defectors the colony of higher cooperators grows. The growth is a very slow process and takes tens of thousands of time steps.

The steady state in this case is not critical. The power spectrum of fluctuation of activity has a $`1/f^2`$ form, which is same as the power spectrum of a random walk.

## 7 Summary

We have explored the possibility of survival of cooperative tendency in a model of locally shared common resources. The spatial extension of the model allows the cooperators to form stable colonies (domains). The selfish agents can only survive in the neighborhood of cooperators, therefore they are found on the boundaries between these colonies. Thus, natural selection has provision for coexistence of both type of behaviors.

In addition, this model shows interesting statistical behavior. A moderate mutation rate makes the model critical with fluctuation of all time scales. The critical fluctuations are helpful in evolution of cooperative behavior in the initial stage. However, it work against the evolution of higher cooperative level. In case of slow mutation rate the system becomes noncritical and an interesting process of evolution starts. The higher cooperators form a coat of defectors and slowly advance in the territory of lower cooperators. Thus the average cooperative level slowly increases.