Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.stsci.edu/~tle/papers/fokker_planck.pdf Äàòà èçìåíåíèÿ: Thu Sep 18 07:06:49 2008 Äàòà èíäåêñèðîâàíèÿ: Tue Oct 2 12:33:26 2012 Êîäèðîâêà: koi8-r Ïîèñêîâûå ñëîâà: ï ï ï ï ï ï ï ï ï ï

The Astrophysical Journal, 647:539 ­ 551, 2006 August 10
# 2006. The American Astronomical Society. All rights reserved. Printed in U.S.A.

TIME-DEPENDENT STOCHASTIC PARTICLE ACCELERATION IN ASTROPHYSICAL PLASMAS: EXACT SOLUTIONS INCLUDING MOMENTUM-DEPENDENT ESCAPE
Peter A. Becker and Truong Le and Charles D. Dermer
E. O. Hulburt Center for Space Research, Naval Research Laboratory, Washington, DC 20375; truong.le@nrl.navy.mil, dermer@gamma.nrl.navy.mil Received 2005 November 28; accepted 2006 April 22
1

Center for Earth Observing and Space Research, George Mason University, Fairfax, VA 22030-4444; pbecker@gmu.edu

ABSTRACT Stochastic acceleration of charged particles due to interactions with magnetohydrodynamic ( MHD) plasma waves is the dominant process leading to the formation of the high-energy electron and ion distributions in a variety of astrophysical systems. Collisions with the waves influence both the energization and the spatial transport of the particles, and therefore it is important to treat these two aspects of the problem in a self-consistent manner. We solve the representative Fokker-Planck equation to obtain a new, closed-form solution for the time-dependent Green's ´ function describing the acceleration and escape of relativistic particles interacting with Alfven or fast-mode waves characterized by momentum diffusion coefficient D( p) / pq and mean particle escape timescale tesc ( p) / pqþ2 , where p is the particle momentum and q is the power-law index of the MHD wave spectrum. In particular, we obtain solutions for the momentum distribution of the ions in the plasma and also for the momentum distribution of the escaping particles, which may form an energetic outflow. The general features of the solutions are illustrated via examples based on either a Kolmogorov or Kraichnan wave spectrum. The new expressions complement the results obtained by Park and Petrosian, who presented exact solutions for the hard-sphere scattering case (q ¼ 2) in addition to other scenarios in which the escape timescale has a power-law dependence on the momentum. Our results have direct relevance for models of high-energy radiation and cosmic-ray production in astrophysical environments such as -ray bursts, active galaxies, and magnetized coronae around black holes. In particular, we outline an application of the new results to black hole sources that produce outflows of relativistic hadrons, with associated predictions that can be tested using GLAST. Subject headings: acceleration of particles -- black hole physics -- cosmic rays -- galaxies: jets -- g methods: analytical -- plasmas

1. INTRODUCTION Observations of high-energy radiation from a variety of astrophysical sources imply the presence of significant populations of nonthermal (often relativistic) particles. Much of the observational data is characterized by strong variability on very short timescales. Nonthermal distributions are naturally produced via the Fermi process, in which particles interact with scattering centers moving systematically and /or randomly. The first-order Fermi mechanism treats particle acceleration in converging flows, such as shocks, and is thought to be important at the Earth's bow shock, in certain classes of solar flares, for cosmic-ray acceleration by supernova remnants, and in sources with relativistic outflows, including blazars and -ray bursts ( for reviews, see Blandford & Eichler 1987; Kirk 1994). The second-order process, as it was originally conceived by Fermi, involved the stochastic acceleration of particles scattering with randomly moving magnetized clouds ( Fermi 1949). Later refinements of this idea replaced the magnetized clouds with magnetohydrodynamic ( MHD) waves (e.g., Melrose 1980). The second-order, stochastic Fermi process now finds application in a wide range of astrophysical settings, including solar flares ( Miller et al. 1990; Liu et al. 2004a), clusters of galaxies (Schlickeiser et al. 1987; Petrosian 2001; Brunetti
Also Department of Physics and Astronomy, George Mason University, Fairfax, VA 22030-4444.
1

et al. 2004), the Galactic center ( Liu et al. 2004b; Atoyan & Dermer 2004), and -ray bursts ( Waxman 1995; Dermer & Humi 2001). The standard approach to modeling the acceleration of nonthermal particles via interactions with MHD waves involves obtaining the solution to a steady-state transport equation that incorporates treatments of systematic and stochastic Fermi acceleration, radiative losses, and particle escape (e.g., Schlickeiser 1984; Schlickeiser & Steinacker 1989; Liu et al. 2006). However, the prevalence of variability on short timescales in many sources calls into question the validity of the steady-state interpretation. Park & Petrosian (1995) provided a comprehensive review of the various time-dependent solutions that have been derived in the past 30 years. In most cases, it is assumed that the momentum diffusion coefficient, D( p), and the mean particle escape timescale, tesc ( p), have power-law dependences on the particle momentum p. Although the set of existing solutions covers a broad range of possible values for the associated power-law indices, there are certain physically interesting cases for which no analytical solution is currently available. For example, Dermer et al. (1996) have shown that for the stochastic acceleration of relativistic particles due to resonant interactions with plasma waves in black hole magnetospheres, one obtains D( p) / pq and tesc ( p) / pqþ2 , where q is the index of the wavenumber distribution (see x 2). The analytical solution for the time-dependent Green's function in this situation is of particular physical interest, but it has not appeared in the previous literature. This has 539


540

BECKER, LE, & DERMER

Vol. 647

motivated us to reexamine the associated Fokker-Planck transport equation for this case, and to obtain a new family of closedform solutions for the secular Green's function. The resulting expression, describing the evolution of an initially monoenergetic particle distribution, complements the set of solutions discussed by Park & Petrosian (1995). Our primary goal in this paper is to present a detailed derivation of the exact analytical solution and to demonstrate some of its key properties. Detailed applications of our results to the modeling of particle acceleration in black hole accretion coronae and other astrophysical environments will be presented in a separate paper. The remainder of the paper is organized as follows. In x 2we review the fundamental equations governing the acceleration of charged particles interacting with plasma waves. The transport equation is solved in x 3 to obtain the time-dependent Green's function, and illustrative results are presented in x 4. The astrophysical implications of our work are summarized in x 5, and additional mathematical details are provided in the Appendix. 2. FUNDAMENTAL EQUATIONS Charged particles in turbulent astrophysical plasmas are expected to be accelerated via interactions with whistler, fast-mode, ´ and Alfven waves propagating in the magnetized gas. Here we consider a simplified isotropic description of the wave energy distribution, denoted by W (k ), where W (k )dk represents the energy density due to waves with wavenumber between k and k × dk . The transport formalism we consider assumes a powerlaw distribution for the wave-turbulence spectrum , which implies definite relations between the momentum diffusion coefficient and the momentum-dependent escape timescale (e.g., Melrose 1974; Schlickeiser 1989a, 1989b; Dermer et al. 1996). MHD waves injected over a narrow range of wavenumber cascade to larger wavenumbers, forming a Kolmogorov or Kraichnan power spectrum over the inertial range with W (k ) / k þq , where q ¼ 5/3 and q ¼ 3/2 for the Kolmogorov and Kraichnan cases, respectively (e.g., Zhou & Matthaeus 1990). The specific forms we adopt for the transport coefficients in x 2.2 are based on the physics of the resonant scattering processes governing the wave-particle interactions (see Dermer et al. 1996). We assume that the nonthermal particle distribution is isotropic, and focus on a detailed treatment of the propagation of particles in momentum space due to wave-particle interactions. The spatial aspects of the transport (i.e., the confinement of the particles in the acceleration region) are treated in an approximate manner using a momentum-dependent escape term. 2.1. Transport Equation The fundamental transport equation describing the propagation of particles in the momentum space can be written in the flux-conservation form (e.g., Becker 1992; Schlickeiser 1989a, 1989b) & !' @f 1@ @f f S ( p; t ) 2 ¼þ 2 p A( p) f þ D( p) × ; þ @t p @p @p tesc ( p) 4 p 2 Ï 1÷ where p is the particle momentum, f ( p; t ) is the particle distribution function, D( p) denotes the momentum diffusion coefficient, tesc ( p) is the mean escape time, S ( p; t ) represents particle sources, and A( p) describes any additional, systematic acceleration or loss processes, such as shock acceleration or synchrotron /inverseCompton emission. The quantity in square brackets in equation (1) describes the flux of particles through the momentum

space ( Tademaru et al. 1971), and the source term is defined so that S ( p; t ) dp dt gives the number of particles injected into the plasma per unit volume between times t and t × dt with momenta between p and p × dp. The total particle number density n(t ) and energy density U (t ) of the distribution f ( p; t )are computed using Z1 Z1 4p 2 f ( p; t )dp; U (t ) ¼ 4p 2 f ( p; t )dp; Ï2÷ n(t ) ¼
0 0

where the particle kinetic energy is related to the Lorentz factor and the particle momentum p by ¼ ( þ 1)mc 2 ; ¼ p2 ×1 m2 c 2 1
=2

;

Ï 3÷

and m and c denote the particle rest mass and the speed of light, respectively. Rather than solve equation (1) directly to determine f ( p; t ) for a given source term S ( p; t ), it is more instructive to first solve for the Green's function, fG ( p0 ; p; t0 ; t ), which satisfies the equation & !' @ fG 1@ @ fG 2 ¼þ 2 p A( p) fG þ D( p) p @p @t @p fG ( p þ p0 )(t þ t0 ) × þ ; Ï 4÷ 2 tesc ( p) 4 p0 where p0 is the initial momentum and t0 is the initial time. The source term in this equation corresponds to the injection of a single particle per unit volume with momentum p0 at time t0 . The particle number and energy densities associated with the Green's function are given by Z1 4p 2 fG ( p0 ; p; t0 ; t )dp; nG ( t ) ¼ 0 Z1 U G (t ) ¼ 4p 2 fG ( p0 ; p; t0 ; t )dp: Ï 5÷
0

Once the Green's function solution has been determined, the particular solution associated with an arbitrary source distribution S ( p; t ) can be computed using the integral convolution (e.g., Becker 2003) Z tZ 1 f ( p; t ) ¼ fG ( p0 ; p; t0 ; t )S ( p0 ; t0 )dp0 dt0 ; Ï 6÷
0 0

where we have assumed that the particle injection begins at time t ¼ 0 and no particles are present in the plasma for t < 0. 2.2. Transport Coefficients In the Appendix (x A1), we demonstrate that for arbitrary particle energies, the mean rate of change of the particle momentum due to stochastic acceleration is related to the momentum diffusion coefficient D( p) via D dp E 1 dþ 2 à pD: ¼2 Ï 7÷ dt stoch p dp The corresponding result for the mean rate of change of the kinetic energy is (see x A2 of the Appendix and Miller & Ramaty 1989) D d E 1 dþ 2 à p vD ; ¼2 Ï 8÷ dt stoch p dp


No. 1, 2006

STOCHASTIC PARTICLE ACCELERATION and we have introduced the dimensionless constants a Aö ; Dö 1 : D ö tö

541

where v is the particle speed. If the MHD wave spectrum has the ´ power-law form W / k þq associated with Alfven and fast-mode waves, then the momentum dependences of the diffusion coefficient D( p) and the mean escape time tesc ( p) describing the resonant pitch-angle scattering of relativistic particles are given by (e.g., Dermer et al. 1996; Miller & Ramaty 1989) q qþ2 p p 22 ; tesc ( p) ¼ tö ; Ï 9÷ D ( p) ¼ D ö m c mc mc where Dö / sþ1 and tö / s are constants. We shall focus on transport scenarios with q 2, so that tesc is either a decreasing or constant function of the momentum p. In order to maintain the physical validity of the escape-timescale approach used here, we must require that tesc exceed the light-crossing time L/c for a source with size L. This implies a fundamental upper limit to the particle momentum when q < 2. By combining equations (7) and (9), we find that the mean rate of change of the momentum for relativistic particles accelerated stochastically by MHD waves is given by qþ1 D dp E p ¼ (q × 2)Dö mc : Ï10÷ dt stoch mc For simplicity, we assume that the momentum dependence of the additional, systematic loss/acceleration processes appearing in the transport equation (1), described by the coefficient A( p), mimics that of the stochastic acceleration (eq. [10]). We therefore write qþ1 p ; Ï11÷ A( p) ¼ Aö mc mc where the constant Aö / sþ1 determines the positive (negative) rate of systematic acceleration ( losses). Note that this formulation precludes the treatment of loss processes with a quadratic energy dependence (e.g., inverse Compton or synchrotron) since that would imply q ¼ 3, which is outside the range considered here. However, first-order Fermi acceleration at a shock or energy losses due to Coulomb collisions can be treated by setting q ¼ 2 with Aö either positive or negative, respectively. This suggests that the results obtained here are relevant primarily for the transport of energetic ions. However, even in this application one needs to bear in mind that synchrotron and inverse Compton losses will become dominant at sufficiently high energies (e.g., Schlickeiser 1984; Schlickeiser & Steinacker 1989; Liu et al. 2006). It is convenient to transform to the dimensionless momentum and time variables x and y, defined by x p ; mc y Dö t ; Ï12÷

Ï15÷

Note that x equals the particle Lorentz factor in applications involving ultrarelativistic particles. The constant a describes the relative importance of systematic gains or losses compared with the stochastic process. 2.3. Fokker-Planck Equation Additional physical insight can be obtained by reorganizing equation (13) in the Fokker-Planck form ö @ NG @2 @á ¼ 2 Ïx q NG ÷ þ ( q × 2 × a) x q þ 1 N G @x @y @x þ x 2þq NG × (x þ x0 )( y þ y0 );

Ï16÷

where we have defined the Green's function number distribution NG using NG (x0 ; x; y0 ; y) 4m3 c 3 x 2 fG (x0 ; x; y0 ; y): Ï17÷

The Fokker-Planck coefficients appearing in equation (16), which describe the evolution of the particle distribution due to the influence of stochastic and systematic processes, are given by ( Reif 1965) 1 d 2 ¼ x q; 2 dy D dx E ¼ ( q × 2 × a) x dy
q þ1

;

Ï18÷

in which case the transport equation (4) for the Green's function becomes @ fG 1@ a @ þ 1×q à 2 ×q @ f G x x fG ¼2 þ2 x @x x @x @y @x (x þ x0 )( y þ y0 ) þ x 2 þq f G × ; Ï13÷ 2 4 m 3 c 3 x 0 where p0 ; x0 mc y0 Dö t0 ; Ï14÷

where the first coefficient describes the ``broadening'' of the distribution due to momentum space diffusion, and the second represents the mean ``drift'' of the particles (i.e., the average acceleration rate). Equation (16) is equivalent to equation (49) from Park & Petrosian (1995) if we set their parameters app ¼ 2 × a and spp ¼ q þ 2, where spp denotes the power-law index describing the momentum dependence of the escape timescale. Our particular choice for spp reflects the physics of the resonant wave-particle interactions, as represented by equations (9). The Fokker-Planck form of equation (16) clearly reveals the fundamental nature of the transport process. In particular, we note that in the limit a ! 0, the drift coefficient hdx/dyi reduces to the purely stochastic result (eq. [10]), as expected when systematic gains / losses are excluded from the problem. The total number and energy densities are computed in terms of NG using (cf. eq. [5]) Z1 nG ( y ) ¼ NG (x0 ; x; y0 ; y)dx; Z 10 UG ( y) ¼ xNG (x0 ; x; y0 ; y)dx; Ï19÷
0

where (see eq. [3]) " ¼ mc
2

þ

x ×1

2

à

#
1=2

þ1 :

Ï20÷

Since the physical situation considered here corresponds to the injection of a single particle per unit volume at ``time'' y ¼ y0 ,it follows that nG ( y0 ) ¼ 1.


542

BECKER, LE, & DERMER 3. SOLUTION FOR THE GREEN'S FUNCTION

Vol. 647

The result obtained for the Green's function number distribution NG has two important special cases depending on whether q ¼ 2or q < 2. Park & Petrosian (1995) obtained the exact solution to equation (16) for the hard-sphere case ( Ramaty 1979) with q ¼ 2 and their parameter spp ¼ 0, corresponding to a momentumindependent escape timescale. In this section we derive the exact solution to the time-dependent problem with q < 2 and spp ¼ q þ 2, which describes the physics of the wave-particle interactions (see eqs. [9]). 3.1. Laplace Transformation We define the Laplace transformation of NG using Z L(x0 ; x; s)
0 1

where M and U denote the confluent hypergeometric functions (Abramowitz & Stegun 1970) and pffiffiffi s × (a × 3) pffiffiffi ; 2(2 þ q)  a×3 : 2þq Ï28÷

The constants A and B appearing in equation (27) are determined by ensuring that the function L is continuous at z ¼ z0 , and that it also satisfies the derivative jump condition dL lim !0 dz
z0 × z0 þ

c0 eþsy0 ¼þ (2 þ q)z

2 0

(3 z0 c0

þ2q)=(2þq)

¼

þe 2x

þsy0

0

pffiffiffi ; Ï29÷

e

þsy

NG (x0 ; x; y0 ; y)dy: R

Ï21÷
1 þsy 0e

By operating on the Fokker-Planck equation (16) with we obtain d2 q dá (q × 2 × a)x Ï x L÷ þ dx dx 2
q þ1

dy,

ö L þ x
þsy

obtained by integrating the transport equation (25) with respect to z in a small range around the source momentum z0 . The constant B can be eliminated by combining the continuity and derivative jump conditions. After some algebra, the solution obtained for A is e
þsy
0

2þq

L þ sL Ï22÷ A¼þ e
z0 =2 (a×2)=(qþ2) z0

¼ þe or, equivalently, d 2L × dx 2

0

(x þ x0 );

M (; ; z0 ) pffiffiffi ; 2x 0 W ( z 0 )

Ï30÷

! q þ 2 þ a dL (1 þ q)(2 × a) s × þ 2qþ2 þ q L x dx x2 x x eþsy0 (x þ x0 ) : ¼þ xq Ï23÷

where W (z) denotes the Wronskian, defined by W (z) M (; ; z) d d U (; ; z) þ U (; ; z) M (; ; z): dz dz Ï31÷ Using equation (13.1.22) from Abramowitz & Stegun (1970), we find that W is given by the exact expression W (z) ¼ þ þ(  ) z þ  e z ; þ( ) Ï32÷

This equation can be transformed into standard form by introducing the new momentum variables pffiffiffi 2 x z(x) 2þq
2 þq

;

pffiffiffi 2 x z 0 ( x0 ) 2þq

2þq 0

:

Ï24÷

After some algebra, we find that the equation for L now becomes z
2

which can be combined with equation (30) and the continuity relation to obtain A¼ þ( )z0
 ×(a×2)=(qþ2) þsy

! d 2 L a × 1 dL (1 þ q)(2 × a) z 2 sz z × × þþ L dz 2 q þ 2 dz 4 (2 þ q)2 c0 (2 þ q)2 c0 eþsy0 (z þ z0 ) z (3þ2q)=(2þq) ; ¼þ 2þq c0 Ï25÷

e 0 eþ z pffiffiffi 2x0 þ( )

0

=2

M (; ; z0 )

;

Ï33÷

B¼

þ( )z

 ×(a×2)=(qþ2) þsy e 0

0 þz e pffiffiffi 2x0 þ(  )

0

=2

U (; ; z0 )

:

Ï34÷

where pffiffiffi 2 : c0 2þq The solutions to equation (25) obtained for z 6¼ high- and low-energy boundary conditions are & þz=2 (a×2)=(2þq) AU (; ; z); z L(z0 ; z; s) ¼ e BM (; ; z); Ï26÷ z0 that satisfy the given by z ! z0 ; z z0 ;

The final solution for the Laplace transformation L can therefore be written as L(z0 ; z; s) ¼
 ( þ( ) z 0 z pffiffiffi þ(  )2x0 z0 a×2)=(2þq)

;e where

þsy

0

e

þ(z×z0 )=2

M (; ; z

min

)U (; ; z

max

);

Ï35÷

Ï27÷

z

min

min(z; z0 );

z

max

max (z; z0 ):

Ï36÷


No. 1, 2006

STOCHASTIC PARTICLE ACCELERATION

543

Since the poles are associated with the function þ( ) in equation (35), we need to make use of the identity lim (s þ sn )þ( ) ¼ pffiffiffi (þ1)n 2(2 þ q) ; n! Ï42÷

s!sn

which follows from equations (28) and (38). Combining equations (35), (41), and (42), we find that the residues are given by Res(sn ) ¼ (þ1)n (2 þ q)e sn ( n!þ( )x0 ;e
þ(z×z0 )=2 yþy0 )  z0

( z z0

a×2)=(2þq)

M (þ n;; z

min

)U (þ n;; z

max

):

Ï43÷

Based on equations (13.6.9) and (13.6.27) from Abramowitz & Stegun (1970), we note that the confluent hypergeometric functions appearing in this expression reduce to Laguerre polynomials, and therefore our result for the residue can be rewritten after some simplification as
Fig. 1.--Integration contour C used to evaluate eq. (37 ).

Res(sn ) ¼

n! e ;e

sn ( yþy0 )  z0

(2 þ q) z ( þ(  × n) x 0 z0
( (z)Pn þ1)

a×2)=(2þq)

3.2. Inverse Transformation The solution for the Green's function NG can be found using the complex Mellin inversion integral, which states that Z ×i1 1 e sy L(z0 ; z; s)ds; Ï37÷ NG (z0 ; z; y0 ; y) ¼ 2i þi1 where is chosen so that the line Re s ¼ lies to the right of any singularities in the integrand. The singularities are simple poles located where jþ( )j ! 1, which corresponds to ¼ þn, with n ¼ 0, 1, 2, : : : Equation (28) therefore implies that there are an infinite number of poles situated along the real axis at s ¼ sn ,where pffiffiffi sn ¼ ½2(q þ 2)n þ a þ 3 ; n ¼ 0; 1; 2; : : : Ï38÷ We can therefore employ the residue theorem to write I 1 X e sy L(z0 ; z; s)ds ¼ 2i Res(sn );
n ¼0

þ(z×z0 )=2 (  þ1) Pn

(z0 );

Ï44÷

where P

(  þ1) n

(z) denotes the Laguerre polynomial.

3.4. Closed-Form Expression for the Green's Function The result for the Green's function number distribution NG is obtained by summing the residues (see eq. [40]), which yields N G ( z 0 ; z ; y0 ; y) ¼
1 X n! e s n ( yþy0 )  z0

;e

n ¼0 þ(z×z0 )=2 (  þ1) Pn

(2 þ q) z þ(  × n)x0 z0
( (z)Pn þ1)

(a×2)=(2þq)

(z0 ):

Ï45÷

Ï39÷

This convergent sum is a useful expression for the Green's function. However, further progress can be made by employing the bilinear generating function for the Laguerre polynomials, given by equation (8.976) from Gradshteyn & Ryzhik (1980). After some algebra, we find that the general closed-form solution can be written in the form pffiffiffiffiffiffiffiffi 2 þ q x (a×1)=2 zz0 NG (x0 ; x; y0 ; y) ¼ x0 x0 1þ ! pffiffiffiffiffiffiffiffi (z × z0 )(1 × ) 2 zz0 ; exp þ Iþ1 ; 2(1 þ ) 1þ Ï46÷

where C is the closed integration contour indicated in Figure 1 and Res(sn ) denotes the residue associated with the pole at s ¼ sn . Based on asymptotic analysis, we conclude that the contribution to the integral due to the curved portion of the contour vanishes in the limit R ! 1, and consequently we can combine equations (37 ) and(39)toobtain N G ( z 0 ; z ; y0 ; y) ¼
1 X n¼0

Res(sn ):

Ï40÷ where pffiffiffi p 2 2 þq 2 x ; z0 (x0 ) z(x) 2þq 2þ a×3 ; ( y; y0 ) e 2(qþ2)(  2þq ffiffiffi x q
2þq 0

Hence we need only evaluate the residues in order to determine the solution for the Green's function. 3.3. Evaluation of the Residues The residues associated with the simple poles at s ¼ sn are evaluated using the formula (e.g., Butkov 1968) Res(sn ) ¼ lim (s þ sn )e sy L(z0 ; z; s):
s!s
n

; ; Ï47÷

pffiffi yþy0 )

Ï41÷

which is the main result of the paper. Note that the solution for NG depends on the time parameters y and y0 only through the ``age'' of the injected particles, y þ y0 , as indicated by the form


544

BECKER, LE, & DERMER

Vol. 647

for . In the limit ! 0, corresponding to infinite escape time, the Green's function number distribution reduces to a=2 (xx0 )(3þq)=2 x ¼ N G ( x0 ; x; y0 ; y) 2( y þ y ) x ¼0 (2 þ q) x0 0 0 " # " # 2 þq 2 þq (x × x0 ) 2(xx0 )(2þq)=2 ; exp þ I  þ1 : Ï48÷ (2 þ q)2 ( y þ y0 ) (2 þ q)2 ( y þ y0 ) The exact solution for the time-dependent Green's function describing the evolution of a monoenergetic initial spectrum with q < 2 given by equation (46) represents an interesting new contribution to particle transport theory. The corresponding solution for the hard-sphere case with q ¼ 2, given by equation (43) of Park & Petrosian (1995), can be stated in our notation as N G ( x0 ; x; y0 ; y) ¼ ( eþk( yþy0 ) x pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 x0 ( y þ y0 ) x0
a×1)=2

tion for the number distribution associated with continual monoenergetic particle injection can be written as N ( x ; y ) 4 m 3 c 3 x 2 f ( x ; y ) ¼ Z N0 y NG (x0 ; x; y0 ; y)dy0 ; Dö 0 Ï53÷

where we have used equation (14) to make the substitution dt0 ¼ dy0 /Dö . Since NG depends on the time parameters y and y0 only through the combination y þ y0 (see eqs. [46] and [49]), it follows that Z N0 y N ( x; y) ¼ NG (x0 ; x; 0; y 0 )dy 0 : Ï54÷ Dö 0 This is a more convenient form for N (x; y) because y now appears only in the upper integration bound. For general, finite values of y, the time-dependent particular solution for N (x; y) given by equation (54) must be computed numerically by substituting for NG using either equation (46) or (49), depending on the value of q. However, as y ! 1, the solution rapidly approaches a stationary result representing a balance between injection, acceleration, and particle escape. The form of the stationary solution, called the ``steady-state Green's G function,'' Nss , can be obtained by directly solving the transport equation (1) with @ f /@ t ¼ 0 and S ( p; t ) ¼ N0 ( p þ p0 ). Alternatively, the steady-state Green's function can also be computed by taking the limit of the time-dependent solution, which yields Z N0 1 G NG (x0 ; x; 0; y 0 )dy 0 : Ï55÷ Nss (x) lim N (x; y) ¼ y!1 Dö 0 In the q < 2 case, we can substitute for NG using equation (46) and evaluate the resulting integral using equation (6.669.4) from Gradshteyn & Ryzhik (1980). After some algebra, we obtain (a×1)=2 N0 x G (xx0 )(2þq)=2 Nss (x) ¼ (2 þ q)x0 Dö x0 pffiffiffi 2þq ! pffiffiffi 2þq ! xmin xmax ; IÏþ1÷=2 KÏþ1÷=2 ; Ï56÷ 2þq 2þq where  ¼ (a × 3)/(2 þ q) and x
min

q ¼2

! þ(ln x þ ln x0 )2 ; exp ; 4( y þ y0 )

Ï49÷

where k (a × 1)2 × 2 × a × : 4 Ï50÷

We note that the general solution for NG given by equation (46) agrees with equation (49) in the limit q ! 2 as required. Furthermore, equation (46) allows q to take on negative values if desired, and it is also applicable over a broad range of both positive and negative values for a. Recall that when a ¼ 0, there are no systematic acceleration or loss processes included in the model. Positive (negative) values for a imply additional systematic acceleration ( losses). 3.5. Transition to the Stationary Solution The analytical solutions we have obtained for the Green's function provide a complete description of the response of the system to the impulsive injection of monoenergetic particles at any desired time. The generality of these expressions allows one to obtain the particular solution for the distribution function f associated with any arbitrary momentum-time source function S using the convolution given by equation (6). One case of special interest is the spectrum resulting from the continual injection of monoenergetic particles beginning at time t ¼ 0, described by the source term ( N 0 ( p þ p0 ) ; t ! 0; Ï51÷ S ( p; t ) ¼ 0; t < 0; where N0 denotes the rate of injection of particles with momentum p0 per unit volume. We assume that no particles are present for t < 0. Combining equations (6) and (51), we find that the time-dependent distribution function resulting from monoenergetic particle injection is given by f ( p; t ) ¼ N0 Z
0 t

min(x; x0 );

x

max

max (x; x0 ):

Ï57÷

Likewise, for the case with q ¼ 2, we can substitute for NG using equation (49) and then utilize equation (3.471.9) from Gradshteyn & Ryzhik (1980) to conclude that the steady-state solution is given by (a×1)=2 pffiffi N0 x x max þ k G pffiffiffi ¼ ; Ï58÷ Nss (x) q ¼2 x min 2x0 Dö k x0 where k is defined by equation (50). The steady-state solutions given by equations (56) and (58) agree with the results obtained by directly solving the transport equation. Due to the asymptotic behavior of the Bessel K (z) function for large z, equation (56) G indicates that Nss exhibits an exponential cutoff at high energies when q < 2 (Abramowitz & Stegun 1970). This corresponds physically to the fact that the escape timescale decreases with increasing particle momentum in this case. However, when q ¼ 2 (eq. [58]), the spectrum displays a pure power-law behavior at

fG ( p0 ; p; t0 ; t )dt0 :

Ï52÷

By transforming to the dimensionless variables x and y and employing equation (17), we conclude that the particular solu-


No. 1, 2006

STOCHASTIC PARTICLE ACCELERATION

545

high energies because the escape timescale is independent of the particle momentum. Specific examples illustrating these behaviors will be presented in x 4. 3.6. Escaping Particle Distribution The various expressions we have obtained for the Green's function NG describe the momentum distribution of the particles remaining in the plasma at time t after injection occurring at time t0 . However, since our model incorporates a physically realistic, momentum-dependent escape timescale tesc ( p) given by equation (9), it is also quite interesting to compute the spectrum of the escaping particles, which may form an energetic outflow capable of producing observable radiation. In general, the number dis tribution of the escaping particles, N esc (x; y), is related to the current distribution of particles in the plasma, N (x; y), via N
esc

nario can be analyzed by using equation (54) to substitute for N (x; y) in equation (59), which yields Zy NG (x0 ; x; 0; y 0 ) dy 0 : Ï64÷ N esc (x; y) ¼ N0 x 2þq
0

In the limit y ! 1, the escaping particle spectrum approaches the steady-state result Z1 esc NG (x0 ; x; 0; y 0 ) dy 0 : Nss (x) lim N esc (x; y) ¼ N0 x 2þq
y!1 0

Ï65÷ By comparing equations (55), (62), and (65), we conclude that esc Nss (x) ¼ Dö x
2 þq G Nss (x) ¼ N0 à N esc G

( x) :

Ï66÷

( x; y)

þ1 tesc

N (x; y) ¼ Dö x

2 þq

N ( x; y) ;

Ï59÷

where we have used equations (9) and (15) to obtain the final result. The quantity N esc (x; y)dx represents the number of particles escaping per unit volume per unit time with dimensionless momenta between x and x × dx. An important special case is the evolution of the escaping particle spectrum resulting from impulsive monoenergetic injection at dimensionless time y ¼ y0 . In this application, equation (59) gives the Green's function number spectrum for the escaping particles, defined by esc NG (x0 ; x; y0 ; y) Dö x
2þq

The final result can be combined with equation (63) to show that the escaping spectrum satisfies the normalization relation Z1 esc Nss (x)dx ¼ N0 ; Ï67÷
0

as expected for the case of continual steady-state injection. Note that analytical expressions for the steady-state escaping G esc spectrum Nss (x) can be obtained by substituting for Nss (x) in equation (66) using either equation (56) or (58), depending on the value of q. We obtain N 0 x 2 þq esc Nss (x) ¼ (2 þ q) x0 ;I for q < 2 and esc Nss for q ¼ 2. 4. NUMERICAL RESULTS The new result we have derived for the secular Green's function (eq. [46]) displays a rich behavior through its complex dependence on momentum, time, and the dimensionless parameters q, a,and , which are related to the fundamental physical transport coefficients Dö , Aö , and tö via equations (15). Here we present several example calculations in order to illustrate the utility of the new solution. Detailed applications to astrophysical situations, including active galaxies and -ray bursts, will be presented in subsequent papers. The panels on the left-hand side of Figure 2 depict the timedependent Green's function solution, NG , describing the evolution of the particle distribution in the plasma resulting from impulsive monoenergetic injection at y ¼ 0. Results are presented for the hard-sphere scattering case (q ¼ 2), computed using equation (49), and for the Kolmogorov (q ¼ 5/3) and Kraichnan (q ¼ 3/2) cases, evaluated using equation (46). The only acceleration mechanism considered here is the stochastic acceleration associated with the second-order Fermi process, and therefore we set a ¼ 0. The escape parameter is set equal to unity, so that the N 0 x 2 þq pffiffiffi (x) ¼ 2x 0 k ( x x0
a×1)=2 Ï þ1÷=2

NG (x0 ; x; y0 ; y);

Ï60÷

where NG is evaluated using either equation (46) or (49) de esc pending on the value of q. We can use this expression for NG to compute the total escaping number distribution resulting from an impulsive flare occurring at t ¼ 0 by integrating the escaping distribution with respect to time, which yields à
esc NG

(a×1)=2 x (xx0 )(2þq)=2 x0 pffiffiffi 2þq ! pffiffiffi 2þq ! xmin xmax KÏþ1÷=2 2þq 2þq

Ï68÷

Z ( x)
0

1

N

esc G

dt ¼ x

2 þq

Z
0

1

NG (x0 ; x; 0; y )dy ; Ï61÷

0

0



x x

max min

þpffiffi k Ï69÷

where we have used equation (60) and taken advantage of the fact that NG depends on y and y0 only through the difference y þ y0 (cf. eq. [54]). By comparing equations (55) and (61), we deduce that
esc à NG

Dö ( x) ¼ x N0

2þq

G Nss (x);

Ï62÷

so that the total escaping spectrum is simply proportional to the steady-state Green's function resulting from continual injection, as expected. The process considered here corresponds to the injection of a single particle at time t ¼ 0, and therefore we find esc that the normalization of à NG is given by Z
0 1 esc à NG (x)dx ¼ 1;

Ï63÷

which provides a useful check on the numerical results. Another interesting example is the case of continual monoenergetic particle injection commencing at time t ¼ 0, described by the source term given by equation (51). The time-dependent buildup of the escaping particle spectrum N esc (x; y) in this sce-


546

BECKER, LE, & DERMER

Vol. 647

Fig. 2.-- Green's function solutions to the time-dependent stochastic particle acceleration equation for x0 ¼ 1. The left panels depict the impulsive-injection solution, NG (eqs. [46] and [49]), and the right panels illustrate the response to uniform, continuous injection, N (eq. [54] ), with N0 ¼ Dö ¼ 1. Note that N approaches the corresponding steady-state solution (eqs. [56] and [58]) as y increases. The indices of the wave turbulence spectra are indicated, with q ¼ 2 for hard-sphere scattering, q ¼ 5/3 for a Kolmogorov cascade, and q ¼ 3/2 for a Kraichnan cascade.

timescale for escape is comparable to the diffusion timescale. As the wave turbulence spectrum becomes steeper (i.e., as q increases), a larger fraction of the turbulence energy is contained in waves with long wavelengths, which interact resonantly with higher energy particles. Steeper turbulence spectra therefore produce harder particle distributions, as can be confirmed in the plots. Consequently we conclude that an ensemble of hard-sphere scattering centers is more effective at accelerating nonthermal particles compared with a Kolmogorov wave spectrum, which in turn is more effective than a Kraichnan spectrum. The panels on the right-hand side of Figure 2 illustrate the buildup of the particle spectrum in the plasma, N (x; y), due to

continual monoenergetic injection beginning at y ¼ 0, computed by evaluating numerically the integral in equation (54). We have set a ¼ 0 and therefore the acceleration is purely stochastic. As y ! 1, the spectrum approaches the steady-state form given by equation (56) for q < 2 or by equation (58) for q ¼ 2. In the hard-sphere case (q ¼ 2), the particle spectrum displays a power-law shape at high energies in agreement with equation (58). However, when q < 2, particle escape dominates over acceleration at high energies, and therefore the steady-state distribution is truncated, even in the absence of radiative losses. This effect produces the quasi-exponential turnovers exhibited by the stationary spectra when q ¼ 5/3 and q ¼ 3/2. Particle escape


No. 1, 2006

STOCHASTIC PARTICLE ACCELERATION

547

Fig. 3.--Evolution of the particle distribution in the plasma resulting from monoenergetic injection with q ¼ 5/3, a ¼ 0, ¼ 0:1, and x0 ¼ 1. Panel a depicts the Green's function, NG , resulting from impulsive monoenergetic injection (eq. [46]), and panel b illustrates the response to continual monoenergetic injection (eq. [54]) and the corresponding steady-state solution (eq. [56]) for N0 ¼ Dö ¼ 1.

in these cases mimics energy losses due to, for example, synchrotron emission from electrons. Figures 3 and 4 illustrate the effects of varying the values of the escape parameter and the systematic acceleration/loss parameter a for the q ¼ 5/3 case. In Figure 3, we set a ¼ 0 and ¼ 0:1, which represents a particle escape timescale that is an order of magnitude larger than the corresponding case depicted in Figure 2. The longer escape time allows the particles to be accelerated to higher mean energies before diffusing out of the plasma, and the decay of the particle density at late times takes place much more slowly than for larger values of . The hardening of the spectra causes the quasi-exponential cutoffs to move to higher energies, and the same effect can also be noted in the corresponding stationary solutions. The calculations represented in Figure 4 include additional systematic acceleration processes that are modeled by setting a ¼ 2. The enhanced particle acceleration further hardens the spectrum and shifts the cutoff to even higher energies. Although the slope of the low-energy particle distribution (x < x0 ) is not altered much when only is varied (see Figs. 2 and 3), this slope becomes significantly steeper when a is increased, as indicated in Figure 4.

In the left-hand panels of Figures 5 and 6 we plot the timedependent solution for the Green's function describing the es esc caping particles, NG , resulting from impulsive monoenergetic particle injection (eq. [60]) when q ¼ 5/3. The right-hand panels illustrate the buildup of the escaping spectrum, N esc (eq. [64]), resulting from continual injection beginning at y ¼ 0. Note the esc transition to the steady-state solution, Nss (eq. [65]), as y ! 1. In Figure 5 we set ¼ 0:1 and a ¼ 0, and in Figure 6 we set ¼ 0:1 and a ¼ 2. Included for comparison are the corresponding spectra describing the particle distributions in the plasma at the same values of y. We point out that the escaping particle spectra are significantly harder than the in situ distributions, which reflects the preferential escape of the high-energy particles resulting from the momentum dependence of the escape timescale when q < 2 (see eq. [9]). 5. DISCUSSION AND SUMMARY We have derived new, closed-form solutions (eqs. [46] and [68]) for the time-dependent Green's function representing the stochastic acceleration of relativistic ions interacting with MHD waves. The analytical results we have obtained describe the

Fig. 4.--Same as Fig. 3, except for a ¼ 2.


548

BECKER, LE, & DERMER

Vol. 647

Fig. 5.--Evolution of the particle distribution resulting from monoenergetic injection with q ¼ 5/3, a ¼ 0, ¼ 0:1, x0 ¼ 1, N0 ¼ 1, and Dö ¼ 1. Panel a treats the case of impulsive injection, with the thin lines representing the particle distribution in the plasma (eq. [46]) and the thick lines denoting the escaping particle spectrum (eq. [60]) in units of . Panel b illustrates the response to continual monoenergetic injection for the particle distribution in the plasma (eqs. [54] and [56]; thin lines) and for the escaping particle spectrum (eqs. [64] and [66]; thick lines). See the discussion in the text.

time-dependent distributions for both the accelerated (in situ) and the escaping particles. The Fokker-Planck transport equation considered here includes momentum diffusion with coefficient D( p) / pq , particle escape with mean timescale tesc / pqþ2 , and additional systematic acceleration/losses with a rate proportional to pqþ1 , where p is the particle momentum and q is the index of the wave turbulence spectrum. This specific scenario describes ´ the resonant interaction of particles with fast-mode and Alfven waves, which is one of the fundamental acceleration processes in high-energy astrophysics (e.g., Dermer et al. 1996). The new analytical result complements the work of Park & Petrosian (1995) since it is applicable for any q < 2 provided spp ¼ q þ 2, where spp is the power-law index used by these authors to describe the momentum dependence of the escape timescale. The most closely related previous result is given by their equation (59), which treats the case with q 6¼ 2and spp ¼ 0, corresponding to a momentum-independent escape timescale. Our analytical solution (eq. [46]) agrees with theirs in the limit q ! 2, as expected. The general features of our new solution were dis-

cussed in x 4, where it was demonstrated that increasingly hard particle spectra result from larger values of the wave index q,smaller values of the escape parameter , and larger values of the systematic acceleration parameter a. The rich behavior of the solution as a function of momentum, time, and the parameters q, , and a provides useful physical insight into the nature of the coupled energetic/spatial particle transport in astrophysical plasmas. The solutions presented here can be used to describe the acceleration and transport of relativistic ions in astrophysical environments in which the turbulence spectrum is very poorly known and can be approximated by a power law, such as -ray bursts, active galaxies, magnetized coronae around black holes, and the intergalactic medium in clusters of galaxies. For example, the hard X-ray emission from black hole jet sources such as Cygnus X-1 ( Malzac et al. 2006) and the microquasar LS 5039 (Aharonian et al. 2005) could be powered by the stochastic acceleration of ions in a black hole accretion disk corona that subsequently escape and interact with surrounding material. In this scenario, persistent acceleration of monoenergetic particles

Fig. 6.--Same as Fig. 5, except for a ¼ 2.


No. 1, 2006

STOCHASTIC PARTICLE ACCELERATION

549

injected into the corona, or flaring episodes averaged over a sufficiently long time, would produce a time-averaged escaping distribution of relativistic protons given by equation (68) for q < 2. Assuming only stochastic acceleration, so that a ¼ 0, the number distribution of escaping particles with x ! x0 takes the form pffiffiffi 2þq ! x esc (7þ3q)=2 Nss (x) / x KÏþ1÷=2 2þq 8 > 3þ2q 2 þ q 1=(2þq) >x > pffiffiffi ; xT ; > < pffiffiffi 2þq ! / > 2 þ q 1=(2þq) x > (5þ2q)=2 >x > ; x 3 pffiffiffi exp þ : : 2þq Ï70÷ When the escaping hadrons collide with ambient gas or stellar wind material, they would generate X-rays and -rays via a pion production cascade with a very hard spectrum leading up to a quasiexponential cutoff. The Gamma-ray Large Area Space Telescope (GLAST ) will be able to provide detailed spectra from Galactic black hole sources and unidentified EGRET -ray sources to test for the existence of this hard component. Our new analytical solution can also be used to model the variability of radiation from ions accelerated in the accretion-disk coronae of Seyfert

galaxies by changing the level of turbulence. Flaring $100 MeV to GeV radiation could be weakly detected by GLAST as a consequence of this process. Additional applications of our work include studies of the stochastic acceleration of relativistic cosmic rays in -ray bursts ( Dermer & Humi 2001). The results we have obtained describe both the momentum distribution of the particles in the plasma, and the momentum distribution of the particles that escape to form energetic outflows. Since the solutions do not include inverse Compton or synchrotron losses, which are usually important for energetic electrons, they are primarily applicable to cases involving ion acceleration. However, the new solution can be used to treat relativistic electrons if the escape timescale is shorter than the synchrotron /inverseCompton timescale. In order to treat the acceleration of relativistic electrons, a generalized calculation including both stochastic particle acceleration and radiative losses is desirable, and we are currently working to extend the analytical model presented here to incorporate these effects. Beyond the direct utility of the new analytical solutions for probing the nature of particle acceleration in astrophysical sources, we note that they are also useful for benchmarking numerical simulations.

T. L. is funded through NASA GLAST Science Investigation DPR-S-1563-Y. The work of C. D. D. and P. A. B. is supported by the Office of Naval Research. The authors also acknowledge the useful comments provided by the anonymous referee.

APPENDIX In this appendix we explore the relationship between the momentum diffusion coefficient D( p) and the mean rate of change of the particle momentum and kinetic energy, using the integral method outlined by Subramanian et al. (1999). A1. MEAN RATE OF CHANGE OF MOMENTUM DUE TO STOCHASTIC ACCELERATION In the case of pure stochastic acceleration, the transport equation (1) reduces to ! @f 1@ @f ¼2 p 2 D( p) : @t p @p @p The mean momentum, hpi, is defined as a function of t by 1 h pi n Z
0 1

ÏA1÷

4p3 f ( p; t )dp;

ÏA2÷

where the number density n is related to f via equation (2). The associated mean rate of change of the momentum is given by D dp E 1 Z 1 @f dp: ÏA3÷ 4 p3 ¼ dt n0 @t Combining this relation with the transport equation (A1) yields D dp E 1 Z 1 @ @f p2D 4 p ¼ dp: dt n0 @p @p Integrating by parts once gives Z D dp E 1 ¼þ dt n0 and integrating by parts a second time yields D dp E 1 Z ¼ dt n0
1 1

ÏA4÷

4 p 2 D

@f dp; @p

ÏA5÷

4 f

dþ 2 à p D dp: dp

ÏA6÷


550

BECKER, LE, & DERMER

Vol. 647

It is sufficient for our purposes to consider the evolution of the particle distribution f satisfying the monoenergetic initial condition ¼ A0 ( p þ p0 ): ÏA7÷ f ( p; t )
t ¼t0

Combining this relation with equation (A6), we find that at the initial time D dp E 1 dþ 2 pD ¼2 dt p dp
t ¼t0

t ¼ t0 , à :
p ¼p
0

ÏA8÷

Dropping the subscript ``0'' without loss of generality, we conclude that the mean rate of change of the momentum for particles with momentum p at time t is given by D dp E 1 dþ 2 à pD; ÏA9÷ ¼2 dt p dp which agrees with equation (7). A2. MEAN RATE OF CHANGE OF KINETIC ENERGY DUE TO STOCHASTIC ACCELERATION The mean kinetic energy, hi, is defined by hi 1 n Z
0 1

4 p 2 f d p;

ÏA10÷

where is given by equation (3). The associated mean rate of change is given by D d E 1 Z 1 @f dp; 4 p 2 ¼ dt n0 @t which can be combined with the transport equation (A1) to obtain D d E 1 Z 1 @ @f p 2D 4 ¼ dp: dt n0 @p @p As in the previous section, we integrate by parts to find that Z D d E 1 ¼þ dt n Next we employ the derivative relation d ¼v dp and integrate by parts again to obtain D d E 1 Z ¼ dt n0
1

ÏA11÷

ÏA12÷

1

4
0

d 2 @ f pD dp: dp @p

ÏA13÷

ÏA14÷

4 f

dþ 2 à vp D dp: dp

ÏA15÷

Applying the initial condition given by equation (A7), we find that at the initial time t ¼ t0 , D d E 1 d þ 2 à vp D ¼2 : dt p dp
t ¼t0 p¼p
0

ÏA16÷

Dropping the subscript ``0'' without loss of generality, we conclude that the mean rate of change of the kinetic energy for particles with momentum p at time t is given by D d E 1 dþ 2 à p vD ; ¼2 dt p dp which agrees with our equation (8) and also with equation (13) from Miller & Ramaty (1989). ÏA17÷


No. 1, 2006

STOCHASTIC PARTICLE ACCELERATION

551

REFERENCES Abramowitz, M., & Stegun, I. A. 1970, Handbook of Mathematical Functions Melrose, D. B. 1974, Sol. Phys., 37, 353 ( New York: Dover) ------. 1980, Plasma Astrophysics ( New York: Gordon and Breach) Aharonian, F., et al. 2005, Science, 309, 746 Miller, J. A., Guessoum, N., & Ramaty, R. 1990, ApJ, 361, 701 Atoyan, A., & Dermer, C. D. 2004, ApJ, 617, L123 Miller, J. A., & Ramaty, R. 1989, ApJ, 344, 973 Becker, P. A. 1992, ApJ, 397, 88 Park, B. T., & Petrosian, V. 1995, ApJ, 446, 699 ------. 2003, MNRAS, 343, 215 Petrosian, V. 2001, ApJ, 557, 560 Blandford, R., & Eichler, D. 1987, Phys. Rep., 154, 1 Ramaty, R. 1979, in AIP Conf. Proc. 56, Particle Acceleration Mechanisms in Brunetti, G., Blasi, P., Cassano, R., & Gabici, S. 2004, MNRAS, 350, 1174 Astrophysics, ed. J. Arons, C. Max, & C. McKee ( New York: AIP), 135 Butkov, E. 1968, Mathematical Physics ( London: Addison-Wesley) Reif, F. 1965, Fundamentals of Statistical and Thermal Physics ( New York: Dermer, C. D., & Humi, M. 2001, ApJ, 556, 479 McGraw-Hill ) Dermer, C. D., Miller, J. A., & Li, H. 1996, ApJ, 456, 106 Schlickeiser, R. 1984, A&A, 136, 227 Fermi, E. 1949, Phys. Rev., 75, 1169 ------. 1989a, ApJ, 336, 243 Gradshteyn, I. S., & Ryzhik, I. M. 1980, Table of Integrals, Series, and Products ------. 1989b, ApJ, 336, 264 ( New York: Academic Press) Schlickeiser, R., Sievers, A., & Thiemann, H. 1987, A&A, 182, 21 Kirk, J. G. 1994, Plasma Astrophysics, ed. A. O. Benz & T. J.-L. Courvoisier Schlickeiser, R., & Steinacker, J. 1989, Sol. Phys., 122, 29 ( Berlin: Springer), 225 Subramanian, P., Becker, P. A., & Kazanas, D. 1999, ApJ, 523, 203 Liu, S., Melia, F., & Petrosian, V. 2006, ApJ, 636, 798 Tademaru, E., Newman, C. E., & Jones, F. C. 1971, Ap&SS, 10, 453 Liu, S., Petrosian, V., & Mason, G. M. 2004a, ApJ, 613, L81 Waxman, E. 1995, Phys. Rev. Lett., 75, 386 Liu, S., Petrosian, V., & Melia, F. 2004b, ApJ, 611, L101 Zhou, Y., & Matthaeus, W. H. 1990, J. Geophys. Res., 95, 14881 Malzac, J., et al. 2006, A&A, 448, 1125