Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.stsci.edu/~tle/papers/invisic_paper2.pdf Äàòà èçìåíåíèÿ: Tue Jun 10 22:30:22 2008 Äàòà èíäåêñèðîâàíèÿ: Tue Oct 2 13:14:24 2012 Êîäèðîâêà: Windows-1251 Ïîèñêîâûå ñëîâà: m 106

The Astrophysical Journal, 632:476 - 498, 2005 October 10
# 2005. The American Astronomical Society. All rights reserved. Printed in U.S.A.

PARTICLE ACCELERATION AND THE PRODUCTION OF RELATIVISTIC OUTFLOWS IN ADVECTION-DOMINATED ACCRETION DISKS WITH SHOCKS
Truong Le
E. O. Hulburt Center for Space Research, Naval Research Laboratory, Washington, DC 20375; truong.le@nrl.navy.mil

and Peter A. Becker
1

Center for Earth Observing and Space Research, George Mason University, Fairfax, VA 22030-4444; pbecker@gmu.edu Received 2005 March 9; accepted 2005 June 22

ABSTRACT Relativistic outflows ( jets) of matter are commonly observed from systems containing black holes. The strongest outflows occur in the radio-loud systems, in which the accretion disk is likely to have an advection-dominated structure. In these systems, it is clear that the binding energy of the accreting gas is emitted primarily in the form of particles rather than radiation. However, no comprehensive model for the disk structure and the associated outflows has yet been produced. In particular, none of the existing models establish a direct physical connection between the presence of the outflows and the action of a microphysical acceleration mechanism operating in the disk. In this paper we explore the possibility that the relativistic protons powering the jet are accelerated at a standing, centrifugally supported shock in the underlying accretion disk via the first-order Fermi mechanism. The theoretical analysis employed here parallels the early studies of cosmic-ray acceleration in supernova shock waves, and the particle acceleration and disk structure are treated in a coupled, self-consistent manner based on a rigorous mathematical approach. We find that first-order Fermi acceleration at standing shocks in advection-dominated disks proves to be a very efficient means for accelerating the jet particles. Using physical parameters appropriate for M87 and Sgr AÃ , we verify that the jet kinetic luminosities computed using our model agree with estimates based on observations of the sources. Subject headingg accretion, accretion disks -- black hole physics -- galaxies: jets -- hydrodynamics s:

1. INTRODUCTION A large body of observational evidence has established that extragalactic relativistic jets are commonly associated with radioloud active galactic nuclei (AGNs), which may contain hot, advection-dominated accretion disks. However, the precise nature of the mechanism responsible for transferring the gravitational potential energy from the infalling matter to the small population of nonthermal particles that escape to form the jet is not yet clear (see, e.g., Livio 1999). The most promising jet acceleration scenarios proposed so far are the Blandford-Znajek mechanism ( Blandford & Znajek 1977) and the electromagnetic cocoon model ( Lovelace 1976; Blandford & Payne 1982), which involve the extraction of energy from the rotation of the black hole or the accretion disk in order to power the outflow. While conceptually attractive, one finds that the complex physics involved in these models tends to obscure the nature of the fundamental microphysical processes. In particular, the introduction of the relativistic particles that escape to form the jet is usually made in an ad hoc manner without any reference to the dynamics of the associated accretion disk, although recent magnetohydrodynamic simulations carried out by De Villiers et al. (2005) and McKinney & Gammie (2004) have achieved a higher level of self-consistency. Given the relative complexity of the electromagnetic models, it is natural to ask whether the outflows can be explained in terms of well-understood microphysical processes operating in the hot, tenuous disk, such as the possible acceleration of the jet particles at a standing accretion shock.
Also at: Department of Physics and Astronomy, George Mason University, Fairfax, VA 22030- 4444.
1

It has been known for some time that inviscid accretion disks can display both shocked and shock-free (i.e., smooth) solutions depending on the values of the energy and angular momentum per unit mass in the gas supplied at a large radius (e.g., Chakrabarti 1989a; Chakrabarti & Molteni 1993; Kafatos & Yang 1994; Lu & Yuan 1997; Das et al. 2001). Shocks can also exist in viscous disks if the viscosity is relatively low (Chakrabarti 1996; Lu et al. 1999), although smooth solutions are always possible for the same set of upstream parameters ( Narayan et al. 1997; Chen et al. 1997). Hawley et al. (1984a, 1984b) have shown through general relativistic simulations that if the gas is falling with some rotation, then the centrifugal force can act as a ``wall,'' triggering the formation of a shock. Furthermore, the possibility that shock instabilities may generate the quasi-periodic oscillations (QPOs) observed in some sources containing black holes has been pointed out by Chakrabarti et al. (2004), Lanzafame et al. (1998), Molteni et al. (1996), and Chakrabarti & Molteni (1995). Nevertheless, shocks are ``optional'' even when they are allowed , and one is always free to construct models that avoid them. However, in general the shock solution possesses a higher entropy content than the shock-free solution, and therefore we argue based on the second law of thermodynamics that when possible, the shocked solution represents the preferred mode of accretion ( Becker & Kazanas 2001; Chakrabarti & Molteni 1993). Our primary objective in this paper is to explore the consequences of the presence of a shock in an advection-dominated accretion flow (ADAF ) disk for the acceleration of the nonthermal particles in the observed jets. The question of whether viscosity needs to be included in the disk model is difficult to answer in general. Several authors have shown that shock solutions are possible in viscous (e.g., Chakrabarti 1990, 1996; Lu et al. 476


OUTFLOWS FROM ADVECTION-DOMINATED DISKS 1999; Chakrabarti & Das 2004) as well as inviscid disks (e.g., Chakrabarti 1989a, 1989b, 1996; Abramowicz & Chakrabarti 1990; Yang & Kafatos 1995; Das et al. 2001). In particular, Chakrabarti (1990) and Chakrabarti & Das (2004) demonstrated that shocks can exist in viscous disks if the angular momentum and the viscosity are relatively low. Since the acceleration of particles in shocked disks has never been investigated before, in this first study we focus on inviscid flows containing isothermal shocks (e.g., Chakrabarti 1989a; Kafatos & Yang 1994; Lu & Yuan 1997), while deferring the treatment of viscous disks to future work. However, it is clearly important to address the potential connection between this idealized, inviscid calculation and the physical properties of real accretion disks, which undoubtedly have nonzero viscosity. We argue that the results presented here should be qualitatively similar to those obtained in a viscous disk, provided a shock is present, in which case efficient firstorder Fermi acceleration is expected to occur. While the possible existence of standing shocks in viscous disks is a controversial issue at the present time, we believe that the work of Chakrabarti (1990, 1996), Lu et al. (1999), and Chakrabarti & Das (2004) provides sufficient support for the possibility to motivate the present investigation. Although the effect of a standing shock in heating the gas in the postshock region has been examined by a number of previous authors for both viscid (Chakrabarti & Das 2004; Lu et al. 1999; Chakrabarti 1990) and inviscid (e.g., Lu & Yuan 1997, 1998; Yang & Kafatos 1995; Abramowicz & Chakrabarti 1990) disks, the implications of the shock for the acceleration of nonthermal particles in the disk have not been considered in detail before. However, a great deal of attention has been focused on particle acceleration in the vicinity of supernova-driven shock waves as a possible explanation for the observed cosmic-ray energy spectrum ( Blandford & Ostriker 1978; Jones & Ellison 1991). In the present paper we consider the analogous process occurring in hot ADAFs around black holes. These disks are ideal sites for first-order Fermi acceleration at shocks, because the plasma is collisionless and therefore a small fraction of the particles can gain a great deal of energy by repeatedly crossing the shock. Shock acceleration in the disk therefore provides an intriguing possible explanation for the powerful outflows of relativistic particles observed in many radio-loud systems ( Le & Becker 2004). The dynamical model for the disk /shock /outflow employed here was discussed by Le & Becker (2004), who demonstrated that the predicted kinetic power in the jets agrees with the observational estimates for M87 and Sgr AÃ . In this paper we present a more detailed development of the dynamical model, including a careful examination of the implications of the shock acceleration process for the evolution of the relativistic particle distribution in the disk and the jet. The number and energy densities of the relativistic particles are determined, along with the hydrodynamic structure of the disk, in a self-consistent manner by solving the fluid dynamical conservation equations and the transport equation simultaneously using a rigorous mathematical approach. In this sense, the model presented here represents a new type of synthesis between studies of accretion dynamics and particle transport. The remainder of the paper is organized as follows. In x 2we discuss the ADAF model assumptions and the possibility of shock acceleration in ADAF disks, and the general structure of the disk /shock model is examined in x 3. The isothermal shock jump conditions and the asymptotic variations of the physical parameters at both large and small radii are discussed in x 4. In x 5 we analyze the steady state transport equation governing the

477

distribution of the relativistic particles in the disk and the jet. Solutions for the number and energy density distributions of the relativistic particles are obtained in x 6, and detailed applications to the disks /outflows in M87 and Sgr AÃ are presented in x 7. The astrophysical implications of our results are discussed in x 8. 2. MODEL BACKGROUND Accretion onto a black hole involves differentially rotating flows in which the viscosity plays an essential role in transporting angular momentum outward , thereby allowing the accreting gas to spiral in toward the central mass ( Pringle 1981). In the ADAF model, it is assumed that the mass accretion rate is much smaller than the Eddington rate, ME cÀ2  À1 LE ? 2:2 ; 10À9 
À1



M M



M yrÀ1;

ð 1Þ

where the radiative efficiency parameter  P 10% and the Eddington luminosity is defined by LE 4GMmp c/T for pure, fully ionized hydrogen, with T , M, mp ,and c denoting the Thomson cross section, the black hole mass, the proton mass, and the speed of light, respectively. Due to the sub-Eddington accretion rates in these systems, the plasma is rather tenuous, and this strongly inhibits the efficiency of two-body radiative processes such as free-free emission. The gas is therefore unable to cool effectively within an accretion time, and consequently the gravitational potential energy dissipated by viscosity is stored in the gas as thermal energy instead of being radiated away (e.g., Narayan et al. 1997). The low density also reduces the level of Coulomb coupling between the ions and the electrons, resulting in a two-temperature configuration with the ion temperature (Ti $ 1012 K) close to the virial value and a much lower electron temperature (Te $ 109 K). In this scenario, most of the energy is advected across the horizon into the black hole, and the resulting X-ray luminosity is far below the Eddington value ( Becker & Le 2003; Becker & Subramanian 2005) with  T1. When the ion temperature is close to the virial temperature, as in ADAFs, the disk is gravitationally unbound (e.g., Narayan et al. 1997; Blandford & Begelman 1999; Becker et al. 2001). It follows that the original ADAF model was not entirely selfconsistent, since it neglected outflows. This motivated Blandford & Begelman (1999) to propose the self-similar advection-dominated inflow-outflow solution (ADIOS) to address the question of selfconsistency by including the possibility of powerful winds that carry away mass, energy, and angular momentum. In this Newtonian, nonrelativistic model, the dynamical solutions are not applicable near the event horizon, and therefore the ADIOS approach cannot be used to obtain a global understanding of the disk structure. This led Becker et al. (2001) to modify the ADIOS scenario to include general relativistic effects by replacing the Newtonian ? potential with the pseudo-Newtonian form ( Paczynski & Wiita 1980) È(r) ÀGM ; r À rS ð 2Þ

where rS ? 2GM /c 2 is the Schwarzschild radius for a black hole of mass M. This modified model is known as the self-similar relativistic advection-dominated inflow-outflow solution ( RADIOS). Despite the success of the self-similar RADIOS model in describing the general features of the disk /outflow structure, it does not


478

LE & BECKER

Vol. 632

provide a comprehensive picture, since no explicit microphysical acceleration mechanism is included. It is therefore natural to explore possible extensions to the ADAF scenario that incorporate a concrete acceleration mechanism capable of powering the outflows. The idea of shock acceleration in the environment of AGNs was first suggested by Blandford & Ostriker (1978). Subsequently, Protheroe & Kazanas (1983) and Kazanas & Ellison (1986) investigated shocks in spherically symmetric accretion flows as a possible explanation for the energetic radiation emitted by many AGNs. However, in these papers the acceleration of the particles was studied without the benefit of a detailed transport equation, and the assumption of spherical symmetry precludes the treatment of acceleration in disks. The state of the theory was advanced by Webb & Bogdan (1987) and Spruit (1987), who employed a transport equation to solve for the distribution of energetic particles in a spherical accretion flow characterized by a self-similar velocity profile terminating at a standing shock. While more quantitative in nature than the earlier models, these solutions are not applicable to disks, since the geometry is spherical and the velocity distribution is inappropriate. Hence, none of these previous models can be used to develop a single, global, self-consistent picture for the acceleration of relativistic particles in an accretion disk containing a shock. The success of the diffusive ( first-order Fermi) shock acceleration model in the cosmic-ray context suggests that the same mechanism may be responsible for powering the outflows commonly observed in radio-loud systems containing black holes. As a preliminary step in evaluating the potential relevance of shock acceleration as a possible explanation for the observed outflows, we need to consider the basic physical properties of the hot plasma in ADAF disks. One of the critical issues for determining the efficiency of shock acceleration in accretion disks is the role of particle-particle collisions in thermalizing the high-energy ions. The mean free path for ion-ion collisions is given in cgs units by (Subramanian et al. 1996) kii ? 1:8 ; 105 T2 i ; Ni ln à ð 3Þ

Fig. 1.--Schematic representation of our disk /shock /outflow model. The filled circles in the disk represent the test particles, and the unfilled circles represent the MHD scattering centers moving with the background gas. The test particles are injected at the shock location.

that can be achieved by the particles. First, at very high energies the particles will tend to lose energy to the waves due to recoil. Second, the mean free path k mag will eventually exceed the thickness of the disk as the particle energy is increased, resulting in escape from the disk without further acceleration. 3. TRANSONIC FLOW STRUCTURE As discussed in x 1, various authors have established that shocks can exist in both viscid and inviscid disks. In this first study of particle acceleration in shocked disks, we focus on the inviscid case, because it is the most straightforward to analyze from a mathematical viewpoint, and also because it serves to illustrate the basic physical principles involved. Moreover, we expect that the results obtained in the viscous case will be qualitatively similar to those presented here, since efficient Fermi acceleration will occur whether or not viscosity is present, provided the flow contains a shock. The equations governing the disk structure can yield solutions that include three possible types of standing shocks, namely, (1) Rankine-Hugoniot shocks, where the effective cooling processes are so inefficient that no energy is lost from the surface of the disk; (2) isentropic shocks, where the entropy generated at the shock is comparable to the amount radiated away; and (3) isothermal shocks, where the cooling processes are so efficient that the postshock sound speed and disk thickness remain the same as the preshock values. In the isothermal case, the shock must radiate away both energy and entropy through the upper and lower surfaces of the disk (e.g., Chakrabarti 1989a, 1989b; Abramowicz & Chakrabarti 1990). This renders the isothermal shock model particularly useful from the point of view of modeling outflows, since the energy lost from the shock can be identified with that powering the jet. On the other hand, RankineHugoniot shocks cannot be used if we are interested in any kind of escape. The isentropic shock is an intermediate case. In this paper, we focus exclusively on the isothermal shock model, since this case provides the strongest potential connection with the observed outflows. The model considered here is depicted schematically in Figure 1. In this scenario, the gas is accelerated gravitationally toward the central mass and experiences a shock transition due to an obstruction near the event horizon. The obstruction is provided by the ``centrifugal barrier,'' which is located between the inner and outer sonic points. Particles from the high-energy tail of the background Maxwellian distribution are accelerated at the shock discontinuity via the first-order Fermi mechanism, resulting in the formation of a nonthermal, relativistic particle

where Ni and Ti denote the thermal ion number density and temperature, respectively, and ln à is the Coulomb logarithm. In ADAF disks, kii greatly exceeds the vertical thickness of the disk, and therefore the shock and the flow in general are collisionless. However, the mean free path k mag for collisions between ions and magnetohydrodynamic ( MHD) waves is much shorter than kii for the thermal particles, and it is much longer than kii for the relativistic particles ( Ellison & Eichler 1984; Subramanian et al. 1996). The increase in k mag with increasing particle energy reflects the fact that the high-energy particles will interact only with the highest energy MHD waves. The low-energy background particles therefore tend to thermalize the energy they gain in crossing the shock due to collisions with magnetic waves. Conversely, the relativistic particles are able to diffuse back and forth across the shock many times, gaining a great deal of energy while avoiding thermalization due to the longer mean free path. The probability of multiple shock crossings decreases exponentially with the number of crossings, and the mean energy of the particles increases exponentially with the number of crossings. This combination of factors naturally gives rise to a powerlaw energy distribution, which is a general characteristic of Fermi processes ( Fermi 1954). Two effects limit the maximum energy


No. 1, 2005

OUTFLOWS FROM ADVECTION-DOMINATED DISKS

479

distribution in the disk. The spatial transport of the energetic particles within the disk is a stochastic process based on a threedimensional random walk through the MHD scattering centers in the accreting background gas. Consequently, some of the accelerated particles diffuse to the disk surface and become unbound, escaping through the upper and lower edges of the cylindrical shock to form the outflow, while others diffuse outward radially through the disk or advect across the event horizon into the black hole. In order to analyze the connection between the disk /shock model and the transport /acceleration of the relativistic particles, we consider the set of physical conservation equations employed by Chakrabarti (1989a) and Abramowicz & Chakrabarti (1990), who investigated the structure of a one-dimensional, steady state, axisymmetric, inviscid accretion flow based on the vertically averaged conservation equations. The effects of general relativity are incorporated in an approximate manner by using the pseudoNewtonian form for the gravitational potential per unit mass given by equation (2). The use of such a potential allows one to investigate the complicated physical processes taking place in the accretion disk within the context of a semiclassical framework while maintaining good agreement with fully relativistic calculations (see, e.g., Narayan et al. 1997; Becker & Subramanian 2005). The pseudo-Newtonian potential correctly reproduces the radius of the event horizon, the marginally bound orbit, and the ? marginally stable orbit ( Paczynski & Wiita 1980). Furthermore, the dynamics of freely falling particles near the event horizon computed using this potential agrees perfectly with the results obtained using the Schwarzschild metric, although time dilation is not included ( Becker & Le 2003). 3.1. Transport Rates Becker & Le (2003) and Becker & Subramanian (2005) demonstrated that three integrals of the flow are conserved in viscous ADAF disks, namely, the mass transport rate M ? 4rH v; the angular momentum transport rate J ? Mr 2 À G; and the energy transport rate ? ÀG þ M 1 v 2 þ 1 v 2 þ P þ U þ È ; E 2 2 ð 6Þ ð 5Þ ð 4Þ

where is the kinematic viscosity. The disk half-thickness H is given by the standard hydrostatic prescription H (r ) ? a ; K ð 8Þ

where a represents the adiabatic sound speed, defined by 1=2 P ; ð 9Þ a( r ) and K denotes the Keplerian angular velocity of matter in a circular orbit at radius r in the pseudo-Newtonian potential (eq. [2]), defined by 2 K GM 1 dÈ : ? 2 r dr r (r À rS ) ð10Þ

The quantities M and J are constant throughout the flow, and therefore they represent the rates at which mass and angular momentum, respectively, enter the black hole. The energy trans port rate E generally remains constant, although it will jump at the location of an isothermal shock if one is present in the disk. We can eliminate the torque G between equations (5) and (6) and combine the result with equation (9) to express the energy transport per unit mass as E 1 1 l 2 l0 l a2 þ È; ? v2 À 2 þ 2 þ 2 2r r À1 M ð11Þ

where l (r) r 2 (r) and l0 J /M denote the specific angular momentum at radius r and the (constant) angular momentum transport per unit mass, respectively. Note that equation (11) is valid for both viscid and inviscid flows. 3.2. Inviscid Flow Equations In the present application, viscosity is neglected, and therefore G ? 0 and the specific angular momentum is given by l (r) ? l0 ? constant ð12Þ

throughout the disk. It follows that the flow is purely adiabatic, except at the location of a possible isothermal shock ( Becker & Le 2003). In the inviscid case, equation (11) reduces to ? 1 2 1 l2 a2 vþ þ È: þ 2 2 2r À1 ð13Þ

where is the mass density, v is the radial velocity (defined to be positive for inflow), is the angular velocity, G is the torque, H is the disk half-thickness, v ? r is the azimuthal velocity, U is the internal energy density, and P ? ( À 1)U is the pressure. Each of the various quantities represents a vertical average over the disk structure. We also assume that the ratio of specific heats maintains a constant value throughout the flow. Note that all of the transport rates M , J , and E are defined to be positive for inflow. The torque G is related to the gradient of the angular velocity via the usual expression (e.g., Frank et al. 2002) G ? À4r 3 H d ; dr ð 7Þ

The resulting disk /shock model depends on three free parameters, namely, the energy transport per unit mass , the specific heat ratio , and the specific angular momentum l.The valueof will jump at the location of an isothermal shock if one exists in the disk, but the value of l remains constant throughout the flow. This implies that the specific angular momentum of the particles escaping through the upper and lower surfaces of the cylindrical shock must be equal to the average value of the specific angular momentum for the particles remaining in the disk, and therefore the outflow exerts no torque on the disk ( Becker et al. 2001). Since the flow is purely adiabatic in the absence of viscosity, the pressure and density variations are coupled according to P ? D0 ; ð14Þ


480

LE & BECKER 3.4. Critical Point Analysis

Vol. 632

where D0 is a parameter related to the specific entropy that remains constant except at the location of the isothermal shock, if one is present. By combining equations (4), (8), (9), (10), and (14), we find that the quantity Kr
3=2

( r À rS ) v a

( þ1)=( À1)

ð15Þ

is conserved throughout an adiabatic disk, except at the location of an isothermal shock. Following Becker & Le (2003), we refer to K as the ``entropy parameter,'' and we note that the entropy per particle S is related to K via S ? k ln K þ c0 ; ð16Þ

Equations (21) and (22) can be solved simultaneously to express vc and ac as explicit functions of the critical radius rc , which yields ! GMr 3 À l 2 (rc À rS )2 2 c vc ? 2 ; ð23Þ 2 (5rc À 3rS )(rc À rS ) rc ! GMr 3 À l 2 (rc À rS )2 2 c ac ? ( þ 1) : ð24Þ 2 (5rc À 3rS )(rc À rS ) rc The corresponding value of the entropy parameter K at the critical point is given by (see eq. [15]) Kc ? r
3=2 c

where k is the Boltzmann constant and c0 is a constant that depends on the composition of the gas but is independent of its state. 3.3. Critical Conditions By combining equations (13) and (15), one can solve for the flow velocity v as a function of r using a simple root-finding procedure. However, in order to understand the implications of the transonic (critical) nature of the accretion flow, we must also analyze the properties of the ``wind equation,'' which is the firstorder differential equation governing the flow velocity v . By differentiating equation (13) with respect to r, we obtain the steady state radial momentum equation dv l 2 dÈ 2a da À : ð17Þ v ? 3À dr r dr À 1 dr The derivative of the sound speed appearing on the right-hand side of this expression can be evaluated by using equations (10) and (15) to write 1 da À1 1 d K 1 d v 1 ? À À : a dr þ1 K dr v dr r ð18Þ

(rc À rS ) vc a

( þ1)=( À1) c

:

ð25Þ

By using equations (23) and (24) to substitute for v and a in equation (13), we can express the energy transport parameter in terms of rc , l, and , obtaining ! 1 l2 GM 2 GMr 3 À l 2 (rc À rS )2 c ? À þ : ð26Þ 2 2 2 rc rc À rS À 1 rc (rc À rS )(5rc À 3rS ) This expression can be rewritten as a quartic equation for rc of the form
4 2 N rc À Or 3 þ P rc À Qrc þ R ? 0; c

ð27Þ

where O ? 16 À 3 þ 1 5À 2 P ? 12 þ l 2 À1 8 2 Q? l 2; R ? À1 N ? 5; 2 ; À1 À 6; þ6 2 l; À1

ð28Þ

We can now construct the wind equation by combining equations (10), (17), and (18), which yields 1 dv N ?; v dr D ð19Þ

where the numerator and denominator functions N and D are given by ! GM l2 a 2 3r S À 5r 2a 2 N? À v 2: À 3þ ; D? 2 r þ 1 r (r À rS ) þ1 (r À rS ) ð20Þ The simultaneous vanishing of N and D yields the critical conditions ! 2 GM l2 ac 3r S À 5r c Àþ ? 0; ð21Þ ðrc À rS Þ2 r 3 þ 1 rc ðrc À rS Þ c
2 2ac 2 À v c ? 0; þ1

ð22Þ

where vc and ac denote the values of the velocity and the sound speed at the critical radius, r ? rc .

and we have used natural gravitational units with GM ? c ? 1 and rS ? 2. These equations agree with the corresponding results derived by Das et al. (2001). The four solutions for rc in terms of the three fundamental parameters , l, and can be obtained analytically using the standard formulas for quartic equations (e.g., Abramowitz & Stegun 1970). We refer to the roots using the notation rc1, rc 2 , rc 3 , and rc 4 in order of decreasing radius. The critical radius rc 4 always lies inside the event horizon and is therefore not physically relevant, but the other three are located outside the horizon. The type of each critical point is determined by computing the two possible values for the de^ rivative dv/dr at the corresponding location using L'Hopital's rule and then checking to see whether they are real or complex. We find that both values are complex at rc 2 , and therefore this is an O-type critical point, which does not yield a physically acceptable solution. The remaining roots rc1 and rc 3 each possess real derivatives and are therefore physically acceptable sonic points, although the types of accretion flows that can pass through them are different. Specifically, rc 3 is an X-type critical point, and therefore a smooth, global, shock-free solution always exists in which the flow is transonic at rc 3 and then remains supersonic all the way to the event horizon. On the other hand, rc1 is an -type critical point, and therefore any accretion flow originating at a large distance that passes through this


No. 1, 2005

OUTFLOWS FROM ADVECTION-DOMINATED DISKS

481

point must display a shock transition below rc1 (Abramowicz & Chakrabarti 1990). After crossing the shock, the subsonic gas must pass through another (-type) critical point before it enters the black hole, since the flow has to be supersonic at the event horizon ( Weinberg 1972). 3.5. Shock-free Solutions Even when a shock can exist in the flow, it is always possible to construct a globally smooth (shock-free) solution using the same set of parameters. Smooth flows must pass through the inner critical point located at radius rc 3, which is calculated using the quartic equation (27) for given values of , l,and . The corresponding values for the critical velocity, vc 3 , the critical sound speed, ac 3 , and the critical entropy, Kc 3 , are computed using equations (23), (24), and (25), respectively. Since the flows treated here are inviscid, they have a conserved value for the entropy parameter K (eq. [15]) unless a shock is present. Hence, in a smooth, shockfree flow the value of K is everywhere equal to the critical value Kc 3 . The structure of the velocity profile in a shock-free disk can therefore be determined using a simple root-finding procedure as follows. By eliminating the sound speed a between equations (13) and (15), we obtain the equivalent expression ?
2 1 2 1 l2 1 K c3 vþ þ 2 3 (r À r )2 v 2 2r À1 r S

^ stead, the new inner critical radius, which we refer to as rc3 , must be computed using the downstream value of the energy transport parameter + . We point out that the total energy inflow rate across the horizon, including the rest mass contribution, must be positive, since no energy can escape from the black hole, and therefore we require that c 2 þ þ > 0. 4.1. Isothermal Shock Jump Conditions We assume that the escape of the relativistic particles from the disk results in a negligible amount of mass loss, because the Lorentz factor of the escaping particles is much greater than unity (see Table 2 below). This is confirmed ex post facto by comparing the rate of mass loss, Mloss , with the accretion rate M . We find that for the models analyzed here, Mloss /M P 10À3 , and therefore our assumption of negligible mass loss is justified. Hence, the accretion rate M is essentially conserved as the gas crosses the shock, which is represented by the condition Á M lim M (rà À ") À M (rà þ ") ? 0;
"!0

ð30Þ

!
2

( À1)=( þ1)

þ È;

ð29Þ

where the symbol Á is used to denote the difference between post- and preshock quantities. The specific angular momentum l J /M is also conserved throughout the flow, and therefore we find that Á J ? 0: ð31Þ

where we have set K ? Kc3 . In general, at any radius r equation (29) yields one subsonic root and one supersonic root for the velocity. The subsonic solution is chosen for r > rc3,and the supersonic solution is selected for r < rc3. Once the velocity profile v (r) has been computed, we can obtain the corresponding sound speed distribution a(r) by using equation (15) with K ? Kc3 . Note that the velocity and sound speed solutions can also be calculated by integrating the wind equation (19) numerically, and the results obtained using this approach agree with the root-finding method. 4. ISOTHERMAL SHOCK MODEL Our primary goal in this paper is to analyze the acceleration of relativistic particles due to the presence of a standing, isothermal shock in an accretion disk. Hence, we are interested in flows that pass smoothly through the outer critical radius rc1 and then experience a velocity discontinuity at the shock location, which we refer to as rà . In order to form self-consistent global models, we first need to understand how the structure of the disk responds to the presence of a shock. This requires analysis of the shock jump conditions, which are based on the standard fluid dynamical conservation equations. Since shocks are always optional even when they can occur, we compare our results for the relativistic particle acceleration with those obtained when there is no shock and the flow is globally smooth. The values of the energy transport parameter on the upstream and downstream sides of the isothermal shock at r ? rà are denoted by À and + , respectively. Note that À > þ as a consequence of the loss of energy through the upper and lower surfaces of the disk at the shock location. It is important to emphasize that the drop in at the shock has the effect of altering the transonic structure of the flow in the postshock region. Hence, although the postshock flow must pass through another critical point and become supersonic before crossing the event horizon, the new (inner) critical point is not equal to any of the four roots computed using the upstream energy transport parameter À . In-

Furthermore, the radial momentum transport rate, defined by I 4rH (P þ v 2 ); ð32Þ

must remain constant across the shock, and consequently Á I ? 0: ð33Þ

Based on equations (13) and (31), we find that the jump con dition for the energy transport rate E is given by 1 Áv 2 þ 1 Áa 2 : ÁE ? M 2 À1 ð34Þ

Equations (4), (8), and (13) can be combined with equations (30), (33), and (34) to obtain þ v þ aþ ? À v À aÀ ;
2 2 aþ P þ À aÀ P À ? aÀ À v À À aþ þ v þ ;

ð35Þ ð36Þ ð37Þ

þ À À ?

2 2 2 2 v þ À v À aþ À aÀ þ ; 2 À1

where the minus and plus subscripts refer to quantities measured just upstream and just downstream from the shock, respectively. In the case of an isothermal shock, aþ ? aÀ , and therefore the shock jump conditions reduce to þ vþ ? À vÀ ; Pþ À PÀ ? þ À À ?
2 2 À vÀ À þ vþ 2 2 vþ À vÀ

ð38Þ ; ð39Þ ð40Þ

2

:


482

LE & BECKER

Vol. 632

Combining equations (38) and (39) and substituting for the pressure P using equation (9) yields the velocity jump condition vþ ? À1 MÀ2 < 1; À vÀ MÀ vÀ ; aÀ ð41Þ

sound speed ac1, and the critical entropy Kc1 are then calculated using equations (23), (24), and (25), respectively. Note that since the flow is adiabatic everywhere in the preshock region, it follows that KÀ ? Kc1 : ð47Þ

where MÀ is the incident Mach number of the shock. The corresponding result for the shock compression ratio RÃ is RÃ þ ? M 2 > 1: À À ð42Þ

Hence, the gas density increases across the shock, as expected. Based on equations (15) and (41) and the fact that aþ ? aÀ in an isothermal shock, we find that the jump condition for the entropy parameter K is given by Kþ ? À1 MÀ2 < 1; À KÀ ð43Þ

The profiles of the velocity v (r) and the sound speed a(r)inthe preshock region can therefore be calculated using a root-finding procedure based on equation (29), or, alternatively, by integrating numerically the wind equation (19). The next step is the selection of an initial guess for the shock radius rà and the calculation of the associated shock Mach number MÀ vÀ /aÀ using the preshock dynamical solutions for v (r)and a(r). Based on the value of MÀ, we can compute the jump in the entropy parameter K using equation (43), and consequently we find that the entropy in the downstream region is given by Kþ ? Kc1 ; M2 À ð48Þ

which indicates that entropy is lost from the disk at the shock location due to the escape of the particles that form the outflow ( jet). We can also make use of equation (41) to rewrite the jump condition for the energy transport parameter (eq. [40]) as v2 Á þ À À ? À 2 ! 1 À 1 < 0: 2 M4 À ð44Þ

where we have also used equation (47). In order to determine whether the initial guess for rà is selfconsistent, we employ a second, independent procedure for calculating the entropy in the downstream region based on the critical nature of the flow. In this approach, the downstream energy parameter + is calculated using the jump condition given by equation (44), which yields v2 þ ? À þ À 2 ! 1 À1 : 2 M4 À ð49Þ

The associated rate at which energy escapes from the disk at the isothermal shock location (the ``shock luminosity'') is given by (see eq. [13]) Lshock ÀÁ E ? ÀM Á / ergs sÀ1 : ð45Þ

Eliminating Á between equations (44) and (45) yields the alternative result ! 2 Lshock vÀ 1 1À 2 4 : ? ð46Þ 2 M MÀ 4.2. Shock Point Analysis For a given value of , it is known that smooth, shock-free global flow solutions exist only within a restricted region of the (À, l ) parameter space, and isothermal shocks can occur only in a subset of the smooth-flow region. In order for a shock to exist in the flow, it must be located between two critical points, and it must also satisfy the jump conditions given by equations (41), (43), and (44). The procedure for determining the disk /shock structure is summarized below. The process begins with the selection of values for the fundamental parameters À, l, and . The values of À and l are ultimately constrained by the observations of a specific object, as discussed in x 7. Following Narayan et al. (1997), we assume an approximate equipartition between the gas and magnetic pressures, and therefore we set ? 1:5. The first step in the determination of the shock location is the computation of the outer critical point location rc1 using the quartic equation (27). The associated values for the critical velocity vc1, the critical

We use this value to compute the downstream critical point ra^ dius rc3 based on the quartic equation (27). The associated values ^ for the critical velocity vc3 , the critical sound speed ac3 , and the ^ ^ critical entropy Kc3 are then calculated using equations (23), (24), and (25), respectively. The final step is to compare the value of ^ Kc3 with that obtained for K+ using equation (48). If these two quantities are equal, then the shock radius rà is correct and the disk /shock model is therefore dynamically self-consistent. Otherwise, the value for rà must be iterated, and the search continued. Roots for rà can be found only in certain regions of the (À, l, ) parameter space. By combining the analysis of the shock location discussed above with the critical conditions developed in x 3, we are able to compute the structure of shocked and shock-free (smooth) disk solutions for a given set of parameters À, l,and .The resulting topology of the parameter space is depicted in Figure 2 for the case with ? 1:5, which is the main focus of this paper. Within region I, only smooth flows are possible, and in regions II and III both smooth and shocked solutions are available. No global flow solutions (either smooth or shocked) exist in region IV, with l > lmax . Inside region II, one root for the shock radius rà can be found, and in region III two shock solutions are available, although only one actual shock can occur in a given flow. It is unclear which of the two roots for rà in region III is preferred, since the stability properties of the shocks are not completely understood (e.g., Chakrabarti 1989a, 1989b; Abramowicz & Chakrabarti 1990). However, it is worth noting that the inner shock is always the stronger of the two, because the Mach number diverges as the gas approaches the horizon. The larger compression ratio associated with the inner location leads to more


No. 1, 2005

OUTFLOWS FROM ADVECTION-DOMINATED DISKS

483

Fig. 2.--Plot of the (À, l ) parameter space for an ADAF disk with ? 1:5. Only smooth flows exist in region I, and both shocked and smooth solutions are possible in regions II and III. When l > lmax (region IV ), no steady state dynamical solutions can be obtained. The parameters corresponding to models 2 and 5 are indicated. See the text for a complete discussion.

efficient particle acceleration and enhanced entropy generation, and therefore we expect that the inner solution is preferred in nature. Based on this argument, we focus on the inner shock location in our subsequent analysis. Before we proceed to examine the transport equation for the relativistic particles, it is important to analyze the asymptotic solutions obtained for the dynamical variables near the event horizon and also at large radii. This is a crucial issue, since the nature of the global solutions to the transport equation depends sensitively on the boundary conditions imposed at large and small radii. The asymptotic solutions for the dynamical variables near the event horizon were fully discussed by Becker & Le (2003). We briefly review their results and then perform a similar analysis in order to determine the asymptotic properties of the solutions as r ! 1. 4.3. Asymptotic Behavior Near the Horizon Becker & Le (2003) demonstrated that the variation of the global solutions in a viscous ADAF disk becomes purely adiabatic close to the event horizon, and therefore the asymptotic solutions that they obtain can be directly applied to our inviscid model. Using their equations (47) and (51), we find that the asymptotic variations of the radial velocity v and the sound speed a near the horizon are given by v 2 ( r ) / ( r À rS ) À1 ; a 2 (r) / (r À rS )(1 r ! rS :
À )=(1þ )

The divergence of v as r ! rS implies that it cannot represent the standard velocity in the region near the horizon. However, our dynamical model is consistent with relativity if we interpret v as the radial component of the four-velocity in this region ( Becker & Le 2003; Becker & Subramanian 2005). By combining these relations with equations (4), (8), and (10), we find that the corresponding results for the asymptotic variations of the disk half-thickness H and the density become H (r ) / (r À rS )
( þ3)=(2 þ2)

; (r) / (r À rS ) r ! rS :

À1=( þ1)

; ð51Þ

4.4. Asymptotic Behavior at Infinity We can use the energy transport equation (13) and the entropy equation (15) to obtain the asymptotic solutions for v and a at infinity as follows. In the limit r ! 1, the two dominant terms in equation (13) are and a 2 /( À 1), and therefore we find that a 2 (r) ! ( À 1); r ! 1: ð52Þ

Recalling that K is constant in the adiabatic upstream flow, we can combine equations (15) and (52) to conclude that the asymptotic variation of the inflow velocity is given by v/r
À5=2

; ð50Þ ; r ! 1: ð53Þ


484

LE & BECKER

Vol. 632

Finally, based on equations (8), (10), and (52), we find that the disk half-height varies as H /r
3=2

;

r ! 1:

ð54Þ

We can also combine equations (9), (14), and (52) to conclude that the asymptotic behavior of the density is given by ! const; r ! 1: ð55Þ

5. STEADY STATE PARTICLE ACCELERATION Our goal in this paper is to analyze the transport and acceleration of relativistic particles (ions) in a disk governed by the dynamical model developed in xx 3 and 4. For fixed values of the theory parameters À, l, and , we study the transport of particles in disks with and without shocks. The particle transport model used here includes advection, spatial diffusion, Fermi energization, and particle escape. In order to maintain consistency with the dynamics of the disk, we need to equate the energy escape rate for the relativistic particles with the shock luminosity Lshock given by equation (45). Our treatment of Fermi energization includes both the general compression related to the overall convergence of the accretion flow and also the enhanced compression that occurs at the shock. In the scenario under consideration here, the escape of particles from the disk occurs via vertical spatial diffusion in the tangled magnetic field, as depicted in Figure 1. To avoid unnecessary complexity, we use a simplified model in which only the radial (r) component of the spatial particle transport is treated in detail. In this approach, the diffusion and escape of the particles in the vertical (z) direction is modeled using an escape probability formalism. We treat the relativistic ions as test particles, meaning that their contribution to the pressure in the flow is neglected. This assumption is valid provided the pressure of the relativistic particles turns out to be a small fraction of the thermal pressure, as discussed in x 8. The ions accelerated at the shock are energized via collisions with MHD scattering centers (waves) advected along with the background (thermal) flow, and therefore the shock width is expected to be comparable to the magnetic coherence length, k mag. This approximation is used to determine the rate at which particles escape from the disk in the vicinity of the shock. 5.1. Transport Equation The Green's function, fG (E0 , E, rà , r), represents the particle distribution resulting from the continual injection of N0 particles per second, each with energy E0 , from a source located at the shock radius, r ? rà . In a steady state situation, the Green's function satisfies the transport equation ( Becker 1992) Á @ fG 1 @À 3 E v = :fG þ fsource À fesc ; ? 0 ? À: = F À 2 @E 3E @t ð56Þ where the specific flux F is evaluated using F ? À:fG À vE @ fG ; 3 @E ð57Þ

The quantities E, , v , and Hà H (rà ) represent the particle energy, the spatial diffusion coefficient, the vector velocity, and the disk half-thickness at the shock location, respectively, and the dimensionless parameter A0 determines the rate of particle escape through the surface of the disk at the shock location. The ^ ^ ^ vector velocity v has components given by v ? vr r þ vz z þ v f, and v ? Àvr is the positive inflow speed. The total number and energy densities of the relativistic particles, denoted by nr and Ur , respectively, are related to the Green's function via Z1 Z1 4E 2 fG dE ; Ur (r) ? 4E 3 fG dE ; ð59Þ nr ( r ) ?
0 0

which determine the normalization of fG . Equations (56) and (57) can be combined to obtain the alternative form v = :fG ? : = v @ fG E þ : = ð:fG Þ þ fso 3 @E
urce

À fesc ;

ð60Þ

where the left-hand side represents the comoving (advective) time derivative and the terms on the right-hand side describe first-order Fermi acceleration, spatial diffusion, the particle source, and the escape of particles from the disk at the shock location, respectively. Note that escape and particle injection are localized to the shock radius due to the presence of the -functions in equations (58). Our focus here is on the first-order Fermi acceleration of relativistic particles at a standing shock in an accretion disk, and therefore equation (60) does not include second-order Fermi processes that may also occur in the flow due to MHD turbulence (e.g., Schlickeiser 1989a, 1989b; Subramanian et al. 1999). Under the assumption of cylindrical symmetry, equations (58) and (60) can be rewritten as ! @ fG @ fG 1 1 @ d vz @ fG 1 @ @ fG r þ vz À À ðr v r Þ þ E vr 3 r @r @r @z dz @E r @r @r 0 (E À E0 )(r À rà ) N À A0 c(r À rà ) fG ; ? (4E0 )2 rà Hà ð61Þ where the escape of particles from the disk is described by the final term on the right-hand side. In Appendix A, we demonstrate that the vertically integrated transport equation is given by (see eq. [A9]) @ fG 1@ @ fG 1 @ @ fG rH ? þ ðrH vr ÞE H vr 3r @ r @r @E r @r @r N0 (E À E0 )(r À rà ) À A0 cHà fG (r À rà ); ð62Þ þ (4E0 )2 rà where the symbols fG , vr , and represent vertically averaged quantities. We establish in Appendix B that within the context of our one-dimensional spatial model, the dimensionless escape parameter A0 appearing in equation (62) is given by (see eqs. [ B8] and [ B10]) 3 à 2 < 1; ð63Þ A0 ? cHà where à (À þ þ )/2 denotes the mean value of the diffusion coefficient at the shock location. The condition A0 < 1 is required for the validity of the diffusive picture we have employed in our model for the vertical transport.

and the source and escape terms are given by N0 (E À E0 )(r À rà ) fsource ? ; fesc ? A0 c(r À rà ) fG : (4E0 )2 rà Hà ð58Þ


No. 1, 2005

OUTFLOWS FROM ADVECTION-DOMINATED DISKS 5.2. Number and Energy Densities

485

The energy moments of the Green's function, In(r), are defined by Z1 4E n fG dE; ð64Þ In (r)
0

so that (cf. eqs. [59]) nr (r) ? I2 (r); Ur (r) ? I3 (r): ð65Þ

R1 By operating on equation (62) with 0 4E n dE and integrating by parts once, we find that the function In satisfies the differential equation dIn n þ 1 In d 1d dIn rH ?À H vr ðrH vr Þ þ 3 r dr dr r dr dr n À2 N0 E 0 (r À rà ) À A0 cHà In (r À rà ); ð66Þ þ 4 r à which can be expressed in the flux conservation form " d 2 À n dIn N0 E nÀ2 (r À rà ) 0 (4rHFn ) ? 4rH þ v dr 3 4 r à H à dr # À A0 c(r À rà )In ; ð67Þ

where 0 and are dimensionless constants. Due to the appearance of the inflow speed v in equation (70), we note that exhibits a jump at the shock. This is expected on physical grounds, since the MHD waves that scatter the ions are swept along with the thermal background flow, and therefore they should also experience a density compression at the shock. As discussed above, close to the event horizon inward advection at the speed of light must dominate over outward diffusion. Conversely, in the outer region, we expect that diffusion will dominate over advection as r ! 1. By combining equation (70) with the asymptotic velocity variations expressed by equations (50) and (53), we find that the conditions given by equations (69) are satisfied if > 1, and in our work we set ? 2. Note that the escape parameter A0 is related to 0 via equation (63), which can be combined with equation (70) to write 30 và r A0 ? cHÃ
S

2

rà À1 rS

4 < 1; ð71Þ

where và (vÀ þ vþ )/2 represents the mean velocity at the shock location r ? rà . The value of the diffusion parameter 0 is constrained by the inequality in equation (71). In x 7 we demonstrate that 0 can be computed for a given source based on energy conservation considerations. With the introduction of equations (70) and (71), we have completely defined all of the quantities in the transport equation, and we can now solve for the number and energy densities of the relativistic particles. The particle distribution Green's function fG and its applications will be discussed in a separate paper. 6. SOLUTIONS FOR THE ENERGY MOMENTS Once the disk /shock dynamics have been computed based on the selected values for the free parameters À, l,and using the results in xx 3 and 4, the associated solutions for the number and energy densities of the relativistic particles in the disk can be obtained by solving equation (66). In the case of the number density, nr ? I2 , an exact solution can be derived based on the linear first-order differential equation describing the conservation of particle flux. However, in order to understand the variation of the energy density, Ur ? I3 , we must numerically integrate a second-order equation. The solutions obtained below are applied in x 7 to model the outflows observed in M87 and Sgr Aà . 6.1. Relativistic Particle Number Density The equation governing the transport of the particle number density nr is obtained by setting n ? 2 in equation (67), which yields d Nr ? N0 (r À rà ) À 4rà Hà A0 cnr (r À rà ); dr ð72Þ

where 4rHFn represents the rate of transport of the nth moment, and the flux Fn is defined by n þ1 dIn ; ð68Þ Fn À v In À 3 dr with v ? Àvr denoting the positive inflow speed. In order to close the system of equations and solve for the relativistic particle number and energy densities I2(r) and I3(r) using equation (66), we must also specify the radial variation of the diffusion coefficient . The behavior of can be constrained by considering the fundamental physical principles governing accretion onto a black hole. First, we note that near the event horizon, particles are swept into the black hole at the speed of light, and therefore advection must dominate over diffusion. This condition applies to both the thermal ( background) and the nonthermal (relativistic) particles. Second, we note that in the outer region (r ! 1), diffusion is expected to dominate over advection. Focusing on the flux equation for the particle number density nr , obtained by setting n ? 2 in equation (68), we can employ scale analysis to conclude based on our physical constraints that
r!r

lim
S

(r) ? 0; ( r À rS ) v ( r )

r!1

lim

r v (r) ? 0: (r )

ð69Þ

The precise functional form for the spatial variation of is not completely understood in the accretion disk environment. In order to obtain a mathematically tractable set of equations with a reasonable physical behavior, we use the general form r À1 ; ð70Þ (r ) ? 0 v (r )rS rS

where the relativistic particle transport rate Nr (r) is defined by (cf. eq. [68]) dnr ð73Þ Nr (r) À4rH v nr þ / s À1 ; dr and Nr > 0 if the transport is in the outward direction. Since the source is located at the shock, there are two spatial domains of interest in our calculation of the particle transport, namely,


486

LE & BECKER

Vol. 632

domain I (r > rà ) and domain II (r < rà ). Note that the number density nr(r) must be continuous at r ? rà in order to avoid generating an infinite diffusive flux according to equation (73). Away from the shock location, r 6? rà , and therefore equation (72) reduces to Nr ? const, which reflects the fact that particle injection and escape are localized at the shock. We can therefore write ( N I ; r > rà ; Nr (r) ? ð74Þ NII ; r < rà ; where the constant NI > 0 denotes the rate at which particles are transported outward, radially, from the source location, and the constant NII < 0 represents the rate at which particles are transported inward toward the event horizon. The magnitude of the jump in the particle transport rate at the shock is obtained by integrating equation (72) with respect to radius in a very small region around r ? rà , which yields NI À NII ? N0 À Nesc ; where N
esc

Furthermore, in order to avoid exponential divergence of nr as r ! rS in domain II, we also require that nà ? ÀNII CII ; where CII 1 4 Z
r r
S Ã

ð82Þ

e J(r ) 0 dr : r 0H

0

ð83Þ

By combining equations (75), (76), (80), and (82), we can develop explicit expressions for the quantities nà , NI , NII , and Nesc based on the values of rà and N0 and the profiles of the inflow velocity v (r) and the diffusion coefficient (r). The results obtained are nà ? N0 ; þ C þ 4rà Hà A0 c nà nà NI ? ; NII ? À ; CI CII Nesc ? 4rà Hà A0 cnà : CIÀ1
À1 II

ð75Þ

ð84Þ

4rà Hà A0 cn

Ã

ð76Þ

represents the ( positive) rate at which particles escape from the disk at the shock location to form the outflow ( jet), and nà nr (rà ). If no shock is present in the flow, then A0 ? 0, and there fore Nesc ? 0. Note that the discontinuity in Nr at the shock produces a jump in the derivative dnr /dr via equation (73). We can rewrite equation (73) for the number density in the form Nr dnr v þ nr ? À ; dr 4rH ð77Þ

These relations, along with equation (78), complete the formal solution for the relativistic particle number density nr(r). The solution is valid in both shocked and shock-free disks (the shockfree case is treated by setting A0 ? 0). When a shock is present, the particle escape rate Nesc is proportional to N0 but is independent of E0 by virtue of equations (84). It is interesting to examine the asymptotic variation of nr near the event horizon and also at large distances from the black hole. Far from the hole, advection is negligible, and the particle transport in the disk is dominated by outward-bound diffusion. In this case we can use equation (77) to conclude that NI dnr !À ; dr 4rH r ! 1; ð85Þ

which is a linear, first-order differential equation for nr (r). Using the standard integrating factor technique and employing equation (70) for yields the exact solution nr ( r ) ? e
ÀJ ( r )

where we have used the fact that Nr ? NI for r > rà . By combining equations (70) and (85) with the asymptotic relations given by equations (53) and (54), we find upon integration that nr ( r ) / 1 ; r r ! 1: ð86Þ

! Z 0 Nr (r) r e J ( r ) 0 dr ; nà À 4 r à r 0 H

ð78Þ

where Nr (r) is given by equation (74) and the function J (r) is defined by " 1 1# Zr À À rà r À1 0 À1 J (r ) v dr ? 0 À1 À À1 : ð79Þ rS rS rà According to equation (78), nr (r) is continuous at the shock / source location, as required. Far from the black hole, diffusion dominates the particle transport, and therefore nr should vanish as r ! 1. In order to ensure this behavior, we must have nà ? NI CI ; where CI 1 4 Z
r 1
Ã

In order to study the behavior of nr near the event horizon, we take the limit as r ! rS in equation (78), obtaining after some algebra nr ( r ) ! À NII ; 4rH v r ! rS ; ð87Þ

where we have set Nr ? NII . Comparing this relation with equation (4), we find that nr (r) / (r); r ! rS ; ð88Þ

ð80Þ

e J(r ) 0 dr : r 0H

0

where is the density of the background (thermal) gas. Equation (88) is a natural consequence of the fact that the particle transport near the horizon is dominated by inward-bound advection. We can also combine equations (51) and (88) to obtain the explicit asymptotic form nr (r) / (r À rS )À
1=( þ1)

ð81Þ

;

r ! rS :

ð89Þ


No. 1, 2005

OUTFLOWS FROM ADVECTION-DOMINATED DISKS

487

6.2. Relativistic Particle Energy Density The differential equation satisfied by the relativistic particle energy density Ur ? I3 is obtained by setting n ? 3in equation (66), which yields Hv dUr 4 Ur d 1d dUr rH ?À ðrH vr Þ þ 3 r dr r dr dr dr N0 E0 (r À rà ) À A0 cHà Ur (r À rà ): þ 4 r Ã

we find that Ur / nr . We can therefore use equation (86) to conclude that Ur (r) / 1 ; r r ! 1: ð96Þ

r

ð90Þ

Close to the event horizon, the particle transport is dominated by advection, and therefore Ur and nr obey the standard adiabatic relation Ur / n
4=3 r

;

r ! rS :

ð97Þ

By analogy with equations (72) and (73), we can recast this expression in the flux conservation form ! d Er v dUr N0 E0 (r À rà ) ? 4rH À þ À A0 cUr (r À rà ) ; 3 dr 4 r à H à dr ð91Þ where the relativistic particle energy transport rate Er (r) is defined by 4 dUr vUr þ ð92Þ Er (r) À4rH / ergs sÀ1 ; 3 dr and Er > 0 for outwardly directed transport. Note that unlike the number transport rate Nr , the energy transport rate Er does not remain constant within domains I and II due to the appearance of the first term on the right-hand side of equation (91), which expresses the compressional work done on the relativistic particles by the background flow. Although the energy density Ur must be continuous at the shock /source location in order to avoid generating an infinite diffusive flux, the derivative dUr /dr displays a discontinuity at r ? rà , which is related to the jump in the energy transport rate via equation (92). By integrating equation (91) in a very small region around r ? rà , we find that Á Er ? Lesc À N0 E0 ; where Lesc 4rà Hà A0 cUà / ergs s
À1

Combining this result with equations (88) then yields Ur / (r À rS )
À4=(3 þ3)

;

r ! rS :

ð98Þ

The global solution for Ur (r) can now be expressed as & AQI (r); r > rà ; Ur (r) ? ð99Þ BQII (r); r < rà ; where A and B are constants and the functions satisfy the homogeneous differential equation dQ 4Q d 1d ?À rH ðrH vr Þ þ H vr dr 3 r dr r dr QI (r) and QII (r) (see eq. [90]) dQ ; ð100Þ dr

along with the boundary conditions (see eqs. [96] and [98]) rou QI (rout ) ? rS
t

À

1

;

rin QII (rin ) ? À1 rS



À

4=(3 þ3)

;

ð101Þ

where rin and rout denote the radii at which the inner and outer boundary conditions are applied, respectively. The constants A and B are computed by requiring that Ur be continuous at r ? rà and that the derivative dUr /dr satisfy the jump condition given by equation (95). The results obtained are QII A?B ; QI r ? rà N0 E0 4 (vÀ À vþ )QII B? 4 r à H à 3 þ
0 QII þ

ð93Þ

ð102Þ

ð94Þ

denotes the rate of escape of energy from the disk into the outflow ( jet) at the shock location, and Uà Ur (rà ). If no shock is present, then A0 ? 0 and therefore Lesc ? 0. We remind the reader that the symbol Á refers to the difference between postand preshock quantities (see eq. [30]). Equations (92) and (93) can be combined to show that the derivative jump is given by N0 E0 À L dUr ? Á dr 4 r à H Ã
esc

QII QI0 À À þ A0 cQII QI

! À1

;
r?r
Ã

ð103Þ

4 À U Ã Áv: 3

ð95Þ

The differential equation (90) governing the relativistic particle energy density is second order in radius, and therefore we need to establish two boundary conditions in order to solve for Ur (r). These can be obtained by analyzing the behavior of Ur close to the event horizon and at large distances from the black hole. Far from the hole, advection is negligible, and the particle transport in the disk is dominated by outward-bound diffusion. In this regime, Fermi acceleration is negligible, and consequently

where primes denote differentiation with respect to radius. The solutions for the functions QI (r)and QII (r) are obtained by integrating equation (100) numerically subject to the boundary conditions given by equations (101). Once the constants A and B are computed using equations (102) and (103), the global solution for Ur (r) is evaluated using equation (99). The shockfree case is treated by setting A0 ? 0. This completes the solution procedure for the relativistic particle energy density. The results derived in this section are used in x 7 to model the outflows observed in M87 and Sgr AÃ . 7. ASTROPHYSICAL APPLICATIONS Our goal is to determine the properties of the integrated disk / shock /outflow model based on the observed values for the black


488

LE & BECKER
TABLE 1 Disk Structur e P arameters Model 1.................................. 2.................................. 3.................................. 4.................................. 5.................................. 6..................................
À

Vol. 632

l 3.1040 3.1340 3.1765 3.1280 3.1524 3.1756

+ À0.005676 À0.005746 À0.005875 À0.008723 À0.008749 À0.008778

rc1 93.177 98.524 109.781 131.204 131.874 135.192

rc3 5.478 5.379 5.252 5.408 5.329 5.260

rà 19.624 21.654 24.500 14.073 15.583 16.950

rc ^

3

H

Ã

M

À

RÃ 1.795 1.897 2.052 1.857 1.970 2.076

TÃ 1.28 1.16 1.01 1.65 1.49 1.36

0.001600 0.001527 0.001400 0.001240 0.001229 0.001200

5.798 5.659 5.487 5.850 5.723 5.614

10.369 11.544 13.154 6.819 7.672 8.434
11

1.094 1.125 1.170 1.113 1.146 1.177 K.

Note.--All quantities are expressed in gravitational units (GM ? c ? 1) except TÃ , which is in units of 10

hole mass M, the mass accretion rate M , and the jet kinetic power Ljet associated with a given source. The fundamental free parameters for the theoretical model are À, l, and . Since we set ? 1:5 in order to represent an approximate equipartition between the gas and magnetic pressures (e.g., Narayan et al. 1997), only À and l remain to be determined. Here we describe how global energy conservation considerations can be used to solve for the various theoretical parameters in the model based on observations. 7.1. Energy Conservation Conditions Once the values of M, M , and Ljet have been specified for a source based on observations, we select a value for the free parameter À and then compute l by satisfying the relation Lshock ? Ljet ; ð104Þ

particle energy escapes from the disk due to vertical diffusion is equal to the observed jet power. This condition can be written as L
esc

? Ljet ;

ð106Þ

where Lesc is the energy escape rate given by equation (94). The escape constant A0 appearing in the transport equation is independent of the particle energy in our model, and consequently the escaping particles will have exactly the same mean energy as those in the disk at the shock location. The mean energy of the escaping particles is therefore given by E
esc



UÃ ; nÃ

ð107Þ

where Lshock is the shock luminosity given by equation (46). This result ensures that the jump in the energy transport rate at the isothermal shock location is equal to the observed jet kinetic luminosity. The procedure for determining l also includes solving for the shock location and the critical structure using results from xx 3 and 4. The velocity profile v (r) is computed either by numerically integrating equation (19) or by using a root-finding procedure based on equation (29), and the associated solution for the adiabatic sound speed a(r) is obtained using equation (13). After the velocity profile has been determined, we can compute the number and energy density distributions for the relativistic particles in the disk using equations (78) and (99), respectively. This requires the specification of the injection energy of the seed particles E0 as well as their injection rate N0 . We set the injection energy using E0 ? 0:002 ergs, which corresponds to an injected Lorentz factor À0 E0 /(mp c 2 ) $ 1:3, where mp is the proton mass. Particles injected with energy E0 are subsequently accelerated to much higher energies due to repeated shock crossings. We find that the speed of the injected particles, v0 ? c(1 À ÀÀ2 )1/2 , 0 is about 3 - 4 times higher than the mean ion thermal velocity at the shock location, vrms ? (3kTÃ /mp )1/2 , where TÃ is the ion temperature at the shock. The seed particles are therefore picked up from the high-energy tail of the Maxwellian distribution for the thermal ions. With E0 specified, we can compute the particle in jection rate N0 using the energy conservation condition N0 E0 ? Lshock ; ð105Þ

where nà and Uà denote the number and energy densities of the relativistic particles at the shock location, respectively. Hence, Eesc is proportional to E0 , but it is independent of N0 . We note that equations (76), (94), and (107) can be combined to show that L
esc

? Nesc Eesc ;

ð108Þ

where Nesc is the particle escape rate (eq. [76]). By satisfying equations (104), (105), and (106), we ensure that energy is properly conserved in our model. Taken together, these relations allow us to solve for the various theoretical parameters based on obser vational values for M, M , and Ljet , as explained below. 7.2. Model Parameters Our simulations of the disk structure and particle transport in M87 and Sgr Aà are based on various published observational estimates for M, M , and Ljet. For M87, we set M ? 3 ; 109 M (e.g., Ford et al. 1994), M ? 1:3 ; 10À1 M yrÀ1 (e.g., Reynolds et al. 1996), and Ljet ? 5:5 ; 1043 ergs sÀ1 (e.g., Reynolds et al. 1996; Bicknell & Begelman 1996; Owen et al. 2000). For Sgr Aà , Å we use the values M ? 2:6 ; 106 M (e.g., Schodel et al. 2002) and M ? 8:8 ; 10À7 M yrÀ1 (e.g., Yuan et al. 2002; Quataert 2003). Although the kinetic luminosity of the jet in Sgr Aà is rather uncertain (see, e.g., Yuan 2000; Yuan et al. 2002), we adopt the value quoted by Falcke & Biermann (1999) and set Ljet ? 5 ; 1038 ergs sÀ1. We study both shocked and shock-free solutions spanning the computational domain between the inner radius rin ? 2:001 and the outer radius rout ? 5000, where rin and rout are the radii at which the boundary conditions are applied (see eqs. [101]). Six different accretion /shock scenarios are explored in detail, with the values for the various parameters À, l, +, rc1, rc3, rà , ^ rc3 , Hà , MÀ , Rà , and Tà reported in Table 1. Models 1, 2, and 3 are associated with M87, while models 4, 5, and 6 are used to study Sgr Aà . In our numerical examples, we use natural gravitational units (GM ? c ? 1 and rS ? 2), except as noted. Based

which ensures that the rate at which energy is injected into the flow in the form of the relativistic seed particles is equal to the energy-loss rate for the background gas at the isothermal shock location. In order to maintain agreement between the transport model and the observations, we must also require that the rate at which


No. 1, 2005

OUTFLOWS FROM ADVECTION-DOMINATED DISKS

489

Fig. 3.-- Quantities Lesc /Ljet , Àesc ,and Nesc /N0 for a shocked disk plotted as functions of 0 . Panels (a), (b), and (c) correspond to the M87 parameters (models 1, 2, and 3), and panels (d ), (e), and ( f ) correspond to the Sgr AÃ parameters (models 4, 5, and 6). The model numbers are indicated for each curve. See the discussion in the text.

on the observational values for M and Ljet associated with the two sources, we can use equations (45) and (104) to conclude that Á ? À0:007 for M87 and Á ? À0:01 for Sgr AÃ . These results are consistent with the values for À and + reported in Table 1, and therefore Lshock ? Ljet , as required (see eq. [104]). Next, we use the energy conservation condition Lesc ? Ljet (eq. [106]) to determine the value of the diffusion constant 0 (eq. [70]) for a shocked disk. In Figure 3, we plot Lesc /Ljet , Àesc , and Nesc /N0 as functions of 0 for the M87 and Sgr AÃ parameters, where Àesc Eesc /(mp c 2 ) is the mean Lorentz factor of the escaping particles. The treatment of energy conservation in our disk /shock model is self-consistent when the condition Lesc /Ljet ? 1 is satisfied, which corresponds to specific values of 0, as indicated in Figures 3a and 3d. We find that two 0 roots exist for models 3 and 6, one root is possible for models 2 and 5, and no roots exist for models 1 and 4. Hence, the values of À associated with models 2 and 5 represent the maximum possible values for À that yield self-consistent solutions based on the M87 and Sgr AÃ data, respectively. For illustrative purposes, we focus on the details of the disk structure and particle transport obtained in models 2 and 5. 7.3. Disk Structure and Particle Transport In order to illustrate the importance of the shock for the acceleration of high-energy particles, we examine the structure of the accretion disk both with and without a shock based on the values of the upstream parameters À and l used in models 2 and 5 (see Table 1). In Figures 4a and 4b we plot the inflow speed v (r) and the adiabatic sound speed a (r) for the shocked and

smooth (shock-free) solutions associated with models 2 and 5, respectively. Since we are working within the isothermal shock scenario, the sound speed a is continuous at the shock location. In Figures 4c and 4d we plot the specific internal energy U / ? ( À 1)À1 kT /mp for the background (thermal) gas along with the specific gravitational potential ( binding) energy GM /(r À rS ) as functions of radius for ? 1:5. These results demonstrate that the gas is marginally bound in the absence of a shock and strongly bound when a shock is present. The increased binding of the thermal gas in the disk results from the escape of energy in the outflow, which reduces the sound speed compared with the shock-free case. The enhanced cooling allows the accretion to proceed, thereby removing one of the major objections to the original ADAF scenario ( Narayan & Yi 1994, 1995). We emphasize that these new results represent the first fully self-consistent calculations of the structure of an ADAF disk coupled with a shock-driven outflow, hence extending the heuristic work of Blandford & Begelman (1999) and Becker et al. (2001). Next, we study the solutions obtained for the relativistic particle number and energy density distributions in the disk based on the flow structures associated with models 2 and 5. The related transport parameters are listed in Table 2. In Figures 5 and 6 we plot the global number and energy density distributions obtained in a shocked disk using the model 2 and 5 parameters, respectively. We also include the corresponding results obtained in a shock-free (smooth) disk for the same values of the upstream energy transport rate À and the specific angular momentum l. In each case the densities decrease monotonically with increasing radius. The increase near the horizon is a consequence of


490

LE & BECKER

Fig. 4.--Velocity v (r) (solid lines) and sound speed a(r)?2/( þ 1)1 2 (dashed lines) plotted in units of c for (a) model 2 and (b) model 5. The curves cross at the critical points. Also plotted are the specific internal energy U/ (solid lines) and the specific gravitational potential energy GM /(r À rS ) (dashed lines), both in units of c 2, for (c) model 2 and (d ) model 5. The shocked and shock-free solutions are denoted by the thin and thick lines, respectively.
=

advection, while the decline as r ! 1 reflects the fact that the particles injected at the shock have a very small chance of diffusing to large distances from the black hole. Note that the shocked disk has a lower value for the number density nr at all radii as a consequence of particle escape. However, the shocked disk also displays a higher value for the energy density Ur , which reflects the central role of the shock in accelerating the relativistic test particles. The kinks in the energy and number density distributions at the shock radius r ? rà indicated in Figures 5 and 6 reflect the derivative jump conditions given by equations (75) and (95). The values for the ratios NI /NII and Nesc /N0 reported in Table 2 indicate that most of the injected particles are advected into the black hole, with $20% escaping to form the outflow (see Figs. 3c and 3f ). In order to validate the accuracy of the numerical solutions for nr (r) and Ur (r), we also compare the profiles obtained with the asymptotic relations developed in x 6. We demonstrate in Figure 7 (model 2) and Figure 8 (model 5) that the solutions for both nr (r) and Ur (r) agree closely with the asymptotic expressions given by equations (89) and (98) for small radii and by

equations (86) and (96) for large radii. Note that the values reported by Le & Becker (2004) for nà and Uà were expressed in incorrect units and are given correctly in our Table 2. 7.4. Jet Formation in M87 and Sgr A
Ã

The mean energy of the relativistic particles in the disk is given by (cf. eq. [107]) hEi Ur (r) ; nr (r) ð109Þ

so that hEi? Eesc at r ? rà . In Figure 9 we plot the mean energy as a function of radius in shocked and shock-free disks based on the parameters used for models 2 and 5. The results demonstrate that when a shock is present in the flow, the relativistic particle energy is boosted by a factor of $5 - 6 at the shock location. By contrast, we find that in the shock-free models with the same values for À, l,and 0 the energy is boosted by a factor of only $1.4 -1.5. This clearly demonstrates the essential role of the shock in efficiently accelerating particles up to very high energies, far

TABLE 2 Transport Equation P arameters N0 (sÀ1) 2.75 ; 1046 2.51 ; 1041 nà (cmÀ3) 2:01 ; 104 4:33 ; 105 Uà (ergs cmÀ3) 2:39 ; 10 4:71 ; 10
2 3

Model 2.................................. 5..................................

NI /N

II

0 0.02044 0.02819



Ã

A0 0.0124 0.0158

Nesc /N 0.17 0.18

0

E

esc

/E0

À

esc

À0.18 À0.15

0.427877 0.321414

5.95 5.45

7.92 7.26

Note.--All quantities are expressed in gravitational units (GM ? c ? 1) except as noted.


Fig. 5.-- Global solutions for (a) the relativistic number density and (b) the relativistic energy density computed using the model 2 parameters. The solid and dashed curves correspond to disks with and without shocks, respectively. Note that the number density is higher in the smooth (shock-free) disk due to the absence of particle escape. Conversely, the energy density is higher in the shocked disk due to the enhanced particle acceleration occurring at the shock.

Fig. 6.--Same as Fig. 5, but for model 5.


Fig. 7.--Plots of the numerical solutions for nr (r) and Ur (r) (solid lines) computed using the model 2 parameters in a shocked disk compared with the asymptotic expressions (dashed lines) (a, b) close to the event horizon and (c, d ) at large radii. See the discussion in the text.

Fig. 8.--Same as Fig. 7, but for model 5.


OUTFLOWS FROM ADVECTION-DOMINATED DISKS

493

Fig. 9.--Mean energy of the relativistic particles in the disk, hEi Ur (r)/nr (r) (eq. [109]), plotted in units of the injection energy E0 as a function of radius for (a) model 2 and (b) model 5. Results are indicated for both shocked (solid lines) and shock-free (dashed lines) disk structures.

above the energy required to escape from the disk. Note that close to the event horizon, the mean energy of the relativistic particles is further enhanced by the strong compression of the accretion flow, as indicated by the sharp increase in hEi as r ! rS . The material in the outflow is initially ejected from the disk in the vicinity of the shock as a hot plasma that cools as it expands, with its outward acceleration powered by the pressure gradient in the surrounding plasma. Based on our results for models 2 and 5, we find that the shock /jet locations are given by rà $ 22 and rà $ 16 for M87 and Sgr Aà , respectively. The terminal (asymptotic) Lorentz factor of the jet À1 can be estimated by writing À
1

Sgr Aà , our model indicates that the shock forms at rà $ 16, which is fairly close to the value suggested by Yuan (2000). However, future observational work will be needed to test our prediction for the asymptotic Lorentz factor of Sgr Aà , since no reliable observational estimate for that quantity is currently available. 7.5. Radiative Losses from the Jet It is still unclear whether the outflows observed to emanate from many radio-loud systems containing black holes are composed of an electron-proton plasma or electron-positron pairs, or a mixture of both. Whichever is the case, the particles must maintain sufficient energy during their journey from the nucleus in order to power the observed radio emission, unless some form of reacceleration takes place along the way, e.g., due to shocks propagating along the jet (Atoyan & Dermer 2004a). Protonelectron outflows, such as those studied here, have a distinct advantage in this regard, since most of the kinetic power is carried by the ions, which do not radiate much and are not strongly coupled to the electrons under the typical conditions in a jet (e.g., Felten 1968; Felten et al. 1970; Anyakoha et al. 1987; Aharonian 2002). We therefore suggest that if the observed outflows are proton driven, then they may be powered directly by the shock acceleration mechanism operating in the disk, with no requirement for additional in situ reacceleration in the jet. In this section we confirm this conjecture by considering the energy losses experienced by the protons in the outflow. The ions in the jet lose energy via two distinct channels, namely, (1) direct radiative losses due to the production of synchrotron and inverse Compton emission and (2) indirect radiative losses via Coulomb coupling with

?À

esc

?

Eesc ; mp c 2

ð110Þ

which is based on the assumption that the jet starts off ``slow'' and ``hot'' and subsequently expands to become ``fast'' and ``cold.'' Adopting the Àesc values listed in Table 2 for M87 and Sgr Aà , we obtain À1 ? 7:92 (see Fig. 3b) and À1 ? 7:26 (see Fig. 3e), respectively. We can now compare our model predictions for the shock /jet location and the asymptotic Lorentz factor with the observations of M87 and Sgr Aà . According to Biretta et al. (2002), the M87 jet forms in a region no larger than $30 gravitational radii from the black hole, which agrees rather well with our predicted shock /jet location rà $ 22 for this source. Turning now to the asymptotic (terminal) Lorentz factor, we note that Biretta et al. (1999) estimated À1 ! 6 for the M87 jet, which is comparable to the result À1 ? 7:92 obtained using our model. In the case of


494

LE & BECKER

Vol. 632

the electrons. We evaluate these two possibilities by estimating the corresponding cooling timescales for the outflows in M87 and Sgr AÃ and comparing the results with the jet propagation timescales for these sources. The energy-loss rate due to the production of synchrotron and inverse Compton emission by the relativistic protons escaping from the disk with mean energy Àesc mp c 2 is given by (see eqs. [7.17] and [7.18] from Rybicki & Lightman 1979) dE dt 4 T c À ? 3
2 esc



rad

me mp

2

À Á UB þ Uph ;

ð111Þ

where rà , nà , Hà ,and k mag denote the radius, the proton number density, the vertical half-thickness, and the magnetic mean free path inside the disk at the shock location, respectively. The shock is expected to have a width comparable to k mag , and therefore the sum of the upper and lower face areas of the shock annulus is equal to 4rÃk mag . We also note that the flux of the relativistic protons escaping from the disk into the outflow is given by cnp , where np is the proton number density at the base of the jet. Combining these relations, we can write the proton escape rate in terms of np using Nesc ? 4rà k
mag

cnp :

ð116Þ

where Uph and UB ? B2 /(8) denote the energy densities of the soft radiation and the magnetic field with strength B, respectively. The associated energy-loss timescale is therefore t Àesc mp c 2 3m p c ? (dE=dt )jrad 4T Àesc mp me 2 À UB þ Uph Á
À1

By equating the two expressions for Nesc given by equations (115) and (116), we find that np is related to nà via np k mag ? < 1: nà Hà ð117Þ

rad

:

ð112Þ

In our application to M87, we take B $ 0:1 G based on estimates from Biretta et al. (1991), and we set Àesc $ 8 (see Table 2, model 2). Assuming equipartition between the magnetic field and the soft radiation, this yields for the radiative cooling time trad $ 1012 yr, which suggests that the protons can easily maintain their energy for many millions of parsecs without being seriously effected by synchrotron or inverse Compton losses. For Sgr AÃ , we assume equipartition with B $ 10 G (Atoyan & Dermer 2004b) and Àesc $ 7 (see Table 2, model 5). The radiative cooling time for the escaping protons is therefore trad $ 108 yr. Hence, synchrotron and inverse Compton losses have virtually no effect on the energy of the protons in the Sgr AÃ jet. In addition to synchrotron and inverse Compton radiation, the protons in the jet will also lose energy due to Coulomb coupling with the thermal electrons, which radiate much more efficiently than the protons. The energy-loss rate for this process can be estimated using equation (4.16) from Mannheim & Schlickeiser (1994), which yields dE dt ? 30ne T me c 3 ;
Coul

Since the electron-proton jet must be charge neutral, the electron number density at the base of the jet ne is equal to the proton number density np , and therefore we obtain ne ? kmag nà : Hà ð118Þ

ð113Þ

where ne represents the electron number density in the jet. The associated loss timescale for a proton escaping from the disk with mean energy Àesc mpc 2 is t
Coul



Àesc mp c 2 Àesc mp ? : (dE =dt )jCoul 30ne T cme

ð114Þ

1 Using the relation k mag /Hà ? A0=2 (see eq. [ B8]) along with the results for A0 and nà reported in Table 1 for M87 (model 2), we obtain ne ? 0:11nà ? 2:2 ; 103 cmÀ3. Setting Àesc $ 8, we find that equation (114) yields for the electron-proton Coulomb coupling timescale tCoul $ 3:5 ; 105 yr. Note that this is an extremely conservative estimate, since it is based on conditions at the bottom of the jet, and therefore it suggests that Coulomb coupling between the protons and the electrons is insufficient to seriously degrade the energy of the accelerated ions escaping from the disk as they propagate out to the radio lobes via the jet. For Sgr Aà , we use the model 5 data in Table 2 to obtain ne ? 0:13nà ? 5:6 ; 104 cmÀ3. Setting Àesc $ 7 yields for the Coulomb coupling timescale tCoul $ 1:2 ; 104 yr, which implies that the length of the jet can be as large as several thousand parsecs before much energy is drained from the protons, assuming the material in the jet travels at half the speed of light. We emphasize that these numerical estimates of the importance of radiative and Coulomb losses experienced by the relativistic protons are based on the ``worst-case'' assumption that the conditions at the base of the outflow prevail throughout the jet. In reality, the jet density will drop rapidly as the gas expands, and therefore the true values for the proton energy-loss timescales will be much larger than the results obtained above. This strongly suggests that shock acceleration of the protons in the disk, as investigated here, is sufficient to power the observed outflows without requiring any reacceleration in the jets.

The electron number density ne decreases rapidly as the jet expands from the disk into the external medium. Hence, the most conservative estimate ( based on the strongest Coulomb coupling) is obtained by adopting conditions at the base of the jet, where ne has its maximum value. To estimate the electron number density at the base of the outflow, we begin by calculating the rate at which protons escape from the disk at the shock location. By using equation ( B8) to eliminate A0 in equation (76), we find that the proton escape rate is given by Nesc ? 4 r à k
2 mag Ã

7.6. Radiative Losses from the Disk In the ADAF scenario that we have focused on, radiative losses from the disk are ignored. The self-consistency of this approximation can be evaluated by computing the free-free emissivity due to the thermal gas in the disk. The total X-ray luminosity can be estimated by integrating equation (5.15b) from Rybicki & Lightman (1979) over the disk volume to obtain for pure, fully ionized hydrogen Z1 1 1:4 ; 10À27 Te =2 2 mÀ2 dV ; ð119Þ Lrad ? p
r
S

cn

Ã

H

;

ð115Þ


No. 1, 2005

OUTFLOWS FROM ADVECTION-DOMINATED DISKS

495

where dV ? 4rH dr represents the differential (cylindrical) volume element and Te denotes the electron temperature. We can obtain an upper limit on the X-ray luminosity by assuming that the electron temperature is equal to the ion temperature. Based on the detailed disk structures associated with models 2 and 5, we find that Lrad /Ljet $ 10À2 and Lrad /Ljet $ 10À5 , respectively. However, in an actual ADAF disk the X-ray luminosity will of course be substantially smaller than these values, because the electron temperature is roughly 3 orders of magnitude lower than the ion temperature. Hence, our neglect of radiative losses is completely justified, as expected for ADAF disks. 8. CONCLUSION In this paper we have demonstrated that particle acceleration at a standing, isothermal shock in an ADAF accretion disk can energize the relativistic protons that power the jets emanating from radio-loud sources containing black holes. The work presented here represents a new type of synthesis that combines the standard model for a transonic ADAF flow with a self-consistent treatment of the relativistic particle transport occurring in the disk. The energy lost from the background (thermal) gas at the isothermal shock location results in the acceleration of a small fraction of the background particles to relativistic energies. One of the major advantages of our coupled, global model is that it provides a single, coherent explanation for the disk structure and the formation of the outflow based on the well-understood concept of first-order Fermi acceleration in shock waves. The theory employs an exact mathematical approach in order to solve simultaneously the combined hydrodynamic and particle transport equations. The analysis presented here closely parallels the early studies of cosmic-ray shock acceleration. As in those first investigations (e.g., Blandford & Ostriker 1978), we have employed an idealized model in which the pressure of the accelerated particles is assumed to be negligible compared with that of the thermal background gas (the ``test particle'' approximation). In order to check the self-consistency of this assumption, we have confirmed that the total pressure is dominated by the pressure of the background (thermal) gas throughout most of the disk. However, in the vicinity of the shock the two pressures can become comparable, and this suggests that the dynamical results will change slightly if the test particle approximation is relaxed. We plan to consider this question in future work by developing a ``two-fluid'' version of our model that includes the particle pressure, in analogy with the ``cosmic-ray modified shock'' scenario for cosmic-ray acceleration ( Becker & Kazanas 2001; Drury & Volk 1981). Å We have presented detailed results that confirm that the general properties of the jets observed in M87 and Sgr AÃ can be understood within the context of our disk /shock /outflow model.

In particular, our results indicate that the shock acceleration mechanism can produce relativistic outflows with terminal Lorentz factors and total powers comparable to those observed in M87 and Sgr AÃ . However, in principle even higher efficiencies can be achieved by varying the upstream energy transport rate À,which is the fundamental free parameter in our model. The buildup of the particle pressure in such high-efficiency situations would require relaxation of the test particle approximation, as discussed above. In this paper we have focused on inviscid disks, which are the simplest to analyze. While the inviscid model provides useful insight into the importance of shock acceleration in ADAF disks, this restriction clearly must be lifted in the future, since viscosity plays a key role in determining the structure of an actual accretion disk. We are currently developing a self-consistent viscous disk model in order to explore shock formation and particle acceleration in a more rigorous context. However, we do not expect the presence of viscosity to alter any of the basic conclusions reached in this paper, because significant particle acceleration will occur regardless of the viscosity, provided a shock is present. The existence of shocks in viscous disks is a controversial issue, but several studies suggest that shock formation is possible, provided the viscosity is relatively low. In the absence of a consensus regarding the possible presence of shocks in accretion disks, we believe that it is important to study models with shocks in order to develop theoretical predictions that can be tested observationally. The shock acceleration mechanism analyzed in this paper is effective only in rather tenuous, hot disks, and therefore we conclude that our model may help to explain the observational fact that the brightest X-ray AGNs do not possess strong outflows, whereas the sources with low X-ray luminosities but high levels of radio emission do. We suggest that the gas in the luminous X-ray sources is too dense to allow efficient Fermi acceleration of a relativistic particle population, and therefore in these systems the gas simply heats as it crosses the shock. Conversely, in the tenuous ADAF accretion flows studied here the relativistic particles are able to avoid thermalization due to the long collisional mean free path, resulting in the development of a significant nonthermal component in the particle distribution that powers the jets and produces the strong radio emission. We therefore conclude that the coupled, self-consistent theory for the disk structure and the particle acceleration investigated here provides a natural explanation for the outflows observed in many radio-loud systems containing black holes.

The authors are grateful to Lev Titarchuk for providing a number of useful comments on the manuscript, and also to the anonymous referee for several insightful suggestions that significantly improved the paper.

APPENDIX A TREATMENT OF THE VERTICAL STRUCTURE In principle, the pressure P, density , diffusion coefficient , Green's function fG , and velocity components vr and vz in the disk all display significant variations in the vertical (z) direction. Following Abramowicz & Chakrabarti (1990), we use the first five quantities to represent vertical averages over the disk structure at radius r. However, the vertical variation of the velocity component vz must be treated differently. Here, we assume for simplicity that the vertical expansion is homologous, and therefore the vertical velocity variation is given by vz (r; z) ? B(r)z: ðA1Þ


496

LE & BECKER

Vol. 632

It follows that the vertical velocity at the surface of the disk, z ? H (r), can be written as v z (r ; z)j
z?H

? B(r)H (r):

ðA2Þ

In a steady state situation, we can also express the vertical velocity at the disk surface using vz (r; z)j
z ?H

? vr

dH : dr

ðA3Þ

By combining the two previous expressions, we find that the function B(r) is given by B(r) ? v
r

d ln H : dr

ðA4Þ

This result will prove useful when we vertically integrate the transport equation. Note that in terms of B(r), we can write the divergence of the flow velocity v in cylindrical coordinates as :=v? 1@ @ vz 1 @ ? ðr v r Þ þ ðrvr Þ þ B(r); r @r r @r @z ðA5Þ

where we have assumed azimuthal symmetry. Application of equation (A4) now yields := v ? 1@ ðrH vr Þ: Hr @ r ðA6Þ

The steady state transport equation expressed in cylindrical coordinates is (see eq. [61]) ! N0 (E À E0 )(r À rà ) @ fG @ fG 1 1 @ d vz @ fG 1 @ @ fG r þ vz ? þ vr ðr v r Þ þ À A0 c(r À rà ) fG : ðA7Þ E þ 3 r @r @r @z dz @E r @r @r (4E0 )2 rà Hà R1 Operating on equation (A7) with 0 dz and applying equation (A1) yields, after partially integrating the term containing vz on the left-hand side, ! N0 (E À E0 )(r À rà ) @ 1 1d @ fG 1 @ @ fG rH þ ðrvr Þ þ B HE À A0 cHà (r À rà ) fG ; vr (HfG ) À HBfG ? þ @r 3 r dr @E r @r @r (4E0 )2 rà ðA8Þ

where the symbols fG , vr , and now refer to vertically averaged quantities. Using equations (A4), (A5), and (A6), we can rewrite the vertically integrated transport equation as N0 (E À E0 )(r À rà ) @ fG 1@ @ fG 1 @ @ fG rH ? þ H vr ðrH vr ÞE À A0 cHà (r À rà ) fG : ðA9Þ þ 3r @ r @r @E r @r @r (4E0 )2 rà This expression is used in x 5.1 to analyze the transport of the relativistic particles in the disk. APPENDIX B DERIVATION OF THE ESCAPE PARAMETER The dimensionless parameter A0 appearing in equation (58) determines the rate of particle escape through the surface of the disk due to random walks occurring near the shock location. Since the particles are accelerated as a consequence of collisions with magnetic waves, we assume that the thickness of the shock is comparable to the magnetic mean free path k mag . In order to estimate A0 ,wemodel the escape of the particles from the disk using the analogy of ``leakage'' through an opening in a cylindrical pipe with a radius equal to the half-thickness of the disk at the shock location, Hà . The length of the open section of the pipe is set equal to the shock thickness k mag . The particle number density in the open section is governed by the equation vx dnr nr ?À ; dx tesc ðB1Þ

where vx , nr ,and tesc represent the flow velocity, the relativistic particle number density, and the average time for the particles to escape through the open walls of the pipe via diffusion, respectively. Upon integration, the solution to equation ( B1) is given by x nr (x) ? n0 exp À ; vx tesc ðB2Þ


No. 1, 2005

OUTFLOWS FROM ADVECTION-DOMINATED DISKS

497

where n0 is the incident number density as the flow encounters the opening in the pipe, at x ? 0. We can approximate the solution for nr (x) by performing a Taylor expansion around x ? 0, which yields x nr ( x ) % n0 1 À : ðB3Þ vx tesc The fraction of particles that escape from the pipe can therefore be estimated by setting x ? k k mag nr ? : fesc ? 1 À n0 x?k mag vx tesc
mag

to obtain ðB4Þ

In order to make contact with the disk application, we note that according to equations (75) and (76), the fraction of particles that escape as the gas crosses the isothermal shock is given by fesc ? A0 c ; và ðB5Þ

where và (vþ þ vÀ )/2 is the mean velocity at the shock and we have assumed that advection dominates over diffusion. Eliminating fesc between equations ( B4) and ( B5) and setting vx ? và , we find that A0 ? k mag : ctesc
esc

ðB6Þ is related to k
mag

Within the context of our one-dimensional model for the particle transport in the disk, the mean escape time t and the disk half-thickness at the shock HÃ via tesc ?
2 HÃ HÃ ? ; vdiA ck mag

ðB7Þ

where vdiA ? ck mag /HÃ denotes the vertical diffusion velocity of the protons in the tangled magnetic field near the shock, which is valid provided HÃ /k mag > 1. Eliminating tesc between equations ( B6) and ( B7) then yields k mag 2 A0 ? < 1: ðB8Þ HÃ The diffusion coefficient at the shock is related to the magnetic mean free path by the standard expression (e.g., Reif 1965) ? and therefore equation ( B8) can be rewritten as A0 ? 3 Ã cHÃ 2 ; ðB10Þ ck mag ; 3 ðB9Þ

where à (À þ þ )/2 denotes the average of the upstream and downstream values of on either side of the shock.
REFERENCES Abramowicz, M. A., & Chakrabarti, S. K. 1990, ApJ, 350, 281 Chakrabarti, S. K. 1989a, PASJ, 41, 1145 Abramowitz, M., & Stegun, I. A. 1970, Handbook of Mathematical Functions ------. 1989b, ApJ, 347, 365 ( New York: Dover) ------. 1990, MNRAS, 243, 610 Aharonian, F. A. 2002, MNRAS, 332, 215 ------. 1996, ApJ, 464, 664 Anyakoha, M. W., Okeke, P. N., & Okoye, S. E. 1987, Ap&SS, 132, 65 Chakrabarti, S. K., Acharyya, K., & Molteni, D. 2004, A&A, 421, 1 Atoyan, A., & Dermer, C. D. 2004a, ApJ, 613, 151 Chakrabarti, S. K., & Das, S. 2004, MNRAS, 349, 649 ------. 2004b, ApJ, 617, L123 Chakrabarti, S. K., & Molteni, D. 1993, ApJ, 417, 671 Becker, P. A. 1992, ApJ, 397, 88 ------. 1995, MNRAS, 272, 80 Becker, P. A., & Kazanas, D. 2001, ApJ, 546, 429 Chen, X., Abramowicz, M. A., & Lasota, J. P. 1997, ApJ, 476, 61 Becker, P. A., & Le, T. 2003, ApJ, 588, 408 Das, S., Chattopadhyay, I., & Chakrabarti, S. K. 2001, ApJ, 557, 983 Becker, P. A., & Subramanian, P. 2005, ApJ, 622, 520 De Villiers, J. P., Hawley, J. F., Krolik, J. H, & Hirose, S. 2005, ApJ, 620, 878 Becker, P. A., Subramanian, P., & Kazanas, D. 2001, ApJ, 552, 209 Drury, L. O'C., & Volk, H. J. 1981, ApJ, 248, 344 Å Bicknell, G. V., & Begelman, M. C. 1996, ApJ, 467, 597 Ellison, D. C., & Eichler, D. 1984, ApJ, 286, 691 Biretta, J. A., Junor, W., & Livio, M. 2002, NewA Rev., 46, 239 Falcke, H., & Biermann, P. L. 1999, A&A, 342, 49 Biretta, J. A., Sparks, W. B., & Macchetto, F. 1999, ApJ, 520, 621 Felten, J. E. 1968, ApJ, 151, 861 Biretta, J. A., Stern, C. P., & Harris, D. E. 1991, AJ, 101, 1632 Felten, J. E., Arp, H. C., & Lynds, C. R. 1970, ApJ, 159, 415 Blandford, R. D., & Begelman, M. C. 1999, MNRAS, 303, L1 Fermi, E. 1954, ApJ, 119, 1 Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29 Ford, H. C., et al. 1994, ApJ, 435, L27 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 Frank, J., King, A. R., & Raine, D. J. 2002, Accretion Power in Astrophysics Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 (3rd ed.; Cambridge: Cambridge Univ. Press)


498

LE & BECKER
Pringle, J. E. 1981, ARA&A, 19, 137 Protheroe, R. J., & Kazanas, D. 1983, ApJ, 265, 620 Quataert, E. 2003, Astron. Nachr., 324, 435 Reif, F. 1965, Fundamentals of Statistical and Thermal Physics ( New York: McGraw-Hill) Reynolds, C. S., Matteo, T. D., Fabian, A. C., Hwang, U., & Canizares, C. R. 1996, MNRAS, 283, L111 Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics ( New York: Wiley) Schlickeiser, R. 1989a, ApJ, 336, 243 ------. 1989b, ApJ, 336, 264 Schodel, R., et al. 2002, Nature, 419, 694 Å Spruit, H. C. 1987, A&A, 184, 173 Subramanian, P., Becker, P. A., & Kafatos, M. 1996, ApJ, 469, 784 Subramanian, P., Becker, P. A., & Kazanas, D. 1999, ApJ, 523, 203 Webb, G. M., & Bogdan, T. J. 1987, ApJ, 320, 683 Weinberg, S. 1972, Gravitation and Cosmology ( New York: Wiley) Yang, R., & Kafatos, M. 1995, A&A, 295, 238 Yuan, F. 2000, MNRAS, 319, 1178 Yuan, F., Markoff, S., & Falcke, H. 2002, A&A, 383, 854

Hawley, J. F., Smarr, L. L., & Wilson, J. R. 1984a, ApJ, 277, 296 ------. 1984b, ApJS, 55, 211 Jones, F. C., & Ellison, D. C. 1991, Space Sci. Rev., 58, 259 Kafatos, M., & Yang, R. 1994, MNRAS, 268, 925 Kazanas, D., & Ellison, D. C. 1986, ApJ, 304, 178 Lanzafame, G., Molteni, D., & Chakrabarti, S. K. 1998, MNRAS, 299, 799 Le, T., & Becker, P. A. 2004, ApJ, 617, L25 Livio, M. 1999, Phys. Rep., 311, 225 Lovelace, R. V. E. 1976, Nature, 262, 649 Lu, J., Gu, W., & Yuan, F. 1999, ApJ, 523, 340 Lu, J., & Yuan, F. 1997, PASJ, 49, 525 ------. 1998, MNRAS, 295, 66 Mannheim, K., & Schlickeiser, R. 1994, A&A, 286, 983 McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977 Molteni, D., Sponholz, H., & Chakrabarti, S. K. 1996, ApJ, 457, 805 Narayan, R., Kato, S., & Honma, F. 1997, ApJ, 476, 49 Narayan, R., & Yi, I. 1994, ApJ, 428, L13 ------. 1995, ApJ, 444, 231 Owen, F. N., Eilek, J. A., & Kassim, N. E. 2000, ApJ, 543, 611 ? Paczynski, B., & Wiita, P. J. 1980, A&A, 88, 23