A&A 381, 1066-1079 (2002)
DOI: 10.1051/0004-6361:20011533

Stationary equatorial MHD flows in general relativity

F. Daigne - G. Drenkhahn

Max-Planck-Institut für Astrophysik, Postfach 1317, 85741 Garching bei München, Germany

Received 7 May 2001 / Accepted 29 October 2001

Abstract
We derive a new formulation of the fully general relativistic equations describing a stationary equatorial MHD outflow from a rotating central object. The wind solution appears as a level contour of a "Bernoulli'' function fixed by the requirements that it must pass through the slow and fast critical points. This approach is the general relativistic extension to the classical treatment of Sakurai (1985). We discuss in details how the efficiency of the magnetic to kinetic energy conversion depends mainly on the geometry of the flux tubes and show that the magnetic acceleration can work very well under some conditions. We show how this tool can be used for the study of several astrophysical phenomena, among which gamma-ray bursts.

Key words: MHD - black hole physics - relativity - gamma-ray bursts

1 Introduction

Magnetized winds are believed to be present in many astrophysical objects. They were first put forward in the context of the solar wind. More recently they were discovered to probably play a major role in many situations where a relativistic flow is powered by a central rapidly rotating compact object: wind of pulsars, jets in radio galaxies, quasars, Seyfert galaxies and BL Lac objects, microquasars and even possibly gamma-ray bursts.

The first quantitative model of a magnetic stellar wind was developed by Weber & Davis (1967). They were considering the equations of a stationary, axisymmetric, polytropic flow near the equatorial plane in classical MHD. They found that such a wind can carry off most of the angular momentum of the star and are very efficient to accelerate particles up to very high velocities. An important feature of magnetic winds is the existence of three "critical points'' where the velocity of the flow equals the wave velocity of the three MHD wave modes (the slow, Alfvén and fast modes) whereas in comparison non-magnetic winds (Parker 1958) have only one "critical point'' where the velocity of the flow equals the sound speed (sonic waves being the only present wave mode).

The first extension of this theory to relativistic winds is due to Michel (1969) in the context of radio pulsars. This work was considering cold outflows driven by rapidly rotating highly magnetized neutron stars. The main conclusion was that the efficiency of the magnetic to kinetic energy conversion was extremely low compared to the classical case. Goldreich & Julian (1970) studied cool isothermal relativistic winds. Kennel et al. (1983) extended Michel's model to finite temperatures and relativistic injection speeds. All these works were always limited to the equatorial plane and either completely neglected the effect of gravity or adopted an approximative treatment for it. Okamoto (1978) first included an exact general relativistic description of the gravity field. His work was also not limited to the equatorial plane, applying for that the powerful concept of flux tubes. However he restricted his study to pressureless flows only. In a series of papers Camenzind (1986a,b, 1987 derived a complete set of equations describing a stationary axisymmetric relativistic magnetic wind in an arbitrary metric. He then solved these equations in some particular cases (cold flows, jet geometries).

The goal of this paper is to present a formulation of the equations governing a stationary axisymmetric MHD flow in the equatorial plane including an exact treatment of all effects (thermal pressure, gravity and arbitrary shapes of flux tubes) which allows a direct comparison with the classical model of Weber & Davis (1967), so that the relativistic effects can be easily identified. This is done in Sects. 2 and 3, where we worked by analogy with the formulation of the classical case by Sakurai (1985). Then we study in details the efficiency of the magnetic to energy conversion (Sect. 4), in particular the influence of the flux tubes geometry and of the gravity. We confirm and extend the results of Begelman & Li (1994) and we show that a large variety of situations is expected from very inefficient winds like those considered by Michel (1969) to highly efficient cases. Because our model assumes axisymmetry and focus on the equatorial plane, it fully applies only to simple astrophysical objects like isolated neutron stars. On the other hand it can also describe the outer parts of more complex objects, e.g. compact objects with accretion disks or complex magnetospheres, as long as the magnetic field can be approximated as monopole like at these distances from the source. In the particular case of gamma-ray bursts the possibility of Poynting-flux dominated fireballs is briefly discussed in Sect. 5. This is summarized in Sect. 6.

   
2 The wind equations

The conservation laws of the general-relativistic magnetohydrodynamics have been derived by Bekenstein & Oron (1978). Camenzind (1986b) used these results to obtain the equations governing a stationary axisymmetric wind. In this section, we first recall these equations and we apply them to the particular case of a flow occurring in the equatorial plane. Then we derive a new formulation for this problem where the wind solution is a level contour of a Bernoulli-like function. A very similar formulation was earlier studied by Sakurai (1985) for the non-relativistic case. The similarity allows us to compare our results with those of the classical case.

2.1 Assumptions and basic equations

Any stationary and axisymmetric space-time can be represented by the following metric:

$${\rm d}s^2 = g_{tt} {\rm d}t^2 + 2 g_{t\phi} {\rm d}t {\rm ... ... + g_{\phi\phi} {\rm d}\phi^2 + g_{ab} {\rm d}x^a {\rm d}x^b$$ (1)

where the coordinates t and $\phi$ correspond to the two symmetries of the space-time defined by the two Killing fields $k=\partial_{\rm t}$ (stationarity) and $m=\partial_\phi$(axisymmetry). The metric coefficients g_ab depend only on the two remaining coordinates x^a (a=1,2). Note that in the whole paper we will make use of the (-+++) signature for the metric. The electromagnetic field is described by the field tensor  $F^{\mu\nu}$and the dual field tensor $F^{*\mu\nu} = \frac{1}{2}\epsilon^{\mu\nu\rho\sigma} F_{\rho\sigma}$( $\epsilon^{\mu\nu\rho\sigma}$ being the Levi-Civita alternating tensor) which satisfy the Maxwell equations

$$\nabla_\mu F^{*\mu\nu} = 0.$$ (2)

The motion of the plasma is governed by the energy- and momentum conservation equation

$$\nabla_\mu T^{\mu\nu} = 0,$$ (3)

where the energy-momentum tensor $T^{\mu\nu}$ is made up of the fluid part

$$T^{\mu\nu}_{\rm matter} = \rho \frac{{h}}{{c^2}} u^\mu u^\nu + P g^{\mu\nu}$$ (4)

and of the electromagnetic part

$$% T^{\mu\nu}_{\rm em} = \frac{1}{4\pi}\left( \frac{{b^2}}{... ...mu u^\nu + \frac{b^2}{2} g^{\mu\nu} - b^\nu b^\mu \right).$$ (5)

$b^\mu$ is the magnetic field according to a comoving observer, written as

$$b^\mu = F^{*\mu\nu} u_\nu$$ (6)

with $b^2=b_\mu b^\mu$. In the comoving frame $b^\mu$ reduces to the common magnetic field

$$b^\mu_{\rm co} = F^{*\mu0} = (0,\vec B).$$ (7)

All dissipative effects (heat conduction, viscosity, cooling by radiation, etc.) have been neglected so that the flow is adiabatic. Ideal MHD is also assumed, which means that the proper electric field as seen in the plasma frame vanishes

$$F^{\mu\nu} u_\mu = 0.$$ (8)

The following quantities appear in these equations: $u^\mu$ is the 4-velocity, $\rho$ is the comoving mass density, P is the pressure and h is the specific enthalpy, which is given by (assuming a constant adiabatic index $\gamma$)

$$h = {c^2} + \frac{\gamma}{\gamma-1}\frac{P}{\rho}\cdot$$ (9)

The last assumption is that the particle number (or mass) is conserved

$$\nabla_\mu \left(\rho u^\mu\right) = 0.$$ (10)

Before writing the wind equations, it is useful to define the specific angular momentum of the flow

$$l = -\frac{u_\phi}{u_t} {c^2}.$$ (11)

As a consequence of the symmetries, the non vanishing elements of the electromagnetic tensor are given by

$$- F_{1t} = {c} \Omega(P)\frac{\rho}{\eta(P)}\sqrt{-g} u^2,$$ F_t1 = (12)

$$- F_{2t} = -{c} \Omega(P)\frac{\rho}{\eta(P)}\sqrt{-g} u^1,$$ F_t2 = (13)

$$F_{\phi1} -F_{1\phi} = -c \frac{\rho}{\eta(P)}\sqrt{-g} u^2,$$ = (14)

$$F_{\phi2} -F_{1\phi} = -c \frac{\rho}{\eta(P)}\sqrt{-g} u^{1},$$ = (15)

$$F_{21} = c \frac{\rho}{\eta(P)}\sqrt{-g} \left(\Omega(P) {u^t} - {u}^\phi \right),$$ F₁₂ = (16)

where the angular frequency of the streamline $\Omega(P)$ at its footpoint and the mass flux per unit flux tube $\eta(P)$ are constant along each flow line P. The corresponding magnetic field is

$$\frac{\rho}{\eta(P)} \left[ 1+\frac{{u_t u^t}}{{c^2}} \left( 1-\frac{{\Omega(P)l}}{{c^2}} \right) \right] ,$$ b^t = (17)

$$b^\phi \frac{\rho}{\eta(P)} \left[ \Omega(P) + \frac{{u_t u}^\phi}{{c^2}} \left( 1-\frac{{\Omega(P) l}}{{c^2}} \right) \right],$$ = (18)

$$\frac{\rho}{\eta(P)} \frac{{u_t}}{{c^2}} \left(1-\frac{{\Omega(P)l}}{{c^2}}\right) u^1 ,$$ b¹ = (19)

$$\frac{\rho}{\eta(P)} \frac{{u_t}}{{c^2}} \left(1-\frac{{\Omega(P)l}}{{c^2}}\right) u^2.$$ b² = (20)

It is also useful to give the expression of the classical magnetic field in the frame of the central object for the comparison with the classical case. In this frame, it is related to the dual field tensor by $B_\mu=F_{*\mu0}$. Then

$$-\frac{{1}}{{c^2}}\frac{\rho}{\eta(P)} \left(\Omega(P)u^t-u^\phi\right) g_{t\phi},$$ B^t = (21)

$$B^\phi \frac{{1}}{{c^2}}\frac{\rho}{\eta(P)} \left(\Omega(P)u^t-u^\phi\right) g_{tt} ,$$ = (22)

$$-\frac{{1}}{{c^2}}\frac{\rho}{\eta(P)} \left( g_{tt} + \Omega(P) g_{t\phi} \right) u^1,$$ B¹ = (23)

$$-\frac{{1}}{{c^2}}\frac{\rho}{\eta(P)} \left( g_{tt} + \Omega(P) g_{t\phi} \right) u^2.$$ B² = (24)

Because the flow is stationary and axisymmetric, the total angular momentum L(P) and the total energy $E_{\rm tot}(P)$ are also conserved along each flow line P which provides us with two new equations

 

$$L(P) = -\frac{{u_t}}{{c^2}}\cdot \biggl\lbrace \frac{{h}}{{c^2}} ... ...^2}} l +\left( g_{t\phi} + \Omega(P) g_{\phi\phi} \right) \right] \biggr\rbrace$$ (25)

and

$$E_{\rm tot}(P) = -\frac{{u_t}}{{c^2}}\cdot\biggl\lbrace h + \Omeg... ...2}} l + \left( g_{t\phi} + \Omega(P) g_{\phi\phi}\right) \right]\biggr\rbrace .$$ (26)

It is convenient to write the total energy as $E_{\rm tot}(P) = {c^2}+E(P)+\Omega(P)L(P)$ so that the energy conservation can have the simpler form:

$$E(P) + {c^2} = -h \frac{{u_t}}{{c^2}} \left[1-\frac{{\Omega(P)l}}{{c^2}}\right].$$ (27)

Each flow line P is completely determined by the four constants $\Omega(P)$, $\eta(P)$, L(P) and E(P). The light surface is defined by

$$g_{tt} + 2 g_{t\phi}\Omega(P) + g_{\phi\phi}\Omega^2(P) = 0$$ (28)

and the Alfvén point is fixed by two conditions

$$\frac{{1}}{{c^2}} \left[ g_{tt} + 2 g_{t\phi}\Omega(P) + g_{\phi\phi}\Omega^2(P) \right]_{\rm A} = - M_{\rm A}^2,$$ (29)

$$\frac{{1}}{{c^2}} \frac{\left( g_{t\phi} + \Omega(P) g_{\phi\phi}\right)_{\rm A}} {M_{\rm A}^2} = \frac{L(P)}{E(P)+{c^2}},$$ (30)

where the "Mach'' number M is given by

$$M^2 = \frac{4\pi\eta^2(P)\,h}{\rho {c^2}}\cdot$$ (31)

One sees immediately from (28) and (29) that the Alfvén point stays always inside the light surface (because of  $M_{\rm A}>0$).

2.2 The wind equations in the equatorial plane

We use now the spherical coordinates (x¹=r and $x^2=\theta$) and limit our study to the equatorial plane $\theta=\frac{\pi}{2}$. The specification of the flow line P is dropped from here on and $\Omega$, $\eta$, E, and L are used instead of $\Omega(P)$, $\eta(P)$, E(P), and L(P). Because of the symmetry, $u^\theta$and $B^\theta$ vanish in this plane but this is not necessarily the case for their derivatives. Then the conservation of mass (10) can be written

$$\sqrt{-g} s(r) \rho u^r = \dot{m},$$ (32)

where the function s(r) depends on the geometry of the flux tubes. In the simple case where $\partial_\theta u^\theta = 0$ and $\partial_\theta B^\theta=0$, we have $s(r)={\rm const.}$ (constant opening angle). Otherwise we have

$$\left. \partial_\theta u^\theta \right\vert _{\theta=\frac{\pi}{2}} = \frac{s'(r)}{s(r)}u^r .$$ (33)

The conservation of angular momentum (25) and the conservation of energy (27) read

 

$$L = \left[ g_{\phi\phi}\frac{h}{c^2} - \frac{{g}_{{t}\phi}^{{2}} ... ...}} - {g_{tt} g}_{\phi\phi}}{{c^2}} \frac{\Phi^2\rho}{4\pi\dot{m}^2} \right] u^t$$ (34)

and

 

$$-h \left[\frac{{g_{tt}} +{\Omega g}_{{t}\phi}}{{c^2}} u^t + \frac{{g}_{{t}\phi} + {\Omega g}_{\phi\phi}}{{c^2}} u^\phi \right],$$ (35)

where instead of using $\eta$ we have introduced the magnetic flux $\Phi=\dot{m}/\eta$. The Eqs. (32), (34) and (35) are completed by the normalization of the four velocity

$$g_{tt} \left( u^t \right)^2 + 2 g_{t\phi} u^t u^\phi + g_{\... ...left(u^\phi\right)^2 + g_{\rm rr} \left(u^r\right)^2 = -c^2$$ (36)

and the equation of state. We assume here, like in Sakurai (1985), a polytropic relation $P=\kappa \rho^\gamma$ so that the specific enthalpy is given by

$$h = c^2 + \frac{\gamma}{\gamma-1}\kappa\rho^{\gamma-1}.$$ (37)

The system of Eqs. (32), (34)-(37) describes entirely the flow determined by the six constants $\Omega,E,L,\Phi,\dot{m},\kappa$ and the free function s(r). In addition the two supplementary conditions

 

$$-\frac{{1}}{{c^2}} \left( g_{tt} + 2\Omega g_{t\phi} + \Omega^2 g_{\phi\phi} \right)_{\rm A} M_{\rm A}^2$$ =

$$\frac{4\pi\dot{m}^2}{\Phi^2} \frac{h_{\rm A}}{\rho_{\rm A}c^2}$$ = (38)

and

$$- \frac{{1}}{{c^2}} \left( \frac{ g_{t\phi} + \Omega g_{\p... ...+ \Omega^2 g_{\phi\phi}} \right)_{\rm A} = \frac{L}{E+c^2}$$ (39)

(all quantities with index A are computed at the Alfvén point) must be fulfilled, so that the flow remains regular at the Alfvén point.

The classical limit (for a weak gravitational field and for velocities small compared to the speed of light) of this system of equations in the case where $s(r)={\rm const.}=1$ gives exactly the Eqs. (1) to (6) in Sakurai (1985). Notice that E, $\Phi$ and $\Omega$ have the same meaning in both papers whereas we use here different notations for the mass flux $\dot{m}$, the total angular momentum Land the polytropic constant $\kappa$ which are respectively f, $\Omega {r}_{\rm A}^{ 2}$ and K in Sakurai (1985). Notice also that the relation (39) between L and $r_{\rm A}$ tends towards $L=\Omega {r}_{\rm A}^{{2}}$ in the classical limit, so that all notations are fully consistent.

2.3 A Bernoulli-like formulation

Following Sakurai (1985) we use the dimensionless variables $x=r/r_{\rm A}$ and $y=\rho/\rho_{\rm A}$. The metric coefficients are also normalized to become dimensionless $\tilde g_{tt}(x)= {g_{tt}/c^2}$, $\tilde g_{t\phi}(x)= {g}_{{t}\phi}/({c r})$, $\tilde g_{\phi\phi}(x)= g_{\phi\phi}/r^2$, $\tilde{g}_{rr}(x)=g_{rr}$and $\sqrt{-\tilde{g}}(x)=\sqrt{-g}/r^2$. These coefficients may be not only functions of x but also of some parameters defining the metric (they are given for the Minkowski, Schwarzschild and Kerr metric in Appendix A). We define $\tilde{\varpi}^2 =\left({\tilde g}{^2}_{{t}\phi}-{\tilde g_{tt}\tilde g}_{\phi\phi}\right)/{c^2}$. We normalize s by $\tilde{s}(x)=s/s_{\rm A}$ and define a dimensionless specific enthalpy by $\tilde{h}=h/c^2$. Concerning four-vectors like $u^\mu$or $B^\mu$ we will use the definitions $\tilde{A}^t=A^t$, $\tilde{A}^\phi=r\sin{\theta} A^\phi$, $\tilde{A}^r=A^r$ and $\tilde{A}^\theta=r A^\theta$ so that the spatial part of these vectors is now given in the usual basis $\partial_r, \frac{1}{r}\partial_\theta, \frac{1}{r\sin{\theta}}\partial_\phi$, where all components have the same dimension. We introduce four normalized parameters

$$\beta' \frac{1}{c^2}\left( \frac{\dot{m}}{s_{\rm A}r_{\rm A}^2\rho_{\rm A}} \right)^2,$$ = (40)

$$\Theta' \frac{\gamma\kappa\rho_{\rm A}^{\gamma-1}}{c^2},$$ = (41)

$$\omega' \frac{\left(\Omega r_{\rm A}\right)^2}{c^2},$$ = (42)

$$\frac{E}{c^2},$$ E' = (43)

and we are now able to rewrite the system of Eqs. (32) and (34) to (37):

$$\sqrt{\beta'} = \sqrt{-\tilde{g}}\ \tilde{s}\ x^2 y \frac{\tilde{u}^r}{c}$$ (44)

 

$$\frac{K_{\rm A}}{M_{\rm A}^2}\left(E'+1\right) = x\left\lbrace \l... ...tilde{h}(1)}{M_{\rm A}^2} \tilde{\varpi}^2 x y\right] \tilde{u}^t \right\rbrace$$ (45)

 

$$E'+1 = - \tilde{h}\left[\left( \tilde g_{tt} + \sqrt{\omega'}\til... ...i} + \sqrt{\omega'}\tilde g_{\phi\phi} x \right)\frac{\tilde{u}^\phi}{c}\right]$$ (46)

 

$$-1 = \tilde g_{tt} \left(\tilde{u}^t\right)^2 + 2\tilde g_{t\phi}... ...tilde{u}^\phi}{c}\right)^2 + \tilde{g}_{rr}\left(\frac{\tilde{u}^r}{c}\right)^2$$ (47)

$$\tilde{h}(y) = 1+\frac{\Theta}{\gamma-1}y^{\gamma-1}.$$ (48)

In (46) and (45) the constants $L/r_{\rm A}c$ and $\Phi^2\rho_{\rm A}/4\pi\dot{m}^2$ have been eliminated using the conditions (38) and (39) at the Alfvén point, which now read

$$-\left( \tilde g_{tt}(1) + 2\sqrt{\omega'}\tilde g_{t\phi}(1) + \omega'\tilde g_{\phi\phi}(1) \right) M_{\rm A}^2$$ =

$$\frac{4\pi\dot{m}^2}{\Phi^2} \frac{\tilde{h}(1)}{\rho_{\rm A}}$$ = (49)

and

$$- \frac{\tilde g_{t\phi}(1) + \sqrt{\omega'}\tilde g_{\phi\phi}(1)} {M_{\rm A}^2} = \frac{L}{r_{\rm A} c}\frac{1}{E'+1}$$ (50)

so that

$$\frac{\Phi^2\rho_{\rm A}}{4\pi\dot{m}^2} = \frac{\tilde{h}(... ...{\rm A} c} = \frac{K_{\rm A}}{M_{\rm A}^2}\left(E'+1\right)$$ (51)

with $K_{\rm A} = \tilde g_{t\phi}(1) + \sqrt{\omega'}\tilde g_{\phi\phi}(1)$.

The component $\tilde{u}^{r}$ can be expressed from (44) and substituted into (47) to provide a first relation between $\tilde{u}^t$ and $\tilde{u}^\phi$. A second relation between these two components is given by (45) when subtracting $K_{\rm A}/M_{\rm A}^2\times$ (46). It allows us to express all components of the four velocity as functions of only xand y. Then the remaining Eq. (46) becomes a Bernoulli-like equation $\tilde{H}(x,y)={\rm const.}$ as (10) of Sakurai (1985). We do not further elaborate on the different steps which lead to this final expressions:

$$\tilde{h}(y) = 1+\frac{\Theta'}{\gamma-1}y^{\gamma-1},$$ (52)

$$\tilde{u}^t = \sqrt{K(x,y)}\frac{N(x,y)}{\sqrt{\mathcal{D}(x,y)}},$$ (53)

$$\frac{\tilde{u}^\phi}{c} = \sqrt{K(x,y)}\frac{D(x,y)}{\sqrt{\mathcal{D}(x,y)}},$$ (54)

$$\frac{\tilde{u}^r}{c} = \frac{\sqrt{\beta'}}{\sqrt{-\tilde{g}}\tilde{s}x^2y},$$ (55)

 

$$\tilde{H}(x,y)+1 \tilde{\varpi}^2 x \tilde{h}(y)\,\sqrt{K(x,y)} \frac{\left\vert\mathcal{N}(x,y)\right\vert}{\sqrt{\mathcal{D}(x,y)}}$$ =

= E'+1, (56)

where we have introduced the following auxiliary functions:

$$1 + \frac{\tilde{g}_{rr}}{-\tilde{g}} \frac{\beta'}{\tilde{s}^2 x^4 y^2},$$ K(x,y) = (57)

$$\left[ \left(M_{\rm A}^2+\sqrt{\omega'}K_{\rm A}\right) \tilde g_... ...\rm A} \tilde g_{t\phi} \right] \tilde{h}(y) -\tilde{h}(1)\tilde{\varpi}^2 x y,$$ N(x,y) = (58)

$$-\left[ \left(M_{\rm A}^2+\sqrt{\omega'}K_{\rm A}\right) \tilde g... ... g_{tt} \right] \tilde{h}(y) -\sqrt{\omega'}\tilde{h}(1)\tilde{\varpi}^2 x^2 y,$$ D(x,y) =

$$\mathcal{N}(x,y) \left( \tilde g_{tt} + 2\sqrt{\omega'}x\ \tilde g_{t\phi} + \omega' x^2\ \tilde g_{\phi\phi} \right) \tilde{h}(1) y + M_{\rm A}^2\tilde{h}(y),$$ = (59)

$$\mathcal{D}(x,y) -\left[ \tilde g_{\phi\phi} D^2(x,y) + 2\tilde g_{t\phi} D(x,y) N(x,y) +\tilde g_{tt} N^2(x,y)\right].$$ = (60)

Equation (56) is the Bernoulli equation we will now consider. The solution y(x) of the wind equations appears as the level contour E' of the surface H(x,y). Notice that the classical limits of $\beta'$, $\Theta '$ and $\omega '$ are not exactly the corresponding $\beta$, $\Theta$ and $\omega$ parameters used by Sakurai (1985), who made the choice of normalizing these quantities with the value of the gravitational potential at the Alfvén point $GM/r_{\rm A}$ (see Eqs. (11a), (11b) and (11c) of Sakurai 1985) whereas we used c² to make the definition of these parameters more general. However a simple relation applies between Sakurai's and our parameters: $\beta/\beta'=\Theta/\Theta'=\omega/\omega'=r_{\rm A}/r_{\rm g}$ where $r_{\rm g}=GM/c^2$ is the gravitational radius. This is also valid for the classical limit of the Bernoulli function $\tilde{H}(x,y)$ and the definition used by Sakurai (1985).

Before studying $\tilde{H}(x,y)$ in the following section, we have to note that this function is not defined everywhere in the region x>0, y>0 as it is the case in the classical limit. The function $\mathcal{D}(x,y)$ must be strictly positive so that $\tilde{H}(x,y)$is well defined and the velocities are not imaginary. The domain where this condition applies is determined in Appendix C. In the sub-Alfvénic region (y>1) this domain lies always inside the light surface, its location is given more precisely in Appendix B.

   
3 Description of the solutions

Figure 1: Solution plane of the wind equations. The gray region corresponds to the domain where the Bernoulli function is well defined. In the sub-Alfvénic region (y>1), it is limited by the light surface. The thick line (Alfvén mode) separates the sub- and the super-Alfvénic modes. The dashed line indicates the slow (y>1) and fast (y<1) mode Mach curves and the dotted line the gravitational throat curve. The gravitational throat curve is very close to the fast mode Mach curve for this particular case. The slow (S) and fast (F) critical points are the intersections of the Mach and throat curves. The Alfvén point is indicated by A.

Open with DEXTER

Figure 2: A few level contours of the Bernoulli function are shown. The physical solution (thick line) starts in the sub-Alfvénic region, crosses the slow mode Mach curve at the slow critical point (S) then reaches the Alfvén point (A) and enters the super-Alfvénic region where it crosses the fast mode Mach curve at the fast critical point (F). This calculation has been made for a Schwarzschild black hole with m=0.01 and the parameters ( $\gamma =4/3$, $\Theta '=0.04$ and $\omega '=0.5$) have been chosen so that the different points are well separated.

Open with DEXTER

3.1 Properties of the Bernoulli function

As described by Sakurai (1985) in the classical case, the Bernoulli function $\tilde{H}(x,y)$ has the following properties:

1.

At y=1 ( $\rho=\rho_{\rm A}$) $\tilde{H}(x,y)$ diverges if $x \ne 1$ and remains finite if x=1 ( $r=r_{\rm A}$). It means that all solutions going from the sub-Alfvénic region (y>1) to the super-Alfvénic region (y<1) must pass through the Alfvén point x=y=1 which is the only "hole'' in the infinite "wall'' y=1.

2.

Two important curves in the x-y plane are the so called slow/fast mode Mach curve defined by

$$\frac{\partial\tilde{H}}{\partial y}(x,y) = 0$$ (61)

(the slow mode corresponds to the sub-Alfvénic region y>1 and the fast mode to the super-Alfvénic region y<1) and the gravitational throat curve (so called by analogy with the de Laval nozzle) defined by

$$\frac{\partial\tilde{H}}{\partial x}(x,y) = 0.$$ (62)

At the intersections of these two curves, the function $\tilde{H}(x,y)$ is locally flat, corresponding to an X-type critical point (or O-type point).

3.

All level contours of $\tilde{H}(x,y)$ going from $y\to +\infty$ in the sub-Alfvénic region to $y\to 0^+$ in the super-Alfvénic region must cross these critical lines. They have to cross them simultaneously to be not interrupted. This means that the solution must pass through two critical points defined as the slow (respectively fast) critical point $x_{\rm s},y_{\rm s}$ (resp. $x_{\rm f},y_{\rm f}$), intersection of the slow (resp. fast) mode Mach curve and the gravitational throat curve. This imposes two new conditions for the solution:

$$\tilde{H}\left(x_{\rm s},y_{\rm s}\right) = E' \qquad{\rm and}\qquad \tilde{H}\left(x_{\rm f},y_{\rm f}\right) = E'.$$ (63)

Figure 1 shows the x-y plane for a particular choice of the parameters with the slow/fast mode Mach curve, the gravitational throat curve. Different level contours of the function $\tilde{H}(x,y)$ are shown in Fig. 2. The solution is one of the level contours and passes through the slow and fast critical points. These figures are very similar to Fig. 1 of Sakurai (1985) obtained in the classical case. Notice that Begelman & Li (1994) restrict their study to cold flows. A consequence of this assumption is that a purely radial flow is a singular case where the fast critical point is at infinite radius. The fast point moves inward to finite radii only if the flow diverges over-radially, like in the situations we will study in the next section. This is not a necessary condition in our more general model. Because we include gravity and thermal pressure the fast point is always located at finite distances even in a purely radial flow.

3.2 Classification of wind solutions by dimensionless parameters

Figure 3: Parameter space: there are no wind solutions in the lower left part ("Static'': the limits are only approximatively indicated). The upper left part ("Centrifugal'') corresponds to winds where the thermal pressure is negligible whereas in the lower right part ("Thermal'') the centrifugal force plays no essential role. For a Schwarzschild metric with m=0.01 (so $\omega '<0.98$ and $\omega <98$) we have computed several series of solutions with constant terminal Lorentz factor $\Gamma _{\infty }=1.01$, 1.5, 10, 100 and 1000 (solid lines). For $\Gamma _{\infty }=10$ we show also the case where m=0.1 (dashed line: $\omega <8$) and m=0.001 (dotted line: $\omega <998$). The dotted horizontal line corresponds to $\xi =1$ (initially the electromagnetic and the matter energy fluxes are equal) for m=0.01. For m=0.01 and $\Gamma _{\infty }=1.01$ and 1.5 we indicate at different positions the value of the efficiency $\mathit{eff}$ of the magnetic to kinetic energy conversion.

Open with DEXTER

Equation (56) can formally be written as

$$\tilde{H}\left(x,y\,; \beta',\Theta',\omega', m, g, \gamma\right) = E',$$ (64)

where m stands for the parameters defining the metric (see Appendix A for the definition of m in usual cases) and g stands for the parameters fixing the geometry of the flux tubes (i.e. defining $\tilde{s}(x)$). We will for the moment restrict our study to the case where $\tilde{s}={\rm const.}=1$ and we fix the adiabatic index to $\gamma =4/3$ for the rest of this paper. From the two conditions at the slow/fast critical points (63) we can fix two parameters. In practice, following Sakurai (1985), we adjust $\beta'$ so that the value of $\tilde{H}$ is the same at the two critical points and this value gives E'. Therefore all solutions are determined by only three independent parameters: $\Theta '$, $\omega '$ and m. Compared to the classical case studied by Sakurai (1985), there is one supplementary parameter: m. This is due to the presence of a characteristic length scale related to the structure of the space-time (typically the gravitational radius GM/c²) which has no classical counterpart.

The signification of these three parameters is clear: m measures the intensity of the gravitational field, $\Theta '$ gives the strength of the thermal pressure (the pressureless case which is often considered corresponds to $\Theta'=0$) and $\omega '$ measures the strength of the centrifugal force in accelerating the wind, or equivalently the effect of the magnetic field.

In the limit of the flat space-time (Minkowski metric) it is well known that the slow point does not exist anymore and that the solution cannot start at an arbitrary small radius because it cannot cross the slow mode Mach curve without having an infinite derivative . We therefore cannot describe the region near the source but this is clearly because gravity influences the solution significantly at small radii and the Minkowski approximation breaks down. In the case without gravity, as in the classical case, a supplementary parameter ($\beta'$) must be fixed (for instance by fixing the mass flux). In many cases the physical conditions at the basis of the wind are a complex and not well understood question, and some assumptions made everywhere else in the flow are probably not valid here (like the adiabaticity). In this context, just adopting a Minkowski metric and including in the $\beta'$ parameter all the unknown physics fixing the mass flux at the basis of the wind is more elegant than adopting a Schwarzschild (or Kerr) metric and applying the solution up to the source.

Before exploring the three-parameters space we have defined, it is useful to express the relation between m, $\omega '$ and $\Theta '$and more usual physical quantities. If m is the ratio of the gravitational radius of the source over the Alfvén radius, we have

$$m = 0.21\cdot \left(\frac{r_{\rm A}}{10^{6}\,{\rm cm}}\right)^{-1} \left(\frac{M}{1.4~M_\odot}\right).$$ (65)

The angular frequency $\Omega$ can be identified with the rotation rate of the source and fixes the value of $\omega '$:

$$\omega' = 0.11\cdot \left(\frac{\Omega}{10^4\,{\rm Hz}}\right)^2 \left(\frac{r_{\rm A}}{10^6\,{\rm rm}}\right)^2.$$ (66)

The local sound speed is defined by $(c^{\rm s}/c)^2=\gamma P/\rho h$ so that $\Theta '$ can be related to the ratio of the sound speed at the Alfvén point $c^{\rm s}_{\rm A}$ over the speed of light (for $\gamma =4/3$):

$$\Theta' = 1.0\times10^{-2}\cdot \left(\frac{c^{\rm s}_{\rm A... ...ac{1-3\left(c^{\rm s}_{\rm A}/c\right)^2}{0.97} \right)^{-1}.$$ (67)

The three parameters $\Theta '$, $\omega '$ and m being fixed, the solution of the wind equations is simply found by adjusting $\beta'$so that the condition (63) is fulfilled. For each step with $\beta'$ fixed, the slow and fast critical points are determined by a simple Newton-Raphson procedure. The exact expressions of $\tilde{H}$ and its derivatives are used. We have explored in detail the parameter space and the results are presented in Fig. 3. Notice that for a given m, $\omega '$ is limited to the interval

$$\omega'_{\rm min} \le \omega' \le \omega'_{\rm max},$$ (68)

where $\omega'_{\rm max}$ is due to the condition that the Alfvén point lies inside the light surface and $\omega'_{\rm min}$ is non vanishing only in the case of the Kerr metric. In this metric there are no solutions without rotation ($\omega'=0$) because the matter is forced to rotate in the vicinity of the central source. The analytical expressions of $\omega'_{\rm min}$ and $\omega'_{\rm max}$ are given in Appendix A. Here we have considered a Schwarzschild metric with different values of m. We show the results in $\Theta$-$\omega$ coordinates, where $\Theta=\Theta'/m$ and $\omega=\omega'/m$ are the parameters used by Sakurai (1985). Like in the classical case, there are no wind solutions in the lower left part of the plane ("Static'') because neither the centrifugal force (magnetic acceleration) nor the thermal pressure are sufficient to power the wind. The limits of this region are only approximatively indicated. For a pure thermal wind ($\omega=0$) and a Schwarzschild metric, the minimal value of $\Theta$ is given by

$$\Theta_{\rm min} = \left(\gamma-1\right)\frac{\frac{1}{\sqrt{1-2m}}-1}{m},$$ (69)

tending towards $\gamma-1$ for $m\to0$, in agreement with the classical case. In the pressureless case ($\Theta=0$) the minimal value of $\omega$ tends towards (3/2)^3/2 for $m\to0$, which also corresponds to the limit given by Sakurai (1985). The upper left part corresponds to winds where the centrifugal force dominates and the lower right part corresponds to pure thermal winds. In Fig. 3 we have plotted several solutions with constant terminal Lorentz factor $\Gamma _{\infty }=1.01$, 1.5, 10, 100 and 1000 for m=0.01 and in the particular case $\Gamma _{\infty }=10$ we have also plotted the same curves for m=0.1and m=0.01 to show the effect of varying the gravitational field.

   
4 Efficiency of the magnetic to kinetic energy transfer

4.1 Expressions of the energy fluxes

To study the efficiency of the winds computed in the previous section, we need to express the different components of the energy flux along the flow:

$$\dot{E}_{\rm matter} \dot{m}\left(-h \frac{{u_t}}{{c^2}}\right)$$ =

$$\dot{m}\left( c^2 + E + h \frac{\Omega}{c} \frac{u_\phi}{c} \right)$$ =

$$-\left( \tilde g_{tt} \tilde{u}^t + \tilde g_{t\phi}\frac{\tilde{u}^\phi}{c} \right)\ \tilde{h}(y)\ \dot{m} c^2,$$ = (70)

$$\dot{E}_{\rm em} \dot{m}\frac{\Phi^2\rho}{4\pi\dot{m}^2} \Omega\frac{{g}_{{t}\phi}^2-{g_{tt}g}_{\phi\phi}}{{c^2}} \left(\Omega u^t - u^\phi\right)$$ =

$$\dot{m}\left( \Omega L-h\frac{\Omega}{c}\frac{u_\phi}{c} \right)$$ =

$$\sqrt{\omega'}x\tilde{\varpi}^2 \frac{\tilde{h}(1)y}{M_{\rm A}^2}... ...( \sqrt{\omega'}x\tilde{u}^t-\frac{\tilde{u}^{\rm\phi}}{c} \right)\ \dot{m}c^2,$$ = (71)

$$\dot{E}_{\rm total} \dot{E}_{\rm matter} + \dot{E}_{\rm em}$$ =

$$\dot{m}\left(c^2 + E + \Omega L\right)$$ =

$$-\left( \tilde g_{tt}(1) + \sqrt{\omega}'\tilde g_{t\phi}(1) \right)\ \frac{E'+1}{M_{\rm A}^2}\ \dot{m}c^2 .$$ = (72)

Notice that the matter part is made up of the rest-mass, kinetic and internal energy of the matter. We define two parameters: the initial baryonic load $\eta$

$$\frac{1}{\eta} = \left. \frac{\dot{E}_{\rm matter}}{\dot{m}c^2}\right\vert _{x_0}$$ (73)

and the ratio of the initial power injected in the electromagnetic field over the initial power injected in the matter

$$\xi = \left. \frac{\dot{E}_{\rm em}}{\dot{E}_{\rm matter}} \right\vert _{x_{0}},$$ (74)

where x₀ is the radius where the wind starts. The value of $\xi$depends only weakly on x₀, and we have arbitrary chosen x₀=6m(or r₀=6r_g). One sees that $\eta$ will be fixed by $\Theta '$ [via the initial value of $\tilde{h}(y)$] whereas $\xi$ depends strongly on $\omega '$. Along the flow, the internal energy is converted into kinetic energy which accelerates the wind. Therefore if there were no magnetic field, the Lorentz factor at infinity would be $\Gamma_\infty=1/\eta$. However, the magnetic field can also contribute to the acceleration (when coupled with the rotation) and depending on the efficiency of the conversion of electromagnetic into kinetic energy, the terminal Lorentz factor can be larger than $1/\eta$, with a maximum value (complete conversion) given by

$$\Gamma_\infty^{\rm max} = \frac{1+\xi}{\eta}\cdot$$ (75)

In reality the conversion will never be complete and will be estimated by the following fraction

$$\mathit{eff} = 1 - \frac{\left.\dot{E}_{\rm em}\right\vert _... ...em}\right\vert _{0}} = \frac{\eta\Gamma_{\infty}-1}{\xi}\cdot$$ (76)

4.2 Inefficient conversion

In the pressureless case ($\Theta'=0$) it is well known that the magnetic to kinetic energy transfer is very inefficient for high terminal Lorentz factors (Michel 1969)

$$\frac{\Gamma_\infty \dot{m} c^2}{\dot{E}_{\rm total}} = \frac{1}{\Gamma_\infty^2}\cdot$$ (77)

The curves for constant terminal Lorentz factor in Fig. 3 show clearly that for highly relativistic winds ( $\Gamma _{\infty }=10$, 100 and 1000) the terminal Lorentz factor $\Gamma_\infty$ is independent of $\omega$ (or equivalently of $\xi$) which means that there can only be a tiny magnetic to kinetic energy conversion. When $\omega$ is very close to the maximal allowed value and the outflow is Poynting flux dominated (corresponding to the case where the Alfvén point is at the light surface radius) this tendency is not valid anymore. In this region the terminal Lorentz factor $\Gamma_\infty$ depends strongly on $\omega$ and is almost independent of $\Theta$ (or equivalently of $\eta$). However, even in this case only a tiny fraction of the magnetic energy is converted into kinetic energy. The converted energy amount is great compared to the initial energy in the matter part and therefore leads to a greater increase in $\Gamma$ and $\dot E_{\rm matter}$ throughout the flow. The efficiency $\mathit{eff}$ of the conversion is maximal in the pressureless case ($\Theta=0$) but is still rather small. In this case the parameters $\eta$ and $\xi$ are given by $\eta \simeq 1$ and $\xi \simeq \dot E_{\rm total}/\dot mc^2 - 1 = \Gamma^{3}_{\infty}-1$ which corresponds to

$$\mathit{eff} \simeq \frac{1}{1+\Gamma_\infty+\Gamma_\infty^2}\cdot$$ (78)

These tendencies are still present in mildly relativistic winds, as can be seen in Fig. 3 for $\Gamma_\infty=1.5$ where we have indicated the evolution of $\mathit{eff}$ along the curve. It is only within the classical limit that the conversion becomes important. This is shown in Fig. 3 (The case $\Gamma _{\infty }=1.01$ ( $v_\infty=0.14\,c$)). Here the efficiency of conversion reaches $\sim$60% in the pressureless case, which is in agreement with the classical study of Sakurai (1985).

4.3 Efficient conversion

All wind solutions showing (except in the classical limit) a very bad efficiency of the electromagnetic to kinetic energy conversion have been calculated for a particular geometry corresponding to $\tilde{s}={\rm const.}=1$ in our dimensionless units. This corresponds to magnetic flux tubes of constant opening angle. Under the assumption that the velocity is purely radial and constant at infinity it is possible to predict analytically the asymptotic behavior of the flow for any kind of geometry $\tilde{s}(x)$:

$$\tilde{u}^t \to \Gamma_\infty,$$

$$\tilde{u}^r \to \sqrt{\Gamma_\infty^2-1},$$

$$\tilde{u}^\phi \to 0,$$

$$y \simeq \frac{\sqrt{\beta'}}{\sqrt{\Gamma^2_{\infty}-1}} \frac{1}{\tilde{s}x^2},$$

so that the asymptotic expressions of the energy fluxes are

$$\frac{\dot{E}_{\rm matter}}{\dot{m}c^2} \to \Gamma_\infty,$$ (79)

$$\frac{\dot{E}_{\rm em}}{\dot{m} c^2} \to \frac{\omega'\sqrt... ...(1)}{M_{\rm A}^2} \frac{1}{v_\infty}\frac{1}{\tilde{s}} \cdot$$ (80)

From the last equation, one sees that the magnetic to kinetic energy conversion depends strongly on $\tilde{s}$. At infinity $\tilde{s}\to 0$ is unphysical because it would mean that the energy diverges. The case where the opening angle is constant at infinity corresponds to $\dot{E}_{\rm em}\to{\rm const.}> 0$ at infinity so that the conversion is not complete and the case where the opening angle diverges ( $\tilde{s}\to +\infty$) gives $\dot{E}_{\rm em}\to 0$ so that the conversion is complete. These results indicate that all models considered in the previous section are inefficient due to a particular choice of the geometry: $\tilde{s}={\rm const}$. This assumption is certainly correct at very large distance from the source but the opening angle may have variations at smaller radii. Equation (80) indicates that every region where the opening angle increases is a region of efficient magnetic to kinetic energy transfer. This is in agreement with the results of Begelman & Li (1994).

To check that the geometry is really the key parameter governing the efficiency of such winds we have computed some models using various laws for the evolution of the opening angle $\tilde{s}(x)$. The results are shown in Fig. 4 and confirm the previous analysis. We have considered a Schwarzschild metric with m=0.01 and a wind model characterized by $\omega'=0.97$ and $\Theta '=0.1$ (so the energy flux is initially dominated by the electromagnetic energy flux). We plot the different energy fluxes and the "Lorentz factor'' $\tilde{u}^{\rm t}$ as well as the geometrical function $\tilde{s}(x)$ we have used in each case. Figure 4a corresponds to the inefficient case $\tilde{s}={\rm const.}=1$. Figure 4b corresponds to the case where $\tilde{s}$increases in a region located between x₁=10 and x₂=18: the magnetic to kinetic energy conversion is immediately better. The efficiency $\mathit{eff}$ increases also in geometries with different shapes (Figs. 4d, f, g, h) and different locations of the $\tilde s>1$-region, provided that this region lies beyond the fast point as shown by Begelman & Li (1994). In this case $\tilde{s}_{\infty}$ is the only relevant quantity which governs $\mathit{eff}$. A $\tilde s>1$-region within the fast point like in Fig. 4e does not increase $\mathit{eff}$ and is similar to the purely radial case.

Figure 4: Effect of the geometry on the efficiency of the magnetic to kinetic energy conversion: we consider a Schwarzschild metric with m=0.01 and a wind solution characterized by $\Theta '=0.65$ and $\omega '=0.74$. This corresponds to an initial energy flux which is dominated by the electromagnetic part: $\xi =3.1$ for all models so that initially $75\%$ of the energy is magnetic. All solutions presented here have $\eta =2\times 10^{-2}$ except for cases e) ( $\eta =1.7\times 10^{-2}$) and h) ( $\eta =1.8\times 10^{-2}$). The slow and fast critical points are located at $x_{\rm s}=3.0\times 10^{-2}$ and $x_{\rm f}=2.0$ except for cases d) ( $x_{\rm f}=1.6$) and h) ( $x_{\rm f}=1.6$). On each figure ${\rm d}\log\tilde{s}/{\rm d}\log{x}$, $\tilde{s}(x)$ and the different components of the energy flux (matter/em and total) are presented as functions of the radius x. The "Lorentz factor'' $\tilde{u}^{\rm t}$ is also shown (dotted line). Three vertical dotted lines show the location of the slow (s), the Alfvén (A) and the fast point (f). Case a) $\tilde{s}={\rm const.}=1$. For this particular choice of the geometry, the conversion is extremely inefficient ( $\mathit{eff}=3.8\times 10^{-4}$) and the terminal Lorentz factor equals $\Gamma _{\infty }=50=1/\eta$. Case b) $\tilde{s}$ increases between x1=10 and x2=18 reaching a maximal slope ${\rm d}\log{\tilde{s}}/{\rm d}\log{x}=1$. The efficiency improves a lot: $\mathit{eff}=0.26$ and $\Gamma _{\infty }=90$. Case c) same as  b) but $\tilde{s}$ increases between x1=1000 and x2=1800 ( x2/x1 is the same). It changes neither the efficiency nor the terminal Lorentz factor. Case d) same as b) but $\tilde{s}$ increases between x1=1.5 and x2=2.7 ( x2/x1 is the same), i.e. before the position of the fast point in the reference solution a). Again the efficiency $\mathit{eff}=0.24$ and the terminal Lorentz factor $\Gamma _{\infty }=88$ are almost unchanged. Notice that the fast critical point has moved to be almost at x1. Case e) same as b) but $\tilde{s}$ increases between x1=0.1 and x2=0.18 ( x2/x1 is the same), i.e. before the Alfvén point. The efficiency is again very low: $\mathit{eff}=1.2\times 10^{-4}$ and $\Gamma _{\infty }=58\simeq 1/\eta$. Case f) same as b) but with a maximal slope of ${\rm d}\log{\tilde{s}}/{\rm d}\log{x}=4$. The efficiency is better: $\mathit{eff}=0.69$ and $\Gamma _{\infty }=157$. Case g) same as b) but the region where $\tilde{s}$ increases is larger: x1=10 and x2=100. Again the efficiency is better: $\mathit{eff}=0.68$ and $\Gamma _{\infty }=156$. Case h) we have considered a case where $\tilde{s}$ increases from x1=0.1 to x2=104 with a maximal slope ${\rm d}\log{\tilde{s}}/{\rm d}\log{x}=0.4$. Almost 90% of the magnetic energy is converted into kinetic energy ( $\mathit{eff}=0.88$) so that $\Gamma _{\infty }=206$.

Open with DEXTER

Figure 5: The influence of the geometry, gravitation and thermal energy upon the electromagnetic energy conversion. The solid, dotted and dashed lines correspond to different gravitational field strengths, m=0.001, m=0.01 and m=0.1. a) The case for $\omega '=0.95\,\omega '_{\rm max}$, $\Theta '=0.01$ which is a cold Poynting-flux dominated outflow. b) Thermal energy dominated fireball: $\omega '=0.1\,\omega '_{\rm max}$, $\Theta '=10$. The dotted line lies very close to the solid line and is not visible. c) Non-relativistic case where the rest mass dominates with $\omega '=0.1\,\omega '_{\rm max}$, $\Theta '=0.1$.

Open with DEXTER

Begelman & Li (1994) showed that the electromagnetic energy flux decreases like

$$\frac{\left.\dot{E}_{\rm em}\right\vert _{\infty}} {\left.... ...\vert _{x_0}} = \frac{\tilde{s}(x_{\rm f})}{\tilde{s}_\infty}$$ (81)

for a cold flow in Minkowski metric. If the asymptotic regime is already reached in the region where the opening angle increases, Eq. (80) shows that this relation should still be valid in the most general case, independent of the gravitational field or of the initial amount of thermal energy. To check the validity of this result we consider 9 different $\omega',\theta',m$ combinations, illustrating all possible situations and for each of them we compute the evolution of the efficiency when varying $\tilde{s}_{\infty}$. As Fig. 4 shows, the exact shape of the geometry is not important, so we adopt a particular choice where $\tilde{s}$ rises from $\tilde{s}_{0}=1$ to $\tilde{s}_{\infty}$ between x₁=100 and x₂=200. This region lies always in the super-Alfvénic region, which as discussed above is the condition for magnetic to kinetic energy conversion. Figure 5 shows the quantity

$$\frac{\left.\dot E_{\rm em}\right\vert _{\infty}} {\left.\... ...ht\vert _{x_0}} \cdot\frac{\tilde{s}_{\infty}}{\tilde{s}_{0}}$$ (82)

plotted over $\tilde{s}_{\infty}/\tilde{s}_{0}$ for the 9 different cases. Notice that with our choice of geometry $\tilde{s}(x_{\rm f})=\tilde{s}_{0}$. One sees that gravity and pressure changes the simple picture a bit. In the cold cases (Figs. 5a and c) the converted energy fraction decreases for a stronger gravitational field and high values of m. On the other hand the gravitational field increases the energy conversion by a small amount in the hot thermal dominated case as seen in Fig. 5b. But (81) remains still valid within a factor of 2. For the cases of low gravity and low thermal energy (solid lines in Figs. 5a and c) the quantity (82) approaches 1as expected.

We can therefore conclude that the flow geometry always dominates the energy conversion and all other parameters play an only minor role.

   
5 Application to gamma-ray bursts

Since the discovery of their optical afterglows gamma-ray bursts (hereafter GRBs) have been known to be located at cosmological distance. More than ten redshifts have been measured from z=0.43(GRB990712) to z=4.5 (GRB000131). The corresponding radiated energy in the gamma-ray domain (20-20000keV) ranges from $5\times10^{51}\,$erg (GRB970228) to $2\times 10^{54}\,$erg (GRB990123) assuming isotropic emission. Most sources that have been proposed to explain such a huge release of energy in a few seconds involve a rapidly rotating compact stellar-mass core. Among them the two most popular are mergers of compact objects (neutron stars binary or neutron star - black hole systems) or collapses of very massive stars to a black hole (collapsars) (Mészáros & Rees 1992; Narayan et al. 1992; Mochkovitch et al. 1993; Woosley 1993; Paczynski 1998). In both cases, the resulting system is a stellar mass black hole surrounded by a thick torus made of stellar debris or of infalling stellar material partially supported by centrifugal forces. An other interesting proposition (Usov 1992; Kluzniak & Ruderman 1998; Spruit 1999) associates GRBs with highly magnetized millisecond pulsars. The location of the detected optical counterparts, well inside their host galaxy and possibly associated with star-forming regions, seems to favor the collapsar scenario. However the other propositions cannot be ruled out, at least for short bursts, for which no optical counterpart has been detected yet.

Whatever the source is, the released energy must initially be injected in a wind which eventually becomes relativistic. The existence of such a relativistic wind has been directly inferred from the observations of radio scintillation in GRB970508 (Frail et al. 1997) and is also needed to avoid photon-photon annihilation. The absence of signature of this last process in the BATSE spectra of GRBs implies very high Lorentz factor for the wind: $\Gamma \sim 100$-1000(Goodman 1986; Baring 1995). The second step consists in the conversion of a fraction of the wind kinetic energy into gamma-rays, probably via the formation of shocks within the wind itself (Rees & Mészáros 1994; Daigne & Mochkovitch 1998). Such internal shocks are expected if the wind is generated with a highly non uniform distribution of the Lorentz factor so that rapid layers catch up with slower ones. In the last step, the wind is decelerated when it interacts with the environment of the source and the resulting external shock is responsible for the afterglow observed in X-ray, optical and radio bands.

The origin of the relativistic wind is the most complex of the three steps in this scenario. Several proposals have been made but only few calculations have been performed so that none appears to be fully conclusive. However it is suspected that large magnetic fields play an important role. In a previous paper (Spruit et al. 2001) we have considered different possible geometries of magnetic fields in GRB outflows and we have proposed that in many cases, dissipation of magnetic energy by reconnection should occur. The model we have presented in this paper allows us to investigate these questions in more details. In particular we focus on the case where the outflow generated by the central engine is initially Poynting flux dominated (in the following, we assume that only $10\%$ of the energy is initially injected in the matter). To be consistent with the observations showing that at the beginning of the afterglow emission, the matter flow is highly relativistic, we also impose that the terminal Lorentz factor has a large value (in the following, we will adopt $\Gamma_\infty=100$). This implies a reasonable efficiency of the magnetic to kinetic energy conversion. The goal of the study presented in this section is to illustrate that there are geometries allowing such an efficiency and to discuss the possibility of magnetic reconnection in this scenario.

Spruit et al. (2001) have shown that for typical GRB outflows the MHD approximation is valid to very large distance ( $10^{19}\,$cm) which is the main assumption of our calculations. The second main assumption - the stationarity of the flow - is of course less justified in the case of GRBs. However we can estimate the time scale to reach the stationary regime in our wind solutions as the time needed by a particle starting from the basis of a flow line to reach the Alfvén point:

$$t_{\rm stat} = \frac{r_{\rm a}}{c}\int_{x_{0}}^{1} \frac{\tilde{u}^{t}}{\tilde{u}^{r}}{\rm d}x$$ (83)

(in the source frame). Let us estimate this time scale in a particular case. We consider a Poynting flux dominated wind (we adopt $\xi=9.0$so that only 10% of the energy flux is initially injected in the matter) with a moderately low initial baryonic load (we take $\eta=1/50$). We impose that the terminal Lorentz factor is $\Gamma_\infty=100$. If there were no magnetic to kinetic energy conversion, the Lorentz factor at infinity would only be $1/\eta=50$. In order to get a final Lorentz factor of 100, we need to assume that the geometry allows an efficiency $\mathit{eff} = \left(\eta\Gamma_{\infty}-1\right)/\xi = 1/9$. We have shown in Sect. 4 that this implies $\tilde{s}_{\infty}/\tilde{s}_{0} \simeq \xi / \left(1+\xi-\eta\Gamma_{\infty}\right)=1.125$. For a given m, the two other parameters $\Theta '$ and $\omega '$ are fixed by the values of $\xi$ and $\eta$. We find that the following set of parameters: m=0.069, $\Theta'=3.8$ and $\omega'=0.78$ fulfill the requirements and corresponds to a reasonable value of the Alfvén radius $r_{\rm A}$ and the angular frequency $\Omega$ in the case of a millisecond pulsar-like source ( $M=1.4~M_\odot$) which is most likely leading to an equatorial flow as we are considering here: $r_{\rm A}=3.0\times 10^{6}\,$cm and $\Omega = 8.8\times 10^{3}\,$Hz. Figure 6 shows the evolution of the "Lorentz factor'' and the electromagnetic and matter energy fluxes in this case. We have assumed a simple geometry like those in Sect. 4 where $\tilde{s}$ increases in a region located between x₁=300 and x₂=900, well outside the fast critical point radius. The corresponding time scale to reach the stationary regime $t_{\rm stat}$ is between $\sim$ $r_{\rm A}/c$ and $\sim$ $2 r_{\rm A} /c$, depending on the adopted value of the initial radius x₀. As $r_{\rm A}/c = 10^{-4}\,$cms^-1 here, this is compatible with the timescale of the variability observed in GRBs profiles. This means that when the physical conditions at the basis of the flow vary on a time scale ms the flow reacts instantaneously to reach a new stationary state corresponding to the new boundary conditions. Thus our calculation is a good approximation for the relativistic wind of GRBs. If the wind produced by the source lasts for a duration $t_{\rm w}$, our solution is appropriate for the physical quantities within the corresponding shell when it is located at radius r.

On the solution we present in Fig. 6, the acceleration occurs in two phases. First the initial thermal energy is converted into kinetic energy, the magnetic energy remaining unchanged. This phase ends at $r\simeq 10^{9}\,$cm where $\Gamma_{\infty}\simeq 1/\eta \simeq 50$. The second phase occurs in the region where the opening angle increases. Here a magnetic to kinetic energy conversion takes place. We define the acceleration radius $r_{\rm acc}$ as the radius where the flow reaches a Lorentz factor of $\Gamma=0.95\Gamma_\infty$ and the acceleration can be considered as finished. The value of this radius is completely dominated by the unknown flow geometry and equals $r_{\rm acc} \simeq 2.7\times 10^{9}\,$cm in this case. Even if the location of the region where the opening angle diverges would extend to higher radii up to 10¹⁰-$10^{11}\,$cm, this radius is well below two other important radii: the photosphere radius $r_{\rm ph}$ where the wind becomes transparent and the reconnection radius $r_{\rm rec}$ where the reconnection of the magnetic field should occur. These two radii have been estimated in Spruit et al. (2001). The photosphere radius is the solution of

$$1 = \int_{r}^{r+2\Gamma^2 c t_{\rm w}} \frac{\kappa \rho}{2\Gamma^2}{\rm d}r$$ (84)

and is independent of the duration of the burst $t_{\rm w}$ as long as $r_{\rm ph} \ll 2\Gamma^2 c t_{\rm w} = 6\times 10^{14}\,{\rm cm}\cdot\left(\Gamma/100\right)^2 \left(t_{\rm w} / 1\,{\rm s}\right)$. Here we have $r_{\rm ph}=6.2\times 10^{10}\,$cm for $\dot{E}_{\rm total}=10^{51}\,$erg/s and $r_{\rm ph}=5.9\times 10^{11}\,$cm for $\dot{E}_{\rm total}=10^{52}\,$erg/s. This interval is marked by a thick line in Fig. 6. As the remaining thermal energy in the wind at such a large radius is very small, our adiabatic wind solution applies up to the reconnection radius, where magnetic dissipation starts. This radius is given by

$$r_{\rm rec} \simeq \frac{\pi c}{\epsilon \Omega}\Gamma^2 \left(1+\frac{1}{\xi}\right)^{1/2}$$ (85)

where $\epsilon<1$ is a numerical factor of order unity measuring the reconnection speed in unit of the Alfvén speed. In our case we have $r_{\rm rec}\simeq 1.1\times 10^{11}\,{\rm cm}\cdot\epsilon$. As the magnetic energy flux is still 80% of the total energy flux at $r_{\rm rec}$, a very large amount of energy can possibly be dissipated at this large distance. Depending on the value of $\dot{E}_{\rm total}$ and $\epsilon$, such reconnection events may start when the wind is still optically thick (low $\epsilon$, high $\dot{E}_{\rm total}$) or when the wind is already transparent (high $\epsilon$, low $\dot{E}_{\rm total}$). As the dissipated magnetic energy is probably first converted into thermal energy, the consequences for the wind may be very different in these two cases. (i) if the wind is optically thick, this injection of thermal energy should be converted, at least partially (up to the photosphere radius) into kinetic energy, leading to a third phase of acceleration; (ii) on the other hand, if the wind is transparent, reconnection events could directly contribute to the observed emission. Notice that all the radii we have computed are usually small compared to the typical radius where internal shocks occur (with $\Gamma_\infty=100$)

$$r_{\rm IS}\simeq 3\times 10^{14}\,{\rm cm}\cdot \left(\frac{t_{\rm var}}{1\,{\rm s}}\right),$$ (86)

where $t_{\rm var}$ is the typical time scale of the variability in the initial distribution of the Lorentz factor and also small compared to the deceleration radius where the external shock becomes efficient (with $\Gamma_\infty=100$)

$$r_{\rm dec}\simeq 5\times 10^{16}\,{\rm cm}\cdot \left(\fr... ...\right)^{1/3} \left(\frac{n}{1\,{\rm cm}^{-3}}\right)^{-1/3},$$ (87)

where n is the density of the external medium and E the total energy of the wind at this radius. So these two "standard'' mechanisms are not affected by the reconnection events. However the relevant energy flux will be the kinetic energy flux at $r_{\rm acc}$, possibly increased to a larger value if the reconnection starts in the optically thick regime.

	Figure 6: Geometry, "Lorentz factor'' and energy fluxes for our example. The vertical dotted lines mark the radii of the slow-, the Alfvén-, fast point and the acceleration radius.
Open with DEXTER

   
6 Conclusions

We have presented here a new formulation of the equations governing a stationary axisymmetric MHD flow in the equatorial plane. This formulation includes an exact treatment of all effects: thermal pressure, gravity and arbitrary shape of the magnetic flux tubes. The wind solution appears as the level contour of a Bernoulli-function which passes through two particular points: the slow and fast critical points. It allows a direct comparison with the classical model of Weber & Davis (1967), in particular in the formulation given by Sakurai (1985). Thus the specifically relativistic effects are easily identified.

We have used our model to extend the study of the magnetic to kinetic energy conversion made by Begelman & Li (1994). We show that the main parameter which fixes this efficiency is the shape of the magnetic flux tubes. In the case of a constant opening angle, non-relativistic flows have a good efficiency of the magnetic to kinetic energy conversion but as soon as the terminal Lorentz factor is greater than $\sim$1.5, this efficiency decreases rapidly. Such relativistic winds are not able to transfer a large fraction of their magnetic energy to the matter. On the other hand, regions where the opening angle diverges from the constant case are very efficient in converting magnetic into kinetic energy, even in the ultra-relativistic case. This is true as long as such regions are located beyond the fast critical point. Gravity and the thermal pressure play only a minor role.

In Sect. 5, we apply this model in the context of gamma-ray bursts (GRBs hereafter). In the case where the wind produced by the source of GRBs is initially Poynting flux dominated, we have shown that the efficiency of the acceleration strongly depends on the geometry of the magnetic flux tubes. We found that a large variety of situations is expected. If the magnetic tubes have the possibility to diverge strongly from a constant opening angle, it is possible that most of the energy is eventually in kinetic form. On the other hand it is very likely that the magnetic to kinetic energy conversion is incomplete and that the wind is still Poynting flux dominated when it has reached its terminal Lorentz factor. We have demonstrated on one example that such a wind can lead to very promising situations compared to the standard picture: a large amount of the magnetic energy can be dissipated at large radii by reconnection. This reconnection can start when the wind is optically thick or already transparent. So the large magnetic energy reservoir could have two effects: a supplementary acceleration phase increasing the final magnetic to kinetic energy conversion efficiency and/or a direct contribution to the emission. These two possibilities will be investigated in a future work.

Acknowledgements

The authors would like to thank Dr. H. C. Spruit for many stimulating discussions, important suggestions and reading the manuscript.

   
Appendix A: Metric coefficients

The following table gives the metric coefficients in normalized units in three cases: the Minkowski (M), Schwarzschild (S) and Kerr (K) space-times.

$$\begin{array}{\vert c\vert c\vert c\vert c\vert} \hline & {... ...sqrt{1-2m+a^2m^2}}{1+a^2m^2+2 a^2 m^3}\\ \hline \end{array}$$ (A.1)

where $m=GM/r_{\rm a}c^2$ and a=Jc/GM² (where J is the total angular momentum of the black hole and $0\le a\le1$). The radii $x_{\rm h}$ and $x_{\rm e}$ are respectively the radius of the horizon and of the ergosphere. We consider only the case where $0\le m \le \frac{1}{2}$ (the Alfvén point is outside the ergosphere). The minimum value $\omega'_{\rm min}$ of $\omega '$ corresponds to the condition $K_{\rm A} \ge 0$ (positive total angular momentum L) and the maximum value $\omega'_{\rm max}$ corresponds to the condition $M_{\rm A}^2>0$ (the Alfvén point must be inside the light surface).

   
Appendix B: The light surface

The light surface is defined by (28) which limits the region defined by

$$\xi^2(x) a^2 m^2\omega' - 1 + 2\frac{m}{x}\left(1-a m\sqrt{\omega'}\right)^2 + \omega' x^2$$ =

$$\le$$ 0. (B.2)

In the general case (with $\omega'\le\omega'_{\rm max}$, this corresponds to a domain $x_{e} \le x^{-}_{\rm lc} \le x \le x^{+}_{\rm lc}$ including the Alfvén point (x_a=1). In the Minkowski case, $x^{-}_{\rm lc}=0$ and $x^{+}_{\rm lc}=\frac{1}{\sqrt{\omega'}}$.

   
Appendix C: Domain of definition of the Bernoulli function

Figure C.1: Domain of definition of the Bernoulli function: we consider a Schwarzschild black hole and draw this domain for different values of m and $\omega '$ (variations of $\Theta '$ correspond only to a contraction or a dilatation of the domain along the y axis. We adopt here $\Theta '=1$). Cases a) and b): m=0.2<m*. There are three possible configurations depending on the value of $\omega '$. Case a) corresponds to $\omega '=0.3<\omega '_*$ and case b to $\omega '=0.5>\omega '_*$. The case $\omega '=\omega '_*$ looks very similar to case b) and is not plotted. We added the fast/slow mode Mach curves (dotted line) and the gravitational throat curve (dashed line) with the corresponding slow (S) and fast (F) critical point. The thick line is the solution passing through S, F and the Alfvén point A. Case c) m=m*=1/3. This is the critical case where it is still possible to find a solution (here $\omega '=0.3$ and $\Theta '=1$).

The Bernoulli function $H\left(x,y\right)$ is defined for ${\cal D}(x,y)>0$ which gives the following condition:

 

A(x) Y² - 2 B(x) Y + C(x) > 0 (C.1)

where $Y=\tilde{h}(1) y / \tilde{h}(y)$ is a function of y and $\Theta '$ only and A, B and C are functions of x, $\omega '$, a and m only (notice that it is completely independent of the function $\tilde{s}(x)$). The function Y(y) is strictly increasing from Y(0)=0 to $Y(+\infty)=+\infty$ with Y(1)=1. So we can focus to the pressureless case $\Theta'=0$ and Y=y, all other cases $\Theta'>0$ corresponding only to a contraction of the domain along the y-axis. The coefficients A, B and C are given by

$$A(x) = -\left(\tilde{\varpi}^2x\right)^2\xi^2(x)$$ (C.2)

$$B(x) = M_{\rm A}^2\left(\tilde{\varpi}^2 x\right)^2$$ (C.3)

$$\tilde{\varpi}^2\left\lbrace \tilde{\varpi}^2(1) \xi^2(x)\right.$$ C(x) = (C.4)

$$\left. - M_{\rm A}^2 \left( \tilde g_{tt}(1)\tilde g_{\phi\phi} x... ...)\tilde g_{t\phi} x + \tilde g_{\phi\phi}(1)\tilde g_{tt} \right)\right\rbrace.$$

The radius x being fixed, the Eq. (C.1) has 1 or 2 positive roots delimiting the domain where H(x,y) is well defined. Many configurations are possible for a Kerr metric and we do not specify them in details here. We discuss only the case of the Schwarzschild metric. Four different configurations are possible depending on m (with a critical case at m_*=1/3) and $\omega '$(with a critical case at $\omega'=27 m^2 \left(1-2m\right)^2$). This is illustrated in Fig. C.1.

References

 
• Baring, M. 1995, Ap&SS, 231, 169 In the text NASA ADS  
• Begelman, M. C., & Li, Z.-Y. 1994, ApJ, 426, 269 In the text NASA ADS  
• Bekenstein, J. D., & Oron, E. 1978, Phys. Rev. D, 18, 1809 In the text NASA ADS  
• Camenzind, M. 1986a, A&A, 156, 137 In the text NASA ADS  
• Camenzind, M. 1986b, A&A, 162, 32 In the text NASA ADS  
• Camenzind, M. 1987, A&A, 184, 341 In the text NASA ADS  
• Daigne, F., & Mochkovitch, R. 1998, MNRAS, 296, 275 In the text NASA ADS  
• Frail, D., Kulkarni, S., Nicastro, S., Feroci, M., & Taylor, G. 1997, Nature, 389, 261 In the text NASA ADS  
• Goldreich, P., & Julian, W. 1970, ApJ, 160, 971 In the text NASA ADS  
• Goodman, J. 1986, ApJ, 308, L47 In the text NASA ADS  
• Kennel, C. F., Fujimura, F. S., & Okamoto, I. 1983, J. Astrophys. Geophys. Fluid Dyn., 26, 147 In the text  
• Kluzniak, W., & Ruderman, M. 1998, ApJ, 505, L113 In the text NASA ADS  
• Mészáros, P., & Rees, M. 1992, MNRAS, 257, 29 In the text  
• Michel, F. C. 1969, ApJ, 158, 727 In the text NASA ADS  
• Mochkovitch, R., Hernanz, M., Isern, J., & Martin, X. 1993, Nature, 361, 236 In the text NASA ADS  
• Narayan, R., Paczynski, B., & Piran, T. 1992, ApJ, 395, L83 In the text NASA ADS  
• Okamoto, I. 1978, MNRAS, 185, 69 In the text NASA ADS  
• Paczynski, B. 1998, ApJ, 494, L45 In the text NASA ADS  
• Parker, E. N. 1958, ApJ, 128, 664 In the text NASA ADS  
• Rees, M., & Mészáros, P. 1994, ApJ, 430, L93 In the text NASA ADS  
• Sakurai, T. 1985, A&A, 152, 121 In the text NASA ADS  
• Spruit, H. C. 1999, A&A, 341, L1 In the text NASA ADS  
• Spruit, H. C., Daigne, F., & Drenkhahn, G. 2001, A&A, 369, 694 In the text NASA ADS  
• Usov, V. V. 1992, Nature, 357, 472 In the text NASA ADS  
• Weber, E. J., & Davis, L. J. 1967, ApJ, 148, 217 In the text NASA ADS  
• Woosley, S. E. 1993, ApJ, 405, 273 In the text NASA ADS