Dynamic Monte Carlo Simulations of Radiatively Accelerated GRB Fireballs
aa r X i v : . [ a s t r o - ph . H E ] F e b MNRAS , 1–15 (2015) Preprint 5 November 2018 Compiled using MNRAS L A TEX style file v3.0
Dynamic Monte Carlo Simulations of RadiativelyAccelerated GRB Fireballs
Atul Chhotray, ⋆ Davide Lazzati Department of Physics, Weniger Hall, Oregon State University, Corvallis 97331, OR, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present a novel Dynamic Monte Carlo code (DynaMo code) which self-consistentlysimulates the Compton scattering driven dynamic evolution of a plasma. We usethe DynaMo code to investigate the time–dependent expansion and acceleration ofdissipationless GRB fireballs by varying their initial opacities and baryonic content.We study the opacity and energy density evolution of an initially optically thick,radiation–dominated fireball across its entire phase space - in particular during the R ph < R sat regime. Our results reveal new phases of fireball evolution: a transitionphase with a radial extent of several orders of magnitude - the fireball transitionsfrom Γ ∝ R to Γ ∝ R , a post–photopsheric acceleration phase - where fireballsaccelerate beyond the photosphere, and a Thomson–dominated acceleration phase- characterized by slow acceleration of optically thick, matter–dominated fireballsdue to Thomson scattering. We quantify the new phases by providing analyticalexpressions of Lorentz factor evolution, which will be useful for deriving jet parameters. Key words: stars: gamma-ray burst: general – radiation: dynamics – radiation mech-anisms: thermal – radiative transfer – scattering
Gamma–Ray Bursts (GRBs henceforth) are intense burstsof high energy radiation (high energy X-rays and γ rays).First detected serendipitously about five decades ago(Klebesadel et al. 1973), today GRBs are regularly detectedby space-based satellites and are known to be one of thebrightest explosions in the universe ( L ∼ − ergs / s)(see Piran 2004). These large luminosities along with smallvariability time-scales and the emission of high energygamma rays ( ∼ + keV) led to the compactness problemin GRBs. The compactness problem was resolved by in-voking ultra-relativistic motion (Paczynski 1986; Goodman1986) of the emitting source, which has been confirmed byobservations (Frail et al. 1997; Frail et al. 2001)GRBs are thought to be powered by the core-collapseof a massive rotating star (leading to Long GRBs) andthe merger of two neutron stars / a neutron star andblack hole (resulting in Short GRBs) (see Kumar & Zhang2015). Although ultra-relativistic jets are invoked to explainGRBs, the mechanism responsible for the jet productionfrom GRB progenitors is not well understood and is underinvestigation. The two mechanisms that have been proposed ⋆ E-mail: [email protected] (AC) to launch and accelerate jets to relativistic speeds are 1)magnetic fields and 2) radiation. Several prior works havestudied the driving role of magnetic fields in collimatingand accelerating relativistic outflows in astrophysicalenvironments (see for e.g., Blandford & Znajek 1977,McKinney 2006, Komissarov 2011, Tchekhovskoy et al.2011). Consequently, they have been been proposed as thejet production mechanism in GRBs (Drenkhahn & Spruit2002; Lyutikov et al. 2003). In this paper we will explore indetail radiation as an alternative mechanism to accelerate(and possibly launch) relativistic outflows.Radiative acceleration of outflows is a well known as-trophysical phenomenon (e.g., continuum and line–drivenstellar winds – see Castor et al. 1975; Smith & Owocki2006). To produce relativistic outflows associated withGRBs and AGNs, radiative acceleration due to external ra-diation sources has been investigated. Madau & Thompson(2000) studied the radiative acceleration of cold, op-tically thin plasmas due to external radiation sources.Zampieri et al. (2003) studied radiative acceleration of lowdensity ion-electron plasma by incident transient radiationalong with radiation-induced internal electric fields. Theabove mentioned works have primarily focused on thedynamical evolution of a low-density, optically thin plasmaunder the influence of radiation sources external to the © Chhotray & Lazzati
Figure 1.
Diagram illustrating the phase space evolution of aGRB fireball based on the fireball model. The y axis depicts theratio of radiation to matter energy density and the x-axis plotsthe opacity. In the phase space diagram, GRB fireballs begin theirevolution in quadrant 1 ( τ ≫ , e γ / e m ≫ ) and end it in quadrant3 ( τ ≪ , e γ / e m ≪ ). The arrows show possible evolutionarypaths of fireballs. The shaded regions near τ ∼ and e γ / e m ∼ represent transition zones or regions where fireball evolution hasnot been well studied. plasma. In contrast to the earlier studies, one of the firstand foremost theoretical models to understand the physicsof GRBs was the GRB fireball model, based on a hot,spherically expanding, outflow optically thick to radiation.The fireball model assumes an initially optically thickplasma, with radiation and matter in thermal equilibrium,and the radiation energy density exceeding the rest massenergy density significantly (Meszaros et al. 1993; Piran2004). When these conditions are met, radiation pressuredue to the trapped photons dominates the outflow (orfireball - we will use these terms interchangeably) evolutionleading to an accelerating, ultra-relativistic fireball.In this article, we will investigate from a micro-physicalperspective Comptonization driven acceleration by radia-tion advected with the outflow. We present our DynamicMonte Carlo code (DynaMo code) that we have used to self-consistently investigate the scattering induced acceleration,and relativistic expansion of a spherical fireball. This paperis organized as follows: § § § § §
4, wesummarize our work and draw our conclusions.
The idea of the fireball was first advocated byCavallo & Rees (1978) to explain the phenomenon ofGRBs. Created by depositing a large amount of energyonto matter confined to a small volume, a sufficiently dense fireball will be optically thick ( τ ≫ ) to its ownradiation. This radiation can drive the acceleration andultra-relativistic expansion of the fireball if the radiation en-ergy density exceeds the fireball’s rest mass energy density.Thus, the fireball model resolved the GRB compactnessproblem by employing radiation as a mechanism to drivethe radiation-matter mixture to ultra-relativistic speeds(Paczynski 1986; Goodman 1986). Radiation to matterenergy density and opacity are important factors thatgovern the evolution of a fireball, and this is graphicallyrepresented in Fig. 1 by the phase-space diagram. In thisdiagram, GRB fireballs start optically thick and radiation–dominated in the top–right part of the graph. We refer tothis part of the diagram as quadrant 1 (here exponentsof both logarithmic axes are positive, analogous to firstquadrant in the Cartesian coordinate system where x andy are positive). The fireballs eventually evolve to quadrant3 (bottom–left portion of Fig. 1 where exponents of bothlogarithmic axes are negative, similar to the third quadrantin the Cartesian coordinate system) where they becomeoptically thin and matter–dominated. As we discuss laterin this section, the evolutionary trajectory of a fireball asit evolves from quadrant 1 to 3 depends on the fireball’sinitial parameters.To parameterize the fireball the most important physicalquantities are 1) the total energy E (includes energies ofboth matter and radiation), 2) its rest mass M , and 3) theinitial radius of the fireball R (i.e., fireball evolution beginswith bulk Lorentz factor Γ = at R ). Radiation-matterinteractions transform the fireball’s internal energy intobulk kinematic motion. Assuming that the fireball remainshighly opaque ( τ ≫ ) during expansion, the maximum pos-sible Lorentz factor η attained (using energy conservation)is given by - η = EMc . (1)To characterize the evolution (acceleration and expansion)of the fireball plasma, we require the bulk Lorentz factor Γ and the radius R of the fireball. These physical quantities donot evolve independently during the radiation–dominatedacceleration phase (quadrant 1 in Fig. 1). The mathemati-cal relationship between these quantities and the comovingtemperature T ′ can be obtained using energy, momentum,and entropy conservation, as given by (Meszaros et al. 1993;Piran 2004) – Γ ∝ R ∝ T ′ . (2)Piran (1999) has studied an adiabatically expanding fireballwith infinite opacity using hydrodynamic equations. Theconservation equations for energy, momentum, and numberof particles, are as follows: R ρ Γ = c ′ (3) R e Γ = c ′′ (4) R (cid:18) ρ + e (cid:19) Γ = c ′′′ , (5)where R denotes the fireball radius, and c ′ , c ′′ and c ′′′ areconstants (see Piran 1999 for an exact definition of the vari-ables used). The boundary conditions for this system are 1)the fireball starts from rest at an initial radius R and 2) MNRAS , 1–15 (2015) n Radiative Acceleration in GRB Fireballs as R → ∞ the Lorentz factor equals η . Using these bound-ary conditions, the conservation equations can be solved toobtain - R = R Γ ( η − ) / ( η − Γ ) / . (6)The radiation–dominated acceleration phase can end if 1)the optically thick fireball’s radiation energy density be-comes comparable to or less than the rest mass energy den-sity (i.e., radiation energy does not dominate matter energydensity), and/or 2) radiation escapes at the photopshere dueto the expanding plasma’s decreasing opacity. In the formercase, the trapped photons lose their energy by continuallyaccelerating the plasma to the maximum possible Lorentzfactor η (see eq. 1). The fireball, thereafter, coasts at η andthis new phase is thus termed the saturation or the coastingphase. In the GRB phase space diagram (Fig. 1), this evo-lution is represented by the fireball moving from quadrant 1to 3 via quadrant 4. Using eq. 2, the characteristic radius atwhich the saturation phase begins can be computed as - R sat = η Γ R = η R . (7)In the latter case, the photons escape the fireball if theexpanding plasma’s optical depth does not remain largeenough, thereby, bringing an end to the acceleration process.In the phase space diagram, this evolution is represented bythe fireball moving from quadrant 1 to 2 and eventually to3. The characteristic theoretical radius where radiation de-couples and escapes the plasma (assumed to occur when theoptical depth τ ∼ ) is termed the photospheric radius andis denoted by R ph . The absence of radiation causes the driv-ing force to vanish, and as a result the plasma coasts at theLorentz factor achieved at the photospheric radius.In the previous paragraphs we outlined different evolution-ary paths for GRB fireballs. The fireball’s evolution dependson its initial energy density and opacity, and the effect ofboth these parameters can be mathematically representedby the ratio of the characteristic radii, i.e., R ph / R sat . It isimportant to point out that most earlier works have stud-ied extreme regimes of fireball evolution, e.g., the hydro-dynamical evolution of a fluid with infinite opacity (Piran1999). Paczynski (1990) explored the evolution of super–eddington winds assuming throughout, 1) an optically thickflow (similar to fireball evolution from quadrant 1 to 4),and 2) a constant velocity past the photosphere. Anotherwell studied scenario is when a radiation–dominated fireballreaches the photosphere and then suddenly loses its radia-tion (Meszaros et al. 1993; Rees & M´esz´aros 2005). Earlierworks have not studied the realistic scenario where, due tothe expansion and acceleration, the fireball 1) gradually be-comes optically thin (evolves from quadrant 1 to 2) and losesphotons, and 2) is no longer radiation energy dominated(transitions from quadrant 1 to 4). These transition zonesare represented by the highlighted regions in Fig. 1. In addi-tion, analytical calculations encounter difficulties in study-ing fireball evolution for small opacities (e.g., the Lorentzfactor evolution is unknown in quadrant 2 – as indicatedby the question mark in Fig. 1). In §
3, we will show anddiscuss the results of fireball evolution in all regimes andacross the transition zones. These results are obtained usingour DynaMo code, which we detail in the next section.
In this section we detail the methodology behind the Dy-naMo code, which simulates the Compton scattering drivenexpansion and evolution of an outflow. The outflow is a scat-tering dominated plasma composed of photons, leptons (andprotons). Figure 2 shows a diagrammatic cross-sectionalview of the outflow geometry, which is a conical wedge withan opening angle θ c and encapsulates a section of the fire-ball’s spherical shell. Also shown is a fireball shell as it travelsradially outward, expands and becomes optically thin (thedifferently colored arcs represent the expanding shell at dif-ferent radial positions and opacities).For our simulations, we are primarily concerned withthe physical quantities in three distinct reference frames - • Lab Frame - The lab or the laboratory frame is theframe at rest with respect to the GRB progenitor producingthe outflow or the host galaxy. Any observer at rest in thelab frame observes the jet moving with a bulk Lorentz factor Γ . • Jet Frame - Also known as the comoving frame, thisreference frame travels with the GRB outflow at the bulkLorentz factor Γ . From the perspective of an observer in thelab frame, the jet frame would be in motion along the ra-dial direction (which is along the z direction as depicted inFig. 2). Any physical quantity computed in the jet frame willbe primed, e.g., the four–momentum of an electron in the co-moving frame will be denoted as p ′ µ e − . This frame is naturallysuited for calculating local properties of the fluid/plasma(e.g., comoving temperature T ′ ). • Particle Frame - At any given time, the motion of anyparticle in the plasma will be different from the bulk motionof the jet plasma due to the random thermal motion of theparticle. The particle frame is the frame comoving with thelepton (or proton) selected for the scattering (or collision)process (we will discuss the particle interactions in greaterdepth later in this section). The interaction cross-sections(total and differential) are expressed in terms of physicalquantities defined in the frame of the particle involved inthe interaction.All physical quantities in the particle frame will be doubleprimed, for example, the four–momentum of a photon asobserved in the frame of an electron will be denoted by p ′′ µγ . The first step of our simulation is the initialization of theGRB fireball with the appropriate parameters. As we simu-late a section of a fireball shell encapsulated within a conicalwedge, the following parameters need to be specified to ini-tialize the simulation – • Wedge Opening Angle: θ c = π / radians.To prevent causality violation and unphysical scatterings, itis imperative that the condition θ c ≪ / Γ is always satisfiedduring fireball evolution. • Initial Bulk Lorentz factor: Γ = The GRB fireball initially possesses only internal energy andno bulk kinetic energy. • Wedge Outer Radius: r outer , = cm (Pe’er et al.2015) • Wedge Inner Radius: r inner , = × cm MNRAS , 1–15 (2015)
Chhotray & Lazzati
Figure 2.
Diagram depicting the cross-sectional view of the sim-ulated wedge. The colored arcs represent a shell traveling radi-ally outward (along the z direction), expanding and becomingmore transparent. Also shown is an implementation of the peri-odic boundary condition with particle A exiting and re-enteringthe wedge at an angle ϕ with the lateral surface. • Particle Initialization:– Particle Count: The code can simulate interactionsbetween photons, leptons (electrons) and baryons (pro-tons), and thus requires the number N and type of par-ticles to be specified. Unless otherwise stated, the totalphoton count is N γ = , and the electron and protoncounts are N e = N p = .– Spatial Distribution: The code synthesizes the speci-fied number N of particles, and uniformly distributes themwithin the specified inner and outer radii of the wedge.– Distribution Functions: The user can synthesize par-ticles by supplying the appropriate parameters to the in-built energy distribution functions (e.g., comoving tem-perature T ′ for thermal distribution, power law index p fornon-thermal distribution, the energy for mono-energeticdistribution).Though the code can accept any temperature value, thesimulations shown in this paper start at T ′ = . × K.At this temperature the electrons are distributed accord- ing to a Maxwell J¨uttner distribution, the protons adopta Maxwell Boltzmann distribution, and the photons con-form to a Blackbody distribution.
As one of the motivations of our code is to track and evolvethe plasma and its constituents, it is essential to track thefour–momenta and positions of all the particles in the labframe. In order to dynamically evolve the system a time stepis required, and the maximum value of that step is dictatedby the mean free times of the photons in the plasma. Letus first consider photon-lepton scatterings. From a physicalperspective, all the particles travel in straight lines until aphoton-lepton pair scatters, exchanges energy and momenta,and then the particles again travel in straight lines until thenext photon-lepton scattering occurs. This time interval be-tween successive photon-lepton scatterings is the collision-free time experienced by all particles in the plasma (we shallrefer to this interval as the collision time - t coll ). Thus, thedynamic time step in our simulation is this collision time.The reader should refer to Appendix A for a detailed calcu-lation of the collision time.As experienced by a single photon, the infinitesimal opac-ity of a medium (at rest) as given in Rybicki & Lightman(1979) is – d τ = n σ dl , (8)where n denotes the number density of the scatterers, σ is the effective Thomson cross section (we explain this ingreater detail in Appendix E), and dl is the infinitesimalpath length traversed by the photon.Using eq. 8 and as shown in detail in Appendix A, the meanfree time t mf of the photon is (see eq. A2) – t mf = nc σ, (9)where c = × cm / s, is the speed of light. If the scatteringmedium contains N p photons, the mean free time τ Pop for anyphoton within this population would be given by – t Pop = Í N p i = nc σ = Í N p i = (cid:16) t mf (cid:17) = t mf N p . (10)Physically, if the number of photons in a medium increasesso should the likelihood of a photon colliding, thereby de-creasing the mean free time of the entire population.Now we shall discuss how the opacity and mean free quan-tities are modified if the scatterers themselves are in mo-tion. The motion of scatterers introduces a velocity de-pendence in the opacity. As shown in Appendix A2 andAbramowicz et al. (1991), a single photon immersed in amedium of moving N s scatterers, experiences an opacity d τ given by – d τ = N s Õ j = n j σ ( − β j cos θ j ) dl , (11)where β j is the speed of the scatterer traveling in the j th di-rection normalized by the speed of light and θ j is the anglebetween the three momenta of the photon and the scatterer. MNRAS , 1–15 (2015) n Radiative Acceleration in GRB Fireballs Let us now consider a plasma having multiple photons inter-acting with multiple scatterers. As shown in Appendix A2,the mean free time of the entire population is – t Pop = V Í i Í j σ ( − β j cos θ ij ) c = (cid:16)Í i (cid:16)Í j t ij (cid:17)(cid:17) = Í i (cid:16) t i (cid:17) . (12)Thus all the particles experience an average or mean freetime t Pop given by eq. 12. We note that this is an averagevalue and therefore the extraction of the real collision timeor the actual free travel time between successive scatterings(which is different from the mean free time) requires the useof a Monte Carlo acceptance-rejection algorithm. The prob-ability of one photon scattering between the time interval t and t + dt is proportional to – P ( t ) ∝ e − tt Pop , (13)which we use to obtain the collision time t coll for our simu-lation. The collision time t coll determines the travel time betweenphoton–lepton scatterings. As some photon–lepton pairs aremore likely to scatter than others, the code thus needs toaccount for this likelihood to identify and select the scatter-ing pair.Intuitively, the photon more likely to scatter will have asmaller free travel time as compared to the other photons.The mean free time is a measure of the free travel time fora given photon, and we use each photon’s mean free time(eq. A15) and Monte Carlo techniques to identify and selectthe scattering photon. For any given photon, the probabilityof scattering is inversely proportional to its mean free time.Thus by inverting the mean free time t i (see eq. A15) weobtain a number p i proportional to the likelihood for thatphoton to scatter. For p i to be a probability distribution, wenormalize the inverses obtained from all photons, by obtain-ing the norm A – A = N p Õ i = p i = N p Õ i = t i . (14)By normalizing the inverses using A , we create a probabilitydistribution amenable to Monte Carlo techniques. As eachphoton has a mean free time and an associated scatteringprobability, it is uniquely represented along the probabilitydistribution by a certain range of values. We generate a uni-formly distributed random variable in the interval [0,1) toidentify which photon scatters.The code follows a similar procedure to select the scatteringlepton. To generate the probabilities for each lepton, we usethe mean free time of interaction of the selected photon withthe leptons in the plasma (given by eq. A13). As a result,we obtain the inverses p i , j and the new norm A i - A i = N e Õ j = p i , j = N e Õ j = i , j (15)As illustrated by the photon selection process, each normal-ized probability p i , j / A i represents a range of values, thereby,identifying a unique lepton. We generate a uniform random variable in the interval [0,1), and thus identify which leptonamong the N e leptons scatters with the selected photon. Between the scattering/collision events, the particles freelypropagate for the duration of the time step. The position ofany particle at time t is given by- x i ( t ) = x i ( t − t coll ) + c p i p t coll (16)where x i ( t − t coll ) denotes the i th position component of theparticle at time t − t coll , p i and p denote respectively thespatial and the zeroth component of the four–momenta ofthe particle. As explained in § § ϕ with the exiting surface, then a correspondingparticle is inserted through the wedge surface opposite tothe exit location, but at the same radius (see Fig. 2). Theinserted (re-entering) particle makes the same angle ϕ withthe entering surface as the exiting particle made with theexit surface, and possesses the same energy. This schemeis similar to the periodic boundary condition for parallelsurfaces (where the inclination angle between the surfacesis zero) where the particle momentum vectors at exit andentry surfaces are parallel and hence unrotated. As our sim-ulated wedge has inclined surfaces, a lateral surface crossingrequires the particle’s three momenta be rotated by an angleof θ c (or twice the opening angle of the wedge). It shouldbe noted that energy and total momentum magnitude of theparticle remains constant during this scheme’s implementa-tion. However, the direction of three momentum does notremain constant. MNRAS , 1–15 (2015)
Chhotray & Lazzati
The dynamical energy transfer mechanism in our simula-tions is Compton / Inverse Compton scattering. Once theparticle selection step decides which pair of particles will bescattered, we use the Compton scattering algorithm (Klein–Nishina regime included) discussed in Chhotray & Lazzati(2015) for the scattering event. Post scattering, we have anew four–momenta for both the scattered lepton and pho-ton which we use to calculate a new time step for the nextscattering.
In general, GRB fireballs are baryon loaded (Meszaros et al.1993) making them significantly more massive than baryon–deficient fireballs. We wish to explore whether the baryons(protons) in these baryon–loaded fireballs can accelerate andattain relativistic speeds (just like electrons do). We ignorethe photon–proton ( γ p ) scatterings because the photon en-ergies have to be comparable to rest mass energy of protons( ∼ GeV) for the scattering process to transfer an apprecia-ble amount of energy. To accelerate the protons the DynaMocode accounts for energy transfer between electrons and pro-tons via electron–proton ( e − p ) collisions. To simplify thephysics and the collision process, the code does not carryout e − p collisions using the Coulomb cross section, insteadit performs elastic collisions between electrons and protons.We note that our goal is to explore if protons can be ac-celerated to relativistic speeds, and not to simulate the e − p collision process to perfection. As our results demonstrate(see § m eff = m e instead of realprotons that are 1838 times more massive. This is similar tothe reduced ion-masses employed by multi-dimensional PICsimulations (see Spitkovsky 2008).The Coulomb cross section for interaction between chargedparticles (like e − p interactions) is much larger than the crosssection of γ e − interactions. Due to the attractive nature ofthe forces between these charged particles, there will bemany e − p interactions for every γ e − interaction. In otherwords, for every Compton scattering event, there will beseveral elastic e − p collisions. As a consequence, between suc-cessive γ e − scatterings, the code ensures that each electroncollides with a proton (for every γ e − scattering the numberof e − p collisions equal the number of e − p pairs in the sys-tem).To perform this elastic collision, the code begins by ran-domly selecting one e − p pair from the plasma. To simulatean e − p collision and the involved energy transfer, we firstLorentz transform to the proton’s reference frame and com-pute the transformed four–momentum of the incoming elec-tron. To calculate the transferred energies, we move to thecenter of momentum frame (COM frame) of the system. The R (cm) Γ Γ vs R R ph /R sat = 0.073R ph /R sat = 0.157R ph /R sat = 0.338R ph /R sat = 1.15 R ph /R sat = 11.5R ph /R sat = 115.0R Figure 3.
Evolution of Lorentz factor with radius for fireballscharacterized by different initial opacities. The legend displaysthe R ph / R sat value for each fireball and the corresponding colorrepresenting it. advantage of employing this frame is that the collision hereis always head on and the net three momenta of the e − p sys-tem, pre and post collision is always zero. We Lorentz trans-form the four–momenta of the colliding proton and electronto this frame, which completely specifies the pre–collision ge-ometry. To specify the post–collision geometry, we randomlygenerate the polar angle θ ′ and the azimuthal angle φ ′ in theCOM frame. Now, by employing conservation laws, we canobtain the four–momentum of the proton after collision.The four–momenta thus obtained in the COM frame is de-boosted to get the momenta in the proton’s frame and sub-sequently, in the lab frame. This provides us with the post–collision momenta of both the proton and the electron inthe lab frame. As stated earlier, for every γ e − scatteringevent, the code collides all the e − p pairs once, leading to en-ergy and momentum exchange among these particles. As aresult, though the photons and the protons do not interactdirectly, the electrons act as an intermediary for transferringenergy between radiation and baryons. Once the code performs the scattering and collisions, it up-dates the positions and four–momenta of all particles in-volved in the simulation. After each time step the code com-putes the radius of the farthest and nearest lepton from theorigin, thereby, determining the fireball’s outer and innerradii, respectively. To evaluate if a photon has escaped fromthe fireball, the code compares the radii of all photons withthe fireball’s outer radius. If any photon’s radius exceeds theouter radius of the fireball, that photon is deemed to haveescaped the fireball and is no longer a part of our simula-tion (it will not be involved in the time-step calculation northe scattering process). The code runs until all the photonsescape and are unable to further accelerate the fireball.
MNRAS , 1–15 (2015) n Radiative Acceleration in GRB Fireballs -2 -1 R ph /R sat Γ t e r m TheorySimulation
Figure 4.
Theoretical and simulated values of terminal Lorentzfactor Γ term plotted against different R ph / R sat ratios. The blackcurve plots the theoretical value of the terminal bulk Lorentz fac-tor (obtained from the fireball model). The blue points representthe results obtained from our DynaMo simulations. We have performed simulations of Compton scattering in-duced – radiative acceleration and expansion of a fireball,with and without baryon loading. Figure 3 shows the radialevolution of the bulk Lorentz factor of baryon–deficient fire-balls for several ratios of R ph / R sat (which serves as a proxy forthe opacity). A detailed calculation of R ph / R sat can be foundin Appendix E. These results confirm the idea that greaterthe opacity of a plasma, the greater the terminal Lorentz fac-tor achieved by the plasma due to acceleration by embeddedradiation. Alternatively, as opacity indicates the number ofscatterings occurring in a medium, these results confirm thatan increasing number of scatterings is necessary for a plasmato continue accelerating radiatively by converting its inter-nal energy into bulk motion. The evolution of our simulatedfireballs is in strong agreement with the fireball model’s pro-portionality relation Γ ∝ R , during the radiation–dominatedacceleration phase. In addition to Γ ∝ R , the simulated fire-balls also follow the relation T ′ ∝ R − . The radial evolutionof comoving temperature and other code tests can be foundin Appendix D. An interesting and new feature that emerges from Fig. 3is the transition of the simulated Lorentz factors from the Γ ∝ R (radiation dominated) to Γ ∝ R (matter dominated)regime. As we pointed out in § Γ ∝ R regime occurs because fireball’s energy density is no longerdominated by radiation. Further, we note that the magni-tude of curvature or turnover of the Lorentz factor (the mea-sure of how rapidly Γ transitions from Γ ∝ R to Γ ∝ R ) isopacity dependent. This implies that more opaque fireballs(e.g, blue curve in Fig. 3) transition more smoothly andgradually (with larger curvature) as compared to less opaquefireballs (e.g., green curve). A consequence of this gradualtransition is an increase in the radial extent of the fireballs’ transition zone, which spans several orders of magnitude (itis also larger in extent than the Γ ∝ R regime). Figure 4 plots the terminal Lorentz factor of several fireballsagainst their corresponding R ph / R sat values. For a given R ph / R sat ratio, the black curve depicts the theoreticallycalculated terminal Lorentz factor ( Γ Th ). The blue pointsshow the terminal Lorentz factor obtained by DynaMosimulations (let us call them Γ Sim ). For R ph / R sat < , thefireball model predicts that all radiation escapes at thephotoshere and the Lorentz factor saturates, i.e., Γ Th = Γ ph (see § R ph / R sat < simulated Lorentz factors exceed the cor-responding theoretical values, i.e., Γ Sim > Γ Th . Physically,this excess represents an acceleration phase occurring after the theoretical photosphere. Paczynski (1990) was able toidentify this excess acceleration regime but was unable toanalytically study this regime using the obtained solutions.This is due to the fact that to obtain the outflow solutionsthe flow velocity beyond the photosphere was assumed tobe constant. Our results show (and this is also mentionedin Paczynski 1990) that this assumption breaks downfor outflows that remain radiation dominated until theirphotosphere (thereby experiencing acceleration beyondtheir photospheres).The reason underlying this excess acceleration lies inthe definition of opacity. The opacity of a medium is aprobabilistic quantity and can only provide informationregarding the probability of escape of a photon (Pe’er 2008).However, the theoretical values (such as the Γ Th at thephotospheric radius) are calculated assuming that all thephotons in the fireball escape at τ = (Rees & M´esz´aros2005) and do not account for the opacity’s probabilistic na-ture. By using the particle tracking feature of the DynaMocode we find that even for optically thin fireballs (i.e., τ ≤ )a fraction of photons are still trapped, which concurs withthe probabilistic nature of opacity. These trapped photons(that are unaccounted for in the fireball model) continuescattering beyond the photosphere and are the reason why Γ Sim exceeds Γ Th . It is interesting to note that the smallerthe value of R ph / R sat the larger the discrepancy betweensimulated and theoretical values. This can be attributedto the fact that the (trapped) radiation energy density islarger for smaller R ph / R sat fireballs, which leads to greateracceleration beyond the photosphere. As R ph grows, theradiation energy density decreases and it becomes harderfor less energetic, trapped photons to provide that extrapush, causing Γ Sim to converge to Γ Th . On the fireball phasediagram (Fig. 1), the post–photospheric acceleration phaseoccurs during the transition from quadrant 1 to 3 via 2. According to the standard fireball model, for R ph / R sat ≥ thefireball enters the saturation phase and thereafter coasts at η (see § R ph / R sat ≥ and in contrastto the results of the post–photospheric acceleration phase,Fig. 4 shows that Γ Sim ≤ Γ Th . Another interesting observa- MNRAS , 1–15 (2015)
Chhotray & Lazzati R (cm) Γ R ph /R sat = 0.073R ph /R sat = 0.157R ph /R sat = 0.338R ph /R sat = 1.15R ph /R sat = 11.5 R ph /R sat = 115.0 ProtonsElectrons R Figure 5.
Same as Fig. 3 but with fireballs that are baryonloaded. Also plotted is the Lorentz factor evolution of protons(represented by + markers). tion is that the for increasing R ph / R sat values the discrep-ancy between simulated and theoretical values is reduced,and Γ Sim → Γ Th = η asymptotically. These features can beattributed to the fact that 1) at this stage these fireballsare matter dominated and 2) Comptonization is inefficientfor transferring energy from low energy photons to leptons.The matter–dominated fireball phase implies that the av-erage lepton energy is significantly larger than the averagephoton energy. As seen from the comoving frame, the pho-ton energies are significantly lower than the rest mass en-ergy of electrons and scatterings here occur in the Thomsonregime . In order to extract the final remnants of the sig-nificantly smaller photon energies and reach the maximumallowed Lorentz factor, a large number of (Compton) scat-terings are required (Chhotray & Lazzati 2015). Our resultsshow that once fireballs become matter dominated, Comp-tonization is not efficient in completely converting internalenergy into bulk motion. On the phase–space diagram shownin Fig. 1, the fireball encounters this phase as it evolvesfrom quadrant 1 to 3, transitioning via quadrant 4. It isthis Thomson–dominated phase that is responsible for thegradual fireball acceleration, the gradual Γ turnover (whichbecomes increasingly smooth for higher opacities) and thelarge extent of the transition phase. Thus, optically thickfireballs that become matter dominated accelerate gradu-ally, and require extremely large opacities (or an extremelylarge number of scatterings) to approach η . Figure 5 plots the radially evolving Lorentz factors ofbaryon–loaded fireballs, with each colored curve cor-responding to a unique value of R ph / R sat . The fireballevolution is again in strong agreement with the fireballmodel’s proportionality relations (see eq 2) during theradiation–dominated acceleration phase. We observe againthat fireballs with larger R ph / R sat attain larger Lorentzfactors. Thus, both Figs. 3 & 5 confirm the idea that largeropacity leads to larger acceleration, independently of thefireball’s baryonic content. R (cm) Γ R ph /R sat = 0.037R ph /R sat = 0.073R ph /R sat = 0.157R ph /R sat = 0.338 R ph /R sat = 1.15R ph /R sat = 11.5R ph /R sat = 115.0∞ opacity Figure 6.
Lorentz factor evolution with increasing radii of fire-balls characterized by different initial opacities. The differentlyshaped and colored markers plot eq. 17, which is obtained bycurve fitting (see § § The inclusion of baryons increases the effective mass ofthe plasma and in comparison to baryon–free plasmasrequires more photons for acceleration. This leads to anincrease in computational time and memory consumption.Since both baryon–devoid and baryon–loaded plasmas haveextremely similar Lorentz factor evolution and terminalvalues at saturation, we think it is better to simulate justbaryon–free fireballs as they require less computationaltime and resources.
In this section, we present and discuss the analytical ex-pressions obtained by fitting the simulated data depicted inFig. 3. Each simulation in Fig. 3 is characterized by uniqueinitial outflow conditions – value of injection radius R , max-imum possible bulk Lorentz factor η (attained when all in-ternal energy is converted to bulk motion) and the uniqueparameter σ (a measure of opacity of the plasma; its relationto R ph / R sat is described in detail in Appendix E). The ex-pression we use to model the radial evolution of the Lorentzfactor Γ ( R ) is (the expression used has a form similar toBeuermann et al. 1999) – Γ ( R ) = Γ ∞ ( σ, η ) (cid:20) + (cid:16) r acc ( σ,η ) r (cid:17) s ( σ ) (cid:21) s ( σ ) , (17)which depends primarily on three parameters, which in turnare functions of the initial outflow conditions. These threeparameters are – • s ( σ ) = . − . σ • Γ ∞ ( σ, η ) = (cid:16) log η [ − exp (− . log σ )] . (cid:17) • r acc ( σ, η, R ) = . R η MNRAS000
In this section, we present and discuss the analytical ex-pressions obtained by fitting the simulated data depicted inFig. 3. Each simulation in Fig. 3 is characterized by uniqueinitial outflow conditions – value of injection radius R , max-imum possible bulk Lorentz factor η (attained when all in-ternal energy is converted to bulk motion) and the uniqueparameter σ (a measure of opacity of the plasma; its relationto R ph / R sat is described in detail in Appendix E). The ex-pression we use to model the radial evolution of the Lorentzfactor Γ ( R ) is (the expression used has a form similar toBeuermann et al. 1999) – Γ ( R ) = Γ ∞ ( σ, η ) (cid:20) + (cid:16) r acc ( σ,η ) r (cid:17) s ( σ ) (cid:21) s ( σ ) , (17)which depends primarily on three parameters, which in turnare functions of the initial outflow conditions. These threeparameters are – • s ( σ ) = . − . σ • Γ ∞ ( σ, η ) = (cid:16) log η [ − exp (− . log σ )] . (cid:17) • r acc ( σ, η, R ) = . R η MNRAS000 , 1–15 (2015) n Radiative Acceleration in GRB Fireballs Figure 7.
An updated version of Fig. 1 with the DynaMo code’ssimulation results in context. The differently colored curves rep-resent fireballs starting with different initial opacities. Fireballsthat start comparatively optically thin ( R ph / R sat < ) evolve fromquadrant 1 to 3 via 2. Comparatively optically thick fireballs( R ph / R sat ≥ ) evolve from quadrant 1 to 3 via quadrant 4. Allcurves eventually drop because only photons trapped inside thefireball contribute to the radiation energy density e γ . The Lorentzfactor evolution across the entire phase space is given by eq. 17. Fig. 6 compares the Lorentz factors calculated from the nu-merical expression obtained via curve fitting (i.e., eq. 17 anddepicted by markers) and the simulation results (representedby the solid lines). As the markers and the solid lines show,the numerical and simulated results are in good agreementwith each other. The dashed red line is the Lorentz factorcalculated from eq. 6, which was obtained from the hydro-dynamic evolution of an infinite opacity fireball (see § The mechanism that launches and accelerates the ultra–relativistic GRB jets is a mystery. The outflows inGRBs can be magnetically accelerated if the jets arestrongly magnetized, i.e., the fields carry a significantamount of total energy, in the form of Poynting fluxwhich is converted into (bulk) kinetic energy (Lyubarsky2010;Kumar & Zhang 2015). Several works have shown thatthese jets can be accelerated adiabatically via magneticpressure, e.g. (Contopoulos 1995; Granot et al. 2011). Un-der certain magnetic field configurations, magnetic recon-nection and dissipation can also accelerate GRB outflows(Drenkhahn & Spruit 2002; Drenkhahn 2002; Zhang & Yan2011). For magnetically accelerated outflows, the propor-tionality relation between the Lorentz factor and the outflowradius grows as Γ ∝ R / (Kumar & Zhang 2015), which isdistinctly different from the Lorentz factor evolution of ra-diation dominated outflows (eq. 2). As several authors haveinvestigated magnetized outflows, in our work we focusedon an alternate launching and acceleration mechanism foroutflows, i.e., radiation. In this paper, we have studied the bulk kinematic evolutionof a radiation–driven GRB fireball via Monte Carlo simula-tions. We have presented our novel Dynamic Monte Carlocode (DynaMo) code, which we use to self–consistently sim-ulate, the Compton scattering–induced expansion and accel-eration of a GRB fireball. Earlier works have studied particu-lar phases of fireball evolution analytically (Meszaros et al.1993; Piran 1999). The analytical approximations do notproperly address the question of self consistent fireball evo-lution, especially when 1) the radiation energy density be-comes comparable to or less than the fireball’s rest massenergy density and/or 2) the fireball becomes optically thin(the case when R ph < R sat ). Fig. 7 is an updated version ofthe GRB phase space diagram (Fig. 1) with the evolutionof DynaMo code’s simulated fireballs plotted onto the dia-gram. We have studied fireball evolution across all regimes,including the transition regimes that are the highlighted re-gions in phase diagrams (see Figs. 1 & 7). We summarizeour results as follows – • The evolution of a GRB fireball (or any outflow) canbe summarized by Fig. 7. All radiation–dominated and op-tically thick GRB fireballs start in the upper–right quadrant(quadrant 1) and evolve towards quadrant 3 (optically thin,matter dominated). Depending on where the fireball origi-nates within quadrant 1, it will enter quadrant 3 via pathsthrough quadrants 2 or 4. • We have investigated the effect of initial opacity on fire-ball evolution by simulating several fireballs, each startingwith a different initial opacity (as differentiated by the col-ored curves in Figs. 3 and 5). Our results are in agreementwith the fireball model, as optically thick fireballs achievehigher Lorentz factors, and saturate at larger radii than com-paratively thin ones. • We have also investigated the effects of baryon load-ing of fireballs on the Lorentz factor evolution, as shownby Fig. 5. Our results show identical evolutionary behav-ior of baryon–deficient and baryon–loaded fireballs (Figs. 3and 5 evolve similarly). Being baryon loaded, these fireballsrequire more photons to accelerate the additional mass andhence require more computational time and resources. Ourresults strongly suggest that the evolution of baryon–loadedfireballs can be accurately predicted from the significantlyfaster and less memory intensive, baryon–deficient fireballsimulations. This can also be realized with the aid of eq. 1,as the mass term in the denominator does not differentiatebetween baryonic or leptonic mass. • A remarkable result that can be seen from both Figs. 3and 5 is the existence of a transition regime, with anopacity dependent radial extent. Occuring between theradiation–dominated acceleration phase (where Γ ∝ R ) andthe phase where the Lorentz factor flattens ( Γ ∝ R ), thetransition phase has 1) a significant radial extent (whichexceeds the radial extent of Γ ∝ R phase) and 2) a curvature/ turnover radius which gradually increases with increasingopacity. In Fig. 7, this regime begins within the highlightedregions, where opacity approaches unity and/or radiationenergy no longer dominates the fireball energy density. • The simulation results show the existence of a post–photospheric acceleration phase (see Fig. 4) as predictedby Paczynski (1990), during the previously unexplored
MNRAS , 1–15 (2015) Chhotray & Lazzati R ph < R sat regime. This phase is encountered by the fireballas it travels from quadrant 1 to 3 via 2 (see Fig. 7). Inthis regime, simulated fireballs’ terminal Lorentz factorsare found to be larger than the theoretical Lorentz factorsobtained using the fireball model. In other words, thesimulated fireballs continue accelerating beyond the theo-retical photosphere. We attribute this post–photosphericacceleration phase to the energetic photons still trapped inthe plasma. The particle tracking feature informs us that1) not all photons escape at the theoretical photosphere(when τ ∼ ) and 2) photons at the these small radii aresufficiently energetic. These trapped, energetic photonscontinue to scatter and accelerate the fireball beyond themodel’s theoretical values. • Another interesting result is the existence of theThomson–dominated acceleration phase, detected in theregime R ph ≥ R sat (see Fig. 4). Conversely to the post–photospheric acceleration phase, this phase is character-ized by the theoretical Lorentz factor exceeding the simu-lated terminal Lorentz factor. On the phase space diagram(Fig. 7), this phase arises as the plasma evolves from quad-rant 1 to 3 via 4. This phase can be attributed to the factthat in the optically thick and matter–dominated regime,the average energy of matter is greater than the average pho-ton energy. As a result, energy transfer per scattering is low(Thomson scattering) and for Comptonization to transfersignificant energy from radiation to matter, a large numberof scatterings is required. The Thomson–dominated scatter-ing phase is responsible for the gradual acceleration and thelarge transition / gradual turnover phase for fireballs withlarge R ph ≥ R sat .We have shown that radiation can accelerate fireballsbeyond the photosphere (post–photospheric acceleration– see §
3) and to Lorentz factors larger than previouslyestimated. A powerful implication of our DynaMo code isthat the photosphere cannot be assigned a single value or aradius (and not all radiation escapes at this value), insteadit corresponds to a region or volume from where photonsgradually escape (Pe’er 2008; Beloborodov 2011; Lazzati2016). While escaping through this volume, the radiationscatters and accelerates the plasma leading to larger thanexpected acceleration. Another consequence of our results isthe re–definition of saturation radius. The saturation radiusdefined using eq. 7 does not hold true if the energy densityis not dominated by radiation. Our results show that atthe end of the radiation–dominated acceleration phase,an opaque fireball can still accelerate but only gradually(Thomson–dominated acceleration phase). As a result, thefireball can attain η only for extremely large opacities (orequivalently a large number of scatterings) and at radiisignificantly larger than the saturation radius (see Fig. 4).We quantify the radial evolution of the Lorentz factorby fitting the simulated data using χ technique. Theanalytical expression obtained (see eq. 17) is parametrizedusing the fireball’s initial opacity. The advantage of thisexpression is that it can successfully capture the evolutionof the GRB fireball across its entire phase space. Thisincludes the transition regimes as well as the fireballevolution in quadrant 2. As an example, by supplying therelevant input parameters (such as the initial opacity) to eq. 17, the Lorentz factor at any radius can be obtained foran evolving GRB fireball.Our simulated fireballs are hot and dense enough thatelectron–positron pair processes can become important inchanging the photon and lepton counts. At the beginningof the fireball evolution the changing particle count onlyserves to increase the opacity which decreases the escapeprobability of the photons. As a result, radiation andmatter remain in equilibrium as the fireball expands,cools and eventually reaches a temperature where pairprocesses become unimportant and thus can be ignored.In all our simulations, pairs become irrelevant long beforeradiation escapes the plasma and thus pair processes arenot accounted for.The results obtained in this paper will be useful for studyingphotospheric emission from GRBs. Though this paper onlyfocuses on GRB fireballs, the DynaMo code can be usedfor studying radiative acceleration in relativistic outflowsfrom other astrophysical environments, such as AGNs andmicroquasars. The DynaMo code’s particle tracking featureallows the user to not only obtain the spectrum of theescaping radiation, but also determine when and where thephotons escape the outflow. The code can thus producetime–resolved spectra and light curves from GRB fireballs,which will be the subject of future publications. ACKNOWLEDGEMENTS
We thank the anonymous reviewer for suggestions to im-prove the manuscript. This work made use of NumPy(van der Walt et al. 2011) for computation, and Matplotlib(Hunter 2007) for preparing figures. This work was sup-ported in part by NASA Swift Grant NNX15AG96G.
REFERENCES
Abramowicz M. A., Novikov I. D., Paczynski B., 1991, ApJ, 369,175Beloborodov A. M., 2011, ApJ, 737, 68Beuermann K., et al., 1999, A&A, 352, L26Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433Castor J. I., Abbott D. C., Klein R. I., 1975, ApJ, 195, 157Cavallo G., Rees M. J., 1978, MNRAS, 183, 359Chhotray A., Lazzati D., 2015, ApJ, 802, 132Contopoulos J., 1995, ApJ, 450, 616Drenkhahn G., 2002, Astronomy & Astrophysics, 387, 714Drenkhahn G., Spruit H. C., 2002, Astronomy and Astrophysics,391, 1141Frail D. A., Kulkarni S. R., Nicastro L., Feroci M., Taylor G. B.,1997, Nature, 389, 261Frail D. A., et al., 2001, ApJ, 562, L55Goodman J., 1986, ApJ, 308, L47Granot J., Komissarov S. S., Spitkovsky A., 2011,Monthly Notices of the Royal Astronomical Society, 411,1323Hunter J. D., 2007, Computing in Science & Engineering, 9, 90Klebesadel R. W., Strong I. B., Olson R. A., 1973, ApJ, 182, L85Komissarov S. S., 2011, Mem. Soc. Astron. Italiana, 82, 95Kumar P., Zhang B., 2015, Phys. Rep., 561, 1Lazzati D., 2016, ApJ, 829, 76Lyubarsky Y. E., 2010, MNRAS, 402, 353Lyutikov M., Pariev V. I., Blandford R. D., 2003, ApJ, 597, 998MNRAS000
Abramowicz M. A., Novikov I. D., Paczynski B., 1991, ApJ, 369,175Beloborodov A. M., 2011, ApJ, 737, 68Beuermann K., et al., 1999, A&A, 352, L26Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433Castor J. I., Abbott D. C., Klein R. I., 1975, ApJ, 195, 157Cavallo G., Rees M. J., 1978, MNRAS, 183, 359Chhotray A., Lazzati D., 2015, ApJ, 802, 132Contopoulos J., 1995, ApJ, 450, 616Drenkhahn G., 2002, Astronomy & Astrophysics, 387, 714Drenkhahn G., Spruit H. C., 2002, Astronomy and Astrophysics,391, 1141Frail D. A., Kulkarni S. R., Nicastro L., Feroci M., Taylor G. B.,1997, Nature, 389, 261Frail D. A., et al., 2001, ApJ, 562, L55Goodman J., 1986, ApJ, 308, L47Granot J., Komissarov S. S., Spitkovsky A., 2011,Monthly Notices of the Royal Astronomical Society, 411,1323Hunter J. D., 2007, Computing in Science & Engineering, 9, 90Klebesadel R. W., Strong I. B., Olson R. A., 1973, ApJ, 182, L85Komissarov S. S., 2011, Mem. Soc. Astron. Italiana, 82, 95Kumar P., Zhang B., 2015, Phys. Rep., 561, 1Lazzati D., 2016, ApJ, 829, 76Lyubarsky Y. E., 2010, MNRAS, 402, 353Lyutikov M., Pariev V. I., Blandford R. D., 2003, ApJ, 597, 998MNRAS000 , 1–15 (2015) n Radiative Acceleration in GRB Fireballs Madau P., Thompson C., 2000, ApJ, 534, 239McKinney J. C., 2006, MNRAS, 368, 1561Meszaros P., Laguna P., Rees M. J., 1993, ApJ, 415, 181Paczynski B., 1986, ApJ, 308, L43Paczynski B., 1990, ApJ, 363, 218Pe’er A., 2008, ApJ, 682, 463Pe’er A., Barlow H., O’Mahony S., Margutti R., Ryde F., LarssonJ., Lazzati D., Livio M., 2015, ApJ, 813, 127Piran T., 1999, Phys. Rep., 314, 575Piran T., 2004, Reviews of Modern Physics, 76, 1143Rees M. J., M´esz´aros P., 2005, ApJ, 628, 847Rybicki G. B., Lightman A. P., 1979, Radiative processes in as-trophysics. Wiley-InterscienceSmith N., Owocki S. P., 2006, ApJ, 645, L45Spitkovsky A., 2008, ApJ, 673, L39Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS,418, L79Zampieri L., Turolla R., Foschini L., Treves A., 2003, ApJ, 592,368Zhang B., Yan H., 2011, ApJ, 726, 90van der Walt S., Colbert S. C., Varoquaux G., 2011,Computing in Science & Engineering, 13, 22
APPENDIX A: TIME STEP CALCULATIONDETAILS
In this section we detail our calculation of the DynaMocode’s dynamic time step. We begin with a discussion of theopacity experienced by a photon immersed in a collection ofscatterers at rest. We use this opacity to calculate the meanfree path for a single photon. We then extend this analysisto calculate the opacity of a medium containing several pho-tons and the mean free time of this population of photons.We then detail how to incorporate the motion of scatterersinto the opacity and mean free time calculations. As before,we derive the expression for mean free time for a single pho-ton immersed among scatterers in motion. We use this resultto find the mean free time for a population of photons in amoving medium.
A1 For scatterers at rest
Using eq. 8 (which is valid for scatterers at rest) as the defin-ing equation of opacity, the mean free path ( l mf ) traversedby the photon can be computed to be – l mf = n σ . (A1)We take advantage of the constancy of the speed of light toobtain the mean free time t mf as – t mf = nc σ . (A2)The reader should note that this is the mean or average timethat the photon travels between scatterings. For a mediumcontaining N p photons, the total infinitesimal opacity can beobtained by summing the individual contribution from eachphoton – d τ Pop = N p Õ i = d τ i = N p Õ i = n σ dl i = n σ N p Õ i = dl i = n σ l mf , Pop N p , (A3)where the subscript Pop is used for population and l mf , Pop de-notes the mean free path for the photon population. Similar to eq. A1, we can express the mean free path of the entirepopulation as – l mf , Pop = n σ N p = lN p . (A4)The mean free time for this population can be computed tobe – t Pop = l mf , Pop c = N p cn σ = t mf N p . (A5) A2 For scatterers in motion
In this section we derive the opacity and mean quantitieswhen the scatterers in the medium are moving (i.e., themedium itself is moving) and compare them with the samequantities obtained when the scatterers are at rest. We be-gin with the case of just a single photon and then extendour calculation to a system containing multiple photons.Consider the case where each of the scatterers travel at aunique velocity and hence contribute uniquely to the opac-ity. The differential opacity seen by a photon due to scat-terers traveling at a distinct velocity given by ® β k is (seeAbramowicz et al. 1991) – d τ k = n k σ ( − ® β k . ˆ r ) dl (A6) d τ k = n k σ ( − β k cos θ k ) dl (A7)where n k denotes the number density of scatterers travelingin the k th direction with a speed β k , ˆ r is the unit vectoralong the direction of propagation of the photon and θ k isthe angle between the photon’s direction of propagation and k th direction.The total opacity as observed by the photon will be the sumof the opacities due to each individual scatterer and can bewritten as – d τ = Õ k τ k = Õ k n k σ ( − ® β k . ˆ r ) dl (A8) d τ = Õ k τ k = Õ k n k σ ( − β k cos θ k ) dl (A9)where Í k denotes the sum over the directions and speeds as-sociated with the scatterers. These expressions can be usedto write the mean free time experienced by a photon travel-ing through a moving medium as – t mf = Í k n k c σ ( − β k cos θ k ) . (A10)The reader should note that on comparison with a mediumat rest, the photon observes a decrease in opacity (and anincrease in the mean free time) if it moves along the direc-tion of motion of the scatterer. For motion anti-parallel tothe scatterer’s motion, an increase in opacity is experiencedand as a consequence, the mean free time is reduced.If the total number of scatterers (say leptons) is fixed anddenoted by the number N e , we can re-write eq. A9 by elim-inating n k as – d τ = Í N e k σ ( − β k cos θ k ) V dl , (A11)where V denotes the volume occupied by the scatterers. Re-writing the mean free time obtained in eq. A10, we obtain MNRAS , 1–15 (2015) Chhotray & Lazzati – t mf = l mf c = V σ c Í N e k = ( − β k cos θ k ) . (A12)Let us now extend the above analysis to the most generalscenario – a plasma containing N p photons and N e leptons(which play the role of scatterers), and each lepton canbe in motion. Similar to the calculation performed in Ap-pendix A1, the opacity experienced by the i th photon due to any lepton traveling along the j th direction is – d τ i , j = σ ( − β j cos θ ij ) V dl . (A13)Using eqs. A12 and A13, the mean free time for the interac-tion between i th photon – j th lepton pair becomes – t i , j = l ij c = V σ c ( − β j cos θ ij ) . (A14)Using eq. A12 and A14, the mean free time experienced bythe i th photon as it interacts with all N e leptons is – t i = V Í N e j = σ ( − β j cos θ ij ) c = V σ c Í j ( − β j cos θ ij ) = Í j (cid:16) t ij (cid:17) . (A15)Most generally, the net infinitesimal opacity of the photonpopulation can be written as – d τ net = N p Õ i = d τ i = N p Õ i = N e Õ j = σ ( − β j cos θ ij ) dlV . (A16)The use of eqs. A12 and A16 provides us with the mean freetime for the entire population, which is – t mf , Pop = V σ c Í i Í j ( − β j cos θ ij ) = Í i Í j (cid:16) t ij (cid:17) = Í i (cid:16) t i (cid:17) . (A17) APPENDIX B: RADIUS AND BULK LORENTZFACTOR CALCULATION
As stated in § R , bulk Lorentz factor Γ and thetemperature T ′ measured in the frame comoving with thefireball. In this section we illustrate how we compute thefireball’s radius and bulk Lorentz factor from the individualpositions and momenta of the matter particles (leptons). Inother words, we show how to find the radius and velocity ofthe Center of Momentum (COM) frame associated with thesystem of particles.As the simulated fireball is a thin–shell, its constituent par-ticles are distributed within the shell’s finite width. Let usdenote the i th particle’s position as r i and velocity v i . Gen-erally, the particles do not travel along the radial direction(which lies along the ˆ z direction as explained in §
2) becauseour simulation’s origin is not coincident with each particle’sposition. As a consequence, we can resolve each particle’snet velocity parallel and perpendicular to the radial vector at the particle’s location. The radial components will con-tribute to the bulk motion whereas the non–radial compo-nents will provide us with a means to measure the temper-ature of the particles. The radial component of the velocitycan be written as – v rad , i = ® v i . ˆ z . (B1)In general, the radial direction at each particle’s position ® r i is different. However, if the angular size of the shell (in otherwords, the opening angle of our simulated fireball) is smallthen the positions of all particles lie within a narrow cone,and their radial vectors are approximately aligned. Usingeq. B1, we can compute the radial Lorentz factor as – γ rad , i = r − v , i c . (B2)The radius of the fireball R (or the radius of the COM) canbe computed from the positions of the constituent particlesas – R = Í i γ rad , i ® r i . ˆ r i Í i γ rad , i . (B3)To calculate the bulk Lorentz factor Γ , we start with theequation for total momentum of the system which gives us– Γ MV = Õ i γ rad , i m i v rad , i (B4)where V = c p − / Γ . This equation can be simplified to – Γ β = Í i γ rad , i m i β rad , i M = Í i γ rad , i β rad , i N , (B5)where we β = V / c , β rad , i = v rad , i / c and M = Í Ni m i = Nm .Using the relation ( Γ β ) = Γ − , we can find the bulk Lorentzfactor from the radial Lorentz factors of individual particlesas – Γ = s(cid:18) Í i γ rad , i β rad , i N (cid:19) + . (B6) APPENDIX C: COMOVING TEMPERATURECALCULATION
As the fireball accelerates and expands, it uses internal en-ergy to fuel its bulk kinetic motion. The temperature of anobject is a measure of the internal energy (attributed to ran-dom kinetic energy) content of that object. The net motionof the particles forming the fireball, in particular, the mat-ter (e.g., leptons) – is a composite of bulk motion (due tooutward directed radiative acceleration) and random motion(due to temperature). Thus, determination of particles’ tem-perature depends upon the particles’ random speeds whichcan be correctly calculated only after accounting for theirbulk velocity. As a result, we need to calculate the veloci-ties of the particles as observed in the comoving frame andcompute the temperature in that frame.As discussed in section 2.2.1, the fireball is a thin sphericalshell that is expanding radially outward. The constituentparticles of this shell are distributed randomly within the
MNRAS000
MNRAS000 , 1–15 (2015) n Radiative Acceleration in GRB Fireballs shell and interact with other particles in the shell. The ra-dial speed of the center of mass can be computed as – V COM = Í i γ i m i v i , rad Í i γ i m i , (C1)where v i , rad , γ i and m i are the radial speed, Lorentz factor,and mass respectively, of the i th lepton and are defined inAppendix B. Due to the small opening angle of our simulatedfireball, the COM velocity (bulk velocity) can be written as– ® V COM = V COM ˆ z . (C2)Using this bulk velocity, we can Lorentz transform the four–momenta of the particle to the COM frame. This transfor-mation removes the radial–bulk motion component from theparticle momentum, leaving behind the momentum only dueto random motion. Let us denote the random velocity of the i th lepton (as measured in the COM frame) by v ′ i / COM . TheLorentz factor associated with this random motion is – γ ′ i / COM = r − v ′ / COM c . (C3)The average random kinetic energy or internal energy (de-noted by < KE ′ > ) of the system can be computed as – < KE ′ > = Í N e i = KE ′ i N e = N e Õ i = (cid:16) < γ ′ i / COM > − (cid:17) mc N e , (C4)where KE ′ i is the random kinetic energy of the i th leptonand N e denotes the total number of leptons present in thefireball. As the average kinetic energy and temperature arerelated by – < KE ′ > ≃ k B T ′ , (C5)the average temperature T ′ of the leptons can be computedas – T ′ ≃ Í N e i = (cid:16) < γ ′ i / COM > − (cid:17) mc k B N e , (C6)where k B denotes the Boltzmann’s constant. APPENDIX D: CODE TESTSD1 Temperature Evolution
This section explores the evolution of comoving tempera-ture T ′ (as measured in the jet frame) with the radius of thefireball. Figure D1 depicts the comoving temperature evolu-tion of baryon–free fireballs, parameterized by their differentinitial opacities (as discussed earlier in §
3, baryon–loadedfireball evolution is similar to the evolution of baryon–freefireballs, and hence the results of one hold true for the other).Fig. D1 shows that during the radiation–dominated acceler-ation phase the comoving temperature of all fireballs growsproportionally to R − , thus verifying the fireball model’s pro-portionality relations (see eq. 2). The fireballs deviate from T ′ ∝ R − path when the radiation dominance ends (as dis-cussed in § R (cm) T ′ ( K ) R ph /R sat = 0.073R ph /R sat = 0.157R ph /R sat = 0.338R ph /R sat = 1.15 R ph /R sat = 11.5R ph /R sat = 115.0R −1 Figure D1.
Evolution of the comoving temperature with thefireball radius. Each colored curve represents a fireball that beginsevolution with a unique initial opacity. Also shown (by the dashedline) are the fireball model’s theoretical predictions during theradiation–dominated acceleration phase. R (cm) T ′ ( K ) R ph /R sat = 0.157R ph /R sat = 1.15R ph /R sat = 115.0 R −1 R −4/3 R −2/3 Figure D2.
Comoving temperature evolution with radius for fire-balls with R ph / R sat = . , .
15 & 115 . The black circles depict aline that decays proportionally to R − / , the red squares depicta decay proportional to R − / , and the black dashed line repre-sents a curve evolving proportionally to R − . The deviations ofthe evolutionary trajectories from the theoretical predictions arediscussed in § D1. Also plotted are the error bars for the simula-tions computed from their standard deviation.
Figure D2 investigates these deviations in greater detail.Similar to Fig. D1, Fig. D2 plots the comoving temperaturewith the fireball radius, but only for three fireballs (insteadof six) characterized by R ph / R sat = . , .
15 & 115 . . Onlythree fireballs are plotted to improve clarity and develop abetter understanding of the comoving temperature devia-tions. For each fireball displayed in Fig. D2, the transitionfrom the radiation–dominated phase to the transient phaseproduces different evolutionary relationships between T ′ and R . If the fireball transitions from quadrant 1 to 3 via 2 (seeFig. 7, during the post–photospheric acceleration phase) thereduced optical depth causes radiation to decouple. As a re- MNRAS , 1–15 (2015) Chhotray & Lazzati R (cm) ∆ ( c m ) R ph /R sat = 115. 0R R Figure D3.
Evolution of the width of a fireball shell (for R ph / R sat = ) as the fireball expands. As predicted by theory, theshell maintains a constant width during the acceleration phase. sult, the cyan curve (the least optically thick of the plottedfireballs) drops below the T ′ ∝ R − dashed line, and begins toevolve along the black dotted circles representing T ′ ∝ R − / path.The blue curve follows a different evolutionary path, a con-sequence of the blue fireball transitioning from quadrant 1to 3 via quadrant 4 (see Fig. 7). The blue curve represents ahighly opaque fireball that achieves matter domination be-fore transparency, and as a result, undergoes the Thomson–dominated acceleration phase (see § T ′ ∝ R − / . During the matter–dominated phase, the bluecurve in Fig. D2 jumps above the T ′ ∝ R − dashed line, andfollows the red squares. However, when the fireball even-tually enters the optically thin regime (transitioning fromquadrant 4 to 3) it begins following the black circles.The yellow curve represents a fireball ( R ph / R sat = . ) thatis comparatively optically thick than the cyan one. As aresult, its initial deviation (from the radiation–dominatedphase) is similar to blue curve than the cyan curve. Fig. 7shows that the yellow curve becomes matter–dominated be-fore becoming optically thin, and transitions from quadrant1 to 3 via 4. As a result, the yellow curve jumps over the T ′ ∝ R − (dashed) line just like the blue curve, but onlyslightly. As the fireball represented by the yellow curve isnot as optically thick as the one represented by the bluecurve, the yellow curve becomes optically thin at a com-paratively smaller radius. As a result, it drops below the T ′ ∝ R − line earlier than the blue fireball, and begins evolv-ing along the black circles. Thus, we find that the fireball’stransient phases can explain the temperature evolution andits deviation from T ′ ∝ R − proportionality relation. D2 Shell Width Evolution
Meszaros et al. (1993) found that the width of a fireball shellin the lab–frame remains constant during the radiation–dominated acceleration phase. The shell starts expandingproportionally to R during the matter–dominated phase. R (cm) Γ OriginalRestart R Figure D4.
Evolution of Lorentz factor with radius for fireballswith R ph / R sat = . . The blue curves represents the original sim-ulation, which starts with no bulk motion. The red curve repre-sents the simulation ‘restarted’ from the original simulation withinitial conditions extracted at Γ ∼ . Also plotted are the errorbars for both the simulations computed from their standard de-viation. Figure D3 depicts the evolution of the width of the shell ∆ in the lab frame and its evolution with the fireball radius,and it shows that the model is in agreement with the theorybarring the transition phase. D3 Restarted Simulations
As shown in §
3, the DynaMo code reproduces the propor-tionality relations of the fireball model (see eq. 2). A power-ful feature of the code is the ability to start (or restart) thesimulations at a specific point along the simulation (posi-tions and four–momenta of particles are all that is requiredto initialize the simulation). In this section we discuss theresults obtained by restarting a completed simulation at aspecific point along its evolution and compare them with theresults of the original simulation.The simulations we are considering in this section are speci-fied by R ph / R sat = . (see Fig. 3). The specific data pointalong the evolution which we use for the starting conditionsfor the ‘restarted simulation’ refers to the point when thebulk Lorentz factor of the original simulation equals 3. Asshown in Figs. D4 and D5, the two simulations evolve con-sistently as expected, and both simulations are in agreementwith the proportionality relations of the fireball model. APPENDIX E: CALCULATION OF THE R PH AND R SAT
A relativistically expanding fireball achieves saturationwhen all its internal energy is converted into bulk motion,and as a result its bulk Lorentz factor achieves the maximumpossible value. Using the fireball initialization parametersdefined in § η ∼ . . In this section wecalculate the photospheric radius R ph in terms of the param-eters defined in § MNRAS000
A relativistically expanding fireball achieves saturationwhen all its internal energy is converted into bulk motion,and as a result its bulk Lorentz factor achieves the maximumpossible value. Using the fireball initialization parametersdefined in § η ∼ . . In this section wecalculate the photospheric radius R ph in terms of the param-eters defined in § MNRAS000 , 1–15 (2015) n Radiative Acceleration in GRB Fireballs R (cm) T ′ OriginalRestart R −1 R −4/3 Figure D5.
Evolution of comoving temperature with radius forfireballs with R ph / R sat = . . The blue curves represents theoriginal simulation which started with no bulk motion. The redcurve is started with initial conditions extracted from the dataof the original simulation (when Γ ∼ ). Also plotted are theerror bars for the two simulations (computed from their standarddeviation). differently during the acceleration and saturation phases ofthe fireball evolution, therefore we consider two cases: 1)the fireball attains the photospheric radius before satura-tion, i.e., R ph < R sat and 2) when fireball achieves saturationbefore reaching the photospheric radius R ph ≥ R sat . E1 Case I: R ph / R sat < The first case considers the scenario that the photosphericradius is reached before the saturation radius, i.e., R ph < R sat .The saturation radius R sat (the radius at which the fireballattains η ) has already been obtained in eq. 7 as – R sat = η R Γ = EMc R Γ . (E1)The calculation of the photospheric radius R ph begins withthe calculation of opacity of the evolving fireball. As statedin Abramowicz et al. (1991), the lab frame opacity τ is givenby – τ = ∫ ∞ R ph n σ ( − β cos θ ) dR , (E2)where n is the number density of scattering particles andthe ( − β cos θ ) term accounts for velocity dependence ofopacity ( β is the bulk speed of the outflow normalized bythe speed of light). Assuming nR = n R = constant (where n is the initial particle density), and for cos θ ∼ (smallangle approximation), eq. E2 can be written as – τ = ∫ ∞ R ph n R σ R Γ dR . (E3)If the fireball is within the acceleration regime (see eq. 2),we can write – Γ ph R = Γ R . (E4) Combining eqs. E3 and E4, the opacity equation becomes – τ = n R σ ∫ ∞ R ph R . (E5)The photospheric radius R ph is attained when τ ∼ , whichimplies – R ph = n R σ ! / . (E6)The ratio of the photospheric and the saturation radii canbe written as – R ph R sat = n R σ ! / η R . (E7)To vary the results of eq. E7 we use σ as an effective cross–section. The connection between the effective cross–sectionalparameter σ and the Thomson cross–section σ T is – σ ≡ σ T σ coeff V c , (E8)where V is the initial volume of the fireball shell and σ coeff is the parameter we employ to scale the cross–section (andregulate the number of scatterings). Re–writing eq. E7 usingeq. E8 – R ph R sat = n R σ ! / η R = N R σ eff c ! / η R . (E9)In our simulations we have kept the injection radius of, andthe number of scatterers N (and therefore the initial numberdensity n V = N ) in the fireball as fixed quantities. In orderto explore fireballs with different opacity, we have to varythe quantity σ by varying the parameter σ coeff as they areproportional. The physical effect of increasing (decreasing)this parameter is to increase (decrease) the opacity (andas a consequence also the number of scatterings). As anexample, for σ coeff = the ratio R ph / R sat = . , and R ph / R sat = . results from σ coeff = (see fig. 3). E2 Case II: R ph / R sat ≥ We now explore the case when the R ph ≥ R sat . At the satu-ration radius, as all the internal energy has been convertedinto bulk kinetic energy, the acceleration of the fireball ceasesand the fireball coasts beyond that radius at η . If the sat-uration radius is smaller than the photospheric radius, thisimplies radiation is still trapped within the fireball. As thefireball is no longer in the acceleration regime (see eq. 2), Γ ph = η , and the fireball opacity evolves differently from thecase R ph < R sat . Using eq. E2, we can write – τ = = ∫ ∞ R ph n R σ R η dR . (E10)The photospheric radius obtained by integration can be ex-pressed as – R ph = n σ R η . (E11) MNRAS , 1–15 (2015) Chhotray & Lazzati
Therefore, the ratio of the photospheric to the saturationradii is – R ph R sat = n σ R η . (E12)Using eq. E8, we can rewrite eq. E12 as – R ph R sat = n σ R η = N σ eff V R cV η = N σ eff R c η . (E13)As discussed in Appendix E1, by tuning the parameter σ coeff we can regulate the fireball opacity. For example, by using σ coeff = we produce the ratio R ph / R sat = . . Similarly,by increasing σ coeff to we obtain an increased ratio R ph / R sat = . This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000