Table Of ContentMon.Not.R.Astron.Soc.000,1–11(2010) Printed1February2008 (MNLATEXstylefilev1.4)
Relativistic blastwaves and synchrotron emission
⋆
T. P. Downes1 , P. Duffy2 and S.S. Komissarov3
1School of Mathematical Sciences, Dublin City University,Glasnevin, Dublin9, Ireland
2Department of Mathematical Physics, National University of Ireland, Dublin, Dublin4, Ireland
3Department of Applied Mathematics, Universityof Leeds, Leeds LS2 9JT, UnitedKingdom
Accepted —.Received—;inoriginalform—
2
0
0
ABSTRACT
2
n Relativistic shocks can accelerate particles, by the first order Fermi mechanism,
a which then emit synchrotron emission in the post shock gas. This process is of par-
J ticular interest in the models used for the afterglow of gamma ray bursts. In this
2 paper we use recent results in the theory of particle accelerationat highly relativistic
2 shocks to model the synchrotron emission in an evolving, inhomogeneous and highly
relativistic flow.We have developeda numericalcode which integratesthe relativistic
2
Euler equations for fluid dynamics with a general equation of state, together with a
v
simple transport equation for the accelerated particles. We present tests of this code
1
and, in addition, we use it to study the gamma ray burst afterglow predicted by the
7
fireball model, along with the hydrodynamics of a spherically symmetric relativistic
2
2 blastwave.
0 We find that, while, broadly speaking, the behaviour of the emission is similar to
1 thatalreadypredictedwithsemi-analyticapproaches,thedetailedbehaviourissome-
0 whatdifferent.The“breaks”inthesynchrotronspectrumbehavedifferentlywithtime,
/ and the spectrum above the final break is harder than had previously been expected.
h
p Theseeffectsareduetotheincorporationofthegeometryofthe(spherical)blastwave,
- alongwithrelativisticbeamingandadiabaticcoolingoftheenergeticparticlesleading
o to a mix, in the observed spectrum, between recently injected ”uncooled” particles
r
and the older ”cooled” population in different parts of the evolving, inhomogeneous
t
s flow.
a
: Key words: gamma rays: bursts – hydrodynamics – radiation mechanisms: non-
v
thermal – relativity
i
X
r
a
1 INTRODUCTION tering the particle distribution is almost isotropic both up-
stream and downstream. This results in a power law index
Relativistic shock fronts arise in astrophysics when a rel-
for the particle spectrum which only depends on the shock
ativistic flow propagates into ambient material such as in
compression ratio, r = U /U . In relativistic flows the sit-
thejetsofactivegalacticnuclei(AGNs)andfireballmodels 1 2
uation is more complicated. When v is comparable to U
of gamma ray bursts (GRBs). These shock fronts are also 1,2
thetimescaleforisotropisationoftheenergeticparticlepop-
sites of energetic particle acceleration which in turn leads
ulation by the magnetic turbulence (which determines the
to non-thermal emission processes such as synchrotron ra-
degree of anisotropy of the particle distribution upstream
diation and inverse Compton scattering. Particle accelera-
and downstream) becomes important. Moreover even if the
tion at relativistic shock fronts can occur via the first or-
particledistribution wereisotropic inonefluidframeit will
der Fermi mechanism. In this picture, energetic particles of
beanisotropicwhenviewedintheshockorotherfluidframe.
speed v are scattered back and forward across the shock
Asaresult,unlikethecasefornonrelativisticflows,thepar-
front by magnetic turbulence. Each crossing involves the
ticle spectrum depends not only on the shock strength but
particle meeting the flow head-on leading to an increase of
also on the details of how the particles are scattered (Kirk
energy. At nonrelativistic shocks, v >> U where U and
1,2 1 &Duffy1999).Thiscreatesdifficultieswhenstudyingshock
U are the upstream and downstream flow speeds in the
2 accelerationinflowssuchasspherical,relativisticblastwaves
shock rest frame respectively. With diffusive particle scat-
(Blandford & McKee 1976), the energetic particle distribu-
tion is the solution of a kinetic equation which describes
⋆ scattering in thelocal fluid rest frame.
AffiliatedtotheDublinInstitute forAdvancedStudies.
(cid:13)c 2010RAS
2 T. P. Downes et al.
However, there has been recent progress with ultra- inSect.5wegiveadiscussionofthehydrodynamicsofblast-
relativisticshocks,Γ>10,forwhichithasbeenshownthat waveswithnon-zeroinitialradius.Theresultsarepresented
the index of the power-law distribution of energetic parti- and discussed in Sect.6.
cles, which undergo isotropic scattering, is q = 2.23 where
N(E) ∝ E−q is the differential number of particles with
energy E (Kirk et al. 2000). This result has potentially im- 2 HYDRODYNAMICS
portant implications for fireball models of GRBs (Cavallo
Wewishtosimulateasphericallysymmetricexplosionwhich
& Rees 1978; Goodman 1986; Paczyn´ski 1986) in which a
givesrisetorelativisticvelocities.WeusetherelativisticEu-
large amount of energy is suddenly released in the form of
lersystemofequationsinordertocalculatetheevolutionof
an optically thick electron-positron plasma. The presence
thissystem.Thecodeusedisasecond-orderGodunov-type
of baryonic matter (Rees & M´esz´aros 1992) results in the
scheme (e.g. van Leer 1977) which uses a linear Riemann
conversion of the initial energy into a relativistically mov-
solver for shocks, and a non-linear one for strong rarefac-
ing bulk flow. The forward and reverse shocks will accel-
tions. Appendix A gives the details of the code, along with
erate particles by the above mechanism in the presence of
sometestcalculations.Inthissectionwebrieflydescribethe
magneticturbulence.Similarly,internalshockscanoccurbe-
equations used.
causeofvariabilityinthecentralengineofaGRBandthese
The conservation equations for inviscid relativistic hy-
arealsositesofrelativisticshockacceleration. Inthispaper
drodynamics in spherical symmetry are
we present a model for the synchrotron emission of shock
acceleratedparticlesdownstreamoftheforwardandreverse ∂ 1 ∂
(Γρ)+ r2Γρβ = 0 (1)
shocks in a spherically symmetric relativistic blastwave. In ∂t r2∂r
contrast to the strict fireball model our initial conditions ∂ wΓ2β + 1 ∂ r2 wΓ2β(cid:0)2+p (cid:1) = 2p (2)
are for a gas at rest with a large amount of internal ther- ∂t r2∂r r
mal energy confined to a small volume. As in the fireball (cid:0) ∂ w(cid:1)Γ2−p +(cid:2) 1(cid:0)∂ r2wΓ2β(cid:1)(cid:3) = 0 (3)
case the gas will accelerate to the point where all of the ∂t r2∂r
initialenergyisconvertedintobulkkineticenergyiftheex-
where Γ(cid:0)is the flu(cid:1)id Lorentz(cid:0)factor, ρ(cid:1) is the proper density,
ternal medium is sufficiently tenuous. Ultimately the flow
β is the velocity in c = 1 units, w is the enthalpy and p is
will become self-similar when theamount of swept upmass
the proper pressure. Time, t, and distance, r, refer to the
becomes significant and theflow will decelerate thereafter.
coordinatesmeasuredintheobserver’sframe.Wecanrelate
Ourprincipleaiminthispaperistostudytheinfluence
theenthalpy,density and pressure by
of an evolving, inhomogeneous flow on the energetic parti-
γ∗
cles produced by the first order Fermi mechanism at rela- w=ρ+ p (4)
γ∗−1
tivistic shocks. This is of central importance in the theory
ofGRBsandmanypapershaverecentlyapproachedrelated If the gas is composed of N species which each have mass
problemsonanumberofdifferentlevels.Theadiabaticevo- mi and number density ni then, for a Synge gas (Falle &
lution of a GRB interacting with an external medium, and Komissarov 1996),
theassociated emission, hasbeenconsidered (Panaitescu & N N −1
Kumar2000) whereananalyticmodelisintroducedforthe γ∗ = m n K3(ξi) −1 mini (5)
Lorentz factor of the shock front and relativistic remnant. γ∗−1 "i=1 i i(cid:18)K2(ξi) (cid:19)#"i=1 ξi #
Asimilar hydrodynamicalmodelhasalsobeen employedin X X
Moderski, Sikora & Bulik (2000) for thedeceleration of the where ξi = mkTi and Kα are the modified Bessel functions.
blastwavebutforbeamed ejecta.Suchageometry hasalso In this work we assume two species - electrons and protons
beendiscussedinarecentpaper(Granotetal.2001) where - which are in thermal equilibrium. The variable γ∗ defined
a full hydrodynamical simulation is carried out. In this pa- above is different to the ratio of specific heats, γ which is
perwesolvethefullhydrodynamicalequationsforaspheri- used in thecalculation of, for example, sound-speeds.With
callysymmetricsystemwhichstartsfromrestbutbeginsto this relation we should be able to solve the evolution equa-
expand relativistically, and ultimately becomes self-similar, tions 1 to 3. However, it should be noted that in relativis-
as a result of the enormous amount of internal energy in tichydrodynamicsthereare no explicit relations giving the
thecentralengine.Moreover,particlesareacceleratedatall primitive variables ρ, β and p in terms of the conserved
shockwaveswhichariseinoursimulation,whichinthiscase variables Γρ, wΓ2β and wΓ2−p. Therefore an iterative al-
are the forward and reverse shocks, with a spectrum deter- gorithmmustbedevelopedtodothis.Theprecisenatureof
minedbyaccelerationtheory.Thismodelwillbegeneralised this algorithm can affect how high a Lorentz factor can be
to jet typegeometries in futurework. simulated by the code. We have found that the code used
In order to achieve this we have developed a hydro- here will simulate flows up to Lorentz factors of at least
dynamical code, described in Sect. 2 and Appendix A, several hundredreliably.
which integratestherelativistic Eulerequationsreliably for
Lorentz factors of up to several hundred. This code is ap-
plicable to a gas with a general equation of state in one 3 PARTICLE ACCELERATION AND
dimension and incorporates adaptive hierarchical mesh re- SYNCHROTRON EMISSION
finement.InSect. A1we present tests of our hydrodynami-
3.1 The particle distribution
calcode.Ourmodelfortheinjectionofacceleratedparticles
into the flow, and their subsequent cooling, is presented in Relativisticshocks,inthepresenceofmagneticfluctuations
Sect.3.TheinitialconditionsarepresentedinSect.4,while which enable multiple shock crossings, can accelerate par-
(cid:13)c 2010RAS,MNRAS000,1–11
Relativistic blastwaves and synchrotron emission 3
ticles by the diffusive shock mechanism (Kirk & Schneider E )E−p (H being the Heviside function) which describes
min
1989). As described in the introduction it has been shown injectionofapowerlawspectrum,startingattimetˆ witha
0
that at ultra-relativistic shocks this produces a power law minimumparticleenergyE .Wecansolveequation(7)by
min
distributionofparticleswithindex2.23,i.e.N(E)∝E−2.23, finding the characteristic curves along which N(E,tˆ)dE =
so that theenergy densitydoes not divergeat high particle N(E ,tˆ)dE whereaparticlewithenergyE attˆ coolsto
0 0 0 0 0
energies. Theacceleration timescale inrelativisticflows will an energy E at time tˆ. From equation (6) we have
scale roughly with the particle mean free path in the tur-
bulent magnetic field divided by the shock speed; although
there is no published derivation of the precise acceleration d E−1ρ31 =αρ31B2 (8)
dtˆ
timescale. In this paper we are interested in what happens (cid:16) (cid:17)
tothepowerlawpopulation producedbyhighlyrelativistic whichcanbesolvedtoshowthatalongacharacteristiccurve
shocks over the much longer timescale which is the size of
theblastwavedividedbyc.Therefore,overhydrodynamical 1
ρ3E
timescales, we can treat acceleration as an impulsive injec- E = 0 (9)
tion of energetic particles, with thisuniversalpower law up 0 ρ13 −αE tˆρ13B2dtˆ
tˆ0
to arbitrarily high energies, into each fluid element which
passesthrougharelativistic shock.Wewill addresstherole N is conservedRso that
of an intrinsic cut-off in the shock accelerated population,
eitherasaresultofafiniteaccelerationtimescaleorlosspro- 1 1
cesses, in a future paper. The particles subsequently suffer dE0 = ρ03ρ3 (10)
synchrotron and adiabatic losses from which we can calcu- dE ρ31 −αE tˆρ13B2dtˆ 2
late the emission. We need to introduce three free param- tˆ0
(cid:16) (cid:17)
eters which, when taken together with the results of our and thesolution toR(7) becomes
hydrodynamicalsimulations, will giveacompletepictureof
theblastwave’sevolutionandtheassociated emission ofen-
ergetic particles. These parameters are N(E,tˆ)=Q ρ1−3pρ31E−p ρ31 −αE tˆρ13B2dtˆ p−2 (11)
0 0
• theratio,ǫb,betweenthemagneticfieldenergydensity Ztˆ0 !
and the thermal energy density,
In order to calculate the particle spectrum we need to
• thefraction,ǫ ,ofthedownstreamthermalenergyden-
e evaluatetheintegralin equation (11) foreach fluidelement
sity which is converted into high energy electrons as a fluid
whichhaspassedthroughashockatsomepoint.Nofluidel-
cell is shocked and,
ementpassesthroughtwoormoreshocksduringthesesimu-
• theenergy,E ,of thelowest energyparticle which is
min lations;ourmethodwouldhavetobeadaptedinsuchacase
produced at the shock. The differential energy spectrum is
(e.g. for overtaking internal shocks). We introduce a conti-
a power law in energy of index 2.23 from E .
min nuity equation for the integral quantity in equation (11).
Onceanelectronisinjectedattheshockitwillbescat- DefineI ≡ t0+t ρ31B2dtwherethefactorofΓ−1 boostsob-
teredinthelocalfluidframeandsufferbothsynchrotronand t0 Γ
server time to the comoving fluid time. Then the equation
adiabaticlosses. Ifthelengthscale overwhichtheelectrons R
for thedistribution of I is simply
areisotropisedismuchshorterthananyotherlengthscaleof
interest then therelativistic electrons will respond adiabat- ∂ (ΓρI)+ 1 ∂ r2ΓρIβ =B2ρ34 (12)
ically totheexpansionorcontractionoftheflow.Adiabatic ∂t r2∂r
losses, or indeed gains, are then described by the fact that Thus the quantit(cid:0)y I incr(cid:1)eases by ρ1/3B2∆t/Γ with each
p/ρ1/3 is constant where p and ρ are the particle momen-
time-step∆t,following thefluidelement,asrequired.Ifthe
tumand fluiddensityin thelocal fluid frame. With E =pc
gradientofΓβiseverlessthanacertainnegativevalue,cho-
forultra-relativisticparticlesthecombinedsynchrotronand
sen to ensure that the fluid element has passed through a
adiabatic losses are described in thecomoving frame by
shock,then I is set to0 at that point.The aboveconserva-
tion equation can easily be solved in conjunction with the
1ρ˙
E˙ =−αB2E2+ E (6) hydrodynamical equations using the same numerical meth-
3ρ
ods. Note that I is not kept at zero in unshocked fluid,but
where α is a constant and B the magnetic field strength in this does not affect the synchrotron emission since, in such
thelocal fluidframe. Consider nowa fluidelement which is fluid,thereare noenergetic particles.
shockedattimetˆ whenapowerlawdistributionofenergetic
0
particles is injected and where tˆis time measured in the
comoving frame. This population then evolves according to 3.2 Synchrotron emission
theequation
We approximate the emission of a single particle to be a
delta function in frequency with a single particle emissivity
∂N ∂
∂tˆ + ∂E E˙N =Q(E,tˆ) (7) in the fluid frame given by jνˆ(γ) = a0γ2B2δ(νˆ−a1γ2B)
where γ is the Lorentz factor of the electron in the fluid
where N(E(cid:0),tˆ) is(cid:1)the differential number of particles of en- frame, a =σ c/6π and a =e/2πm c with σ the Thom-
0 T 1 e T
ergy E at time tˆ. The losses, E˙, are given by equation (6) soncrosssection.Withthelocalelectronspectrumgivenby
while the injection term is Q(E,tˆ) = Q δ(tˆ− tˆ)H(E − 11, the emissivity in the local fluid frame is then
0 0
(cid:13)c 2010RAS,MNRAS000,1–11
4 T. P. Downes et al.
simplecasetothetreatmentusedinModerskietal.(2000).
In particular from those particles which have not radiated
away a significant fraction of their energy we get the un-
cooled spectrum of F ∝ ν−(p−1)/2 while the higher energy
ν
s particles have a steeper, cooled spectrum of F ∝ ν−p/2.
r ν
Below the minimum observed frequency the flux scales as
. θ φ . Fν ∝ ν1/3. At even lower frequencies the effect of syn-
z chrotronself-absorptionwillbecomeimportant,althoughwe
L have not explicitly calculated this, leading to a ν2 part of
thespectrum.
4 INITIAL CONDITIONS
Instudiesofthefireballmodelofgamma-raybursts,thepa-
Figure 1. Diagramofthetechniqueusedtocalculatethespec-
rametersusedarethetotalmass, M,andtheinitial radius,
trum of photons arriving at an observer at a given time t(obs).
The fluidelement is located at r inthe simulation, and the cur- R0. In general, E, the total energy, is thought to lie some-
renttimeist.This,inconjunctionwithL,thedistance between where in the region 1051-1054 ergs. The ratio between the
theobserverandtheblastwave, definesθ andφ.Seetext. rest-mass energy in the blast (i.e. the mass of the baryonic
component), and the total energy determines two critical
radii. These are the radius at which the baryons have been
Eˆνˆ= jνˆ(γ)N(γ,tˆ)dγ (13) accelerated up to their maximum velocity, Rc, and the ra-
Z diusatwhichtheshellofejectedbaryonshavesweptuptheir
Theemissivityintheobserver’sframeisrelatedtothat own mass in interstellar material, R . This latter radius is
d
inthefluidframebyEν =D2Eˆνˆ.TheDopplerfactor,D,for the radius at which the baryonic shell begins to decelerate.
thefluidelement with a threevelocity β makingan angle θ The former also depends on R , the initial radius of the
0
tothelineofsightisD=[Γ(1−µβ)]−1whereµ=cosθ.The blast. This initial radius is thought to be quite small, with
observedfrequencyisrelatedtothephotonfrequencyinthe R ∼1012 cm.
0
fluid frame by ν =Dνˆ. The observed intensity of radiation In order to get all the initial energy in the blast con-
ataninclination φtothelinebetweentheobserverandthe vertedintokinetic energy of thebaryons,it is also required
centreof theexplosion is thatthefireballbeoptically thick topaircreation untilthe
∞ radiusoftheblastexceedsR .Ifthiswerenotthecase,then
c
Iν(φ,t0)= Eν(s,φ,t)ds (14) thephotonswouldescapeintospacebeforeaccelerating the
Z0 baryons to their maximum velocity. In our case we put the
where t0 is the time of observation and t = t0−s/c is the initialenergyoftheblastintothermalpressure.Thismeans
time of emission in the observer’s frame a distance s away thatitisalreadyinthekineticenergyofthebaryons.Hence
from the observer (figure 1). The observed flux density is the latter consideration need not be taken into account in
obtained byintegrating theintensity over theentire source ourinitial conditions.
Theresultsoftwosimulationsarepresentedinthispa-
Fν(t0)= Iν(φ,t0)dΩ (15) per. Each of thesimulations used thefollowing properties:
ZΩ
• E =1051 ergs
In cylindrical coordinates (z,r) centred on the explosion,
• Ratio of energy to mass: η = E = 580. This value
s≈L−z since φ≪1 so that Mc2
being chosen so that the Newtonian value of R is twice
d
2π rc rc R . In reality, R is slightly less than R due to relativistic
F (t )= dz E (z,r,t)rdr (16) c d c
ν 0 L2 ν effects.
Z−rc Z0
• R =1.2×1014 cm.Thisismuchlargerthanthepostu-
with L the distance from the observer to the centre of the 0
latedphysicalvalue.However,reducingR to1012cmwould
explosionandr thesizeofthecomputationaldomain.Ifr 0
c 0 impose extremely severe computational costs on the simu-
is the initial radius of the explosion then the interval of t
0 lation. This increase in the value of R is not expected to
since the arrival of the first photon is t = t −(L−r )/c 0
0 0 0 changethepropertiesofthesynchrotronemissioncalculated
so that the flux can be written as F (t ). By the change
ν o as we are interested in the afterglow in this work, and not
of variable z = r +c(t−t ) the integration over z can be
0 0 theemission from theacceleration phase.
convertedtooneovert,thetimeofemissionintheobserver’s
• ǫ =0.01-thisvalueisconsistentwithobservedvalues
frame, e
found by,e.g. Wijers & Galama (1999)
2πc tend rc • ǫ =0.01 in one simulation, and ǫ =0.1 in the other.
F (t )= dt E (r +c(t−t ),r,t)rdr (17) b b
ν 0 L2 ν 0 0 This is the only way in which thetwo simulations differ.
Z0 Z0
In the simple case where the fluid is assumed to move We must also choose the injection energy E for the
min
directly towards the observer with a constant velocity we shockaccelerationprocess.Physicallythiswillbesetbynon-
have successfully compared our method with the exact re- linear plasma processes. The highest energy thermal parti-
sults.Theaboveanalysis,whichisessentiallythesameasin clesdownstreamofashockcanleakintotheupstreamform-
Komissarov & Falle (1997), also corresponds for the latter ingabeamofparticleswhichcanexciteplasmaturbulence.
(cid:13)c 2010RAS,MNRAS000,1–11
Relativistic blastwaves and synchrotron emission 5
This turbulencethen scatters particles back into the down-
stream fluid and multiple shock crossings become possible. 10
4-velocity
Theproblemofdeterminingthedetailsofthisprocessisan log(proper density)
unsolved one in the theory of particle acceleration at rela- 8
tivisticshocks.However,itisclearthatthevalueE must Reverse
min
correspond to the high energy tail of the downstream ther- 6 Shock
Ingoing
malpopulation.Thisishardlysurprisingsincethethickness
Rarefaction
of a shock wave will be set by the gyroradius of a thermal 4
particle in a collisionless plasma and first order Fermi ac-
Forward
celeration willonlybeapplicabletoparticleswithgyroradii 2
Shock
muchlargerthanthis.Consequentlywehavefixedthevalue
of Emin = 10mec2 which is chosen so that it exceeds the 0
thermal energy for all shocked fluid elements. Setting E Contact
min
Discontinuity
as a factor of a few times thethermal energy at theinstant -2
a piece of fluid is shocked does not change our results ap- 0 0.5 1 1.5 2 2.5 3 3.5 4
preciably. r
Initially,then,thereisasphereofradiusR whichisat
0
restattheorigin.Thisspherecontainsalltheinitialenergy
10
and mass of the blast. Outside this sphere is a gas with a 4-velocity
density of 1 cm−3. 8 log(proper density)
Ingoing
6 Reflected Rarefaction
5 HYDRODYNAMIC EVOLUTION OF A Rarefaction
4
SPHERICAL BLASTWAVE
Reverse Shock
2
Inthissectionwediscussthequalitativeevolutionofaspher-
ical blastwave with finite initial radius. This discussion ap- Forward Shock
0
pliestobothrelativisticandnon-relativisticcases.Thepur-
Contact Discontinuity
pose is to elucidate, in a qualitative fashion, how the for- -2
ward and the reverse shocks are formed, and how they be-
have,from ahydrodynamicperspective.Figure2showsthe -4
0 0.5 1 1.5 2 2.5 3 3.5 4
early-timeevolutionofthesystem,computedusingthecode
r
presented here. The initial conditions used were as follows:
Figure2. Plotsofthespatialcomponentofthe4-velocity,along
3 on r≤1 with the proper density for times of t = 1.5 (top) and t = 3
ρ =
1 on r>1 (bottom). Allthewaves inthesystem arelabelled.Seetext.
(cid:26)
u = 0
103 on r≤1
p = 3×10−3 on r>1 fast as it would in the planar case, and the reverse shock is
(cid:26) formed to decelerate material behind the forward shock to
These initial conditions are similar to those in Sect. theappropriate speed.
A1.3,andarechosensothattherelevantfeaturesareclearly When the rarefaction reaches r = 0 it is reflected and
visible in theresults. enhanced.This“reflectedrarefaction”canbeseenclearlyat
Initiallyashockisdrivenintotheambientmedium(the t = 3 in figure 2. This rarefaction will eventually catch up
“forward shock”infigure2),whileararefaction propagates with the reverse shock/blastwave system (see figure 4 and
from r = R back to r = 0 (the “ingoing rarefaction” in Sect.6.1).Whentherarefactioncatchesupwiththereverse
0
figure 2). Note that, at t=1.5, the ingoing rarefaction has shock,theram-pressureofthematerialenteringthisshockis
stillnotreachedtheorigin.Theforwardshockgetsprogres- drastically reduced, because it has been rarefied and decel-
sively stronger as blast material gets accelerated down the erated.Henceitcannolongersupportthethermalpressure
pressuregradientoftherarefaction.Thereisalsoanentropy between the forward and reverse shocks. As a result, the
wave(the“contactdiscontinuity”)whichmovesoutfromR material between these two shocks expands.Since theram-
0
with a speed slightly lower than that of the blastwave. pressure of thematerial entering theforward shock has not
In addition, a reverse shock is created almost imme- changed (assuming constant external density), this neces-
diately between the entropy wave and the rarefaction (see sitates the reverse shock slowing down with respect to the
figure 2). If the system were in planar symmetry then one rest-frame of the initial blast. Indeed, it will actually begin
would not expect this reverse shock since there are only propagating towards theorigin. The speed at which it does
three characteristics in hydrodynamics(see, e.g. Landau & this depends on the initial value of η, since this determines
Lifschitz 1966). We get this fourth wave purely as a result thestrengthof theinitial(and sothereflected)rarefaction.
ofthenon-planarnatureofthesystem.Itformsbecause,as For example, if η is very small, then the force due to the
the forward shock moves outwards it sweeps up an increas- pressure gradient will produce a relatively small increase in
ing amount of interstellar matter per unit distance. This velocity of thematerial in theinitial sphere.Hencetherar-
means that the forward shock does not move outwards as efaction will be weak. If, on the other hand, η is large then
(cid:13)c 2010RAS,MNRAS000,1–11
6 T. P. Downes et al.
the force due to the pressure gradient will accelerate the
blast material to high velocities very quickly, producing a 100000
very strong rarefaction. For the conditions in Sect. 4 it be-
comes sufficiently strong to cause ejected material to flow 3) 10000
^
m
back towards theorigin also. c
When thereverse shockreaches r=0 it reboundsand, es/ 1000
cl
asitpropagatesoutwardsagain,itweakens.Oncethisshock arti
is sufficiently weak to be ignored, we have reached the self- y (p 100
similar stage of evolution. Prior to this point the finite ra- sit
n
dius of the initial blast plays a part, and hence there is a e 10
D
significant length-scale in the system. If the blast was ini- er
p
tially ofzeroradiusthenthewholeprocess describedabove Pro 1
wouldhappeninfinitelyquickly(or,equivalently,would not
happen at all). 0.1
Wecanrelatetheabove“hydrodynamical”picturewith 0 50000 100000 150000 200000
the one used by, for example, Rees & M´esz´aros (1992), as Distance (light-seconds)
follows. The coasting radius, R , is reached when all ma-
c 70
terial has been accelerated down the pressure gradient of
the ingoing rarefaction. This, of course, never occurs, but 60
one can define R as being the radius when “most” of the
c
material has been accelerated. c) 50
of
the rTeflheecdteedcerlearraetfaiocntiorandciuast,chReds,uispthweinthththeerardeviuesrsaetswhhoicckh. units 40
This is when the material between the forward and reverse y (
shock expands back towards the origin, necessarily deceler- ocit 30
el
ating significantly as it does so. 4-v 20
10
6 RESULTS
0
0 50000 100000 150000 200000
Wediscusstheresultsofthesimulationsintwosections.The
Distance (light-seconds)
first deals with the detailed hydrodynamic evolution of the
blastwave. The second concerns the observed spectra and 100
light-curves.
10
2) 1
6.1 Hydrodynamic results ^
m
c 0.1
Hereweshowthehydrodynamicaspectoftheresultsofour es/
n
simulations. Figure 3 shows the distribution of proper den- dy 0.01
sity,velocityandpressureafter1×105 seconds.Wecansee e (
ur 0.001
the two rarefactions mentioned in Sect. 5. The head of the ss
e
reflectedrarefactionliesatthelocationofthepeakvelocity, Pr 0.0001
with its tail at r = 0. The original rarefaction lies between
1e-05
thepeakofvelocityandtheriseinthedensityandpressure
plots at ∼1×105 light-seconds. 1e-06
0 50000 100000 150000 200000
The forward shock can be identified as the right-most
Distance (light-seconds)
rise in pressure and density, while the reverse shock is just
to theleft of this.
Figure 3. Plotsofdensity, 4-velocityandpressure(top tobot-
Figure 4 shows the same plots as figure 3, but for a tom)foratimeof1×105 secondsaftertheinitialblast.Wecan
time of 7×106 seconds. Here we can see the reverse shock see the original and reflected rarefaction, as well as the forward
has begun to separate significantly from theforward shock. andreverseshocks.Seetext.
Thereflectedrarefactionisnowtheonlyrarefaction present
in the system. The reverse shock is expanding back into
this rarefaction due to the insufficient ram-pressure of the the influence of these errors on the overall solution, or on
material inthereflectedrarefaction tosupportthepressure thecalculatedsynchrotronemission,isnegligibleduetothe
of the material between the forward and reverse shocks. low density of theregion containing theerrors.
It is possible to see oscillations in the density plot just Figure 5 shows the latter stages of the non-self-similar
behindthereverseshock.Theseoscillations ariseout ofnu- evolution of thesystem.Atthis stage thereverseshock has
mericalerrorsduringthetimethatthereflectedrarefaction propagatedallthewaytor=0andhasrebounded.Wecan
comescloseto,butisnotincontactwith,thereverseshock. seethatithasweakenedgreatly,andwillcontinuetodoso.
The rarefaction is extremely strong here, and it is not sur- Once this shock is sufficiently weak to be ignored the sys-
prising that numerical errors become significant. However, temwillevolveinaself-similarway.Wecanseetheentropy
(cid:13)c 2010RAS,MNRAS000,1–11
Relativistic blastwaves and synchrotron emission 7
100 100
10
3) 1 3) 10
^ ^
m m
s/c 0.1 s/c
e e
cl cl
arti 0.01 arti 1
p p
y ( 0.001 y (
sit sit
n n
e 0.0001 e 0.1
D D
1e-05
1e-06 0.01
0 2e+06 4e+06 6e+06 8e+06 1e+07 0 2e+06 4e+06 6e+06 8e+06 1e+071.2e+071.4e+07
Distance (light-seconds) Distance (light-seconds)
12 6
10 5
8
c) c) 4
nits of 6 nits of 3
u 4 u
city ( 2 city ( 2
o o
Vel Vel 1
0
-2 0
-4 -1
0 2e+06 4e+06 6e+06 8e+06 1e+07 0 2e+06 4e+06 6e+06 8e+06 1e+071.2e+071.4e+07
Distance (light-seconds) Distance (light-seconds)
1 100
0.01 10
2) 2)
^ ^
m 0.0001 m 1
c c
s/ s/
e e
n n
y 1e-06 y 0.1
d d
e ( e (
ur ur
s 1e-08 s 0.01
s s
e e
Pr Pr
1e-10 0.001
1e-12 0.0001
0 2e+06 4e+06 6e+06 8e+06 1e+07 0 2e+06 4e+06 6e+06 8e+06 1e+071.2e+071.4e+07
Distance (light-seconds) Distance (light-seconds)
Figure4. Asforfigure3,butforatimeof7×106 secondsafter Figure 5. As for figure 3, but for a time of 1.4×107 seconds
theinitialblast. aftertheinitialblast.
6.2 Spectra and light curves
errors generated bythevery strongreflected rarefaction re-
maining just ahead of the reverse shock now. These errors, In this section we present the results from the calculations
although apparently significant here, do not affect the syn- of thesynchrotron emission.
chrotron results dueto the low densities in this region, and
the low relativistic boosting suffered by emission from this
material. Both these effects make the influence of emission
6.2.1 Spectra
from this material negligible.
Figure 6 shows a plot of the spectrum observed 24 hours
after the initial blast could have been observed for thecase
whereǫ =0.1(hereafterreferredtoasthehighǫ case)and
b b
(cid:13)c 2010RAS,MNRAS000,1–11
8 T. P. Downes et al.
case gives a slightly steeper spectrum at these frequencies
1e+10 than the low ǫ case. For high magnetic fields the contri-
High B-fields b
1e+09 Low B-fields bution from cooled electrons will be stronger closer to the
shock than in the low magnetic field case. Therefore, while
1e+08
we still expect the uncooled electrons to be preferentially
y) 1e+07 boosted in the spectrum, the effect of the cooled electrons
J
micro- 1e+06 will bIteisstwroonrgtherdfiosrcuhsisgihnegrtmheagbneehtaicvifiouelrdos.fthe‘critical’ fre-
x ( 100000 quencies in the spectra with time - i.e. where the breaks
u
Fl 10000 occur. The lower break, which arises due to thepresence of
electrons injected with the minimum energy E behaves
1000 min
roughly as t−0.73. This is the same for both the low and
100 high ǫ cases. This is a much less dramatic decrease with
b
10 timethanthatpredicted bySari& Piran (1997) andis due
100001e+061e+081e+101e+121e+141e+161e+181e+20 tothefactthat,here,wedefinetheminimumenergyofinjec-
Frequency (Hz) tionapriori.Thisenergyischosensothatthegyroradiusof
anelectronatthisenergywouldbeacceleratedbytheFermi
Figure6. Thespectrumobservedat24hrsaftertheinitialblast
mechanism.Thelatter authorschoose todefinethenumber
observationforthecaseofǫb=0.01(“LowB-fields”)andǫb=0.1 of electrons accelerated, and this, in conjunction with the
(“HighB-fields”).
energy in energetic particles, defines E . This produces
min
the different behaviours of the lower break in the spectra
for ǫ =0.01 (hereafter referred to as the low ǫ case). It is with time.
b b
clearthatinbothcasesthespectrumisabrokenpower-law. Wefindthatthefrequencyofthesecond(upper)break
In the first section of the spectra, the flux goes as goes like t−1.2 for the low ǫb case and t−1.75 for the high ǫb
ν1/3, as expected for synchrotron radiation below the peak case.Thisisamuchfasterdecreasewithtimethanthatpre-
frequency of the lowest energy electron. At lower frequen- dictedinSari&Piran(1997).Thisdifferencecannotbedue
cies,whicharenotplottedhere,synchrotronself-absorption to the different choice of definition of Emin, but it could be
would give a flux of ν1/3. The second part of the spectrum hypothesised that it is due to the effects of adiabatic cool-
is, in both cases, the power law ν−0.615 indicating that this ing. This cooling will cause thesecond critical frequency to
part of the spectrum is dominated by emission from elec- decrease faster where the flow is divergent. However, simu-
trons,withadistributionofE−2.23,whichhavenotsuffered lations of this system without adiabatic cooling show that,
significant adiabatic or synchrotron cooling. whileadiabaticcoolingdoeshavethiseffectontheobserved
The break from this part of the spectrum to the fi- spectrum, it is too small to explain the discrepancy. It is
nal, steeper, part occurs in different places in the spectra notclearwhatspecificallycausesthisdifferencebetweenour
plotted. For the high ǫ case, the break occurs at a lower resultsand those of Sari& Piran (1997), butit mustbere-
b
frequency than for the low ǫ case as would be expected membered that here we perform a much fuller treatment of
b
since, in higher magnetic fields, the losses suffered by the the dynamics and so it is not surprising that we get some
energetic population are correspondingly higher. However, significant differences with semi-analytic approaches.
while in the final part of the spectrum the exponent pre-
dictedfromsimpletheorywouldbe−p/2=−1.115,wefind
the spectrum to be slightly harder. In the high ǫb case, we 6.2.2 Light curves
have F ∝ ν−1.019, and in the low ǫ case, F ∝ ν−0.981.
ν b ν
The reason for the slightly harder spectra at high frequen- Inthissectionwediscussthelightcurvesresultingfromthe
ciesisthenon-uniformvelocitydistributioninthe“shell”of simulationsdescribedabove.Figure7showsplotsofthelight
ejectedmaterial(see,e.g.,figure4).Materialmovingathigh curvesatlow,mediumandhighfrequencies.Thesefrequen-
velocitytowardstheobserverwillbemoreheavilyweighted cies are chosen so that the points are always in the regime
in the spectra than material moving at lower velocity, due whereFν ∝ν1/3 (lowfrequencies),Fν ∝ν1−2p (mediumfre-
to relativistic beaming. Since such high velocity material quencies) and, for high frequencies, the frequency is always
occursimmediately behindtheforward shock,thismaterial in the third section of thespectrum.
will contribute more to the spectrum than material further It should be noted in the following discussion that the
back from the shock. The fluid just behind the shock con- range of times covered by these light curves is 1 to 24 hrs
tains an electron population which has only recently been after the first signal of the blast reaches the observer. Ini-
accelerated andsoemission from herewillbedominatedby tially we see an increase in all the light curves shown, with
uncooled electrons to a very high frequency. Material from theexception of thehigh frequency,high ǫ case, where the
b
further behind the shock will have emission dominated by emission is roughly constant over the first couple of hours.
cooled electrons down to a lower frequency. Hence, above a Thisincreaseismuchfasterthanthatpredictedbyprevious
certain cut-off, we expect to get a mixture between emis- analytic work. For example, for the low ǫ case at low fre-
b
sion from cooled and un-cooled electrons, with un-cooled quencies,thefluxgoesasroughlyt1.4 atearlytimes(cfSari
electronsbeingpreferentially weighted.Thisleadstoanex- & Piran, 1997, where the maximum increase is t1/2). This
ponent for the spectrum lying between the uncooled value effect is not due to any problems with numerical resolution
(1−p = −0.615) and the cooled value (−p = −1.115). This since the emission in the 1st hour is dominated by mate-
2 2
conclusion is further supported by thefact that thehigh ǫ rial between the forward and reverse shocks approximately
b
(cid:13)c 2010RAS,MNRAS000,1–11
Relativistic blastwaves and synchrotron emission 9
The medium frequency light curves initially increase,
1e+09 asalreadystated,andthenfalloffdramatically.Again,itis
High B-fields
Low B-fields difficulttoseeanypointatwhichthebehaviourofthelight
curve is a true power-law. The high frequency light curves
alsobegintofalloffverysteeply,andagain,itisdifficultto
y) see any sign of a power-law behaviourin thecurves.
J
o- The lack of a clear power-law behaviour in the light
micr 1e+08 curves may well be due to the restricted time-scale over
x ( whichthecurvesarecalculated(from1hourto1day).How-
u
Fl ever,itisclearthat,ifthereisabrokenpower-lawbehaviour
thenthebreaks aresmeared out, and,in addition, thefinal
fall-offofthelight-curvewouldseemtobemuchfasterthan
previously predicted.
1e+07 The unexpected behaviour of the low and medium fre-
1000 10000 100000 quency light curves is more likely to do with the different
FreqTuiemnec y(s )(Hz) treatmentoftheminimumenergyintheacceleratedelectron
population.Inourcase,sincewefixE a priori,weallow
1e+09 min
High B-fields the available energy to define the number of energetic elec-
Low B-fields
trons.Previously it has beencommon todefinethenumber
ofenergeticelectronsandhencetheavailableenergydefines
E .Thiswillleadtoquitedifferentbehavioursofthelower
y) 1e+08 min
J break frequency with time.
o-
cr
mi
x (
u
Fl 1e+07 7 CONCLUSIONS
We have presented a model for the synchrotron emission
of energetic particles downstream of relativistic, spherical
1e+06 shock waves. The hydrodynamical part of the problem has
1000 10000 100000 beensolvednumericallywiththesimplersimulationsagree-
FreqTuiemnec y(s )(Hz) ingwithresultspublishedelsewhere.Ontheotherhandthe
particleacceleration aspecthasbeentreatedasaninjection
1000
High B-fields processwithrelativisticshocksleavingapopulationofener-
Low B-fields geticparticlesimmediatelydownstream.Thespectralindex
isknownfromsemi-analyticworksothatweneedonlyspec-
ify the fraction of the downstream thermal energy which is
Jy) converted into energetic particles and a lower cut-off en-
o-
micr 100 ecrogoyli.nTghaenpdaratdicilaebsastuicbsleoqsuseesn.tlTyhloesehyenderorgdyynbaymsiycnaclhrreosturlotns
ux ( havecapturedtheevolutionofboththeforwardandreverse
Fl
shocks, which are of principal interest for particle accelera-
tionandradiativeemission,aswellastherarefactionwaves.
Thespectraemittedfromoursystemlargelyagreewiththe
simple predictions for the low energy, uncooled part of the
10
spectrum.However,at thehigherfrequency,cooled, partof
1000 10000 100000
theelectronpopulationtherelativisticboostingofthemate-
FreqTueimncey ( s()Hz)
rial coming straight towards theobserverhardensthespec-
Figure7. Lightcurvescalculatedatlow,intermediate,andhigh trumfromthepurecooledvaluewhichcomesfrommaterial
frequencies.Seetext. furtherdownstream and from material which isnot moving
directlytowardstheobserver.Further,non-trivialbehaviour
hasbeenfoundforthetemporalvariationofboththebreak
6×106 seconds after the initial blast, when the two shocks frequencies and the light curves. Adiabatic losses coupled
are well-resolved and separated by thecode. with integrating the spectrum over a spherical system tend
Thebehaviourofthelowfrequencylightcurvesisqual- tosmearoutthebreakfrequencieswhicharepredictedfrom
itatively different to the other two. It keeps increasing for simple scaling arguments.
the duration of the simulation. This is unsurprising as it It is clear that we need to include several other effects
will only begin todecrease when thecritical frequencyν before bringing our results into contact with observations.
min
corresponding to Γ passes through the frequency of the Ageneralisation toa2-D hydrodynamicalcode hasalready
min
light curve, and we have chosen this frequency so that this begun.Theroleofinternalshockswillalsobeincludedalong
does not happen. It is very difficult to say that the light withtheireffectonacceleration.Theenergeticparticlesobey
curve behaves anything like a power-law with breaks over a transport equation, the solution of which determines the
this time-scale based on our results. spectralindexandweareworkingtowardsincludingsucha
(cid:13)c 2010RAS,MNRAS000,1–11
10 T. P. Downes et al.
solutioninourmodel.Morecomplicatedradiativeprocesses −ri2−1/2F~in−+11//22 +S~in+1/2∆t (A1)
can also be included such as the inverse-Compton and syn-
i
chrotron self-Compton processes. However, what we have where superscripts refer to the time index and subscripts
done in this paper is to combine a realistic hydrodynami- refer to thespatial index.Also,
calmodelwithresultsfromparticleacceleration theoryina U~n = Γρ,Γρβ,wΓ2−p T n (A2)
relativistic flow. The computed hydrodynamical evolution, i
i
swpiethctroabsaenrdvaltigiohntsciunrvfuestuwriellwaollrokw. ustocomparethismodel F~in±1/2 = h(cid:0)Γρβ,wΓ2β2+p,w(cid:1)Γi2β T N (A3)
i±1/2
h(cid:0) 2p T n (cid:1) i
S~n = 0, ,0 (A4)
i r
ACKNOWLEDGMENTS (cid:20)(cid:16) (cid:17) (cid:21)i
WewouldliketothankLukeDruryforusefuldiscussionson ItshouldbenotedthatthesourcetermS~ shouldbevolume-
theincorporationofadiabaticlossesintotheelectrondistri- averagedover[ri−1/2,ri+1/2](Falle&Komissarov1996).So,
bution. We are grateful to an anonymous referee for com- for spherical symmetry, and using a first order scheme we
ments which led to improvements in the paper. This work getthatthesourceterm fortheradialmomentumequation
was supported by the TMR programme of the European is
Union undercontract FMRX-CT98-0168. Sn =3pn ri+1/2+ri−1/2 (A5)
i i ri2+1/2+ri+1/2ri−1/2+ri2−1/2
For a second order scheme, which allows p to vary linearly
REFERENCES in a given cell the correct form for Sn is
i
BlandfordR.D.,McKeeC.F.1976,Phys.Fluids,19,1130 Sn = 3pn ri+1/2+ri−1/2
CDaovwanlleosGT..,PR.,eResayMT.J.P.,.,1917989,8,MAN&RAA,S3,3118,31,133509 i i ri2+1/2+ri+1/2ri−1/2+ri2−1/2
FalleS.A.E.G.,1991, MNRAS,250,581 +g n 9(ri+1/2+ri−1/2)2(ri2+1/2+ri2−1/2) (A6)
FGaololedmS.aAn.EJ..G,.1,9K86o,mAispsJa,ro3v08S,.SL.4,71996,MNRAS,278,586 pi (cid:20) 4(ri2+1/2+ri+1/2ri−1/2+ri2−1/2)2 (cid:21)
Granot J., Miller M., Piran T., Suen W.M., Hughes P.A., 2001, wheregpni isthegradientofthepressureincelliattime-step
astro-ph/010038 n.
KhokhlovA.M.,1998,J.Comp.Phys. 143,519 ThefluxesF~ arecalculatedfromtheprimitivevariables
KirkJ.G.,DuffyP.,1999, J.Phys.G,25,R163 extrapolated to the cell edges using non-linear averaging of
KirkJ.G.,SchlickeiserR.,Schneider P.,1988,ApJ,328,269
thegradients(e.g.vanLeer1977;Downes&Ray1998),mak-
KirkJ.G.,SchneiderP.,1989,A&A,225,559
ing the code second order in space. A non-linear Riemann
KirkJ.G.,Guthmann A.W.,GallantY.A.,Achterberg A.,2000,
solver is used in the presence of strong rarefactions, and
ApJ,542,235
otherwise a linear Riemann solver is used in order to save
KobayashiS.,PiranT.,SariR.,1999,ApJ,513,669
KomissarovS.S.,FalleS.A.E.G.,1997,MNRAS,288,833 computational time. The flux terms F~in±+11//22 are calculated
Landau L.D., Lifschitz E.M.,1966, Fluid Mechanics. Pergamon, using a first order Godunov-scheme from the conditions at
Oxford time n. This makes thecode second order in time.
M´esza´rosP.,LagunaP.,ReesM.J.,1993, ApJ,415,181 ThedetailsofthelinearandnonlinearRiemannsolvers
ModerskiR.,SikoraM.,BulikT.,2000,ApJ,529,151 follow closely the discussions of Falle (1991) and Falle &
Paczyn´ski B.,ApJ,308,L43 Komissarov (1996).
PanaitescuA.,KumarP.,2000,ApJ,543,66
The code also contains an option for hierarchical grid
PiranT.,ShemiA.,NarayanR.,1993,MNRAS,263,861
refinement. The algorithm used for this is loosely based on
ReesM.J.,M´esza´rosP.,1992,MNRAS,258,41P
thatofKhokhlov(1998).However,thisoptionisnotusedin
SariR.,1997,ApJ,489,L37
thework presentedhere,sowewill not discussthisfurther.
SariR.,PiranT.,1995, ApJ,455,L143
SariR.,PiranT.,1997, ApJ,485,270
vanLeerB.,1977,J.Comp.Phys.,23,276
A1 Tests of the hydrodynamical code
Wijers,R.A.M.J.,Galama,T.J.,1999,ApJ,523,177
The code used in this work has been rigorously tested
against, in particular, results presented in Falle & Komis-
APPENDIX A: NUMERICAL DETAILS sarov (1996). The code reproduces their published results.
Below we show three examples of these tests for complete-
In this appendix we discuss the details of thecode used for
ness.
thecalculationspresentedinthispaper.Wealsopresentthe
results of threetest calculations.
Weemploy asecond orderfinitevolumeGodunov-type A1.1 Test 1: Shock-tube problem with γ = 4
3
scheme to solve the equations 1 to 3 in spherical symme-
In thistest we use thefollowing initial conditions:
try. Assuming that the grid cells are defined so that cell i
occupies the space [ri−1/2,ri+1/2] then the scheme can be ρ = 1
written as u = 0
U~in+1 =U~in− r3 3−∆tr3 ri2+1/2F~in++11//22 p = 1100−32 oonn xx≤>00..55
i+1/2 i−1/2 h (cid:26)
(cid:13)c 2010RAS,MNRAS000,1–11
