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

The Astrophysical Journal, 661:416 Y 429, 2007 May 20
# 2007. The American Astronomical Society. All rights reserved. Printed in U.S.A.

PARTICLE ACCELERATION IN ADVECTION-DOMINATED ACCRETION DISKS WITH SHOCKS: GREEN'S FUNCTION ENERGY DISTRIBUTION
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 2006 November 16; accepted 2007 January 6

ABSTRACT The distribution function describing the acceleration of relativistic particles in an advection-dominated accretion disk is analyzed using a transport formalism that includes first-order Fermi acceleration, advection, spatial diffusion, and the escape of particles through the upper and lower surfaces of the disk. When a centrifugally supported shock is present in the disk, the concentrated particle acceleration occurring in the vicinity of the shock channels a significant fraction of the binding energy of the accreting gas into a population of relativistic particles. These high-energy particles diffuse vertically through the disk and escape, carrying away both energy and entropy and allowing the remaining gas to accrete. The dynamical structure of the disk /shock system is computed self-consistently using a model previously developed by the authors that successfully accounts for the production of the observed relativistic outflows ( jets) in M87 and Sgr Aö . This ensures that the rate at which energy is carried away from the disk by the escaping relativistic particles is equal to the drop in the radial energy flux at the shock location, as is required for energy conservation. We investigate the influence of advection, diffusion, and acceleration on the particle distribution by computing the nonthermal Green's function, which displays a relatively flat power-law tail at high energies. We also obtain the energy distribution for the particles escaping from the disk, and we conclude by discussing the spectrum of the observable secondary radiation produced by the escaping particles. Subject headingg accretion, accretion disks -- black hole physics -- galaxies: jets -- hydrodynamics s: 1. INTRODUCTION Accretion flows onto supermassive black holes with mass M k 108 M are believed to power the high-energy emission observed from quasars and active galactic nuclei (AGNs), as discussed by Lynden-Bell (1969; see also Novikov & Thorne 1973; Rees 1984; Ford et al. 1994). When the accretion rate is relatively low compared with the Eddington value, the flow is radiatively inefficient, and the gas temperature approaches the virial value. In this case the disk has an advection-dominated structure, and most of the binding energy is swallowed by the black hole (e.g., Narayan et al. 1997; Blandford & Begelman 1999; Becker et al. 2001; Becker & Le 2003). Despite the inefficiency of X-ray production in these sources, advection-dominated disks are luminous radio and -ray emitters, and they frequently display powerful bipolar jets of matter escaping from the central mass, presumably containing highenergy particles (e.g., Sambruna et al. 2004; Di Matteo et al. 2000; Allen et al. 2000; Urry & Padovani 1995; Owen et al. 2000). 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
Also at: Department of Physics and Astronomy, George Mason University, Fairfax, VA 22030- 4444.
1

1991). In the present paper we consider the analogous process occurring in hot, advection-dominated accretion flows (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 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 radio and -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 precluded the treatment of acceleration in disks. The state of the theory was advanced significantly 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 directly 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. Accretion disks around black holes are not spherical, except perhaps in the innermost region, and therefore a transport equation written in cylindrical geometry is required in order to describe the acceleration of energetic particles in the accreting plasma. In this 416


PARTICLE ACCELERATION IN SHOCKED DISKS paper we develop a self-consistent cylindrically symmetric model describing the acceleration of protons and /or electrons via firstorder Fermi acceleration processes operating in a hot, advectiondominated accretion disk containing an isothermal shock. This is the third and final paper in a series studying the production of relativistic particles in ADAF disks with shocks. In Le & Becker (2004, hereafter Paper I ), we successfully applied our inflow/ outflow model to explain the observed kinetic power of the outflows in M87 and Sgr Aö , and in Le & Becker (2005, hereafter Paper II ), we discussed the detailed dynamics of ADAF disks containing isothermal shocks. On the basis of this dynamical model, we can self-consistently compute the structure of the shocked disk and also the rate of escape of particles and energy from the disk at the shock location. In this paper we provide further details regarding the distribution in space and energy of the accelerated particles in the disk, and we also describe the energy distribution of the relativistic particles escaping from the disk at the shock location. Our focus here is on the computation of the Green's function for monoenergetic particle injection , which describes in detail how the particles are advected, diffused, and accelerated within the disk, resulting in a nonthermal particle distribution with a high-energy power-law tail. We also discuss the spectrum of the observable secondary radiation produced by the escaping particles, including radio (synchrotron) as well as inverse Compton X-ray and -ray emission, which provides a quantitative basis for developing observational tests of the model. The first-order Fermi acceleration scenario that we focus on here (and in Papers I and II ) parallels the early studies of cosmicray acceleration in supernova shock waves. In analogy with the supernova case, we expect that shock acceleration in the disk will naturally produce a power-law energy distribution for the accelerated particles. Our earlier work in Paper II established that the first-order Fermi mechanism provides a very efficient means for accelerating the jet particles. In this paper, the solution for the Green's function, fG (E ; r), describing the evolution of monoenergetic seed particles injected into the disk, is obtained by solving the transport equation using the method of separation of variables. The eigenvalues and the corresponding spatial eigenfunctions are determined using a numerical approach because the differential equation obtained for the spatial function cannot be solved in closed form. The resulting solution for the Green's function provides useful physical insights into how the relativistic particles are distributed in space and energy. We can also use the solution for fG (E ; r) to check the self-consistency of the theory by comparing the numerical values obtained for the number and energy densities by integrating the Green's function with those obtained in Paper II by directly solving the associated differential equations for these quantities. The remainder of the paper is organized as follows. In x 2we briefly review the inviscid dynamical model used to describe the structure of the ADAF disk /shock system. In x 3we develop and solve the associated steady state particle transport equation to obtain the solution for the Green's function, fG (E ; r), describing the evolution of monoenergetic seed particles injected from a source located at the shock radius, and in x 4 we evaluate the Green's function using dynamical parameters appropriate for M87 and Sgr Aö . In x 5 we discuss the observational implications of our results for the production of secondary radiation by the outflowing relativistic particles, and in x 6 we review our main conclusions. 2. ACCRETION DYNAMICS It has been known for some time that inviscid accretion disks can display both shocked and shock-free (i.e., smooth) solutions,

417

Fig. 1.-- Schematic representation of the disk /shock /outflow model. The filled circles in the disk represent the accelerated particles, and the open circles represent the MHD scattering centers, which advect with the background flow velocity. The seed particles are injected at the shock location.

depending on the values of the energy and angular momentum per unit mass in the gas supplied at a large radius (e.g., Chakrabarti 1989; 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 shockfree 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). The possible presence of shocks in hot ADAF disks with low viscosity motivated us to explore in Paper II the relationship between the dynamical structure of the disk /shock system and the possible acceleration of relativistic particles that can escape from the disk to power the outflows commonly observed from radioloud accretion flows around both stellar-mass and supermassive black holes. The disk /shock /outflow model 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 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 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.


418 2.1. Transport Rates

LE & BECKER

Vol. 661

Becker & Le (2003) and Becker & Subramanian (2005) demonstrated that three integrals of the flow are conserved in viscous ADAF disks when outflows are not occurring. The three integrals correspond to the mass transport rate, M ¼ 4rH u; the angular momentum transport rate, J ¼ Mr 2 þ G; and the energy transport rate, 1 v2 × 1 u2 × P × U × õ ; E ¼ þG × M 2 2 Ï 2÷ Ï 1÷

analytically, while retaining reasonable agreement with the dynamics of actual disks with low viscosities. We utilize the set of physical conservation equations employed by Chakrabarti (1989) and Abramowicz & Chakrabarti (1990), who investigated the structure of one-dimensional, inviscid, steady state, axisymmetric accretion flows. In the inviscid limit, the torque G vanishes, and the angular momentum per unit mass transported through the disk is given by l (r) J ¼ r 2 (r) ¼ constant: M Ï 9÷

We can combine equations (3) and (7) to show that the energy transported per unit mass in the inviscid case can be written as Ï 3÷ E 1 1 l2 c2 ¼ u2 × × s × õ; 2 2 r2 þ1 M Ï10÷

where is the mass density, is the angular velocity, G is the torque, H is the disk half-thickness, u is the positive radial inflow speed, v ¼ r is the azimuthal velocity, U is the internal energy density, P ¼ ( þ 1)U is the pressure, and õ (r ) þGM r þ rS Ï 4÷

denotes the pseudo-Newtonian gravitational potential for a black hole with mass M and Schwarzschild radius rS ¼ 2GM /c 2 ´ ( Paczynsky & Wiita 1980). Each of the various densities and velocities 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 quantities M and J remain constant throughout the entire flow, and therefore they represent the rates at which mass and angular momentum, respectively, enter the black hole. The energy transport rate E also remains constant, except at the location of an isothermal shock if one is present in the disk. 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 Ï 5÷

wherewehavealsousedthe relation l ¼ rv . The resulting disk / shock model depends on three free parameters, namely, the energy transport per unit mass, , the specific heats ratio, , and the specific angular momentum, l . The value of 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). 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. The condition þ > × must be satisfied as a consequence of the loss of energy through the upper and lower surfaces of the disk at the shock location. In order to satisfy global energy conservation, the rate at which energy escapes from the gas as it crosses the shock must be equal to the jet luminosity, Ljet , and therefore we have L
jet

þà E ¼ þM à / ergs sþ1 ;

Ï11÷

where is the kinematic viscosity. The disk half-thickness H is given by the standard hydrostatic prescription H (r ) ¼ cs ; K Ï 6÷

where cs represents the adiabatic sound speed, defined by =2 P 1 cs ( r ) ; Ï 7÷ and K denotes the Keplerian angular velocity of matter in a circular orbit at radius r in the pseudo-Newtonian potential (eq. [4]), defined by
2 K (r)

where à denotes the difference between quantities on the downstream and upstream sides of the shock, so that à ¼ × þ þ < 0. The drop in at the shock has the effect of altering the transonic structure of the flow in the postshock region, compared with the shock-free configuration (see Paper II ). In the inviscid case, the flow is purely adiabatic except at the location of an isothermal shock (if one is present), and therefore the pressure is related to the density profile via P ¼ D0 ; Ï12÷

GM 1 dõ : ¼ 2 r dr r ( r þ rS )

Ï 8÷

2.2. Inviscid Flow Equations Following the approach taken in Paper II, we focus here on the inviscid case, since that allows the disk structure to be computed

where D0 is a parameter related to the specific entropy that remains constant except at the location of a shock. The process of constructing a dynamical model for a specific disk /shock system begins with the selection of values for the fundamental parameters þ , l , and . For given values of these three quantities, a unique flow structure can be obtained using the procedure discussed in xx 3 and 4 of Paper II, which yields a numerical solution for the inflow speed u (r). Following Narayan et al. (1997 ), we assume an approximate equipartition between the gas and magnetic pressures, and therefore we set ¼ 1:5. The range of acceptable values for þ and l is constrained by the observations of a specific object, as discussed in x 7 of Paper II.


No. 1, 2007

PARTICLE ACCELERATION IN SHOCKED DISKS

419

3. STEADY STATE PARTICLE TRANSPORT The gas in ADAF disks is hot and tenuous and therefore essentially collisionless. In this situation, the relativistic particles diffuse mainly via collisions with magnetohydrodynamic ( MHD) ´ scattering centers (e.g., Alfven waves), which control both the spatial transport and the acceleration of the particles. The magnetic scattering centers are advected along with the background flow, and therefore the shock width is expected to be comparable to the magnetic coherence length, k mag . Particles from the highenergy 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 distribution in the disk. 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 power-law energy distribution, which is a general characteristic of Fermi processes ( Fermi 1954). Two effects limit the maximum energy 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 will eventually exceed the thickness of the disk as the particle energy is increased, resulting in escape from the disk without further acceleration. As a consequence of the random walk, 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 (see Fig. 1). Our goal in this paper is to analyze the Green's function fG (E ; r) in the disk on the basis of the transport equation and the dynamical results discussed in Paper II. We focus specifically on models 2 and 5 from Paper II, which describe the disk structures for M87 and Sgr Aö , respectively. 3.1. Transport Equation The particle transport formalism we use includes advection, spatial diffusion, first-order Fermi acceleration, and particle escape. Our treatment of the Fermi process includes both the general compression related to the overall convergence of the accretion flow and the enhanced compression that occurs at the shock. To avoid unnecessary mathematical complexity, we use a simplified model for the spatial transport in which only the radial (r)component 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 will adopt the test particle approximation used in Paper II, meaning that the dynamical effect of the relativistic particle pressure on the flow structure is assumed to be negligible. This assumption is valid provided that the pressure of the relativistic particles is much smaller than the background (thermal ) pressure, which is verified for models 2 and 5 in x 8 of Paper II. We also assume that the injection of the seed particles and the escape of the accelerated particles from the disk are both localized at the isothermal shock radius, which is consistent with the strong concentration of first-order Fermi acceleration in the vicinity of the shock. The localization of the particle injection and escape at the shock radius allows us to make a direct connection between the jump in the energy flux in the disk at the shock location and the energy carried away by the escaping particles. This connection is essential in order to maintain self-consistency between the dynamical and transport calculations, as discussed in Paper II. In a steady state situation, the Green's function, fG (E; r), representing the particle distribution resulting from the continual injection of mono-

energetic seed particles from a source located at the shock radius (r ¼ rö ) satisfies the transport equation ( Becker 1992) à @ fG 1 @þ3 E v =:fG × fsource þ fesc ; ¼ 0 ¼þ:= F þ 3E 2 @ E @t where the specific flux F is evaluated using F ¼ þ:fG þ v E @ fG ; 3 @E Ï14÷ Ï13÷

and the quantities E, , and v represent the particle energy, the spatial diffusion coefficient, and the vector velocity, respectively. ^ ^ ^ The velocity has components given by v ¼ vr r × vz z × v f, with vr ¼ þu < 0, since we are dealing with inflow. The injection of the monoenergetic seed particles is described by the function fso ¼ N0 (E þ E0 )(r þ rö ) ; (4E0 ) 2 rö Hö Ï15÷

urce

where E0 is the energy of the injected particles, N0 is the particle injection rate, and Hö H (rö ) represents the half-thickness of the disk at the shock location. The escape of the particles through the upper and lower surfaces of the disk is represented by the term fesc ¼ A0 c(r þ rö ) fG ; Ï16÷

where the dimensionless parameter A0 is computed using (see Appendix B of Paper II ) 3 ö A0 ¼ cHö 2 < 1; Ï17÷

and ö (þ × × ) /2 denotes the mean value of the spatial diffusion coefficient at the shock radius. The subscripts ``þ'' and ``+'' will be used to denote quantities measured just upstream and just downstream from the shock, respectively. The condition A0 < 1 is required for the validity of the diffusive treatment of the vertical escape employed in our approach. In order to compute the Green's function fG (E; r), we must specify the injection energy of the seed particles, E0 ,aswellas their injection rate, N0 . Following the approach taken in Paper II, 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 , is several times larger than the mean ion ther0 = mal 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 local Maxwellian distribution. With E0 specified, we can compute the particle injection rate N0 using the energy conservation condition N0 E0 ¼ þM à / ergs sþ1 ; Ï18÷

where à < 0 represents the jump in the energy transport per unit mass as the plasma crosses the shock (see eqs. [10] and [11]). This relation ensures that the energy injection rate for the seed particles is equal to the energy loss rate for the background gas at the isothermal shock location.


420

LE & BECKER

Vol. 661

The total number and energy densities of the relativistic particles, denoted by nr and Ur , respectively, are related to the Green's function via Z nr ( r ) ¼ Z Ur (r) ¼
0 0 1 1

where uþ and u× denote the positive inflow speeds just upstream and downstream from the shock, respectively. The velocity jump condition appropriate for an isothermal shock is (see Paper II ) u× 1 ¼ ; uþ M2 þ Ï24÷

4E 2 fG (E ; r) dE ; 4E3 fG (E; r) dE ; Ï19÷

which determine the normalization of fG . Equations (13) and (14) can be combined to obtain the alternative form v = :fG ¼ E @ fG := v × := Ï:fG ÷ × fsource þ fesc ; 3 @E Ï20÷

where the left-hand side represents the comoving (advective) time derivative and the terms on the right-hand side describe firstorder Fermi acceleration, spatial diffusion, the particle source, and the escape of particles from the disk at the shock location, respectively. Note that the escape and particle injection are localized to the shock radius due to the presence of the -functions in equations (15) and (16). Our focus here is on the first-order Fermi acceleration of relativistic particles at a standing shock in an accretion disk, and therefore equation (20) does not include secondorder Fermi processes that may also occur in the flow due to MHD turbulence (e.g., Schlickeiser 1989a, 1989b; Subramanian et al. 1999; Becker et al. 2006). Since no energy loss mechanisms are included in the model, all of the particles injected with energy E0 will experience an increase in energy. Under the assumption of cylindrical symmetry, equations (15), (16), and (20) can be combined to obtain v @ fG @ fG E @ fG × vz þ @r @z 3 @E N0 (E ¼ (4 ! 1d d vz 1@ @ fG (r v r ) × r þ r dr r @r dz @r þ E 0 ) (r þ rö ) þ A0 c(r þ rö ) fG ; E0 )2 rö Hö Ï21÷ where the escape of particles from the disk is described by the final term on the right-hand side. We focus here on the vertically integrated form of the transport equation, which can be written as (see Paper II, Appendix A) @ fG 1d E @ fG 1 @ @ fG (rHu) rH ¼þ × þHu r dr 3 @E r @r @r @r 0 (E þ E0 ) (r þ rö ) N þ A0 cHö (r þ rö ) fG ; × (4E0 )2 rö Ï22÷ where u ¼ þvr and fG and are vertically averaged quantities. Due to the presence of the velocity derivative on the right-hand side of equation (22), the first-order Fermi acceleration of the particles is most pronounced near the shock, where u is discontinuous. In the vicinity of the shock, we find that du ! (uþ þ u× ) (r þ rö ); dr r ! rö ; Ï23÷

where Mþ is the incident Mach number on the upstream side of the shock. Although the Fermi acceleration of the particles is concentrated at the shock, the rest of the disk also contributes to the particle acceleration because of the general convergence of the MHD waves in the accretion flow. In order to close the system of equations and solve for the Green's function using equation (22), 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, which we have demonstrated in Paper II. 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 utilize the general form (see Paper II, x 5) 2 r þ1 ; Ï25÷ ( r ) ¼ 0 u( r ) r S rS where 0 > 0 is a dimensionless constant that can be computed for a given source on the basis of energy considerations (see x 4). Due to the appearance of the radial speed u in equation (25), 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. We also point out that the dimensionless escape parameter A0 can be evaluated by combining equations (17) and (25), which yields 4 30 uö rS 2 rö þ 1 < 1; Ï26÷ A0 ¼ cHö rS where uö (uþ × u× ) /2 denotes the mean value of the inflow speed at the shock location. Note that the value of the diffusion parameter 0 is constrained by the inequality in equation (26). 3.2. Separation of Variables For values of the particle energy E > E0, the source term in equation (22) vanishes, and the remaining equation is separable in energy and space using the functions kn Eþ f n (E ; r ) ¼ Yn (r); Ï27÷ E0 where kn are the eigenvalues, and the spatial eigenfunctions Yn (r) satisfy the second-order ordinary differential equation þ Hu dYn ¼ dr

r

kn d 1d dYn (rHu)Yn × rH þ A0 cHö (r þ rö )Yn : r dr 3r dr dr

Ï28÷

The eigenvalues kn are determined by applying suitable boundary conditions to the spatial eigenfunctions, as discussed below. Once a numerical solution for the inflow speed u(r) has been obtained


No. 1, 2007

PARTICLE ACCELERATION IN SHOCKED DISKS

421

following the procedures described in Paper II, we can compute H (r)and (r) using equations (6) and (25), respectively. Since the coefficients in equation (28) cannot be represented in closed form, analytical solutions cannot be obtained, and therefore the eigenfunctions Yn (r) must be determined numerically. Away from the shock (r ¼ rö ), equation (28) reduces to 6 ! d 2 Yn d ln (rH ) u dYn kn uYn d ln (rHu) × ¼ 0; × × dr dr dr dr 2 3 which can be rewritten as ! d 2 Yn d ln (rHu) 2 rS dYn × × × dr r þ rS 0 (r þ rS ) 2 dr dr 2 kn rS Yn d ln (rHu) ¼ 0; × dr 3 0 ( r þ r S ) 2 Ï29÷

Fig. 2.-- Eigenvalues for model 2 ( M87; squares) and model 5 (Sgr Aö ; circles).

Ï30÷ terms in equation (30), we find that the global solutions for the spatial eigenfunctions can be written as ( Yn (r) ¼
in Gn (r);

where we have used equation (25) to substitute for . Once a solution for the inflow speed u(r) has been determined using the procedure described in x 2, the variation of the disk halfthickness H (r) can be computed by combining equations (1), (6), (7), and (12). The result obtained is 1 r K (r )
3=2

r

rö ;

an G

out n

( r ) ; r ! rö ;

Ï34÷

H (r) ¼

( r þ r S ) u( r ) K

!

Ï1þ ÷=Ï1× ÷

;

Ï31÷

in ou where an is a matching constant and Gn (r)and Gn t (r)denote the fundamental solutions to equation (30) possessing the asymptotic behaviors

where K is given by equation (8) and K is the ``entropy parameter,'' which remains constant except at the location of an isothermal shock if one is present in the flow (see Paper II; Becker & Le 2003). By using the solutions obtained for u (r)and H (r), we are able to compute all of the coefficients in equation (30), which allows us to solve numerically for the spatial eigenfunctions Yn (r). 3.3. Eigenvalues and Eigenfunctions The global solution for the eigenfunction Yn (r) must satisfy the continuity and derivative jump conditions associated with the presence of the shock /source at radius r ¼ rö . By integrating equation (28) with respect to the radius in the vicinity of the shock, we find that à(Yn ) ¼ 0; kn dYn à uYn × ¼þ cA0 Yn (rö ); 3 dr Ï32÷ Ï33÷

in in Gn (r) ! gn (r) ou Gn t

r þ1 rS 1 rþ out (r) ! gn (r) ; rS



þ

kn =Ï3 ×3÷

; r ! 1;

r ! rS ; Ï35÷

at small and large radii, respectively. The value of an is determined using the continuity condition (see eq. [32]), which yields an ¼
in Gn (rö ) : out (r ) Gn ö

Ï36÷

where the symbol à denotes the difference between postshock and preshock quantities, and we have made use of the fact that H is continuous across the isothermal shock. Equations (32) and (33) establish that Yn (r) must be continuous at the shock location, and its derivative must display a jump there. Since the spatial eigenfunctions Yn (r) satisfy the second-order linear differential equation given in equation (30), we must also impose two boundary conditions in order to determine the global solutions and the associated eigenvalues kn . The required conditions can be developed on the basis of our knowledge of the dynamical structure of the disk near the event horizon (r ¼ rS ) and at large radii (see Becker & Le 2003 and x 4of Paper II ). By using the Frobenius expansion method to identify the dominant

The validity of the asymptotic forms in equations (35) is verified in x 4.1 via comparison with the numerical solutions obtained for the spatial eigenfunctions. We use a bidirectional integration scheme to solve for the eigenvalues kn on the basis of the boundary conditions in equations (35). The value of kn is iterated using a bisection method until the Wronskian of the inner and outer solutions vanishes at a matching radius located in the postshock region. Once a particular eigenvalue has been determined to the required accuracy, the matching constant an is computed using equation (36 ). This procedure is repeated until the desired number of eigenvalues and eigenfunctions have been obtained. The sequences of eigenvalues associated with the model 2 ( M87 ) and model 5 (Sgr Aö ) parameters are plotted in Figure 2. Note that k1 $ 4 in all cases, indicating that the acceleration is quite efficient, in analogy with the case of cosmicray acceleration (e.g., Blandford & Ostriker 1978). The first two eigenvalues in models 2 and 5 are k1 ¼ 4:165 and k2 ¼ 6:415 and k1 ¼ 4:180 and k2 ¼ 6:344, respectively. It follows that at high energies the power-law slope of the energy distribution is dominated by the first eigenvalue, k1 , since all of the subsequent eigenvalues are much larger (see Fig. 2).


422 3.4. Eigenfunction Orthogonality

LE & BECKER

Vol. 661

We can establish several useful general properties of the eigenfunctions by rewriting equation (30) in the Sturm-Liouville form, ! d dYn S (r ) × kn ! (r)Yn (r) ¼ 0; dr dr where ( " 1 1 #) þ þ rH 1 rö r S (r ) exp þ1 þ þ1 rö H ö ö 0 rS rS Ï38÷ and the weight function ! (r) is defined by ! (r ) uS d ln (rHu) : 3 dr Ï39÷ Ï37÷

Multiplying each side of this expression by the product Yn (r)!(r) and integrating from r ¼ rS to r ¼ 1 yields Z
1

fG (E0 ; r)Yn (r) ! (r) dr ¼
rS

Nmax X m ¼1

Z b
m r

1

Ym (r)Yn (r) ! (r) dr:
S

Ï44÷

On the basis of the orthogonality of the eigenfunctions, as expressed by equation (41), we note that only the m ¼ n term survives on the right-hand side of equation (44), and we therefore obtain Z1 Z1 2 fG (E0 ; r)Yn (r) ! (r) dr ¼ bn Yn (r) ! (r) dr: Ï45÷
rS rS

Solving this expression for the expansion coefficient bn yields R1 r fG (E0 ; r )Yn (r)!(r ) dr bn ¼ S ; Ï46÷ In where the quadratic normalization integrals, I n , are defined by Z1 2 In Yn (r) ! (r) dr: Ï47÷
r
S

Note that ! (r) displays a -function discontinuity at the shock through its dependence on the derivative of u (r). In the vicinity of the shock, we can combine equations (23), (38), and (39) to show that uþ þ u !( r ) ! 3 ö
×

(r þ rö );

r ! rö :

Ï40÷

On the basis of the boundary conditions satisfied by the spatial eigenfunctions Yn (r) (see eqs. [35]), we demonstrate in the Appendix that the eigenfunctions satisfy the orthogonality relation Z
r 1

Yn (r)Ym (r)!(r) dr ¼ 0;
S

m ¼ n: 6

Ï41÷

To complete the calculation of the expansion coefficients, we need to evaluate the distribution function at the source energy, fG (E0 ; r). This can be accomplished by using equation (23) to substitute for the velocity derivative in equation (22) and then integrating with respect to E in a small range around the injection energy E0, which yields 8 > 3N0 < ; r ¼ rö ; 3 Ï48÷ fG (E0 ; r) ¼ (4)2 E0 rö Hö (uþ þ u× ) > : 0; r 6¼ rö : Although the source term in the transport equation (21) is a -function, we note that the Green's function evaluated at the injection energy E0 is not a -function. Instead, as indicated by equation (48), we find that the Green's function at the injection energy is zero for radii away from the shock, and it is equal to a finite value at the shock (r ¼ rö ). From a physical point of view, this reflects the fact that the particle injection occurs at the shock location in our model. The strong acceleration experienced by the particles at the shock eliminates the appearance of a -function behavior in the resulting Green's function. It is interesting to compare this with the behavior observed in the case of particle acceleration in a plane-parallel supernova shock wave, which was analyzed in detail by Blandford & Ostriker (1978). Our results are similar to theirs, except that the Blandford & Ostriker Green's function evaluated at the injection energy E0 has a finite value at all locations in the flow, because there is no convergence in the plasma surrounding the plane-parallel supernova shock wave. However, in our model, the plasma converges at all radii in the accretion flow (not just at the shock), and this additional particle acceleration causes the Green's function to vanish at the injection energy for all r 6¼ rö . Utilizing equation (48) to substitute for fG (E0 ; r) in equation (46) and carrying out the integration, we obtain bn ¼ N0 Yn (rö ) 23 E0 rö Hö ö I n ; Ï49÷

3.5. Eigenfunction Expansion We have established that Yn (r) is the solution to a standard Sturm-Liouville problem, and therefore it follows that the set of eigenfunctions is complete. Consequently, the Green's function fG (E; r) can be represented using the eigenfunction expansion (see eq. [27]) fG (E; r) ¼
Nmax X n¼1

bn Yn (r)

E E0

þ

k

n

;

E ! E0 ;

Ï42÷

where bn are the expansion coefficients (with either positive or negative signs) and the value of Nmax is established by analyzing the term-by-term convergence of the series. We use the value Nmax ¼ 10 in our numerical examples, which generally yields an accuracy of at least three decimal digits. Note that fG (E ; r) ¼ 0 for all E < E0 because no deceleration processes are included in the particle transport model. The expansion coefficients b1, b2, b3, : : : can be calculated by utilizing the orthogonality of the spatial eigenfunctions. At the source energy, E ¼ E0 , equation (42) reduces to fG (E0 ; r) ¼
Nmax X m ¼1

bm Ym (r):

Ï43÷

(4)


No. 1, 2007

PARTICLE ACCELERATION IN SHOCKED DISKS
TABLE 1 Disk Structur e P arameters M (M) 3.0 ; 10 2.6 ; 106
9

423

Source M87 ................... Sgr Aö ...............

Model 2 5

M (M yrþ1) 1.3 ; 10 8.8 ; 10
þ1 þ7

Ljet (ergs sþ1) 5.5 ; 10 5.0 ; 1038
43

þ 1.527 ; 10 1.229 ; 10þ
þ3 3

l 3.1340 3.1524

+ þ5.746 ; 10 þ8.749 ; 10þ
þ3 3

rö 21.654 15.583

uö =c 0.108 0.124

Note.--All quantities are expressed in gravitational units (GM ¼ c ¼ 1), except as indicated.

where we have made use of the fact that the weight function ! (r)has a -function behavior close to the shock (see eq. [40]). The singular nature of the weight function also needs to be considered when computing the normalization integrals I n defined in equation (47 ). By combining equations (40) and (47 ), we find that Z r ö þ I n ¼ lim ! (r)Yn2 (r) dr !0 rS Z1 uþ þ u× 2 ! (r)Yn (r) dr × × Yn2 (rö ); Ï50÷ 3 ö r ö × and we use this expression to evaluate the normalization integrals. 4. ASTROPHYSICAL APPLICATIONS In Paper II we investigated the acceleration of particles in an inviscid ADAF disk containing an isothermal shock. For a given source with a known black hole mass M,accretion rate M , and jet kinetic power Ljet , we found that a family of flow solutions can be obtained, with each solution corresponding to a different value of the diffusion parameter 0 (see eq. [25]). We concluded that the upstream energy flux, þ , is maximized for a certain value of 0 , and we adopted that particular value for 0 in our models for M87 and Sgr Aö in Paper II. In Table 1 we list the values of the pa rameters M, M , Ljet , þ , l, × , rö ,and uö associated with models 2 and 5, which describe M87 and Sgr Aö , respectively. 4.1. Numerical Solutions for the Eigenfunctions In order to compute the Green's function fG (E ; r), we must first solve for the eigenvalues kn and eigenfunctions Yn (r) following the procedures described in xx 3.3 and 3.4, based on equation (30).

Once the eigenvalues and eigenfunctions have been determined, we can compute the expansion coefficients bn using equation (49), after which fG (E ; r) is evaluated using equation (42). In Figure 3 in we compare the fundamental asymptotic solutions Gn (r) and out in out Gn (r) with the asymptotic solutions gn (r)and gn (r)for model 2 with n ¼ 1 (see eqs. [35]). The corresponding results for model 5 are plotted in Figure 4. The excellent agreement between the G(r)and g (r) functions confirms the validity of the asymptotic relations we have employed near the event horizon and also at large radii. Figures 3 and 4 also include plots of the global solutions for the first eigenfunction Y1 (r) associated with models 2 and 5, respectively (see eq. [34] ). Note that the global solutions display a derivative jump at the shock location, as is required according to equation (33). 4.2. Green's Function Particle Distribution We can combine our results for the eigenvalues, eigenfunctions, and expansion coefficients to evaluate the Green's function fG (E ; r) using equation (42). In Figures 5a and 5b we plot fG (E; r) as a function of the particle energy E at various radii r in the disk for models 2 and 5, respectively. Note that at the injection energy (E ¼ E0 ), the Green's function is equal to zero except at the source/shock radius, r ¼ rö , in agreement with equation (48). This reflects the fact that the particles are rapidly accelerated after they are injected, and therefore all particles have energy E > E0 once they propagate away from the shock location. The positive slopes observed at low energies when r 6¼ rö result from the summation of terms with either positive or negative expansion coefficients bn (see x 3.5). At energies above the turnover, the spectrum has a power-law shape determined by the first eigenvalue, k1 $ 4 (see Fig. 2). The relatively flat slope of the Green's function above the turnover reflects the high efficiency of the particle acceleration in the shocked disk. By contrast,

in ou in Fig. 3.-- Fundamental solutions to eq. (28), G1 (r)and G1 t (r) (solid lines), for model 2 ( M87 ), compared with the corresponding asymptotic solutions g1 (r)and gout (r)(dashed lines)in panels a and b, respectively (see eqs. [35]). The associated global solution for the first eigenfunction Y1 (r) (eq. [34]) is plotted in panel c.The 1 shock location at r ¼ rö is indicated.


424

LE & BECKER

Vol. 661

Fig. 4.-- Same as Fig. 3, but for model 5, which describes Sgr Aö .

in situations involving weak acceleration, the Green's function has a strong peak at the injection energy, surrounded by steep wings (e.g., Titarchuk & Zannias 1998). In Paper II we established that most of the injected particles are advected inward by the accreting MHD waves, eventually crossing the event horizon into the black hole. Hence, only a small fraction of the particles are able to diffuse upstream to larger radii, which is confirmed by the strong attenuation of the particle spectrum with increasing r that is apparent in Figure 5. The particles in the inner region (r < rö ) exhibit the greatest overall energy gain because they are able to cross the shock multiple times, and they also experience the strong compression of the flow near the event horizon. This is consistent with the results for the mean energy distribution, h E i Ur /nr , plotted in Figure 9 of Paper II. 4.3. Number and Energy Density Distributions In our numerical examples, we have generated results for the Green's function fG (E ; r) by summing the first 10 eigenfunctions

using equation (42). Once the Green's function energy distribution has been determined at a given radius, we can integrate it with respect to the particle energy E to obtain the corresponding values for the number and energy densities. By combining equations (19) and (42) and integrating term by term, we obtain nG ( r ) 4 E r UrG (r) 4E
3 0 Nmax X bn Yn (r) ; k þ3 n ¼1 n Nmax X bn Y n ( r ) : k þ4 n ¼1 n

4 0

Ï51÷

Equations (51) provide an important tool for checking the selfconsistency of our entire formalism. This is accomplished by comparing the results for the number and energy densities computed using equations (51) with the corresponding values obtained by directly solving the differential equations for nr and Ur using the method developed in Paper II. In principle, the two sets of results should agree closely if our procedure for computing

Fig. 5.-- Results for the relativistic particle Green's function fG (E; r) in units of ergsþ3 cmþ3, computed using eq. (42) for (a) M87 (model 2) and (b) Sgr Aö (model 5). The value of the radius r in units of GM /c 2 is indicated for each curve, with r ¼ rö denoting the shock location.


No. 1, 2007

PARTICLE ACCELERATION IN SHOCKED DISKS

425

Fig. 6.-- Plots of (a) the relativistic particle number density and (b) the relativistic particle energy density in cgs units for model 2, which describes M87. The solid lines represent the solutions computed by numerically integrating the associated differential equations, and the filled circles represent the corresponding results obtained by integrating the Green's function using eqs. (51). The location of the shock at r ¼ rö is indicated. Note the close agreement between the results, which confirms the accuracy of our computational method.

fG (E; r) is robust. We compare the model 2 results obtained for nr and Ur in Paper II with those computed using equations (51) in Figure 6. The corresponding comparison for model 5 is carried out in Figure 7. Note the close agreement between the two sets of results, which confirms the validity of the analysis involved in searching for the eigenvalues, computing the eigenfunctions, solving for the expansion coefficients, and calculating the Green's function. Moreover, the results in Figures 6 and 7 indicate that equation (42) converges successfully using the first 10 terms (Nmax ¼ 10) in the series. 4.4. Escaping Particle Distribution The Green's function fG (E; r) computed using equation (42) represents the energy distribution of the relativistic particles inside

the disk at radius r. However, our primary goal in this paper is the calculation of the energy distribution of the relativistic particles escaping from the disk at the shock radius, rö . In order to compute the energy spectrum of the escaping particles, we need to integrate equation (16) over the volume of the disk, which yields N
esc E

¼ (4E )2 rö Hö cA0 fG (E ; rö );

Ï52÷

esc where NE dE represents the number of particles escaping from the disk per unit time with energies between E and E × dE.At the highest particle energies, the escaping particle number spectrum es has a power-law shape, with NE c / E þ ,where ¼ k1 þ 2and k1 is the first eigenvalue (see eqs. [42] and [52]). The total number

Fig. 7.-- Same as Fig. 6, but for model 5, which describes Sgr Aö .


426

LE & BECKER

Vol. 661

esc esc Fig. 8.-- Plots of (a) the number distribution NE and (b) the energy distribution E NE for the relativistic particles escaping at the shock location, r ¼ rö (see eq. [52]). The solid curves represent model 2 ( M87 ), and the dashed curves represent model 5 (Sgr Aö ).

of particles escaping from the disk per second, Nesc , is related to es the spectrum NE c via Z1 es NE c dE ¼ 4rö Hö cA0 nö ; Nesc ¼ Ï53÷
0

integral in equation (53) are in excellent agreement with the values listed in Table 2, which were taken from Paper II. These tests confirm the validity of the analytical approach we have utilized in our derivation of the Green's function, given by equation (42). 5. OBSERVATIONAL IMPLICATIONS The jets of relativistic particles produced by the first-order Fermi shock acceleration mechanism considered here will flow outward from the disk with initially high pressure and low velocity. The rapid expansion of the plasma will lead to strong acceleration as the gas cools and the random internal energy is converted into ordered kinetic energy (see Paper II ). The terminal Lorentz factor, þ1 , associated with the M87 (model 2) and Sgr Aö (model 5) jets considered here is found to be 7.92 and 7.26, respectively, on the basis of comparison with the observed jet kinetic luminosities. However, robust confirmation of the model requires a detailed comparison between the predicted emission and the observed multifrequency spectra for a variety of sources. The jets envisioned here are expected to be charge neutral, with the protons ``dragging'' along an equal number of electrons. The electrons are expected to cool rapidly via inverse Compton and bremsstrahlung emission produced while the particles are still close to the active nucleus. On the other hand, the relativistic protons will cool inefficiently via direct radiation, and they also lose very little energy to the electrons via Coulomb coupling due to the low gas density. Hence, it is expected that the ions will maintain their high energies until they collide with particles or radiation in the ambient medium. Detailed consideration of this

where nö nr (rö ) is the relativistic particle number density at the shock location. The corresponding result for the total energy escape rate is given by Z1 esc NE EdE ¼ 4rö Hö cA0 Uö ; Ï54÷ Lesc
0

where Uö Ur (rö ) is the relativistic particle energy density at the shock location. The escaping particle distributions for M87 (model 2) and Sgr Aö (model 5) are plotted in Figure 8. The results obtained in the case of M87 are about 5 orders of magnitude larger than those associated with Sgr Aö due to the corresponding difference in the reported jet luminosity Ljet for these two sources (see Table 1). In Table 2 we list the values obtained for N0, Nesc , nö , Uö , ö , 0, A0 , Eesc ,and þ1 in models 2 and 5, where Eesc Uö /nö denotes the mean energy of the particles escaping from the disk to power the jet and þ1 represents the terminal Lorentz factor of the jet. We note that the values obtained for Lesc using the integral in equation (54) are in excellent agreement with the values for Ljet reported in Paper II and listed in Table 1. The agreement between Lesc and Ljet confirms that our model satisfies global energy con servation. We also find that the values for Nesc obtained using the

TABLE 2 Transport Equation P arameters N0 (sþ1) 2.75 ; 10 2.51 ; 1041
46

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

N

esc

/N0

nö (cmþ3) 2.01 ; 10 4.33 ; 105
4

Uö (ergs cmþ3) 2.39 ; 10 4.71 ; 103
2

ö 0.427877 0.321414

0 0.02044 0.02819

A

0

E

esc

/E

0

þ

1

0.17 0.18

0.0124 0.0158

5.95 5.45

7.92 7.26

Note.--All quantities are expressed in gravitational units (GM ¼ c ¼ 1), except as indicated.


No. 1, 2007

PARTICLE ACCELERATION IN SHOCKED DISKS

427

issue is beyond the scope of the present paper and will be pursued in future work. However, we briefly discuss a few of the probable scenarios for the production of the observed radio, X-ray, and -ray emission by the jet particles below. Scenario 1.--Anyakoha et al. (1987 ) and Morrison et al. (1984) have suggested that high-energy protons in the jets may undergo proton-proton collisions with the ambient medium, resulting in the production of charged and neutral pions. The neutral pions subsequently decay into -rays, and the charged pions decay into electron-positron pairs, which produce further -rays via inverse Compton and annihilation radiation. These processes can explain the observed radio and -ray emission in blazar sources such as 3C 273. Anyakoha et al. (1987 ) demonstrated that the production of the observed -rays requires protons with Lorentz factors in the range 1:4 P P 86, and the production of the observed radio emission requires protons with Lorentz factors in the range 7 P P 300. On the basis of our model, we find that the terminal proton Lorentz factor, þ1 , for the jets in M87 and Sgr Aö is given by þ1 $ 8and þ1 $ 7, respectively (see Table 2). We therefore conclude that the escaping protons accelerated in our model reach energies high enough to produce the observed radio and -ray emission in M87 and Sgr Aö . In addition to the radiation emitted by the escaping, outflowing protons, the relativistic protons remaining inside the disk will produce additional -rays (e.g., Eilek & Kafatos 1983; Bhattacharyya et al. 2003, 2006 ), although the escape of these photons from the disk may be problematic due to - attenuation ( Becker & Kafatos 1995). Scenario 2.--Usov & Smolsky (1998) suggested that the nonthermal -ray and X-ray emission from AGN jets may be produced by energetic electrons accelerated due to collisions between relativistically moving jet clouds and the ambient medium. The collisions accelerate the electrons up to energies $mp c 2 þ1 ,where þ1 is the terminal Lorentz factor of the jet in our model. Gamma rays are produced when the high-energy electrons scatter external radiation. On the basis of the results presented by Usov & Smolsky (1998), we find that the mean energyofthe scattered photons emerging from the jet in this scenario is " P 400þ2 keV, where 1 þ1 $ 7 Y 8 (see Table 2). Additional X-rays and -rays may also be produced as a result of direct upscattering of background photons by beamed electrons if the Lorentz factor þ1 $ 10, which is comparable to the range predicted by our model (e.g., Dermer & Schlickeiser 1993). Scenario 3.--This last scenario is quite different from the first two. Sikora & Madejski (2000) argue that the jets in optically violently variable (OVV ) quasars contain a mixture of compositions, since the pure e× -eþ jets overproduce soft X-rays and the pure e-p jets underproduce nonthermal X-rays. To resolve these issues, they suggest that initially e-p plasma is injected at the base of the jets, with the protons providing the inertia to account for the kinetic luminosity of the jet. This is consistent with our model of proton /electron escape at the shock. Subsequently, e× -eþ pairs are produced via interactions with the hard X-rays/soft -rays emitted from the hot accretion disk corona. This leads to velocity and density perturbations in the jet and also to the formation of shocks, where the particles are further accelerated. The consequence of the Sikora & Madejski (2000) model is that pair production in the accretion-disk corona can be responsible for the fast ($day) variability observed in OVV quasars at very small angles relative to the jet axis, which is also detected in the MeV-GeV -ray regime. We can expect a similar effect in our theoretical model, since the escape of the relativistic protons/electrons also occurs at the base of the jet, at the shock location. Our model therefore provides a natural connection between the particle acceleration

mechanism operating in the disk and the formation of the energetic outflows considered by Sikora & Madejski (2000). 6. CONCLUSIONS In this paper we have demonstrated that particle acceleration at a standing, isothermal shock in an ADAF 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 outflows based on the well-understood concept of firstorder Fermi acceleration at shock waves. The theory employs an exact mathematical approach in order to solve simultaneously the coupled hydrodynamical and particle transport equations. 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. The transport model we have utilized includes the effects of spatial diffusion, first-order Fermi acceleration, bulk advection, and vertical escape through the upper and lower surfaces of the disk. All of the theoretical parameters can be tied down for a source with a given mass, accretion rate, and jet luminosity. The diffusion coefficient model employed in this work is consistent with the behavior expected as particles approach the event horizon and also at large radii. Close to the event horizon, inward advection at the speed of light dominates over outward diffusion, as expected, while at large r, diffusion dominates over advection. The vertical escape of the energetic particles through the upper and lower surfaces of the disk at the shock location occurs via spatial diffusion in the tangled magnetic field. The approach taken here closely parallels the early studies of cosmic-ray shock acceleration. As in those first investigations (e.g., Blandford & Ostriker 1978), we have employed the test particle approximation, in which the pressure of the accelerated particles is assumed to be negligible compared with that of the thermal background gas. This approximation was shown to be valid for our specific applications to M87 and Sgr Aö in Paper II. We have computed the Green's function describing the particle distribution inside the disk, as well as the energy distribution for the relativistic particles that escape at the shock location to power the jets. The escaping particle distributions for model 2 ( M87) and model 5 (Sgr Aö ) each possess high-energy power-law tails, as expected for a Fermi mechanism. The relatively flat slope of the energy distribution also implies the potential for an interesting connection between the acceleration process explored here and the formation of the universal cosmic-ray background. This is


428

LE & BECKER

Vol. 661

consistent with the results of Szabo & Protheroe (1994), who demonstrated that AGNs may be an important source of cosmic rays in the region of the ``knee.'' The coupled, self-consistent model for the disk structure, the acceleration of the relativistic particles, and the production of the outflows presented in this paper represents the first step toward a fully integrated understanding of the relationship between accretion dynamics and the

formation of the relativistic jets commonly observed around radioloud compact sources. The authors would like to thank the referee, Lev Titarchuk, for providing several insightful comments that led to significant improvements in the manuscript.

APPENDIX ORTHOGONALITY OF THE SPATIAL EIGENFUNCTIONS We can establish the orthogonality of the spatial eigenfunctions Yn (r) by starting with the Sturm-Liouville form (see eq. [37]): ! d dYn S (r ) ÏA1÷ × kn !(r)Yn (r) ¼ 0; dr dr where !(r)and S (r) are given by equations (39) and (38), respectively. Let us suppose that kn and km denote two distinct eigenvalues (kn 6¼ km ) with associated spatial eigenfunctions Yn (r) and Ym (r). Since Yn and Ym each satisfy equation (A1) for their respective values of k, we can write & ' ! d dYm S (r) ÏA2÷ Yn (r) × km !(r)Ym (r) ¼ 0; dr dr & ' ! d dYn S (r ) ÏA3÷ Ym (r) × kn !(r)Yn (r) ¼ 0: dr dr Subtracting equation (A3) from equation (A2) yields ! ! d dYm d dYn S (r ) S (r ) Yn (r) þ Ym (r) ¼ (kn þ km ) ! (r)Yn (r)Ym (r): dr dr dr dr We can integrate equation (A4) by parts from r ¼ rS to r ¼ 1 to obtain 1 Z 1 dYm dYm dYn þ Y n (r )S (r ) dr S (r ) dr rS dr dr rS Z1 dYn 1 dYn dYm × þYm (r)S (r) dr S (r ) dr rS dr dr rS Z1 ¼ ( kn þ km ) !(r)Yn (r)Ym (r) dr:
r
S

ÏA4÷

ÏA5÷

Upon cancellation of like terms, this expression reduces to dYm dYn þ Ym ( r ) S (r) Yn (r) dr dr !
1

Z ¼ ( kn þ km )

1

!(r)Yn (r)Ym (r) dr:
rS

ÏA6÷

r

S

On the basis of the asymptotic behaviors of Yn (r)as r ! rS and r ! 1 given by equations (35), we conclude that the left-hand side of equation (A6) vanishes, leaving Z1 !(r)Ym (r)Yn (r) dr ¼ 0; m 6¼ n: ÏA7÷
r
S

This result establishes that Ym and Yn are orthogonal eigenfunctions relative to the weight function ! (r) defined in equation (39).
REFERENCES Becker, P. A., Subramanian, P., & Kazanas, D. 2001, ApJ, 552, 209 Bhattacharyya, S., Bhatt, N., & Misra, R. 2006, MNRAS, 371, 245 Bhattacharyya, S., Bhatt, N., Misra, R., & Kaul, C. L. 2003, ApJ, 595, 317 Blandford, R. D., & Begelman, M. C. 1999, MNRAS, 303, L1 Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29 Chakrabarti, S. K. 1989, PASJ, 41, 1145 ------. 1990, MNRAS, 243, 610 ------. 1996, ApJ, 464, 664 Chakrabarti, S. K., Acharyya, K., & Molteni, D. 2004, A&A, 421, 1

Abramowicz, M. A., & Chakrabarti, S. K. 1990, ApJ, 350, 281 Allen, S. W., Di Matteo, T., & Fabian, A. C. 2000, MNRAS, 311, 493 Anyakoha, M. W., Okeke, P. N., & Okoye, S. E. 1987, Ap&SS, 132, 65 Becker, P. A. 1992, ApJ, 397, 88 Becker, P. A., & Kafatos, M. 1995, ApJ, 453, 83 Becker, P. A., & Kazanas, D. 2001, ApJ, 546, 429 Becker, P. A., & Le, T. 2003, ApJ, 588, 408 Becker, P. A., Le, T., & Dermer, C. D. 2006, ApJ, 647, 539 Becker, P. A., & Subramanian, P. 2005, ApJ, 622, 520


No. 1, 2007

PARTICLE ACCELERATION IN SHOCKED DISKS

429

Chakrabarti, S. K., & Das, S. 2004, MNRAS, 349, 649 Chakrabarti, S. K., & Molteni, D. 1993, ApJ, 417, 671 ------. 1995, MNRAS, 272, 80 Chen, X., Abramowicz, M. A., & Lasota, J. P. 1997, ApJ, 476, 61 Das, S., Chattopadhyay, I., & Chakrabarti, S. K. 2001, ApJ, 557, 983 Dermer, C. D., & Schlickeiser, R. 1993, ApJ, 416, 458 Di Matteo, T., Quataert, E., Allen, S. W., Narayan, R., & Fabian, A. C. 2000, MNRAS, 311, 507 Eilek, J. A., & Kafatos, M. 1983, ApJ, 271, 804 Fermi, E. 1954, ApJ, 119, 1 Ford, H. C., et al. 1994, ApJ, 435, L27 Frank, J., King, A. R., & Raine, D. J. 2002, Accretion Power in Astrophysics (Cambridge: Cambridge Univ. Press) 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 ( Paper I ) ------. 2005, ApJ, 632, 476 ( Paper II ) Lu, J., Gu, W., & Yuan, F. 1999, ApJ, 523, 340 Lu, J., & Yuan, F. 1997, PASJ, 49, 525 ------. 1998, MNRAS, 295, 66

Lynden-Bell, D. 1969, Nature, 223, 690 Molteni, D., Sponholz, H., & Chakrabarti, S. K. 1996, ApJ, 457, 805 Morrison, P., Roberts, D., & Sadun, A. 1984, ApJ, 280, 483 Narayan, R., Kato, S., & Honma, F. 1997, ApJ, 476, 49 Novikov, I., & Thorne, K. S. 1973, in Black Holes, Les Houches 1972, ed. C. DeWitt & B. S. DeWitt ( New York: Gordon & Breach), 343 Owen, F. N., Eilek, J. A., & Kassim, N. E. 2000, ApJ, 543, 611 ´ Paczynsky, B., & Wiita, P. J. 1980, A&A, 88, 23 Protheroe, R. J., & Kazanas, D. 1983, ApJ, 265, 620 Rees, M. J. 1984, ARA&A, 22, 471 Sambruna, R. M., Gambill, J. K., Maraschi, L., Tavecchio, F., Cerutti, R., Cheung, C. C., Urry, C. M., & Chartas, G. 2004, ApJ, 608, 698 Schlickeiser, R. 1989a, ApJ, 336, 243 ------. 1989b, ApJ, 336, 264 Sikora, M., & Madejski, G. 2000, ApJ, 534, 109 Spruit, H. C. 1987, A&A, 184, 173 Subramanian, P., Becker, P. A., & Kazanas, D. 1999, ApJ, 523, 203 Szabo, A. P., & Protheroe, R. J. 1994, Astropart. Phys., 2, 375 Titarchuk, L., & Zannias, T. 1998, ApJ, 493, 863 Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 Usov, V. V., & Smolsky, M. V. 1998, Ap&SS, 259, 331 Webb, G. M., & Bogdan, T. J. 1987, ApJ, 320, 683 Yang, R., & Kafatos, M. 1995, A&A, 295, 238