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

Nathan Walker Oregon State University
Department of Physics, 301 Weniger Hall, Oregon State University
Corvallis, OR 97331, USA  Tyler Parsotan Astrophysics Science Division, NASA
Goddard Space Flight Center,Greenbelt, MD 20771, USA  Davide Lazzati Oregon
State University
Department of Physics, 301 Weniger Hall, Oregon State University
Corvallis, OR 97331, USA

###### Abstract

Long-duration gamma-ray bursts (LGRBs), thought to be produced during core-
collapse supernovæ, 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 post-processing using the Monte Carlo Radiation Transfer
(MCRaT) 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 $\alpha$
and $E_{0}$, while the picture with the $\beta$ 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.

Gamma-ray bursts(629)) — Radiative transfer simulations(1967) — Hydrodynamical
simulations (767)

††facilities: 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.

##  1 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 supernovæ, a process in which stars more massive than about
$8M_{\odot}$ 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 relativistic 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
supernovæ 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.

##  2 Computational Methods

###  2.1 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 $\delta\theta$, defined with respect to the jet axis. All
cells within the interval $\delta\theta$ and with a radius between
$R_{inj}\pm\frac{1}{2}\,c\delta t$ are selected, where $c$ is the speed of
light and $\delta 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,

$T^{{}^{\prime}}_{i}=\left(\frac{3p_{i}}{a}\right)^{\frac{1}{4}},$ (1)

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)

$dN_{i}=\frac{\xi\,T_{i}^{{}^{\prime}3}\,\Gamma_{i}}{w}\,dV_{i},$ (2)

where $dN_{i}$ is the expected number of photons in the $i^{th}$ RHD cell,
$\xi$ is the photon number density coefficient from $n_{\gamma}=\xi T^{3}$
($\xi=20.29$ for a Planck spectrum and $\xi=8.44$ for a Wein spectrum),
$T_{i}^{{}^{\prime}}$ is the comoving fluid temperature, $\Gamma_{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)

$\lambda_{i}=\frac{dr}{d\tau_{T}}=\frac{1}{\sigma_{T}\,n_{i}^{{}^{\prime}}\,\Gamma_{i}\left(1-\beta_{i}\cos\theta_{fl,i}\right)},$
(3)

where $\sigma_{T}$ is the Thomson cross section, $n_{i}^{{}^{\prime}}$ is the
comoving lepton number density, $\beta_{i}$ is the fluid velocity in units of
$c$, and $\theta_{fl,i}$ is the angle between the fluid velocity and the
photon velocity. A random scattering time for each photon is drawn from the
distribution

$P_{i}(t)\propto e^{-\frac{c}{\lambda_{i}}\,t},$ (4)

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.

###  2.2 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
four-momentum at the end of the simulation, and the detection time is
calculated as

$t_{d}=t_{p}+t_{real}-t_{j},$ (5)

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\times 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,

$L_{t}=\frac{1}{\Delta\Omega\Delta t_{bin}}\sum_{i}w_{i}E_{i},$ (6)

where $E_{i}$ is the energy of each photon, $\Delta t_{bin}$ is the time bin,
and
$\Delta\Omega=2\pi[\cos(\theta_{v}-\Delta\theta/2)-\cos(\theta_{v}+\Delta\theta/2)]$
is the solid angle the detector occupies, with $\Delta\theta$ being the
angular acceptance range centered around $\theta_{v}$. We also bin the photon
energies to calculate spectra via

$\frac{dN_{e}(E)}{dE\,dt}=\frac{1}{\Delta E_{bin}\Delta\Omega\Delta
t}\sum_{i}w_{i},$ (7)

where all the terms are the same as in Equation 6 except $\Delta E_{bin}$ and
$\Delta 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),

$\displaystyle\begin{split}&N_{E}(E)=\\\
&A\left(\frac{E}{100\text{keV}}\right)^{\alpha}\exp(-E/E_{0}),\\\
&\qquad\qquad\qquad\qquad\qquad E\leq(\alpha-\beta)E_{0}\\\
&A\left[\frac{(\alpha-\beta)E_{0}}{100\text{keV}}\right]^{\alpha-\beta}\exp(\beta-\alpha)\left(\frac{E}{100\text{keV}}\right)^{\beta}\\\
&\qquad\qquad\qquad\qquad\qquad E\geq(\alpha-\beta)E_{0}.\end{split}$ (8)

to spectra obtained from equation 7. In equation 8, $\alpha$ and $\beta$ 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 $E_{pk}=(2+\alpha)\,E_{0}$.

In order for the calculated spectra and light curves to correspond to what an
observer would see, the optical depth must reach a value $\tau\sim 1$. We
calculate the optical depth (Parsotan et al., 2018) as:

$\tau_{i}^{n}=\sum_{j=i}^{L}\left<N\right>_{j}^{n},$ (9)

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 $\left<N\right>_{j}^{n}$ 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 $\lesssim 1$,
corresponding to the photosphere condition of $\tau\sim 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.

###  2.3 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_{-}-n_{+}=Y_{e}\,\frac{\rho}{m_{p}},$ (10)

where $n_{\pm}$ are the $e^{\pm}$ number densities. In the absence of
$e^{\pm}$ 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 $\rho$ 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,

$e^{-}+p\rightarrow n+\nu,\,\,\,\,\,e^{+}+n\rightarrow p+\bar{\nu},$ (11)

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}$: $\rho\rightarrow Y_{e}\,\rho$.

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^{\circ}$ 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.

##  3 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
$\sim 10^{13}$ cm. The 16TI simulation in Lazzati et al. (2013) was performed
on an adaptive mesh grid with a maximum resolution of $4\times 10^{6}$ cm and
output files were saved every $\delta 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\delta t/\delta r$ which, with
the spatial and temporal resolutions from the 16TI simulation used here,
results in a light crossing ratio as large as $\sim 1500$.

Our methods are similar to Parsotan et al. (2018), with a key difference being
that we inject $\sim 2\times 10^{5}$ photons for $\sim 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
$\theta_{v}<3$°when conducting mock observations. For the electron-to-nucleon
ratio we use the values $Y_{e}=1,0.7,0.4,\,\text{and}\,\,0.1$ to cover the
cases of a small to large neutron component.

Figure 1 shows lightcurves obtained at a viewing angle of
$\theta_{v}=1^{\circ}$ alongside the time-resolved best fit parameters
$\alpha$ and $\beta$ for the Band function (equation 8), in addition to the
peak energy $E_{pk}=(2+\alpha)\,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 $\sim 8$ s, with a
brighter peak at $\sim 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 ProcessMCRaT, 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 $\alpha$, while the other pairs of parameters
have no notable correlations. This strong correlation between $\alpha$ and
$E_{0}$ plays a part in the evolution of Band function parameters for all four
of the MCRaT simulations in this work.

Figure 1: Light curves and time resolved best fit parameters of the 4 MCRaT
simulations: (a.) $Y_{e}=0.1$, (b.) $Y_{e}=0.4$, (c.) $Y_{e}=0.7$, (d.)
$Y_{e}=1$. The best fit parameter $\alpha$ is shown in red, $\beta$ is shown
in blue, and $E_{pk}$ is shown in green. $\beta$ is not shown when a
comptonized function provides a better fit than the Band function  Figure 2:
Time-integrated spectra for MCRaT simulations with $Y_{e}=1$, shown in red,
and $Y_{e}=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, $E_{0}$, for both simulations. Both spectra were
calculated using photons collected over the the first 40 s of the lightcurves
in figure 1. Figure 3: Corner plot resulting from fitting the Band function to
the spectrum from the $Y_{e}=0.1$ simulation with a Markov-Chain Monte Carlo
algorithm. The four parameters are the low energy slope $\alpha$, the high
energy slope $\beta$, the break energy $E_{0}$, and the normalization
parameter $N$. This clearly shows a tight correlation between $E_{0}$ and
$\alpha$, with less prominent correlations between all other parameters.

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 shows how the Band function
parameters $\alpha$, $\beta$, 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 $\alpha$ are correlated
(e.g. Figure 3), the low energy slope $\alpha$ mirrors this behavior, with
simulations having larger neutron components displaying smaller values for
$\alpha$. Panel (b), however, shows no clear trend.

Figure 4: Evolution of the Band function parameters for spectra computed from
each value of $Y_{e}$: (a) the low energy slope $\alpha$; (b) the high energy
slope $\beta$; and (c) the break energy $E_{0}$. 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 $\alpha$ and $E_{0}$, respectively. Spectra obtained for
all four values of $Y_{e}$ start off hot, having a high $E_{0}$, and gradually
cool as the MCRaT simulations progress. $E_{0}$ obtained from the $Y_{e}=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 $\alpha$
reaching a smaller value for $Y_{e}=0.1$ than for other simulations. Panel (b)
shows no discernible pattern for $\beta$.

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 $\sim 1.8\times 10^{13}$ cm. While the photons in all of our MCRaT
simulations reach $\tau=1$, some only do so at this sharp drop. This rapid
decay is due to the sum in Equation 9 only 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.

Figure 5: 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 $Y_{e}$ begin to decouple from the
jet sooner, resulting in higher energies for those simulations. The energy
from the $Y_{e}=0.1$ is nearly constant after $R\sim 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\sim 10^{13}$ cm.  Figure 6: Optical depth
(Equation 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, $\tau\sim
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 $\tau$ at $\sim 1.8\times
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 $\tau\sim 1$, and the MCRaT simulations with $Y_{e}=0.1$ and $Y_{e}=0.4$
reach this value before the drop. The MCRaT simulations with $Y_{e}=0.7$ and
$Y_{e}=1$, however, reach this value right at the edge of the drop, indicating
that these simulations are still somewhat coupled to 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 $\theta_{v}=1^{\circ},\,2^{\circ},\,3^{\circ}$. 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 (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 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. Qualitatively, 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 $\alpha$ and $E_{0}$ depend on $Y_{e}$, with
best-fit power laws shown as dashed and dash-dotted lines. $\beta$ is not
shown due to the lack of a clear pattern in Figure 4. Neither $\alpha$ 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 $\alpha$ 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.

Figure 7: a.) Amati and b.) Yonetoku correlations for all four values of
$Y_{e}$. 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
$Y_{e}$. In a.), the solid gray line shows the Amati Relationship from Tsutsui
et al. (2009), with the dotted gray line showing the $1\sigma$ 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 $Y_{e}$ corresponding to
higher $E_{pk}$, $E_{iso}$, and $L_{iso}$. 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:
Golenetskii relation for all values of $Y_{e}$ over the first 40 s of each
burst. Each value of $Y_{e}$ 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\sigma$ intervals. Every simulation shows good agreement
with the Golenetskii relation, with smaller values of $Y_{e}$ corresponding to
higher values of $E_{pk}$ and Luminosity, similar to Figure 7 Figure 9: The
$Y_{e}$ effect on (a) the Band function parameters $\alpha$ and $E_{0}$, and
on (b) total radiated energy. The x-axis in both panels shows $Y_{e}^{-1}$ so
the size of the neutron component increases to the right. In (a), the red
squares show the break energy $E_{0}$ and the black triangles show the low
energy slope $\alpha$, with the dashed and dashed-dotted lines showing the
best fit trend lines for $E_{0}$ and $\alpha$, respectively. The break energy
$E_{0}$ clearly increases as the neutron component gets larger, and $\alpha$
clearly decreases nearly symmetrically as evidenced by the $E_{0}$ slope of
-0.26 and the $\alpha$ slope of 0.297. The low energy slope $\beta$ 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.

##  4 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+\alpha]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_{e}^{-1}$ ($E_{0}\propto Y_{e}^{\zeta}$) yields an index
$\zeta=-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 $\alpha$ 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 $\alpha$. This is supported corresponding power law for
$\alpha(Y_{e}^{-1})$, which is 0.297. In contrast to the other parameters,
$\beta$ has no clear trend, possibly due to the fact that the high energy tail
of the spectrum forms relatively close to the photosphere compared to the
lower frequency parts of the spectrum, which are characterized by $\alpha$ and
$E_{0}$.

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 $\alpha$
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^{\circ}$, 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 supernovæ 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.

N. W. and D.L. acknowledge support from NSF grant AST-1907955

## References

  * Abramowicz et al. (1991) Abramowicz, M. A., Novikov, I. D., & Paczynski, B. 1991, The Astrophysical Journal, 369, 175, doi: 10.1086/169748
  * Amati et al. (2002) Amati, L., Frontera, F., Tavani, M., et al. 2002, Astronomy and Astrophysics, 390, 81, doi: 10.1051/0004-6361:20020722
  * Arita-Escalante et al. (2023) Arita-Escalante, J., Parsotan, T., & Cenko, S. B. 2023, The Astrophysical Journal, 954, 208, doi: 10.3847/1538-4357/ace7d3
  * Band et al. (1993) Band, D., Matteson, J., Ford, L., et al. 1993, The Astrophysical Journal, 413, 281, doi: 10.1086/172995
  * Beloborodov (2003) Beloborodov, A. M. 2003, The Astrophysical Journal, 588, 931, doi: 10.1086/374217
  * Beloborodov (2010a) —. 2010a, Monthly Notices of the Royal Astronomical Society, 407, 1033, doi: 10.1111/j.1365-2966.2010.16770.x
  * Beloborodov (2010b) —. 2010b, Monthly Notices of the Royal Astronomical Society, 407, 1033, doi: 10.1111/j.1365-2966.2010.16770.x
  * Beloborodov (2013) —. 2013, The Astrophysical Journal, 764, 157, doi: 10.1088/0004-637X/764/2/157
  * Bloom et al. (1999) Bloom, J. S., Kulkarni, S. R., Djorgovski, S. G., et al. 1999, Nature, 401, 453, doi: 10.1038/46744
  * Bucciantini et al. (2012) Bucciantini, N., Metzger, B. D., Thompson, T. A., & Quataert, E. 2012, MNRAS, 419, 1537, doi: 10.1111/j.1365-2966.2011.19810.x
  * Chhotray & Lazzati (2015) Chhotray, A., & Lazzati, D. 2015, The Astrophysical Journal, 802, 132, doi: 10.1088/0004-637X/802/2/132
  * Costa et al. (1997) Costa, E., Frontera, F., Heise, J., et al. 1997, Nature, 387, 783, doi: 10.1038/42885
  * Derishev et al. (1999) Derishev, E. V., Kocharovsky, V. V., & Kocharovsky, V. V. 1999, The Astrophysical Journal, 521, 640, doi: 10.1086/307574
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306, doi: 10.1086/670067
  * Galama et al. (1998) Galama, T. J., Vreeswijk, P. M., van Paradijs, J., et al. 1998, Nature, 395, 670, doi: 10.1038/27150
  * Ghisellini et al. (2000) Ghisellini, G., Celotti, A., & Lazzati, D. 2000, Monthly Notices of the Royal Astronomical Society, 313, L1, doi: 10.1046/j.1365-8711.2000.03354.x
  * Giannios & Spruit (2007) Giannios, D., & Spruit, H. C. 2007, Astronomy & Astrophysics, 469, 1, doi: 10.1051/0004-6361:20066739
  * Golenetskii et al. (1983) Golenetskii, S. V., Mazets, E. P., Aptekar, R. L., & Ilyinskii, V. N. 1983, Nature, 306, 451, doi: 10.1038/306451a0
  * Groot et al. (1998) Groot, P. J., Galama, T., van Paradijs, J., & Kouveliotou, C. 1998, The Astrophysical Journal, 493, L27, doi: 10.1086/311125
  * Hjorth et al. (2003) Hjorth, J., Sollerman, J., Møller, P., et al. 2003, Nature, 423, 847, doi: 10.1038/nature01750
  * Ito et al. (2018) Ito, H., Levinson, A., Stern, B. E., & Nagataki, S. 2018, Monthly Notices of the Royal Astronomical Society, 474, 2828, doi: 10.1093/mnras/stx2722
  * Klebesadel et al. (1973) Klebesadel, R. W., Strong, I. B., & Olson, R. A. 1973, The Astrophysical Journal, 182, L85, doi: 10.1086/181225
  * Lazzati (2016) Lazzati, D. 2016, The Astrophysical Journal, 829, 76, doi: 10.3847/0004-637X/829/2/76
  * Lazzati et al. (2009) Lazzati, D., Morsony, B. J., & Begelman, M. C. 2009, The Astrophysical Journal, 700, L47, doi: 10.1088/0004-637X/700/1/L47
  * Lazzati et al. (2013) Lazzati, D., Morsony, B. J., Margutti, R., & Begelman, M. C. 2013, The Astrophysical Journal, 765, 103, doi: 10.1088/0004-637X/765/2/103
  * Lu et al. (2012) Lu, R.-J., Wei, J.-J., Liang, E.-W., et al. 2012, The Astrophysical Journal, 756, 112, doi: 10.1088/0004-637X/756/2/112
  * MacFadyen et al. (2001) MacFadyen, A. I., Woosley, S. E., & Heger, A. 2001, The Astrophysical Journal, 550, 410, doi: 10.1086/319698
  * Morsony et al. (2007) Morsony, B. J., Lazzati, D., & Begelman, M. C. 2007, ApJ, 665, 569, doi: 10.1086/519483
  * Narayan et al. (2001) Narayan, R., Piran, T., & Kumar, P. 2001, The Astrophysical Journal, 557, 949, doi: 10.1086/322267
  * Nava et al. (2012) Nava, L., Salvaterra, R., Ghirlanda, G., et al. 2012, Monthly Notices of the Royal Astronomical Society, 421, 1256, doi: 10.1111/j.1365-2966.2011.20394.x
  * Parsotan & Lazzati (2018) Parsotan, T., & Lazzati, D. 2018, The Astrophysical Journal, 853, 8, doi: 10.3847/1538-4357/aaa087
  * Parsotan & Lazzati (2021) —. 2021, The Astrophysical Journal, 922, 257, doi: 10.3847/1538-4357/ac2428
  * Parsotan & Lazzati (2022) —. 2022, The Astrophysical Journal, 926, 104, doi: 10.3847/1538-4357/ac4093
  * Parsotan et al. (2018) Parsotan, T., Lopez-Camara, D., & Lazzati, D. 2018, The Astrophysical Journal, 869, 103, doi: 10.3847/1538-4357/aaeed1
  * Pe’er & Ryde (2011) Pe’er, A., & Ryde, F. 2011, The Astrophysical Journal, 732, 49, doi: 10.1088/0004-637X/732/1/49
  * Pe’er et al. (2006) Pe’er, A., Mészáros, P., & Rees, M. J. 2006, The Astrophysical Journal, 642, 995, doi: 10.1086/501424
  * Preece et al. (1998) Preece, R. D., Briggs, M. S., Mallozzi, R. S., et al. 1998, The Astrophysical Journal, 506, L23, doi: 10.1086/311644
  * Rees & Meszaros (1994) Rees, M. J., & Meszaros, P. 1994, The Astrophysical Journal, 430, L93, doi: 10.1086/187446
  * Rossi et al. (2004) Rossi, E. M., Beloborodov, A. M., & Rees, M. J. 2004, in AIP Conference Proceedings, Vol. 727, 198–202, doi: 10.1063/1.1810830
  * Ryde et al. (2011) Ryde, F., Pe’er, A., Nymark, T., et al. 2011, Monthly Notices of the Royal Astronomical Society, 415, 3693, doi: 10.1111/j.1365-2966.2011.18985.x
  * Tsutsui et al. (2009) Tsutsui, R., Nakamura, T., Yonetoku, D., et al. 2009, Journal of Cosmology and Astroparticle Physics, 2009, 015, doi: 10.1088/1475-7516/2009/08/015
  * Vurm & Beloborodov (2016) Vurm, I., & Beloborodov, A. M. 2016, The Astrophysical Journal, 831, 175, doi: 10.3847/0004-637X/831/2/175
  * Woosley & Bloom (2006) Woosley, S., & Bloom, J. 2006, Annual Review of Astronomy and Astrophysics, 44, 507, doi: 10.1146/annurev.astro.43.072103.150558
  * Woosley & Janka (2005) Woosley, S., & Janka, H.-T. 2005, Nature Physics, 1, 147, doi: 10.1038/nphys172
  * Woosley & Heger (2006) Woosley, S. E., & Heger, A. 2006, The Astrophysical Journal, 637, 914, doi: 10.1086/498500
  * Yonetoku et al. (2004) Yonetoku, D., Murakami, T., Nakamura, T., et al. 2004, The Astrophysical Journal, 609, 935, doi: 10.1086/421285