text:
# Energy and Momentum of a Class of Rotating Gravitational Waves

## I Introduction

The notion of energy has been one of the most thorny and important problems in Einstein’s theory of General Relativity (GR). There have been many attempts  to get a well defined expression for local or quasi-local energy and momentum. However, there is still no generally accepted definition known. As a result, different people have different points of view. Cooperstock  argued that in GR, energy and momentum are localized in regions of the non-vanishing energy and momentum tensor and consequently gravitational waves are not carriers of energy and momentum in vacuum. The gravitational waves, by definition, have zero stress-energy tensor. Thus the existence of these waves was questioned. However, the theory of GR indicates the existence of gravitational waves as solutions of Einstein’s field equations . Infact this problem arises because energy is not well defined in GR.

The problem for gravitational waves was resolved by Ehlers and Kundt , Pirani  and Weber and Wheeler  by considering a sphere of test particles in the path of the waves. They showed that these particles acquired a constant momentum from the waves. Qadir and Sharif  presented an operational procedure, embodying the same principle, to show that gravitational waves impart a momentum. Rosen  investigated whether or not cylindrical gravitational waves have energy and momentum. He used the energy-momentum pseudo tensors of Einstein and Landau Lifshitz and carried out calculations in cylindrical polar coordinates. However, he arrived at the conclusion that the energy and momentum density components vanish. These results supported the conjecture of Scheidegger  that physical system cannot radiate gravitational energy. Later, he pointed out  that the energy and momentum densities turn out to be non-vanishing and reasonable if the calculations are performed in Cartesian coordinates. Rosen and Virbhadra  explicitly evaluated these quantities in the Einstein’s prescription by using Cartesian coordinates and found them finite and well defined. Virbhadra  then used Tolman, Landau-Lifshitz and Papapetrou’s prescriptions to evaluate the energy and momentum densities and found that the same results turn out in all these prescriptions.

In this paper we use Einstein and Papapetrou’s prescriptions to evaluate energy and momentum densities to a class of cylindrical rotating gravitational waves. As we shall see from the analysis given in the paper, when rotation is included, the problem is considerably complicated. We find that the results obtained by these two prescriptions are not exactly the same. However, it is shown that they both reduce to the same result for a special case of cylindrical gravitational waves. In the next section, we shall describe the class of rotating cylindrical gravitational waves. In sections three and four, we evaluate energy and momentum using Einstein and Papapetrou’s prescriptions respectively. Finally, we shall discuss the results.

## II A Class of Rotating Cylindrical Gravitational Waves

A class of solutions of the gravitational field equations describing vacuum spacetimes outside rotating cylindrical sources is given by the line element of the form

$$ds^2=e^{2\gamma 2\psi }(dt^2d\rho ^2)\mu ^2e^{2\psi }(\omega dt+d\varphi )^2e^{2\psi }dz^2,$$
(1)
in the cylindrical coordinates $`(\rho ,\varphi ,z)`$. Here the metric functions $`\gamma ,\mu ,\psi `$ and $`\omega `$ depend on the coordinates $`t`$ and $`\rho `$ only. In general the spacetimes have two Killing vectors, one is associated with the invariant translations along the symmetry axis, $`\xi _{(z)}=z`$, and the other is associated with the invariant rotations about the axis, $`\xi _{(\varphi )}=\varphi `$. Obviously, these Killing vectors are orthogonal. When $`\omega =0`$, the metric represents spacetimes without rotation, in which the polarization of gravitational waves has only one degree of freedom and the direction of polarization is fixed . It is to be noticed that if we take $`\omega =0`$ and $`\mu =\rho `$, the above metric reduces to a special case of cylindrical gravitational waves . Einstein’s vacuum field equations for the metric form (1) are given by

$$(\mu \psi _v)_u+(\mu \psi _u)_v=0,$$
(2)
$$\mu _{uv}\frac{l^2}{8}\mu ^3e^{2\gamma }=0,$$
(3)
$$\omega _v\omega _u=l\mu ^3e^{2\gamma },$$
(4)
$$\gamma _u=\frac{1}{2\mu _u}(\mu _{uu}+2\mu \psi _u^2),$$
(5)
$$\gamma _v=\frac{1}{2\mu _v}(\mu _{vv}+2\mu \psi _v^2),$$
(6)
where $`\psi _u=\frac{\psi }{u}`$, etc. The subscripts $`u=t\rho `$ and $`v=t+\rho `$ are retarded and advanced times respectively. Here $`l`$ is a constant length characteristic of the rotation of the system which is positive and is specifically attributed with rotating gravitational waves. For $`l=0`$ we have $`\omega =\omega (t)`$ from Eq.(4) and $`\mu _{tt}=\mu _{\rho \rho }`$ from Eq.(3). A simple transformation to a rotating frame reduces the waves to non-rotating generalized Beck spacetimes which have been studied by many authors .

In order to have meaningful results in the prescription of Einstein and Papapetrou, it is necessary to transform the metric in Cartesian coordinates. Let us transform the metric in Cartesian coordinates by using

$$x=\rho \mathrm{cos}\varphi ,y=\rho \mathrm{sin}\varphi .$$
(7)
The corresponding metric in these coordinates will become

$$ds^2=e^{2\gamma 2\psi }(dt^2\frac{1}{\rho ^2}(xdx+ydy)^2)\mu ^2e^{2\psi }(\omega dt+\frac{1}{\rho ^2}(xdyydx))^2e^{2\psi }dz^2.$$
(8)
## III Energy and Momentum in Einstein’s Prescription

The energy-momentum complex of Einstein  is given by

$$\mathrm{\Theta }_a^b=\frac{1}{16\pi }H_{a,c}^{bc},$$
(9)
where

$$H_a^{bc}=\frac{g_{ad}}{\sqrt{g}}[g(g^{bd}g^{ce}g^{cd}g^{be})]_{,e},a,b,c,d,e=0,1,2,3.$$
(10)
$`\mathrm{\Theta }_0^0`$ is the energy density, $`\mathrm{\Theta }_0^a`$ are the momentum density components, and $`\mathrm{\Theta }_a^0`$ are the components of energy current density. The Einstein energy-momentum satisfies the local conservation laws

$$\frac{\mathrm{\Theta }_a^b}{x^b}=0.$$
(11)
The required non-vanishing components of $`H_a^{bc}`$ are

$$H_0^{01}=\frac{1}{\mu \rho ^3}(\mu ^2x2\mu \mu _\rho \rho x+\rho ^2xe^{2\gamma }\mu ^2\omega ^2\rho ^2x+2\mu ^2\omega \dot{\gamma }\rho ^2y\mu ^4\omega \omega _\rho \rho xe^{2\gamma }),$$
(12)
$$H_0^{02}=\frac{1}{\mu \rho ^3}(\mu ^2y2\mu \mu _\rho \rho y+\rho ^2ye^{2\gamma }\mu ^2\omega ^2\rho ^2y2\mu ^2\omega \dot{\gamma }\rho ^2x\mu ^4\omega \omega _\rho \rho ye^{2\gamma }),$$
(13)
$$H_1^{01}=\frac{1}{\rho ^4}(2\mu \dot{\gamma }\rho y^2+2\dot{\mu }\rho x^2\mu ^3\omega _\rho xye^{2\gamma }),$$
(14)
$$H_1^{02}=\frac{1}{\rho ^4}(\mu \omega \rho ^3+2\mu \dot{\gamma }\rho xy+\mu ^3\omega _\rho y^2e^{2\gamma }),$$
(15)
$$H_2^{01}=\frac{1}{\rho ^4}(2\mu \dot{\gamma }\rho xy2\dot{\mu }\rho xy\mu \omega \rho ^3\mu ^3\omega _\rho x^2e^{2\gamma }),$$
(16)
$$H_2^{02}=\frac{1}{\rho ^4}(2\mu \dot{\gamma }\rho x^2+2\dot{\mu }\rho y^2+\mu ^3\omega _\rho xye^{2\gamma }),$$
(17)
$$H_3^{03}=\frac{2}{\rho }(\mu \dot{\gamma }2\mu \dot{\psi }+\dot{\mu }),$$
(18)
$$H_0^{12}=\mu \omega _\rho +2\mu \omega \gamma _\rho 2\mu _\rho \omega ,$$
(19)
$$H_0^{03},H_1^{03},H_2^{03},H_3^{01},H_3^{02},H_0^{23},H_0^{31}=0.$$
(20)
Using Eqs.(12-20) in Eq.(9), we obtain energy and momentum densities in Einstein’s prescription

$`\mathrm{\Theta }_0^0={\displaystyle \frac{1}{16\pi \mu ^2\rho ^3}}[\mu ^2(\mu +\mu _\rho \rho 2\mu _{\rho \rho }\rho ^2\mu \omega ^2\rho ^22\mu \omega \omega _\rho \rho ^3\mu _\rho \omega ^2\rho ^3)`$ (21)
$`+\rho ^2(\mu +2\mu \gamma _\rho \rho \mu _\rho \rho )e^{2\gamma }\mu ^4\rho ^2(\mu \omega _\rho ^2+\mu \omega \omega _{\rho \rho }`$ (22)
$`2\mu \omega \omega _\rho \gamma _\rho +3\mu _\rho \omega \omega _\rho )e^{2\gamma }],`$ (23)
$`\mathrm{\Theta }_1^0={\displaystyle \frac{1}{16\pi \rho ^5}}[2\mu \rho ^2\dot{\gamma }x6\dot{\mu }x^32\dot{\mu }\rho ^2x2\dot{\mu _\rho }\rho x^3\mu \omega ^2\rho y`$ (24)
$`\mu ^2\rho y(\mu \omega _\rho \mu \omega _{\rho \rho }\rho +2\mu \omega _\rho \gamma _\rho \rho 3\mu \mu _\rho \omega _\rho \rho )e^{2\gamma }],`$ (25)
$`\mathrm{\Theta }_2^0={\displaystyle \frac{1}{16\pi \rho ^4}}[2\mu \rho \dot{\gamma }y2\dot{\mu _\rho }\rho ^2y+\mu \omega \rho x\mu \omega _\rho \rho ^2x\mu _\rho \omega \rho ^2x`$ (26)
$`+\mu ^2x(\mu \omega _\rho \mu \omega _{\rho \rho }\rho +2\mu \omega _\rho \gamma _\rho \rho 3\mu _\rho \omega _\rho \rho )e^{2\gamma }],`$ (27)
$`\mathrm{\Theta }_0^1={\displaystyle \frac{1}{16\pi \mu ^2\rho ^3}}[\mu ^2(\dot{\mu }x+2\dot{\mu _\rho }\rho x+2\mu \omega \dot{\omega }\rho ^2x2\mu \omega \ddot{\gamma }\rho ^2y+2\mu \omega _\rho \gamma _\rho \rho ^2y`$ (28)
$`\mu \omega _{\rho \rho }\rho ^2y+2\mu \omega \gamma _{\rho \rho }\rho ^2y+\dot{\mu }\omega ^2\rho ^2x+2\mu _\rho \omega \gamma _\rho \rho ^2y3\mu _\rho \omega _\rho \rho ^2y`$ (29)
$`2\mu _{\rho \rho }\omega \rho ^2y)\rho ^2x(2\mu \dot{\gamma }\dot{\mu })e^{2\gamma }+\mu ^4\rho x(\mu \dot{\omega }\omega _\rho `$ (30)
$`+\mu \omega \dot{\omega _\rho }3\mu \omega \omega _\rho \dot{\gamma }+3\dot{\mu }\omega \omega _\rho )e^{2\gamma }],`$ (31)
$`\mathrm{\Theta }_0^2={\displaystyle \frac{1}{16\pi \mu ^2\rho ^3}}[\mu ^2(2\dot{\mu _\rho }\rho y+\mu \omega _{\rho \rho }\rho ^2y2\mu \omega _\rho \gamma _\rho \rho ^2y2\mu \omega \gamma _{\rho \rho }\rho ^2y+2\mu \omega \dot{\omega }\rho ^2y`$ (32)
$`+2\mu \dot{\omega }\dot{\gamma }\rho ^2x+2\mu \omega \ddot{\gamma }\rho ^2x+\dot{\mu }\omega ^2\rho ^2y+2\dot{\mu }\omega \dot{\gamma }\rho ^2x+3\mu _\rho \omega _\rho \rho ^2y`$ (33)
$`2\mu _\rho \omega \gamma _\rho \rho ^2y+2\mu _{\rho \rho }\omega \rho ^2y)\rho ^2y(2\mu \dot{\gamma }\dot{\mu })e^{2\gamma }+\mu ^4\rho y(\mu \dot{\omega }\omega _\rho `$ (34)
$`+\mu \omega \dot{\omega _\rho }2\mu \omega \omega _\rho \dot{\gamma }+3\dot{\mu }\omega \omega _\rho )e^{2\gamma }],`$ (35)
$$\mathrm{\Theta }_0^3=\mathrm{\Theta }_3^0=0.$$
(36)
Now for $`\omega =0`$ and $`\mu =\rho `$, Eqs.(21)-(26) become

$$\mathrm{\Theta }_0^0=\frac{1}{8\pi }e^{2\gamma }(\psi _\rho ^2+\psi _t^2),$$
(37)
$$\mathrm{\Theta }_1^0=\frac{1}{4\pi \rho }x\psi _\rho \psi _t,$$
(38)
$$\mathrm{\Theta }_2^0=\frac{1}{4\pi \rho }y\psi _\rho \psi _t,$$
(39)
$$\mathrm{\Theta }_0^1=e^{2\gamma }\mathrm{\Theta }_1^0,$$
(40)
$$\mathrm{\Theta }_0^2=e^{2\gamma }\mathrm{\Theta }_2^0,$$
(41)
$$\mathrm{\Theta }_0^3=\mathrm{\Theta }_3^0=0.$$
(42)
These are the energy and momentum densities of cylindrical gravitational waves given by Rosen and Virbhadra .

## IV Energy and Momentum in Papapetrou’s Prescription

The symmetric energy-momentum complex of Papapetrou  is given by

$$\mathrm{\Omega }^{ab}=\frac{1}{16\pi }N_{,cd}^{abcd},$$
(43)
where

$$N^{abcd}=\sqrt{g}(g^{ab}\eta ^{cd}g^{ac}\eta ^{bd}+g^{cd}\eta ^{ab}g^{bd}\eta ^{ac}),$$
(44)
and $`\eta ^{ab}`$ is the Minkowski spacetime. The energy-momentum complex satisfies the local conservation laws

$$\frac{\mathrm{\Omega }^{ab}}{x^b}=0.$$
(45)
The locally conserved energy-momentum complex $`\mathrm{\Omega }^{ab}`$ contains contributions from the matter, non-gravitational and gravitational fields. $`\mathrm{\Omega }^{00}`$ and $`\mathrm{\Omega }^{0a}`$ are the energy and momentum (energy current) density components. The required non-vanishing components of $`N^{abcd}`$ are given as

$$N^{0101}=\frac{1}{\mu \rho ^3}(\mu ^2\rho ^2+\mu ^2x^2\mu ^2\rho ^2\omega ^2y^2+\rho ^2y^2e^{2\gamma }),$$
(46)
$$N^{0102}=\frac{1}{\mu \rho ^3}(\mu ^2xy+\mu ^2\omega ^2\rho ^2xy\rho ^2xye^{2\gamma }),$$
(47)
$$N^{0202}=\frac{1}{\mu \rho ^3}(\mu ^2\rho ^2+\mu ^2y^2\mu ^2\rho ^2\omega ^2x^2\rho ^2x^2e^{2\gamma }),$$
(48)
$$N^{0303}=\frac{\mu }{\rho }(1+e^{2\gamma 4\psi }),$$
(49)
$$N^{0121}=\frac{1}{\rho }(\mu \omega x),$$
(50)
$$N^{0122}=\frac{1}{\rho }(\mu \omega y),$$
(51)
Substituting Eqs.(36-41) in Eq.(33), we have energy and momentum density components in Papapetrou’s prescription

$`\mathrm{\Omega }^{00}={\displaystyle \frac{1}{16\pi \mu ^2\rho ^3}}[\mu ^2(\mu +\mu _\rho \rho 2\mu _{\rho \rho }\rho ^2\mu \omega ^2\rho ^22\mu \omega \omega _\rho \rho ^3\mu _\rho \omega ^2\rho ^3)`$ (52)
$`+(\mu \rho ^2+2\mu \gamma _\rho \rho ^3\mu _\rho \rho ^3)e^{2\gamma }],`$ (53)
$`\mathrm{\Omega }^{01}={\displaystyle \frac{1}{16\pi \mu ^2\rho ^3}}[\mu ^2(\dot{\mu }x+2\dot{\mu _\rho }\rho x+\mu \omega y\mu \omega _\rho \rho y\mu \omega _{\rho \rho }\rho ^2y`$ (54)
$`+2\mu \omega \dot{\omega }\rho ^2x+\dot{\mu }\omega ^2\rho ^2x\mu _\rho \omega \rho y2\mu _\rho \omega _\rho \rho ^2y`$ (55)
$`\mu _{\rho \rho }\omega \rho ^2y)\rho ^2x(2\mu \dot{\gamma }\dot{\mu })e^{2\gamma }],`$ (56)
$`\mathrm{\Omega }^{02}={\displaystyle \frac{1}{16\pi \mu ^2\rho ^3}}[\mu ^2(\dot{\mu }y+2\dot{\mu _\rho }\rho y+2\mu \omega \dot{\omega }\rho ^2y+\dot{\mu }\omega ^2\rho ^2y\mu \omega x`$ (57)
$`+\mu \omega _\rho \rho x+\mu \omega _{\rho \rho }\rho ^2x+\mu _\rho \omega \rho x+2\mu _\rho \omega _\rho \rho ^2x`$ (58)
$`+\mu _{\rho \rho }\omega \rho ^2x)\rho ^2y(2\mu \dot{\gamma }\dot{\mu })e^{2\gamma }],`$ (59)
$$\mathrm{\Omega }^{03}=\mathrm{\Omega }^{30}=0.$$
(60)
We see that for $`\omega =0`$ and $`\mu =\rho `$, Eqs.(42)-(45) yield

$$\mathrm{\Omega }^{00}=\frac{1}{8\pi }e^{2\gamma }(\psi _\rho ^2+\psi _t^2),$$
(61)
$$\mathrm{\Omega }^{01}=\frac{1}{4\pi \rho }x\psi _\rho \psi _te^{2\gamma },$$
(62)
$$\mathrm{\Omega }^{02}=\frac{1}{4\pi \rho }y\psi _\rho \psi _te^{2\gamma },$$
(63)
$$\mathrm{\Omega }^{03}=0.$$
(64)
These turn out to be energy and momentum density components for cylindrical gravitational waves given by Virbhadra .

## V Discussion

We have evaluated energy and momentum density components for a class of rotating cylindrical gravitational waves by using prescriptions of Einstein and Papapetrou. It can be seen that the energy and momentum densities for a class of rotating gravitational waves are finite and well defined in both the prescriptions. It follows from Eqs.(21-26) and (42-45) that though the energy-momentum complexes of Einstein and Papapetrou are not exactly the same but are similar upto certain terms. However, it is interesting to note from Eqs.(27-32) and (46-49) that both the results reduce to the same energy and momentum densities of a special case of cylindrical gravitational waves as given in .

Acknowledgment

The author would like to thank Prof. Chul H. Lee for his hospitality at the Department of Physics and Korea Scientific and Engineering Foundation (KOSEF) for postdoc fellowship at Hanyang University Seoul, KOREA.

References

 Penrose, R.: Proc. Roy. Soc. London A381(1982)53;
Misner, C.W. and Sharp, D.H. Phys. Rev. 136(1964)B571.

 Bondi, H.: Proc. Roy. Soc. London A427(1990)249.

 Cooperstock, F.I.: Found. Phys. 22(1992)1011; in Topics in Quantum Gravity and Beyond: Papers in Honor of L. Witten eds. Mansouri, F. and Scanio, J.J. (World Scientific, Singapore, 1993) 201; in Relativistic Astrophysics and Cosmology, eds. Buitrago et al (World Scientific, Singapore, 1997)61; Annals Phys. 282(2000)115.

 Misner, C.W., Thorne, K.S. and Wheeler, J.A.: Gravitation (W.H. Freeman, San Francisco, 1973).

 Ehlers, J. and Kundt, W.: Gravitation: An Introduction to Current Research, ed. L. Witten (Wiley, New York, 1962)49.

 Pirani, F.A.E.: Gravitation: An Introduction to Current Research, ed. L. Witten (Wiley, New York, 1962)199.

 Weber, J. and Wheeler, J.A.: Rev. Mod. Phy. 29(1957)509;
Weber, J.: General Relativity and Gravitational Waves, (Interscience, New York, 1961).

 Qadir, A. and Sharif, M.: Physics Letters A 167(1992)331.

 Rosen, N: Helv. Phys. Acta. Suppl. 4(1956)171.

 Scheidegger, A.E.: Rev. Mod. Phys. 25(1953)451.

 Rosen, N.: Phys. Rev. 291(1958)110.

 Rosen, N. and Virbhadra, K.S.: Gen. Rel. and Grav. 26(1993)429.

 Virbhadra, K.S.: Pramana-J. Phys. 45(1995)215.

 Mashhoon, Bahram , McClune, James C. and Quevedo, Hernando: Class. Quantum Grav.17(2000)533.

 Kramer, D., Stephani, H. , MacCallum, M. and Herlt, E.: Exact Solutions of Einstein’s Field Equations (Cambridge University Press, Cambridge, 1980).

 Verdaguer, E.: Phys. Rep. 229(1993)1;
Griffths, J.B.: Colliding Plane Waves in General Relativity (Oxford University Press, Oxford, 1991);
Alekseev, G.A. and Griffths, J.B.: Class. Quantum Grav. 13(1996)2191.

 M$`\ddot{o}`$ller, C.: Ann. Phys. (NY) 4(1958)347; 12(1961)118.

 Papapetrou, A.: Proc. R. Irish. Acad. A52(1948)11.


format: text

text:
# Acceleration and collimation of relativistic plasmas ejected by fast rotators

## 1 Introduction

Radio pulsars (Michel michel91 (1991)), superluminal galactic objects (Mirabel & Rodriguez mirabel94 (1994), mirabel99 (1998)) and AGNs (Begelman, Blandford & Rees blandford (1984), Pelletier et al. pelletier (1996)) are observed to be associated with relativistic outflows. AGNs and superluminal galactic objects produce jets, while radio pulsars demonstrate a surprisingly high efficiency of transforming the pulsar rotational energy into kinetic energy of the relativistic wind (Kennel & Coroniti kennel (1984)). Hence, the problem of acceleration and collimation of relativistic plasmas is among the most important in the physics of relativistic outflows.

Nevertheless, radio pulsars, AGNs and superluminal galactic objects are objects of rather different nature and scales. It may well be that the mechanisms of acceleration of the outflows in these sources are not the same. However, there seems to exist a mechanism which operates very well for all these objects, namely, the so-called magnetocentrifugal acceleration where the acceleration of the plasma outflow occurs in the magnetic field of the rotating object. Such conditions are realized in the magnetospheres of radio pulsars, in the winds from accretion disks (Blandford & Payne bp (1982)) and apparently in superluminal galactic objects (Mirabel & Rodriguez mirabel99 (1998)). However, it has not been investigated whether the magnetocentrifugal mechanism could provide the acceleration of the relativistic plasma in these objects to the observed bulk motion velocities.

The analysis of problem of the relativistic plasma acceleration and collimation in a rotating magnetosphere started with Michel (michel (1969)), where it was shown that the magnetocentrifugal acceleration of the plasma is not effective for a prescribed monopole-like magnetic field. Subsequently, the relativistic plasma outflow has been investigated in several other studies and for different astrophysical conditions (Okamoto okamoto78 (1978), Ardavan ardavan (1979), Beskin et al. bgi (1983), Camenzind max86 (1986)). In particular, the ultrarelativistic plasma outflow from an accretion disk was investigated numerically in Camenzind (max87 (1987)) while the structure of electromagnetically driven relativistic jets was discussed in Li (li93 (1993)) under the assumption of jet confinement by an ambient medium. Collimation of stellar magnetospheres to jets in the force-free approximation was studied in Sulkanen & Lovelace (ml (1990)) and Fendt et al. (fendt95 (1995)). Later, Fendt & Camenzind (fendt96 (1996)) investigated the plasma dynamics in the MHD approximation in magnetic configurations obtained via the force-free approximation. A class of self-similar solutions of the relativistic plasma outflow has been obtained in Contopoulos (ianos94 (1994)). It was found that the efficiency of the plasma acceleration can be higher than in Michel’s solution due to the effect of the magnetic collimation of the wind to the axis of rotation (Begelman & Li (li (1994)), Fendt & Camenzind fendt96 (1996), Takahashi & Shibata (ts (1998))). Thus, the problem of the relativistic plasma acceleration in a rotating magnetic field appeared to be closely connected to the problem of plasma collimation.

Strong evidence has been accumulated for a correlation between collimated astrophysical outflows (jets) and accretion disks. Direct evidence for a disk-jet connection has been demonstrated by HST observations of the Herbig-Haro object HH30 (Burrows et al. burrows (1996)), although on the theoretical front this connection was discussed earlier (Camenzind max90 (1990), Ferreira & Pelletier ferreira (1993), Falke & Biermann falke (1995)). A disk-jet connection is also supported by observations of the active galaxy M87 (Ford et al. ford (1994)) and superluminal galactic objects (Mirabel & Rodriguez mirabel98 (1998)). These observations seem to lead to the conclusion that jets require the existence of accretion disks for their formation (e.g. Livio livio (1999)) and therefore accretion disks should be a necessary element of the collimation mechanism. However, we have by now reliable evidence that there is at least one source producing jets without the presence of an accretion disk. Namely, recent observations of the Crab Nebula by the Chandra observatory leave no doubts that the Crab pulsar indeed produces jets (Weiskopff et al. crab (2000)). Therefore, it would be important to understand the mechanism of plasma collimation in all known jet sources, including the Crab pulsar.

Such a mechanism is provided by magnetic collimation of outflows from accretion disks and it was first suggested by Lovelace (lovelace (1976)), Blandford (blandford76 (1976)) and Bisnovatyi-Kogan & Ruzmaikin (bkgr (1976)). Subsequently, it was shown from an analysis of the full set of the MHD equations for nonrelativistic (Heayvaerts & Norman heyvaerts (1989)) and relativistic outflows (Chiueh et al. chiueh (1991), Bogovalov bog95 (1995)) that collimation is a rather general property of magnetised winds from a rotating object. A magnetised fast-mode supersonic wind inevitably will be partially cylindrically collimated along the axis of rotation at distances large compared to all characteristic dimensions of the source, such as the size of the accretion disk, provided that: the flow is non-dissipative, the magnetised central object is rotating, the total magnetic flux of one polarity reaching infinity (open field lines) is finite, the polytropic index of the plasma exceeds unity and the pressure of ambient matter is low enough not to be able to terminate the wind before the formation of the jet. This result is independent of the specific properties of the central source. Collimation occurs spontaneously due to the internal Lorentz force arising in the magnetised plasma which is ejected from the rotating source and therefore can be regarded as a process of self-collimation (hereafter SC process). Note that the assumption of a polytropic gas with a constant polytropic index $`\gamma `$ does not seem really restrictive, since several exact MHD solutions for collimated outflows have been constructed for the case of a variable index $`\gamma `$, Sauty & Tsinganos (1994), Vlahakis & Tsinganos (1998).

Although the SC process has been questioned (e.g., Okamoto okamoto99 (1999)), a numerical investigation of this process and the characteristics of jets produced via this mechanism in the cold (Bogovalov & Tsinganos bogts (1999)) or hot (Tsinganos & Bogovalov tsbog (2000), Sauty et al (1999)) nonrelativistic winds leave no doubt that this mechanism works. This mechanism provides the formation of cylindrically collimated outflows at a large distance from the central source while collimation is accompanied by essentially more efficient acceleration of the plasma at the equatorial region as compared with the original study of Michel (michel (1969)). Evidence for the focusing of MHD outflows from accretion disks to the axis of rotation at distances of the order of the dimension of accretion disks has been obtained in other studies as well, such as in Ouyed & Pudritz (op (1997)), Ustyugova et.al. (ustyug (1999)) and Krasnopolskii et. al. (krasno (1999)).

In the analysis of Bogovalov & Tsinganos (1999), the relaxation method was used to obtain the solution of the steady state outflow problem, i.e., the time-dependent problem was solved in order to get a steady-state solution. This method was first employed by Washimi & Shibata (shibata (1993)) for outflows from stellar objects. Later this method was applied in Bogovalov (1996) for the investigation of the SC process in an initially monopole-like magnetic field. It seems that this numerical method allows us to solve accurately all problems connected with the critical surfaces (Bogovalov 1997b , Tsinganos et al. tsincr (1996)).

In our previous study (Bogovalov 1997a ), the relaxation method has been applied to the investigation of the relativistic plasma outflow and it was found that the relativistic plasma is very weakly collimated and accelerated even if the Poynting flux dominates the matter energy flux at the equator. However, an analytical study performed in this work has shown that an efficient SC of the relativistic plasma is still possible for very fast rotation of the central source. Recently it has been argued that a relativistic plasma is poorly accelerated even under this condition (Beskin et al. beskin98 (1998)). However, the SC process was totally neglected in this work. Therefore, the question concerning the efficiency of the acceleration and SC of the relativistic winds from fast rotators remains open.

## 2 Basic assumptions

Outflows of magnetised plasmas are described by the system of the familiar nonlinear MHD equations. It is not surprising that this system has many solutions corresponding to different physical conditions. Then, a reasonable choice of the appropriate approximation and a sound model is crucially important to arrive at some physically meaningful solution.

To formulate the basic demands of the model let us consider the general properties of winds and jets. According to observations, the length of the jets greatly exceeds all characteristic scales of the central source. Therefore, for an investigation of the SC process and the characteristics of the jets, the theory should allow us to obtain the solution on scales greatly exceeding all sizes of the central source.

Observations of winds and jets from several classes of astrophysical objects such as our Sun (Parker parker (1963)), or, YSO’s (Burrows et al. burrows (1996)), pulsars (Kennel & Coroniti kennel (1984)) and AGNs (Ferrari et al. ferrari (1996)) show that they are super-magnetosonic at relatively large distances from the source. At the same time all winds must be subsonic near the central source where the basic parameters of the outflows are controlled. Therefore, the theory must describe flows of mixed type, containing the small scale regions of a subsonic flow and the large scale regions of the super magnetosonic flow. The scales of these regions may differ by several orders of magnitude. These properties of the outflows dictate the choice of the approximations and mathematical models.

### 2.1 Approximations

The force-free approximation has been widely used for the analysis of axisymmetric relativistic outflows (Scharlemann & Wagoner sw (1973), Mestel & Wang mw (1979), Beskin et al. bgi (1983), Lyubarskii lyubarskii (1990), Fendt et al. fendt95 (1995)). For simplification reasons, the axisymmetric stellar magnetosphere can be assumed initially to have a dipolar magnetic field (Contopoulos et al.kazanas (1999)). However, this approximation is not appropriate for the investigation of the SC process for several reasons. For example, it describes only subsonic flows of plasma since the force-free transfield (Grad-Shafranov) equation is elliptic everywhere (Fendt et al. fendt95 (1995)). Therefore, this approximation cannot be used for the investigation of plasma outflow in the supersonic regions. Moreover, the collimating and accelerating forces are equal to zero in this approximation. For these reasons, the force-free approximation gives ambiguous results. For example, the structure of the flow of the relativistic plasma at large distances from the source obtained in Sulkanen & Lovelace ml (1990) and in Fendt et al. fendt95 (1995) differs from the structure obtained in Contopoulos et al. (kazanas (1999)). Although Sulkanen & Lovelace (ml (1990)) and Contopoulos et al. (kazanas (1999)) used similar models, Contopoulos et al. (kazanas (1999)) obtained a radially outflowing plasma while Sulkanen & Lovelace (ml (1990)) obtained collimated jets. Jets were obtained as well in Fendt et al. (fendt95 (1995)). Actually all these works are mathematically correct. The difference between them is only in the boundary conditions at infinity or at the outer boundary of the integration domain. These works clearly demonstrate that the flow of plasma at large distances from the source in the force-free approximation basically depends on the boundary conditions at infinity. But there is no such dependence in real supersonic flows. This is the most important qualitative difference between the subsonic and super fast magnetosonic flows.

In the steady state subsonic flow the MHD signals propagate down and upstream of the flow. That is why the plasma flow is defined by the forces affecting the plasma at a given point, upstream and down stream of the flow. This is not the case in the super fast magnetosonic flow. No MHD signals can propagate upstream of the supersonic flow. The supersonic flow is defined only by the forces which affect the plasma upstream and at the given point and does not depend on the forces which affect the plasma downstream of the flow. For this reason the supersonic flow does not depend on the boundary conditions located downstream of the flow (actually the situation is more complicated; the flow should cross the fast mode critical surface, Bogovalov 1997b ). The simplest approximation which allows us to describe the mixed type outflows is the ideal MHD approximation. This approximation is used in this paper.

### 2.2 Model

However, even the use of the MHD approximation does not always allow us to investigate the SC process in the supersonic region. There is a wide class of self-similar solutions for the nonrelativistic (Blandford & Payne bp (1982), Vlahakis & Tsinganos tv (1999) e.g.) and relativistic plasma (Contopoulos ianos94 (1994)). All these solutions terminate not so far from the source. Independent of why this termination happens, the self-similar solutions cannot be used for the verification of the predictions of the SC theory since the central source has infinitely large size with infinite magnetic flux in the upper hemisphere. At any distance from the source, the dimension of the source exceeds the distance to the source. However, we are interested in the opposite case, when the distance is large compared with the dimension of the central source. Therefore, we demand that the appropriate model should have finite geometrical dimensions, finite magnetic flux of the open field lines of every polarity and initially uncollimated plasma outflow. The last condition is important for us since our objective is to study the process which transforms the initially uncollimated outflow into the outflow containing the jet. The word “initially” means that if one of the conditions for the SC process is not fulfilled (for example if the central source is not rotating or the magnetic field is absent) the wind expands radially at large distances from the source. Actually, a lot of models satisfy all these conditions.

The SC is a general property of magnetised winds from rotating sources. It must take place in a wind from any source provided that the conditions specified by Heayvaerts & Norman (heyvaerts (1989)), Chiueh et al. (chiueh (1991)) and Bogovalov (bog95 (1995)) are fulfilled. Therefore “convenience” and “simplicity” become reasonable arguments for the choice of the model for the investigation of this process. The model of a star with a wind in an initially monopole-like (below m-l) magnetic field firstly proposed by Sakurai (sakurai (1985)) is the most simple and convenient. This model allows us to avoid the mathematically complicated problem connected with the coexistence of open and closed field lines near the source which actually has no direct relation to the SC process and usually is considered particular (Keppens & Goedbloed keppens (1999)).

The model with the initially monopole-like magnetic field describes more a wider class of outflows, including the time-dependent outflows from oblique rotators (Bogovalov bog99 (1999)) and outflows from objects with arbitrary distribution of the polarity of the magnetic field at the base (Tsinganos & Bogovalov tsbog (2000)).

In this paper we consider the outflow of an initially cold relativistic plasma and neglect only gravitation since we are interested in the effect of the electromagnetic forces only on the plasma.

## 3 Basic equations and methods of solution.

The plasma outflow from astrophysical objects consists of two zones where the properties of the flow differ qualitatively. The nearest zone is located near the surface of the central source and includes all the critical surfaces of the steady state flow (Bogovalov 1997b ). The outer boundary of the nearest zone can be located rather arbitrary but always downstream of the fast mode critical surface. The flow of the plasma in the nearest zone is of the mixed type. It is subsonic near the source and is super fast magnetosonic at the outer boundary. The far zone is located downstream of the nearest zone. The plasma flow in the far zone is super fast magnetosonic. Different methods should be applied for the solution of the problem in each of these zones.

The solution of the problem of the steady state plasma outflow can be divided in two steps (Bogovalov & Tsinganos bogts (1999)). In the first step, the problem is solved in the nearest zone by the relaxation method and in the second, the steady state solution in the far zone is obtained by solving directly the transfield equation as a Cauchy problem. The solution in the nearest zone is used to specify the initial values for the flow in the far zone.

### 3.1 Solution in the nearest zone.

The system of equations defining the dynamics of the relativistic plasma in the ideal MHD approximation is as follows,

$$\rho c(\dot{𝐔}+(𝐯)𝐔)=q𝐄+\frac{1}{c}𝐣\times 𝐁,$$
(1)
$$𝐄+\frac{1}{c}𝐯\times 𝐁=0,$$
(2)
$$\times 𝐁=\frac{4\pi }{c}𝐣+\frac{1}{c}\dot{𝐄},$$
(3)
$$\frac{1}{c}\dot{𝐁}=\times 𝐄,$$
(4)
$$\dot{\rho }+(𝐯\rho )=0.$$
(5)
Here $`q=(𝐄)/4\pi `$ is the induced electric charge density, $`𝐣`$ is the electric current density, $`𝐄`$ is the electric field, $`𝐔=𝐯\gamma /c`$ is the spatial component of the plasma four-velocity and $`\rho `$ is the mass density of the plasma.

The problem was solved in the cylindrical system of coordinates. Some details of the solution method are given in appendix A. The numerical simulation was performed only in the quarter of the total box of simulation due to the symmetry of the flow in relation to the equatorial plane and the axis of rotation. The boundaries of the nearest zone include the surface of the star, the axis of rotation, the equator and the outer boundaries, which are a cylinder (right hand side of the box) and a plane perpendicular the axis of rotation at some distance from the source (upper side of the box). The outer boundaries are located in the supersonic domain of the flow beyond the fast mode critical surface. This guarantees the independence of the solution in the nearest zone of the position and conditions of the outer boundaries.

The number of independent boundary conditions on the surface of the star is equal to the number of MHD waves leaving the surface to satisfy the causality in the flow (Bogovalov 1997b ). They are as follows:
1. A constant plasma density $`\rho _0`$.
2. A constant Lorentz-factor $`\gamma _0`$ of the plasma in the corotating frame system.
3. A constant normal component of the poloidal magnetic field $`B_0`$.
4. The continuity of the tangential component of the poloidal electric field.
5. A constant temperature of the plasma equal to zero.

The boundary conditions at the axis of rotation correspond to the axial symmetry of the flow. The boundary conditions at the equator follow from the symmetry of the flow in relation to the equatorial plane. This means that $`\rho (z)=\rho (z)`$, $`v_r(z)=v_r(z)`$, $`v_z(z)=v_z(z)`$, $`B_r(z)=B_r(z)`$, $`B_z(z)=B_z(z)`$, $`v_\phi (z)=v_\phi (z)`$ and $`B_\phi (z)=B_\phi (z)`$.

At the outer boundaries only one side internal derivatives were used for calculation of the derivatives in the finite differences, since the solution does not depend on the conditions down the flow from these boundaries.

### 3.2 Solution in the far zone

The flow in the far zone is fast mode supersonic. The system of equations describing the steady state flow is hyperbolic. Therefore, a Cauchy problem can be formulated for this flow.

It is convenient to introduce the flux function $`\psi `$ as follows

$$𝐁_p=\frac{\psi \times \widehat{\phi }}{r},$$
(6)
Where $`B_p`$ is the poloidal magnetic field of the axisymmetric outflow and $`\widehat{\phi }`$ is the unit azimuthal vector.

The stationary MHD equations admit well known four integrals (Tsinganos tsin82 (1982)). They are:

The ratio of the poloidal magnetic and mass fluxes, $`F(\psi )`$

$$F(\psi )=\frac{B_p}{4\pi \rho v_pc}.$$
(7)
The total angular momentum per unit mass $`L(\psi )`$,

$$rU_\phi FrB_\phi =L(\psi ).$$
(8)
The corotation frequency $`\mathrm{\Omega }(\psi )`$ in the frozen-in condition

$$U_\phi B_pU_pB_\phi =r\gamma B_p\mathrm{\Omega }(\psi )/c.$$
(9)
The total energy $`W(\psi )`$ per particle in the equation for total energy conservation,

$$\gamma F(\psi )r\mathrm{\Omega }(\psi )B_\phi /c=W(\psi ).$$
(10)
This system of equations should be supplemented by a relativistic relationship between the components of the four-velocity

$$\gamma ^2=1+U_p^2+U_\phi ^2.$$
(11)
Momentum balance across the poloidal field lines is expressed by the transfield equation. For analysing the behavior of the plasma at large distances, it is convenient to deal with this equation in an orthogonal curvilinear coordinate system ($`\psi ,\eta `$) formed by the poloidal magnetic field lines and by the lines of the poloidal electric field. $`\psi `$ varies with the motion across the magnetic field lines, while $`\eta `$ varies with the motion along the magnetic field lines. A geometrical interval in these coordinates can be expressed as

$$(d𝐫)^2=g_\psi ^2d\psi ^2+g_\eta ^2d\eta ^2+r^2d\phi ^2,$$
(12)
where $`g_\psi ,g_\eta `$ are the corresponding line elements (components of the metric tensor).

According to Landau & Lifshitz (1975) the equation $`T_{\psi ;k}^k=0`$, where $`T^{ij}`$ is the energy-momentum tensor and ”; k” denotes covariant differentiation, will have the following form in these coordinates,

$$\begin{array}{c}\frac{}{\psi }\left[\frac{B^2E^2}{8\pi }\right]\frac{1}{r}\frac{r}{\psi }\left[U_\phi v_\phi c\rho \frac{B_\phi ^2E^2}{4\pi }\right]\\ \\ \frac{1}{g_\eta }\frac{g_\eta }{\psi }\left[U_pv_pc\rho \frac{B_p^2E^2}{4\pi }\right]=0\end{array}$$
(13)
with the electric field $`E=\frac{r\mathrm{\Omega }}{c}B_p`$.

The unknown variables here are $`z(\eta ,\psi )`$ and $`r(\eta ,\psi )`$. The metric coefficient $`g_\eta `$ can be obtained from the transfield equation (13),

$$g_\eta =\mathrm{exp}(\underset{0}{\overset{\psi }{}}G(\eta ,\psi )𝑑\psi ),$$
(14)
where

$$G(\eta ,\psi )=\frac{\frac{}{\psi }\left[\frac{B^2E^2}{8\pi }\right]\frac{1}{r}\frac{r}{\psi }\left[U_\phi v_\phi c\rho \frac{B_\phi ^2E^2}{4\pi }\right]}{\left[U_pv_pc\rho \frac{B_p^2E^2}{4\pi }\right]}.$$
(15)
The lower limit of the integration in (14) is chosen to be 0 such that the coordinate $`\eta `$ is uniquely defined. In this way $`\eta `$ coincides with the coordinate $`z`$ where the surface of constant $`\eta `$ crosses the axis of rotation.

The metric coefficient $`g_\psi `$ can be obtained from definition (6) in terms of the magnitude of the poloidal magnetic field as follows

$$g_\psi =\frac{1}{rB_p}.$$
(16)
The orthogonality condition

$$r_\eta r_\psi +z_\eta z_\psi =0,$$
(17)
and the relationships

$$g_\eta ^2=r_\eta ^2+z_\eta ^2,g_\psi ^2=r_\psi ^2+z_\psi ^2$$
(18)
give us that

$$r_\eta =\frac{z_\psi g_\eta }{g_\psi },z_\eta =\frac{r_\psi g_\eta }{g_\psi },$$
(19)
with $`g_\eta `$ calculated by the expression (14). Here $`r_\eta =r/\eta `$, $`z_\eta =z/\eta `$, $`r_\psi =r/\psi `$, $`z_\psi =z/\psi `$. For the numerical solution of the system of equations (19) the two step Lax-Wendroff method on the lattice with a dimension equal to 1000 is used.

Equations (19) should be supplemented by appropriate boundary conditions and initial values on some initial surface of constant $`\eta `$ located in the nearest zone, but down the flow from all the critical surfaces. The form of the initial surface of constant $`\eta `$ was obtained numerically through integration of the following equations,

$$\frac{r}{\psi }=\frac{B_z}{rB_p^2},\frac{z}{\psi }=\frac{B_r}{rB_p^2}.$$
(20)
The values $`B_z`$, $`B_r`$ and integrals $`W(\psi ),L(\psi ),\mathrm{\Omega }(\psi ),F(\psi )`$ were taken from the solution of the problem in the nearest zone.

Boundary conditions on the axis of rotation and the equatorial plane are the same as the conditions in the nearest zone. No conditions at infinity were specified.

## 4 Parameterisation of the solution

A stationary solution of the relativistic plasma flow from the star with an initially m-l magnetic field is defined by three dimensionless parameters (Bogovalov bog99 (1999)). One of them is the ratio of the radius of the star $`R_{star}`$ to the radius of the initial fast mode magnetosonic surface $`R_f=\sqrt{(B_0R_{star})^2/4\pi \rho _0v_0cU_0}`$, where $`U_0=v_0\gamma _0/c`$. The word ”initial” here refers to the radius of the fast mode surface when the star is not rotating. But this dependence can be neglected under the condition $`R_{star}R_f`$ and $`R_{star}R_l`$, where $`R_l=c/\mathrm{\Omega }`$ is the radius of the light cylinder, and $`\mathrm{\Omega }`$ is the angular velocity of the uniformly rotating star. Therefore, the solution crucially depends only on two dimensionless parameters which can be chosen rather arbitrary. It is convenient to choose one of them as follows

$$\sigma =(\frac{R_f}{R_l})^2=\frac{B_0^2}{4\pi \rho _0\gamma _0v_0^2}(\frac{R_{star}\mathrm{\Omega }}{c})^2.$$
(21)
$`\sigma `$ is the squared ratio of the initial radius of the fast mode surface to the light cylinder radius. According to (21) the radius of the light cylinder in units of $`R_f`$ is

$$\frac{R_l}{R_f}=\sqrt{1/\sigma }.$$
(22)
The parameter $`\sigma `$ is closely related to the so-called magnetisation parameter usually used in the physics of the relativistic winds from radio pulsars (Arons arons (1996)). The magnetisation parameter is the ratio of the Poynting flux density to the flux density of the mass energy (including rest mass energy). The parameter $`\sigma `$ used here is equal to the ratio of the Poynting flux to the mass energy flux at the equator of the rotator with the initially m-l magnetic field provided that $`\sigma /U_0^2<<1`$ (Bogovalov bog99 (1999)). In our previous work (Bogovalov 1997a ), we used another definition of the $`\sigma `$ parameter introduced initially by Michel (michel (1969)). This parameter is denoted by $`\sigma _M`$. There is a simple relationship between them $`\sigma _M=\sigma U_0`$. It follows from this relationship that $`\sigma _M`$ is equal to the Poynting flux density per particle on the equator provided that $`\sigma /U_0^2<<1`$.

It is convenient to take the second parameter as follows

$$\alpha =\frac{R_f}{R_lU_0}=\frac{\sqrt{\sigma }}{U_0}.$$
(23)
The physical sense of this parameter becomes clear if we pay attention to the fact that in the model of Michel (michel (1969)) the terminating four-velocity of the plasma on the equator is $`U_{\mathrm{}}=(\sigma _M^{1/3})`$ and $`\alpha `$ can be presented as follows

$$\alpha =(U_{\mathrm{}}/U_0)^{3/2}.$$
(24)
According to Eq. (24), the parameter $`\alpha `$ shows how effectively the magnetic forces arising in the rotating magnetic field accelerate the plasma. If $`\alpha <1`$ acceleration is not important compared to the initial velocity of the plasma. It is natural to call the rotators with $`\alpha <1`$ as slow rotators. Correspondingly, the rotators with $`\alpha 1`$ should be fast rotators which could in principle (but not necessarily) accelerate the plasma. The parameter $`\alpha `$ introduced here is the relativistic generalization of the $`\alpha `$ parameter introduced earlier by Bogovalov & Tsinganos (bogts (1999)) for the nonrelativistic plasma. Numerical simulations of the nonrelativistic plasma outflow show that this parameter indeed divides rotators into two groups with different efficiency of acceleration and collimation of the plasma. Below we will see that this division remains valid also in the relativistic limit.

Representation of the solution in the parameters $`\alpha `$, $`\sigma `$ is convenient since the steady state outflow in the initially m-l magnetic field has a simple analytical solution in the limit $`\alpha 0`$ independent of the value of $`\sigma `$ (Bogovalov 1997a ) and the transition to the nonrelativistic limit corresponds simply to the limit $`\sigma 0`$ and $`v_0/c1`$.

For comparison with observation, Eq. (21) should be expressed in variables general for the m-l model and real astrophysical objects. An astrophysical object can be considered as equivalent to the object with the m-l magnetic field if they both produce winds with the same magnetic flux, the same mass flux, the same initial Lorentz factor and finally both rotate with the same angular velocity. The magnetisation parameter in these variables takes the form

$$\sigma =\frac{\psi _t^2}{(2\pi )^2m\dot{N}\gamma _0v_0}(\frac{\mathrm{\Omega }}{c})^2,$$
(25)
where $`\psi _t`$ is the open field lines flux of one polarity in the wind, $`\dot{N}`$ is the rate of ejection of particles with mass $`m`$. For radio pulsars, $`\psi _t=\pi B_0R_p^2`$ where $`R_p=R_{star}\sqrt{R_{star}\mathrm{\Omega }/c}`$ is the radius of the polar cap on the surface of the pulsar ( Ruderman & Sutherland rs (1975)).

The position of the radio pulsars in the parameter space $`\alpha `$, $`\sigma `$ was discussed in detail by Bogovalov (bog99 (1999)). All radio pulsars eject Poynting-dominated winds ($`\sigma >1`$) but basically all of them are slow rotators ($`\alpha <1`$). One of the youngest and fastest rotating pulsars with a strong magnetic field on the surface is the pulsar in the Crab Nebula. This pulsar rotates with $`\mathrm{\Omega }200`$ and ejects $`e^\pm `$ pairs with the rate $`\dot{N}410^{38}`$ per second and average $`\gamma _0200`$ ( Daugherty & Harding dh (1982)). At $`B_0410^{12}`$ G the magnetisation parameter of the wind $`\sigma 1.310^4`$ and $`\alpha 0.57`$ The light cylinder radius of the Crab pulsar is located at the distance $`R_L=1.510^8`$ cm from the axis of rotation. The radius of the fast magnetosonic surface $`R_f`$ exceeds the light cylinder radius by a factor of 114. Marginally, the Crab pulsar could be the fast rotator provided that $`\gamma _0<200`$. The uncertainties in estimates of this parameter do not rule out this possibility. But all other, more slowly rotating pulsars are certainly the slow rotators.

In AGNs and superluminal sources $`\gamma _0`$ is defined from observations of jets and lies in the range $`110`$. If the collimation of the plasma is really due to the SC process we should expect that $`\alpha >1`$ and $`\sigma >\gamma _0^2`$. So, the winds from AGNs and superluminal sources should be mainly Poynting flux dominated.

## 5 Basic results

In Fig. 1 we show the result of the numerical simulation of the plasma outflow in the nearest zone of the star, for a time interval $`t=2.62R_f/v_0`$ and for the set of the values of the two key parameters, $`\alpha =1.0`$ and $`\sigma =51`$. This time interval is larger by a factor of 2.62 compared to the typical time the plasma and all MHD perturbations take to travel from the surface of the star to the fast mode surface. Note that during this time interval, a steady state solution is reached, <sup>1</sup><sup>1</sup>1This solution is stable. In the previous work (Bogovalov 1997a ) the time was expressed in the units $`1/\mathrm{\Omega }`$. Therefore, the ratio of the simulation time in the present work over the time in the previous work is $`(t_{present}/t_{previous})(R_f\mathrm{\Omega }/v_0)(t_{present}/t_{previous})\sqrt{\sigma }`$ provided that $`v_0c`$. Here $`t_{present}`$ and $`t_{previous}`$ are the dimensionless times in the present and previous works. In the present work the solution is stable up to the time $`t_{present}=2.62`$ at $`\alpha =1`$ and $`\sigma =51`$. In the previous work the solution for $`\alpha =0.0014`$ and $`\sigma =3.3`$ was destroyed after $`t_{previous}=2.02`$. Thus, in the present work much stronger Poynting-dominated wind was simulated, for 9.26 times longer than in the previous one. The simulation shown in the Fig. 4 lasted even longer. This confirms the stability of the solution. since all perturbations which could be excited at the start of the simulation have left the box of simulation. The initial fast mode surface is circular with a radius equal to 1 in our units. Note that the final fast mode surface shown in Fig. 1 has a small amplitude kink structure at a polar angle around 45 degrees. Such a structure is formed in all simulations performed with the present method and in the present geometry (Washimi & Shibata shibata (1993), Bogovalov & Tsinganos bogts (1999)). Basically it is due to the shape of the surface of the star which is actually approximated by a sequence of step functions in the cylindrical system of coordinates.

The SC theory predicts that a wind of relativistic plasma consists of a cylindrically collimated core surrounded by a radially expanding flow (Heyvaerts & Norman heyvaerts (1989), Chiueh et al. chiueh (1991) and Bogovalov bog95 (1995)). An inspection of Fig. 1 shows that the collimation of the plasma along the axis of rotation is practically invisible in the nearest zone. The results of the same simulation in the far zone are presented in Fig. 2. At a first glance, in the far zone there is no cylindrical collimation. However, a closer examination of the results shows that a cylindrically collimated core is actually formed in the wind.

Fig. 3 gives the distribution of the poloidal magnetic field across the axis of rotation for two values of $`z`$. This distribution differs from that expected for a m-l wind. The m-l poloidal magnetic field is described by the equation $`B_p=B_0/(r^2+z^2)`$. And, at a fixed $`z`$ the poloidal field normalized to its value at the axis of rotation ($`B_{axis}`$) depends on $`r`$ as $`B_p/B_{axis}=1/(1+(r/z)^2)`$. Then, at $`z=1121R_f`$ the poloidal magnetic field remains practically constant at $`r5z=1121`$. This dependence is represented by the horizontal line $`B_p/B_{axis}=1`$ in Fig. 3. However, our numerical simulation gives a totally different picture. Namely, the poloidal magnetic field has clearly a maximum at the rotation axis and then drops down with the distance $`r`$ from the axis, for example at $`rR_f`$ to half its value from the axis. It is important that this distribution is practically independent of $`z`$. This result implies that the wind near the axis of rotation has achieved an asymptotic state and the flow is collimated exactly along the axis of rotation here.

The m-l model allows us to perform a quantitative comparison of the SC theory and our numerical experiment. The structure of the cylindrically collimated flow at large distances from the source is especially simple provided that the central source rotates uniformly while the plasma velocity and the ratio of the poloidal magnetic and mass fluxes (function $`cF`$ in Eq. (7)) are constant across the jet. In this case the dependence of $`B_p`$ on $`r`$ at the fixed $`z`$ for a cold plasma is given by the simple analytical expression (Bogovalov bog95 (1995))

$$\frac{B_p}{B_0}=\frac{1}{1+(r/R_j)^2},$$
(26)
where $`R_j=\gamma _jv_j/\mathrm{\Omega }=R_fU_j/U_0\alpha `$ is the characteristic radius of the collimated flow of the cold plasma, while $`\gamma _j`$ and $`v_j`$ are the Lorentz factor and the velocity of the plasma in this flow, respectively. Eq. (26) is a direct consequence of the SC theory and is valid at $`rR_l`$. The dependence of the normalised plasma density on $`r`$ is identical to the distribution (26) since $`F`$ and the plasma velocity are assumed to be constant across the cylindrically collimated flow. In this case the plasma density is simply proportional to the poloidal magnetic field, since $`\rho =B_p/4\pi v_pcF(\psi )`$ (see Eq. (7)).

All these conditions are indeed fulfilled in our case. The outflow speed at the vicinity of the axis of rotation is close to the speed of light ( $`v_p=0.87c`$ at $`\gamma _0=2`$ and $`v_p=0.98c`$ at $`\gamma _0=7`$) and is independent of $`r`$. The function $`F`$ is assumed independent of $`\psi `$ for all field lines and the star rotates uniformly. A comparison of the calculated dependence of $`B_p`$ on $`r`$ with the analytical prediction, Eq. (26), shows that the SC theory correctly predicts the asymptotic properties of the wind structure.

The cylindrically collimated flow in Fig. (2) is not visible since its transversal scale is extremely small. The flux of the poloidal magnetic field contained in this flow is small as well. The dependence of the magnetic flux on $`r`$ is shown in Fig. (3) by the dashed-dotted line. For plotting convenience, this flux is multiplied by 100 and is also normalised to the total magnetic flux in the upper hemisphere, $`\psi _{eq}`$, which coincides with the magnetic flux at the equator. The mass flux normalised to the total mass flux in the upper hemisphere is presented by the same curve since the mass flux is proportional to the poloidal magnetic field in our case. The fraction of the magnetic and mass flux in the cylindrically collimated flow is of the order $`0.1\%`$. The rest of the magnetic and mass flux expands radially, in accordance with the SC theory.

The same physical picture emerges for a larger value of $`\alpha `$. In Fig. (4) we show in the nearest zone the results of the simulation of the flow of plasma for $`\alpha =4`$.

The increase $`\alpha `$ does not change significantly the structure of the flow, although collimation now becomes visible even in the nearest zone. The same flow in the far zone is shown in Fig. (5). It may be seen from this figure that a very narrow collimated beam is formed in the wind. The structure of the cylindrically collimated flow is shown in Fig. (6) for two values of $`z`$.

The efficiency of the acceleration of the plasma by a rotating object is important in astrophysical applications. This acceleration in the model under consideration is expected to be most efficient near the equatorial plane since at the axis of rotation there is no accelerating force at all. The accelerating force in the poloidal direction is proportional to the product $`j_{}B_\phi `$, where $`j_{}`$ is the component of the electric current perpendicular to a field line. The accelerating force in the azimuthal direction is proportional to the product $`j_{}B_p`$. Thus, all accelerating forces are proportional to $`j_{}`$. This implies that the acceleration occurs only if the electric currents cross the poloidal field lines. The electric currents shown by the dashed line in Figs. (1) and (4) indeed cross the poloidal field lines. This is very well seen especially in Fig. (4) near the equator.

The dependence of the four-velocity of the plasma on the distance from the centre of the source is shown in Fig. (7). The distance is given in a logarithmic scale. The energy of the particles increases slower than logarithmically and reaches a constant value. It follows from energy conservation (10) that the maximum possible terminal four-velocity $`U_{\mathrm{}}`$ at the equator is $`U_{\mathrm{}}=U_0(v_{\mathrm{}}/v_0)\sigma `$, since $`W`$ at the equator can be estimated as $`\gamma _0\sigma `$, provided that $`\sigma 1`$ and the toroidal magnetic field is estimated as $`B_\phi r\mathrm{\Omega }B_p/v_0`$ (Bogovalov 1997a ). All Poynting flux should be transformed into kinetic energy of the plasma in this case. But this value of the velocity is never achieved. Michel (michel (1969)) calculated for $`U_{\mathrm{}}`$ the value $`U_{michel}=(\sigma U_0)^{1/3}`$ for the plasma flow in a prescribed and fixed m-l poloidal magnetic field. Therefore, $`U_{michel}/U_0=\alpha ^{2/3}`$ (see Eq. 24). It is reasonable to estimate the efficiency of the plasma acceleration as the ratio $`\xi =U_{\mathrm{}}/U_0\sigma `$ which is approximately the ratio of the part of the Poynting flux transformed into the kinetic energy of the plasma (this is valid, provided that $`U_{\mathrm{}}U_0`$) to the initial value of the Poynting flux at the equator in the relativistic limit. For the Michel solution (michel (1969)) $`\xi =\alpha ^{2/3}/\sigma `$. The same value is given by Beskin et al. (beskin98 (1998)) for the analytical solution of the same problem under the assumption that the poloidal magnetic field goes to a monopole-like field at $`\alpha \mathrm{}`$. It follows from Fig. (7) that $`U_{\mathrm{}}/U_0`$ increases with $`\alpha `$ essentially faster than $`\alpha ^{2/3}`$. For the flow with $`\alpha =4`$ Michel’s estimate gives $`U_{p,\mathrm{}}/U_0=2.5`$, while the numerical solution gives more than 8. Correspondingly, the efficiency of the acceleration remarkably exceeds the efficiency of the acceleration given by Michel and achieves $`16.7\%`$ for the flow at $`\alpha =4`$ and $`\sigma =48`$. It follows from Fig. 7 that the remarkable part of the acceleration occurs at large distances from the star. It happens at the collimation of the magnetic flux to the axis of rotation due to disbalance arising between the tension and pressure gradient of the toroidal magnetic field (Begelman, & Li li (1994), Fendt & Camenzind fendt96 (1996)). This mechanism was not taken into account in the analysis of Michel (michel (1969)) and Beskin et al. (beskin98 (1998)) since they assumed a fixed m-l magnetic field.

It follows from Fig. 7 that the largest energy is achieved not at the equator but at $`r=1000R_f`$ for $`\alpha =4`$. To understand why this happens it is convenient to look on the dependence of $`U_p/U_0`$ on $`\psi `$ across field lines. This dependence is presented in Fig. (8) by two upper curves 1 ( for $`z=1000R_f`$) and 2 (for $`z=10350R_f`$). These curves have small amplitude irregular structure which is purely of numerical origin. It arises at large distances where the program deals with large numbers. The maximum of the acceleration is achieved not at the equator, but at $`\psi /\psi _{eq}0.3`$. This is explained by the fact that the asymptotical regime predicted by the SC theory is actually not achieved at these $`z`$ for $`\psi /\psi _{eq}>0.05`$. It appears that the larger $`\psi `$ the larger the distance which is necessary to reach the asymptotic regime of the outflow. This property of the SC process has already been found in the nonrelativistic outflow simulation (Bogovalov & Tsinganos bogts (1999)).

The physics of this phenomena can be understood as follows. It can be shown from simple geometrical consideration that

$$\frac{1}{g_\eta }\frac{g_\eta }{\psi }=\frac{1}{rR_cB_p},$$
(27)
where $`R_c`$ is the curvature radius of the poloidal magnetic field lines, with $`R_c`$ positive if the centre of curvature is in the domain between the line and the axis of rotation and negative in the opposite case. With this expression, Eq. (13) takes the form

$$\begin{array}{c}\frac{}{\psi }\left[\frac{B_p^2}{8\pi }\right]+\frac{1}{8\pi r^2}\frac{}{\psi }\left[r^2(B_\phi ^2E^2)\right]\frac{1}{4\pi r}\left[\frac{r}{\psi }\right]\frac{U_\phi ^2B_p}{U_pF(\psi )}\\ \\ \frac{[U_pF(\psi )(1(r\mathrm{\Omega }/c)^2)B_p]}{4\pi rR_cF(\psi )}=0.\end{array}$$
(28)
Let us consider the flow in the region of radial expansion, as in the case under consideration, where the SC is small. For simplicity, we neglect the acceleration of the plasma. The first term in this equation decreases with $`r`$ as $`r^4`$. The second one decreases as $`r^2`$, since $`B_\phi ErB_p`$ (Bogovalov 1997a ). The third term is proportional to the product $`rU_\phi ^2B_p`$. This term decreases as $`r^3`$, since $`rB_p1/r`$ and $`U_\phi 1/r`$ due to angular momentum conservation. Thus, only the second and fourth terms survive in Eq. (28) at these large distances. Therefore, this equation takes the form

$$\frac{1}{8\pi x^2}\frac{}{\psi }(x^2(B_\phi ^2E^2))\frac{(U_p+F(\psi )x^2B_p)}{4\pi xR_cF(\psi )}=0,$$
(29)
where $`x=r\mathrm{\Omega }/c`$. From the frozen-in condition (2) we may estimate the toroidal magnetic field at large distances when $`v_\phi 0`$, where $`v_\varphi `$ is the azimuthal velocity.

$$B_\phi =\frac{Ec}{v}.$$
(30)
In the asymptotic regime plasma moves along straight lines, when $`R_c\mathrm{}`$. The structure of the flow in this regime is defined by the trivial equation

$$\frac{}{\psi }(x^2(B_\phi ^2E^2))=0.$$
(31)
Together with Eq. (30) this equation gives that the relativistic generalization of the Heyvaerts & Norman integral $`I_{HN}=|xB_\phi |/\gamma `$ (Chiueh et al. chiueh (1991)) should be constant in the region of radial expansion. This integral, shown in Fig. (8) by the two lower curves for essentially different distances, is actually constant in the small region $`0.01<\psi /\psi _{eq}<0.05`$ (Heyvaerts & Norman plateau in the figure). But at $`\psi /\psi _{eq}>0.05`$, $`I_{HN}`$ is not constant. This contradicts the expected behavior of $`I_{NH}`$, but, as we show below, this implies that the asymptotic regime is still not reached in this region of $`\psi `$.

The characteristic scale of collimation $`D_{coll}`$ which defines the distance where the flow reaches the asymptotic regime can be estimated as follows. $`d\theta /dR=1/R_c`$, where $`\theta `$ is the polar angle of the plasma velocity and $`R`$ is the distance to the center of the source. In the radially expanding wind $`d\psi =B_pr^2\mathrm{sin}(\theta )d\theta `$. Taking into account also that at large distances $`B_\phi =(r\mathrm{\Omega }/v_p)B_p`$ (Bogovalov 1997a ), the rate of turning of the trajectory of the plasma follows from Eq. (29) and is defined by the equation

$$|\frac{d\theta }{dR}|\frac{\mathrm{sin}2\theta }{RU_0^2(\frac{1}{\sigma }+\mathrm{sin}^2\theta )}.$$
(32)
Integration of this equation gives the estimate of the distance at which the line of the flow turns through the angle $`\mathrm{\Delta }\theta `$. The scale of the collimation of the flow line $`D_{coll}`$ is defined by the distance at which the flow line turns through the angle $`\mathrm{\Delta }\theta \theta `$. The integration of equation (32) can be performed from the fast mode surface since this consideration is valid only in the supersonic region. The scale of the collimation is estimated as

$$D_{coll}R_f\mathrm{exp}\left[\frac{\theta U_0^2(\frac{1}{\sigma }+\mathrm{sin}^2\theta )}{\mathrm{sin}2\theta }\right].$$
(33)
According to this equation there are two regions of the flow with the essentially different dependence of $`D_{coll}`$ on the parameters of the flow. The smallest $`D_{coll}`$ takes place for small angles limited by the value $`\theta ^2<1/\sigma `$. For these angles the scale of the collimation is (see also Bogovalov & Tsinganos bogts (1999))

$$D_{coll}R_f\mathrm{exp}\frac{1}{2\alpha ^2}.$$
(34)
Similar to the nonrelativistic outflow the collimation of the plasma close to the rotational axis occurs near the fast mode surface at $`\alpha >1`$. At the initially m-l outflow the magnetic flux depends on the angle $`\theta `$ as $`\psi =\psi _{eq}(1\mathrm{cos}\theta )`$. Therefore, the part of the magnetic flux where the collimation occurs as in the nonrelativistic limit is

$$\frac{\psi }{\psi _{eq}}\frac{1}{2\sigma }.$$
(35)
The collimation scale of the other part of the flow ( at low latitudes) is defined by the equation

$$D_{coll}R_f\mathrm{exp}\frac{\theta \mathrm{sin}\theta U_0^2}{2\mathrm{cos}\theta }.$$
(36)
and at the angles $`\theta ^22/U_0^2`$ the collimation scale becomes large compared with $`R_f`$ independent of $`\sigma `$. It follows from this result that the remarkable bend of the flow lines at low latitudes is possible only for the plasma with relatively low $`U_0`$. The collimation scale increases exponentially with the growth of the Lorentz factor of the plasma. These properties of the flow explain the numerical results.

In the simulation we obtain such a behavior of the flow, with three regimes: a core where the current increases with $`\psi `$, the Heyvaerts & Norman regime and the regime where the enclosed current increases (Fig. 8). Thus, an extremely narrow cylindrical core is formed which contains the magnetic flux $`\psi <(1/2\sigma )\psi _{eq}`$. This core is formed at a distance comparable to $`R_f`$ for $`\alpha >1`$, similarly to the nonrelativistic limit. For larger $`\psi `$ the collimation scale is defined by Eq. (36). For $`\gamma _0=7`$ the characteristic scale of collimation is of the order of $`R_f10^{10.6\theta \mathrm{tan}\theta }`$. It implies that the smaller $`\theta `$, the earlier the asymptotic regime will be achieved. This explains why the Heyvaerts & Norman plateau is formed only at the relatively small $`\theta `$ (small $`\psi `$). At larger $`\theta `$ the scale of collimation exponentially increases and finally diverges at $`\theta =\pi /2`$. Thus, the asymptotic regime cannot be achieved at low latitudes. At the final distance from the source there should exist a region at low latitudes (large $`\psi `$) where $`I_{NH}`$ is not constant. It is this behavior which is found in our calculations.

For $`\alpha =4`$ the collimation is more efficient and the effect of collimation appears already in the nearest zone where the plasma is not yet accelerated. At the estimation of the collimation scale it is necessary to take into account acceleration of the plasma to a Lorentz factor greater than 10 near the equator. This could be done by integration of Eq. (32) from a radius $`R_{in}>R_f`$ where the energy of the plasma practically reaches the terminal value. If we take into account the acceleration, the asymptotic regime (in which $`I_{NH}`$ does not depend on $`\psi `$ in a radially expanding wind) will be achieved at a distance of the order $`D_{coll}>R_f10^{22\theta \mathrm{tan}\theta }`$. It follows from this expression that at $`\alpha =4`$ the asymptotic regime at $`R10^4R_f`$ will be achieved only for $`\psi /\psi _{eq}<0.1`$. This value agrees well with the numerical results presented in Fig. 8 in the limits of uncertainties.

Comparison of the present work with our previous work (Bogovalov & Tsinganos bogts (1999) ) shows that the relativistic plasma is collimated less efficiently than the nonrelativistic one. In particular, the fraction of the collimated mass flux in the nonrelativistic flow was close to $`1\%`$, while in the relativistic limit it is of the order $`0.1\%`$. This difference can be understood from simple qualitative consideration. The dynamics of the cold plasma is basically controlled by two forces. One of them is the volumetric Lorenz force, $`(𝐣\times 𝐁)/c`$ and the other is the volumetric Coulomb force, $`q𝐄`$ (see Eq.(1)). The collimation of the plasma to the axis of rotation is defined by the components of these forces perpendicular to the stream lines. The perpendicular component of the Lorenz force is $`F_L=\frac{1}{c}j_{}B_\phi `$, where $`j_{}`$ is the component of the poloidal electric current parallel to the poloidal magnetic field and $`B_\phi `$ is the toroidal magnetic field. The perpendicular component of the Coulomb force is $`F_C=qE_{}`$. These forces are shown schematically in Fig. (4). It is important to note that these two forces are directed along opposite directions. The Lorenz force collimates the plasma while the Coulomb force decollimates it. In the nonrelativistic limit the Coulomb force is neglected, since it is negligible in comparison to the Lorenz force. The ratio of the Coulomb and Lorenz forces is of the order $`(v/c)^2`$ in this limit. In the relativistic limit the Coulomb force becomes practically equal to the Lorenz force and it strongly decreases the SC of the plasma. That is why the relativistic generalization of the parameter $`\alpha `$ introduced earlier for the nonrelativistic winds (Bogovalov & Tsinganos bogts (1999)) has additionally the Lorentz factor in the denominator (see Eq. (23)).

Although the Lorentz factor in the denominator of Eq. (23) decreases remarkably the value of $`\alpha `$, it was expected that the collimation of the relativistic plasma will be strong enough at the condition $`\alpha >1`$. In the nonrelativistic case the collimation becomes strong already in the subsonic region under this condition. The most surprising result is that the relativistic plasma is practically not collimated, even under this condition. The increase of the collimating force does not provide more efficient collimation of the relativistic plasma! This happens since another relativistic effect comes into play.

The curvature $`\kappa _c=1/R_c`$ of the stream lines of the nonrelativistic plasma is defined by the simple equation

$$\kappa _c=\frac{F_{}}{mnv^2},$$
(37)
where $`m`$ is mass of particles, $`F_{}`$ is the sum of all forces affecting the plasma perpendicular to the stream line. It is seen from this equation that an increase of the force $`F_{}`$ results in an increase of $`\kappa _c`$, such that collimation becomes stronger. In other words, as it is evident from Eq. (37), $`\kappa _c`$ is proportional to the collimating force and inversly proportional to the mass density of the outflowing plasma. The situation changes essentially in the case of the relativistic, Poynting dominated plasma ($`\sigma >1`$) where the energy flux is dominated by the electromagnetic field. Then, according to the familiar relationship between energy and mass ($`E=mc^2`$), the flux of the energy is equivalent to some inertial mass flux. It follows from Eq.(29) that eq. (37) for the relativistic Poynting-dominated plasma is modified as follows

$$\kappa _c=\frac{F_{}}{nv^2\gamma (m+\frac{B_\phi ^2}{4\pi nc^2})}.$$
(38)
This equation demonstrates that the effective mass of the Poynting- dominated plasma is dominated by the toroidal magnetic field. It is important that now the effective mass increases by a term proportional to $`B_\phi ^2`$. The collimating force is also proportional to the same term (see eq. (29)). Therefore, the increase of the collimating force is accompanied by a corresponding increase of the effective mass density of the plasma. This in turn results in the fact that the curvature $`\kappa _c`$ of the stream lines does not increase with the increase of the collimating force. This is the reason why the collimation scale (36) at low latitudes does not depend on the magnetic field and angular velocity of the source.

## 6 Conclusion

This work shows that the self-collimation (SC) process occurs in relativistic winds as predicted by the SC theory. The wind consists of a radially expanding part which surrounds a cylindrically collimated flow. The radius of the cylindrically collimated flow is $`R_jR_f/\alpha `$, where $`R_f`$ is the initial radius of the fast magnetosonic surface. The mass flux in the collimated part is small compared with the total mass flux in the wind ($`10^3`$), even for relatively small Lorentz factors. SC essentially decreases for larger Lorentz factors. Acceleration of the plasma occurs in the radially expanding wind while the plasma in the jet is not accelerated in this cold outflow model.

The relatively small part of the mass flux in the cylindrically collimated flow has been also found in our simulations of nonrelativistic outflows (Bogovalov & Tsinganos bogts (1999)), although the magnetic flux contained in the relativistic case appears smaller by a factor of order 10. In addition, the collimation scale at low latitudes does not depend on the magnetic field and the angular velocity of the central source for $`\alpha >1`$, in contrast with our results for the nonrelativistic winds in the same model. This means that the weak collimation of the relativistic plasma is due to the intrinsic property of the plasma and not to a property of the model. Collimation decreases due to the affect of the electric force which becomes comparable with the collimating Lorenz force in the relativistic

limit and in the high $`\sigma `$ limit the collimation is saturated due to the contribution of the electromagnetic field into the inertial mass of plasma.

The most striking feature of the collimated flows obtained in the model under consideration is that a very small fraction ($`10^3`$) of the mass flux exists in the collimated flow in comparison to the mass flux in the wind. This result is apparently in conflict with observations of AGN. The source of the wind energy is the gravitational energy released by the accreted matter. Therefore, the energy flux in the wind should be less than the luminosity of the accretion disk where about half of the gravitational energy is released (Shakura & Sunyaev ss (1973)). The total luminosity of AGN $`L_{tot}`$ in soft electromagnetic emission ( below X-rays) should be comparable with the energy released at the accretion. According to our results the energy flux in the jet $`\dot{E}_{jet}`$ should be by the factor $`10^3`$ less than the accretion energy, $`\dot{E}_{jet}<10^3L_{tot}`$. This disagrees with observations. X-ray and gamma-ray emission is usually attributed to the emission from jets (Urry & Padovani urry (1995)). Sometimes, the luminosity of the jets is less only by a factor of 10 than the total luminosity of AGN at other wavelengths (Elvis et al. elvis (1994)). The energy flux in this emission can be considered as the lower limit on the energy flux in the jet (however, we do not take into account beaming here). If so, we have in AGN $`\dot{E}_{jet}>10^1L_{total}`$.

This disagreement between theory and observations implies that either the present model has no observed counterpart or the SC mechanism is not able in principle to provide collimation of the major fraction of the outflow into the jet. It is important in this regard to understand from the theoretical point of view what modifications should be made in order to provide magnetic collimation of the major fraction of the relativistic wind. This problem has been already investigated by Bogovalov & Tsinganos (bogts00 (2001)) in relation to nonrelativistic winds. It was shown in this study that in a model consisting of a central source ejecting a high-speed wind and an accretion disk ejecting a nonrelativistic magnetised outflow, all the mass flux from the central source can be cylindrically collimated indirectly by the surrounding disk-wind. Apparently, a collimation of the major part of a relativistic outflow from a central source can be provided by the same mechanism in the context of the present model as well.

###### Acknowledgements.

I acknowledge support from the University Joseph Fourier and the Institut Universitaire de France as well as the warm hospitality of Prof. G. Pelletier and Dr. J. Ferreira during my stay in the Laboratoire d’Astrophysique de l’Observatoire de Grenoble, where part of this work was done. This work was partially supported by the Ministry of Education of Russia in the framework of the program “Universities of Russia - fundamental research”, project N 990479 and by collaborative INTAS-ESA grant N 99-120. The author is sincerely grateful to the unknown second referee and K.Tsinganos for the large amount of work they put in reading of this manuscript and their addition of useful comments. The author thanks V.B. Komberg for useful discussion of the observational properties of AGN.

## Appendix A Details of the numerical simulation in the nearest zone.

For the numerical solution of the time dependent problem, the system of equations should be presented as

$$\frac{}{t}\widehat{\mathrm{\Phi }}=A(\widehat{\mathrm{\Phi }},\widehat{\mathrm{\Phi }}_i,\widehat{\mathrm{\Phi }}_{ik}),$$
(39)
where $`\widehat{\mathrm{\Phi }}`$ is the vector of the plasma state which includes the density, components of the velocities and magnetic fields. $`\widehat{\mathrm{\Phi }}_i`$ and $`\widehat{\mathrm{\Phi }}_{ik}`$ are the first and the second order derivatives of $`\widehat{\mathrm{\Phi }}`$ on the spatial variables. Eq. (1) contains the time derivative $`\dot{𝐄}`$ in the right hand part generated by the displacement current in Eq. (3). Therefore, some special transformation of Eq. (1) is necessary to get form (39). The specific form of this transformation is not unique. In our previous work (Bogovalov 1997a ) unsuccessful choice of this transformation resulted in the instability of the numerical code in the Poynting flux-dominated regime. Here we use another form of the equation of motion which provides stable numerical code in all regimes.

This equation can be obtained as follows. Time derivative of the frozen-in condition (2) gives that

$$\dot{𝐄}=\frac{1}{\gamma }(\dot{𝐔}\times 𝐁)\frac{1}{\gamma }(𝐔\times \dot{𝐁})\frac{\dot{\gamma }}{\gamma }𝐄.$$
(40)
Substitution of Eq. (40) and Eq. (3) together with the relativistic relationship $`\dot{\gamma }=v_i\dot{U}_i/c`$ in Eq. (1) gives the equation of motion in the form

$`M_{ik}\dot{U}_k=\rho (𝐯)U_i+qE_i/c+{\displaystyle \frac{1}{4\pi c}}((\times 𝐁)\times 𝐁)_i`$ (41)
$`{\displaystyle \frac{1}{4\pi c^3}}(v_i(𝐁\dot{𝐁})\dot{B}_i(𝐁𝐯)),`$
where

$$M_{ik}=(\rho +\frac{B^2}{4\pi c^2\gamma })\delta _{ik}\frac{1}{4\pi c^2\gamma }(B_iB_k+S_iv_k/c)$$

and $`S_i=(𝐄\times 𝐁)_i`$. Terms with $`\dot{𝐁}`$ in the right hand part of Eq. (41) are eliminated with help of Eq. (4). After that, multiplication of equation of motion (41) on the inverse matrix

$`M_{ik}^1={\displaystyle \frac{1}{(\rho +\frac{B^2}{4\pi c^2\gamma })}}[\delta _{ik}+`$ (42)
$`{\displaystyle \frac{B_iB_k}{4\pi \rho c^2\gamma }}+{\displaystyle \frac{(𝐯𝐁/c)S_iB_k/4\pi \rho c^2\gamma +S_iv_k/c}{(4\pi \rho c^2\gamma +B^2(𝐒𝐯/c)}}]`$
gives the equation of motion in the form (39) appropriate for the numerical solution.

A two-step Lax-Wendroff numerical scheme was used for the numerical solution of the problem (Press et al. press (1990)). The analytical solution obtained for the slow rotators (Bogovalov 1997a ) was used as the initial state of the plasma.

The testing of the numerical code was performed by comparison of the results with the results obtained with the nonrelativistic version of the code (Bogovalov & Tsinganos bogts (1999)), by comparison with predictions of the analytical theory (Bogovalov bog92 (1992), Bogovalov 1997a ) and by simulation of the perturbation propagation in the wind. One of these tests is presented in Fig. 9. The beginning of the star rotation starts the shock wave propagating with the Alfvenic velocity in the isotropic wind. This velocity is very close to 1 for the taken parameters of the wind. The shock has an oscillating structure typical for the collisionless shocks. The front is smoothed due to the numerical diffusion. The first (largest) maximum can be accepted as the boundary of the shock. It follows from the figure that the wave propagates with the velocity very close to the light velocity which is equal $`c/v_0=1.01`$ in the dimensionless variables normalized to $`v_0`$ at $`\gamma _0=7`$.


text:
# Phase Diagram of a Two-Dimensional Neutral Classical Coulomb Gas from a Non-perturbative sine-Gordon Expansion

## ACKNOWLEDGMENTS

Lu acknowledges Dr. S. Y. Cho, Prof. M. H. Chung and Mr. M. D. Kim for useful discussions. This project was supported by the Korea Research Foundation (99-005-D00011). Lu’s work was also supported in part by the National Natural Science Foundation of China under the grant No. 19875034.


text:
# Tractable approximations for probabilistic models: The adaptive TAP mean field approach

## Abstract

We develop an advanced mean field method for approximating averages in probabilistic data models that is based on the TAP approach of disorder physics. In contrast to conventional TAP, where the knowledge of the distribution of couplings between the random variables is required, our method adapts to the concrete couplings. We demonstrate the validity of our approach, which is sofar restricted to models with non-glassy behaviour, by replica calculations for a wide class of models as well as by simulations for a real data set.

Probabilistic models (for a review see e.g. ) find widespread applications in many areas of data modeling. Their goal is to explain complex observed data by a set of unobserved, hidden random variables based on the joint distribution of both sets of variables. The price that a modeler has to pay for the high degree of flexibility of these models is the vast increase in computational complexity when the number of hidden variables is large.

Both statistical inference about hidden variables and training usually require computation of marginal distributions of the hidden variables which for exact calculation demands infeasible high dimensional sums or integrals. Since similar types of calculations are ubiquitous in the computations of thermal averages, there is a great deal of interest in adopting approximation techniques from statistical physics. For a variety of cases, when a standard tool, the Monte Carlo sampling technique reaches its limits, a simple mean field (MF) method, which neglects correlations of random variables has yielded good results in a variety of probabilistic data models. The MF approximation yields a closed set of nonlinear equations for the approximate expectation values of random variables which usually can be solved in a time that only grows polynomially in the number of variables. At present, there is a growing research activity trying to overcome the limitations of the simple MF method by partly including the dependencies of variables but still keeping the approximation tractable (for a review see ).

Various researchers  have discussed applications of the so-called TAP MF approach, originating in the statistical physics of disordered systems, first introduced by Thouless, Anderson and Palmer (TAP)  to treat the Sherrington-Kirkpatrick (SK) model of disordered magnetic materials . Under the assumption that the couplings (interactions) between random variables are themselves drawn at random from certain classes of distributions, the TAP equations become exact in the thermodynamic limit of infinitely many variables. Unfortunately, the Onsager correction to the simple, naive MF theory will explicitly depend on the distribution of these couplings. Two models with the same connectivities but different distributions for the couplings, like e.g. the SK model and the Hopfield model  have different expressions for the Onsager corrections (see e.g. , chapter XIII).

In order to use the TAP method as a good approximation for models of real data, the lack of knowledge of the underlying distribution of the couplings (which are usually functions of the observed data) should be compensated by an algorithm which adapts the Onsager correction to the concrete set of couplings. Simply taking the correction from a theory that assumes a specific distribution may lead to suboptimal performance. This letter presents a solution to this problem for an important class of probabilistic models. As a check of the validity of the approach, we show that our method leads to the exact results in the thermodynamic limit for large classes of probability distributions over the couplings.

We will consider probabilistic models of the type

$$P(𝐒)=\frac{\rho (𝐒)}{Z(𝜽,𝐉)}\mathrm{exp}\left[\underset{i<j}{}S_iJ_{ij}S_j+\underset{i}{}S_i\theta _i\right]$$
(1)
where the set $`𝐒=(S_1,\mathrm{},S_N)`$ denotes the (hidden) random variables of the model. Any observed (i.e. fixed) quantities are assumed to be encoded in the matrix $`𝐉`$ and the fields $`𝜽`$. The term $`\rho (𝐒)_j\rho _j(S_j)`$ is a product distribution which also contains all constraints of the $`S_i`$ (the range, discreteness, etc). In its simplest version, when $`𝐒`$ is a real variable with positive measure $`\rho `$, the class of models (1) contains Ising models (such the SK and Hopfield models), Gaussian process models , probabilistic independent component analysis  and combinatorial optimization problems . If we lift the restrictions that all variables must be real random variables, we can treat a variety of important models with dependencies between the $`S_i`$ that are defined through a set of fields $`_{i=1}^Nx_{ij}S_i`$. We will give two examples. Bayesian learning in single layer neural networks is described by a Gibbs distribution $`P(𝐒)P_0(𝐒)_{j=1}^mF(_{i=1}^Nx_{ij}S_i),`$ where $`𝐒`$ is a weight vector of the network being trained on a number of $`m`$ data vectors with components $`x_{ij}`$ in a $`N`$ dimensional space. $`P_0`$ is a prior distribution of the weights and $`F`$ is the Likelihood quantifying the goodness of fit to the data . A second example is given by the class of Bayesian belief networks on a directed graph which are promising models for adaptive expert systems. They are defined by $`P(𝐒)=_iP(S_i|\text{pa}(S_i))`$ where $`S_i\{0,1\}`$ and $`pa`$ denotes the parents of $`S_i`$, i.e. the variables in the graph that feed their information into $`S_i`$ via directed bonds. A specific type is the sigmoid belief networks , where $`P(S_i|\text{pa}(S_i))=\frac{e^{S_ih_i}}{1+e^{h_i}}`$ with $`h_i=_{jpa(S_i)}x_{ij}S_j`$. The latter two models can be easily brought into the form (1) by the standard ‘field-theoretic’ trick of introducing Dirac $`\delta `$-functions and their exponential representations using purely imaginary conjugate variables $`\widehat{𝐒}=(\widehat{S}_1,\mathrm{},\widehat{S}_m)`$. This leads to an augmentation of the space of variables to the set $`(𝐒,\widehat{𝐒})`$. The hatted variables have the complex single variable distributions $`\widehat{\rho }(\widehat{S})=\frac{dh}{2\pi i}e^{\widehat{S}h}e^{H(h)}`$ in case of the neural network model and $`\widehat{\rho }(\widehat{S})=\frac{dh}{2\pi i}e^{\widehat{S}h}/(1+e^h)`$ for the belief network (where $`m=N`$). The augmented coupling matrix is of the form $`𝐉=\left(\begin{array}{cc}𝐀& 𝐁\\ 𝐁^T& \mathrm{𝟎}\end{array}\right),`$ where $`B_{ij}=x_{ij}`$ and $`𝐀=0`$ for the neural network and $`A_{ij}=B_{ij}=x_{ij}`$ for the belief net.

We will derive both an adaptive TAP-like approximation for the marginal distribution $`P_i(S)_{ji}dS_jP(𝐒)`$ and the free energy $`F(𝐉,𝜽)=\mathrm{ln}Z(𝐉,𝜽)`$. The free energy corresponds to the negative log probability of the observed data which can be used as a yardstick for deciding which model best fits to the data.

Our derivation will be based on the cavity approach introduced by . We will assume that we are not dealing with a glassy system with its many ergodic components, but that all averages are for a single state. This is (as shown for many of the teacher-student scenarios studied in the statistical mechanics of neural networks) usually expected to hold when the probabilistic model is well matched to the data. Defining the field $`h_i=_jJ_{ij}S_j`$, the marginal distribution of $`S_i`$ can be written as

$`P_i(S)={\displaystyle \underset{ji}{}dS_jP(𝐒)}={\displaystyle \frac{\rho _i(S)}{Z_i}}e^{H_i(S)},`$ (2)
where we have introduced an effective single variable Hamiltonian $`H_i(S)`$ with corresponding partition function $`Z_i`$. Defining an auxiliary average over the distribution of the system with variable $`S_i`$ left out by $`\mathrm{}_{\backslash i}`$, we get

$$H_i(S)=\mathrm{ln}e^{Sh_i}_{\backslash i}=\underset{k}{}\frac{\kappa _k^{(i)}}{k!}S^k$$
(3)
where $`\kappa _k^{(i)}`$ are the cumulants of this cavity distribution, i.e. $`\kappa _1^{(i)}=h_i_{\backslash i}`$ and $`\kappa _2^{(i)}=h_i^2_{\backslash i}h_i_{\backslash i}^2`$ etc.

The basic physical assumption, which is the major ingredient of all cavity derivations of the TAP mean field theory , is that all variables $`S_j`$ have only weak mutual dependencies. Mathematically expressed within the so–called clustering hypothesis  this becomes equivalent to the vanishing of all cumulants $`\kappa _k^{(i)}`$ with $`k>2`$ for fully connected systems. In the case, where the $`S_j`$ are real variables with positive measure, this corresponds to a central limit theorem for the cavity fields. Under this assumption, setting $`V_i=\kappa _2^{(i)}`$, we get

$`h_i`$ $`=`$ $`{\displaystyle \frac{1}{Z_i}}{\displaystyle 𝑑S\rho _i(S)\frac{}{S}e^{H_i(S)}}=h_i_{\backslash i}+V_iS_i`$ (4)
$`P_i(S)`$ $`=`$ $`{\displaystyle \frac{\rho _i(S)}{Z_i}}e^{\left(_jJ_{ij}S_jV_iS_i+\theta _i\right)S+\frac{1}{2}V_iS^2}`$ (5)
for $`i=1,\mathrm{},N`$. So far, the approach is well known. The new aspect of our paper is in the way we compute the $`V_i`$’s. Since these reaction terms account for the weak influence between random variables, they can be computed self-consistently from the matrix of susceptibilities $`\chi _{ij}\frac{S_i}{\theta _j}`$. We make the approximation that upon differentiation, the $`V_i`$’s are held constant which is consistent with the fact that the $`V_i`$’s are expected to be selfaveraging quantities in the thermodynamic limit. Under this assumption we get from eq. (5) $`\chi _{ij}=\chi _{ii}\left(\delta _{ij}+_k(J_{ik}V_k\delta _{ik})\chi _{kj}\right)`$ which can be solved with respect to $`𝝌`$ and yields $`𝝌=(𝚲𝐉)^1,`$ where $`𝚲=\text{diag}\{V_i+1/\chi _{ii}\}`$ is a diagonal matrix. The Fluctuation Dissipation Theorem (again assuming that we deal with a single state), shows that $`𝝌`$ also equals the matrix of correlations $`C_{ij}=S_jS_kS_jS_k`$. By specializing to the diagonal elements, we can compute $`V_i`$ as a function of $`S_i^2S_i^2`$ by solving

$$S_i^2S_i^2=\frac{^2\mathrm{ln}Z_i}{\theta _i^2}=\left[\left(𝚲𝐉\right)^1\right]_{ii}$$
(6)
for $`i=1,\mathrm{},N`$. The sets of equations (5) together with (6) constitute the first main result of this letter. They yield closed sets of equations for the first and second moments of $`S_i`$ which in turn enables us to approximate the full marginal distribution of $`S_i`$ and the correlation functions. For comparison we note that the naive mean field approximation (for real random variables) is obtained by setting $`V_i=0`$. Selfinteractions $`V_iS_i`$ determined by the the linear response method have also been introduced in  as a heuristics to correct the naive MF equations for Boltzmann machines. A sanity check of the internal consistency of our approach is obtained by the fact that the matrix $`𝝌`$ must be positive definite. (If a group of the variables are complex, this has to hold for the submatrix of the real random variables).

The next task is to compute the adaptive TAP approximation to the free energy $`F(𝐉,𝜽)=\mathrm{ln}Z(𝐉,𝜽)`$. It is useful to generalize our model eq. (1) to a one parameter class of models where the interaction $`𝐉`$ is replaced by $`s𝐉`$ with $`0s1`$ and to define the Legendre transform (Gibbs free energy) by

$$\mathrm{\Phi }_s(𝐦,𝐌)=F(s𝐉+𝝀,𝜽+𝜸)+\underset{i}{}\gamma _im_i+\underset{i}{}\frac{\lambda _i}{2}M_i,$$

where $`\gamma _i`$ and $`\lambda _i`$ are external fields conjugate to $`S_i`$ and $`S_i^2`$ which must be chosen to extremize the right hand side and $`𝝀`$ is a diagonal matrix with entries $`\lambda _i`$. The solutions $`𝐦^e`$ and $`𝐌^e`$ of the sets of equations $`_{m_i}\mathrm{\Phi }_s=_{M_i}\mathrm{\Phi }=0,`$ determine the correct equilibrium expectation values $`S_i_s=m_i^e`$ and $`S_i^2_s=M_i^e`$ (the index indicates that the expectation is taken with parameter $`s`$). Our desired approximation to the free energy is finally obtained as $`F(𝐉,𝜽)=\mathrm{\Phi }_1(𝐦^e,𝐌^e)`$. To compute $`\mathrm{\Phi }_1`$ we differentiate $`\mathrm{\Phi }_s`$ with respect to $`s`$, to show that

$$\mathrm{\Phi }_1=\mathrm{\Phi }_0\frac{1}{2}_0^1𝑑s\left\{\underset{i,j}{}m_iJ_{ij}m_j+Tr(𝝌_s𝐉)\right\}$$
(7)
with $`\chi _{s,ij}=S_iS_j_sS_i_sS_j_s`$. Inserting our TAP approximation $`𝝌_s=(𝚲_ss𝐉)^1`$ and integrating, we obtain

$`\mathrm{\Phi }_1`$ $`=`$ $`\mathrm{\Phi }_0{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}m_iJ_{ij}m_j+\mathrm{\Delta }\mathrm{\Phi }`$ (8)
$`\mathrm{\Delta }\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}det(𝚲𝐉){\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}V_i\chi _{ii}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\mathrm{ln}\chi _{ii}`$ (9)
with $`\chi _{ii}=M_im_i^2`$. The first two terms constitute the naive mean field approximation to $`\mathrm{\Phi }`$ and the last term $`\mathrm{\Delta }\mathrm{\Phi }`$ is the Onsager correction. Note, that this result is not equivalent to a truncation of a power series expansion of $`\mathrm{\Phi }`$ to second order in $`s`$ (a Plefka expansion ) but contains terms of all orders. A different way to derive this result is obtained from the observation that the functional form of the Onsager term $`V_i`$ in the TAP equations does not depend on the specific single variable densities $`\rho (𝐒)`$. Hence, we may compute this universal form by calculating $`\mathrm{\Phi }`$ for an exactly solvable model, i.e. for a Gaussian $`\rho `$ and subtract the naive mean field part. This is related to the strategy used by Parisi and Potters  in order to derive the TAP equations for a spin glass model with orthogonal random matrix $`𝐉`$.

To check the significance of our approach, we will next show that it will give the correct results for the statistical mechanics in the thermodynamic limit $`N\mathrm{}`$ for a large class of distributions of the random matrix $`𝐉`$. For simplicity, we specialize to models with only one type of single variable distribution $`\rho _i(S)=\rho (S)`$. Selfaveraging properties of the models can be computed within the replica framework by averaging the free energy over the distribution of the random matrix $`𝐉`$. This requires the calculation of the asymptotic scaling of the function $`K_N(𝐀)\frac{1}{N}\mathrm{ln}\left[e^{\frac{1}{2}Tr(\mathrm{𝐀𝐉})}\right]_𝐉`$ for the matrix $`A_{ij}=_{a=1}^nS_{ia}S_{ja}`$, where the $`𝐒_a`$, are $`n`$ replicas of the variables. Following Ref.  and assuming the scaling $`K_N(𝐀)TrG(𝐀/N)`$ as $`N\mathrm{}`$ where the function $`G`$ characterizes the random matrix ensemble, the averaged free energy will depend only on the single set of orderparameters given by $`q_{ab}\frac{1}{N}_iS_{ia}S_{ib}`$. This is characteristic for models with matrices $`𝐉`$ of extensive connectivity. E.g., the SK-model with coupling matrix of independent components of variance $`\frac{\beta }{N}`$ has $`G(r)=\frac{(\beta r)^2}{4}`$ and the Hopfield model with $`J_{ij}=_{\mu =1}^{\alpha N}x_i^\mu x_j^\mu `$ and independent $`x_i^\mu `$ with variance $`\frac{\beta }{N}`$ leads to $`G(r)=\frac{\alpha }{2}(\mathrm{ln}(1\beta r)+\beta r)`$. Under the assumption of replica symmetry, the averaged free energy $`f=\frac{1}{N}[\mathrm{ln}Z]_𝐉`$ is obtained by extremizing

$`f(q,\mathrm{\Delta })`$ $`=`$ $`G(\mathrm{\Delta })+\mathrm{\Delta }(qG{}_{}{}^{\prime \prime }(\mathrm{\Delta })+G^{}(\mathrm{\Delta }))`$ (10)
$`{\displaystyle Dz\mathrm{ln}}`$ $`dS`$ $`\rho (S)\mathrm{exp}\left[\sqrt{2qG{}_{}{}^{\prime \prime }(\mathrm{\Delta })}zS+G^{}(\mathrm{\Delta })S^2\right]`$ (11)
with respect to the off-diagonal orderparameter $`q=q_{ab}`$ and to $`\mathrm{\Delta }=q_{aa}q`$, where $`Dz=\frac{dz}{\sqrt{2\pi }}e^{z^2/2}`$.

We can show the correspondence for $`N\mathrm{}`$ of the adaptive TAP method and replica theory. A disorder average gives the conventional TAP result for the Onsager coefficients: $`V_i=V=2G^{}(\overline{\chi }),`$ with $`\overline{\chi }=\frac{1}{N}_i[\chi _{ii}]_𝐉`$. To compare the TAP Gibbs free energy eq. (8) with the replica symmetric free energy (10), we compute $`\widehat{f}=lim_{N,\gamma \mathrm{}}\frac{1}{\gamma N}\left[\mathrm{ln}𝑑𝐦𝑑𝐌\mathrm{exp}(\gamma \mathrm{\Phi })\right]_𝐉,`$ where the paths of integration must be chosen such that the integral converges. The integral will be dominated by the values for $`𝐦`$ and $`𝐌`$ which fulfill the TAP equations. Evaluating this expression using the replica method shows that both free energies coincide, i.e. $`\widehat{f}=f`$. It is also possible to translate the condition of positive definiteness of the susceptibility matrix $`𝝌`$ into the thermodynamic limit. We can show that this stability is satisfied for $`12G^{\prime \prime }(\overline{\chi })\frac{1}{N}_i[\chi _{ii}]_𝐉^2>0,`$ which coincides with the well known AT stability condition of replica theory .

We have performed two types of simulations of the TAP approaches on Bayesian neural network learning problems. For the first case (Fig. 1) we test the self–consistency of our method on a real data set, ‘Sonar – Mines versus Rocks’  of size $`m=104`$ with with binary class labels $`y_j=\pm 1`$ and a $`N=60`$ dimensional input space. The prior is $`P_0(𝐒)\mathrm{exp}(𝐒𝐒/2)`$ and the likelihood $`F(\widehat{h}_j)=\varphi (y_i\widehat{h}_j/\sigma )`$, with $`\widehat{h}_j=_ix_{ij}S_i`$, $`\varphi (t)=_t^{\mathrm{}}Dz`$ and $`\sigma ^2=0.5`$. We compute the prediction for the average (conjugate) cavity field $`\widehat{h}_j_{\backslash j}=\widehat{h}_j\widehat{V}_j\widehat{S}_j`$, using eq. (4). The fraction of negative terms $`y_j\widehat{h}_j_{\backslash j}`$ equals the ‘leave-one-out’ estimate $`ϵ_{\mathrm{loo}}`$ which provides an important practical estimator for the generalization error of the network. If our theory takes the reaction of the remaining variables correctly into account, this prediction should be close to the ‘exact’ average cavity field obtained by leaving one example out and solving the TAP eqs. for the remaining $`m1`$ examples. Fig. 1 shows excellent agreement between the two computations and we find $`ϵ_{\mathrm{loo}}=ϵ_{\mathrm{loo}}^{\mathrm{exact}}=33/104`$. For comparison, the conventional TAP approach , which assumes a distribution of input data vectors with independent components, leads to a wrong result, $`ϵ_{\mathrm{loo}}=41/104`$ and $`ϵ_{\mathrm{loo}}^{\mathrm{exact}}=33/104`$.

In the second set of simulations Fig. 2 we demonstrate that the adaptive TAP method yields the correct statistical physics for the case of the linear Ising perceptron . This has prior distribution $`P(S)=\frac{1}{2}\delta (S1)+\frac{1}{2}\delta (S+1)`$ and likelihood $`F(\widehat{h}_j)\mathrm{exp}((y_i\widehat{h}_j)^2/2\sigma ^2)`$, where we have chosen $`\sigma ^2=0.2`$ and $`N=60`$ in the simulations. See Ref.  for a discussion of this model in the context of demodulation in communications systems. To compare with the replica results , we have generated inputs at random and outputs using a noise free teacher perceptron sampled from the same prior. The small deviations between theory and TAP simulations close to the first order transition are attributed to hysteresis effects.

It will be interesting to extend our adaptive TAP method to glassy systems (generalizing the ideas of chapter V in ) where the present approach would fail, e.g. indicated by the appearance of negative eigenvalues in the susceptibility matrix. However, one may speculate that in such cases solving the TAP equations may be highly nontrivial. Acknowledgments. We thank T. Tanaka for providing us with his preprint . This research is supported by the Swedish Foundation for Strategic Research as well as the Danish Research Councils through the Computational Neural Network Center (CONNECT) and the THOR Center for Neuroinformatics.