Cosmic ray acceleration at blast waves from type Ia supernovae

0.4 MB·English

by  Hyesung Kang

Journal of The Korean Astronomical Society
39: 1 11, 2006
∼
COSMIC RAY ACCELERATION AT BLAST WAVES FROM TYPE Ia
SUPERNOVAE
HYESUNG KANG
DepartmentofEarthSciences, PusanNationalUniversity,Pusan609-735,Korea
E-mail: [email protected]
(Received November 10, 2006; Accepted December 6, 2006)
ABSTRACT
7
We have calculated the cosmic ray (CR) acceleration at young remnants from Type Ia supernovae expanding
0
into a uniform interstellar medium (ISM). Adopting quasi-parallel magnetic fields, gasdynamic equations and the
0
diffusion convection equation for the particle distribution function are solved in a comoving spherical grid which
2
expands with the shock. Bohm-type diffusiondue to self-excitedAlfv´enwaves,drift and dissipationofthese waves

n in the precursorandthermal leakage injection were included. With magnetic fields amplified by the CR streaming
a instability, the particle energycan reachup to 1016Z eV atyoung supernovaremnants (SNRs) of severalthousand
J
years old. The fraction of the explosion energy transferred to the CR component asymptotes to 40-50 % by that

2 time. For a typical SNR in a warm ISM, the accelerated CR energy spectrum should exhibit a concave curvature
with the power-law slope flattening from 2 to 1.6 at E > 0.1 TeV.
1 ∼
v Key Words : cosmic ray acceleration – supernova remnants – hydrodynamics – methods:numerical
7
2
0
1 I. INTRODUCTION equationfortheparticlemomentumdistribution,f(p),
0 issolvedalongwithsuitablymodifiedgasdynamicequa-
7 Diffusive shock acceleration (DSA) is widely ac- tions (e.g., Kang & Jones1991). Accurate solutionsto
0 cepted as the primary mechanism through which cos- the CR diffusion-convection equation require a com-
/ mic raysaregeneratedinawide rangeofastrophysical putational grid spacing significantly smaller than the
h
p shocks (Drury, 1983; Malkov & Drury 2001 and ref- particle diffusion length, ∆x xd(p) = κ(p)/us,
- erences therein; Kang & Jones 2002; Kang 2003). It where κ(p) is diffusion coefficien≪t and us is the shock
o is well known that the CR energy density is compara- speed. In a realistic diffusion transportmodel, the dif-
r ble to the gas thermal energy density in the ISM and fusion coefficient has a steep momentum dependence,
t
s playsimportantdynamicalrolesintheevolutionofour κ(p) ps, with s 1 2. So a wide range of length
a Galaxy. Most of galactic cosmic rays, at least up to scales∝is required∼to b−e resolved in order to follow
v: 1014 eV of the proton energy, are believed to be accel- the CR acceleration from the injection energy (typ-
i eratedbySNRswithinourGalaxyviaFermifirstorder ically p /m c 10−2) to highly relativistic energy
X inj p
process(Blanford&Ostriker1978;Lagage&Cesarsky (p/m c 1). T∼his constitutes an extremely challeng-
p
r 1983; Blandford & Eichler 1987; Drury et al. 2001). ing num≫erical task, requiring rather extensive compu-
a
Time-dependent, kinetic simulations of the CR ac- tational resources.
celerationatSNRshaveshownthatupto50%ofexplo- To overcome this numerical problem, Berezhko and
sion energy can be converted to CRs, when a fraction collaborators (e.g., Berezhko et al. 1995) introduced a
10−4 10−3 of incoming thermal particles are injected “change of variablestechnique” in which the radialco-
−
intotheCRpopulationatthesubshock(e.g.,Berezhko, ordinate is normalized by the diffusion length, x (p),
d
Ksenofontov, & Yelshin 1995; Berezhko, & Vo¨lk 1997; at each particle momentum for the upstream region.
Vo¨lk & Berezhko 2005; Kang & Jones 2006). This This allowed them to solve the coupled system of gas-
shouldbeenoughtoreplenishthegalacticCRsescaping dynamicequationsandtheCRtransportequationwith
fromourGalaxywithLCR 1041erg s−1. X-rayobser- κ(p) p. Theirschemewasdesignedforsimulationsof
vationsofyoungSNRssuch∼asSN1006andCasAindi- super∝novaremnants,whichwererepresentedbypiston-
catethepresenceofTeVelectronsemittingnonthermal drivensphericalshocks inone-dimensional(1D) spher-
synchrotronemissionimmediatelyinsidetheouterSNR ical geometry.
shock (Koyama et al. 1995; Bamba et al. 2003). They
On the other hand, Kang and collaborators have
provide clear evidences for the efficient acceleration of
taken an alternative approach that is based on a more
the CR electrons at SNR shocks. Also recent HESS
conventional Eulerian formalism. Adaptive Mesh Re-
observationof SNR RXJ1713.7-3946indicates possible
finement (AMR) technique andsubgrid shock tracking
detections of π0 decay γ rays from the hadronic CRs
technique were combined to build CRASH (Cosmic-
accelerated by the SNR shock (Aharonian et al. 2004;
Ray Amr SHock) code in 1D plane-parallel geometry
Berezhko & Vo¨lk 2006).
(Kang et al. 2001) and in 1D spherical symmetric ge-
Inthe kinetic equationapproachto numericalstudy ometry (Kang & Jones 2006). In order to implement
of CR acceleration at shocks, the diffusion-convection the shock tracking and AMR techniques effectively in
– 1 –
2 KANG
a spherical geometry, we solve the fluid and diffusion- the upstream region is represented by the term
convection equations in a frame comoving with the
outer spherical shock. Adopting a comoving frame ∂Pc
W(r,t)= v , (5)
A
turns out to be a great numerical success, since we − ∂r
can achieve numerical convergence at a grid resolution
where v = B/√4πρ is the Alfv´en speed. This term
muchcoarserthanthatrequiredinanEuleriangrid. In A
derives froma simple modelin whichAlfv´en wavesare
thecomovinggrid,theshockremainsatthesameloca-
amplified by streaming CRs and dissipated locally as
tion, so the compression rate is applied consistently to
heat in the precursor region (e.g., Jones 1993).
the CRdistributionatthe subshock,resultinginmuch
moreaccurateandefficientlowenergyCRacceleration. TheCRpopulationisevolvedbysolvingthediffusion-
convection equation,
Inthepresentpaper,weapplythesphericalCRASH
code for the CR acceleration at remnant shocks from
∂g ∂g 1 ∂ ∂g
Type Ia SNe expanding into a uniform interstellar +(u+u ) = r2(u+u ) 4g
medium, assuming a quasi-parallelfield geometry. De- ∂t w ∂r 3r2∂r w (cid:18)∂y − (cid:19)
(cid:2) (cid:3)
tails ofthe numericalmethod aredescribedin II. The 1 ∂ ∂g
§ + r2κ(r,y) , (6)
simulation results are presented and discussed in III, r2∂r (cid:20) ∂r(cid:21)
§
followed by a summary in IV.
§
where g =p4f, with f(p,r,t) the pitch angle averaged
II. NUMERICAL METHOD CR distribution, and y = ln(p), while κ(r,y) is the
diffusion coefficient parallel to the field lines (Skilling
(a) BASIC EQUATION
1975). For simplicity we express the particle momen-
Here we consider the CR acceleration at a quasi- tum, p in units mpc hereafter and consider only the
parallel shock where the magnetic field lines are par- proton CR component. The wave speed is set to be
allel to the shock normal. So we solve the standard uw = vA in the upstream region, while we use uw = 0
gasdynamic equations with CR pressure terms added in the downstream region. This term reflects the fact
in the Eulerian formulation for one dimensional spher- that the scatteringby Alfv´enwavestends to isotropize
ical symmetric geometry. the CR distribution in the wave frame rather than the
gas frame.
∂ρ ∂ 2
+ (ρu)= ρu, (1) (b) Spherical CRASH code
∂t ∂r −r
Details of the CRASH code in 1D spherical sym-
∂(ρu) + ∂ (ρu2+P +P )= 2ρu2, (2) metricgeometrycanbefoundinKang&Jones(2006),
∂t ∂r g c −r so we briefly describethe basic featureshere. We solve
theequations(1)-(6)inacomovingframethatexpands
with the instantaneous shock speed, since a spherical
∂(ρeg) ∂ ∂Pc 2 shockcanbemadetobestationarybyadoptingcomov-
+ (ρe u+P u)= u (ρe u+P u),
∂t ∂r g g − ∂r − r g g ing variables which factor out a uniform expansion or
(3) contraction. Because the shock is at rest and tracked
accurately as a true discontinuity, we can refine the
region around the gas subshock at an arbitrarily fine
∂S ∂ 2 (γg 1) level. The AMR technique allows us to “zoom in” in-
+ (Su)= Su+ − [W(r,t) L(r,t)],
∂t ∂r −r ργg−1 − side the precursor structure with a hierarchy of small,
(4) refined grid levels applied around the shock. The re-
where P and P are the gas and the CR pressure, sultisanenormoussavingsinbothcomputationaltime
g c
respectively,e =P /[ρ(γ 1)]+u2/2is the totalen- and data storage over what would be required to solve
g g g
− theproblemusingmoretraditionalmethodsonasingle
ergyofthe gasper unit mass. The evolutionofa mod-
ified entropy, S = Pg/ργg−1, is followed everywhere fine grid.
except across the subshock, since for strongly shocked
(c) The Thermal Leakage Injection Model
flows numerical errors in computing the gas pressure
fromthetotalenergycanleadtospuriousentropygen- Theinjectionratewithwhichsuprathermalparticles
erationwith standardmethods, especially in the shock are injected into CRs at the subshock depends in gen-
precursor(Kang,Jones,&Gieseler2002). Totalenergy eralupontheshockMachnumber,fieldobliquityangle,
conservation is applied only across the subshock. The and strength of Alfv´en turbulence responsible for scat-
remaining variables, except for L and W, have stan- tering. The CRASH codes treat this process naturally
dardmeanings. The injectionenergylossterm,L(r,t), and self-consistently via “thermal leakage” through
accounts for the energy carried by the suprathermal lowest momentum bins. The thermal leakage injec-
particles injected into the CR component at the sub- tionmodelemulatestheprocessbywhichsuprathermal
shock. Gas heating due to Alfv´en wave dissipation in particles wellinto the tailofthe postshockMaxwellian
Cosmic Ray Acceleration at SNRs 3
Fig. 1.— Top two panels show the evolution of a typical SNR in a gasdynamic simulation. Bottom two panels show the
positions of the outer forward shock (Rout) and inner reverse shock (Rin) and the speed of the outer shock (Uout) as a
function of time. The Sedov Taylor similarity solutions for the outershock, RST and UST are also shown for comparison.
distributionleakupstreamacrossaquasi-parallelshock ǫ = 0.16 0.2 in this study. Dependence of the CR
B
−
(Malkov&Vo¨lk1998;Malkov&Drury2001). Thisfil- injectionandaccelerationonthisparameterwillbedis-
tering process is implemented numerically by adopting cussed below.
a “transparency function”, τ (ǫ ,υ), that expresses
esc B
the probability of supra-thermal particles at a given (d) A Bohm-like Diffusion Model
velocity, υ, leaking upstream through the postshock
Self-excitationofAlfv´enwavesbytheCRstreaming
MHD waves (Kang et al. 2002). One free parameter
instability in the upstreamregionis anintegralpartof
cisonthtreolisnvtehrissefruantciotioonf;thneamamelpy,litǫuBde=ofBt0h/eBp⊥o,stwshhoicchk the DAS at SNRs (Bell 1987; Vo¨lk et al. 1988; Lucek
MHDwaveturbulenceB⊥ tothegeneralmagneticfield &Bell2001). Theparticlesareresonantlyscatteredby
aligned with the shock normal, B (Malkov & Vo¨lk those waves, diffuse across the shock, and get injected
0
1998). Plasmahybridsimulationsandtheorybothsug- into the Fermi first-order process. These complex in-
gest that 0.25 < ǫ < 0.35, so that the model is well teractions are represented by the diffusion coefficient,
B
whichisexpressedintermsofameanscatteringlength,
constrained. H∼oweve∼r, such large values of ǫ lead to
B
λ, and the particle speed, υ, as κ(x,p) = λυ/3. The
veryefficientinitialinjectionandmostoftheshocken-
Bohm diffusion model is commonly used to represent
ergy is quickly transferred to the CR component for
a saturated wave spectrum (i.e., λ = r , where r is
strong shocks considered here (50 < M < 300 ini- g g
s
tially), causing a numerical problem at the very early the gyro-radius), κB(p) = κnp2/ (p2 + 1)1/2. Here
stage of simulations. So we adopted smaller values, κn = mc3/(3eB) = 3.13 × 1022cm2s−1Bµ−1, and Bµ
is the magnetic field strength in units of microgauss.
4 KANG
Table 1. Model Parameters
Model n E B ǫ r t u P
ISM o µ B o o o o
(cm−3) (1051 ergs) µG (pc) (years) (104 km s−1) (10−6erg cm−3)
S1 0.3 1 30 0.16 3.19 255. 1.22 1.05
S2 0.3 4 30 0.16 3.19 127. 2.45 4.20
S3 0.3 1 5 0.16 3.19 255 1.22 1.05
S4 0.003 1 30 0.16 14.8 1182. 1.22 1.05 (-2)
S5 0.3 1 5 0.2 3.19 255 1.22 1.05
Because of the steep momentum dependence for non- of the outer shock speed, U (t) can be approximated
out
relativistic particles (p 1), simulations with a Bohm by the Sedov solution U (t/t )−2/5 for t/t > 0.2.
diffusion model require≪extremely fine grid resolution In order to take accountSTfo∝r the CoR acceleratioon from
around the shock where freshly injected CRs are con- free expansion stage through ST stage, we begin the
centrated. Instead we adopt a Bohm-like diffusion co- calculations with the ST similarity solution at t/t =
o
efficientthatincludesaweakernon-relativisticmomen- 0.5. In principle we could start the simulations from
tum dependence, t/t 0.1. Such simulations, however, require very
o
∼
long computational time, since we need to include the
ρ
0 extremely hot region with fast sound speeds near the
κ(r,p)=κ p . (7)
n
· (cid:18)ρ(r)(cid:19) explosioncenter. Wecarriedoutonemodelfromt/t =
o
0.2 and compared it with the case started from t/t =
o
Previous studies showed that simulations using these 0.5 in the next section. Since the total CRenergy gain
two types of diffusion coefficient produced very simi- is proportional to the kinetic energy passed through
lar results (Kang et al. 2001). The assumed density the shock, E = 2π ρ U3 R2 dt, we expect our
sw 0 out out
dependence for κ accounts for compression of the per- calculations should capRture the key aspects of the CR
pendicular component of the wave magnetic field and acceleration at the outer SNR shock.
also inhibits the acoustic instability in the precursor
We terminated the simulation at t/t = 10, while
of highly modified CR shocks (Kang, Jones, & Ryu, o
theSTstageendstypicallywhentheshockslowsdown
1992). Hereafter we use the subscripts ’0’, ’1’, and ’2’ to U <300 km s−1 around t/t 3000. However, it
to denote conditions far upstream of the shock, imme- ST o ∼
has been shown that the highest momentum, p , is
max
diately upstream of the gas subshock and immediately
achievedbytheendofthefreeexpansionstageandthe
downstream of the subshock, respectively.
transfer of explosion energy to the CR component oc-
curs mostly during the early ST stage (e.g., Berezhko
III. Simulations of Sedov-Taylor Blast Waves
et al. 1997). We will show in the next section that
the shockpropertiesandtheCRaccelerationefficiency
Forasupernovaremnantpropagatingintoauniform
reachroughlytimeasymptoticvaluesfort/t >1. Also
ISM, the CR acceleration takes place mostly during o
itwassuggestedthatnon-linearwavedampingandthe
freeexpansionandSedov-Taylor(ST hereafter)stages,
wave dissipation due to ion-neutralcollisions may sup-
sincetheshockslowsdownsignificantlyafterward. Fig.
presstheMHDwavessignificantlyatthelateSTstage,
1 shows the evolution of a typical SNR calculated by
leading to fast particlediffusion andinefficient acceler-
a hydrodynamicscode without the CR pressureterms.
ation. (V¨olk et al. 1998; Ptuskin & Zirakashvili 2003)
This demonstrates that the ST solution is established
only after the inner reverse shock is reflected at the
IV. Results
center at t/t 7 (see below for the definition of nor-
o
∼
malization constants). Before that time, the reverse
(a) SNR Model Parameters
shock is strong and dynamically important. In our
simulations, however, we will ignore the reverse shock, WeconsideraTypeIasupernovaexplosionwiththe
because the current version of CRASH code can treat ejectamass,M =1.4M ,inawarmorhotISMwith
ej
only one shock. Applicationof ourAMR algorithmfor a uniform density. Mode⊙l parameters are summarized
multiplesphericalshocksisnotsimple,sinceitrequires in Table 1. The fiducial model, labeled S1 in Table 1,
multiple, comoving grids. The CR acceleration at the hastheexplosionenergy,E =1051 ergs,andtheback-
o
reverse shock is thought to be not important, because ground density, n =0.3 cm−3. The pressure of the
ISM
the kinetic energy passed through the reverse shock is background gas is set to be P 10−12 erg cm−3,
ISM
≈
relatively small. Also adiabatic losses by CRs acceler- which determines the soundspeed of the upstreamgas
atedearlyonintheinteriorandthenadvectedoutward and so the Mach number of the SNR shock. Recent
through the ST phase would generally be very large. X-ray observationsof young SNRs indicate a magnetic
On the other hand, Fig. 1 shows that the evolution field strengthmuchgreaterthanthe meanISM fieldof
Cosmic Ray Acceleration at SNRs 5
Fig. 2.— Evolution of model S1 SNR expanding into a uniform ISM at t/t =0.5, 1., 6., and 10. The model parameters
o
are Mej =1.4M⊙, Eo =1051 ergs, nISM = 0.3cm−3, and Bµ =30. The injection parameter for thermal leakage injection,
ǫ =0.2. The lower left panel shows the volume integrated CR spectrum, G(p)= f(r,p)p4r2dr. The initial condition at
B
t/to =0.5 (solid line) is set by theSedov-Taylorsimilarity solution. R
5µG, or values expected by compression of that field rate.
(e.g., Berezhko et al. 2003; Vo¨lk et al. 2005). It is be- The physical quantities are normalized, both in the
lievedthatthe magnetic fieldupstreamfromthe shock numericalcodeandintheplotsbelow,bythefollowing
isamplifiedbytheCRstreaminginstabilityinthepre- constants:
cursorregion(Bell 1978;Lucek & Bell 2001). Thus, to
represent this effect we take B = 30µG as the fiducial 3M 1/3
ej
field strength. The strength of magnetic field deter- ro = ,
(cid:18)4πρ (cid:19)
mines the size of diffusion coefficient, κ , and the drift o
n
speed of Alfv´en waves relative to the bulk flow. The ρ r5 1/2
AvAlf,0v´en=sp(1e.e8dkims sg−iv1e)nBµb/y√vnAISM=. vAT,0h(eρ/sρe0c)o−n1d/2mwohdeerle, to =(cid:18)uEo=oor(cid:19)/t ,,
S2, assumes a higher explosion energy, while the third o o o
model, S3, assumes the ISM magnetic field, rather ρo =(2.34 10−24gcm−3)nISM,
×
than the amplified field. Model S4 assumes a hot ISM P =ρ u2.
(T 106 K), while all other models assume a warm o o o
ISM≈(T 104 K). For models S1-S4 ǫB = 0.16 is These values are also given in Table 1 for reference.
≈
adopted. Model S5 has the same parameters as model
S3 except ǫ = 0.2, which allows a higher injection
B
6 KANG
Fig. 3.— Immediatepre-subshockdensity,ρ1, post-subshock density,ρ2,shock Mach number,Ms, post-subshock CR and
gaspressureinunitsoftherampressureoftheunmodifiedSedov-Taylorsolution,ρ0US2T ∝(t/to)−6/5,andtheCRinjection
parameter, ξ,are plotted for models S1 and S2.
(b) Remnant Evolution celeration during the early free expansion stage would
affect the results at late ST stage. The total CR en-
Fig. 2 shows the evolution of model S1 for t/t =
o ergy accelerated by t/t = 10 differs about 3 % and
0.5 10. The solid lines are for the initial condition o
− theCRspectrumextendstoslightlyhigherpmax inthe
which is assumed be the Sedov Taylor similarity solu-
simulationstartedfrom t/t =0.2,as one wouldantic-
tion (R /r = ξ (t/t )0.6, U /u = 0.6ξ (t/t )−0.4, o
ST o s o ST o s o ipate from the longer acceleration interval. The small
where ξ = 1.15167) extrapolated to t/t = 0.5. Ini-
s o differences indicate that omitting the evolution before
tially there is no pre-existing CRs and so all CR par-
t/t =0.5doesnotaffecttheoverallconclusionsofthis
ticles are freshly injected at the shock. The CR pres- o
work.
sure becomes dominated over the gas pressure and the
density compression across the total shock becomes
(c) CR Injection and Acceleration
ρ /ρ 7 after t/t 1. As the shock slowsdownand
2 0 o
the CR≈pressure incr∼eases, the subshock Mach number The efficiency of the particle injection is quantified
decreases and the thermal population cools down, re- by the fraction of particles swept through the shock
sultinginlessefficientparticleinjectionatlowenergies. that have been injected into the CR distribution:
Wenotethattheinnerboundaryofthesimulationgrid 4πr2dr 4πf (p,r,t)p2dp
moves out with the expanding comoving grid. ξ(t)= CR , (8)
R R4πr2n u dt
We repeated the same calculation starting from an s 0 s
earlier time, t/t = 0.2, to explore how the CR ac- where f is the CR dRistribution function, n is the
o CR 0
particle number density far upstream and r is the
s
Cosmic Ray Acceleration at SNRs 7
Fig. 4.—Pre-subshockdensity,ρ1,post-subshockdensity,ρ2,shockMachnumber,Ms,post-subshockCRandgaspressure
in unitsof theram pressureof theunmodifiedSedov-Taylorsolution, ρ0US2T ∝(t/to)−6/5, and theCR injection parameter,
ξ,areplotted for modelsS3and S4. Theleft panelsalso showmodel S5with ǫ =0.2 (longdashed linesand dottedlines).
B
shock radius. B =30µG. While the mode with higher B has smaller
Figs. 3-4showtheevolutionofshockpropertiessuch diffusioncoefficientandsotheCRspectrumextendsto
as the compression factors, subshock Mach number, higher pmax (see Fig. 5), the reduction of the injection
postshock pressures, and the injection fraction for all and accelerationdue to the wavedrift leads to less CR
fivemodels. Mostofthesequantitiesapproachtotime- energy.
asymptotic values. The shock Mach number is the key Model S5 has a larger value of ǫ =0.2 than model
B
parameter that determines the CR injection and ac- S3 does, but otherwise they have the same model pa-
celeration efficiency. Model S4 has much lower initial rameters. Weakerturbulence(largerǫ )leadtohigher
B
Machnumber,soitshowslessinjectionrate(ξ 10−5) thermal leakage, so the injection rate ξ is about 15
compared to the other models (ξ > 10−4). T∼he com- times larger than that of models S1-S3. This in turn
pression factor in the precursor, ρ∼/ρ 2 3, while results in slightly higher CR pressure and more signifi-
1 0
≈ −
the compressionfactor acrossthe totalshockstructure cant precursor compression. However, it demonstrates
varies somewhat, ρ /ρ 7.0, 8.2, 10., and 4.8 for thattheCRaccelerationefficiencydependsonlyweakly
2 0
≈
models S1-S4, respectively. The postshock CR pres- on the injection rate.
sure relative to the ram pressure of the Sedov solution Fig. 5 shows the volume integrated CR spectrum,
is Pc,2/ρ0US2T ≈ 0.35, 0.39, 0.42, and 0.25 for models G(p) = 4πg(p)r2dr in code units and the slope,
S1-S4, respectively. q = d(lnRGp)/dlnp 2. The maximum momentum
− −
We note that model S4 with B = 5µG has the depends on the SNR parameters as p u2t/κ
higher CR pressure than the fiducial model S1 with max ∝ s n ∝
8 KANG
Fig. 5.—IntegratedCRnumber,G = g(p)r2dr,inarbitraryunitsandtheslope,q=−d(lnG )/dlnp−2,areshownat
p p
t/to = 1, 3, 6, and 10. SeeTable 1 for moRdel parameters for models S1-4.
B E ρ−1/3. The values of log(p ) 6.5, 6.8, 5.7, for all models. The total CR energy accelerated up to
µ o ISM max ≈
and7.0att/t =10for modelsS1-S4,respectively,are t/t =10isE /E =0.55,0.57,0.62and0.39formod-
o o C o
roughly consistent with this relation. els S1-S4, respectively. So models S2 (larger E ) and
o
The CR spectra for models S1-S3 with high initial S3 (smaller Bµ) are a bit more efficient in transferring
Mach numbers show the canonical nonlinear concave the SN explosion energy to the CR energy, compared
curvature for p 1, which comes from the large den- to the fiducial model. As mentioned earlier, compari-
sity compression≫factor (ρ /ρ 4) due to the domi- son between models S1 and S3 indicates that strongly
2 0
nantCRpressure. Forthesemo≫delstheslopeq isclose amplifiedmagneticfieldcarriesfasterAlfv´enwavesand
to 2 atlow energies,but flattens to 1.6 at highener- leads to less efficient CR acceleration,even though the
gies below the upper momentum cu∼toff. For model S4 particles are accelerated to higher energies. Also com-
with lower initial Mach number the density compres- parisonbetweenmodelsS3andS5showsthatthetotal
sion is only ρ /ρ 4.8, so the flattening of G(p) at CR energy increases only slightly, less than 3 %, even
2 0
high momenta is n≈ot so prominent. Instead, the de- thoughthe injectionrateincreasesbyafactorofabout
crease of the CR injection rate due to the weakening 15 with a largervalue of ǫB for model S5. In model S4
subshock reduces the particle numbers at low momen- with a hotter ISM, the shock is weaker and so the CR
tum, leading to actually flatter spectra there. injection and acceleration are less efficient.
Fig. 6 shows the integrated energies, E /E = Oursimulationsimplythatonaverageabout10−4
4π e r2dr, where e , e , and e are theidenosity 10−3 of the incoming particles are injected to the C−R
i th kin C
oftRhermal,kinetic andcosmicrayenergy,respectively. population at the shock front and up to 50 % of the
The kinetic energy reduces only slightly and is similar SN explosion energy can be transferred to CRs during
Cosmic Ray Acceleration at SNRs 9
Fig. 6.— Integrated thermal, kinetic and CR energies inside the simulation volume as a function of time. See Table 1 for
model parameters for models S1-S5. The long dashed line in the lower left panel shows the volume integrated CR energy
for model S5.
the ST stage, if the magnetic field direction is radial perpendicular shocks may not be significantly smaller,
(i.e., quasi-parallel field). This last result is consis- compared to parallel shocks.
tentwiththe calculationspreviouslydonebyBerezhko
and collaborators using a different numerical scheme V. SUMMARY
(e.g., Berezhko & Vo¨lk 1997). The CR injection rate,
however, probably depends strongly on the angle be- The evolution of cosmic ray modified shocks de-
tween the magnetic field and the shock normal direc- pends on complex interactions between the particles,
tion. Inamorerealisticmagneticfieldgeometry,where waves in the magnetic field, and underlying plasma
a uniform ISM field is swept by the spherical shock, flow. We have developed numerical tools that can em-
only 10-20 % of the shock surface has a quasi-parallel ulatesomeofthoseinteractionsandincorporatedthem
field geometry (V¨olk et al. 2003). In the shock surface intoastandardnumericalschemeforgasdynamicprob-
region where the field is perpendicular, the injection lems. Specifically, diffusive shock acceleration can be
rate is expected to be reduced and so the CR accel- followed by a kinetic approach (i.e., CRASH code) in
eration efficiency would be smaller. Thus the CR en- whichadiffusionconvectionequationfortheCRdistri-
ergy conversion factor averaged over the entire shock bution function is solved with an appropriate diffusion
surface could be significantly smaller than 50 %, per- coefficient(Kang&Jones1991;Kanget al. 2001). The
haps about 10 %. On the other hand, Giacolne (2005) injection of CR particles can be treated by a thermal
showedthattheprotonscanbeinjectedefficientlyeven leakage injection model with a transparency function
at perpendicular shocks in fully turbulent fields due to τesc(ǫB,υ) (Kang et al. 2002). Drifts of resonantly-
fieldlinemeandering. Insuchcasetheinjectionrateat scattering Alfv´en waves and heating of the thermal
10 KANG
plasma due to the dissipation of those waves in the Berezhko, E. G., Ksenofontov, L. T., & Yelshin, V., 1995,
precursor region has been included through a simple EfficiencyofCRaccelerationinsupernovaremnants,Nu-
model (Jones 1993;Kang & Jones 2006). clear Physic B., 39, 171.
Inthepresentpaper,weappliedthesphericalCRASH Berezhko, E. G., & V¨olk, H. J., 1997, Kinetic theory of
codetothe problemoftheCRaccelerationattherem- cosmic rays and gamma rays in supernova remnants. I.
nant shock from typical Type Ia SNe propagating into Uniform interstellar medium,Astropart. Phys. 7, 183
auniforminterstellarmedium. Themainresultsofour Berezhko, E. G., Ksenofontov, L. T., & V¨olk, H. J., 2003,
simulations can be summarized as follows: Confirmation of strong magnetic field amplification and
1) Since the CR injection and acceleration depend nuclearcosmic ray acceleration in SN1006, A&Ap,423,
L11
primarilyon the shockMachnumber, the temperature
of the ISM is an important factor. SNRs in a warm Berezhko,E.G.,&V¨olk,H.J.,2006, Theoryofcosmicray
ISM inject about 10−4 10−3 of the particles passed production in the supernova remnant RX J1713.7-3946,
through the shock and−transfer about 50% of the ex- A&Ap,451, 981
plosion energy to the CR component. SNRs in a hot Blandford, R. D., & Ostriker, J. P., 1978, Particle acceler-
ISMinject10timeslessparticles,buttheystilltransfer ation by astrophysical shocks, ApJ,221, L29
about 40 % of the explosion energy to CRs. Blandford, R.D.,&Eichler, D.,1987, Particleacceleration
2) At sources, the hadronic CR spectrum should at astrophysical shocks - a theory of cosmic-ray origin,
be consistent with a power-law function, N(E)dE Phys. Rept.,154, 1
E−2dE up to p/mpc 102, but it may gradually fla∝t- Drury, L. O’C., 1983, An Introduction to the Theory of
ten to E−1.6 toward ∼the Knee energy and beyond up Shock Acceleration of Energetic Particles in Tenuous
to E 1016Z eV (p/m c 107). Plasmas, Rept. Prog. Phys., 46, 973
p
∼ ∼
3) About 40-50 % of the explosion energy can be Drury, L. O’C., Ellison, D. E., Aharonian, F. A. et al. ,
transferedtoCRs,ifthemagneticfieldisradial(i.e.,quasi- 2001,Testofgalacticcosmic-raysourcemodels-Working
parallel shocks). However,more realistic consideration Group Report, Space ScienceReviews, 99, 329
ofthefieldgeometryrelativetothesphericalshocksur- Giacalone, J., 2005, The Efficient Acceleration of Thermal
facecouldleadtomuchsmallerenergyconversionrate. Protons by Perpendicular ShocksApJ, 628, L37
So a conservative estimate could be order of 10 %.
Jones, T. W., 1993, Alfv´en wave transport effects in the
4) If the magnetic field is amplified in the precursor time evolution of parallel cosmic-ray-modified shocks
regiondue to the streaminginstability, as indicated by ApJ,413, 619
recentX-rayobservationsofyoungSNRs, the particles Kang, H., & Jones, T. W., 1991, Numerical studies of dif-
are accelerated to higher energies due to smaller dif- fusive particle acceleration in supernova remnants, MN-
fusion coefficient. However, faster Alfv´en waves in the RAS,249, 439
precursortend to reduce the CR injection and acceler-
Kang, H., Jones, T. W., & Ryu, D., 1992, Acoustic insta-
ation efficiencies.
bility in cosmic ray mediated shocks, ApJ, 385, 193
In conclusion, the galactic cosmic rays possibly up
Kang, H., Jones, T. W., LeVeque, R. J., & Shyue, K. M.,
to 1016Z eV could be originated from supernova rem-
2001, Time Evolution of Cosmic-Ray Modified Plane
nants. This result is quite robust, regardless of SNR Shocks,ApJ, 550, 737
model parameters and details of microphysics involved
Kang,H.,Jones,T.W.,&Gieseler,U.D.J.,2002,Numer-
in the injection process,providedthat the Bohmdiffu-
ical Studies of Cosmic-Ray Injection and Acceleration,
sion is valid at young SNRs of several thousand years
ApJ,579, 337
old.
Kang, H., & Jones, T. W., 2002, Acceleration of Cosmic
RaysatLargeScaleCosmicShocksintheUniverse,Jour-
ACKNOWLEDGEMENTS
nal of Korean Astronomical Society, 35, 159
This work was supported for two years by Pusan Kang, H., 2003, Cosmic Ray Acceleration at Cosmological
NationalUniversityResearchGrant. Theauthorwould Shocks: Numerical Simulations of CR Modified Plane-
like to thank T.W. Jones for helpful comments on the ParallelShocksJournalofKoreanAstronomicalSociety,
paper. 36, 1
Kang, H., & Jones, T. W., 2006, Numerical studies of dif-
REFERENCES
fusive shock acceleration at spherical shocks Astropart.
Phys. 25, 246
Aharonianetal.,H.E.S.S.collaboration,2004,High-energy
Koyama, K., Petre, R., Gotthelf, E. V. et al. , 1995, Ev-
particleaccelerationintheshellofasupernovaremnant,
idence for Shock Acceleration of High-Energy Electrons
Nature, 432, 75
in the SupernovaRemnant SN:1006, Nature, 378, 255
Bamba, A., Yamazaki, R, Ueno, M. & Koyama, K., 2003,
Lagage, P. O., & Cesarsky, C. J., 1983, The maximum
Small-Scale Structureof the SN 1006 Shock with Chan-
energy of cosmic rays accelerated by supernova shocks,
dra Observations, ApJ, 589, 827
A&Ap,118, 223
Bell, A. R., 1978, The acceleration of cosmic rays in shock
fronts. I, MNRAS,182, 147