Cosmic ray acceleration and escape from supernova remnants

0.53 MB·English

by  AR Bell

Mon.Not.R.Astron.Soc.000,1–17(2013) Printed31January2013 (MNLATEXstylefilev2.2)
Cosmic ray acceleration and escape from supernova
remnants
⋆
A. R. Bell , K.M. Schure, B. Reville and G. Giacinti
3
Clarendon Laboratory, Universityof Oxford, Parks Road, Oxford OX1 3PU, UK
1
0
2

n
a
J
ABSTRACT
0
3 Galacticcosmicray(CR)accelerationtothekneeinthespectrumatafewPeVisonly
possible if the magnetic field ahead of a supernova remnant (SNR) shock is strongly

] amplifiedbyCRescapingtheSNR.Amodelformulatedintermsoftheelectriccharge
E
carriedbyescapingCRpredictsthe maximumCRenergyandtheenergyspectrumof
H CR released into the surrounding medium. We find that historical SNR such as Cas
. A, Tycho and Kepler may be expanding too slowly to accelerate CR to the knee at
h
the present time.
p
- Key words: cosmicrays,accelerationofparticles,shockwaves,magneticfield,ISM:
o
supernova remnants
r
t
s
a
[

1 1 INTRODUCTION forBohm diffusion(Stageetal 2006, Uchiyamaet al2007).
v Furthermoreif themean free pathwere much larger than a
Duringdiffusiveshockacceleration(DSA)cosmicrays(CR)
4 LarmorradiusaccelerationbySNRtothekneeintheGalac-
6 gain energy by repeatedly passing back and forth between ticCRspectrumwouldbeverydifficult.ThemaximumCR
2 theupstreamanddownstreamplasmas(Krymskii1977,Ax-
energyisdeterminedbytheconditionthatt cannotex-
7 ford et al 1977, Bell 1978, Blandford & Ostriker 1978). CR accel
ceedtheaget ofanSNR(Lagage&Cesarsky1983a,b).
. diffuseaheadoftheshocktoformaprecursorwithanexpo- SNR
1 Assuming D = cr , the maximum CR energy T in eV
nentialscaleheightD /u whereu istheshockvelocityand u g max
0 u s s is given by T = 0.03B u2t PeV where B is the
3 DuistheCRdiffusioncoefficientupstreamoftheshock.The upstream magmnaextic field inµmGic7ro1G00a0uss, u is theµshGock ve-
1 averagedwell-timespentupstreamoftheshockbetweensuc- locity in units of 10,000 km s−1 and t 7 is the SNR age
v: cessive shock crossings is 4Du/cus for relativistic particles in 1000’s of years. For a typical young S10N00R expandinginto
(Bell2012).This,alongwiththecorrespondingdownstream
i theinterstellarmediumwithoutmagneticfieldamplification
X dwell-time,determinestherateatwhichCRareaccelerated.
B =3, u =0.5 and t =0.4, giving T =0.01PeV
Lagage&Cesarsky(1983a,b)showedthatthetimetakenfor µG 7 1000 max
r whichfallsafactorof 100shortofthatrequiredtoexplain
a CRaccelerationistaccel =4Du/u2s+4Dd/u2dwhereDdisthe the Galactic CR spect∼rum. SNR in the Sedov phase do not
downstream diffusion coefficient and u is the downstream
d fare better.Their shock velocities decrease in proportion to
fluid velocity in the shock rest frame. Since u = u /4 for
d s time−3/5 andradius−3/2,solittlebenefitaccruesfromtheir
a strong shock it might appear that thedownstream dwell-
larger age and radius. This posed a serious problem for dif-
time determines the acceleration rate, but we can expect
fusive shock acceleration as an explanation of the Galactic
D D partly because the magnetic field is increased by
d ≪ u spectrumuntilitwasshownthataplasmainstabilitydriven
compression by the shock, partly because the compressed
by streaming CR in the upstream precursor could amplify
downstream magnetic field is more closely perpendicularto
themagneticfieldaheadoftheshockandfacilitaterapidac-
the shock normal, and partly because thedownstream field
celeration tohigherenergies (Lucek &Bell 2000, Bell 2004,
is disturbed and more irregular after passing through the
2005).
shock.If thedownstream and upstream dwell-times are the
samet =8D /u2.Afurthercommonassumptionisthat
accel u s
The phenomenon of magnetic field amplification pro-
Bohm diffusion applies: D = cr or D = cr /3 where r
u g u g g
vides a mechanism by which the CR energy can be raised
is the CR Larmor radius. This assumes a diffusion model
significantly beyond 0.01PeV, but there remains the ques-
in which CR are scattered by irregularities in the magnetic
tion of why the fields are amplified to the observed magni-
fieldsuchthatthescatteringmeanfreepathisoftheorderof
tude,upto100’sµGinthehistoricalSNR(Vink&Laming
theCRLarmorradius.Thereissomeobservationalevidence
2003,Berezhkoetal2003,V¨olketal2005),andwhyCRare
accelerated toafewPeVratherthan0.1or10PeV.Wetry
⋆ E-mail:[email protected] to answer these questions by examining the self-consistent
(cid:13)c 2013RAS
2 A.R. Bell, K.M. Schure, B. Reville and G. Giacinti
interactionbetweenstreamingCR,behavingkinetically,and Falle1986,Drury&Downes2012,Malkov&Diamond2009,
theupstream plasma behavingmagnetohydrodynamically. Rogachevskii et al 2012, Schure & Bell 2011) tend to have
IthasbeenknownformanyyearsthattheescapeofCR longer growth times but turbulent amplification may in-
upstreamoftheshockisanimportantpartoftheoverallac- creasethegrowthrate(Bykovetal2011).Anextendeddis-
celerationprocessasdiscussedbelowinthefinalparagraphs cussion of instabilities driven by cosmic ray streaming can
of section 3. We find that the combined CR-MHD system befound in Schureet al (2012).
organises itself to allow a suitable number of CR to escape From the Lagage & Cesarsky (1983a,b) limit as de-
upstream.The CR drivemagnetic field amplification which scribed in section 1 it is clear that proton acceleration to
in turn regulates the number of escaping CR. If a smaller a few PeV is only possible if the magnetic field is strongly
numberof CR escaped, the magnetic field would beinsuffi- amplified above its characteristic interstellar value of a few
cientlyamplifiedtoconfineandacceleratetheCR.Ifalarger µG. The following argument places an upper limit on the
number of CR escaped, the magnetic field would grow too maximumenergyofaCRthatcanbestronglyscatteredby
rapidlytoallowtheirescape.Henceaself-regulatingsystem amagneticfieldgrowingonthescaleofaCRLarmorradius.
isset upthatdeterminesthenumberand maximumenergy Thebackgroundthermalplasmaishighlymagnetisedonthe
of escaping CR. scale of a PeV CR Larmor radius. Consequently the mag-
The paper is organised as follows. Sections 2 and 3 neticfieldis‘frozenin’tothethermalplasma.Magneticfield
present approximate calculations showing how a limit on amplificationoccursastheplasmamovesandstretchesfield
the CR energy is placed by the need for CR to drive mag- lines. Assuming the perturbed magnetic field does not far
netic field amplification by escaping upstream. Sections 4 exceedtheinitialfield,thej Bforcedisplacesaplasma
CR
to7describeVlasov-Fokker-Planck(VFP)simulations that elementamaximumdistances× (j Bt2)/2ρintimet
max CR
∼
support theargumentsof sections 2 and 3. Sections 8to 11 whereB is theinitial seed field.For thestretched magnetic
apply the results to supernova remnants and the Galactic field to strongly scatter CR as required for Bohm diffusion
cosmic rayspectrum.Readersunfamiliar with VFPsimula- the displacement must grow to the order of a CR Larmor
tions may wish to read sections 1 to 3 and 8 to 11 before radius, that is s T/cB where T is the CR energy in
max
∼
returningtothecomputationalvalidationandillustrationin eV.TheCRcurrentdensityj carriesanenergyfluxj T
CR CR
sections 4 to 7. in theupstream plasma rest frame which cannot far exceed
the energy flux ρu3/2 carried by the upstream plasma into
s
the shock, where ρ is the upstream mass density. We de-
fine a CR acceleration efficiency η such that j T = ηρu3.
2 CONDITIONS FOR STRONG MAGNETIC CR s
We then have two equations: s =(j Bt2)/2ρ T/cB
FIELD AMPLIFICATION max CR ∼
and j = ηρu3/T. When combined, they yield a CR en-
CR s
Wreseonaasnsutmhyebtrhidat(NmRaHgn)eitnisctafibeilldityisdegsecnreibraetdedbybByeltlh(e20n0o4n)-. teorgTy T∼∼1.5((ηηcuu373s))11//22BBµtG.tT10h0i0sPeexVprweshseiorenufo7risTthiseesqhuoicvkalevnet-
This is one of a class of plasma instabilities driven by CR locity in units of 10,000km s−1, BµG is the seed magnetic
streaming. In its simplest form, CR have a Larmor radius fieldinµGfromwhichamplificationbegins,andt1000 isthe
much greater than the wavelength of spiral perturbations time in 1000’s of years. Characteristically for young SNR,
in a zeroth order uniform magnetic field. Because of their u7 = 0.5, t1000 = 0.4, η = 0.03 (see Appendix A) and
large Larmor radius the streaming CR, carrying an elec- BµG = 3 where the magnetic field of a few µG represents
tric current density j , are essentially undeflected by the the seed field from which the instability grows. With these
CR
perturbedfield butthej Bforce actstowards thecen- estimates, T 0.1PeV, which is an order of magnitudeless
CR× ∼
tre of thespiral. A corresponding reactive force acts on the than the energy of the knee. This estimate is independent
backgroundplasmatoexpandthespiral.Thisstretchesand of anydetailed instability theoryexcept for theassumption
increases the magnitude of the perturbed magnetic field, thatmagneticfieldamplificationtakesplacethroughplasma
thereby increasing the jCR × B force in a positive feed- motionsgenerated bythejCR×B force actingon thescale
back loop that drives the instability. NRH appears to be ofaCR Larmorradius.Ithighlightsthedifficultyofampli-
the most rapidly growing instability driven by CR stream- fying magnetic field by orders of magnitude on the scale of
ing on a relevant scalelength. Other instabilities that have the Larmor radius of a PeV proton and the need for an in-
attracted significant attention are the resonant Alfven in- stabilitythatcanprovidenon-lineargrowthofthemagnetic
stability(Kulsrud&Pearce1969,Wentzel1974) thatgrows field.
withspatialwavelengthsspatiallyresonantwiththeCRLar- The NRH instability offers a way round this diffi-
mor radius and the Weibel instability (Weibel 1959) that culty by initially growing very rapidly on a scale much less
grows quickly on the spatial scale of an electron or pro- than the CR Larmor radius. Since the small-scale mag-
ton collisionless skin depth, c/ω = 5.3(n /cm−3)−1/2km, netic field grows unstably the j B force grows expo-
pe e CR
×
c/ω = m /m c/ω .TheWeibelinstabilityisimportant nentially in time. The scale-size increases to the CR Lar-
pi p e pe
for interactions engaging thermal or mildly suprathermal mor radius during non-linear growth. The NRH growth
p
electronsandions,andmaybeeffectiveforCRacceleration rate is γ = (kj B /ρ)1/2 where B is the zeroth or-
CR 0 0
to low energies, but it grows on a scale too small to scat- der magnetic field and k is the wavenumber on which
ter PeV ions which have a Larmor radius of 109c/ω . the instability grows. The maximum NRH growth rate is
pe
∼
The Alfven instability grows on the desired spatial scale γ = 0.5j (µ /ρ)1/2 which occurs at a wavenumber
max CR 0
butgrowslessquicklythantheNRHinstabilityintheSNR k =0.5µ j /B .γ is(k r )1/2 timeslargerthan
max 0 CR 0 max max g
shock environment. Instabilities that grow on scales larger the NRH growth rate for kr = 1, thus allowing the mag-
g
than the CR Larmor radius (Bykov et al 2011, Drury & netic field and the j B force to increase rapidly. The
CR
×
(cid:13)c 2013RAS,MNRAS000,1–17
Cosmic ray acceleration and escape from supernova remnants 3
NRH and Alfven growth rates are similar at kr = 1 (see is carried by the highest energy CR being accelerated. The
g
AppendixEandBell 2004).Usingtheabovenomenclature, Larmor radius of a PeV proton in the interstellar medium
k r =ηu3/cv2 = 2 104η u3n B−2 where v is the is comparable with the radius of the historical SNR and
max g s A × 0.03 7 e µG A
Alfven speed, η = 0.03η and n is the electron density its scattering mean free path is very much greater than the
0.03 e
in cm−3. For a discussion of the value of η see Appendix SNR radius. These escaping CR pass into the interstellar
A and Bell (2004). The growth time of the fastest growing medium with low probability of further encounter with the
mode is then γm−1ax = 50η0.03−1n−e1/2u−73TPeV years where SNRshock.StronglyscatteredCRatlowerenergiesarecon-
T istheenergyinPeVoftheCRdrivingtheinstability. finedbytheshockandareadvectedawaydownstreamafter
PeV
γ−1 = 400 years for η = 1, n = 1, u = 0.5 which acceleration.
max 0.03 e 7
implies that even the NRH instability struggles to amplify Suppose that the CR distribution follows a p−4 power
the magnetic field sufficiently to accelerate CR to PeV en- law upto amomentum pmax, thenin steady statetheelec-
ergies in the historical SNR. For growth by many e-folding trical current of CR escaping in a small band of energies
theNRHinstabilitymustbedrivenbyCRwithenergiesless abovetheenergy eTmax =cpmax is
than 1PeV. j =eu πp3 f (p ) (2)
Asremarkedabove,thefastest growing modegrows on CR s max 0 max
ascalekm−1ax whereasefficientCRscatteringrequiresfluctu- wheref0 istheisotropic part of theCR distribution inmo-
ationsinthemagneticfieldonascaler .Theaboveanalysis mentum space at the shock. Equation (2) is derived from
g
showed that k r = 2 104η u3n B−2. The instabil- equation(11a)belowbyintegratinginspaceacrosstheshock
max g × 0.03 7 e µG
ity initially grows on a scale too small to effectively scatter and deriving the rate at which CR reach the momentum
PeVCR.Twofactorssavethesituation.Firstly,asthemag- pmax. The derivation can befound in AppendixB. In com-
netic field grows from a few µG to a few hundred µG the parison the CR pressure at the shock, where the CR distri-
CR Larmor radius decreases by thecorresponding factor of bution extends down to p = mc and m is the proton mass
100. Secondly,thecharacteristic scalelength of the struc- is
∼
tured magnetic field increases during non-linear growth of 4π p
P = cp4 f (p )ln max (3)
the instability. Figure 1 is drawn from figures 3 and 4 of CR 3 max 0 max mc
Bell (2004). The graph of energy densities shows how the (cid:16) (cid:17)
giving
magnitude of the magnetic field grows to many times its
sleenedgthvagluroewisndaurtiinmgeth∼is1t0iγmm−e1axa.sTshhoewcnhabryacttheerisgtriecys-sccaalele- jCR =0.05uTsmPaCxR (4)
images.Bythetimet=10thescalelengthhasgrowntothe where we have assumed ln(p /mc) = 14 for a CR dis-
max
sizeofthecomputationalboxfromitsinitiallymuchsmaller tribution extending from 1GeV to 1PeV. The condition for
scale. Both k and rg decrease by a large factor by the time magnetic field amplification that jCRdt = 10 ρ/µ0 can
C10Rγm−c1auxr.reTnhtesosiimtuwlaatsioimnproepssriobdleucteoddienmfiognusrtera1tehsatdroangfixCeRd berearranged to give a valuefor TRmax: p
sficgautrteerpinrogvibdyesthsteromnaggenveitdiecnficeelfdorarftaepridti=nit1ia0lγgm−r1aoxw,tbhuotftthhee Tmax =0.005PρCuR2s ρu3strµρ0 (5)
fastest growing modes into large scale magnetic structures which is equivalent to
able to scatter CR on thescale of a Larmor radius.
P
T =8 n1/2u3t CRPeV (6)
max e 7 1000 ρu2
s
For characteristic values (u = 0.5, n = 1, t = 0.4,
3 CR ESCAPE UPSTREAM OF THE SHOCK 7 e
P /ρu2 = 0.3) T 100TeV. As expected from the
CR s max ∼
The previous section does not provide an estimate of the above discussion, this falls short of the few PeV required
magnitude of the amplified magnetic field, but it does de- to explain Galactic CR. This will be discussed further in
fineconditionsunderwhichstrongmagneticfieldamplifica- section 8.
tionoccurs.CRaccelerationtoPeVenergiesrequiresstrong Zirakashvilietal(2008a,b)developedarelatedanalytic
magneticfieldamplificationwhichfixesthenumberofinsta- modeloftheexcitationoftheNRHinstabilityandCRcon-
bilitye-foldingstotherange5-10.ThisinturnfixestheCR finement upstream of the shock. They derived an estimate
current.Fromfigure1wetaketheconditionforstrongmag- for the maximum CR energy with a similar dependency on
netic field amplification in a particular volume of upstream n1/2u3t in the numerator of their equation (19) but with
e
plasmatobethat γmaxdt 5inthetimebeforetheshock an additional denominator that depends on the magnitude
∼
overtakes it. Since γ = 0.5j µ /ρ, the condition for oftheamplifiedmagneticfield.Theiranalysisoperatesat a
R max CR 0
field amplification is moredetailedlevelthanourssincetheyconsiderthegrowth
p
in amplitude of the amplified magnetic field over a range
Q = j dt=10 ρ/µ (1) of scales and small angle CR scattering by small scale field.
CR CR 0
Z However, the dominant physics is similar in their analysis
p
whereQ isthetotalelectricchargeofCRpassingthrough and ours, and their estimate of the CR energy (equation
CR
unit surface area upstream of the shock before the shock (19)) is similar to that in ourequation (6).
arrives. (It makes no difference whether the CR passing Notethatthemagneticfielddoesnotenterintothees-
through the plasma have high or low energies. It is only timation of the maximum energy CR in equation (6). The
the number of CR escaping upstream that counts.) Since magnetic field is assumed to grow to whatever magnitude
only the highest energy CR escape, the areal charge Q and spatial scale required to confine CR at energies less
CR
(cid:13)c 2013RAS,MNRAS000,1–17
4 A.R. Bell, K.M. Schure, B. Reville and G. Giacinti
Figure 1. Plots of the magnitude of a magnetic field driven by CR streaming into the page (2D slices of a 3D simulation). The field
growsfromnoiseonasmallspatialscaleatt=0.Byt=10thespatialscalehasgrowntothesizeofthecomputationalbox.Thegraph
showshowthemeanmagnetic, kineticandthermalenergydensitiesincreasewithtime.Thisfigureiscomposedfromfigures3and4of
Bell(2004)wherefurtherdetailscanbefound.Theunitsoftimearesuchthatγmax=1.26.
than T . This assumption is justified on the basis that theshockstructureisincludedsincetherelativisticCRpres-
max
themagnitudeofthefieldincreasesbyasubstantialnumer- sure increases the density jump at the shock resulting in
ical factor in time γ−1 even in the non-linear regime and a spectrum flatter than T−2 at high CR energy. However,
max
that this is accompanied by rapid growth in characteristic the conditions for a steady state solution do not present a
spatialscale.Forexample,themagneticfieldanditsspatial compellingargumentforafreeescapeboundarysincetime-
scale are very much larger in figure 1 when γ dt 8 dependent solutions offer an acceptable option. Indeed the
max
≈
than when γ dt 5 even though the time is different Lagage & Cesarsky limit on CR energy is based on the as-
max ≈ R
by a factor of only 1.6. Allowing the field to grow for just sumptionthatthemaximumCRenergyincreaseswithtime.
R
a little longer gives a much larger magnetic field and spa- ThecaseforCRescapewasstrengthenedbytherecognition
tial scale. Consequently the primary parameter is not the thatmagneticfieldamplificationisessentialsincetherewill
magnitude of the magnetic field but the charge carried by always bea distance upstream at which amplification is in-
escaping CR. operativeandCRmustescape.Similarly,theremustalways
Itisclearthatmagneticfieldamplificationisonlypossi- beanupperlimittothespatialscaletowhichfieldisampli-
bleifapopulationofhighenergyCRescapeupstreamofthe fied and a corresponding limit, through the Larmor radius,
shock. A magnetic field capable of stopping their escape is on the momentum to which CR can be accelerated. This
onlyproducedafterthearealchargeQ hasescaped.Thus ledtotheimposition of afreeescapeboundaryeitheratan
CR
CR escape upstream of the shock is an essential aspect of imposed distance upstream of the shock (eg Caprioli et al
diffusive shock acceleration when magnetic field amplifica- 2010b,Ohiraetal2010)oratanimposedCRmomentum(eg
tionisoperational. IfachargeQ hasnotpassedthrough Ellison & Bykov 2011). Whether the free escape boundary
CR
aparticularpointupstreamoftheshockthenCRcannotbe is imposed in momentum or configuration space, the mo-
confinedat that distance from theshock. mentumatwhichCRescapedependsonthemagneticfield.
Previous authors have also concluded that CR escape Someauthors(egPtuskinetal2010)assumethatthefieldis
upstream is inevitable, but generally for different reasons. amplified until themagnetic energy density reaches a given
One possibility is that CR escape through filamentary cav- fraction of ρu2s. The fractional magnetic energy density can
ities in the magnetic field (Reville & Bell 2012, 2013) but bechosentomatchobservation(Ptuskinetal2010,V¨olket
usuallyitisassumedthatCRescapeatafreeescapebound- al2005). Alternativelythefield canbechosentomatch the
ary in position or momentum. Free escape boundaries were saturation value estimated by Bell (2004). Others (Caprioli
introduced at an early stage in the development of shock etal2010a,Drury2011)chooseamathematicalformforthe
acceleration theory (Ellison et al 1981). It was appreciated amplifiedmagneticfieldthatallowsformultiplepossibilites.
thatsteadystatemodelswereimpossiblewithoutCRescape Ptsukin & Zirakashvili (2003) and Caprioli et al (2009a,b,
(Ellison & Eichler, 1984, Berezhko & Ellison 1999, Malkov 2010b)includemagneticfieldgenerationduetotheresonant
etal2000).ThestandardT−2 testparticleenergyspectrum instability as it develops ahead of the shock. Vladimirov et
atstrongshocksdivergesifintegratedtoinfiniteenergy.The al (2006) include magnetic field generation in response to
divergence is even stronger if non-linear CR feedback onto gradientsintheCRpressure.Ellison et al(2012) choosean
(cid:13)c 2013RAS,MNRAS000,1–17
Cosmic ray acceleration and escape from supernova remnants 5
amplified magnetic field suitable for SNR expansion into a of theCR to theMHD timescale to a minimum while stay-
circumstellar wind. ingclosetotherangeofconceivableSNRexpansionspeeds.
In contrast to the above, CR escape in our model is Our greatest approximations are made in the momentum
determined by the escaping CR electric charge rather than space representation of the CR distribution function since
the magnetic field generated in the upstream plasma, al- this is the aspect of the calculation in which thenumberof
thoughmagneticfieldamplificationisimplicitlyrequiredby dimensions can be reduced.
ourmodel.Revilleetal(2009)alsousedtheescapingfluxto The Vlasov-Fokker-Planck (VFP) equation (equation
calculategrowthratesinanefforttomotivatearealisticfree- 8) is important in the physics of laser-produced plasmas
escape boundary location in their steady-state non-linear whereitissolved infinitedifferenceform tomodelelectron
model. Their discussion on self-regulation of CR precursors transport. The successful use of VFP simulation to model
howeverdid not extendto themaximum CR energy. electron transport in laser-produced plasmas stretchesback
more than 30 years (Bell et al, 1981) so it is natural to ap-
plythetechniquestoCRwhichobeythesameequation.The
distribution function f(r,p,t) of charged particles (cosmic
4 A NUMERICAL MODEL
rays or energetic electrons) is usually represented in spher-
We now set out to test the above conclusions as far as ical co-ordinates (p,θ,φ) in momentum space. A common
we are able with a numerical model that includes the representation of the distribution function is as a sum of
self-consistent interaction of CR modelled kinetically with spherical harmonics:
a background plasma modelled magnetohydrodynamically.
Standard MHD equations describe the background plasma f(r,p,t)= flm(r,p,t)Pl|m|(cosθ)eimφ
except that a −jCR ×B force is added to the momentum Xl,m
equation: l=0, m= l,l f−m =(fm)∗ (9)
∞ − l l
du 1
ρdt =−∇P − µ0B×(∇×B)−jCR×B (7) whahremreonfilcm.(fr,mp(,rt,)p,ist)tihsefucnocetffiiocnieonfttfiomre,thpeos(itl,imon)asnpdhemriacga-l
l
asdescribedinLucek&Bell(2000)andBell(2004).TheCR nitudeofmomentump.Thesphericalharmonicsdecribethe
distributionfunctionf(r,p,t)atposition randmomentum angularstructureonshellsofconstantmagnitudeofmomen-
pisdefinedinthelocalfluidrestframeandevolvesaccord- tum.Reviewsof theVFPtechnique,orpaperscontaininga
ing totheVlasov-Fokker-Planck(VFP) equation significant review element, are Bell et al (2006), Tzoufras
et al (2011) and Thomas et al (2012). For information on
df ∂f ∂u ∂f ∂f
dt =−vi∂ri +pi∂rij ∂pj −ǫijkeviBj∂pk +C(f) (8) trehaedaeprpislicreafteirornedoftoVRFePvitlleech&niBqeulels(2t0o13C)Rwahciccheleinracltuiodnestahne
where C(f) is an optional collision term included to repre- appendix setting out the full VFP equations for a spheri-
sent scattering by magnetic fluctuations on a small scale. cal harmonic expansion. Bell, Schure & Reville (2012) also
The electric field is zero in the local fluid rest frame. applyVFPtechniquestoCRacceleration.Theuseofanex-
Quadratic terms in the local fluid velocity u are neglected pansionintensors(usedhere)asanalternativetospherical
on the assumption that u c. The CR current j , re- harmonics is discussed by Schure& Bell (2012).
CR
≪
quired for insertion into the MHD momentum equation, is VFPsimulationwasusedbyBelletal(2011)forthecal-
calculated by integration over f in momentum space. The culation of CR acceleration by oblique shocks. They found
magnetic field in theCR kineticequation is taken from the that an expansion to 15th harmonic could be needed for
MHD calculation. oblique shocks because of the abrupt change in magnetic
As in Lucek & Bell (2000) and Bell (2004) three spa- field direction at the shock, but that only a few harmon-
tial dimensions are needed to represent the turbulence ad- ics are needed for quasi-parallel shocks. Here we reduce
equately. Since our aim is to investigate the mutual inter- thecomputationalsizeoftheproblembymodellingparallel
actionsofmagneticfieldamplification,CRacceleration and shocks in which the zeroth order magnetic field is parallel
CR escape upstream of the shock we need to model the to the shock normal. There are good reasons to suppose
complete system including the shock and the complete CR that the first few terms in the expansion capture the es-
precursor. Because we model the detailed interaction be- sential physics. Firstly, as shown by Bell (2004) the NRH
tween the CR and the distorted magnetic field we need to instability is driven by the CR current density j . Higher
CR
resolvetheCRLarmorradiusinconfigurationspaceandthe order anisotropies do not directly contribute to the insta-
rotation of CR trajectories in momentum space. As a con- bility. Secondly, the CR precursor scaleheight is c/u times
s
sequence the numerical model should be 6-dimensional in the CR scattering mean free path in standard DSA theory.
momentum-configuration space with spatial scales extend- DSAtheoryisbasedonthediffusiveapproximationinwhich
ing from the CR Larmor radius of the lowest energy CR only the first order anisotropy is needed and only the first
to theprecursor scaleheight of thehighest energy CR. This two terms (isotropy and drift anisotropy) in the harmonic
would be impossible without extraordinary computational expansion are retained. In diffusion theory the higher or-
resourcessoourstrategyistodesignacomputationalmodel der anisotropies are damped by scattering. Hence it might
thatincludesalltheimportantprocessesataminimallevel. be supposed that only the first two terms are needed and
We retain the three dimensions in configuration space but an adequate representation of the CR distribution func-
limittherangeofCRmomentumtoafactorof10sothatwe tion might be f(r,p,t) = f0(r,p,t) + f0(r,p,t)cosθ +
0 1
donothavetoresolvetheLarmorradiusofverylowenergy f1(r,p,t)sinθeiφ .This‘f +f ’expansionallowsfreeCR
ℜ{ 1 } 0 1
CR. We choose a shock velocity u = c/5 to keep the ratio propagation along magnetic field lines but restricts trans-
s
(cid:13)c 2013RAS,MNRAS000,1–17
6 A.R. Bell, K.M. Schure, B. Reville and G. Giacinti
port across the magnetic field because the direction of the The notation is simplified by introducing the vector g
i
anisotropic drift term is rotated by the field. However, this relatedtoshearandvorticity( f )inCRmotiontorep-
1
∇×
expansionomitsanessentialfeatureofCRtransport.Itdoes resenttheoff-diagonalcomponentsofthestresstensor.The
notallowCRtogyrateastheytravelalongamagneticfield on-diagonal components of f (i = j) are approximately
ij
line. The f +f expansion allows CR to propagate along accounted for by multiplying the c∂f /∂r term in equa-
0 1 0 i
−
fieldlinesandseparatelyitallowsCRtogyratearoundfield tion(11b) by9/5 toallow forthestresstensorcontribution
lines,butitdoesnotallowCRtodobothatthesametime. to compressional waves and to allow freely streaming CR
Spiral trajectories require the inclusion of the off-diagonal to propagate at 3/5c instead of 1/3c, as shown in Ap-
components of the stress tensor f . Without the stress ten- pendix D. The equations to be solved (their derivation is
2 p p
sor, CR cannot resonantly interact in space with magnetic given in AppendixC) are then
perturbations on the scale of a Larmor radius. Clearly this
∂f ∂(f u ) c∂f ∂u 1 ∂(p3f )
would rule out the Alfven instability, and more surprising 0 + 0 i = i + i 0
it also rules out the NRH instability. The role of the stress ∂t ∂ri −3∂ri ∂ri 3p2 ∂p
tensorisdiscussedbySchureetal(2011)inwhichthelinear
∂f ∂(f u ) 9 ∂f 1 ∂g ceB
NRH dispersion relation is derived with the perturbed CR i + i j = c 0 cǫ k ǫ jf
distribution expressed as a tensor expansion. We proceed ∂t ∂rj −5 ∂ri − 5 ijk∂rj − ijk p k
on the basis that theCR distribution function may be ade-
quatelyrepresentedbyanisotropicpart(0thorderinus/c) ∂gi + ∂(giuj) =cǫ ∂fk ν g (12)
plus a drift component (1st order in us/c) plus a term rep- ∂t ∂rj ijk∂rj − B i
resenting thestress tensor (2nd order in us/c). This second Equivalently,expressed in vectornotation,
orderexpansionismoreeasily representedintheequivalent
tensor notation instead of spherical harmonics: ∂∂ft0 +∇.(uf0)=−3c∇.f1+ ∇3p.u2 ∂(p∂3pf0)
p p p
f(r,p,t)=f (r,p,t)+f (r,p,t) i +f (r,p,t) i j (10)
0 i p ij p2
∂f 9c c
wherethetraceoffij iszerobecauseitisalreadyaccounted ∂t1 +∇.(uf1)=− 5 ∇f0− 5∇×g1−Ω×f1
for in the isotropic term f . Following Johnston (1960) the
0
reducedVlasovequationfortheCRdistributionfunctionis ∂g
1 + .(ug )=c f ν g (13)
then: ∂t ∇ 1 ∇× 1− B 1
∂f ∂(f u ) c∂f ∂u 1 ∂(p3f ) where Ω = ecB/p is the vector CR Larmor frequency and
0 + 0 i = i + i 0 (11a)
∂t ∂ri −3∂ri ∂ri 3p2 ∂p we take νB = ecB/p. The presence of the curls of f1 and
g facilitates thepropagation of transversemodesin f and
1 1
∂∂fti + ∂(∂friujj) =−c∂∂fr0i − 25c∂∂frijj −ǫijkcepBjfk (11b) gC1Rapsrnoepeadgeadtifoonrhatelaicnalamngolteiotnotahloenwgamveavgencettoirckfi.eldlinesor
Theterminationoftheharmonicexpansionatthestress
tensor makes the computation tractable with available re-
∂f ∂(f u ) c ∂f ∂f c ∂f
∂tij + ∂irjkk =−2(cid:18)∂rji + ∂rji(cid:19)+ 3δij∂rkk spoaunrscieosn.pArpovpiednedsiacensaDdeqanudateEreshporewsetnhtaattiothneotfrtuhneceastseednteixa-l
physicsof CR propagation and CR-driveninstability.
ceB
k (ǫ f +ǫ f ) (11c)
− p kil lj kjl li
wherequadraticordertermsinthevelocityuhavebeenne- 5 THE SIMULATION
glectedandwehaveomittedtermsinvolving∂f/∂ptimesa
gradient of u apart from the term in equation (11a) for the A full 3D simulation with realistic parameters is not possi-
evolution of f . This amounts to the neglect of second or- ble because of the large ratio of the largest distances (the
0
derFermiaccelerationandaccelerationresultingfromshear CR precursor scaleheight and the CR free escape distance)
to the smallest distance (the shortest wavelength on which
motionsinthebackgroundhydrodynamics.Inourproblem,
theseprocessesaresmallincomparisonwithaccelerationat the NRH instability grows). The correspondingly large ra-
theshock. tio of the largest timescale (the SNR expansion time) to
theshortesttimescale(theshortestNRHgrowthtime)simi-
We have chosen to terminate the set of equations at
equation(11c)bysettingf tozero.Additionally,wesoften larlymakessubstantialdemandsoncomputerresources.The
ijk
computational constraintsare discussed in AppendixF.
the termination of the tensor expansion by replacing the
We artificially increase the magnetic field and stretch
magneticrotationinequation(11c)byadampingtermwith
a damping rate equal to the magnetic gyration frequency. otherparameters in a favourable direction, choosing:
The logic of this approximation is that random rotation of n =0.1cm−3 u =60,000km s−1
e s
the stress tensor anisotropy leads to its damping. A similar
assumption underlies Bohm diffusion which replaces rota- B =47µG v =3 105m s−1 T =100TeV (14)
tionoffi byadampingrate,therebyterminatingthetensor 0 A × inject
expansionatf whereasweterminateitatf .Ourapprox- where T is the energy at which CR are injected at the
i ij inject
imation makesintuitivesense, and we support it in Appen- shock. The zeroth order magnetic field is aligned along the
dicesDandEbyexaminingitseffectonpropagatingmodes shocknormal.Theinstabilityisseededbyimposingrandom
and on the NRHinstability. fluctuationsonthemagneticfieldwithameanmagnitudeof
(cid:13)c 2013RAS,MNRAS000,1–17
Cosmic ray acceleration and escape from supernova remnants 7
9µG.CRareinjectedattheshockintothelowestmomentum simulation to initialise the upstream background plasma at
bin (width ∆p) according to therule rest relative to thecomputational grid. The more usual ap-
proachofinitiatingthesimulationbysettingthebackground
∂f
4πp2∆p cp ∂t0 =−constant×∇.u min(ρ0u2,U) (15) plasmainmotiontowardsareflectiveboundarywouldexac-
erbate the effect of numerical diffusion on small structures
where ρ0 is the density upstream of the shock and u is the duringadvection.Instead,adensepiston isinitialised mov-
local background fluid velocity relative to the initially sta- ingleftwardfromtherightboundary,pushingashockbefore
tionaryupstreamplasma.U isthelocalthermalenergyden-
it into thestationary backgroundplasma.
sity.Thisprescriptionisdesignedtoinjectasuitableenergy
density of CR at the shock, dependent on the choice of the
constant,whilstavoidinganegativevalueofU duetoexces-
sive transfer of energy to CR from cold plasma at the foot
oftheshock.TheresultingCRcurrentdensityaheadofthe
shock is displayed in panel (e) of figures 2 and 3. The peak
CR current density in the population freely escaping ahead
of theshock in figure2 is 6 SIMULATION RESULTS IN 3D
j =1.1 10−14Ampm−2 (16) Results from the 3D simulation are presented in figures 2
CR
× and 3 at t=4.1 107sec and t=6.2 107sec respectively.
This current densitycorresponds to η=0.026, giving Note that the ho×rizontal spatial scale×s, in the direction of
γ−1 =2.3 106sec cγ−1 =7 1014m (17) the shock normal, are artificially compressed by a factor
max × max × of 24. The actual aspect ratio of the computational box is
k−1 =7 1011m r =7 1013m 177:1. The plot of the background plasma density in panel
max × g × (a)showsthepositionofthedenseplasmapistonpropagat-
ing leftwards and pushing the shock ahead of it. The high
r k =100 M =200 (18)
g max A pressureregion in panel(b)isduetoplasma heatingat the
We use a spatial grid with ∆x = ∆y = ∆z = 1.4 1012m. shock. The CR pressure is plotted in panel (c). The maxi-
×
Tencellsinmomentumcoveranenergyrangefrom 100TeV mumCR pressure occursat theshock. Panel(c) of figure2
to1PeVwithlogarithmicspacing.Thereare32cellsineach shows the CR population dividing into a population freely
of x and y with periodic boundary conditions. In contrast propagating ahead of the shock and a population confined
5676 cells are used in z with reflective boundaryconditions by magnetic field (panel (d)) near the shock. CR confined
for the MHD part of the code and for CR at the right neartheshockcontinuetobeaccelerated.Thebandedstruc-
hand boundary. CR reaching the left hand boundary are ture in the CR pressure ahead of the shock in panel (c) is
disposed of on theassumption that theyescape freely. Cor- exaggeratedbythecolourcoding.Itrepresentsonlyavaria-
respondingly,thecomputationalboxextends7.9 1015mby tionofafewpercentintheCRpressureandprobablyarises
4.5 1013mby4.5 1013m.Thebox-sizeinxan×dy iscom- from a small oscillation in the injection rate at the shock.
× ×
parable with the initial CR Larmor radius, but the Larmor Also, the horizontal spatial compression of figure 2 distorts
radiuscontractssignificantly asthemagnetic fieldisampli- theaspect ratio of thebands.
fied.Thecell-sizeisπ−1 timesthewavelengthofthefastest TheescapingCRareanessentialaspectoftheprocessof
growing mode, which is barely sufficient to represent the CR confinement by the self-generated magnetic field. They
initialgrowthoftheinstability,butthecharacteristicscale- excite the instability and carry the escaping areal charge
length of the instability increases rapidly as its amplitude QCR = jCRdt identified in equation (1) as necessary for
grows. These parameters are only marginally sufficient to substantial field amplification. In figure 2 the escaping CR
R
representthephysics,butforexampleahalvingofthecom- are only just reaching the left hand of the grid where the
putational cell-size would increase the computational cost integral jCRdtissmallandfieldamplificationisnegligible.
by a factor of 16. The simulation in figure 3 was run for In contrast, immediately in front of the shock in figure 2
R
6.2 107sec = 2yr using 128 processors for 75 hours on jCRdthasreachedthecriticalvalueneededforstrongfield
×
the SCARF-LEXICON cluster at the UK Rutherford Ap- amplification causing theonset of CR confinement.
R
pleton Laboratory. The MHD part of the code is the same By the time of figure 3 most of the escaping popula-
as that used by Lucek & Bell (2000) and Bell (2004). The tionhaspassedthroughthefreeescapeboundaryattheleft
VFPpartofthecodethatmodelstheCRusesasecond or- hand end of the grid. By thistime a large magnetic field in
derRunge-Kuttaschemeinconfigurationspace,adonor-cell the upstream plasma has switched off CR escape creating
schemeinmagnitudeofmomentum,andtheBorisalgorithm an expanseof low CR densitybetween theconfinedand es-
for rotation in momentum due to magnetic field as used in caping CR populations. Panel (d) shows that the magnetic
particle-in-cell codes (Birdsall & Langdon 1985). The VFP field is amplified by an order of magnitude by the escaping
equations are formulated in tensor notation rather than in CR over alarge distance ahead of theshock.
sphericalharmonics,butthetwoformulationsaresimilarfor The separation of CR into escaping and confined pop-
the low order expansion used here, and their numerical so- ulationsmatchesexpectationsfromtheargumentinsection
lution isinformed byexperiencewith theKALOSspherical 3. The CR charge per unit shock area in the escaping pop-
harmoniccode(Belletal,2006) whereRunge-Kuttaadvec- ulation is 1 10−7 Coulomb m−2 in figure 2 in good
tion is found toberobust, sufficientlyaccurate, and a good agreement∼with×the estimate of 1.3 10−7 Coulomb m−2
×
fit to the form of the equations. Numerical diffusion due fromequation(1).Panels(e)infigures2and3plottheCR
tolimited spatialresolution isameliorated bydesigningthe currentdensity ahead of theshock.
(cid:13)c 2013RAS,MNRAS000,1–17
8 A.R. Bell, K.M. Schure, B. Reville and G. Giacinti
Figure2. 2Dslicesofa3Dsimulationwithpeakvaluesinbrackets:electrondensity(9.7cm−3),pressure(340nPa),CRpressure(2.3nPa),
magnitudeofthemagneticfieldperpendiculartotheshocknormal(610µG),CRcurrentparalleltotheshocknormal(1.1×10−14Am−2).
t=4.12×107sec(1.3years).Thedimensionsofthecomputational boxare7.9×1015mby4.5×1013m.Thehorizontalspatialscale,in
thedirectionoftheshocknormal,isartificiallycompressedbyafactorof24.
Figure3. 2Dslicesofa3Dsimulationwithpeakvaluesinbrackets:electrondensity(6.9cm−3),pressure(320nPa),CRpressure(1.9nPa),
magnitudeofthemagneticfieldperpendiculartotheshocknormal(590µG),CRcurrentparalleltotheshocknormal(4.2×10−15Am−2).
t=6.18×107sec(2.0years).Thedimensionsofthecomputational boxare7.9×1015mby4.5×1013m.Thehorizontalspatialscale,in
thedirectionoftheshocknormal,isartificiallycompressedbyafactorof24.
(cid:13)c 2013RAS,MNRAS000,1–17
Cosmic ray acceleration and escape from supernova remnants 9
Figure4. 2Dsimulations.Plotsoff0p3atCRenergies100,167,278and464TeVattimesbetween4.1×107sec(toprow)&16×107sec
(bottom row),withcomputational boxlengthsLproportionaltotbetween 7.9×1015m&3.2×1016mchosensuchthattheCRtravel
thelengthoftheboxduringthesimulationtime.Thewidthofthecomputational boxis4.5×1013mineachplot.Theverticalaxesare
indimensionlessunitswithpeakvaluesgivenagainsttheaxis.
7 SIMULATION RESULTS IN 2D subsequent rows of figure 4 the spatial dimensions normal
to the shock are expanded by factors of 2, 3 and 4 such
Figure 4 follows the calculations of figures 2 and 3 to later
thatCRtravelthelengthofthegridintherespectivetimes
times with thesame parametersbut in 2D instead of 3D to (t=8.2 107, 1.2 108 & 1.6 108sec).
reduce the demand on computer resources. The top row of × × ×
plots in figure 4 is the 2D equivalent of the 3D results in In 3D at t = 4.1 107sec the dip in the CR density
×
figure 2 at t = 4.1 107sec. The dimensions of the spatial separating the escaping and confined CR has only recently
×
grid are thesame in figure2and thetop row of figure4. In formedasseeninfigure2(c).Atthesametimein2D(figure
(cid:13)c 2013RAS,MNRAS000,1–17
10 A.R. Bell, K.M. Schure, B. Reville and G. Giacinti
4 top row) the separation between the confined and escap- R,CR are confinedif
ingpopulations is relatively more developed.Otherwise the
resultsaresimilarin2Dand3D.Thereasonfortheslightly R ηρ(r)u2s(r) r2dr=10R2 ρ(R) (19)
earlier development of CR confinement in figure 4 is that Z0 T(r) s µ0
magnetic structures can expand more freely in 2D. In 2D,
Differentiating this equation with respect to R and assum-
magnetic structures expand in the ignored dimension with-
ing a power law dependence of density on radius, ρ(R) =
out running into stationary or counterpropagating plasma. ρ (R/R )−m,
0 0
Nonlineardevelopmentin 2D istherefore less inhibitedand
themagneticfieldgrowsmorerapidly.Themagneticfieldin T(R)= η√µ0 u2R√ρ (20)
2D reaches a steady value similar to that found in the 3D 5(4 m) s
−
calculation.Theeventuallimitonthemagnitudeofthemag- Defining u = u /10,000km s−1, R = R/parsec, η =
7 s pc 0.03
netic field may be set by magnetic tension which operates η/0.03, n =ρ/2 10−21kg cm−3 (such that n is approxi-
e e
equally in 2D and 3D. mately theelectro×n densityin cm−3),andtakingm=0for
Comparison of the CR distributions at different times
expansion intoa uniform medium,
in figure4shows thatat early times(t=4.1 107sec) only
the low energy CR are confined. At t = 8.2× 107sec CR T =230 η0.03n1e/2u27Rpc TeV (21)
×
areconfinedat100and167TeVwithconfinementbeginning ASNRwith u =0.6, n =1andR =1.7, representative
7 e pc
to occur at 278TeV and no evidence of the development of of Cas A, would then accelerate CR to 140TeV. A SNR
two separate populations in the few CR reaching 464TeV. with u = 0.5, n = 0.1 and R = 10,∼representative of
7 e pc
By t = 1.2 108sec a larger number of CR have reached SN1006, would accelerate CR to 180TeV. These energies
464TeV with×a local minimum in CR density f0p3 evident areafactorof10lower thanthee∼nergyoftheknee.Abbasi
at all energies up to 464TeV in front of the shock. By t = et al (2012) place the knee at 4-5PeV, although data from
1.2 108sec CR at 100 and 167 TeV are strongly confined otherexperimentsindicatea lower energy andtheturnover
×
withashortprecursorscaleheightaheadoftheshock,which in the spectrum is not well defined as shown in their figure
is unsurprisingsince theLarmor radiusof a 100TeV proton 15.
ina700µG magneticfield is5 1012m.By t=1.6 108sec Our model places a considerable question mark over
× ×
CR are confined at all energies upto 464TeV. the ability of the well-known historical SNR to accelerate
At all four times in figure 4 the CR spectrum is much CR to the knee. Acceleration by SNR such as Cas A and
steeper than the steady-state test particle spectrum (f0 SN1006fails toreach thekneein ouranalysisbecausetheir
∝
p−4)becauseof CRloss intothedownstream plasma. Mag- expansion is already significantly decelerated. Zirakashvili
netic field amplification takes place ahead of the shock and &Ptuskin(2008a) also reach theconclusion (theirTable2)
the downstream field only becomes large when the shock that the historical SNR do not accelerate CR beyond 100-
overtakes the field amplified in the upstream. Hence at 200TeV,butnotethattheirparameterη differsfromour
esc
early times CR escape downstream. This causes a reduc- parameter η by a factor of two (η =2η).
esc
tion in thenumberof CR accelerated to high energy at the InitialexpansionvelocitiesofveryyoungSNRcanreach
shock, thereby steepening the spectrum. Downstream con- 30,000 km s−1 (Manchester et al 2002) or possibly higher
finement at the shock improves at later times as shown by for some typesof SN (Chevalier& Fransson 2006). Accord-
the downstream gradients in f0 at t = 1.2 108sec and ing to equation (21) expansion into a density of 1 cm−3 at
×
t= 1.6 108sec. However, unphysical numerical relaxation 30,000 km s−1 for 16 years to a shock radius of 0.5 parsec
×
of the spatially compressed downstream magnetic field due accelerates CR to 1PeV. Moreover, the energy processed
tolimitedspatialresolutionmaybeplayingapartinallow- throughtheshocki∼scomparabletothatforexpansionto1.5
ing CR to escape downstream. parsec at 6,000 km s−1, therebycontributinga comparable
There is slight evidence of pulsed acceleration as seen CR energy content to the Galactic energy budget. A com-
in the three separate peaks in the 278TeV CR density at plementaryperspectiveonthesameproblemisobtainedby
t = 1.2 108sec and t = 1.6 108sec but overall the CR expressing the maximum CR energy in terms of the mass
× ×
profiles develop in an orderly manner. M =4πρR3/3sweptupbytheshockandthecharacteristic
energy of the blast wave E =Mu2/2. If M is the mass in
s ⊙
unitsofasolarmassandE istheenergyinunitsof1044J,
44
themaximum CR energy is
T =0.5η n1/6E M−2/3 PeV (22)
0.03 e 44 ⊙
8 APPLICATION TO SNR
which indicates that a typical SNR should be able to ac-
Our simulations support the model developed in section 3. celerate CR to 1PeV and that the maximum CR energy
∼
According to our model CR are confined and accelerated if is greater for the same energy E given to a smaller mass
the electrical charge of CR escaping upstream of the shock M.ThemaximumCRenergyisnearlyindependentofden-
reachesQ =10 ρ/µ Coulomb m−2.Wenowapplythis sity n because a shock expands to a larger radius in a low
CR 0 e
tosphericalSNRshocks.TheCRcurrentdensityataradius densitymedium before deceleration (Hillas 2006).
p
Risj =ηρu3r2/R2T duetoCRacceleratedtoenergyeT Intheself-similarSedovphase,E isconstant,M ispro-
CR s
when the shock radius was r. Sinceonly thehighest energy portiontoR3andhenceT R−2.ThemaximumCRenergy
∝
CR escape upstream we assume that the CR reaching the decreases with radius during the Sedov phase until mag-
radius R are monoenergetic with energy eT. T evolves as netic field amplification ceases, at which point our analysis
theshockexpands.WhentheSNRshockreachestheradius intermsofthechargecarriedbyescapingCRisinapplicable.
(cid:13)c 2013RAS,MNRAS000,1–17