Kinetic approaches to particle acceleration at cosmic ray modified shocks

0.4 MB·English

by  Elena Amato

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed1February2008 (MNLATEXstylefilev2.2)
Kinetic approaches to particle acceleration at cosmic
ray modified shocks
8
0 Elena Amato1⋆, Pasquale Blasi1 , and Stefano Gabici2
0 † ‡
2 1INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi, 5, 50125, Firenze, Italy

n
a 2Max-Planck-Institut fur Kernphysik, Heidelberg, Germany
J

9

Accepted —-.Received—–
]
h
p
-
o
r
t ABSTRACT
s
a
Kinetic approaches provide an effective description of the process of particle
[

2 accelerationat shock fronts and allow to take into account the dynamical re-
v
actionoftheacceleratedparticlesaswellastheamplificationoftheturbulent
3
2
magneticfieldasduetostreaminginstability.Thelatterdoesinturnaffectthe
7
3 maximum achievable momentum and thereby the acceleration process itself,
.
5
inachainofcausalitywhichis typicalofnon-linearsystems.Differentkinetic
0
7 approaches are characterized by different levels and types of approximations
0
: that also imply different computational times. Here we present the results
v
Xi of two such approaches: one which is mathematically rigorous but rather de-
r manding from the point of view of computational time, and the other which
a
is computationally very fast but based on an ansatz that, while physically
justified, is not rigorous.The identification of possible differences canbe cru-
cial in assessing the possibility of implementation of one such calculation in
hydrodynamical codes for supernova explosions. Special emphasis is given to
a discussion of the appearance of multiple solutions in both approaches.
Key words: acceleration of particles - shock waves
⋆ E-mail:[email protected]
† E-mail:[email protected]
‡ E-mail:[email protected]
2 E. Amato, P. Blasi and S. Gabici
1 INTRODUCTION
Recent observations of non-thermal radiation from shell-type supernova remnants (SNRs)
are showing that effective particle acceleration takes place in these astrophysical objects (see
Funk(2007)forarecent review).Inparticularthedetectionofspatiallyresolvednon-thermal
X-raysfromthinregionsclose totheforwardshock inseveral SNRshasshowed thefirst clear
evidence forstrongmagneticfieldamplificationascouldbeexpected iftheshockisefficiently
accelerating cosmic rays. The recent detection of TeV gamma rays (Aharonian et al. (2004,
2006a,b, 2005, 2006c)) is likely to provide further information on the magnetic field amplifi-
cation by adding to multifrequency observations and by possibly allowing the measurement
of the shape of the TeV gamma ray spectrum.
In the scenario with large magnetic field, low fluxes of gamma radiation through inverse
Compton scattering (ICS) of electrons are expected, thereby leaving a hadronic origin of the
gamma ray emission as a more likely possibility. This latter conclusion is however of more
limited strength and requires further confirmation through a continuous effort to detect
SNRs in gamma rays, not only in the TeV region but also in the GeV energy range, as
should become possible after the upcoming launch of the GLAST gamma ray telescope.
X-rays observed in the 1-10 keV energy range from SNRs are generated by synchrotron
emission of relativistic electrons and the thickness of the brightness profiles allows one to
estimate the total field in the shock vicinity to be of order 100-500 µG in basically all cases
in which measurements exist (V¨olk, Berezhko & Ksenofontov (2005)). This appears to be
the strongest evidence for efficient cosmic ray acceleration at SNR shocks, where here we
use the word efficient to indicate that an appreciable fraction of the kinetic pressure at the
shock is transformed into accelerated particles, which in turn change the dynamics of the
shock. In this regime the standard test particle theory (TPT) fails and a non-linear theory
is required (see Malkov & Drury (2001) for a review). From the phenomenological point of
view the introduction of the non-linear effects is crucial.
The reaction of the accelerated particles onto the shock structure has been calculated
within different approaches, both semi–analytical (see Malkov & Drury (2001) for a recent
review) and numerical (e.g. Ellison, Baring & Jones (1996); Kang & Jones (2006) and ref-
erences therein). Semi-analytical kinetic approaches provide a detailed description of the
process of particle acceleration at cosmic ray modified shocks, including particle spectra,
Kinetic approaches to shock acceleration 3
modified shock dynamics and thermodynamics, and have recently been generalized to in-
clude magnetic field amplification (Amato & Blasi (2005, 2006)).
The main limitation of these approaches is the fact that they all solve the stationary
problem of acceleration. This problem is simply ill defined in the case of acceleration of
protons (or nuclei) in sources such as supernova remnants: in fact in these cases there is
no appreciable energy loss process, and stationarity cannot be reached. This is in principle
true even in test particle approaches, but there the problem is limited to momenta close to
the maximum momentum, while the shape of the spectrum at lower momenta is rather well
defined and stationary. In the non-linear regime, the entire particle spectrum changes as a
consequence of the time dependence of the maximum momentum, therefore the assumption
of stationarityis inprinciple not well justified, thoughit istacitly assumed thatthe situation
to describe is one of quasi-stationarity.
Several approaches to time-dependent shock acceleration at cosmic ray modified shocks
have been put forward (Bell (1987); Falle and Giddings (1987); Kang & Jones (2006)), but
allofthemarebasedonnumericalapproaches(e.g.finitedifferencesolutionoftheequations).
In particular, Falle and Giddings (1987) found that quasi-stationary solutions are rather
good approximations to the full time-dependent solutions. Despite this, a generalization
of the semi-analytical approaches to include time dependence would be highly desirable,
and would probably help to shed some light on the problem of the appearance of multiple
solutions, as discussed below.
The first kinetic model that provided a semi-analytical (stationary) solution of the prob-
lem was developed by Malkov (1997) and Malkov, Diamond & V¨olk (2000) for strongly
modified shocks with a given, spatially constant, diffusion coefficient. A simple, approxi-
mate kinetic approach was later proposed by Blasi (2002, 2004). This approach provides
an accurate description of the shock modification and its effects on particle acceleration
provided the diffusion coefficient has a sufficiently strong dependence on momentum p (the
underlying assumption is somewhat similar to the starting assumption adopted by Eichler
(1979)).
It also allows to obtain the solution to the problem in a very short computational time,
which is particularly important when the particle acceleration process must be described
in the context of complex and time consuming hydrodynamical simulations (for instance
in the case of SNRs). Blasi, Gabici & Vannoni (2005) showed, among other things, that
the solutions obtained with Malkov’s method and those of Blasi’s method are in very good
4 E. Amato, P. Blasi and S. Gabici
agreement for Bohm diffusion coefficient, while the agreement becomes gradually worse
for Kraichnan (D(E) E1/2) and Kolmogorov (D(E) E1/3) diffusion coefficients, as
∝ ∝
expected. Even in these two cases, however, the differences were not very large and the
model of Blasi (2002, 2004) still provided an acceptable description of the main physical
aspects of the problem.
Morerecentlyasolutionofthesystemofequationsdescribingthetransportofaccelerated
particles and the dynamics of the shock was found by Amato & Blasi (2005), where the
diffusion coefficient could formally be taken as an arbitrary function of both momentum
and spatial coordinate (although solutions were actually computed for a few specific choices
of D(x,p)). Contrary to the simple approach of Blasi (2002) where a physically reasonable
ansatz was used to simplify the problem, this approach is mathematically rigorous but more
expensive from the computational point of view.
In a second paper by Amato & Blasi (2006) the self-generation of the waves and the
determination of the resulting diffusion coefficient, though in the context of quasi-linear
theory, were introduced. The method for taking into account magnetic field amplifica-
tion during the acceleration process at a modified shock, as assessed in this work, also
allowed for the determination of the maximum momentum in these complex circumstances
(Blasi, Amato & Caprioli (2007)).
Since the two methods introduced above aim at describing the same problem but using
a different set of approximations and, very important, with a quite different computational
cost, we wish to assess here the advantages and disadvantages of using one or the other.
While doing so, we also discuss an important feature which is found in both approaches,
and more generally in semi-analytical kinetic approaches (as well as in stationary two-fluid
models Drury & V¨olk (1980, 1981)), namely the appearance of multiple solutions.
Weconcludethatformostapplicationsrelatedtothedescriptionofthephenomenologyof
particleaccelerationinsupernovaremnantsthemodelofBlasi(2002,2004)canbesuccesfully
applied despite its simplifications, provided the diffusion coefficient scales rapidly enough
with momentum. This is certainly the case for Bohm diffusion. Moreover we complete this
simple model with recipes that allow us to determine the level of magnetic field amplification
and the maximum energy of the accelerated particles.
The present paper is structured as follows: in 2 we briefly introduce the two kinetic
§
approaches of Blasi (2002), model A hereafter, and Amato & Blasi (2005), model B. In
the subsections we discuss in detail the results of the two theoretical approaches. In 3 we
§
Kinetic approaches to shock acceleration 5
describe the introduction in the first approach of a recipe for the acceleration time which
allows us to estimate the maximum momentum of the accelerated particles. We conclude in
4.
§
2 TWO KINETIC APPROACHES
Inthis section webriefly discuss two kinetic approaches, the oneofBlasi(2002,2004) (Model
A) that is expected to work for diffusion coefficients D(p) that change rapidly enough with
the particle momentum, and the one proposed by Amato & Blasi (2005) (Model B) and
valid for any choice of D(p,x). The main point of the section is to illustrate the results
obtained with the two approaches, identifying the advantages and possible shortcomings of
the first, which is being implemented in some computationally heavy hydrodynamical codes,
in which the formal solution of Model B would be hardly possible to introduce due to the
much longer computational time required.
The basic equations are the transport equation:
∂ ∂ ∂f(x,p) 1 du ∂f(x,p)
D(x,p) f(x,p) u + p +Q(x,p) = 0, (1)
∂x " ∂x #− ∂x 3 dx! ∂p
and the conservation equation for the total momentum:
1 1
ξ (x) = 1+ U(x) U(x)−γg. (2)
c γ M2 − − γ M2
g 0 g 0
Here f(x,p) is the particle distribution function in the shock frame, u(x) is the velocity of
the background fluid, which equals u downstream and changes continuously upstream, from
2
u immediately upstream of the subshock to u at upstream infinity. The quantity U(x) =
1 0
u(x)/u is the normalized velocity, bound to equal unity at x . In Eq. 2, ξ (x) =
0 c
→ −∞
P (x)/ρ u2 is the cosmic ray pressure at the position x, normalized to the kinetic pressure
CR 0 0
ρ u2 at upstream infinity. It is customary to introduce the compression factor R = u /u
0 0 sub 1 2
at the subshock and the total compression factor R = u /u . Assuming homogenisation
tot 0 2
of the cosmic ray plasma in the downstream section (df/dx− = 0, a consequence of the
assumption of stationarity) one easily obtains (Blasi (2002)) that the distribution function
of the particles at the shock location is
3R ηn p dp′ 3R U (p′)
tot 0 tot p
f (p) = exp , (3)
0 RtotUp(p)−1! 4πp3inj (−Zpinj p′ RtotUp(p′)−1)
where we introduced the function U (p) = u /u , with
p p 0
1 0
u = u dx(du/dx)f(x,p) . (4)
p 1
− f (p)
0 Z−∞
6 E. Amato, P. Blasi and S. Gabici
This result is very general and is not related to the assumptions of the model of Blasi (2002).
Here we assumed, as usual, that the injection is a delta function in space (at the shock) and
inmomentum (attheinjectionmomentum p )intheformQ(x,p) = ηngas,1u1δ(p p )δ(x),
inj 4πp2inj − inj
with n = n R /R the gas density immediately upstream (x = 0−) and η the fraction
gas,1 0 tot sub
of the particles crossing the shock which are going to take part in the acceleration process.
The compression factors R and R , if the gas in the upstream frame evolves adiabat-
sub tot
ically, are related through the following expression:
1
2 (γ +1)Rγg (γ 1)Rγg+1 γg+1
R = Mγg+1 g sub − g − sub . (5)
tot 0 " 2 #
Model A
This model was proposed by Blasi (2002) and put in a form to include the reacceleration
of seed particles by Blasi (2004). It is characterized by a simple implementation that makes
it fast in terms of numerical computation while keeping all the main physical ingredients of
the problem. For these reasons it has been implemented in hydrodynamical codes in order
to calculate the instantaneous spectrum of accelerated particles in an expanding supernova
shell (Ellison et al. (2007)).
The essence of the approach is based on the observation that the function f(x,p) in Eq. 4
isexpected tosufferanexponential suppression atadistancex suchthat 0 dxu(x)/D(p)
p xp ∼
1. The spatial position where this suppression is found in general dependRs on momentum if
the diffusion coefficient depends on momentum. These spatial locations for different values
of p are well defined if D(p) is a strong function of momentum, as was first pointed out
by Eichler (1979), and become more ill defined for D(p) weakly dependent on p. As a first
approximation we may assume that the distribution function is f(x,p) = f (p) at x 6 x
0 p
| | | |
and vanishes at x > x . With this assumption Eq. 4 gives u u(x ). The distance x has
p p p p
| | | | ≈
the meaning of the typical distance to which particles with momentum p can diffuse away
from the shock in the upstream fluid, and the fluid velocity at that point is u(x ) u . It
p p
≈
follows that Eq. 2 can be transformed into an equation which depends only on the particle
momentum p:
1 1
ξ (p) = 1+ U U−γg, (6)
c γ M2 − p − γ M2 p
g 0 g 0
where
Kinetic approaches to shock acceleration 7
4π pmax
ξ (p) dpp3v(p)f (p), (7)
c ≈ 3ρ u2 0
0 0 Zp
and v(p) is the velocity of particles with momentum p. Differentiating Eq. 6 with respect to
momentum we finally have:
dU 1 4π
p p 1 U−(γg+1) = p4v(p)f (p). (8)
dp " − M02 p # 3ρ0u20 0
For p p we have U u /u = R /R , while for p p one has U 1. The
inj p 1 0 sub tot max p
→ → → →
procedure tocalculatethesolutionistherefore straightforward:foragiven guess ofR (and
sub
therefore of R through Eq. 5) we can solve the differential equation Eq. 8 coupled with the
tot
expression for f (p) as a function of U (Eq. 3) in an iterative way. In general we end up with
0 p
U (p ) = 1, which implies that the given value of R is not a solution of our problem.
p max sub
6
The solution corresponds to that value of R (and R ) for which U (p ) = 1. This value
sub tot p max
also identifies completely the function U and the distribution function at the shock f (p),
p 0
which is the sought after result. It is important to notice that this procedure does not lead
only to the determination of f (p) but also to the values of the thermodynamical quantities
0
at the shock, such as the temperature of the gas upstream and downstream. What Model
A cannot provide is the spatial dependence of the quantities in the precursor, although we
provide below a physical ansatz that may indirectly make this piece of information available.
We later check the results by comparing them with the outcome of Model B.
In both approaches that we describe here we adopt the injection recipe known as thermal
leakage (Blasi, Gabici & Vannoni (2005); Gieseler et al. (2000)) which connects the value of
η to that of R in a unique way through the relation:
sub
4
η = (R 1)ξ3e−ξ2. (9)
3π1/2 sub −
Here ξ is a parameter that identifies theinjection momentum asa multiple ofthe momentum
of the thermal particles in the downstream section (p = ξp ). The latter is an output of
inj th,2
the non-linear calculation, since we solve exactly the modified Rankine-Hugoniot relations
together with the cosmic rays’ transport equation. For the numerical calculations that follow
we typically use values of ξ between 3 and 4, which correspond to fractions of order
∼
10−4 10−5 of the particles crossing the shock to be injected in the accelerator.
−
Model B
A formal solution of the transport equation and the conservation equations for an ar-
8 E. Amato, P. Blasi and S. Gabici
bitrary choice of the diffusion coefficient in both its dependence on momentum and spatial
coordinate was found by Amato & Blasi (2005).
Malkov (1997) and Amato & Blasi (2005) showed that an excellent approximation to the
solution f(x,p) has the form
q(p) u 0 u(x′)
f(x,p) = f (p)exp (1 1) dx′ , (10)
0 "− 3 − u2 Zx D(x′,p)#
where q(p) = dlnf0(p) is the local slope of f (p) in momentum space. In fact the factor
− dlnp 0
1 u1 in the argument of the exponential was introduced by Blasi, Amato & Caprioli (2007)
−u2
in order to satisfy exactly the boundary condition at the shock, even for weakly modified
shocks.
In terms of the distribution function (Eq. 10), we can also write the normalized pressure
in accelerated particles as:
4π pmax 0 U(x′)
ξ (x) = dp p3v(p)f (p)exp dx′ , (11)
c 3ρ0u20 Zpinj 0 "−Zx xp(x′,p)#
where for simplicity we introduced x (x,p) = 3D(p,x).
p q(p)u0
By differentiating Eq. 11 with respect to x we obtain
dξ
c
= λ(x)ξ (x)U(x), (12)
c
dx
where
pmaxdp p3 1 v(p)f (p)exp 0dx′ U(x′)
λ(x) =< 1/x > = pinj xp(x,p) 0 − x xp(x′,p) , (13)
p ξc R pmax dp p3v(p)f (p)exp h 0dRx′ U(x′) i
pinj 0 − x xp(x′,p)
h i
and U(x) is expressed aRs a function of ξ (x) througRh Eq. 2.
c
Finally, after integration by parts of Eq. 4, one is able to express U (p) in terms of an
p
integral involving U(x) alone:
0 1 0 U(x′)
U (p) = dx U(x)2 exp dx′ , (14)
p Z−∞ xp(x,p) "−Zx xp(x′,p)#
which allows one to easily calculate f (p) through Eq. 3.
0
Eqs. 2 and 12 can be solved by iteration in the following way: for a fixed value of the
compression factorat the subshock, R ,the valueof thedimensionless velocity at theshock
sub
is calculated as U(0) = R /R . The corresponding pressure in the form of accelerated
sub tot
particles is given by Eq. 2 as ξ (0) = 1 + 1 Rsub 1 Rsub −γg. This is used as a
c γgM02 − Rtot − γgM02 Rtot
(cid:16) (cid:17)
boundary condition for Eq. 12, where the functions U(x) and λ(x) (and therefore f (p)) on
0
the right hand side at the kth step of iteration are taken as the functions at the step (k 1).
−
In this way the solution of Eq. 2 at the step k is simply
Kinetic approaches to shock acceleration 9
0
ξ(k)(x) = ξ (0)exp dx′λ(k−1)(x′)U(k−1)(x′) , (15)
c c −
(cid:20) Zx (cid:21)
with the correct limits when x 0 and x . At each step of iteration the functions
→ → −∞
U(x), f (p), λ(x) are recalculated (through Eq. 2, Eqs. 14 and 3, and Eq. 13, respectively),
0
until convergence is reached. The solution of this set of equations, however, is also a solution
of our physical problem only if the pressure in the form of accelerated particles as given
by Eq. 2 coincides with that calculated by using the final f (p) in Eq. 6. This occurs for a
0
value of R , which fully determines the solution of our problem for an arbitrary diffusion
sub
coefficient as a function of location and momentum.
2.1 Spectra and velocity profiles
All approaches to particle acceleration at modified shocks predict the formation of a pre-
cursor in the upstream region, resulting in a gradient of the velocity profile of the fluid.
Since qualitatively the spectrum of the accelerated particles is still determined by an effec-
tive compression factor felt by the particles of given momentum, and the velocity in the
precursor increases with the distance from the shock, it is easy to infer that the spectrum
of the accelerated particles is not expected to be a power law and more precisely that it
should be concave (steeper at low energies and flatter at high energies). Here we discuss the
detailed shape of the spectrum at the shock as obtained through the two kinetic approaches
described above. In Fig. 1 (left panel) we plot the spectra as a function of the momentum of
particles for model B (solid lines) and for model A (dashed lines). The curves are obtained
for p = 105m c, u = 5 108cm s−1, ξ = 3.5 and for the values of the Mach number
max p 0
×
at upstream infinity M = 10, 100, 1000 (curves labeled as 1, 2 and 3, respectively). Bohm
0
diffusion is assumed. The agreement between the two sets of curves is excellent for relatively
low Mach numbers (M 10) and remains good even up to much larger Mach numbers, and
0
∼
in fact for all values we have tried. The largest discrepancies between the two methods are
at the level of 20%. The reason for such discrepancies is to be found in the assumption
∼
that U (p ) is required to equal unity in the approach of Model A.
p max
The velocity profile of the fluid in the precursor is plotted in Fig. 1 (right panel) for the
two models (again solid and dashed lines respectively).
On the x-axis we plot the distance x from the shock in the upstream region in units of
10 E. Amato, P. Blasi and S. Gabici
Figure 1. Left Panel: Spectra of accelerated particles for pmax = 105mpc, u0 = 5×108cm s−1, ξ = 3.5 and for M0 =
10,100,1000 (curves labeled as 1, 2 and 3 respectively). Bohm diffusion coefficient is adopted. The solid lines are obtained
with the calculation of Amato&Blasi (2005) (model B), the dashed ones with that of Blasi (2002) (model A). Right Panel:
Velocity profiles in the precursor for the cases in the left panel. Again the solid curves are computed with Model B and the
dashedoneswithModelA.Thespatialcoordinateisinunitsofxp∗,with|xp∗|=D(pmax)/(u0Up(pmax))asdiscussedinthe
text.
x , which is defined as
p∗
D(p )
max
x = .
p∗
| | u U (p )
0 p max
Some comments are required on the calculation of U(x) for Model A. As discussed in the
previous section, this model does not keep any information about the spatial dependence
of the quantities in the precursor, although such information is somehow contained in the
relation between a momentum p and the mean diffusion length of particles with such mo-
mentum, x(p) D(p)/u (p). The dashed lines in Fig. 1 (right panel) are obtained in the
p
| | ≈
following way: for a given location x upstream, the equation x = D(p)/u (p) is inverted
p
and a corresponding value p of the minimum momentum of particles that may have diffused
to the point x is obtained. At this point the velocity U(x) (in units of u ) is by definition
0
U (p) for the value of p corresponding to x. By definition the fluid velocity in the simple
p
model is bound to be unity at x/x = 1 because no particles are supposed to be able to
p
reach farther regions. In the exact solution there is a spread in the distances that can be
diffusively reached at given momentum and the transition to U(x) = 1 is smoother. This
difference in the velocity profile affects mainly the results for the spectrum at p p , but
max
∼
since these particles carry an appreciable amount of energy in the case of modified shocks,
the whole spectral shape is somewhat affected (at the level of at most 20% in the strongly
∼
modified cases).
As discussed above, the ingredient of model A that makes it much faster in terms of com-
putational time consists in the assumption that particles with given momentum all travel to