AMRVAC and Relativistic Hydrodynamic simulations for GRB afterglow phases

0.9 MB·English

by  Z. Meliani

Mon.Not.R.Astron.Soc.000,??–??(2007) Printed3February2008 (MNLATEXstylefilev2.2)
AMRVAC and Relativistic Hydrodynamic simulations for GRB
afterglow phases
Zakaria Meliani1,2⋆, Rony Keppens1,3,4, Fabien Casse5 and Dimitrios Giannios2
1FOM-InstituteforPlasmaPhysicsRijnhuizenP.O.Box12073430BENieuwegein,Netherlands
2MaxPlanckInstituteforAstrophysics,Box1317,D-85741GarchingGermany
3CentreforPlasmaAstrophysics,K.U.Leuven,Belgium
4SterrenkundigInstituut,Utrecht,Netherlands
5AstroParticule&Cosmologie(APC)
7
UMR7164CNRS-Universite´Paris7,10rueAliceDomonetLe´onieDuquet75025ParisCedex13,France
0
0
2

n
Accepted;Received
a
J

5
ABSTRACT
1
We applya noveladaptivemesh refinementcode, AMRVAC, to numericallyinvestigatethe

variousevolutionaryphasesintheinteractionofarelativisticshellwithitssurroundingcold
1
Interstellar Medium (ISM). We do this for both 1D isotropic as well as full 2D jetlike fire-
v
4 ballmodels.ThisisrelevantforGammaRayBursts,andwedemonstratethat,thankstothe
3 AMRstrategy,weresolvetheinternalstructureoftheshockedshell-ISMmatter,whichwill
4 leaveitsimprintontheGRBafterglow.WedeterminethedecelerationfromaninitialLorentz
1 factorγ = 100uptothealmostNewtonianγ (2)phaseoftheflow.We presentaxisym-
0 metric2Dshellevolutions,withthe2Dextent∼chOaracterizedbytheirinitialopeningangle.In
7 suchjetlikeGRBmodels,wediscussthedifferenceswiththe1DisotropicGRBequivalents.
0
These are mainly due to thermally induced sideways expansions of both the shocked shell
/
h and shocked ISM regions.We found that the propagating2D ultrarelativistic shell does not
p accreteallthesurroundingmediumlocatedwithinitsinitialopeningangle.PartofthisISM
- matter gets pushed away laterally and forms a wide bow-shockconfigurationwith swirling
o
flow patterns trailing the thin shell. The resulting shell deceleration is quite different from
r
t thatfoundin isotropicGRB models.Aslongasthe lateralshellexpansionis merelydueto
s
ballisticspreadingoftheshell,isotropicand2Dmodelsagreeperfectly.Asthermallyinduced
a
: expansionseventuallylead to significantly higherlateral speeds, the 2D shell interacts with
v
comparablymoreISMmatteranddeceleratesearlierthanitsisotropiccounterpart.
i
X
Key words: GammaRays: Afterglow,Hydrodynamics,Theory– ISM:jets and outflows–
r
Galaxies:jets,ISM–methods:numerical,relativity,AMR
a
1 INTRODUCTION Marti&Mu˜ller2003).Theenormoustimeandlengthscaleranges
associatedwithviolentastrophysicalphenomenainrelativistichy-
Many high energy astrophysical phenomena involve relativistic
drodynamics(RHD),makeAdaptiveMeshRefinement(AMR)an
flows and shocks. For example, relativistic flows are invoked to
important algorithmic ingredient for computationally affordable
explain the observed properties of various compact astrophysical
simulations. RHD numerical simulations, particularlywhen com-
objects(Arons2004;Ferrari1998;Corbel2004).Astrophysicalrel-
binedwithAMRcapabilities,caninvestigatemanydetailsofrela-
ativisticflowscanreachaLorentzfactorof 2 10inassociation
− tivisticflowregimesrelevantforastrophysics.
withjetsfromSeyfertandradioloudgalaxies(Pineretal.2003),
orevengouptoLorentzfactors102 103 forGammaRayBurst
In this paper, we concentrate on relativistic dynamics in the
−
(GRB) scenarios (Sari&Piran 1999; Soderberg&Ramirez-Ruiz
fireball model for the afterglow phases of GRBs, in one and two
2001;Me´sza´ros2006).Inthelastdecade,continueddevelopment
dimensional simulations. Since the follow-up detection of GRBs
ofnumericalalgorithmsandtheincreaseincomputerpowerhave
in X-ray (Costaetal. 1997) and their afterglows at longer wave-
allowedtosignificantlyprogressinhigh-resolution,hydrodynamic
lengths (Sahuetal. 1997; VanParadijsetal. 1997; Galamaetal.
numerical simulations in both special and general relativity (see
1997; Frailetal. 1997; Piroetal. 1998), the cosmological origin
ofGRBshasbeenestablished(Metzger1997;Wijers1997).These
detections confirmed the predictions from the fireball theoretical
⋆ E-mail:[email protected] model (Rhoads 1993; Katz 1994; Me´sza´ros&Rees 1997; Vietri
c 2007RAS
(cid:13)
2 Z. Meliani,R. Keppens,F. Casseand D. Giannios
1997). In this model, a compact source releases a large amount theevolutionofanequivalent1Dsphericalshellisenlighteningin
of energy in a very short timescale, producing a fireball expand- thisrespect.
ing with relativistic velocity. Its internal energy gets fully con- Thispaperisorganisedasfollows.Westartbyreviewingthe
verted to kinetic energy, leading to a shell expanding with very relativistichydrodynamicequations.InSection3,weincludesev-
high Lorentz factor. This cold shell continues to expand and in- eralteststodemonstratetheAMRVACcodepotentialforrealistic
teractwiththecircumburstmedium,producingarelativisticshock- RHDcomputations.InSection4,wepresentourmainastrophysi-
dominatedevolution. Astheshell sweepsupthematter,itbegins calapplicationtoGRBflowsin1Dand2Dmodels.
todecelerate.Here,weinvestigatethedetailsofsuchpropagating
relativisticshells withthe relativistichydrodynamics code AMR-
VAC (Bergmansetal. 2004). TheAMRVAC code (Keppensetal.
2 RELATIVISTICHYDRODYNAMICEQUATIONS
2003)ishereforthefirsttimeappliedtothenumericallychalleng-
ingregimeofhighLorentzfactor,andwethereforeincludeava- The special relativistichydrodynamic evolution of a perfect fluid
rietyoftestproblems,demonstratingtherobustnessaswellasthe is governed by the conservation of the number of particles, and
limitationsofourcomputationalstrategy. energy-momentumconservation.Thesetwoconservationlawscan
bewrittenas
Up till recently, analysis of GRB flows have largely
been done analytically (Shemi&Piran 1990; Sari&Piran (ρuµ) =0, (Tµν) =0. (1)
µ µ
1995; Me´sza´ros&Rees 1997; Chiang&Dermer 1999),
whereρ,~u = γ,γ~v ,andTµν = ρhuµuν + pgµν define, respec-
combined with numerical approaches usually employing a (cid:0) (cid:1)
tively, the proper density, the four-velocity and the stress-energy
Lagrangian code. These latter works mainly investigate spher-
tensor of the perfect fluid. Their definition involves the Lorentz
ically symmetric GRB scenarios for obvious computational
factorγ,thefluidpressure p,andtherelativisticspecificenthalpy
convenience (Panaitescuetal. 1997; Kobayashietal. 1999;
h = 1+e+ p/ρ where e is the specific internal energy. For the
Kobayashi&Sari 2000). Recently, some analytical works
(inverse)metricgµν,wetaketheMinkowskimetric.Unitsaretaken
started to investigate the multidimensional jet structure in GRBs
wherethelightspeedequalsunity.
(Rhoads 1997, 1999; Panaitescu&Me´sza´ros 1999; Sarietal.
These equations can be writtenin conservative form involv-
1999; Kumar&Panaitescu 2000; Panaitescu&Kumar 2003;
ingtheCartesiancoordinateaxesandthetimeaxisofafixed‘lab’
Chengetal. 2001; Orenetal. 2004; Kumar&Granot 2003),
Lorentzianreferenceframeas
and some numerical simulations emerged as well, but restricted
to relatively low (order 25) Lorentz factor (Granotetal 2001; ∂U 3 ∂Fj
Cannizzoetal.2004).Higherspeedswereobtainedinthenumer- ∂t +X ∂xj =0. (2)
j=1
icalsimulationofthepropagationofanaxisymmetricjetthrough
a collapsing rotating massive star, as investigated by Aloyetal. Theconservedvariablescanbetakenas
(2000)toanalysethefirstphaseofGRBs.Inthesesimulations,the U = D=γρ,S~ =γ2ρh~v,τ=γ2ρh p γρ T , (3)
jetisfurtherfollowedafterbreakouttoamaximumLorentzfactor h − − i
ofγ 44,whichisstillrelativelysmalltothevaluesrequired andthefluxesarethengivenby
max ∼
for GRBs by the fireball model. Therefore, an important area of
currentinvestigationsinGRBcontextistomodelthedynamicsof F=hργ~v,γ2ρh~v~v+pI,γ2ρh~v γρ~viT , (4)
−
narrowjetsofultra-relativisticallyflowingejecta.Thisismotivated
whereIisthe3 3identitymatrix.Toclosethissystemofequa-
bytheneedtoreducethetotalamountofenergyreleasedinGRBs, ×
tions,weusetheequationofstate(EOS)foranidealgas,whichis
byassumingthesejetstopointtowardstheobserver,ascompared
thepolytropicequationwiththepolytropicindexΓ,
to fully isotropic equivalents. This need is particularly clear for
theexemplarycasesof GRB990123 (Kulkarnietal.1999),GRB p = (Γ 1) ρe. (5)
−
050820A(Cenkoetal.2006),andGRB050904(Frailetal.2006).
At each timestep inthenumerical integration, theprimitive
Thedetectionofpolarization(Covinoetal.1999;Wijersetal. variables (ρ,~v,p) involved in flux expressions should be derived
1999; Greineretal. 2003; Lazzatietal. 2004) gave further sup- fromtheconservativevariablesUresultinginasystemofnonlinear
porttothejetlikemodel.Evidencefornarrowcollimatedoutflows equations.Onecanbringthissystemintoasingleequationforthe
in GRBs is sustained also by the achromatic breaks in the after- pressurep,
glowlightcurveswhichwaspredictedanalytically(Rhoads(1997,
p+Γp(γ(p)2 1)
1999); Sari et al. (1999)) and then observed in a large number τ+D γ(p)D − =0, (6)
− − Γ 1
of GRBs (Staneketal. 1999; Sarietal. 1999; Bergeretal. 2000; −
Panaitescu 2005; 2006; Barthelmyetal. 2005). In Bloometal. which, once solved for pyields~v = S~ . Thisnonlinear equa-
τ+p+D
(2003),variousGRBswereanalysedandin16ofthem,thecom- tion(6)issolvedusingaNewton-Raphsonalgorithm.
bination of these breaks in the spectrum and the jet-like model
was used to deduce their effective energy, which was about E
∼
1051ergs.ThehalfopeningangleofsuchjetsinGRBsisinferred
3 TESTINGAMRVAC
to be of order few degrees. As a result, the afterglow produc-
ing shocked region is collimated too, with a similar initial open- In view of the challenges in the numerical investigation of rela-
ing angle (Frailetal 2001; Bergeretal. 2003; Cenkoetal. 2006; tivistic fluids, we include here several substantial test results for
Panaitescu 2005; Dar&DeRu´jula 2004). In our 2D simulations, codevalidation.Weperformedalargeseriesoftests,someofthem
wewillconcentrateontheafterglowphasesintheGRBevolution showninthissection.Animportantsubclassoftestcasesisformed
startingfromcollimatedejecta,anddiscussthosedynamicaleffects byRiemannproblems,whosenumericalsolutioncanbecompared
causing opening angle changes in detail. Direct comparison with toanalyticalsolutions.Other,2Dtestsshownherehavenoknown
c 2007RAS,MNRAS000,??–??
(cid:13)
AMRVACand RelativisticHydrodynamicsimulationsforGRB afterglowphases 3
Table1.L1errorsofthedensityforthe1DRiemannproblem1withuni-
formgridshownatt=0.36
Numberofgridpoints L1
200 1.15 10 1
400 6.4 ×10−2
800 3.2×10−2
1600 1.9×10−2
3200 1.06× 10−2
× −
nonvanishing tangential velocities, for two ideal gases withpoly-
tropic index Γ = 5/3. We separate two different constant states
p = 103, ρ = 1.0 (left) and p = 10 2, ρ = 1.0 (right). For
Figure1.One-dimensional relativistic shockprobleminplanargeometry L L R − R
att=0.36.Thesolidlinesaretheanalyticalsolution. thetransversevelocityweformninecombinationsofthepairvy,L
and v . Asin Ponsetal.(2000); Mignoneetal. (2005), we take
y,R
v =(0.0,0.9,0.99)cincombinationwithv =(0.0,0.9,0.99)c.
y,L y,R
Thespatialseparationbetweenthetwostatesisinitiallyatx=0.5.
analyticalsolution.Therefore,wecomparetheresultsofoursim- Theresultsatt=0.4areshowninFig. 2,wherewealsooverplot
ulationswithsimilarresultspreviouslyobtainedbyothercodesas theexactsolutionusingthecodeinMarti&Mu˜ller(2003).
documentedintheastrophysicalliterature. Therelativisticeffectsinthesetestsaremainlythermodynam-
TheAdaptiveMesh Refinement versionof theVersatileAd- icalinthefirstmildtest,andareduetocouplingbetweenthether-
vection Code (AMRVAC) is specifically designed for simulating modynamics(throughspecificenthalpy)andkineticproperties(by
dynamicsgovernedbyasystemof(near-)conservationlaws.Avail- theinitialtangentialvelocities).Forsmalltangentialvelocitycases,
able equations are the Euler and magnetohydrodynamic systems, weuseonlyaresolutionof200cellsonthebaseleveland4levels.
in both classical and special relativistic versions. The discretiza- However,forahightangentialvelocitycase,weusehighbaseres-
tionisfinitevolumebased,andvariousshock-capturingalgorithms olution400with10levelstoresolvethecontactdiscontinuityand
can be used. The automated AMR strategies implemented vary thetailoftherarefactionwave.Infact,forahightangentialveloc-
from the original patch-based to a novel hybrid block-based ap- ityatleft(inthehighpressurestate),theeffectiveinertiaoftheleft
proach(VanderHolst&Keppens2006).Theseproceduresgener- stateincreases. Thismakes theoccurring shock moveslowerand
ateor destruct hierarchicallynested gridswithsubsequently finer decreasesthedistancebetweenthetailoftherarefactionwaveand
meshspacing.TherefinementcriterionusedinAMRVACisbased the contact discontinuity. As also found in Zhang&MacFadyen
onaRichardson typeerrorestimator(Keppensetal.2003).Inall (2006),itremainsanumericalchallengetocapturethecontactdis-
followingtests,weusearefinementratioof2betweenconsecutive continuity properly, which we only managed here by allowing a
levels,unlessstatedotherwise. veryhigheffectiveresolution.
3.1 One-dimensionaltestproblems 3.1.2 ShockHeatingTest
3.1.1 Riemannproblems Inanother 1Dtestcase,acoldfluidhitsawallandashockfront
propagates back into the fluid, compressing and heating it as the
In1DRiemannproblems,wefollowtheevolutionofaninitialdis- kineticenergyisconvertedintointernalenergy.Behindtheshock,
continuity between two constant thermodynamical states. In 1D the fluid becomes at rest. This test has an analytical solution in
RHD, we then typically find the appearance of up to three non- planarsymmetryasconsideredbyBlandford&McKee(1976),and
linear waves. Generally, one findsashock wavepropagating into thejumpconditionsare
the lower density/pressure medium, a rarefaction wave propagat-
p = ρ (γ 1)(γ Γ+1) ,
ingatsoundspeedintothedensermedium,andbetweenthesetwo 2 1 1− 1
states,therecanbeacontactdiscontinuity.Intheteststhatfollow, ρ = ρ γ1Γ+1,
calculations are done in Cartesian geometry on a spatial domain 2 1 Γ 1
− γ v
06x61.TheexactsolutionsforRiemannproblemsinrelativistic v = (Γ 1) 1 1 . (7)
hydrodynamicsarediscussedforvanishingtangentialspeed(i.e.y sh − γ1+1
or/andzcomponentsforvelocities)inMarti&Mu˜ller(1994)and Thesegivethepostshockpressurep anddensityρ valuesinterms
2 2
forarbitrarytangentialflowvelocityinPonsetal.(2000). oftheincomingdensityandLorentzfactor,togetherwiththeshock
In a first, mild test, we assume an ideal gas with polytropic propagationvelocityv .
sh
index 5/3 and initial constant states characterized by p = 13.3, Inourtestwetakethesameinitialconditionsasintherecent
L
ρ = 10.0(left)and p = 0.66 10 6,ρ = 1.0(right),separated paperbyZhang&MacFadyen(2006),whereacoldfluidp=10 4
L R × − R −
atthelocationx=0.5.Theresultsatt =0.36areshowninFig.1 withadensityρ=1.0hasanimpactvelocityofv = 1.0 10 10 .
1 (cid:16) − − (cid:17)
witharesolutionof100cellsonthebaseleveland4levels,where This corresponds to a Lorentz factor γ = 70710.675. The tem-
wealsooverplot theexact solution.Inthetable1,wepresentthe perature after the shock becomes relativistic, and therefore we
L =Σ ∆x ρ ρ(x) normerrorsofthedensityρ,whereρ(x)is take the polytropic index Γ = 4/3. Hence the shock velocity is
1 (cid:16) j(cid:17)| j− j | j
theexactsolution.Theaccuracyofourresultiscomparabletothat v = 0.33332862. TheAMRsimulationisdonewith20cellson
sh
ofLucas-Serranoetal.(2004);Zhang&MacFadyen(2006). the base level and 4 levels on the spatial range 0 < x < 1. The
In a second test, we look particularly into effects due to resultatt=2,withthereflectivewallatx=1,isshowninFig.3.
c 2007RAS,MNRAS000,??–??
(cid:13)
4 Z. Meliani,R. Keppens,F. Casseand D. Giannios
Figure2.One-dimensionalrelativisticshockproblemsinplanargeometrywithtangentialvelocitiesvy.Theresultspresentedcorrespondtot=0.4.Thesolid
linesaretheanalyticalsolution(Ponsetal.2000);fromlefttorightvy,R=(0,0.9,0.99)c,andfromtoptobottomvy,L=(0,0.9,0.99)c.
The exact solution is overplotted as well. In this test, because of
v/c
theconstantstatebehindtheshock,themaximumimpactLorentz
factorthatcanbeachievednumericallyislimitedonlybythepre- P/8e9
cisionoftheNewton-Raphsonsubroutine.Thistestisimportantto
demonstrateitsaccuratetreatment,inviewoftheintendedsimula- ρ/5e5
tionsaimedatafterglowsinGRBs.Indeed,intheshell-frame,the
circumburstmediumhitsthedenseshellwithahighLorentzfactor.
Inaprocesssimilartowhatisfoundintheabovetest,thekinetic
energyoftheimpactingmediumisconvertedtothermalenergyof
theexternalmedium.Viewedinthelabframe,thesweptupcircum-
burstmediumwillhavesimilarlyhighLorentzfactorandwillform
ahot shocked layer ahead of thecontact interface. Notealsothat
Fig.3indicatesthatourdiscretizationandwalltreatmentdoesnot
sufferfromthevisibledensityerrorsseeninZhang&MacFadyen
Figure 3. One-dimensional shock heating problem in planar geometry,
(2006).
whereacoldfluidhitsawalllocatedatx=1.Theresultspresentedcorre-
spondtot=2.Thecomputationalgridconsistsof20zoneswith4levelsof
refinement.Thesolidlinesaretheanalyticalsolution.
3.2 Two-dimensionaltests
3.2.1 Arelativistic2DRiemannproblem
SW
ρ,v /c,v /c,p = (0.5,0.0,0.0,1.0),
(cid:16) x y (cid:17)
A twodimensional square region isdivided intofour equal areas
with a constant state each. We fix the polytropic index Γ = 5/3 ρ,v /c,v /c,p SE = (0.1,0.0,0.99,1.0). (8)
(cid:16) x y (cid:17)
andassumefreeoutflowboundaryconditions.Therelativisticver-
Thesimulation isdone with48 48cellsat thelowest grid
sionofthistestwasproposedbyDelZanna&Bucciantini(2002) ×
level,andweallowfor4levels.TheresultisshowninFig.4.Our
and subsequently reproduced by Lucas-Serranoetal. (2004);
resultisinqualitativeagreementwiththoseresultspublished,and
Zhang&MacFadyen (2006) and under slightly improved initial
shows the stationary contact discontinuities between SW-SE and
conditions by Mignoneetal. (2005). We repeat this simulation
SW-NWwith a jump in the transverse velocity. These are some-
withthesameinitialconfigurationfromDelZanna&Bucciantini
whatdiffusedbytheemployed TotalVariationDiminishingLax-
(2002),namely
Friedrichs(TVDLF)discretization(To´th&Odstrcˇil2006).Asim-
ρ,v /c,v /c,p NE = (0.1,0.0,0.0,0.01), ple and easily affordable remedy for improvement is to activate
(cid:16) x y (cid:17)
many more grid levels. Shocks feature across the interfaces NW-
NW
ρ,v /c,v /c,p = (0.1,0.99,0.0,1.0), NEand SE-NE,propagating diagonally totheNEregion, and an
(cid:16) x y (cid:17)
c 2007RAS,MNRAS000,??–??
(cid:13)
AMRVACand RelativisticHydrodynamicsimulationsforGRB afterglowphases 5
Figure4.Densitydistributionforthetwodimensionalshocktubeproblem
att=0.4.WithapolytropicindexΓ=5/3,abaseresolutionof48 48and
×
4AMRlevels.
elongateddiagonalshockstructureformsastheNEsectorrecedes
intotheRHStopcorner.IntheSWcorner,anobliquejet-likestruc-
tureformswithabowshock.
3.2.2 Relativisticjetin2Dcylindricalgeometry
Sinceit is relevant for our 2D GRB simulations, we also present
a two-dimensional simulation of an axisymmetric relativistic jet
propagating ina uniformmedium. Wesimulate theC2jet model
fromMartietal.(1997),butwithanenlargeddomainandathigher
effective resolution. Our computational domain covers the region
0 < r < 15 and 0 < z < 50 jet radii. Initially, the relativis-
tic jet beam occupies the region r 6 1,z 6 1, with v = 0.99, Figure5.Densitydistribution fortheaxisymmetric relativistic jetatt =
jet
130.Atleft,weshowthelabframedensity,atright,weshowtheproper
ρ = 0.01anditsclassicalMachnumber M =6.Inthiscase,the
jet densityinalogarithmicscale.Thecomputationalbasegridconsistsof90
jetissuper-sonicbutitstemperatureisstillclassical,sowecantake ×
300zoneswith5levelsofrefinementandthedomainsizeis15 50.
thepolytropicindexΓ = 5/3.Thedensityoftheexternalmedium ×
isρ = 1.0.Wefollowtheevolutionuntilt = 130,andthisend
ext
resultisshowninFig.5.Weperformedthesimulationwithareso-
lutionatthelowestlevelofthegridsetto90 300,andallowedfor
×
atotal of 5 levelsof refinement eventually achieving an effective
resolutionof1440 4800.
×
In this simulation, the relativistic motion of the flow domi-
nates,thethermalenergyisweakcomparedtothekineticenergy.
Asaresulttheexternalmediuminfluencesonlyweaklythejetand
the Lorentz factor γ 7 flow produces a cocoon structure from
∼
thetipofthejet.Onealsofindsaweaktransverseexpansionofthe
outflowinaccordwithwhatisreportedbyMartietal.(1997).This
transverseexpansionofthejetisinducedbythepressurebuild-up
insidethecocoonBegelman&Cioffi(1989).Inoursimulation,the
averagetransverseexpansionobtainedisoccurringatanestimated
speed of v = 0.11c. Moreover, at the contact interface between
T
shockedexternalmediumandjetmaterial,complexvorticalstruc-
tures form. These originate from Kelvin-Helmoltz type instabili- Figure6.Densitydistribution,inlogarithmicscale,fortheforwardfacing
ties,asaconsequenceofcoldfastjetoutflowmeetingamorestatic stepproblem,att=4.26.
medium.Theaveragepropagationspeedofthejetheadisfoundto
be0.414c,whichisinagreementwiththeone-dimensionalanalyt-
icalestimateof0.42casgivenbyMartietal.(1997).
c 2007RAS,MNRAS000,??–??
(cid:13)
6 Z. Meliani,R. Keppens,F. Casseand D. Giannios
3.2.3 Windtunnelwithstep 4.1 1Disotropicshellevolution
We reproduce here a standard test in the hydrodynamic litera- Inthissimulation,weconsideranISMwithuniformnumberden-
ture, namely the forward facing step test from Emery (1968); sityn = 1cm 3.ManyGRBafterglows(morethan25%)seem
ISM −
Woodward&Colella (1984), but adjusted to the relativistic hy- to be produced in such constant density medium (Chevalier&Li
droregimeasinLucas-Serranoetal.(2004);Zhang&MacFadyen 2000;Panaitescu&Kumar2002;Chevalieretal.2004).Thiscon-
(2006). A horizontal relativistic supersonic flow enters a tunnel stant density medium can be the resultant of a Wolf-Rayet star
withaflatforwardfacingstep.Thetestwasdonewitharesolution progenitor, with its surroundings shaped by a weak stellar wind
50 100zoneswith4levels.Thesizeofthetunnelis06x63and (VanMarleetal.2006).Initiallywesetauniformrelativisticshell
×
06y61.Thestepis0.2inheightanditspositionisatx=0.6.It at R = 1016cm from the central engine, since according to
0
istreatedasareflectingboundary.Theupperboundaryandlower Woods&Loeb(1995)theinteractionoftheshellwiththeISMbe-
boundary for x < 0.6 are also both reflecting. However, the left comesappreciableatthisdistance.TheshellhasaninitialLorentz
boundaryisfixedatthegiveninflowandtherightonehasfreeout- factor of γ = 100 (a γ > 100 is in accord with a shell which is
flow.Initially,thewholecomputationaldomainisfilledwithideal optically thin to gamma-rays (Woods&Loeb 1995; Sari&Piran
gaswith Γ = 7/5 withadensity ρ = 1.4moving at v = 0.999, 1995)),andenergy
x
i.e.,withaLorentzfactorγ = 22.37.TheNewtonianMachnum-
ber isset to 3. Theresult of our simulation isshown inFig. 6 at E=1054ergs=4πγ2R20δρshellc2, (9)
time t = 4.26. In this test, the relativistic flow collides with the
where δ stands for the lab-frame thickness of the shell set to
step, asaresult areverse shock propagates back against theflow
5 1012cm is of the order of the expected value for a fireball
directionandthisshockreflectsfromtheupperboundary.AMach ×
δ max(c∆t,R /γ2), where ∆t is the duration of the GRB. The
stemformsandremainsstationary.Theresultofoursimulationis ∼ 0
ISMandtheshellarecold,andtheinitialpressureissetto p =
comparabletowhatisreportedinLucas-Serranoetal.(2004). ISM
10 3n m c2 and p = 10 3n m c2 respectively. Notethat
− ISM p shell − shell p
thisimpliesahuge initialcontrast inthedensity measured inthe
lab-frame between the shell and the ISM D /D 109, and
shell ISM ∼
thispresents anextremechallenge fromacomputational point of
view.Initially,theenergyoftheshellisthenmainlykinetic.Weuse
4 GRBSANDMODELSFORTHEIRAFTERGLOW
aconstant polytropicindexΓ = 4/3,astheinteractionshell-ISM
PHASE
willbedominatedbytheforwardshock,wherethetemperatureof
ApopularmodelforGRBflowsisknownasthefireballmodel.In theshockedISMbecomesrelativistic.
thismodel,GRBsareproducedbyarelativisticoutflowfollowing Inthissimulation weuseaneffective resolution of 1536000
a violent event near a compact object. A large amount of energy cells corresponding to the highest grid level 10 allowed. We use
ispromptlyreleasedbythecompactsourceinaregionwithsmall thefullAMRcapabilitiesinthissimulation,sincewesimulateon
baryonloading(forareviewseePiran(2005)).Initially,mostofthe adomainofsize[0.3,300] 1016cm,with30000gridpointsonthe
×
energyoftheflowisintheformofinternal(thermal)energy.The lowest level. At t = 0, the shell itself isthen only resolved from
shellexpandsrapidlyconvertingitsinternalenergytokinetic.After gridlevel6onwards,whenweusearefinementratioof2between
theaccelerationphaseiscomplete,theshelliscoldandmoveswith consecutivelevels.Theinitialshellisresolvedbyabout25cellsin
relativisticspeeds. gridlevel10(laterinthedynamicalevolutionthismeansthatthere
This cold shell interacts with the circumburst medium, pro- aremanymoregridpointsthroughoutthewideningstructure).We
ducingstrongshocks.Oursimulationswillconsiderthedynamics usethisveryhigheffectiveresolutiontoavoidanynumericaldiffu-
from this phase onwards. As the shell sweeps up mass from the sionwhichmaycauseanartificialspreadingoftheshell.Weensure
externalmedium,thekineticenergyintherelativisticshellisgrad- thatthroughouttheentiresimulation,gridlevel10isactivatedand
uallytransferredtokineticandinternalenergyintheshockedam- concentratesfinegridsonboththeforwardshockandreverseshock
bientmedium.Moreover,theshellitselfgetstraversedbyareverse regions.Bothareveryimportanttodeterminetheprecisetimingof
shock,whichinturnconvertsthekineticenergyoftheshelltoin- thedeceleration.
ternalenergy. In Fig. 7, we show snapshots taken at lab-frame time t
≃
TheobservedafterglowemissionthatfollowsthepromptGRB 2.2 106scorrespondingtoanearlytimeintheentiresimulation,
×
emissionisbelieved tocomefromsynchrotron (withpossiblein- andinFig.9,weshowsnapshotstakenattimet 1.5 107scor-
≃ ×
verseComptoncontribution)emittingelectronsthatareaccelerated responding to a time when the shock is fully developed, we will
in the forward and reverse shocks (Sarietal. 1998; Galamaetal. concentrate our discussion mainly on this figure 9. These figures
1998).Intheinitialphases oftheshell-ISMinteraction,theelec- demonstrate that we resolve all four regions that characterise the
tronscanbeinthefastcoolingregime(i.e.theircoolingtimescale interaction between an outward moving relativistic shell and the
isshorterthantheexpansiontimescale)and,therefore,radiateeffi- cold ISM. From right to left, we recognize (Fig. 9) (1) the ISM
cientlymostoftheenergyinjectedtothem.Furthermore,ifmostof at rest, (2) the shocked ISM that has passed through the forward
theenergydissipatedintheshocksacceleratestheelectrons,then shock, with its Lorentz factor raised to γ(2) ∼ 30. This swept-up
onehastoconsiderradiativeshocks.Ifeitherofthepreviouscon- ISM gets compressed at the front shock and its number density
ditionsdoesnothold,theradiativelossesintheshocksaresmall. reaches n 75cm 3 Γγ(2)+1 (Sari&Piran 1995). These two
(2) ∼ − ∼ Γ 1
Here, we assume that the radiative losses are dynamically unim- valuescorrespondtotheanaly−ticalestimategivenbyeq.(7)forthe
portant,i,e.,theshocksareadiabaticthroughoutthesesimulations. front shock propagation. Region (3) represents part of the initial
Accordingtothemagnetizationoftheshell,theinteractionshell- shellmaterialwhichisshockedbythereverseshock.Thereverse
ISMandthe spectrum couldchange asisshown inMimicaetal. shockpropagatesbackintothecoldshell,reducingitsLorentzfac-
(2006). Here, we assume that the magnetic field is dynamically torandconvertingitskinetictothermalenergy.Transferofenergy
unimportant. fromtheinitialcoldshellthusoccursbothattheforwardandthere-
c 2007RAS,MNRAS000,??–??
(cid:13)
AMRVACand RelativisticHydrodynamicsimulationsforGRB afterglowphases 7
Figure7.ThefiveregionscharacterizingtheinteractionoftherelativisticshellwiththeISMatt=2.2 106s.Panel(a):Lorentzfactor,(b):logarithmofthe
×
pressure,(c):logarithmofdensity.
Figure8. Ratio ofthermal energy tomassenergy, when theshell reaches distance: R 6.5 1016cm (t 2.2 106s)(left), andR 4.3 1017cm
(t 1.5 107s)(right). ≃ × ≃ × ≃ ×
≃ ×
verseshock.Regions(2)and(3)areseparatedbyacontactdiscon- decreasesenoughduetosphericalexpansion,theLorentzfactorof
tinuity(Me´sza´ros&Rees1992).Atthissphericalcontactsurface, the reverse shock increases and the last part of the shocked shell
thelongitudinalvelocity(Lorentzfactorin1Dcase)andpressure becomeshote ρ .
(3)∼ (3)
remain constant, but there is a jump in density. Furthermore, re- Thereisanotherregion(5)indicatedinthefiguresbehindthe
gion(4)istheunshocked coldmaterialof theshell,moving with shell.Thedensityandthepressureintheregion(5)areverysmall
aLorentzfactorγ(4) = 100.Theweakthermalenergyoftheshell with n(5),min < 10−6cm−3 and p(5),min < 10−6mpc2. Therefore, the
interioritselfdidinduceaslightexpansioninthethicknessofthis region(5)isinthenumericalpointofviewavacuum.Thisregion
partoftheshell. (5) is not of strong interest for the physics of the afterglow, but
The reverse shock separating region (3) and (4) prop- itiscomputationally challenging toresolvetheinterfacebetween
agates into the cold shell with a Lorentz factor γ = the regions (4) and (5) where the ratio of the lab frame density
RS
γin(3e)ffiγ(c4i)e(cid:16)n1t−inv(r3a)vis(4in)/gc2t(cid:17)he∼th2er.5m.aTlheinserregvyercsoenstehnotckasisisNsehwotwonniainn pbreotwpaegeantitohneotwfothreeashchelels(DR(4)/1D0(156)c∼m1),0w14hiilnetthheefiexrsptapnhsiaosneoofftthhee
∼
Fig. 8 (left panel), where we draw the specific thermal energy shellremainsweak.BythetimeshowninFig.7,thiscontrasthas
in the shocked ISM and shell when the shell reaches a distance droppedtoavalueofatmost1012.
R 6.5 1016cm. Thereverseshock remainsNewtonianuntil it Thenear-totaldecelerationoftheshellonlytakesplacewhen
≃ ×
reachesadistancefromtheGRBsourceofR 3.8 1017cm.Then thetwoshocksintheshell-ISMinteractionmanagetoconvert an
∼ ×
it becomes mildly relativistic until R 4.3 1017cm where the important fraction of the kinetic energy of the shell to thermal
∼ ×
reverseshockbecomesveryefficienttoconvertthekineticenergy energy (andtheefficiency of thisconversion depends on whether
tothermalenergy(seeFig.8atright).Beyondthislatterdistance, thereverseshockisrelativisticorNewtonian,asdiscussedabove),
thedensityoftheunshockedshellpartρ hasdecreasedinaccord whiletherestistransferredtotheswept-upISMintheformofki-
(4)
withthesphericalexpansionoftheshell,toρ γ2 ρ .Asa neticandthermalenergy.Inthefirstphaseofthedeceleration,the
(4) ≪ (4) ISM
result,thereverseshockbecomesrelativistic.Thisbehaviorischar- maximumLorentzfactoroftheshelldecreasesgentlyfrom100ata
acteristicforaninitialthincoldrelativisticshelldeceleratingina distanceofR 2.5 1017cmto80atadistanceofR 4.3 1017cm.
∼ × ∼ ×
constantdensityexternalmedium.Infact,untiltheoutwardprop- However,onlyatthelatterdistanceof4.3 1017cm,asuddende-
×
agatingshell reachesR 3.8 1017cm,theshocked ISMmatter creaseof themaximumLorentzfactor of theentireconfiguration
∼ ×
ishote 30.0ρ (whereγ = 30correspondstotheanalyti- from γ = 80 to γ = 30 takes place. This fast drop of the maxi-
(2) ∼ (2) (2)
calsolutionfortherelativisticforwardshocke = γ 1 ρ ), mumLorentzfactorasseeninFig.10coincideswiththemoment
(2) (cid:0) (2)− (cid:1) (2)
whiletheshockedshellmaterialwhichhase =e iscold,since atwhichthereverseshockreachesthebackendofthecoldshell,
(3) (2)
e 0.01ρ .Whenthedensityinthenon-shockedshell(i.e.ρ ) thereby converting its kinetic to thermal energy (Fig. 9). In fact,
(3)∼ (3) (4)
c 2007RAS,MNRAS000,??–??
(cid:13)
8 Z. Meliani,R. Keppens,F. Casseand D. Giannios
Figure9.ThefivezonespresentwhentherelativisticshellinteractswiththeISMatt 1.5 107s.(a)Lorentzfactor,(b)logofpressure,(c)logofdensity.
≃ ×
110 Afterthisphase,theshellstructurecontinuestodecelerateby
γ
max
transferringitskineticenergytoshockedISMmatterattheforward
100
shock.However,asseeninFig.10,itstilltakesacertaintimebe-
90
forethevariation ofmaximum Lorentzfactor now characterizing
80 theshockedISMmatterfollowstheself-similaranalyticalsolution
70 forblast-wavedecelerationasputforwardbyBlandford&McKee
(1976).Fromaboutadistanceof1.2 1018cm,ournumericalsolu-
x 60
a ×
m tion starts to follow the analytical solution precisely. In fact, af-
γ 50
ter the reverse shock traversed the entire initial shell, a forward
40 travelingrarefactionwavepropagates through theentirestructure
30 thereby slowing it down while transferring most of the energy
to shocked ISM regions. This structure does not follow the self-
20
similarprescriptionandcausestheinitialdifference.Intheend,the
10
distancebetweentheforwardshockandthecontact discontinuity
0 0 50 100 150 200 250 300 increasedsufficientlyandtheresultingradialthermodynamicpro-
R [1016 cm] filesinbetweenbecomefullydescribedbytheBlandford-McKee
analyticalsolution.TheLorentzfactorpredictedbytheBlandford-
McKeesolutionbehindtheforwardshock(foranadiabaticshock)
16 Numerical simγu αlaRtio-3n/2 γ -1 α1 R6-2 isγBM=(E/ρISMc2R3)1/2 ∝R−3/2.ThepredictionoftheBlandford-
McKeesolutionfortheLorentzfactoroftheflowisalsoplottedin
BM phase near Sedov phase Fig.10andtheagreementwiththeresultsofthesimulationatthese
8 8 laterstagesofthedecelationisgood.
x Eventually,weenterintothemildlyrelativisticregimeforthe
a
m
blastwaveevolution.ThetransitiontotheSedov-Taylorphaseoc-
γ
4 4 cursbeyondthesimulateddistanceR>300 1016cm,sincewestill
×
haveaLorentzfactorofabout3attheendofthesimulation.The
Sedov-Taylor distance we find is close to the analytical estimate
2 2
givenbyl (3E/4πρ c2)1/3 5 1018cm.
∼ ISM ∼ ×
1 1
140 160 180 200 220 240 260 280 290 300
R [1016 cm] 4.2 2Dmodelingofdirectedejecta
Figure 10. The variation of the maximal Lorentz factor in the propa- Preciseanalysisoftheafterglowphasesrequirestoevolvenumer-
gating shell-ISM structure with time a) when the shell propagates from
ically confined ejecta in more than 1D, propagating in a jet-like
R0 = 1016cm toR = 300×R0 b)when the shell decelerates following fashion into the ISM. Wenow present axisymmetric, 2D simula-
theBlandford-McKeeprofile.
tionsof arelativisticcoldshell propagating inuniform ISMwith
anumber densityn = 1cm 3.Inthiswork,weinvestigatethe
ISM −
uniformmodeljet(Rhoads1999).Theshelldensityandenergyis
set constant throughout theshell, and wetake ittocorrespond to
whenplottingthemaximalLorentzfactorasafunctionofdistance, anisotropic spherical shell containing anequivalent isotropic en-
initially we always observe the Lorentz factor of the unshocked ergy E = 1051ergsand a Lorentzfactor γ = 100. Tomake the
iso
shell matter. As soon as the reverse shock has crossed the entire 2Dcomputationfeasible,wenowstartthesimulationwithashell
initialshell,westarttofollowtheevolutionoftheLorentzfactor thickness δ = 1014cm at a distance R = 1016cm from the cen-
0
of the shocked ISM (at the forward shock) where the maximum tralengine.Intheinitialsetuptheshelloccupiesanannularregion,
Lorentzfactoris30atthatparticularmoment. withhalfopeningangleoftheshellequaltoθ = 1 .Thisangleis
◦
c 2007RAS,MNRAS000,??–??
(cid:13)
AMRVACand RelativisticHydrodynamicsimulationsforGRB afterglowphases 9
rather small withrespect to those typically deduced from model- tothephasewhereE isrelevantinthedynamicstakesplacewhen
jet
ingoftheopticallight-curvebreaksbutstillinagreementwiththe theshockedISMandshellstarttospreadlaterallymuchfasterthan
mostcollimatedGRBflows(Bloometal.2003;Panaitescu2005). whatcorrespondstopureradial(ballistic)flow.
Withthechoiceofarathernarrowjet,weexpectthe2Deffectsto Inthelastpartoftheshell-ISMdecelerationphase,whenthe
appearearlierandtobemorepronouncedwithrespecttoaspher- reverseshock hascrossed theentireinitial shell material,the lat-
icalshellwiththesameisotropicequivalentenergy.Notethatasa eral velocityof shocked shell material reaches acomoving speed
result,thisjetlikeoutflowhasadecreasedeffectiveenergyinthe of v 0.7c. This means that we do find that the lateral ve-
shellE =(θ2/2)E . locitTy,coca∼n be bigger than the sound speed in the medium which
jet iso
TheAMRrunuses4gridlevels,taking400 6000atlevel1, isinaccordwiththeanalyticalresultofSarietal.(1999).Thisis
×
butwithrefinementratiosof2,4and2betweenconsecutivelevels at odds withnumerical findings asthose found inCannizzoetal.
eventuallyachievinganeffectiveresolutionof6400 96000.The (2004),whichemploy amuch reducedresolution ascompared to
×
domain is [0,4] [0.3,30] in R units, as we specifically intend our AMR results (at low resolution, we do obtain a reduced lat-
× 0
to model in detail the most dramatic phase of deceleration prior eral spreading velocity). As a result of this fast lateral spread of
totheBlanford-McKeeevolution.Notethatweinitiallyhaveγ > theshockedmaterial,distinctdifferencesoccurinthedeceleration
1/θ,whichappearstobethecaseforaGRBjet.Beamingeffects stagesascomparedtotheisotropiccase.Thisresultisveryimpor-
areinvoked toexplain theobserved steepening inthedecay light tant,asitshowsthatthelateralspreadingoftheshellisnotrelated
curveresultantfromthetransitionγ>1/θ(indistinguishablefrom onlytotheLorentzfactoroftheshellbuttothetypeofthereverse
anisotropicexplosion) toγ < 1/θ (Panaitescu&Me´sza´ros1999; shock. Inour computation, inanearlyphase thereverse shockis
Panaitescu&Kumar2002;Panaitescu2005).Withthesesetupwe Newtonianandtheexpansionoftheshockedshellpartismodest.
canverifyfromourhighresolutionsimulationwhetheruptotimes However,inalaterphasethereverseshockbecomesrelativisticand
correspondingtothetransitionγ 1/θ,onlytheisotropicenergy thisleadstofasterlateralspreads.However,astheforwardshock
∼
E isrelevantforthedynamicsandtheresultingemission(Piran is always relativistic, already in an early stage the shocked ISM
iso
2000;Granot2005). spreads with high velocity. The overall spreading of the shocked
Theinitialvelocityoftheshellispurelyradial.Notethat,com- ISMand shocked shell configuration can, thus, bequite complex
paredtothe1Disotropiccasepresentedintheprevioussection,the andrathermoreevolvedthantheonethatsemi-analyticalmodels
2D simulation starts with an initial condition containing less en- (Rhoads1999;Panaitescu&Me´sza´ros1999;Sarietal.1999;Piran
ergy.Thisisdoneformerepracticalreasons:wewishtokeepthe 2000)predict.
computationfeasiblewithintwoweek’sexecutiontimeonasingle OnlythatpartoftheISMfoundwithinthesolidangleofthe
processor. Duetothislower energy content, thedeceleration dis- expanding shell is swept up, and this opening angle changes due
tance willbe smaller by about one order of magnitude sinceless to the spreading effects just discussed. In our 2D simulation, we
swept-up ISM mass is sufficient to decelerate the shell. This re- foundinanalogywiththe1D(higherenergy)casefromabove,that
ducestheneedforresolvingmanydecadesofpropagationdistance thereverseshockisinitiallyNewtonian, sothethermalenergyin
as measured in units of the initial shell thickness. In fact the 1D theshockedpartoftheshelldoesnotincreasealotanditslateral
equivalentisotropiccasewiththesameenergyshowsasuddende- expansionremainsweakforawhile.Lateron,itslateralexpansion
creaseofthemaximumLorentzfactorofthedeceleratingconfigu- speed goes up to the 0.7c mentioned above, as the reverse shock
rationthatcorrespondstothereverseshockcrossingoftheshellat becomes relativistic and the material through which the shocked
R 9 1017cm.Thishappenswellbeforethesimulated3 1017cm. shellexpandslaterallyhasalreadybeenbroughttolowerdensities
∼ × ×
Fig. 11shows at the(top),thesound speed contour, and the bytheshockedISMinteraction.
density distribution in a logarithmic scale in the (left), and the The variation of the maximum Lorentz factor is less sud-
Lorentzfactor inthe(right).Atthebottom,azoomintheregion denthanintheequivalent1Dsphericalexplosion,asquantifiedin
around the shell is shown in the (left), the sound speed and the Fig.12.Thepartoftheshellmostdistantfromthesymmetryaxis
lateral velocity, in the (right) the Lorentz factor, and in the (cen- deceleratesbeforethemoreinternalpart.Theshellsweepsupmore
ter) a zoom only on the shell, the density and Lorentz factor. As matterthaninthecorrespondingisotropiccaseintheexternalparts
inthe1Dcase,atfirsttheshellpropagates withaconstant maxi- duetothelateralspreadingeffectsdiscussed.Infact,inthissimu-
mumLorentzfactor,andthisisaccompaniedbyaweakspreadof lationwemaydrawtheanalogybetweentheshellinteractionwith
the shell. Part of this radial shell widening in the bottom part of the ISM and simulations of relativistic AGN jet propagation into
theshellisaffectedbythecreationofaverylowpressureandden- an external medium. As in those cases, the energy is transferred
sityregionbelowtheshell(alsooccurringinthe1Dscenario).This totheISMthroughabowshockstructure. However,ourmodeled
near-vacuumstateremainsattherearpartastheshellmovesaway ejected shell representative of a burst in GRBsis not continually
atthespecifiedLorentzfactor.Inthis2Dsimulation,theunshocked supportedbyinjectionofenergyatthebottom.Asaresult,inthe
shellalsospreadslaterallywithaninitialtransverse(horizontal)ve- (smallopening angle) shell thereisnotreallyevidence of aclear
locity,sincetheshellislaunched withapureradialvelocity.The jet beam as inan AGN jet. In thiscase theinteraction shell-ISM
corresponding maximum initial lateral velocity of the unshocked isdominatedprimarilybyforward,reverseshockpair,andcontact
shellisv =0.0175c.However,theshocked,sweptupISMmatter discontinuityinbetween.Thechanging2Dstructureofthisshock
T
spreadslaterallymuchfaster,due toitshigh thermalenergy con- leadstodifferencesintheshockedISMmassloadedontotheshell.
tent.Initially,thatshockedISMpartspreadswithacomovingve- Asstatedearlier,onlyafractionoftheISMwithintheopeningan-
locityofv 0.4c,whichislessthanthemaximumsoundspeed gleof theshellisswept up,and animportant partof ISMmatter
T,co ∼
allowedbythepolytropicequationofstatec/√3.Duetothisfast gets pushed away laterally as the thin radially confined shell ad-
sidewaysexpansionofshockedshellandISM,themassoftheISM vances. The resulting behavior is clearly influenced by these 2D
hitbytheshellgrowsfasterthanr2.Therefore,thedecelerationof effects,andisthereasonwhythemaximumLorentzfactorofthe
theshellstartsearlierthanintheisotropiccase,seeFig12.Thisre- configurationstartstodecreaseearlierthaninthe1Dscenario.In
sultimpliesthatthetransitionfromthephasewhereE isrelevant, anisotropicscenarioallsweptupmassoftheISMremainsinfront
iso
c 2007RAS,MNRAS000,??–??
(cid:13)
10 Z. Meliani,R. Keppens,F.Casseand D.Giannios
Figure11.Onthe(top),thesoundspeed,logarithmofdensity(right),andLorentzfactorcontoursforthe2Dsimulation(left).Onthe(bottom),azoom
aroundontherelativisticshell,thesoundspeedandthelateralvelocity(left),theLorentzfactor(right),andthezoomontheshell,thedensityandLorentz
factor(center)att=5 106s.
×
c 2007RAS,MNRAS000,??–??
(cid:13)