Nonlinear particle acceleration and thermal particles in GRB afterglows

0.75 MB·

by  Donald C. Warren

Draft version January 17, 2017
PreprinttypesetusingLATEXstyleemulateapjv.5/25/10
NONLINEAR PARTICLE ACCELERATION AND THERMAL PARTICLES IN GRB AFTERGLOWS
Donald C. Warren,1 Donald C. Ellison,2 Maxim V. Barkov,1,3,4 Shigehiro Nagataki1,5
Draft version January 17, 2017
ABSTRACT
7
The standard model for GRB afterglow emission treats the accelerated electron population as
1
0 a simple power law, N(E) ∝ E−p for p (cid:38) 2. However, in standard Fermi shock acceleration a
2 substantialfractionoftheswept-upparticlesdonotentertheaccelerationprocessatall. Additionally,
ifaccelerationisefficientthenthenonlinearbackreactionofacceleratedparticlesontheshockstructure
n
modifies the shape of the non-thermal tail of the particle spectra. Both of these modifications to
a
the standard synchrotron afterglow impact the luminosity, spectra, and temporal variation of the
J
afterglow. Toexaminetheeffectsofincludingthermalparticlesandnonlinearparticleaccelerationon
6
afterglowemission,wefollowahydrodynamicalmodelforanafterglowjetandsimulateaccelerationat
1
numerouspointsduringtheevolution. Whenthermalparticlesareincluded,wefindthattheelectron

population is at no time well-fitted by a single power law, though the highest-energy electrons are;
]
E if the acceleration is efficient, then the power law region is even smaller. Our model predicts hard-
H soft-hard spectral evolution at X-ray energies, as well as an uncoupled X-ray and optical light curve.
Additionally, we show that including emission from thermal particles has drastic effects (factors of
.
h 100and30,respectively)ontheobservedfluxatopticalandGeVenergies. ThisenhancementofGeV
p emission makes afterglow detections by future γ-ray observatories, such as CTA, very likely.
- Keywords: acceleration of particles — ISM: cosmic rays — gamma-ray bursts — magnetohydrody-
o
namics (MHD) — shock waves — turbulence
r
t
s
a
1. INTRODUCTION on in them, and that a great deal of work has been done
[
Gamma-ray bursts (GRBs) are brief flashes of in service of this goal.
1 gamma-ray photons occurring a few times a day across The traditional picture of GRB afterglows is that
v they are synchrotron emission from a population of
the entire sky. After the high-energy prompt emission
0 energetic electrons accelerated by a relativistic exter-
fades, broadband (usually radio to X-ray, but occasion-
7 nal shock. This model for afterglow emission predates
ally GeV gamma-ray) afterglow emission is sometimes
1 the first afterglow detection (Paczyn´ski & Rhoads 1993;
observed. These afterglows are less numerous, though
4 M´esz´aros & Rees 1997). Early analytical models treat-
increasingly frequent as technology improves: the Swift
0 ing the time-dependence of the afterglow invoked syn-
mission (Gehrels et al. 2004) alone reports over 300
1. GRBs with afterglows detected both in optical/UV and chrotron radiation from a single power law of electrons
0 in X-ray6, and of the ∼ 1800 GRBs reported between (Sari et al. 1998; Chevalier & Li 2000). Though very
7 1997and20157 morethan1000weredetectedinatleast efficient in terms of computation time, these models are
1 heavily dependent on assumptions about, e.g., the mag-
one of radio, optical, or X-ray—with 84 detected in all
: netohydrodynamicstateofthefluidinthevicinityofthe
v three8. Afterglowobservationsofferawealthofinforma-
blast wave, the presence of a reverse shock, whether the
i tionabouttheprogenitorsofGRBs, theirenvironments,
X environment of the GRB was wind-like or of constant
the microphysics of particle acceleration, and more. A
density, and many more. The numerous permutations,
r fullaccountingofthephysicsofafterglowsisbeyondthe
a and their consequences for both single-time spectra and
scope of this work, but see Kumar & Zhang (2015) for
long-termlightcurves,arecollectedinGaoetal.(2013).
an extended review of the current state of the field. Suf-
Numericalstudieshavealsobeenundertaken. Numerical
fice it to say that correctly interpreting observations of
hydrodynamics allows more accurate handling of emis-
afterglows is critical to understanding the physics going
sion(1)duringthetrans-relativisticphaseofshockdecel-
eration(vanEertenetal.2010),(2)fromareverseshock
1Astrophysical Big Bang Laboratory, RIKEN, Saitama 351-
(Uhm et al. 2012), (3) due to the multi-dimensional na-
0198,Japan;[email protected]
2Physics Department, North Carolina State University, Box tureoftheblast-wave(vanEertenetal.2012;vanEerten
8202,Raleigh,NC27695,U.S.A. & MacFadyen 2012), and in many other scenarios.
3InstituteforPhysicsandAstronomy,UniversityofPotsdam, Despite the multitude of theoretical and numerical
11476Potsdam,Germany
4DESY,Platanenallee6,15378Zeuthen,Germany analyses and the breadth of situations considered, virtu-
5Interdisicplinary Theoretical Science Research Group ally all work on afterglow emission shares one common
(iTHES),RIKEN,Saitama351-0198,Japan assumption: thattheelectronsresponsibleforthebroad-
6 http://swift.gsfc.nasa.gov/archive/grb table/stats/ bandemissionformapowerlawwithasinglespectralin-
7 http://www.astro.caltech.edu/grbox/grbox.php
dexp(thisletterwillalsobeusedthroughoutthispaper
8 http://www.mpe.mpg.de/jcg/grbgen.html
2 Warren et al.
to refer to momentum, but context should make it clear experience a larger velocity difference between the up-
whichdefinitionappliesatanyparticularinstance). The stream and downstream limits of their scattering, which
rare exceptions considered a Maxwellian thermal distri- in turn causes a harder spectrum at higher particle en-
bution instead of (Pennanen et al. 2014) or in addition ergies (see, e.g., discussion in Berezhko & Ellison 1999;
to (Giannios & Spitkovsky 2009) a power law tail.9 Blasi et al. 2005).
The value for the electron spectral index p, starting WhilemostworkonnonlinearFermiaccelerationhas
from basic shock acceleration principles, is 2.23 (Kirk been done for ions accelerated in nonrelativistic shocks,
et al. 2000; Keshet & Waxman 2005), in agreement with it is important to emphasize that nonlinear effects must
the results of Monte Carlo simulations (Ellison et al. occur in relativistic shocks as well, if acceleration is ef-
1990;Bednarz&Ostrowski1998;Achterbergetal.2001; ficient. Furthermore, if no mechanism acts to transfer
Ellison et al. 2013). Strictly speaking, this result applies energy from ions to electrons, Fermi acceleration places
onlyforultra-relativisticshocksinthecasewheretheup- a large majority of the shock energy into ions. Ellison
streammagneticfieldandshocknormalareparallel,and et al. (2013) have shown that the nonlinear effects re-
then only in the limit where accelerated particles carry quired to ensure flux conservation in relativistic shocks
negligible energy (the so-called “test particle” regime). may lead to orders of magnitude depression of electron
However, particle-in-cell (PIC) simulations show that acceleration efficiency (Warren et al. 2015). Handling
self-generated magnetic field turbulence can overwhelm the interaction between shocks and the particles they
weak pre-existing fields (those for which the magneti- accelerate is of paramount importance to properly un-
zation parameter of the inflowing plasma is σ (cid:46) 10−5, derstanding GRB afterglows.
with σ ≡ B2/(4πγ n m c2) for a shock moving with In addition to the electron spectral index p, GRB
0 0 0 p
Lorentz factor γ into plasma with mean magnetic field models use two parameters to handle the microphysics
0
B and proton number density n ) at relativistic shocks present at the shock: the fraction of energy density in
0 0
(Sironi et al. 2013), and the results are very similar to electrons((cid:15)e;seediscussionsurroundingEquation5)and
those at parallel or entirely unmagnetized shocks (e.g., magnetic field ((cid:15)B). Applying the standard synchrotron
Spitkovsky 2008). model for afterglows shows that the value of (cid:15)e in most
ObservationsofGRBafterglowssuggestsofterspec- bursts is 0.1−0.3 (Santana et al. 2014). The value of
tra than predicted by test-particle theory. Curran et al. (cid:15)B is much less constrained and may vary over several
(2010)attemptedtofitX-rayobservationsbothtoasin- orders of magnitude, from 10−7 to 10−3 (Santana et al.
gle universal power law and to a distribution of p val- 2014; Wang et al. 2015). If a particularly hard electron
ues. Their most likely single value for p was 2.25, in spectral index p (cid:46) 2.4 is chosen, (cid:15)B may even approach
close agreement with theory, but statistical tests reject- 10−2. Theserangesaresignificantlylargerthanthestan-
ed a universal p value. When a Gaussian distribution dard interstellar medium (ISM) value of (cid:15)B ∼10−9, and
of spectral indices was used instead, the best fit was require amplification of the magnetic field by means of
p=2.36±0.59. The large uncertainty is almost entirely plasma instabilities.
duetoalongtailofsofter—ratherthanharder—particle All of the preceding discussion of (cid:15)e and (cid:15)B rests on
spectra. Asearchoveralargerparameterspace(varying, the assumption that the two parameters are constant
e.g., explosion energy, electron spectral index, jet open- throughout the life of the afterglow. This is the de-
ing angle, and energy density in electrons and magnetic fault state for afterglow modeling, though alternatives
field)foundthatGRBafterglowsspannedtheentirecon- have been explored in a handful of cases. Rossi & Rees
sidered range of spectral indices, as soft as p=5 (Ryan (2003) explored a basic two-zone model for the exter-
etal.2015). Forlackofauniversalvalueforpstemming nalshockandafterglow,allowingfordifferentconditions
from observations, the traditional value used in studies at and downstream from the shock. Letting (cid:15)e and (cid:15)B
of GRB afterglows is p=2.4−2.5. vary with time has been proposed to explain both the
It is well known that at non-relativistic speeds, ac- plateau phase of many X-ray light curves (Ioka et al.
celeratedparticlesmaynotformasimplepowerlaw. Ef- 2006) and the chromatic spectral breaks seen in many
ficient acceleration of particles10 moves such shocks be- afterglows (Panaitescu et al. 2006). The two parame-
yondthetest-particleregime,astheacceleratedparticles ters have been fitted to decaying power laws in order to
influence the structure of the shock doing the accelerat- explore rebrightenings associated with a wind-ISM in-
ing. This nonlinear feedback loop results in a shock pre- terface (Kong et al. 2010). The possibility that (cid:15)e and
cursor, which acts to slightly slow the inflowing plasma (cid:15)B might increase with time was considered in Granot
beforeitcrossestheshock. Moreenergeticparticlesthus etal.(2006),asameanstoreducetheγ-rayefficiencyof
the prompt emission below the worryingly large values
9 The work of Eichler & Waxman (2005) also treated the pos- of 0.5−0.9. In each of these cases, however, the values
sibility of separate thermal and accelerated electron populations, for (cid:15)e and (cid:15)B were either chosen to match observations
but used unreasonably low values for the thermal electron energy or given a specific functional form.
becausecrucialnumericalresultshadnotyetbeenpublished.
No work of which we are aware has calculated (cid:15)
10 We must be clear, here: if even 5% of the kinetic energy of e
or (cid:15) from first principles over a large portion of an
inflowingparticles(viewedintheframewheretheshockisatrest) B
istransferredtoahigh-energydistributionofacceleratedparticles, afterglow.11 In Lemoine (2013), though, a functional
then the shock is beyond the test-particle regime. This is lower
thanthe∼10%observedinlong-termelectron-ionPICsimulations
(Sironietal.2013). 11 PIC simulations such as Sironi et al. (2013) and Ardaneh
NLDSA model for GRB afterglows 3
form for decaying microturbulence was assumed based limit. This assumption fails in the immediate vicinity
on PIC results. The synchrotron spectrum of particles of a relativistic shock, as the magnetic field is turbu-
downstream from the shock was calculated using that lent down to length scales of the ion skin depth, which
functional form in a variety of conditions (e.g., gradu- may be many orders of magnitude smaller than the par-
al/rapiddecayofturbulence,weak/stronginverseComp- ticle gyroradius (e.g., Sironi et al. 2013; Plotnikov et al.
tonlosses,etc.). Thelongitudinalextentoftheenhanced 2013). In these conditions the mean free path approach-
turbulence was some 102 −104 ion skin depths. While es λ ∝ p2 rather than ∝ p1 (e.g., Plotnikov et al.
mfp
this is a large distance for PIC simulations, it is ignor- 2011). If, on the other hand, the power spectrum of the
ably small compared to the volumes being simulated in magnetic turbulence extends to length scales compara-
this work. The length scale of our Monte Carlo code is ble to r , then scattering proceeds as p1 at all momenta.
g
r = γ u m c/(eB ), which is σ−1/2 times larger than We do not attempt to address these complexities at this
g0 0 0 p 0
the plasma skin depth; the emission shells are individu- time; previous work on this topic has shown that this
allythousandsofr insize. Here,u istheshockspeed, simple model for particle scattering can adequately re-
g0 0
m is the proton mass, c, and e is the electronic charge. produce some important PIC results for injection and
p
The aim of this paper is to build a physically- early acceleration (Ellison et al. 2013). Ellison et al.
motivated model of a GRB afterglow, from early in the (2016) discussed the effect of a modified diffusion coeffi-
afterglow until very nearly the transrelativistic stage of cient on both CR spectra and nonlinear smoothing, and
deceleration. We pay particular attention to the nonlin- we refer the reader to that paper for more information.
ear effects, from the simultaneous acceleration of elec- For now we simply note that the region of intense
trons and protons, linking the shock structure and the turbulence is contained to a small region around the
cosmic rays being accelerated. Where needed, we in- shock,andsoweassumethattheeffectsofthisregionon
corporate results from PIC simulations to guide our theoverallmodelareminimal. Intheupstreamdirection
parametrization of key microphysics. We calculate pho- the small-scale Weibel instability is contained to a pre-
ton spectra for various times, which we compare against cursorextendingashortdistanceaheadoftheshock(see
observationsintheoptical,X-ray,andgamma-raybands. Equation 6 in Sironi et al. 2013). Downstream from the
The paper is organized as follows. In §2 we discuss shock the turbulence decays rapidly, inverse-cascading
the Monte Carlo code used to generate the particle and tolargerwavelengthsanddecreasinginstrengthast−1/2
photon spectra to be presented later. The Monte Carlo (Lemoineetal.2013,andreferencestherein). Bothscales
code is applied post-process to a hydroynamical base to aredramaticallysmallerthanthedistancesconsideredin
createasequenceofafterglowsnapshots,inamannerwe this paper.
describe in §3. We present our example afterglows and We assume pitch-angle scattering for the particles
compare against observations in §4, and conclude in §5. as they move in the turbulent magnetic field. Under the
pitch-anglescatteringmodel, aparticlewithvelocityv
pt
2. THEMONTECARLOCODE travelsforatimeδt(cid:28)t betweenscatterings. Thequan-
c
The work presented here relies on a Monte Carlo tity t =λ /v is the “collision time,” the fluid-frame
c mfp pt
codetreatingfirst-orderFermiaccelerationinthevicini- time to accumulate a cumulative deflection on the order
tyofaninfinite,steady-state,parallel,collisionlessshock. of 90◦. Our implementation of pitch-angle scattering re-
Thecodefollowsthescatteringandpropagationofparti- quires a parameter, N , which controls the fineness of
g
cles in a turbulent magnetic field around a velocity gra- scattering. Each gyroperiod τ = 2πr /v is divided
g g pt
dient. Momentum and energy fluxes are tracked dur- into N time steps, so δt = τ /N . The maximum de-
g g g
ing propagation, which may then be used to modify the flection a particle may experience at a single scattering
shock structure in an attempt to satisfy the Rankine– event is also controlled by N (Ellison et al. 1990),12
g
Hugoniot conditions at all locations. The process has
(cid:113)
beendescribedindetailin, e.g.,Ellisonetal.(2013)and δθ = 12π/N . (2)
max g
Warren et al. (2015), and the reader is referred to those
papers for a more complete description of the code. In Large values of N imply that only modest deflections
g
this section we briefly discuss key features, and depar- occur at scattering events, while small values of N sug-
g
tures from previous work. gest large-angle scattering (Ellison et al. 1990; Stecker
All particles—regardless of energy or species—are et al. 2007; Summerlin & Baring 2012).
assumed to scatter elastically off turbulence in the lo- The code does not attempt to calculate the wave–
cal fluid frame with a mean free path given by particleinteractionsthatgovernparticlemotioninphys-
ical shocks. Equations 1 and 2 control the transport
λ =r , (1)
mfp g of particles within the shock structure. As such, we
wherer =pc/(ZeB)isthegyroradiusofaparticlewith assume that sufficient turbulence exists on the correct
g
charge number Z in the local magnetic field B. That is, length scales to drive that transport. Further, note
we assume that particle scattering occurs in the Bohm that both equations are position-independent: though
etal.(2015),ofcourse,allowfordirectcalculationofboth(cid:15)e and 12 This equation is specific to the choice of λmfp made in
(cid:15)B. However,thesesimulationsarelimitedbycomputationalpow- Equation 1. More generally, if a particle’s mean free path is a
er to studying time periods vastly shorter than the timescales of function f(rg) of its gyroradius, Equation 2 may be rewritten
(cid:112)
afterglows. δθmax= 12πrg/(Ngf(rg)).
4 Warren et al.
position-dependentdiffusioncoefficientshavebeenincor- of electrons in addition to protons. Since Monte Car-
porated in Monte Carlo calculations in the nonrelativis- lo transport allows for arbitrary particle anisotropy, it
tic regime (Vladimirov et al. 2009; Bykov et al. 2014), has unique advantages over semi-analytic methods that
we defer that particular extension of our code to future avoid direcly modeling thermal particles. Furthermore,
work. the simple parametrization of complex plasma processes
The assumption of position-independent scattering, allows modeling Fermi acceleration, with internally con-
even in the upstream region, is somewhat at odds with sistent injection, to energies not yet accessible with PIC
theoretical predictions (Lemoine & Pelletier 2010) and simulations.
PIC simulations of relativistic shocks propagating into For the basic scattering assumptions, i.e. Equa-
moderatelystrongupstreammagneticfields(Sironietal. tions 1 and 2, this injection model generally predicts ef-
2013). (In the simulations of the latter work, compu- ficient injection and acceleration incompatible with the
tational costs precluded testing the result at the very test particle limit (e.g., Warren et al. 2015). If a sharp,
low upstream magnetic field energy fractions expected unmodified shock is assumed, then energy and momen-
in GRB afterglow shocks.) In both of those works, the tumfluxesarenotconservedwithintheshockstructure:
longitudinal extent of the turbulent upstream precursor the flux downstream from the shock will always exceed
is a few times r (≡ γ m u c/(eB ), the gyroradius of the incoming flux from far upstream.
g0 0 p 0 0
a proton with speed u in a field of strength B ). When In order to discuss test-particle shocks, we have in-
0 0
the upstream magnetic field is ordered rather than tur- troduced a new parameter, P ∈ [0,1]. When a ther-
inj
bulent, the upstream dwelltime of particles increasesby mal particle downstream from the shock would recross
a factor of at least a few (Achterberg et al. 2001). For the shock, P is the chance that it will be injected into
inj
CR acceleration happening over a given (and limited) acceleration process. If the particle fails that check, it is
timeperiod,then,theturbulent-fieldcaseallowsforfew- reflected back downstream to continue scattering. Note
er shock crossing cycles, as particles can scatter further that P =0 corresponds to disabling all diffusive shock
inj
upstreamfromtheshockeverytimetheycrossbackinto acceleration (e.g. Figure 1 in Ellison et al. 2013), while
the upstream region. Protons, which do not experience P =1isthethermalleakagemodelalreadymentioned.
inj
any sort of loss process in our simulation, are limited For any value of P , the code can determine the
inj
in maximum energy entirely by acceleration time. Our self-consistent, modified shock structure, including the
assumptionofaturbulentupstreamregionthereforesup- back-reactionofacceleratedparticles. Typically,themo-
pressesthemaximumenergyoftheprotonspectrum,and mentumandenergyfluxesareconservedtowithin5%ev-
there is good reason to suspect that the CR tail would erywhere in the shock structure. As was shown in, e.g.,
extend to higher energies with a more exact model of Ellisonetal.(2013),andaswillbepresentedbelow,these
the magnetic field around the shock. Electrons, on the modified shocks result in particle distribution functions
otherhand,arelimitedbysynchrotronlossesincurredin thatdeviatesignificantlyfromthetraditionally-assumed
the downstream region; their spectrum is unaffected by power law. The efficient injection we see is nonetheless
turbulence upstream of the shock. consistent with PIC results.
Particle injection is handled via a thermal leakage Inrelativisticshocks,theWeibelinstabilitymediates
model. As unshocked particles first cross the shock and energytransferfromionstoelectronsintheupstreamre-
interact with the downstream plasma, they will scatter gion (Sironi et al. 2013; Ardaneh et al. 2015). The effect
inthedownstreamplasmaframeandbeenergized. Since of this energy transfer is that electrons cross the shock
the downstream flow is subsonic, this “thermalization” carrying as much as 40% of the kinetic energy flux. A
resultsinasignificantfractionofparticleshavingaplas- basic scheme for approximating this energy transfer was
maframespeedgreaterthanthebulkfluidspeed, which implementedin(Warrenetal.2015),uponwhichweiter-
is only mildly relativistic regardless of γ . The Monte ate here. Previously kinetic energy was exchanged only
0
Carlo code follows these particles as they scatter back at the shock itself, simplifying away that the transfer
across the subshock into the upstream region. A par- happens during passage through the Weibel-dominated
ticle that re-enters the upstream region will be further regions of the precursor. In this work thermal kinetic
accelerated as it scatters back downstream, which is our energy is subtracted from protons as they pass through
definition of injection. Implicit in this definition is the the precursor, and added as thermal energy to electrons
assumption that the subshock is transparent; in other atthesameposition. Althoughthecodeonlypropagates
words, mildly superthermal particles can cross without one particle at a time through the shock structure, this
experiencing forces from large-amplitude magnetic tur- implementation mimics simultaneous transfer.
bulence,cross-shockpotentials,orotherphenomenathat The degree of energy transfer is controlled by a pa-
may influence the injection rate. rameter f , which is the approximate fraction of pro-
trans
We do not pretend that this injection captures all ton kinetic energy that is donated to electrons (called
of the complicated plasma physics seen by PIC simula- f inWarrenetal.2015). PICsimulationssuggestthat
ion
tions of the shock layer. Instead, we parametrize the f = 0.4 for shocks with a Lorentz factor above 10,
trans
injection probability in a straightforward fashion and which value we use for all shocks simulated in this work.
connect the injection self-consistently to Fermi accelera- Below a Lorentz factor of 10, the instabilities responsi-
tion and the nonlinear shock structure. It is also trivial bleforenergytransfershouldbequenchedbythehigher
to include the simultaneous injection and acceleration transverse thermal velocities of the inflowing particles
NLDSA model for GRB afterglows 5
relative to their speed of approach (Lemoine & Pelleti- ulations,and(3)apost-processthatdeterminesboththe
er 2011), and so energy transfer should diminish. Early evolutionofdownstreammaterialandtheoverallphoton
testingshowedthatthemannerinwhichf decreases production.
trans
leavesanobservableimprintontheafterglowlightcurve,
and this may offer an important diagnostic connecting 3.1. The hydrodynamical base
observations with the plasma microphysics at the shock. The hydrodynamical picture of a GRB jet is akin to
For now, however, since we know of no PIC simulations thenonrelativisticsupernovaremnantscenario: apiston
ortheoreticalworkdescribingtheturbulentprecursorin of ejecta is driven into the circumburst medium (CBM).
the mildly relativistic and trans-relativistic regime, we This creates a forward shock propagating ahead of the
avoid this complication by taking ftrans to be constant contact discontinuity (the ejecta–CBM interface) and a
over the entire time covered by our simulations. reverse shock propagating down through the ejecta. At
Tomodelobservationsofafterglowsandotherastro- present we focus our attention on the forward shock
physicalshocks, theMonteCarlocodecalculatesphoton alone; particle acceleration may happen at the reverse
production due to the accelerated population. Previ- shock as it passes through the jet, but we defer this top-
ous work (Warren et al. 2015) included three emission ic to future endeavors.
mechanisms, which are used here also: pion decay due As long as shock-accelerated particles do not escape
to hadron-hadron collisions, inverse Compton radiation, theshockinsignificantquantities, thelocationoftheaf-
and synchrotron radiation. However, the implementa- terglow’sforwardshockismostlyindependentofacceler-
tionoftheseprocesseshaschangedbetweenWarrenetal. ation efficiency. Ambient particles crossing a relativistic
(2015) and this work; we provide further details in Sec- shock are heated to relativistic speeds and achieve an
tion 3.3. adiabatic index Γ ≈ 4/3. If particles were able to es-
2
Asignificantdifferencebetweenpreviousworkunder cape the shock in large quantities, the remaining fluid
this Monte Carlo framework and the results presented wouldbemorecompressible(loweringtheeffectiveadia-
here is that we assume that the magnetic field is turbu- baticindex)andsotheforwardshockwouldstaycloseto
lentandcompressedbytheshock. Theupstreamfieldis thepositionofthecontactdiscontinuity(Blondin&Elli-
assumedtobeturbulent,withcomponentsB0(cid:107) andB0⊥. son2001;Warren&Blondin2013;Duffell&MacFadyen
ItisastandardresultofmagnetohydrodynamicsthatB(cid:107) 2014). Another possibility is that the shock could radi-
is unchanged across a shock, while B (x) = B ·z(x), ate away a substantial portion of its kinetic energy. We
⊥ 0⊥
where z(x) = γ β /(γ(x)β(x)) is the density ratio be- will show later that both loss mechanisms are negligible.
0 0
tween the far upstream and local conditions. In com- It is therefore appropriate to use a purely hydrody-
bination, then, the compressed magnetic field strength namical model for the location of the forward shock at
within the shock structure is given by any given time. We have chosen to use the Blandford–
(cid:114) McKeesolutionforanadiabaticblastwave(Blandford&
1 2
B(x)=B + z(x)2, (3) McKee1976),withslightmodificationsfortheinitialac-
0 3 3 celerationperiodasdescribedinKobayashietal.(1999).
The location/velocity/age of the analytical Blandford–
with z(x) defined as above (e.g., Reynolds 1998). As
McKee solution matches very well with numerical re-
mentioned before, we neglect amplification due to the
sultsuntiltheforwardshockhasdeceleratedtoaLorentz
Weibel instability, and assume that the field has de-
factor of 5, and reasonably well down to the trans-
cayed to this value prior to dilution by expansion (see
relativisticregimeofγ ≈2. Downstreamfromtheshock
Section 3.2). 0
the analytical solution diverges from the numerical one
TheMonteCarloapproachasoutlinedaboveisboth
because of the approximations made during the deriva-
complex and versatile, but it does have several limita-
tion;wewilldiscussthesedepartureslater,inSection4.2,
tions. We do not make any assumption about the spec-
as they relate to the results presented in Figure 5.
trum of magnetic turbulence responsible for our parti-
TheBlandford–McKeesolutionscaleswithupstream
cle scattering, beyond that it exists and causes the pre-
density and isotropic energy, and the modifications pro-
scribed scattering. At the moment our scattering mean
posedinKobayashiet al.(1999)introduceanadditional
free paths are position-independent, and further do not
dependenceonthemaximumLorentzfactorachievedby
takeintoaccounttheexpecteddeparturefromBohmdif-
the shock. For the example afterglow we describe here,
fusion at high energies (Plotnikov et al. 2011, 2013; El-
weuseaconstant-densityambientmediumofoneproton
lison et al. 2016). In the context of GRB afterglows,
percubiccentimeter,13 anisotropic-equivalentenergyre-
however, the most severe limitation is the assumption of
leaseof1053 ergs,andamaximumLorentzfactorof400.
asteady-stateshock. Inthefollowingsectionwewillde-
Where needed, the jet is assumed to have an opening
scribe how multiple steady-state simulations are joined
angle of 5◦. The upstream temperature is taken to be
to provide a complete picture of an evolving system.
108 K to speed computation. As long as the Mach num-
3. ANEVOLUTIONARYMODEL ber of the shocks remains much greater than unity, the
Inthissectionwedescribehowasteady-stateMonte
13Despiteoccurringinthevicinityofamassivestar,asignificant
Carlo model is applied to the dynamic environment of a
fraction of GRB afterglows are more compatible with a constant-
GRB afterglow. The process has three key parts: (1) densityCBMthanwithawind-likeCBMprofile(e.g.,Panaitescu
thehydrodynamicpre-process,(2)theMonteCarlosim- &Kumar2002;Schulzeetal.2011).
6 Warren et al.
TABLE 1
Afterglow model, common parameters.
Timestepa tengb tobsc RFS γ0 RRH Ng,TP Ng,NL Pinj Colorinfiguresa
(sec) (sec) (cm)
1 7.5×106 3.0×102 2.3×1017 62.8 3.001 8×104 5×104 0.06 Magenta
2 1.0×107 9.7×102 3.0×1017 40.1 3.002 5×104 2×104 0.06 Purple
3 1.4×107 3.2×103 4.1×1017 25.6 3.003 4×104 1.5×104 0.05 Blue
4 1.8×107 1.1×104 5.5×1017 16.4 3.006 3×104 1×104 0.05 Cyan
5 2.5×107 3.5×104 7.4×1017 10.5 3.014 2×104 8000 0.04 Green
6 3.4×107 1.2×105 1.0×1018 6.7 3.033 1.4×104 6000 0.04 Yellow
7 4.5×107 3.8×105 1.4×1018 4.3 3.08 9000 4000 0.03 Orange
8 6.2×107 1.3×106 1.8×1018 2.7 3.18 5000 2000 0.03 Red
aForfiguresthatcomparedatafromdifferenttimestepsagainsteachother, eachtimestephasacolorthatwillbeconsistentlyusedto
representit.
bElapsedtimeintherestframeofthecentralengine;alsotheISMframe.
cElapsedtimesincethestartoftheGRBforanobserveronEarth,assumingadistanceof1Gpcandcosmologicalparametersasdefined
inPlanckCollaborationetal.(2014).
upstream temperature has minimal effects on the Monte
UpS
Carlo results. Even at time step 8 the sonic Mach num-
region Regions shocked at
ber is approximately 170; all shocks modeled are strong
t
enough to justify our unusually high upstream tempera- i t
taouruerren.MoWnorneetleantoCitvaeirstltohicattirneaatltlhmepiearnrlttoiccolaeflssphslwoascemkptaaucfcpraelmbeyrea,tthaionendshdtohocaekst BTn000 dNdn(pip)i dNn( ip –i ) 1–i 1– 1 dNn( ipt –i ) 2– 2 n it –i 3– 3
nthoetFreeqrmuiireprpoacretsisc.les to be relativistic in order to enter BBii BBdi –p 1 Bddppi – 2 ddNNd((ppp))ii –– 33
Our evolutionary model for GRB afterglows selects i – 2 B
various points during the Blandford–McKee solution for i – 3
detailed modeling with the Monte Carlo code described
in Section 2. The specific times, and their associated
physical properties such as radius and shock speed, are
Shock at t
shown in Table 1. The Table also shows select Monte i
Precursor at t
Carlo parameters used to model the shocks. The only i
free parameter in Table 1 (i.e. not determined by either
Fig. 1.— Schematic showing the multi-shell model used in this
the Blandford–McKee solution or basic shock accelera-
work. At each time step, particle spectra and fluid parameters
tionprinciples)isNg. Convergencetestingoftheshocks aretrackedforboththepopulationofparticlescurrentlyinteract-
showed that larger values than presented in Table 1 do ing with the shock (green region and precursor at left of image).
notsignificantlyaffecttheresultspresentedinthispaper For each shell at each time, we record the density n, the particle
distribution function dN/dp, and the magnetic field strength B.
(e.g.,doublingN inducedachangeinthespectralindex
g Betweentimestepsallshockedpopulationsaresweptdownstream
p of a few hundredths). fromtheshockandevolvedasdescribedinthetext.
At each time step, we model the interaction of par-
is negligible.)
ticles with a shock whose Lorentz factor is given in Ta-
Aftereachtimestep,shockedplasmaissweptdown-
ble 1. We use two means to limit the maximum particle
stream as a single shell. We assume that the fluid in the
energy. First, we enforce a maximum acceleration time
shell has decoupled from the shock, and that no further
that any particle may have, which we set equal to the
acceleration will take place. The result is illustrated in
age of the shock. Second, we place a free escape bound-
Figure1. Atanygiventimesteptheshockisinteracting
ary (FEB) upstream from the shock. Any particle that
with the precursor and the closest downstream region.
reaches this boundary is allowed to escape the shock in
Shells that were shocked at previous time steps are ar-
the upstream direction, and is assumed to freely stream
ranged in an onion-like structure moving down the jet.
awayfromtheshockintothegalaxyatlarge. Wechoose
Thisapproachhasbeenusedtostudythenon-relativistic
2% of R as our FEB location; within this distance the
FS shocks of supernova remnants (e.g., Ellison et al. 2004;
physical curvature of the shock is ignorably small, and
Slane et al. 2014).
ourassumptionofaparallelplaneshockholds. (Inprac-
tice, the FEB location is irrelevant. No particles escape
3.2. Between time steps
upstream before the final two time steps, and even then
the fraction of energy/momentum flux they carry away As illustrated by Figure 1, we must keep track of
three quantities to model the evolution and downstream
NLDSA model for GRB afterglows 7
emissionofpreviously-shockedparticles: thetotaldensi- ence frames, however, and additional relativistic consid-
ty,theparticledistribution,andthestrengthofthelocal erations affect the flux observed at Earth. We elaborate
magnetic field. here, and note that each emission shell is handled inde-
The total density in any particular shell comes from pendently at each time step. Additionally, we assume
the Blandford–McKee solution. Both the upstream den- that all emitting populations are isotropic in their local
sity and forward shock position are known at all times, plasma frame; this will not be the case very near the
so (with the assumption that the jet is conical in shape) shock, but is true for a large majority of the particles
it is simple to count the total number of particles swept being considered.
up at each time step. We then integrate over the down- Both synchrotron and pion decay occur in the fluid
streamregionuntilthetotalnumberofparticlesmatches rest frame. For synchrotron radiation (Rybicki & Light-
the total swept up during that time step; the radius at man 1979), we assume that both the electron distribu-
which these two particle counts are equal separates ad- tionandthemagneticfieldturbulenceareisotropic. This
jacent shells. clearly is not the case upstream from the shock (Weibel
The particle distribution function has a fixed nor- turbulencecreatesanorderedfield,andcosmicraysnec-
malization, which assumes negligible particle transport essarily have a highly anisotropic distribution in order
between shells. The shape of dN/dp once the particles to have scattered ahead of a relativistic shock front); for
decouple from the shock is influenced by two processes: the downstream regions under consideration presently,
adiabatic cooling and radiative losses. As the plasma though,thisisareasonableassumption. Thekernelused
flows away from the shock and expands, on average the tocomputethepionproductioncrosssectioninp−pcol-
particlesarescatteringoffofmagneticturbulencethatis lisions (Kafexhiu et al. 2014) assumes that target parti-
recedingfromthem. Theaveragescatteringeventthere- clesareatrestinthelocalplasmaframe;foranassumed
fore results in a loss of energy, and over time a given isotropic distribution of cosmic ray protons this is true
particle’s energy varies as on average. The photon distributions produced by these
two processes may be transformed into the engine (i.e.,
dE E dV
=−(Γ−1) , (4) ISM)frameusingstandardspecialrelativisticequations.
dt V dt The case of inverse Compton emission is somewhat
whereE istheparticle’senergy, Γistheadiabaticindex moreinvolved, asmultiplephotonfieldsarepresent: the
of a monoenergetic population of particles at energy E, local synchrotron photons, the cosmic microwave back-
andV isthevolumeoftheshell(Reynolds1998). Inad- ground (CMB), and the interstellar radiation field (IS-
ditiontoadiabaticcooling,electronsradiateawayenergy RF). Synchotron self-Compton (SSC) emission happens
via synchrotron and inverse Compton losses; despite po- in the local plasma frame and is transformed alongside
tentiallyhighmagneticfieldsinthevicinityoftheshock, synchrotron and pion decay emission. The photon field
proton synchrotron is never energetically relevant. usedforSSCisthesynchrotronfieldofthelocalemission
As will be seen later (e.g., Figures 3 & 5), the shell: at a given location the number density of locally-
highest-energy electrons radiate away a significant frac- produced synchrotron photons is much greater than the
tion of their energy during the interval between time density of those produced in adjoining shells, so we ig-
steps. To accurately track their energies, we divide the nore photon transport between shells.
first time step after the particle population decouples InverseComptonusingtheCMBorISRFtakesplace
from the shock into a large number of substeps (only in the engine frame, so the electron distributions must
the first time step is treated in this manner; later cool- betransformedintothatframebeforephotonproduction
ing steps for the same population are handled in one canbecomputed. TheCMBisathermalphotonfieldat
interval). The duration of each substep is set so that a temperature of 3.35 K and energy density ≈ 0.60 eV
the highest-energy electrons lose no more than 10−3 of cm−3 (appropriate for a GRB at redshift z ≈ 0.23).
their energy during the substep. The intermediate elec- The ISRF is highly dependent on the environment of
tron spectra are stored for photon production, ensuring the progenitor star, but the work presented here uses
thatweareadequatelycapturingtherapidcoolingofthe the R = 8 kpc ISRF presented in Figure 1a of Porter
electrons in the post-shock environment. et al. (2008). This field is appropriate for a relatively
Finally,themagneticfieldstrengthisdeterminedby isolated progenitor star (i.e., not part of a compact clus-
assuming that the total magnetic energy in each shell is ter) in a Milky Way-like host galaxy. As will be shown
constant in time. The field strength therefore decays as later, however, the IC-ISRF luminosity is many orders
V−1/2. Note that this is a slower decay rate than that of magnitude lower than that due to SSC emission. The
exactshapeandnormalizationoftheISRFisunlikelyto
experienced by particles in the plasma, so the fraction
be observationally important.
of energy in magnetic fields increases with time (albeit
All of the preceding calculations are performed to
slowly).
arrive at a luminosity for the entire emitting shell. Ad-
3.3. Photon production ditional calculations are necessary to convert this lumi-
nosity into a flux as observed at Earth. When shifting
As mentioned in Section 2, we use three process-
from the plasma frame to the engine frame, relativistic
es to determine photon emission: synchrotron, inverse
aberration beams half of the emission into a cone whose
Compton, and pion decay due to proton-proton colli-
half-openingangleisgivenbycos(θ )=β ,whereβ is
sions. These three processes happen in different refer- h rel rel
8 Warren et al.
therelativevelocitybetweenthetwoframes(inthelimit higher than the typical value of 3 µG assumed for the
of β →1, this formula reduces to the more traditional Milky Way. There are locations where such high fields
rel
θ ∼ 1/γ ). This is also the half-opening angle for an mayexistinthecircumburstenvironment,14 butthisval-
h rel
observerlookingatthejet,so(atearlytimes)anobserv- uewaschosenforadifferentreason: onceitiscompressed
er does not see emission from the entire jet surface. We accordingtoEquation3,thedownstreamfieldapproach-
compute the fraction of the jet that lies within the an- es an energy fraction of (cid:15) ≈10−3, when both rest mass
B
gle θ and modify the total luminosity accordingly. We energy and pressure are considered. This is lower than
h
then compute the area of the θ cone at the appropriate thetraditionalvalue,(cid:15) ∼10−2,butinlinewithcurrent
h B
luminosity distance (1 Gpc) to convert from luminosity expectations (it is potentially even on the high side; see
to observed flux. Santana et al. 2014; Beniamini et al. 2015, and see also
Our photon spectra are additionally processed to the discussion of (cid:15) in (Lemoine et al. 2013)). Because
B−
model absorption at high energies, but not at low en- weassumeonlymagneticfieldcompressionandnoaddi-
ergies. Since the simulated afterglow occurs at cosmo- tional amplification, the value of (cid:15) just downstream of
B
logical distances, pair production off of the extragalac- the shock (that is, prior to the expansion described pre-
tic background light (EBL) absorbs energy flux at GeV viously) decays with time as the shock decelerates—see
or higher energies. We use the opacity tables present- Equation3—ratherthanmaintainingaconstantvalueat
ed in Franceschini et al. (2008) to calculate this absorp- all times. This temporal decay at the shock front occurs
tion. Atthelow-energyendofthespectrum,synchrotron in addition to the expected dilution of magnetic ener-
self-absorption (SSA) is the dominant absorption mech- gy density as the plasma expands downstream from the
anism. This process occurs locally within the shocked shock.
material,wherethevelocitydifferencefrombacktofront
may be great enough to induce significant anisotropies 4.1. Particle spectra and energy densities
in the absorbing electron population. We do not model
InFigure2weshowthespectraofcosmicraysforall
this more complicated absorption process here, but nei-
three models, in the local plasma frame. For the sake of
ther will we draw conclusions from the radio to infrared
clarity, thecooledpopulationsthathavebeensweptfur-
parts of our photon spectra.
ther downstream (the interior shells of Figure 1) are not
given. These spectra were taken just downstream from
4. EXAMPLEAFTERGLOWS thesubshockatx=0,tomostcloselymatchthespectral
Having explained the Monte Carlo code (Section 2) indicespredictedbyanalyticalwork. Thoughweassume
and the manner in which it is coupled to a hydrody- an isotropic distribution of particles for all photon pro-
namicalsimulationtocreateanevolvingafterglowmodel duction processes, this is obviously not the case at the
(Section 3), we will now present the results from three subshock. All particles require several diffusion lengths
testcases. Thefirsttwomodelsassumethatacceleration of transport downstream before their angular distribu-
takes place in the test-particle limit, with few enough tionrelaxesintoisotropy. Sincethediffusionlengthisan
particles accelerated that the structure of the shock is increasing function of momentum (Equation 1) higher-
not significantly modified. The particular values of P energy particles remain anisotropic for longer distances
inj
are listed in Table 1. In one model, thermal (uninject- downstream. Theeffectsofanisotropicparticledistribu-
ed) particles are explicitly excluded from the radiating tionswillbeaddressedinfuturework;fornowwemerely
population, though they contribute to the pressure and note the assumption being made.
energydensitydownstream. Werefertothisasthe“CR- The spectra presented in Figure 2 do not use arbi-
only” test case. In the second model, the “TP” model, trary units. Since upstream escape is negligible (no par-
all particles are included in the photon production cal- ticlesescapedduringthefirstsixtimesteps,andevenin
culations. Inthethirdmodel,“NL,”injectionisefficient the last two steps the escaping flux was less than 10−6
enough that the nonlinear backreaction of accelerated of the incoming energy flux), and since the difficulty of
CRs on the shock structure cannot be ignored. scattering upstream against a relativistic shock means
We must comment here on the physical reasonable- most particles are downstream, the total number of par-
nessofallthreemodels. TheCR-onlycaseisnot intend- ticles in each spectrum is the number of particles swept
ed to be physically plausible. There is no way to con- up at a given time step. As such, the spectra are not
serve number, momentum, and energy fluxes in a shock only absolutely normalized within each model: they are
interacting with the ISM unless a thermal population normalized between models as well.
is present. The CR-only case is included only because Both the CR-only case and TP model result in very
it mimics the standard synchrotron model for afterglow simple spectra. This is unsurprising, as they were de-
emission,whichignoresthethermalpopulation. Itserves signed to reproduce the simple power laws used in the
as a control against which our other two models may be standard synchrotron model for afterglows. In the TP
compared; we are not presenting it as a valid alternative
to either the TP or the NL model. Both the TP and NL 14 Specifically, Robishaw et al. (2008) report magnetic field
models are physically possible, and operate at different strengths upwards of 10 mG in starburst galaxies from measure-
ments of Zeeman splitting of the 1667 MHz OH megamaser line.
ends of the injection efficiency spectrum.
Additionally,ambientfieldstrengthswithin100-200pcofthecen-
Allthreemodelafterglowsuseanupstreamambient teroftheMilkyWaymaylocallyreachorexceedmGlevels(Fer-
magnetic field B = 3 mG. This is orders of magnitude ri`ere2009).
0
NLDSA model for GRB afterglows 9
model’s spectra, there is a thermal peak15 of particles
that crossed the shock once (but did not enter the ac-
celeration process), and a high energy tail of shock-
accelerated particles out to a maximum energy. The
thermal peak has been excised from the CR-only case,
while keeping the normalization of accelerated particles
the same. In all cases the spectral index of the proton
distributions decreases with time, as the shock Lorentz
factordrops(seeKeshet&Waxman2005;Ellison2005).
Since protons experience no losses while they in-
teract with the shock, the high-energy turnover of the
CR tail is due to the limit on acceleration time. The
electrons,bycontrast,experiencesignificantsynchrotron
losses in the intense magnetic fields downstream of the
shock. The slow increase in the maximum energy of the
electron spectrum is due to decreased compression also
lowering the magnetic field and cooling rate.
Our Monte Carlo model predicts that protons are
easily accelerated beyond the knee in the cosmic ray
spectrum;thehighmagneticfields,combinedwithpitch-
angle scattering, confine the protons to a small volume
around the shock, and allow for many shock crossing
cycles in the time-limited scenario presented here. The
proton spectrum extends above 1017 eV even at the ear-
liest time step shown in Figure 2, and the maximum en-
ergy slowly increases in successive time steps. This is a
consequenceofourmagneticfieldstructure,whichisuni-
formlystrongintheentirescatteringregiondownstream
from the shock. Both analytical studies and PIC sim-
ulations predict that the magnetic field strength should
decaywithdistancefromtheshock;thiswillincreasethe
scattering length of particles, allow them to spend more
time far from the shock, reduce the number of shock
crossingcyclestheycancomplete,andsolowerthemax-
imumenergy. Amoredetailedstudyofthismodel’spre-
dictions for ultra-high energy cosmic rays, with a more
realistic magnetic field structure, will be published else-
where.
Figure2alsoshowstheparticlespectraproducedby
the NL shocks. While much of the preceding discussion
(regarding,e.g.,spectralindicesandmaximumenergies)
applies to the NL model as well, there are a few impor-
tant points of comparison.
First, the thermal peaks of the two sets of spectra
occuratroughlythesamemomentaforagiventimestep.
The NL model’s peaks are always at slightly lower mo-
mentumthantheTPmodel’speaks;thisisawell-known
result that aids in reducing the downstream fluxes and
ensuring conservation. At the high-energy end, the TP
and NL models are almost identical. Both models’ elec-
tronspectraturnoveratthesameenergy(towithinone
or two bins, a factor of 1.6 or less), despite the spectra
extending over 4-7 decades.
The change from unmodified to NL shocks affects
the shape of the particle spectra. Where the unmodi-
fied shocks generate simple power laws above the ther-
Fig. 2.— Particle spectra for each time step, prior to any cool-
ing,fortheCR-onlycase(toptwopanels),TPmodel(middletwo 15 The term “thermal” is used only to specify those shocked
panels),andNLmodel(bottomtwopanels). Ineachpairofpan- particles that don’t enter the Fermi acceleration process. The re-
els,thetoppanelshowstheprotonspectraandthebottompanel sults we obtain in no way require or presume that these particles
shows the electron spectra. All momenta are in the local plasma equilibrateinathermodynamicsense.
frame,andcolorsrefertotimestepsasdescribedinTable1.
10 Warren et al.
i at time t. The value (cid:15) is defined analogously, with
p
(cid:112)
E(p)=m c2 1+(p/m c)2.
p p
The above definition of (cid:15) is an extension in two
e
important respects of the value commonly used in the
GRB literature. First, values of (cid:15) reported in the lit-
e
erature are assumed to be the fractional energy density
in electrons just after crossing the shock. In the one-
zone approach used in the standard synchrotron model,
the electron energy density is fixed to this value. In our
multi-zoneapproachwefollowthedownstreamevolution
of (cid:15) . The second point of interest about Equations 5
e
and6isthattheyexplicitlyincluderestmass-energy. In
the standard synchrotron model all particles are ultra-
relativistic, so contributions of rest mass-energy to the
energy density are minimal. As the shock slows, the
restmass-energyoftheprotonscannolongerbeignored
compared to their kinetic energy, and a decision must
be made whether to include it. We choose to include
rest mass-energy in our calculation of energy density,
and note that (cid:15) is also defined relative to the total,
B
not simply the kinetic, energy density. The inclusion of
rest-massenergyislargelyirrelevantatearlytimes,when
all particles are highly relativistic. At later times, when
Fig. 3.— Fractional energy density in the outermost shell as a
the thermal particles have plasma-frame Lorentz factors
functionoftime. TheTPmodelismarkedwitharedline,andthe
NL model with a black line. The colored points mark the values of just a few, including rest-mass energy enhances the
at time steps described in Table 1. Top panel: (cid:15)p, the fractional contribution from protons.
energydensityinprotons. Bottompanel: (cid:15)e,thefractionalenergy Thevaluesof(cid:15) and(cid:15) areplottedinFigure3forthe
densityinelectrons. e p
TPandNLmodels,fortheoutermostshell,notingagain
mal peak, the velocity gradients in the NL shocks in- that the simulation follows their evolution throughout
ducemorestructure. Thethermalpeaksarestillpresent the simulation. We do not distinguish between cosmic-
(albeit at slightly lower energies thanks to momentum ray and thermal particles for either model here. (The
and energy conservation), but there is now a substan- distinction between these two populations is obvious in
tial steepening of each spectrum before it levels out into the case of the TP model. As shown in Figure 2, the
somethingapproachingapowerlaw. Aboveacertainen- NL model’s thermal peaks are smoothly joined to the
ergy all particles have long enough mean free paths that accelerated population, making it unphysical to sharply
they can effectively ignore the velocity gradient of the distinguish where the former group ends and the latter
shockprecursorduringtheirscattering. Atthispointall begins.16) TheCR-onlycase,where(cid:15)e isnearlyconstant
bothmodels’spectrasettleintopowerlawswithidentical at ∼ 0.4 for the entire time modeled, is not drawn in
spectral indices (albeit different normalizations relative Figure 3.
to the thermal peak). The difference between the TP and NL models is
Bothparticlespectraandmagneticfieldstrengthare a few percent, but the TP model consistently places a
known everywhere in the jet at all times. This allows us larger fraction of the bulk kinetic energy into electrons
tocalculatedirectly,ratherthaninferfromobservations, thandoestheNLmodel. Thiseffectoccursinspiteofthe
two important parameters in the standard synchrotron fact that the electron-to-proton injection ratio is higher
model: (cid:15) and (cid:15) , the fractional energy density in elec- for the NL model than for the TP model. Finally, we
e B
trons and the magnetic field, respectively. Additionally, note that the effect of including the rest mass energy in
we calculate the time dependence these quantities have, Equation 6 is plainly visible in Figure 3: the value of
rather than fixing them to be constant throughout the (cid:15)e drops from nearly 0.4 at the start of the afterglow to
afterglow. Specifically, we define roughly 0.2 at the end.
(cid:15) ≡ Es(i,t) , (5) 4.2. Photon spectra
s Ee(i,t)+Ep(i,t) In Figure 4 we show photon spectra, all scaled to
a source distance of 1 Gpc, produced by the CR-only
where s could be either e (for electrons) or p (for pro-
case and the TP & NL models. The spectral shapes
tons), and
and normalizations vary substantially between models
(cid:90) dN (p,i,t)
Ee(i,t)= edp E(p)dp (6) 16 The Monte Carlo code, of course, allows us to separate CR
andthermalparticlesbasedonwhetherornottheyhavecompleted
(cid:112) at least one up–down–up cycle. While this is useful for tracking
is the total energy (i.e. E(p) = mec2 1+(p/mec)2), injectionrates,itisunlikelytobephysicallymeaningfulgivenhow
of all electrons (both thermal and accelerated) in cohort smoothlyconnectedthethermalandacceleratedpopulationsare.