The Role of a Neutron Component in the Photospheric Emission of Long-duration Gamma-Ray Burst Jets

Long-duration gamma-ray bursts (LGRBs), thought to be produced during core-collapse supernovae, may have a prominent neutron component in the outflow material. If present, neutrons can change how photons scatter in the outflow by reducing its opacity, thereby allowing the photons to decouple sooner than if there were no neutrons present. Understanding the details of this process could therefore allow us to probe the central engine of LGRBs, which is otherwise hidden. Here, we present results of the photospheric emission from an LGRB jet, using a combination of relativistic hydrodynamic simulations and radiative transfer postprocessing using Monte Carlo radiation transfer code. We control the size of the neutron component in the jet material by varying the equilibrium electron fraction Y e , and we find that the presence of neutrons in the GRB fireball affects the Band parameters α and E 0, while the picture with the β parameter is less clear. In particular, the break energy E 0 is shifted to higher energies. Additionally, we find that increasing the size of the neutron component also increases the total radiated energy of the outflow across multiple viewing angles. Our results not only shed light on LGRBs but are also relevant to short-duration gamma-ray bursts associated with binary neutron star mergers due to the likelihood of a prominent neutron component in such systems.


INTRODUCTION
Our understanding of Gamma-Ray Bursts (GRBs) has evolved dramatically since their discovery in the late 1960's.First detected as short transient bursts of high energy photons (Klebesadel et al. 1973), observations of afterglows (Groot et al. 1998;Costa et al. 1997) and supernova counterparts (Galama et al. 1998;Hjorth et al. 2003;Woosley & Bloom 2006;Bloom et al. 1999) have facilitated a deeper understanding of these otherwise mysterious events.Long duration gamma-ray bursts (LGRBs) are now thought to occur during core-collapse supernovae, a process in which stars more massive than about 8M ⊙ end their lives in a violent explosion, resulting in the formation of either a Black Hole (BH) or a Neutron Star (NS) (Woosley & Janka 2005).After the formation of either a BH or a NS, material from the preceding collapse can accrete around the compact object, providing a possible power source for an ensuing LGRB (e.g.Narayan et al. (2001)).Alternatively, a highly magnetized, fast spinning NS could power a rel-ativistic outflow by tapping into its rotational energy (e.g., Bucciantini et al. 2012).Given the possibility of a NS as either an intermediate or a terminal stage of the supernova, there is a strong possibility of a neutron component in the accreting material, which can then be collimated into a relativistic jet and produce a LGRB.
In spite of this progress, one aspect of GRBs that still remains in contention is the nature of the prompt emission.In LGRBs, the prompt emission can last anywhere from a few seconds to a few minutes (Bloom et al. 1999;MacFadyen et al. 2001) and is characterized by bright, non-thermal spectra (Band et al. 1993).A leading model that explains this emission is the Synchrotron Shock Model (SSM).In this model, the jet expands and reaches the photosphere without producing noticeable radiation.After passing the photosphere, electrons in colliding internal shocks produce non-thermal radiation (Rees & Meszaros 1994).While this model naturally explains the characteristic non-thermal emission of GRBs and is able to fit the spectra of a number of bursts, it may have difficulties in reproducing the peak width of bursts (Beloborodov 2013).In addition, some burst have spectra that are inconsistent with a simple model in which electrons are accelerated impulsively and either do not cool (the line-of-death problem, Preece et al. 1998) or cool radiatively (Ghisellini et al. 2000).Finally, the SSM model has difficulty reproducing the ensamble correlations between properties of different bursts, such as the Amati and the Yonetoku correlations (Amati et al. 2002;Yonetoku et al. 2004).
A viable alternative to the SSM is the so-called photospheric model (e.g.Beloborodov (2010a), Giannios & Spruit (2007), Lazzati et al. (2009), Ryde et al. (2011), Pe'er et al. (2006)).In this model, thermal radiation is produced when the jet is hot and dense near the central engine.As the jet propagates and expands the radiation is shaped through its interaction with the expanding outflow.Effects such as sub-photospheric dissipation (Chhotray & Lazzati 2015;Parsotan et al. 2018;Ito et al. 2018) and multi-color blackbody emission (Pe'er & Ryde 2011) enable this model to account for a non-thermal spectrum.Additionally, as the radiation scatters and propagates with the outflow, it is imprinted with a signature of the history of the outflow that survives until the radiation escapes at the photosphere (Vurm & Beloborodov 2016).Because of this, the composition and dynamics of the jet material are of crucial importance in shaping the observed prompt emission in the photospheric model.An important test of the photospheric model can then be to model the effect that different compositions of the jet material can have on radiation.
Given the possibility of a neutron component in both the compact mergers and supernovae that are thought to produce GRBs, a body of work has been produced that explores the consequences of a neutron component in GRB fireballs.This includes a detailed study on the processes that shape the nuclear composition of the fireball as it expands (Beloborodov 2003), the role neutrons play in heating the jet through collisional processes (Beloborodov 2010b; Rossi et al. 2004), and that of the dynamics of shocks in the explosion (Derishev et al. 1999).However, no work has been done on how neutrons directly shape the observed prompt emission of GRBs.Therefore the role of a neutron component on the photosepheric emission of a LGRB is of particular interest, and a good candidate to further test the photospheric emission model.
In this paper, we use the MCRaT radiative transfer code and the ProcessMCRaT python package (Lazzati 2016;Parsotan & Lazzati 2018;Parsotan et al. 2018;Parsotan & Lazzati 2021) to scatter photons through a 2D relativistic hydrodynamic (RHD) simulation of a LGRB jet (Morsony et al. 2007;Lazzati et al. 2013), and produce mock observables.We control the relative size of the neutron component in the jet material by varying the equilibrium proton-to-nucleon ration Y e .This paper is organized as follows: in Section 2 we summarize how the MCRaT code scatters photons and describe how we take into account a neutron component in the jet; in Section 3 we present results of spectra obtained by varying Y e ; and in Section 4 we discuss our results and their implications.

The MCRaT Code
We use the MCRaT radiative transfer code to individually Compton-scatter a set of photons injected into a RHD simulation of a LGRB jet.In this section we summarize the MCRaT algorithm.Further details on the original algorithm can be found in Lazzati (2016), with improvements found in Parsotan & Lazzati (2018).
MCRaT begins by injecting photons into the output of a RHD simulation.During this injection process, MCRaT selects which RHD cells to inject photons into based on a set of user-defined parameters: the injection radius R inj and the angular interval δθ, defined with respect to the jet axis.All cells within the interval δθ and with a radius between R inj ± 1 2 cδt are selected, where c is the speed of light and δt is the time interval of the selected RHD frame.Once the injection frames are selected, MCRaT determines the four-momentum of the injected photons in each cell by sampling a thermal distribution centered at the local co-moving temperature, where p i is the pressure of the fluid and a is the radiation density constant.The injected photons are then weighted according to (Parsotan et al. 2018) where dN i is the expected number of photons in the i th RHD cell, ξ is the photon number density coefficient from n γ = ξT 3 (ξ = 20.29 for a Planck spectrum and ξ = 8.44 for a Wein spectrum), T ′ i is the comoving fluid temperature, Γ i is the bulk Lorentz factor, dV i is the volume element of the RHD cell, and w is the weight of the injected photons.MCRaT calculates the expected number of photons in each cell via Equation 2, and draws a photon number from a Poisson distribution with a mean given by the expected number of photons.MCRaT then sums over the photon numbers in each cell, and if the total number of injected photons so obtained lies outside the user-specified range, the weights are adjusted and the process of calculating the expected number of photons via Equation 2 repeats.The final weights are those that result in a total number of injected photons that lies within the user-defined range.
Once the injected photon properties are determined, MCRaT scatters each photon according to the properties of the RHD simulation.To begin with, each photon is assigned a mean-free path according to Abramowicz et al. (1991) where σ T is the Thomson cross section, n ′ i is the comoving lepton number density, β i is the fluid velocity in units of c, and θ f l,i is the angle between the fluid velocity and the photon velocity.A random scattering time for each photon is drawn from the distribution and if the smallest of those scattering times is within the time interval of the given hydrodynamical simulation frame, the positions of the photons are all updated by allowing them to travel at the speed of light for the smallest scattering time obtained via Equation 4. Once all photons are updated to a new position in a frame, the photon with the shortest scattering time is scattered with an electron drawn from either a Maxwell-Boltzmann or a Maxwell-Jüttner distribution at the local fluid temperature, with a direction drawn from a random distribution.If the smallest scattering times obtained from equation 4 lies outside the given RHD frame time interval, MCRaT allows the photons to propagate at the speed of light, without scattering, for an amount of time equal to the remaining time in the current RHD frame.Then, MCRaT loads a new simulation frame and the photon mean free paths are all calculated again.This process of calculating photon properties and scattering with electrons is repeated for all the injected photons as they propagate and scatter through all of the provided RHD simulation frames.

Mock Observations
When all the injected photons have been diffused beyond the photosphere we use the ProcessMCRaT package (Parsotan & Lazzati 2022) to conduct mock observations.This software allows for the injected photons to continue propagating unimpeded out to a virtual detector placed at a user-defined radius.To mimic a real observation in which the viewing geometry is unique we count only photons within a given acceptance range around the angle to the observer.The energies of photons are obtained from the time component of the fourmomentum at the end of the simulation, and the detection time is calculated as where t real is the simulation time at the frame used for an observation, t p is the photon detection time, and t j = r d /c is the time it takes for a photon that was emitted at the instant the jet was launched to propagate to the detector.
In the following, all light curves and spectra are obtained by placing the virtual detector at a radius of r d = 2.5 × 10 13 cm, which corresponds to approximately the edge of the RHD simulation.When the photons haven't yet reached the last frame we find the positions of all the photons at the corresponding RHD simulation time and place a detector slightly beyond that point.
After conducting a mock observation, we can bin the photon arrival times to calculate light curves, where E i is the energy of each photon, ∆t bin is the time bin, and is the solid angle the detector occupies, with ∆θ being the angular acceptance range centered around θ v .We also bin the photon energies to calculate spectra via where all the terms are the same as in Equation 6except ∆E bin and ∆t, which are the energy bin width and the time interval over which photons were detected, respectively.
We fit the Band function (Band et al. 1993), to spectra obtained from equation 7.In equation 8, α and β are the low and high energy slopes, respectively, E 0 is the break energy, and A is related to normalization.
The spectral peak is defined with respect to the spectral parameters in equation 8 as In order for the calculated spectra and light curves to correspond to what an observer would see, the optical depth must reach a value τ ∼ 1.We calculate the optical depth (Parsotan et al. 2018) as: where L is the last frame of the RHD simulation and n refers to a group of photons located initially in the i th frame, at some average position R i .The sum over the RHD frame number j goes from the i th frame to the last, with ⟨N ⟩ n j being the average number of scatterings that the n th group of photons experienced in the j th frame.Equation ( 9) essentially counts the number of scatterings each photon undergoes starting from the i th frame and we calculate it by tracking the number of scatterings that individual photons undergo, starting immediately after they are injected.We similarly calculate the average energy of individual photons by tracking their energy throughout the MCRaT simulation.
A group of photons is uncoupled from the jet if the average number of scatterings per photon starting from the i th RHD frame is ≲ 1, corresponding to the photosphere condition of τ ∼ 1.Since this is computed separately for separate groups of photons it allows for the fact that photons in different parts of the jet and cocoon may uncouple at different times.

A Neutron Component in the Fireball
The MCRaT code reads in hydrodynamical data and determines the energy of injected photons via the hydrodynamical pressure (Equation 1), and their mean free paths via the hydrodynamical density and velocity (Equation 3).Normally it is assumed that the total mass of the hydrodynamical simulation is attributed to protons (with a negligible contribution by electrons), and the lepton number density is therefore calculated by dividing the hydrodynamical density by the mass of the proton.This picture changes when we include neutrons in our radiative transfer simulations.
To simulate the role of a neutron component in the fireball, we use the proton-to-nucleon ratio, Y e , defined through the charge neutrality condition (Beloborodov 2003) n where n ± are the e ± number densities.In the absence of e ± pairs, Y e is just the electron-to-nucleon ratio and describes how many electrons there must be in a plasma in order to preserve charge neutrality.The density ρ in Equation 10 can in general consist of both protons and neutrons, and when both are taken into account the result is the equilibrium electron fraction Y e = n p /(n p + n n ).Therefore, increasing the fraction of neutrons in the fireball decreases the electron-to-nucleon ratio, which in turn leaves fewer electrons with which to scatter photons.When calculating photon mean free paths via Equation 3, we can then scale the lepton density by Y e .A larger neutron component reduces the lepton density of the jet.
A neutron component can in principle also change the hydrodynamical behavior of the plasma.When the jet is still near the central engine it is hot and dense enough that the charged current reactions, establish an equilibrium Y e .While these conditions will change as the jet expands, it has been shown that, further from the central engine, neutrons and ions can stay coupled through the acceleration stage as long as the jet has relatively high baryon loading (Beloborodov 2003).
In the same work it was also found that fireballs from neutron rich central engines are likely to remain neutron rich.We therefore do not consider the hydrodynamical effects of neutrons decoupling from protons, and we likewise keep the value of Y e constant throughout our MCRaT simulations.Since the baryons are treated as being in equilibrium we leave the pressure and velocity variables from the RHD simulation unchanged, and we scale the fluid density by the equilibrium electron fraction Y e : ρ → Y e ρ.
While we use a constant value of Y e for each MCRaT simulation we run, the RHD simulation is in fact comprised of material ejected from the central engine, a stellar envelope through which the jet must escape, and a radial power law as the jet propagates into the interstellar medium.All of this materials could, in principle, have a different composition.In light of this, a constant value of Y e applied to the entire RHD domain is just an approximation.To ensure that such approximation gives reliable results, we restrict this study to the region near the jet axis by injecting photons only within the first 3 • relative to the jet axis, where the jet material has a high temperature and Lorentz factor.The role of mixing between materials with different Y e will be explored in a future work.

RESULTS
In this paper we used the FLASH version 2.5 2D RHD simulation in Lazzati et al. (2013) that is based on a 16TI progenitor (Woosley & Heger 2006) in which a jet with initial Lorentz factor of 5 and an opening angle of 10°is injected into the 16TI progenitor for 100 s and propagates out to the photosphere at ∼ 10 13 cm.The 16TI simulation in Lazzati et al. (2013) was performed on an adaptive mesh grid with a maximum resolution of 4 × 10 6 cm and output files were saved every δt = 0.2 s.For our MCRaT simulations to converge according to Arita-Escalante et al. ( 2023), injected photons should travel through multiple RHD cells in each frame.This can be quantified through the light crossing ratio, defined as cδt/δr which, with the spatial and temporal resolutions from the 16TI simulation used here, results in a light crossing ratio as large as ∼ 1500.
Our methods are similar to Parsotan et al. (2018), with a key difference being that we inject ∼ 2 × 10 5 photons for ∼ 50 s of a non-variable jet, which excludes only a constant, low luminosity portion of the lightcurve that is not observed in nature.We also restrict photon injection to the first 3°of the jet as outlined in Section 2. We then adopt a viewing angle of θ v < 3°when conducting mock observations.For the electron-to-nucleon ratio we use the values Y e = 1, 0.7, 0.4, and 0.1 to cover the cases of a small to large neutron component.
Figure 1 shows lightcurves obtained at a viewing angle of θ v = 1 • alongside the time-resolved best fit parameters α and β for the Band function (equation 8), in addition to the peak energy E pk = (2 + α) E 0 , for all 4 values of Y e .Our lightcurve from the Y e = 1 simulation agrees well with past MCRaT results based on similar 16TI simulations (Lazzati 2016;Parsotan & Lazzati 2021), and all lightcurves show a characteristic small peak at ∼ 8 s, with a brighter peak at ∼ 30 s.As Y e is increased, the second peak dims noticeably as evident in panel (d) of Figure 1, where the first peak is brighter than the second.
In Figure 2 we show time-integrated spectra obtained from photons in the Y e = 0.1 and Y e = 1 MCRaT simulations that have reached the final RHD frame.Both spectra in figure 2 were integrated from 0 to 40 s, corresponding to the first two peaks seen in figure 1.As with Figure 1, our spectra with Y e = 1 agrees well with past results.Here, as Y e is decreased, the peak energy shifts to higher frequencies as seen by the dotted lines in Figure 2. We will look at how Y e affects other spectral parameters below.
Figure 3 shows a corner plot for a Band function fit to the Y e = 0.1 spectrum.While spectral parameters in figures 1 and 2 where obtained via a non-linear least squares fitting algorigthm available in ProcessM-CRaT, those in Figure 3 were obtained by fitting a Band function to our MCRaT data with a Markov Chain Monte Carlo algorithm via emcee (Foreman-Mackey et al. 2013).The parameters in Figure 3 are different from those seen in Figure 2 due to the different methodologies used to obtain them.Figure 3 shows a clear correlation between E 0 and α, while the other pairs of parameters have no notable correlations.This strong correlation between α and E 0 plays a part in the evolution of Band function parameters for all four of the MCRaT simulations in this work.
It is also illuminating to analyze the behaviour of the Band Function parameters as the radiation propagates with and through the outflow material.We do this by conducting a mock observation and calculating spectra for multiple intermediate times throughout the MCRaT simulation.At each of these times, the injected photons have scattered through only a portion of the RHD simulation, and thus have some average distance from the central engine.This distance increases as the photons propagate with the outflow until they near the photosphere.For these observations, the position of the detector is determined by the positions of the photons at a given frame.The Band function is fit to the spectrum at each time, and Figure 4  Figure 5 shows the average photon energy as a function of their distance from the central engine.Figure 6 similarly shows how the optical depth (equation 9) of the injected photons.In figures 5 and 6, the photon energy and number of scatterings for each photon are, respectively, calculated for every individual photons starting immediately when they're injected near the central engine.
As stated in the Methods section, for the spectra and lightcurves from MCRaT to correspond what an observer would see from an actual burst, the photons have to decouple from the jet material.Figure 6 shows this directly.All four MCRaT simulations considered here start off with photons that have an optical depth of 10 3 − 10 4 .As the photons scatter and propagate with    the outflow, their optical depth slowly decreases until it reaches a sharp decay at ∼ 1.8 × 10 13 cm.While the photons in all of our MCRaT simulations reach τ = 1, some only do so at this sharp drop.This rapid decay is due to the sum in Equation 9only going to the last RHD simulation frame, instead of all the way out to infinity.The fact that our MCRaT simulations with Y e = 0.7 and Y e = 1 only reach an optical depth of 1 when this artificial drop occurs is indicative of the fact that the photons in these simulations are still relatively coupled to the outflow.A proxy for this can be seen in Figure 5, which shows the same cooling behavior as panel (c) in Figure 4, with photon energies beginning to level off as they approach the photosphere.In particular, it also shows that the photon energy for the simulation with Y e = 0.1 has nearly leveled off while the energies for the other three simulations are still actively decreasing, indicating that the photons are still scattering with the outflow.
In past works, MCRaT has had successes in reproducing various observational correlations of GRBs (Parsotan et al. 2018).Figure 7 shows the Amati and Yonetoku correlations for the four simulations considered here, with viewing angles of θ v = 1 • , 2 • , 3 • .The Amati relation in panel (a) shows two sets of points, one set corresponding to calculations using photons from the first 20 s of the lightcurves in Figure 1   Injected photons in all four of our simulations start of with similar energy and as the photons propagate further from the central engine photons in simulations with lower values of Ye begin to decouple from the jet sooner, resulting in higher energies for those simulations.The energy from the Ye = 0.1 is nearly constant after R ∼ 10 13 cm, while the rest appear to be somewhat coupled to the jet by the time the photons reach the last RHD simulation frame at R ∼ 10 13 cm.
of the lightcurve we use.With the Amati relation, we find that there is some strain when using photons from all 40 s, which is similar to results from Parsotan & Lazzati (2018).However, we can recover the relation if we restrict ourselves to photons from the first 20 s.This is not an entirely new result, since MCRaT analysis of a similar simulation with a short-lived engine (Parsotan et al. 2018) yielded analogous results.Quali-  9) for all four of our MCRaT simulations.Scatterings for each photon are counted, starting when they're injected near the central engine, and accumulate as they propagate out to the photosphere.Initially, τ ∼ 10 3 − 10 4 which is high enough to ensure that the photons are described by a Planck spectrum.There is a significant drop in τ at ∼ 1.8 × 10 13 cm, which corresponds to the average photon position in the last RHD frame.Photons that are fully decoupled from the outflow have an optical depth of τ ∼ 1, and the MCRaT simulations with Ye = 0.1 and Ye = 0.4 reach this value before the drop.The MCRaT simulations with Ye = 0.7 and Ye = 1, however, reach this value right at the edge of the drop, indicating that these simulations are still somewhat coupled to the outflow.
tatively, it is also expected that shortening the duration of the engine reduces the total burst energy (moving points to the left in the Amati plane) with only a relatively small effect on the peak photon energy, likely in the upward direction since bursts tend to have harder spectra in their early phases.
Figure 8 shows the Golenetskii relation for all values of Y e .Each point is calculated by finding the luminosity and time resolved E pk over 1 s time bins for the first 20 s of the lightcurves in Figure 1.As with the Yonetoku relation, we find good agreement with observations without any restrictions on Y e or photons.Moreover, we find that simulations with a larger neutron component tend to push peak energies and luminosities into better agreement with all three observational correlations.
The role of the neutron component in our simulations can be summarized by plotting spectral parameters as a function of Y e .Panel (a) in Figure 9 shows how the Band parameters α and E 0 depend on Y e , with best-fit power laws shown as dashed and dash-dotted lines.β is not shown due to the lack of a clear pattern in Figure 4. Neither α nor E 0 change very much when Y e is near 1.However, as the size of the neutron component increases, corresponding to our simulations with Y e = 0.4 and Y e = 0.1, the spectral parameters begin to change more dramatically.This is consistent with Figures 5 and  6 showing that simulations with a small neutron component are still somewhat coupled to the outflow.Had the injected photons been able to scatter for longer, it is likely changes would be more consistent across the range of Y e considered here.Furthermore, the nearly symmetric slopes of trend lines in panel (a) are consistent with the strong correlation between α and E 0 on display in Figure 3. Additionally, as suggested by Figures 7 and 8, panel (b) in Figure 9 shows that the radiative efficiency increases as the size of the neutron component is increased, and that this effect isn't dependent on viewing angle for the range considered here.

SUMMARY AND DISCUSSION
In this paper we present results from a series of MCRaT radiative transfer simulations that probe the role that a neutron component in the outflow has on the radiation produced in a LGRB.Varying the density of the input RHD simulation controls the size of the neutron component via the lepton density in Equation 3, which in turn changes how the photons interact with the outflow until they reach the photosphere.
Observables, such as spectra and lightcurves, can be produced with the results of our MCRaT simulations.
Our Y e = 1 lightcurve, and the associated time-resolved spectral parameters, show good agreement with past works using similar 16TI RHD simulations (e.g.Parsotan & Lazzati (2021)).We likewise find good agreement between our Y e = 1 time-integrated spectra and those seen in the same paper.
We find clear patterns in the spectral parameters as we vary Y e .In particular, the break energy E 0 (and thus the corresponding peak energy E pk = [2 + α]E 0 ) is shifted to higher energies as Y e decreases (and the size of the neutron component increases).A power-law fit to E 0 as a function of Y −1 e (E 0 ∝ Y ζ e ) yields an index ζ = −0.26.This behavior is consistent with how the radiation in each of our MCRaT simulations decouple from the outflow.Our simulations with Y e = 1 and Y e = 0.7 are still relatively coupled to electrons in the outflow and so the photons are still appreciably cooling when they reach the last frame of the RHD simulation, resulting in a relatively weak power-law index.We also find that α obtained from simulations with a smaller Y e is consistent with a less thermal spectrum than when Y e is larger, and that this behavior is likely due to a strong correlation between E 0 and α.This is supported corresponding power law for α(Y −1 e ), which is 0.297.In contrast to the other parameters, β has no clear trend, possibly due to the fact that the high en-  2012), with the solid gray line showing the line of best fit.All MCRaT simulations follow the Yonetoku relation, with lower values of Ye corresponding to higher E pk , Eiso, and Liso.Similarly to past work with MCRaT, there is some strain with the Amati relation, but this strain is removed when only considering photons from the first 20 s of the jet, when it is experiencing more shocks.Simulations with more neutrons fit both relations better, regardless of which portion of the lightcurve we consider.We also show how radiation evolves from the injected blackbody to the observed Band-type spectra by conducting mock observations, and calculating spectra, using photons before they have finished scattering through the final RHD simulation frame.This shows that all parameters start off more or less equal across all our simulations, and at some point they begin to diverge until they settle to their final values near the photosphere.In particular, E 0 starts off relatively high and decreases gradually as the injected photons propagate through and with the outflow.The low frequency index α mirrors this behavior, probably due to their strong correlation.
Similar behaviour is observed when we track the optical depth and average energy of the injected photons, beginning immediately after injection until they finish scattering.Both quantities start out high, indicating that that the photons are injected into a hot and dense outflow, and so are well-described by the blackbody spectrum.We see a gradual decoupling of the photons from the outflow, which mirrors the behaviour of the spectral parameters.
Finally, we check our simulations against the observational correlations of Amati, Yontetoku, and Golenetskii (Amati et al. 2002;Yonetoku et al. 2004; Golenetskii et al. 1983), and find good agreement with all three, regardless of Y e .This agrees well with past work with MCRaT (Parsotan et al. 2018).However, given the maximum injection angle of 3 • , we are limited to the number of observations we can make.Interestingly, while all of our simulations fit these correlations nicely, those with a larger neutron component tend to lie closer to the trend lines than those with a smaller neutron component.
Generally, these results are very promising as they provide a mechanism for increasing the peak energy predicted by photospheric models of GRB prompt emission.While there is no consensus on the neutron content of GRB outflows, their presence in both core collapse supernovae and binary neutron star mergers suggests that peak energies are at least somewhat higher than seen in past works with MCRaT.The corresponding increase in total radiated energy (which is inevitable since the number of photons is conserved in a pure scattering process) increases radiative efficiency and brings the MCRaT predictions to better agreement with observational correlations.Both of these results can be interpreted by considering a baryon-loaded LGRB outflow: when the outflow is produced near the central engine, it is hot and dense and thus produces blackbody radiation.The outflow is subsequently heated via shocks as it bores its way through the stellar envelope.Eventually the outflow will clear the envelope and begin to cool while its internal energy is converted to bulk kinetic energy.Thus, the initially hot blackbody radiation also cools as it gradually decouples from the matter component of the outflow.When there is a neutron component in the outflow, radiation will decouple sooner and will thus carry with it a signature of the outflow from when it had converted less of its internal energy into bulk kinetic energy, thereby resulting in the observed increase in radiative efficiency.
An important consideration of the material component of GRB outflows, not treated here, is that of mixing.The jet, cocoon, and stellar envelope could all have different neutron components, and mixing between these could thus modify observables.This effect would likely be more prominent at larger viewing angles where mixing is more prominent.Furthermore, the methods discussed here could naturally be extended to sGRB simulations emerging from binary neutron star mergers.Both of these considerations will be explored in future works.
shows how the Band function parameters α, β, and E 0 vary as a function of the photons' average distance from the central engine for all values of Y e we consider.As with Figure 2, all parameters come from time-integrated spectra.Panels (a) and (c) in Figure 4 clearly show the imprint of a neutron component on the spectral parameters of LGRBs.All four of our MCRaT simulations start off hot near the central engine and gradually cool as the photons and outflow propagate.Simulations with a smaller neutron component cool down more, resulting in lower peak energies.Since E 0 and α are correlated (e.g. Figure 3), the low energy slope α mirrors this behavior, with simulations having larger neutron components displaying smaller values for α.Panel (b), however, shows no clear trend.

Figure 1 .
Figure 1.Light curves and time resolved best fit parameters of the 4 MCRaT simulations: (a.) Ye = 0.1, (b.) Ye = 0.4, (c.) Ye = 0.7, (d.) Ye = 1.The best fit parameter α is shown in red, β is shown in blue, and E pk is shown in green.β is not shown when a comptonized function provides a better fit than the Band function

Figure 2 .
Figure 2. Time-integrated spectra for MCRaT simulations with Ye = 1, shown in red, and Ye = 0.1, shown in blue.In both cases circles show data points and the solid lines show the best fit Band functions.The vertical dashed lines show the break energies, E0, for both simulations.Both spectra were calculated using photons collected over the the first 40 s of the lightcurves in figure 1.

Figure 3 .
Figure 3. Corner plot resulting from fitting the Band function to the spectrum from the Ye = 0.1 simulation with a Markov-Chain Monte Carlo algorithm.The four parameters are the low energy slope α, the high energy slope β, the break energy E0, and the normalization parameter N .This clearly shows a tight correlation between E0 and α, with less prominent correlations between all other parameters.

Figure 4 .
Figure 4. Evolution of the Band function parameters for spectra computed from each value of Ye: (a) the low energy slope α; (b) the high energy slope β; and (c) the break energy E0.Each data point represents a parameter obtained from a mock observation conducted at various intermediate steps throughout the MCRaT simulation.At each step, the injected photons have some average position and a detector was placed slightly beyond that point, denoted R det .As the simulation progresses, the position of the detector moves further away from the central engine and the Band function parameters approach their final values near the photosphere.Panels (a) and (c) show clear patterns for α and E0, respectively.Spectra obtained for all four values of Ye start off hot, having a high E0, and gradually cool as the MCRaT simulations progress.E0 obtained from the Ye = 0.1 simulation levels off sooner than for the other simulations, and so maintains a hotter spectra.This behavior is mirrored in panel (a), with α reaching a smaller value for Ye = 0.1 than for other simulations.Panel (b) shows no discernible pattern for β.
(shown in solid colors), while the other set uses photons from the first 40 s (shown in faded colors).Panel (b) shows the Yonetoku relation for photons obtained only during the first 40 s.Here, we see that our simulations agree well with the Yonetoku relation, regardless of Y e or which portion

Figure 5 .
Figure5.Average photon energy computed as a function of distance from the central engine.Injected photons in all four of our simulations start of with similar energy and as the photons propagate further from the central engine photons in simulations with lower values of Ye begin to decouple from the jet sooner, resulting in higher energies for those simulations.The energy from the Ye = 0.1 is nearly constant after R ∼ 10 13 cm, while the rest appear to be somewhat coupled to the jet by the time the photons reach the last RHD simulation frame at R ∼ 10 13 cm.

Figure 6 .
Figure 6.Optical depth (Equation9) for all four of our MCRaT simulations.Scatterings for each photon are counted, starting when they're injected near the central engine, and accumulate as they propagate out to the photosphere.Initially, τ ∼ 10 3 − 10 4 which is high enough to ensure that the photons are described by a Planck spectrum.There is a significant drop in τ at ∼ 1.8 × 10 13 cm, which corresponds to the average photon position in the last RHD frame.Photons that are fully decoupled from the outflow have an optical depth of τ ∼ 1, and the MCRaT simulations with Ye = 0.1 and Ye = 0.4 reach this value before the drop.The MCRaT simulations with Ye = 0.7 and Ye = 1, however, reach this value right at the edge of the drop, indicating that these simulations are still somewhat coupled to the outflow.

Figure 7
Figure 7. a.) Amati and b.) Yonetoku correlations for all four values of Ye.To obtain multiple observations for each simulation, we conduct a mock observation at three different viewing angles.In each figure, different shapes denote viewing angles and different colors denote different values of Ye.In a.), the solid gray line shows the Amati Relationship from Tsutsui et al. (2009), with the dotted gray line showing the 1σ confidence intervals.The faded colors show data obtained from the first 40 s of the lightcurves in Figure 1, while the solid colors show only the first 20 s.In b), the gray dots show observational data of GRBs from Nava et al. (2012), with the solid gray line showing the line of best fit.All MCRaT simulations follow the Yonetoku relation, with lower values of Ye corresponding to higher E pk , Eiso, and Liso.Similarly to past work with MCRaT, there is some strain with the Amati relation, but this strain is removed when only considering photons from the first 20 s of the jet, when it is experiencing more shocks.Simulations with more neutrons fit both relations better, regardless of which portion of the lightcurve we consider.

Figure 8 .
Figure 8. Golenetskii relation for all values of Ye over the first 40 s of each burst.Each value of Ye is denoted by a different color, and each point is calculated by binning the lightcurves shown in Figure 1 into 1 s bins and calculating the time resolved E pk for each bin.The gray solid indicates the Golenetskii relation from Lu et al. (2012), with the dotted gray lines representing the 2σ intervals.Every simulation shows good agreement with the Golenetskii relation, with smaller values of Ye corresponding to higher values of E pk and Luminosity, similar to Figure 7

Figure 9 .
Figure 9.The Ye effect on (a) the Band function parameters α and E0, and on (b) total radiated energy.The x-axis in both panels shows Y −1 e so the size of the neutron component increases to the right.In (a), the red squares show the break energy E0 and the black triangles show the low energy slope α, with the dashed and dashed-dotted lines showing the best fit trend lines for E0 and α, respectively.The break energy E0 clearly increases as the neutron component gets larger, and α clearly decreases nearly symmetrically as evidenced by the E0 slope of -0.26 and the α slope of 0.297.The low energy slope β is not shown due to a lack of a clear pattern in Figure 4.In (b) the different colors show the isotropic energy from mock observations conducted at different viewing angles.As the neutron component is increased, the total radiated energy is increased across multiple viewing angles.

NFacilities:
Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.REFERENCES