Monte Carlo radiation transfer simulations of photospheric emission in long-duration Gamma-Ray Bursts
DDraft version October 13, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
MONTE CARLO RADIATION TRANSFER SIMULATIONS OF PHOTOSPHERIC EMISSION INLONG-DURATION GAMMA-RAY BURSTS
Davide Lazzati Draft version October 13, 2018
ABSTRACTWe present MCRaT, a Monte Carlo Radiation Transfer code for self-consistently computing thelight curves and spectra of the photospheric emission from relativistic, unmagnetized jets. We applyMCRaT to a relativistic hydrodynamic simulation of a long duration gamma-ray burst jet, and presentthe resulting light-curves and time-dependent spectra for observers at various angles from the jet axis.We compare our results to observational results and find that photospheric emission is a viable modelto explain the prompt phase of long-duration gamma-ray bursts at the peak frequency and above,but faces challenges in reproducing the flat spectrum below the peak frequency. We finally discusspossible limitations of these results both in terms of the hydrodynamics and the radiation transferand how these limitations could affect the conclusions that we present.
Subject headings: gamma-ray burst: general — radiation mechanisms: thermal — radiative transfer INTRODUCTION
Photospheric emission is a leading model to explainthe prompt radiation of long-duration gamma-ray bursts(Rees & M´esz´aros 2005; Giannios 2006; Lazzati et al.2009; Lundman et al. 2013; Peer 2015). In the photo-spheric model the radiation is assumed to form deep inthe outflow, where the jet is highly optically thick, and tobe advected out, eventually being released when the jetbecomes transparent. The photospheric model does notspecify the radiation mechanism involved, since the spec-trum is shaped by the interaction with matter, mostlythrough Compton scattering, and not by the emissionprocess itself.The alternative model, the Synchrotron Shock Model(SSM, Rees & Meszaros 1994; Daigne & Mochkovitch1998; Lloyd-Ronning & Petrosian 2002; Daigne et al.2011; Zhang & Yan 2011), assumes instead that the jetcrosses the photosphere virtually deprived of radiation,coasts until the Thomson optical depth is well belowunity, and is reactivated when shells of different speedshock each other and generate non-thermal particles andmagnetic field (Rees & Meszaros 1994). SSM predictsthat the radiation spectrum is directly related to theradiation mechanism and suffers little interaction withthe baryons and leptons of the jet once it is produced.Both models have been successful in explaining some ofthe observed properties of long duration GRBs but arein significant tension with other properties. For exam-ple, SSM easily accounts for the non-thermal, broadbandnature of the prompt spectrum (Piran 1999), but faceschallenges for its intrinsic width (Vurm & Beloborodov2015) and when compared to ensemble correlations suchas the Amati, Yonetoku, Golenetskii, or E-Γ correlations(Golenetskii et al. 1983; Amati et al. 2002; Yonetoku etal. 2004; Liang et al. 2010; Ghirlanda et al. 2012). On theother hand, the photospheric emission has been shown tobe able to reproduce all the correlations (Lazzati et al.2011, 2013, hereafter L13; L´opez-C´amara et al. 2014) but Department of Physics, Oregon State University, 301Weniger Hall, Corvallis, OR 97331, U.S.A. can hardly produce a spectrum that has prominent non-thermal tails at both low- and high-frequencies. Whilesub-photospheric dissipation can cure the high-frequencyproblem producing prominent non-thermal tails (Pe’eret al. 2006; Giannios 2006; Giannios & Spruit 2007; Be-loborodov 2010; Lazzati & Begelman 2010; Vurm et al.2011; Ito et al. 2013, 2014; Chhotray & Lazzati 2015;Santana et al. 2016), it is still unclear whether the low-frequency tails can be reproduced in photon-starved con-ditions without invoking an external radiation mecha-nism such as synchrotron (Pe’er & Ryde 2011; Chho-tray & Lazzati 2015). Encouraging results, however,have been obtained for photospheric emission with con-tributions from an subdominant synchrotron component(Vurm & Beloborodov 2015).In both cases, theoretical models only partially ad-dress the complexity of the GRB phenomenon. Withvery few exceptions, models either simplify the radiationphysics or the geometry of the emitting region. For ex-ample, photospheric radiation computed from jet simula-tions – and therefore fully taking into account the com-plexity of structured and turbulent outflows – assumethat the radiation-matter coupling is infinite below thephotosphere but suddenly drops to zero as soon as thephotosphere is crossed (e.g., Lazzati et al. 2009; Mizutaet al. 2011; Nagakura et al. 2011; L13). On the otherhand, models that consider proper radiation transfer orradiation from magnetized jets with dissipation and par-ticle acceleration are applied to simplified homogeneousplasma regions or to spherical jets with analytic radialstratification and without any polar structure (e.g., Pe’eret al. 2006; Giannios & Spruit 2007; Beloborodov 2010;Lazzati & Begelman 2010; Vurm et al. 2011; Chhotray &Lazzati 2015; Vurm & Beloborodov 2015; Beloborodov2016).In this paper we present the results of a first steptowards merging the two worlds and performing de-tailed radiation transfer calculations on the complex dy-namics of a relativistic jet from high-resolution multi-dimensional numerical simulations (analogous to the re-cent work by Ito et al. 2015). In Section 2 we present a r X i v : . [ a s t r o - ph . H E ] J u l Lazzatithe Monte Carlo Radiation Transfer (MCRaT) code andthe results of several validation tests. In Section 3 wepresent and discuss the results of applying the code to ahydrodynamic (HD) simulation of a GRB jet. Finally, inSection 4 we discuss our results by comparing them withobservations and previous theoretical work, and we dis-cuss limitations and possible improvements of this tech-nique. THE MCRAT CODE
We developed a Monte Carlo Radiation Transfer(MCRaT code) to compute time-resolved light curvesand spectra from hydrodynamic simulations of GRB jets.The approach to the single-photon scattering in the codeis to neglect the energy losses of the electron population,randomly selecting at each scattering an electron from athermal population at a temperature that does not de-pend on previous scatterings. This is a fairly standardapproach that assumes the heat capacity of the fluid to belarge enough that radiation losses can be neglected (e.g.Giannios 2006; Ito et al. 2013). The novelty of the codeis that of performing the radiation transfer on the fastevolving background of a relativistic outflow (Ito et al.2015). The code injects photons in thermal equilibriumwith the fluid well inside the photosphere and performsrepeated Compton and inverse Compton (IC) scatteringoff electrons that are at the local comoving temperatureof the fluid. The present version of MCRaT takes two-dimensional simulations as input, but the photon propa-gation is done in three dimensions to properly take intoaccount the photon propagation angle with respect to thefluid. In other words the photon field includes photonspropagating in directions with non-zero azimuthal angle φ , even if the fluid velocity is constrained to be only inthe radial and polar direction. In more detail, the codeperforms the following steps (see also Figure 1): • Photon injection –
Photons are injected at a userspecified radius R inj and within a user specifiedrange of polar angles from the jet axis. The codedetermines the local temperature, density, and ve-locity vector at the photon’s position by nearest-neighbor interpolation. This can be switched tolinear interpolation at the cost of a longer compu-tational time, but with no visible effect on the re-sults. The photon 4-momentum p γ is generated byrandomly sampling a thermal photon distributionat the local fluid temperature. The distribution canbe chosen between a black body and a Wien spec-trum, the former being appropriate for fireballs inthe acceleration stage, the latter for fireballs in thecoasting phase. The photons’ direction of propa-gation is assumed to be isotropic in the fluid’s restframe at injection. Since it is clearly impossible tosimulate all the photons in a GRB fireball (theirnumber amounting to ∼ for a typical burst),each photon is assigned a weight: w γ,i = 8 π c N inj ξ T (cid:48) i Γ i R i δt δθ sin( θ i ) (1) For the simulations shown in this work, the injection was per-formed at a Thomson optical depth τ T ∼
100 or more. where ξ is the coefficient for the photon numberdensity ( n γ ) in the equation n γ = ξT ; ξ = 20 . ξ = 8 .
44 for a Wien spec-trum; T (cid:48) is the (comoving) fluid temperature atthe photon location; Γ i is the fluid’s bulk Lorentzfactor at the photon location; R is the injection ra-dius (in the Lab frame); δt is the time resolutionof the numerical simulation (time interval betweentwo saved frames); δθ is the polar angular range inwhich the photons are injected (in the Lab frame); θ i is the specific angle at which the i -th photon isinjected (in the Lab frame); and N inj is the numberof photons injected in each frame of the simulation.The subscript i is used for all quantities that arecomputed independently individually for each pho-ton. The weights are chosen so that the sum of theenergies of the injected photons times their weightgives the correct jet luminosity (see the Appendixfor a derivation of the photons weight). • Photon mean free path –
The first step inthe photons propagation is the calculation of theirmean free path in the expanding fluid. This is com-puted as (Abramowicz et al. 1991): λ i = drdτ T = 1 σ T n (cid:48) i Γ i (1 − β i cos θ fl ,i ) (2)where σ T is the Thomson cross-section, n (cid:48) i is the co-moving lepton number density of the fluid at eachphoton position, β i is the fluid’s velocity at the pho-ton position in units of the speed of light c , and θ fl ,i is the angle between the fluid and photon’s veloci-ties. Note that the mean free path should in prin-ciple be integrated over the evolving fluid densityand velocity. However, photons very rarely moveout of their original hydrodynamic resolution ele-ment within the time span δt between a scatteringand the subsequent one, and therefore our approxi-mation that the underlying hydrodynamic does notchange between the times at which the mean freepath is calculated is accurate. We also notice thatwe are assuming that the photon energies are wellbelow the electron’s rest mass energy (511 keV) andwe can ignore Klein-Nishina effects. In most cases,this is a good approximation for a GRB jet in thecomoving frame. • Selection of the interacting photon –
Once allthe mean free paths have been calculated, a randomscattering time for each photon is drawn from thedistributions: p i ( t ) = cλ i e − cλi t (3)The shortest of the obtained times is used as thenext scattering time. • Update of photon positions –
At this point thecode checks whether the selected collision time iswithin the time range of the current HD simulationframe. If so, all photons positions three-vectors R i are updated by propagating the photons at thespeed of light for the selected time interval. If not,a new time interval equal to the time remainingC Radiation Transfer in GRBs 3 Fig. 1.—
Block diagram of the main loop of the MCRaT code. -2 -1 Photon Energy (keV) F ( ν ) T=10 K; n coll = 7T=10 K; n coll = 36T=10 K; n coll = 157T=10 K; n coll = 7T=10 K; n coll = 36T=10 K; n coll = 157Wien Fit
Fig. 2.—
Spectra from a test simulation with a cylindrical outflowof constant velocity and temperature. Two simulations are shown,each at three different stages, for increasing number of scatteringsper photon. for the current HD simulation frame is selected,and all photons positions are updated. No scat-tering takes place in this case, a new HD frame isloaded, and the code returns to the mean free pathstep, computing a new mean free path for the pho-tons in the updated HD and then selecting a newscattering time. Note that, when a new simulationframe is loaded, a new batch of photons is addedat the injection radius following the rules of pho-ton injection as described above in bullet “photoninjection”. • Photon scattering –
When the time and locationof a photon scattering is identified, the code per-forms several steps. First, the photon four momen-tum P i is transformed to the local fluid comovingframe. Second, an electron is randomly drawn froma Maxwell-Boltzmann or Maxwell-J¨uttner distribu-tion at the local fluid temperature T i . The electronis assigned a direction of motion with a probabil-ity distribution p ( ψ ) ∝ sin ψ (1 − β cos ψ ), where ψ is the angle between the photon and electronvelocity vectors. Third, the electron and photonfour momenta are transformed to the electron’scomoving frame, and the reference frame rotatedso that the photon’s propagation is along the z-axis. A scattering angle is generated according to Time (seconds) L M C R a T / L P r e d i c t e d Expectation0.5 deg1.5 deg2.5 deg3.5 deg
Fig. 3.—
Comparison between the expected luminosity and theone computed with MCRaT for a spherical relativistic outflow.For the line of sight very close to the singular polar axis we find a ∼
20 per cent deviation between the expected and obtained results.All the other lines of sight are fully consistent with the predictedluminosity. the Thomson probability distribution p ( θ scatt ) ∝ sin θ scatt (1+cos θ scatt ) and the post-scattering fourmomenta are computed according to momentumand energy conservation. Finally, the photon’s fourmomentum is rotated and doubly boosted back tothe lab frame.The loop described above is performed iteratively untilthe photons reach the edge of the simulation, or the lastsimulation frame is reached (see also Figure 1). Eachtime a new simulation frame is loaded, the code savesthe photon position vectors and momenta four-vectors.In post processing, a script is used to generate a photonevent table that contains the photons detected by a vir-tual detector placed at a certain distance from the burstand along a selected orientation. The event table con-tains photon detection time and energy and the photonweight w γ,i . Code Validation
In order to validate the MCRaT code we ran a suiteof tests for which the expected result is known. All testswere run by imposing an analytic outflow solution on thesame simulation frames that will be used in Section 3 tostudy the light curves and spectra of a long GRB simu-lation. This means that the grid setup (cell centers andsizes) was read from the simulation files, but the valuesof density, pressure, and velocity were substituted withan analytical solution as a function of the cell position.This way, the same spatial resolution of the tests is usedfor the scientific simulation and we can evaluate the in-fluence of the grid size and resolution on the accuracy ofthe results.First we run a test in which the outflow is cylindrical.The velocity is constant and directed along the y axis,the temperature and the density are constant. We runseveral of these simulations with different outflow veloc-ity and/or different temperatures. Photons were injectedwith an initial Wien spectrum. Figure 2 shows the spec-tra for two outflow temperatures at various stages during Lazzati Photon Energy (keV) F ( ν ) Goodman (1986)Black BodyMCRaT
Fig. 4.—
MCRaT spectrum from a spherical outflow comparedwith a black body spectrum at the same temperature and with theprediction of Goodman (1986) for a spherical relativistic fireball.
Photospheric radius Decoupling radius
Fig. 5.—
Radial profile of the average photon energy from a testMCRaT simulation with a spherical flow. Red marks and coloredsolid lines show the radial evolution of the photon mean energy atthe observing angles θ o = 1, 2, and 3 degrees. A dashed line showsthe average energy of photons in equilibrium with the plasma. Thelocation of the photosphere and of the radius at which the radiationand plasma decouple according to Giannios (2012) are shown withvertical arrows. the simulation, i.e., after increasing number of scatter-ings. As expected, the resulting spectra are very well fitby a Wien spectrum (e.g. Chhotray & Lazzati 2015 andreferences therein), both at low temperature (leftmostspectra, non-relativistic electrons) and at high tempera-ture (rightmost spectra, trans-relativistic electrons). Thefigure shows the result from a static configuration with v = 0. Analogous results were obtained for outflows withrelativistic outward velocities (Γ = 10, 100, and 300 weretested).As a second test we run a suite of simulations of spheri-cal outflows, i.e., outflows with only outward radial veloc-ity component. In this case, the Lorentz factor increasesproportionally to the radius until all the internal energyis converted into bulk motion, saturation is reached, andthe fireball transitions from accelerating to coasting atconstant Lorentz factor. In the acceleration stage the electrons’ temperature decreases proportionally to r − ;as the asymptotic Lorentz factor is reached, the fireballcoasts at constant Γ and the electrons’ temperature de-creases as r − / . We stress that this analytical solutionwas imposed, we did not run a hydrodynamic simulationof a spherical outflow, we simply replaced the output ofthe simulation with the analytical prescription.We first compare the photon luminosity of the spheri-cal fireball calculated from the simulations with the an-alytical expectation. The result is shown in Figure 3,where the MCRaT luminosity for four different observersis compared to the analytical expectation. While for ob-servers away from the singular polar axis we find goodagreement, a deviation of ∼
20 per cent is observed at θ = 0 . ◦ . We ascribe this to the artifact of the polaraxis. Since the polar axis is singular in the coordinatesystem, the weight of photons close to the axis can bevery small (cfr. Eq. 1). For infinite number of photons,all small angles are sampled and there is no systematicdeviation in the results. For a finite number of photons,however, significant non-Poissonian deviations can resultfrom under- or oversampling the smallest polar angles.An analogous effect can be seen in Figure 8 of Morsonyet al. 2007 (dashed line) where an unphysical drop ofthe luminosity is observed at small polar angles from thejet axis. In the remainder of this manuscript we willcompute light curves and spectra only for photons prop-agating with a polar angle of at least 0 . ◦ to avoid thisproblematic region.In Figure 4 we show the time-integrated spectrum ac-cumulated between angles θ o = 1 and 2 degrees from thepolar axis . The spectrum is broader than a Planck func-tion, and agrees very well with the prediction of Good-man (1986). Finally, we show in Figure 5 the radial de-pendence of the average photon energy for various off-axis angles compared to the average frequency of pho-tons in equilibrium with the fluid. As described in greatdetail in Beloborodov (2011; see e.g. their Figure 3), thephoton temperature follows the electron temperature atlarge optical depth τ >
10. As the photosphere is ap-proached, the electrons and photons decouple and thephotons start to cool at a slower rate until, at the photo-sphere, the decoupling of the radiation is complete andthe photon temperature remains constant. The MCRaTresult is shown in Figure 5. Red symbols and coloredlines show the average photon energy for viewing angles θ o = 1, 2, and 3 ◦ . A dashed line shows instead the aver-age energy of photons strongly coupled to the electrons.Two vertical arrows show the location of the photosphereand the location at which the Giannios (2012) conditionis attained (the Wien radius, Beloborodov 2013). TheMCRaT result reproduces all the features of the analyt-ical solution (Beloborodov 2011) and confirms that theapproximations adopted so far for computing the peakenergy of the prompt GRB spectrum from simulationsare inaccurate. While it is true that electrons and pho-tons start decoupling at or around the Wien radius, the Note that since each photon propagates along a unique direc-tion, light curves and spectra cannot be calculated at a preciseviewing angle. Rather, photons within a small but non-zero rangeof angles need to be added. In the manuscript, we accumulatespectra and light curves from photons within an acceptance of 1degree, and use the average of these angles to indicate the directionof observation.
C Radiation Transfer in GRBs 5 × × × Distance from engine (cm) A v e r a g e p h o t o n e n e r g y ( K e v ) θ o =1 o θ o =2 o θ o =3 o θ o =4 o θ o =5 o Fig. 6.—
Radial profile of the average photon density for radia-tion propagating at different off-axis angles. Analogous to Figure 5but for the GRB jet simulation. Guidelines are plotted for off-axes1, 2, and 3 degrees, to show the deviation from the self-similarbehaviour as the photosphere is approached (or lack thereof). decoupling is incomplete and the spectrum keeps coolingby interaction with the electrons. This approximationwas adopted in L13 and as a consequence their resultsoverestimate both the light curve luminosity and thepeak frequency of the spectrum. On the other hand, theassumption that the electrons and photons are coupleduntil the photosphere is reached is also incorrect. Thisapproximation, adopted by (Lazzati et al. 2009, 2011;Mizuta et al. 2011; Nagakura et al. 2011) underestimatesthe radiation peak frequency and the light curve luminos-ity. A correct implementation of the radiation transfer innumerical simulations of the photospheric model is there-fore necessary not only to obtain the spectral propertiesbeyond the peak frequency (the low- and high- frequencyspectral indices), but also to correctly estimate the peakfrequency itself. Before showing the MCRaT results fora GRB simulation, we stress that within the viewing an-gles presented in this paper the code is quite accurate,the photon average energies differing only by a few percent for different lines of sight (Figure 5). RESULTS
We have applied the MCRaT code to our fiducial sim-ulation of a long duration gamma-ray burst jet that waspreviously studied in the approximation of a sharp pho-tosphere (Lazzati et al. 2009, 2011; L13). The simula-tion runs for ∼
600 s, at which time the jet head hasreached a distrance of ∼ . × cm from where it waslaunched, the center of the progenitor star. The progeni-tor is modeled as a 16 solar mass Wolf Rayet star, model16TI from Woosley & Heger (2006). The jet is injectedas a boundary condition at r = 10 cm, with luminosity L = 5 . × erg s − , half-opening angle θ = 10 ◦ ,Lorentz factor Γ = 5, and internal energy to allow foran asymptotic jet acceleration to a terminal Lorentz fac-tor Γ ∞ = 400. After t = 100 s, the engine luminosity isdecreased by a factor 1000 to simulate the turn off of theengine. The radial resolution at R = 10 cm, where thephotosphere is expected to lay, is δy = 10 cm, giving δy/y = 10 − , comparable to 1 / Γ . Resolution in the or-thogonal x coordinate is identical since we use square res-olution elements. Shocks at the photospheric radius are × × × × L i s o ( e r g / s ) MCRaT light curveL13 approximation (/10)MCRaT peak energy P e a k F r e q u e n c y ( k e V ) Time since jet launch (s) α βα β Fig. 7.—
Light curve and spectral evolution for photons travelingwithin the acceptance range 0 . ≤ θ ph ≤ . θ o = 1 ◦ ) from the polar axis. In the upper panel, the thicksolid line shows the MCRaT result, the thin line shows the resultfrom L13 (rescaled by a factor 10) and the symbols show the peakenergy in 1 second temporal bins (energy y-axis to the right). Inthe bottom panel, the Band parameters α and β are shown withred and blue symbols, respectively. The α scale is to the left, whilethe β scale is to the right. therefore poorly resolved. Still, this resolution is aboutan order of magnitude better than in the simulations ofIto et al. (2015). As we discuss below, limited resolu-tion is a serious concern for these results and an increaseby at least an order of magnitude is required to ensurethat shocks are resolved and that the simulations capturethe temperature variations that they produce. 300,000photons were injected in the simulation at off-axis angles θ i ≤ ◦ .Under conical expansion, the jet we input in the simu-lations would become optically thin (and therefore haveits photosphere) at R = 6 . × cm (e.g., M´esz´aros& Rees 2000). However, the progenitor-jet interactioninduces a non-conical expansion and prevents the satu-ration of the acceleration at Γ = 400, resulting in a muchlarger photospheric radius. Lazzati et al. (2011) founda photospheric radius of the order of R = 10 cm inthe same simulation, increasing for larger off-axis angles.With MCRaT we do not compute the photospheric ra-dius by integrating the optical depth (as is done in allprevious numerical simulations of photospheric emissionin GRBS); we follow instead the photon propagation.Whether the photons have crossed the photospheric ra-dius or not can be understood by looking at the radialdependence of their average energy. The average photonenergy decreases at radii smaller than the photosphereand attains an asymptotic value beyond it (Beloborodov2011; see also Section 2.1). In Figure 6 we show thetime-averaged MCRaT results for off-axis angles of 1, 2,3, 4, and 5 degrees, each with an acceptance of ± . R ph ∼ . × cm for small off-axis angles Lazzati × × × × L i s o ( e r g / s ) MCRaT light curveMCRaT peak energy P e a k F r e q u e n c y ( k e V ) Time since jet launch (s) α βα β Fig. 8.—
Same as Figure 7 but for an acceptance angle 1 . ≤ θ ph ≤ .
5. L13 results at this angle are not available for compari-son. but is located beyond the outer edge of the simulationbox (which is located at 2 . × cm) for relativelylarge off axis angles θ o > ◦ . In mode detail, the photo-sphere is clearly reached for θ o = 1 ◦ , is almost reachedfor θ o = 2 ◦ , and it is approached for θ o = 3 ◦ . In thefollowing we show results for these three off-axis angles,keeping in mind that the result for θ o = 3 ◦ might beapproximate since the simulations do not appear to beextended enough to completely reach the photosphere.It should be remembered, however, that the photosphereis not at a constant distance throughout the simulation,starting at a large radius and progressively closing in(see, e.g., the blue lines in Figure 4 of L13). Even for thelarger off-axis angles, the photosphere is likely reached atsome time after the initial high-density plug at the jet’shead is overcome. Light Curves and Spectra
Light curves and spectra are computed by binning thephotons within the acceptance angle in time and/or fre-quency, each multiplied times its weight (Eq. 1). Lightcurves for 0 . ≤ θ ph ≤ .
5, 1 . ≤ θ ph ≤ .
5, and2 . ≤ θ ph ≤ . . ≤ θ ph ≤ . × × × × L i s o ( e r g / s ) MCRaT light curveL13 approximation (/4)MCRaT peak energy P e a k F r e q u e n c y ( k e V ) Time since jet launch (s) α βα β Fig. 9.—
Same as Figure 7 but for an acceptance angle 2 . ≤ θ ph ≤ . addition, as shown in Figure 6, the photospheric radiusis contained in our simulations at this off-axis angle. Thethree panels show the photon count spectra for the wholeburst (left panel), a fairly non-thermal time interval (cen-tral panel), and the most luminous time bin of the lightcurve (left panel). In all cases we notice some features inagreement with observed gamma-ray bursts and othersin significant tension with the observations. All spectraare fit with a Band model (Band et al. 1993) and thebest values for the parameters α , β , and hν pk are re-ported. In the Band model, a smoothly joined brokenpower-law function, α is the asymptotic photon index atlow frequency, β is the asymptotic photon index at highfrequency, and hν pk is the photon energy at which the νF ( ν ) spectrum peaks.Going from low- to high-frequencies, we first noticethat the low-energy spectral index α is found to be sig-nificantly larger than in observations. Average observedvalues of α range between -1 and -0.5 (Preece et al. 2000;Goldstein et al. 2013; Gruber et al. 2014; Yu et al. 2016)while we consistently find α >
1. This is not surpris-ing since the difficulty of the photospheric model to re-produce the observations in the low-frequency range iswell known. Our synthetic spectra have power-law in-dices even steeper than a black body, due to the factthat the spectrum was injected with a Wien distributionof frequencies . Injecting a black body spectrum wouldhave made the low-frequency spectrum harder, howeverthe steepening of the spectrum would not have bee sub-stantial enough to reproduce the typical spectra of GRBobservations. We also notice that our simulations do nottest any of the solutions proposed to solve the discrep-ancy. For example, models based on sub photosphericdissipation (e.g. Chhotray & Lazzati 2015) can be testedonly with simulations of higher resolution than the onepresented here, since shocks need to be fully resolvedto affect the radiation spectrum. Models based on theaddition of soft photons from synchrotron (Vurm et al.2011; Vurm & Beloborodov 2015) require instead a radi-ation transfer code that allows for photon emission and The low-frequency photons require more scattering to reachequilibrium and our vary soft tail may also be affected by an in-sufficient number of scattrings (Chhotray & Lazzati 2015).
C Radiation Transfer in GRBs 7absorption, beyond the current capabilities of MCRaT.The peak energy of the light curves is found in the ∼
100 keV range. This is a fairly common value for aburst of intermediate luminosity like the one simulatedhere. Yet, it is on the low tail of the distribution and wewill see in Section 3.2 that there is some strain betweenour results and the Amati and Yonetoku correlations.The best agreement between our results and observa-tions is found at high frequencies. We first notice that wecan exclude both a cutoff power-law and a Black Body fitwiths high confidence. A χ -based F-test yields probabil-ities P < − for the two simpler functions mentionedabove to give a comparable fit with respect to a BandFunction. A high-frequency power-law is therefore re-quired by the data, as in observed bursts. In addition,we find that MCRaT’s β values, between 2 and 3, arefairly typical for long-duration GRBs (Preece et al. 2000;Goldstein et al. 2013; Gruber et al. 2014; Yu et al. 2016).We note that the spectra shown only extend to ∼ M eV due to poor statistics, but there is no indication of a cut-off. In the left panel – the one with the best statistics –we show with smaller symbols and in grey color the datain the spectral bins that have too little photons to beincluded in the fit. We notice that these data points areconsistent with the extrapolation of the high-frequencypower-law and no break or cutoff is required. Very highfrequency photons in the GeV band detected by Fermi(Abdo et al. 2009) are however better explained as for-ward shock emission in the photospheric model (Kumar& Barniol Duran 2009; Ghisellini et al. 2010).Figures 7, 8, and 9 show, besides the light curves, theresults of time-resolved spectroscopy. Photons were ac-cumulated in 1 s bins and a Band spectrum was fit to eachtemporal bin and for each of the three off-axis observers.The top panel of each figure reports the peak energy inkeV (right hand vertical scale), while the bottom panelreports the best fit α and β values. The low-frequency in-dex α refers to the left hand vertical scale, while the high-frequency index β refers to the right hand vertical scale.We first notice that there is quite some diversity in thekind of spectral evolution observed in the three figures.At small off-axis θ o = 1 ◦ , the peak frequency follows ageneral hard-to-soft trend, while the low-frequency index α has a small growing trend. Eventually, both trends arebroken. At times larger than t ∼
40 s the peak energyrises above 100 keV while the spectral index α fluctu-ates with no stable trend. This transition is due to theemergence of the unshocked jet, as already seen in previ-ous simulations (Morsony et al. 2007; Lazzati et al. 2009;L13). The hard-to-soft trend has been observed in somebursts. The α increase, on the other hand is not ob-served but seem irrelevant since, as already noted above,our spectra have α indices in significant disagreementwith observations.Moving away form the axis, at θ o = 2 ◦ we find a lightcurve with very little spectral evolution, except for aninitial hardening and a small softening around the peak.This is quite unusual in observations. Finally, at θ o = 3 ◦ we find an example of tracking behavior, the peak en-ergy increasing at high luminosity and decreasing at lowluminosity. Again, towards the end of the burst a globalhardening is observed, possibly due to the emergence ofthe unshocked jet or to the fact that photons at the endof the burst haven’t yet reached the photosphere. In all three cases, non-thermal high-energy tails are measuredthroughout the burst, with the exception of a few in-tervals in which the high-frequency index was fixed to β = − Correlations
One of the reasons for the success of the photosphericmodel was its ability to reproduce the Amati, Yonetoku,Golenetskii, E iso − Γ, and energy-efficiency correlations(Golenetskii et al. 1983; Amati et al. 2002; Yonetoku etal. 2004; Liang et al. 2010; Ghirlanda et al. 2012). Thiswas demonstrated with the results of numerical simula-tions, and the scaling was ascribed principally to a line ofsight effect (L13). With MCRaT we can re-analyze thesuccess of the photospheric model in reproducing the ob-servational ensemble correlations. We show in Figure 11the result for the Yonetoku correlation between the 1-sbinned peak luminosity and the peak frequency measuredat the time of the peak (again, integrated over 1 s). Giventhe results shown above, it does not come as a surprisethat the MCRaT result are not a perfect match to theobservations. The assumption of a complete decouplingat the Wien radius adopted by L13 (see also Giannios2012) is indeed fairly crude. While it is true that the ra-diation decouples from the leptons at that optical depth,the radiation continues to cool, even though remainingat a temperature that is higher than that of the matter(Beloborodov 2011; see also Figure 5). As a consequence,the MCRaT spectra peak at a frequency that is too lowby a factor ∼ E iso − Γ and energy-efficiency cor-relations, for which no significant discrepancy is found.The agreement, however, is most likely due to the poordefinition of the correlations in the data rather than toa true success of the model. SUMMARY AND DISCUSSION
We have presented a novel numerical code to performMonte Carlo radiation transfer in relativistic jets underthe assumption that the radiation-matter interaction isdominated by Compton scattering off the leptons in thefluid. Our code is analogous to the one recently devel-oped by Ito et al. (2013, 2015) and we find similarresults, even though they applied it to a different set ofsimulations. Neither MCRaT nor the Ito et al. (2013,2015) code cannot take into account the radiation feed-back on the hydrodynamics. However, the hydrodynamic The difference between the two codes is technical (language,dimensionality) but the physics is the same.
Lazzati hν (keV) N ( ν ) ( p h o t o n s / s / k e V ) α =1 . β = − . hν pk =78 keVEntire light curve hν (keV) α =1 . β = − . hν pk =87 keV1s bin at 13s hν (keV) α =1 . β = − hν pk =66 keV1s bin at peak Fig. 10.—
Sample MCRaT photon-count spectra for the 1 . ≤ θ ph ≤ . t = 13 s. Finally, the left panel shows the spectrum at the light curve peak. Each panel shows the best fit Band parameters inan inset and the best fit Band spectrum as a thin solid line. In the left panel bins with some photon counts but not enough statistics tobe included in te fit are shown in gray and with a smaller size. simulations to which the Monte Carlo codes are appliedassume infinite coupling (the equation of state for bothphotons and plasma are the same) and therefore the dy-namics from the HD simulations should be accurate inthe optically thick part of the outflow, where most of theCompton scatterings take place and the spectrum takesits shape. Despite these limitations, both MCRaT andthe Ito et al. code constitute a significant advance overprevious studies, in which the coupling between matterand radiation was assumed to be complete at opticaldepths larger than a critical value and null at smalleroptical depths. Said critical optical depth was assumedto be either unity (Lazzati et al. 2009, 2011; Mizuta etal. 2011; Nagakura et al. 2011), or of several tens (Gian-nios 2012; L13). As analytically predicted for sphericaloutflows by Beloborodov (2011), proper radiation trans-fer shows that the radiation and matter indeed start todecouple at an opacity of several tens. However, the ra-diation keeps cooling significantly until the photosphereis reached, at which point the radiation propagates withconstant average properties . In addition, MC radiationtransfer allows for a full calculation of the emitted spec-trum, while only its peak frequency could be estimatedunder the approximations previously employed.We have applied MCRaT to our fiducial GRB jet simu-lation, used in previous photospheric work. We find thatthe MCRaT results predict light curves that are dim-mer and characterized by lower frequency photons withrespect to the results from L13. Spectral and tempo-ral analysis was performed on the resulting light curvesand spectra. We find that, as predicted by many ana-lytical and semi-analytical previous publications, photo-spheric radiation is characterized by a significant non-thermal high-frequency tail. Fitting the time-resolvedspectra with a Band function we find values of the index β in good agreement with the observations and a signif-icant evolution of the slope over the burst light curve.Because of the continued radiation cooling in the sub- It is envisageable that in the presence of a significant dissipa-tion event, such as an external shock, the radiation properties canbe changed well after the photosphere is crossed. Such events are,however, not present in our simulations. Isotropic equicalent luminosity (erg/s) P e a k e n e r g y ( k e V ) Fig. 11.—
Comparison between the observational Yonetoku cor-relation (black symbols with error bars and best fit solid line) withthe MCRaT results (red symbols). photospheric region, we find that the MCRaT spectraare characterized by photons of lower frequency with re-spect to the L13 results. This deteriorates the agreementof the simulation results with the Amati and Yonetokuensemble correlations. While the MCRaT results are stillwithin the observed dispersion, they do not follow thecorrelation closely and more work is necessary to bet-ter understand the role of photospheric emission in theestablishment of such correlations. In addition, we findthat our simulations do not reach the photospheric radiusfor off-axis angles θ o > ◦ , preventing the analysis of theoff-axis light curves that were deemed to be responsiblefor the correlations (L13). We also find that the low-frequency spectral index is much steeper than observed.This is a well known problem of photospheric emissionin the absence of an efficient emission mechanism or ofstrong sub-photospheric dissipation and the discrepancyfound is not surprising.One possible shortcoming of this work that could beat the basis of the discrepancies discussed above is theresolution of the simulations. In order to study the ef-fect of shock-heating on the spectrum we need to be ableC Radiation Transfer in GRBs 9to resolve the shocks so that their true temperature canbe taken into account. Unfortunately, the resolution ofthe simulations presented here is sufficient only to mod-erately resolve the shocks. As a consequence, continuousheating from recombination, or the kind of shock heatingenvisaged in, e.g., Chhotray and Lazzati (2015), couldnot be studied. We plan to perform higher-resolutionsimulations for both long and short GRB jets in the nearfuture to test the impact of resolution on the results.Some hope comes from the comparison of our results withthose of Ito et al. (2015). Their resolution is lower thanours, but they consider precessing jets that substantiallyincrease dissipation. Even though they do not perform adetailed time-resolved analysis of their spectra, it is clearfrom their Figures 3 and 4 that significant non-thermalfeatures at both ends of the spectrum can be obtained.In addition, we stress that the simulation presented hereinjects a constant luminosity from the inner engine anddoes not have, therefore, internal shocks. Realistic en-gines are likely characterized by intrinsic variability, anda better agreement with the Amati and Yonetoku en-semble correlations could be obtained from a light curvewith the same peak frequency and peak luminosity butwith periods of very weak emission.Additional physics might also be required to obtain abetter match of the synthetic spectra with observations.Shocks likely inject non-thermal particles that are ne-glected by MCRaT and by the Ito et al. (2015) codealike. Magnetic fields and the possibility of injection ofsoft photons from synchrotron are also ignored (Vurm etal. 2011; Vurm & Beloborodov 2015). Additional dissi-pation is also provided by magnetic reconnection (Gian-nios & Spruit 2007; McKinney & Uzdensky 2012). Whileadding such effects would not be trivial, breaking the em-barrassingly parallel structure of the code and makingthe whole calculation longer, it is not impossible and we are currently developing a version of the MCRaT codethat takes into account non-thermal particles and syn-chrotron emission. Since all the missing effects are likelyto increase the non-thermal aspect of the spectrum, ourthe results presented here can be seen as a lower limit onthe non-thermal features of the photospheric spectrumfrom GRB jets.Finally, we would like to comment on the fact that inno case we found the need to add a low-temperature ther-mal component to the spectral model used to fit our spec-tra. Yet, such a component has been detected in severalbursts and is considered one of the proofs of the existenceof a photospheric contribution to the GRB prompt spec-trum (Guiriec et al. 2011; Ghirlanda et al. 2013; Guiriecet al. 2013; Burgess et al. 2014; Nappo et al. 2016). Ourresults, as well as Ito’s results seem to show that the ap-pearance of the photospheric component is always non-thermal and at significantly higher frequencies than thoseat which the thermal components are detected (typicallya few keV to a few tens of keV). The thermal peaks de-tected in BATSE and Fermi observations are thereforeeither arising from a different component of the outflowthat is not present in the simulations or are a signatureof a weak photosphere in a jet dominated by non-thermaland/or non-baryonic energy, such as in a Poynting dom-inated outflow.I would like to thank the anonymous referee for con-structive comments that improved the clarity of this pa-per. I would like to thank Hirotaka Ito and Hiro Na-gataki for useful discussions and for sharing the resultsof their code with me prior to publication. I also wouldlike to thank Andrei Beloborodov, Riccardo Ciolfi, BrunoGiacomazzo, Dimitrios Giannios, Diego L´opez-C´amara,and Brian Morsony, for useful comments and suggestions.This work was supported in part by NASA Swift GI grantNNX15AG96G. APPENDIX
DERIVATION OF THE PHOTON WEIGHT
The luminosity of a uniform jet can be written as (e.g., Morsony et al. 2010): L j = Ω j R (cid:0) p + ρ (cid:48) c (cid:1) Γ c (A1)where p is the pressure of the jet material, and ρ (cid:48) its comoving density. For a jet dominated by radiation ( p = aT (cid:48) / p (cid:29) ρ (cid:48) c ), and considering the solid angle of a one-sided jet, the equation simplifies to: L j = 8 π − cos θ j ) R aT (cid:48) Γ c (A2)Deriving with respect to the polar angle we find the luminosity per unit polar angle (the simulated quantity incylindrical symmetry) and dividing by the average photon energy h ¯ ν = Γ h ¯ ν (cid:48) we obtain the number of photons perunit polar angle as: dN γ dt dθ = 1 h ¯ ν dLdθ = 8 π θ j R aT (cid:48) h ¯ ν (cid:48) Γ c = 8 π θ j R n (cid:48) γ Γ c = 8 π θ j R ξT (cid:48) Γ c (A3)This should be equal to the number of photons we inject in the simulation times the weight assigned to each photon: dN γ dt dθ = w γ,i N inj δt δθ (A4)Solving for w γ,i yields Eq. 1. REFERENCESAbdo, A. A., Ackermann, M., Arimoto, M., et al. 2009, Science,323, 1688 Abramowicz, M. A., Novikov, I. D., & Paczynski, B. 1991, ApJ,369, 1750 Lazzati