Cosmic ray acceleration at relativistic shock waves

# Diffusive cosmic ray acceleration at relativistic shock waves with magnetostatic turbulence

R. Schlickeiser Institut für Theoretische Physik, Lehrstuhl IV: Weltraum- und Astrophysik, Ruhr-Universität Bochum D-44780 Bochum, Germany
###### Abstract

The analytical theory of diffusive cosmic ray acceleration at parallel stationary shock waves with magnetostatic turbulence is generalized to arbitrary shock speeds , including in particular relativistic speeds. This is achieved by applying the diffusion approximation to the relevant Fokker-Planck particle transport equation formulated in the mixed comoving coordinate system. In this coordinate system the particle’s momentum coordinates and are taken in the rest frame of the streaming plasma, whereas the time and space coordinates are taken in the observer’s system. For magnetostatic slab turbulence the diffusion-convection transport equation for the isotropic (in the rest frame of the streaming plasma) part of the particle’s phase space density is derived. For a step-wise shock velocity profile the steady-state diffusion-convection transport equation is solved. For a symmetric pitch-angle scattering Fokker-Planck coefficient the steady-state solution is independent of the microphysical scattering details. For nonrelativistic mono-momentum particle injection at the shock the differential number density of accelerated particles is a Lorentzian-type distribution function which at large momenta approaches a power law distribution function with the spectral index . For nonrelativistic () shock speeds this spectral index agrees with the known result , whereas for ultrarelativistic () shock speeds the spectral index value is close to unity.

acceleration of particles – cosmic rays – relativity – shock waves

## 1. Introduction

One of the most important problems of modern astrophysics is to explain how cosmic ray particles are accelerated to relativistic energies in powerful sources of nonthermal radiation. Diffusive first-order Fermi acceleration at nonrelativistic shock fronts has been regarded as a prime candidate for particle acceleration in astrophysics (for reviews see Drury 1983; Blandford and Eichler 1987). Modern TeV air-Cherenkov telescopes have indeed resolved the shock regions in supernova remnants and identied the shocks as strong emission regions of TeV photons generated by the accelerated particles (Hinton and Hofmann 2009).

It is likely that diffusive shock acceleration also operates efficiently at magnetized shock waves with relativistic speeds. Such relativistic shocks form during the interaction of relativistic supersonic and super-Alfvenic outflows with the ambient ionized interstellar or intergalactic medium (Gerbig and Schlickeiser 2009) producing anisotropic counterstream plasma distribution functions due to shock-reflected charged particles in the upstream medium (Spitkovsky 2008, Sironi and Spitkovsky 2009). Relativistic outflows are a direct consequence of violent explosive events such as in gamma-ray burst sources (Piran 1999) both in the collapsar (Woosley 1993, Paczynski 1998) and supranova (Vietri and Stella 1998) models, but also occur as highly collimated pulsar winds and jets of active galactic nuclei with initial bulk Lorentz factors .

The transport and acceleration of energetic particles in the partially turbulent cosmic magnetic fields associated with shocks is described using the Fokker-Planck equation for the particle distribution function (for a recent derivation, see Schlickeiser 2011). The diffusion approximation for the particle density in the rest frame of the fluid is a well-known simplied form of the Fokker-Planck equation, which results when turbulent pitch-angle scattering is strong enough to ensure that the scale of the particle density variation is signicantly greater than the particle mean free path (Jokipii 1966, Hasselmann and Wibberenz 1968, Earl 1974). Numerical studies confirmed the accuracy of the diffusion approximation in a uniform mean guide magnetic field (Kota et al. 1982).

While for nonrelativistic shock waves the analytic theory of diffusive shock acceleration is well developed (Axford et al. 1977, Krymsky 1987, Blandford and Ostriker 1978, Bell 1978, Drury 1983), for relativistic shock speeds such an analytical theory does not exist sofar even for parallel shock waves, although the underlying Fokker-Planck transport equation (see Eq. (3) below) for the particle dynamics has already been derived (Webb 1985; Kirk, Schneider and Schlickeiser 1988) some years ago. The existing literature concentrated on semi-numerical eigenfunction solutions of this Fokker-Planck transport equation pioneered by Kirk and Schneider (1987), relativistic Monte Carlo simulations (Ellison et al. 1990, Ostrowski 1991, Bednarz and Ostrowski 1998, Summerlin and Baring 2012), and relativistic particle in-cell simulations (Spitkovsky 2008, Sironi and Sptkovsky 2009). It is the purpose of this manuscript to develop an analytical study of cosmic ray acceleration in parallel relativistic magnetized shock waves employing the diffusion approximation in the upstream and downstream regions of the shock wave. The development runs much in parallel with the existing work on nonrelativistic shocks.

## 2. Basic equations

Magnetized space plasmas such as the interstellar medium harbour low-frequency linear () transverse MHD waves (such as shear Alfven and magnetosonic plasma waves) with dispersion relations and , respectively, in the rest frame of the moving plasma. Faradays induction law then indicates for MHD waves that the strength of turbulent electric fields is much smaller than the strength of turbulent magnetic fields. The ordering corresponds to the derivation of cosmic ray transport equations for from the collisionfree Boltzmann equation fir the full phase space distribution to the Fokker-Planck equation for its gyrotropic part , and to the diffusion-convection transport equation for its isotropic part , respectively (for a recent review see Schlickeiser 2011). Accordingly, the cosmic ray anisotropy, defined as the deviation

 g(→X,v,μ,t)=f0(→X,p,μ,t)−F(→X,p,t), (1)

then is small () with respect to . The diffusion approximation applied to the Fokker-Planck transport equation for allows us to relate the cosmic ray anisotropy to the solutions of the diffusion-convection tranport equation for .

Because of the gyrorotation of the cosmic ray particles in the uniform magnetic field, one is not so much interested in their actual position as in the coordinates of the cosmic ray guiding center

 →X=(X,Y,Z)=→x+cqaB20→p×→B0=→x+cqaB0⎛⎜⎝py−px0⎞⎟⎠, (2)

where we orient the large-scale guide magnetic field, which is uniform on the scales of the cosmic ray particles gyradii , along the -axis. and denote the speed and the relativistic gyrofrequency of a cosmic ray particle with mass , charge and energy .

The Larmor-phase averaged Fokker-Planck transport equation in a medium, propagating with the stationary bulk speed with aligned along the magnetic field direction is given by (Webb 1985; Kirk, Schneider and Schlickeiser 1988)

 Γ[1+Uvμc2]∂f0∂t+Γ[U+vμ]∂f0∂z+v(1−μ2)2L∂f0∂μ
 −α(z)(μ+Uv)[μp∂f0∂p+(1−μ2)∂f0∂μ]+Rf0−S(→X,p,t)
 =∂∂μ[Dμμ∂f0∂μ+Dμσ∂f0∂yσ]
 +p−2∂∂yωp2[Dωμ∂f0∂μ+Dωσ∂f0∂yσ], (3)

irrespective of how the Fokker-Planck coefficients are calculated, either by quasilinear (Schlickeiser 2002) or nonlinear (Shalchi 2009) cosmic ray transport theories. denotes the rate of adiabatic deceleration/acceleration in relativistic flows

 α(z)=c2U(z)dΓ(z)dz=dUdzΓ3=d(UΓ)dz (4)

In the Fokker-Planck Eq. (3) the phase space coordinates have to be taken in the mixed comoving coordinate system (time and space coordinates in the laboratory (=observer) system and particle’s momentum coordinates and in the rest frame of the streaming plasma). Moreover, in Eq. (3) we use the Einstein sum convention for indices, and represent the three phase space variables with non-vanishing stochastic fields and . Consequently, the term on the right-hand side generally represents 9 different Fokker-Planck coefficients: but, depending on the turbulent fields considered, not all of them are non-zero and some are much larger than others. accounts for additional sources and sinks of particles.

The focusing length (Roelof 1969)

represents the spatial gradient of the guide magnetic field . It is the only term, resulting from the mirror force in the large scale inhomogeneous guide magnetic field, that we consider (neglecting possible drift effects). for a diverging guide magnetic field; for a converging guide magnetic field.

 (6)

represents continuous () and catastrophic () momentum losses of cosmic ray particles.

For spatially constant flows the rate of adiabatic deceleration/acceleration (4) vanishes, and the remaining flow velocity () dependent terms in Eq. (3) simply result from the Lorentz transformation of special relativity of the laboratory-frame position-time coordinates to the mixed-comoving-frame position-time coordinates . However, for spatially varying flow speeds special relativity no longer applies and has to be replaced by the transformation laws from general realtivity. As noted by Riffert (1986) as well as Kirk, Schneider and Schlickeiser (1988) these introduce connection coefficients or Christoffel symbols of the first kind. In a flat Euclidean space-time metric the terms proportional to in Eq. (3) are exactly these connection coefficients.

## 3. Magnetostatic slab turbulence

As important special case we consider magnetostatic (vanishing turbulent electric fields), isospectral, slab () turbulence. Then the only nonvanishing Fokker-Planck coefficient is . In this case the Fokker-Planck transport equation (3) simplifies to

 Γ[1+Uvμc2]∂f0∂t+Γ[U+vμ]∂f0∂z+v(1−μ2)2L∂f0∂μ
 −α(z)(μ+Uv)[μp∂f0∂p+(1−μ2)∂f0∂μ]+Rf0
 =S(→X,p,t)+∂∂μ[Dμμ∂f0∂μ] (7)

Due to the rapid pitch angle scattering the gyrotropic particle distribution function adjusts very quickly to a distribution function which is close to the isotropic distribution in the rest frame of the moving background plasma. Defining the isotropic part of the phase space density as the -averaged phase space density

 F(→X,p,t)≡12∫1−1dμf0(→X,p,μ,t), (8)

we follow the analysis of Jokipii (1966) and Hasselmann and Wibberenz (1968) to split the total density into the isotropic part and an anisotropic part ,

 f0(→X,p,μ,t)=F(→X,p,t)+g(→X,p,μ,t) (9)

where because of Eq. (8)

 ∫1−1dμg(→X,p,μ,t)=0 (10)

Inserting the ansatz (9) yields for the Fokker-Planck equation (7)

 Γ[1+Uvμc2](∂F∂t+∂g∂t)+Γ[U+vμ](∂F∂z+∂g∂z)
 −α(z)(μ+Uv)[μp∂F∂p+μp∂g∂p+(1−μ2)∂g∂μ]
 +v(1−μ2)2L∂g∂μ+RF+Rg−S(→X,p,t)=∂∂μ[Dμμ∂g∂μ] (11)

Averaging this equation over , using that , leads to the diffusion-convection transport equation

 Γ[∂F∂t+U∂F∂z]−α3p∂F∂p+RF−S(→X,p,t)
 +v2[Γ(∂∂z+Uc2∂∂t)+1L]∫1−1dμμg
 −αU2v[p∂∂p+2]∫1−1dμμg−α2[p∂∂p+3]∫1−1dμμ2g=0, (12)

involving the first and second moment of the anisotropy.

Subtracting Eq. (12) from the Fokker-Planck equation (11) provides

 Γvμ(∂F∂z+Uc2∂F∂t)+Γvμ(∂g∂z+Uc2∂g∂t)+Γ(∂g∂t+U∂g∂z)
 +α3p∂F∂p+Rg−α(μ+Uv)[μp∂F∂p+μp∂g∂p+(1−μ2)∂g∂μ]
 −[vΓ2(∂∂z+Uc2∂∂t)+v2L]∫1−1dμμg
 +αU2v[p∂∂p+2]∫1−1dμμg+α2[p∂∂p+3]∫1−1dμμ2g
 +v(1−μ2)2L∂g∂μ=∂∂μ[Dμμ∂g∂μ] (13)

We note that Eqs. (12) and (13) are still exact.

### 3.1. Diffusion approximation for cosmic shock waves

With we approximate the anisotropy equation (13) to leading order by

 Γvμ(∂F∂z+Uc2∂F∂t)+α[13−μ2−Uμv]p∂F∂p
 ≃∂∂μ[Dμμ∂g∂μ] (14)

In previous nonrelativistic () studies (Jokipii 1966, Hasselmann and Wibberenz 1968, Schlickeiser et al. 2007, Schlickeiser and Shalchi 2008) only the first term on the left-hand side of Eq. (14) has been considered, so that

 vμ∂F∂z≃∂∂μ[Dμμ∂g∂μ] (15)

For gradual flows with non-zero values of for an extended region of space, we have to consider also the term resulting from the momentum gradient . This case will be considered elsewhere.

Here we consider flows with non-zero values of in a very limited region of space such as cosmic shock waves. We chose as laboratory frame the rest frame of the shock wave with the step-wise velocity profile

 U(z)={−U1=const.for 0

with . In this case the rate of adiabatic acceleration (4)

 α=α0δ(z),α0=−(U1Γ1−U2Γ2) (17)

is non-zero only at the position of the shock. We therefore neglect in the anisotropy equation (14) the term resulting from the momentum gradient , and additionally neglect the time derivate of as compared to the spatial gradient of , i.e.

 Uc2∂F∂t≪∂F∂z, (18)

so that Eq. (14) reduces to the streaming anisotropy contribution only

 Γvμ∂F∂z≃∂∂μ[Dμμ∂g∂μ], (19)

which differs from the nonrelativistic anisotropy equation (15) by the additional factor .

Integrating Eq. (19) provides

 Dμμ∂g∂μ=c0+Γvμ22∂F∂z (20)

The integration constant (with respect to ) is determined from the property that the left-hand side of this equation vanishes for , yielding , providing

 ∂g∂μ=−Γv(1−μ2)2Dμμ(μ)∂F∂z (21)

A further integration over together with Eq. (10) provides for the cosmic ray anisotropy

 g(→X,p,μ,t)≃Γv4[∫1−1dμ(1−μ)(1−μ2)Dμμ(μ)
 −2∫μ−1dx(1−x2)Dμμ(x)]∂F∂z, (22)

which is determined by the spatial gradient of the isotropic distribution function with respect to .

### 3.2. Anisotropy moments

With the anisotropy (22) we readily determine the moments needed in the diffusion-convection transport equation (12) as

 ∫1−1dμμg=−ΓvK04∂F∂z (23)

and

 ∫1−1dμμ2g=−ΓvK16∂F∂z (24)

in terms of the two () integrals

 Kn=∫1−1dμμn(1−μ2)2Dμμ(μ) (25)

For the often considered case of symmetric pitch-angle Fokker-Planck coefficients the integral vanishes.

### 3.3. Diffusion-convection transport equation

For further reduction of the diffusion-convection transport equation (12) it is convenient to introduce as new momentum variable

 s=lnp (26)

and to define the convection and diffusion operators

 CF=−α3p∂F∂p+v2L∫1−1dμμg=−α3∂F∂s+v2L∫1−1dμμg (27)

and

 DF=vΓ2∂∂z∫1−1dμμg−αU2v[p∂∂p+2]∫1−1dμμg
 −α2[p∂∂p+3]∫1−1dμμ2g=vΓ2∂∂z∫1−1dμμg
 −αU2v[∂∂s+2]∫1−1dμμg−α2[∂∂s+3]∫1−1dμμ2g, (28)

respectively. With the assumption (18) we then can write the diffusion-convection transport equation (12) in the compact form

 Γ[∂F∂t+U∂F∂z]+CF+DF+RF=S(→X,p,t) (29)

Using for any momentum dependent quantity the identity

 1v∂∂s[A(s)M(p)]=∂∂s[A(s)vM(p)]+A(s)vγ2M(p)
 =∂∂s[A(s)vM(p)]+A(s)(1−β2)vM(p) (30)

with the cosmic ray particle Lorentz factor and the dimensionless particle velocity , the diffusion operator (28) reads

 DF=vΓ2∂∂z∫1−1dμμg−α2∂∂s[Uv∫1−1dμμg+∫1−1dμμ2g]
 −α2[Uv(3−β2)∫1−1dμμg+3∫1−1dμμ2g] (31)

With the moments (21) - (22) we derive after straightforward but tedious algebra

 DF=−Γ∂∂z[Γκzz∂F∂z]+Γ∂∂s[κsz∂F∂z]
 +Γαv4(K1+3−β22UvK0)∂F∂z, (32)

where we introduce the two diffusion coefficients

 κzz=v2K08,
 κsz=κzs=αv12(K1+3U2vK0) (33)

Likewise, the convection operator (27) becomes with the first moment (23)

 CF=−ΓκzzL∂F∂z−α3∂F∂s (34)

Inserting the operators (32) and (34) provides for the diffusion-convection transport equation (29)

 Γ∂F∂t+Γ[U−κzzL+αv4(K1+3−β22UvK0)]∂F∂z+RF
 −α3∂F∂s−Γ∂∂z[Γκzz∂F∂z]+Γ∂∂s[κzs∂F∂z]=S(→X,p,t) (35)

Applying the identity

 ∂∂s[κM]=1p2∂∂p[p3κM]−3κM, (36)

for any momentum dependent quantity to the last term on the left-hand side of Eq. (35) we obtain

 ∂∂s[κsz∂F∂z]=1p2∂∂pp3[κsz∂F∂z]−3κsz∂F∂z
 =1p2∂∂pp2[κszp∂F∂z]−3κsz∂F∂z (37)

Consequently, we find as alternative form for the diffusion-convection transport equation (35)

 Γ∂F∂t+RF−Γ∂∂z[Γκzz∂F∂z]+Γp2∂∂pp2κszp∂F∂z
 +Γ[U−κzzL−3κsz+αv4(K1+3−β22UvK0)]∂F∂z
 −α3∂F∂s=S(→X,p,t) (38)

With the diffusion coefficients (33) we obtain for the difference

 αv4(K1+3−β22UvK0)−3κsz=−αUK0β28
 =−αUc2κzz=−κzzdΓdz, (39)

where we used Eq. (4). The diffusion-convection transport equation (38) then becomes

 Γ∂F∂t−κzzΓL∂F∂z+Γ[U∂F∂z−αp3Γ∂F∂p]−dΓdzΓκzz∂F∂z
 −Γ∂∂z[Γκzz∂F∂z]+Γp2∂∂pp2κpz∂F∂z+RF=S(→X,p,t) (40)

with

 κpz=pκsz=αKpz,
 Kpz=vp12(K1+3U2vK0) (41)

We note the identity

 U∂F∂z−αp3Γ∂F∂p=1Γ∂∂z[ΓUF]−FΓddz(ΓU)
 −1Γp2∂∂p(αp3F3)+αFΓ=1Γ∂∂z[ΓUF]−1Γp2∂∂p(αp3F3), (42)

because

 ddz(ΓU)=ΓdUdz+UdΓdz=ΓdUdz[1+U2c2Γ2]=α (43)

Likewise,

 dΓdzκzzΓ∂F∂z+Γ∂∂z[Γκzz∂F∂z]
 =(dΓdz+Γ∂∂z)[Γκzz∂F∂z]=∂∂z[Γ2κzz∂F∂z] (44)

With these two identities the diffusion-convection transport equation (40) finally reads

 Γ∂F∂t−κzzΓL∂F∂z+∂∂z[Γ(UF−Γκzz∂F∂z)]
 +1p2∂∂p[p2κpzΓ∂F∂z−αp3F3]+RF=S(→X,p,t) (45)

Eq. (45) is the first important new result of this study: the diffusion-convection transport equation of cosmic rays in aligned parallel flows of arbitrary speed containing magnetostatic slab turbulence with the cosmic ray phase space coordinates taken in the mixed comoving coordinate system. It is particularly appropriate to investigate cosmic ray particle acceleration in parallel relativistic flows.

In the limit of nonrelativistic flows so that . the transport equation (45) becomes

 ∂F∂t−κzzL∂F∂z+∂∂z[UF−κzz∂F∂z]
 +1p2∂∂p[p2κpz∂F∂z−αp3F3]+RF=S(→X,p,t), (46)

which differs from the transport theory used in earlier nonrelativistic diffusive shock acceleration theory (Axford et al. 1977, Krymsky 1987, Blandford and Ostriker 1978, Bell 1978, Drury 1983) by the additional third last term on the left-hand side involving , which results from our correct handling of the connection coefficients in Eq. (3). As we will demonstrate below, this additional term provides a major modification of the resulting differential momentum spectrum of accelerated particles in the nonrelativistic flow limit at nonrelativistic particles momenta: instead of a power law distribution of accelerated particles at the shock a Lorentzian distribution function results, which at large momenta then approaches the power law distribution inferred in earlier acceleration theories for nonrelativistic shock speeds.

### 3.4. Slab Alven waves

There are four different magnetostatic slab Alven waves: forward and backward propagating waves which each can be left- or right-handed circularly polarized, respectively. In terms of the cross helicity and the magnetic helicities of forward and backward moving Alfven waves, and power-law type wave intensities with with the same spectral index of all four waves (isospectral turbulence), the two integrals (25) are given by (Dung and Schlickeiser 1990, Schlickeiser 2002)

 Kn=64s−1(RLkmin)2−svkmin(B0δB)2Sn(Hc,σ+,σ−) (47)

for with

 S0=2(2−s)(4−s)G(Hc,σ+,σ−)
 S1=Z[σ++σ−+Hc(σ+−σ−)](3−s)(5−s)G(Hc,σ+,σ−), (48)

where

 G(Hc,σ+,σ−)=4−[σ++σ−+Hc(σ+−σ−)]2 (49)

is always positive. is the cosmic ray gyroradius, the sign of the cosmic ray particle and the smallest wavenumber of the Alfven waves with total magnetic field strength . It is obvious that the two integrals (47) exhibit the same momentum dependence . Moreover, is positive for and negative for .

Defining the maximum cosmic ray momentum with

 pm=|Q|eB0kminc=1.5⋅1016|Q|B(μG)λ100eVc, (50)

where pc denotes the maximum wavelength of the Alfven waves, we obtain for the integrals (47) for

 Kn=KSnv(ppm)2−s (51)

with the constant length

 K=32π(s−1)(B0δB)2λmax=3.1⋅1021λ100(B0δB)2 (52)

For larger momenta the integrals (47) are zero.

Consequently, we find for the diffusion coefficients (33) and (41)

 κzz=vK8S0(ppm)2−s,
 κpz=αKpz,Kpz=pK12(S1+3US02v)(ppm)2−s (53)

## 4. Particle acceleration at relativistic shock waves

We adopt the step-like shock profile (16) and assume particle injection at the position of the shock only with the injection momentum spectrum . Moreover, we assume spatially constant flow velocities and diffusion coefficients in the upstream and downstream region.

In the steady-state case with no losses () and a uniform background magnetic field () the diffusion-convection transport equation (45) in the rest frame of the shock wave reduces to

 ∂∂z[Γ(UF−Γκzz∂F∂z)]−α0δ(z)p2∂∂p[p2KpzΓ∂F∂z−p3F3]
 =S(p)δ(z=0) (54)

We will solve Eq. (54) by the same method as for nonrelativistic parallel step-like shock waves (Axford et al. 1977, Krymsky 1987, Blandford and Ostriker 1978, Bell 1978, Drury 1983).

For the upstream region the steady-state transport equation (54) yields

 F1(z,p)=F0(p)exp[−U1zΓ1κzz,1], (55)

which approaches zero far upstream .

For the downstream region the solution of the steady-state transport equation (54) is given by

 F2(z,p)=F0(p), (56)

which is finite far downstream . At the position of the shock

 F1(z=0,p)=F2(z=0,p)=F0(p) (57)

the distribution function is continuous.

### 4.1. Momentum spectrum of accelerated particles at the shock

The particle momentum spectrum at the position of the shock is obtained by integrating the transport equation (54) from to and considering the limit . This provides the continuity condition for the cosmic ray streaming density at the shock

 Γ1(U1F1−Γ1κzz,1∂F1∂z)0+−Γ2(U2F2−Γ2κzz,2∂F2∂z)0−
 (58)

With the up- and downstream solutions (55) and (56) we obtain

 Γ2U2F0(p)+U1Γ1−U2Γ23p2ddp(p3+3p2Kpz,1U1κzz,1)F0(p)
 =S(p) (59)

where we inserted from Eq. (17).

According to Eqs. (53) we find for the case of slab Alfven waves discussed in Sect. 3.4 for the ratio

 3Kpz,1U1κzz,1=pU1v[2S1S0+3U1v]
 =pU1v[ZR(s,σ+,σ−,Hc)+3U1v] (60)

with the helicity dependent function

 R(s,σ+,σ−,Hc)=(2−s)(4−s)(3−s)(5−s)[σ++σ−+Hc(σ+−σ−)] (61)

The value of the bracket of this function is not greater than for all helicity values, so the function is limited to values

 |Rmax|≤2(2−s)(4−s)(3−s)(5−s), (62)

which for all values of