# High energy neutrino production in gamma-ray bursts: dependence of the
neutrino signal on the jet composition

Valentin De Lia ID and Irene Tamborra ID

###### Abstract

Heavy nuclei can be synthetized or entrained in gamma-ray bursts (GRBs) with
implications on the high-energy neutrino emission. By means of a Monte-Carlo
algorithm, we model nuclear cascades and investigate their impact on the
neutrino production considering kinetic dominated jets (in the internal shock
model, including a dissipative photosphere) as well as Poynting flux dominated
jets (for a jet model invoking internal-collision-induced magnetic
reconnection and turbulence, ICMART). We find that the ICMART model allows for
efficient nuclear cascades leading to an overall larger neutrino fluence than
in the other two jet models. The survival of nuclei and inefficient nuclear
cascades lead to an overall reduction of the neutrino fluence up to one order
of magnitude. However, if nuclei are disintegrated, the neutrino fluence may
be comparable to the one emitted from a jet loaded with protons. Exploring the
parameter space of jet properties, we conclude that the composition and the
bulk Lorentz factor have significant impact on the efficiency of nuclear
cascades as well as the spectral shape of the expected neutrino fluence. On
the other hand, the neutrino spectral distribution is less sensitive to the
power-law index of the accelerated population of protons or heavier nuclei.
For what concerns the diffuse emission of neutrinos from GRBs, we find that
the uncertainty due to the jet composition can be at most comparable to the
one related to the GRB cosmological rate.

##  1 Introduction

Gamma-Ray Bursts (GRBs) are among the most energetic transient astrophysical
phenomena discovered up to date, releasing isotropic energies up to
$E_{\text{iso}}=1.2\times 10^{55}$ erg within a few seconds [1]. The
observation of GRB 170817A in connection with the gravitational wave event GW
170817 has provided the first confirmation that binary neutron star mergers
can harbor ultra-relativistic jets [2, 3, 4, 5, 6, 7, 8], known as short GRBs.
Supernovae have instead been observed in connection with long duration gamma-
ray bursts [9, 10]. However, recent observations challenge the conjectured
different origin of long and short GRBs [11, 12, 13]. Observed long duration
GRBs can be divided in high-luminosity (HL, with isotropic luminosity
$\mathcal{O}(10^{51}$–$10^{53})$ erg/s) and low-luminosity (LL, with isotropic
luminosity smaller than $\mathcal{O}(10^{49})$ erg/s) GRBs [14, 15]. It is
believed that HL and LL GRBs share similar properties, with LL GRBs having
smaller Lorentz boost factor or perhaps corresponding to HL GRBs observed off-
axis. The dynamical evolution of the GRB jet is linked to the properties of
the central engine powering the jet; the jet could be magnetically driven, if
the central engine harbors a strong magnetic field, or it would be better
described by a fireball when the magnetic activity subsides [16, 17].

Gamma-ray bursts are observed as irregular pulses of gamma-rays, with a non-
thermal photon spectral energy distribution. Such non-thermal features suggest
efficient particle acceleration in the jet [18, 19]. Nevertheless, the physics
driving the observed signals is still subject of active debate [20, 21, 17,
22, 23].

Gamma-ray bursts are expected to be sources of ultra-high energy cosmic rays
and high energy neutrinos [24, 25, 26, 27, 28, 29], produced also thanks to
heavier nuclei synthetized or entrained in the jet [30, 31] and eventually
interacting with the photon background. In the aftermath of particle
acceleration, the inelastic interaction of non-thermal high-energy photons
with ultra-relativistic protons or nuclei could lead to the production of
high-energy neutrinos, see e.g. Refs. [32, 33, 24, 34, 35, 36].

It is expected that the GRB jet may entrain nuclei, e.g. loaded at the base of
the jet, captured from the stellar material surrounding the jet and entrained
during propagation, or synthetized in situ [30, 31, 37]. Such nuclei could be
accelerated in the GRB jet and therefore survive photodisintegration against
the intense photon field, only if the dissipation radius and the bulk Lorentz
factor are large enough. Assuming that such conditions are fulfilled, the
resulting neutrino production would be strongly affected.

Previous work [32, 38] investigated the impact of the jet composition on the
high energy neutrino flux, focusing on the internal shock model [18]. In the
light of the strong dependence of the neutrino signal on the GRB emission
mechanism [33], we extend such exploration to the photospheric [39] and the
internal-collision-induced magnetic reconnection and turbulence (ICMART) [40]
models, and explore the neutrino production for a range of jet properties.

This work is organized as follows. Section 2 provides an overview of the
internal shock model, the photospheric one, and the ICMART model. Section 3
offers a summary of the acceleration and cooling timescales adopted to compute
the maximum energies of protons and nuclei. A description of how neutrinos are
produced through photohadronic interactions of protons and neutrons or nuclei
is presented in Sec. 4. Section 5 focuses on the formalism adopted to compute
the neutrino emission from jets loaded with heavier nuclei, including the
modeling of the nuclear cascades using a Monte-Carlo method. Our findings on
the impact of the jet composition on the neutrino emission for our three jet
models are presented in Sec. 6, together with an exploration of the jet
parameter space. We also discuss the consequences that the jet composition
should have on the diffuse emission of high-energy neutrinos in Sec. 7.
Finally, we critically review our findings and conclude in Sec. 8. Details on
the Monte Carlo simulation of nuclear cascades are provided in Appendix A,
while the nuclear photohadronic model is outlined in Appendix B.

##  2 Gamma-ray burst models

This section highlights the main features of the kinetic dominated and
Poynting flux dominated jets considered in this work as well as the adopted
photon spectral energy distributions. We refer the interested reader to Ref.
[33] for a broader overview.

In the following, we distinguish among three reference frames: the comoving
jet frame, the central engine frame, and the observer frame; any jet quantity
$X$ is denoted as $X^{\prime}$, $\tilde{X}$, or $X$ respectively according to
the adopted reference frame. For example, energy or length transform as
$X={\tilde{X}}/{(1+z)}={\Gamma\;X^{\prime}}/{(1+z)}$, while time transforms as
$X=(1+z)\;\tilde{X}={(1+z)\;X^{\prime}}/{\Gamma}$. We consider a relativistic
jet propagating with bulk Lorentz factor $\Gamma$ with respect to the central
engine. During the prompt phase, it is expected that $\theta_{\rm
op}\gg\Gamma^{-1}$, with $\theta_{\rm op}$ being the half opening angle [41].
This justifies the employment of the isotropic equivalent energy
($\tilde{E}_{\text{iso}}$) representing the energy content of the outflow.
Additionally, GRBs are characterized by two timescales. The first one is the
total duration of the event $t_{\text{dur}}$ and the second one is the average
variability timescale of the central engine $t_{\text{var}}$. We thus define
the total isotropic luminosity
$\tilde{L}_{\text{iso}}=\tilde{E}_{\text{iso}}/\tilde{t}_{\text{dur}}$. Any
burst is also characterized by its redshift $z$.

###  2.1 Internal shock model

The internal shock (IS) model has been widely investigated for the GRB prompt
emission [18, 42, 43] because it allows to easily explain the variability of
the GRB lightcurves. Due to the erratic activity of the the central engine,
the outflow results in several shells, each with different Lorentz factor
$\Gamma$. When a faster shell hits a slower shell, dissipation of kinetic
energy occurs. In the following, we focus on a one-zone collision model for
the sake of simplicity, assuming that all collisions are identical. However,
multi-zone collisions models are expected to lead to less efficient neutrino
emission [44, 45, 46, 47, 35]. We assume that, for a GRB of duration
$\tilde{t}_{\text{dur}}$ and variability time scale $\tilde{t}_{\text{var}}$,
$N_{\text{shock}}=\tilde{t}_{\text{dur}}/\tilde{t}_{\text{var}}$ identical
shocks occur at the radius
$R_{\gamma}=R_{\text{IS}}=2\,\Gamma^{2}c\,\tilde{t}_{\text{var}}$, within a
volume $V^{\prime}_{s}=4\pi\,R_{\text{IS}}^{2}\,c\,t^{\prime}_{\text{var}}$.
Unless otherwise stated, the benchmark jet properties adopted for this model
are summarized in Table 1; note that the large jet luminosity considered for
our benchmark jet model most likely does not allow for the survival of all
iron nuclei. While the assessment of the jet conditions leading to a specific
jet composition is beyond the scope of this work, we generalize our findings
in the second part of this paper, exploring variations of some of the jet
characteristic properties.

The fraction of energy dissipated in particle acceleration, and the fractions
of energy that effectively accelerates seed particles (protons, nuclei),
electrons, or contributes to amplify the magnetic field are
$\varepsilon_{\text{d}}=0.2$, $\varepsilon_{A}=0.1$ (where $A$ indicates
protons or nuclei), $\varepsilon_{e}=0.01$, and $\varepsilon_{B}=0.1$ [48,
49]. Hence the magnetic field is

$B^{\prime}=\sqrt{8\pi\frac{\varepsilon_{B}\,\varepsilon_{d}\;E^{\prime}_{\text{iso,s}}}{V^{\prime}_{s}}}\
,$ (2.1)

with $E^{\prime}_{\text{iso,s}}=E^{\prime}_{\text{iso}}/N_{\text{shock}}$.

Shocks accelerate a certain fraction of the charged particles into a power-law
distribution [50]:

$n^{\prime}_{A}(E^{\prime}_{A})=C{E^{\prime}_{A}}^{-k}\exp\left[-\left(\frac{E^{\prime}_{A}}{E^{\prime}_{A,\max}}\right)^{2}\right]\Theta(E^{\prime}_{A}-E^{\prime}_{A,\min})\
,$ (2.2)

with $k\simeq 2.2$; imposing that
$U_{A}=\int_{E^{\prime}_{A,\min}}^{E^{\prime}_{A,\max}}n^{\prime}_{A}(E^{\prime}_{A})\;E^{\prime}_{A}\;dE^{\prime}_{A}$,
we compute the normalization constant $C$, where
$U_{A}={\varepsilon_{d}\;\varepsilon_{A}\;E^{\prime}_{\text{iso},s}}/{V^{\prime}_{s}}$
is the energy density of accelerated seed particles.

Table 1: Summary of the characteristic parameters assumed for our benchmark jet for the internal shock (IS), photospheric (PHOTO) and ICMART models, respectively. Parameter | IS | PHOTO | ICMART
---|---|---|---
$\tilde{E}_{\text{iso}}$ |  $4.5\times 10^{54}$ erg [51]
$z$ |  $2$ [52]
$\Gamma$ |  $300$ [53]
$t_{\text{dur}}$ |  $30$ s [54]
$t_{\text{var}}$ |  $0.5$ s [55]
$\varepsilon_{\text{d}}$ |  $0.2$ [42, 44] |  $0.2$ [42, 44] |  $0.35$ [56]
$\varepsilon_{e}$ |  $0.01$ [48, 49] |  $0.01$ [48, 49] |  $0.5$ [57]
$\varepsilon_{X}$ |  $0.1$ [49] |  $0.1$ [49] |  $0.5$ [57]
$\varepsilon_{B}$ |  $0.1$ [49] |  $0.1$ [49] | n/a
$k$ |  $2.2$ [58, 59] |  $2.2$ [58, 59] |  $2.0$ [60, 61]
$\sigma$ | n/a | n/a |  $45$ [33]

The photon spectrum is modeled through a Band function [62]:

$n^{\prime}_{\gamma}(E^{\prime}_{\gamma})=C\begin{cases}\displaystyle\left(\frac{E^{\prime}_{\gamma}}{E^{\prime}_{0}}\right)^{\alpha}\exp\left(-\frac{(\alpha+2)E^{\prime}_{\gamma}}{E^{\prime}_{\gamma,p}}\right)&\text{if
}E^{\prime}_{\gamma,\text{min}}\leq E^{\prime}_{\gamma}\leq
E^{\prime}_{\gamma,c}\\\
\displaystyle\left(\frac{E^{\prime}_{\gamma}}{E^{\prime}_{0}}\right)^{\beta}\exp(\beta-\alpha)\left(\frac{E^{\prime}_{\gamma,c}}{E^{\prime}_{0}}\right)^{\alpha-\beta}&\text{if
}E^{\prime}_{\gamma,c}\leq E^{\prime}_{\gamma}\leq
E^{\prime}_{\gamma,\text{max}}\\\ \end{cases}$ (2.3)

with $E^{\prime}_{\gamma,p}$, $E^{\prime}_{\gamma,c}$,
$E^{\prime}_{\gamma,\text{min}}$, $E^{\prime}_{\gamma,\text{max}}$ being the
peak, cooling, minimum and maximum energies of the spectrum,
$E^{\prime}_{\gamma,0}=[{(1+z)}/{\Gamma}]\times 100$ keV, respectively.
Relying on GRB observations from Fermi, we choose $\alpha\simeq-1.1$ and
$\beta\simeq-2.2$ [63]. The normalization constant is given by
$U_{\gamma}=\int_{E^{\prime}_{\gamma,\min}}^{E^{\prime}_{\gamma,\max}}n^{\prime}_{\gamma}(E^{\prime}_{\gamma})\;E^{\prime}_{\gamma}\;dE^{\prime}_{\gamma}$
with
$U_{\gamma}={\varepsilon_{d}\;\varepsilon_{e}\;E^{\prime}_{\text{iso},s}}/{V^{\prime}_{s}}$.
The peak energy is computed following the Amati relation [64]:
$\tilde{E}_{\gamma,p}=80\;\left({\tilde{E}_{\gamma,\text{iso}}}/{10^{52}\;\text{erg}}\right)^{0.57}\text{keV}$
with
$\tilde{E}_{\gamma,\text{iso}}=\varepsilon_{d}\varepsilon_{e}\tilde{E}_{\text{iso}}$,
and the cooling energy is
$E^{\prime}_{\gamma,c}=({\alpha-\beta})/({\alpha+2})E^{\prime}_{\gamma,p}$.
Note that our benchmark $\tilde{E}_{\rm{iso}}$ is such that, in the internal
shock model,
$\tilde{L}_{\rm{iso},\gamma}\simeq\varepsilon_{d}\varepsilon_{e}\tilde{E}_{\rm{iso}}/t_{\rm{dur}}$
is slightly lower than what often considered in the literature. Yet, this
value of $\tilde{E}_{\rm{iso}}$ [51] allows to recover the observed
$\tilde{L}_{\rm{iso},\gamma}$ for the photosperic and ICMART models; we choose
to rely on the same jet benchmark properties to allow for a fair comparison
across jet models.

Note that we do not consider neutrino production below the photosphere [65,
36, 66, 67, 68, 69, 23], as neutrinos would be produced with GeV–TeV energies,
while we are interested in comparing different models for neutrino production
in the optically thin region with overall larger energies.

###  2.2 Photospheric model

The photospheric model (PHOTO) is a variation of the internal shock model that
assumes that dissipative processes shape the radiation bulk produced in the
optically thick region below the photosphere [20]. We assume all
characteristic jet parameters as being the same as the ones of the internal
shock model, except for the ones related to the photon spectrum. Protons or
nuclei interact with a photon spectrum described by three components [39, 33]:
a non-thermal photospheric component defined by a Band spectrum and undergoing
dissipation at the internal shock radius, a photospheric Compton up-scattered
spectrum, and a synchrotron component due to electrons cooling in the shock
region. The accelerated particle spectrum is given by Eq. (2.2) with $k=2.2$;
the photon spectrum is [33]:

$n^{\prime}_{\gamma}(E^{\prime}_{\gamma})=\left(\frac{R_{\text{PH}}}{R_{\text{IS}}}\right)^{2}n^{\prime}_{\gamma,\text{PH}}(E^{\prime}_{\gamma})+n^{\prime}_{\gamma,\text{SYNC}}(E^{\prime}_{\gamma})+n^{\prime}_{\gamma,\text{UP}}(E^{\prime}_{\gamma})\
,$ (2.4)

with
$R_{\text{PH}}={\sigma_{T}\;\tilde{L}_{\text{iso}}}/({4\pi\Gamma^{3}m_{p}c^{3}})$
being the radius of the photosphere. Each spectral component is normalized at
the emission radius independently, using the fractions $x_{\text{PH}}$,
$x_{\text{SYNCH}}$ and $x_{\text{UP}}$ of $E^{\prime}_{\text{iso}}$ going into
each part of the spectrum, respectively. We have $x_{\text{PH}}=0.2$,
$x_{\text{SYNCH}}=1.6\times 10^{-4}$, and $x_{\text{UP}}=0.002$ [33]. The
spectral component $n^{\prime}_{\gamma,\text{PH}}(E^{\prime}_{\gamma})$ is
given by Eq. (2.3). The spectral components
$n^{\prime}_{\gamma,\text{SYNC}}(E^{\prime}_{\gamma})$ and
$n^{\prime}_{\gamma,\text{UP}}(E^{\prime}_{\gamma})$ are joint broken power
laws, describing synchrotron self-Compton emission of charged particles in the
fast-cooling regime [70]:

$n^{\prime}_{\gamma}(E^{\prime}_{\gamma})=C\begin{cases}\displaystyle\left(\frac{E^{\prime}_{\gamma}}{E^{\prime}_{\gamma,c}}\right)^{-2/3}&\text{if
}E^{\prime}_{\gamma,\min}\leq E^{\prime}_{\gamma}\leq E^{\prime}_{\gamma,c}\\\
\displaystyle\left(\frac{E^{\prime}_{\gamma}}{E^{\prime}_{\gamma,c}}\right)^{-3/2}&\text{if
}E^{\prime}_{\gamma,c}\leq E^{\prime}_{\gamma}\leq E^{\prime}_{\gamma,p}\ .\\\
\displaystyle\left(\frac{E^{\prime}_{\gamma,p}}{E^{\prime}_{\gamma,c}}\right)^{-3/2}\;\left(\frac{E^{\prime}_{\gamma}}{E^{\prime}_{\gamma,p}}\right)^{-(k+2)/2}&\text{if
}E^{\prime}_{\gamma,p}\leq E^{\prime}_{\gamma}\leq E^{\prime}_{\gamma,\max}\\\
\end{cases}$ (2.5)

For the synchrotron component, $E^{\prime}_{\gamma,i}=({3}/{2})({\hbar
e}/{m_{e}c}){\gamma^{\prime}_{i}}^{2}B^{\prime}$ where $i=p,c$,
$\gamma^{\prime}_{p}=({m_{p}}/{m_{e}})\times({k_{e}-2})/({k_{e}-1})\;\varepsilon_{\text{d}}\,\varepsilon_{e}$
and
$\gamma^{\prime}_{c}=({3\;m_{e}\Gamma})/({4m_{p}n^{\prime}_{b}x_{\text{PH}}\sigma_{T}R_{\text{IS}}})$
[70], $\sigma_{T}=6.65\times 10^{-25}$ cm2 is the Thompson cross-section, and
$n^{\prime}_{b}={\tilde{L}_{\text{iso}}}/({4\pi{R_{\text{IS}}}^{2}m_{p}c^{3}\Gamma^{2}})$
is the total baryon density. For the up-scattered component, we use the same
$E^{\prime}_{\gamma,c}$ and $E^{\prime}_{\gamma,p}$ as for the Band spectrum.
Unless otherwise stated, the benchmark jet properties adopted for this model
are summarized in Table 1.

###  2.3 ICMART model

The ICMART model belongs to the class of Poynting flux dominated jets and
assumes that energy is dissipated at very large distance from the central
engine [40]. In the ICMART model, it is assumed that the jet is composed of
magnetized shells with constant initial magnetization $\sigma=45$ and
different bulk Lorentz factors. Similarly to the internal shock model, these
shells collide with each other. These internal shocks in the optically thick
region alter the ordered magnetic field configuration, triggering magnetic
reconnection and the release of stored magnetic energy at large radii [40].
The radius at which gamma-rays and neutrinos are emitted is independent of the
jet properties ($R_{\gamma}=R_{\text{ICMART}}=10^{15}$ cm).

For the ICMART model, the photon spectrum is also a Band function, see Eq.
(2.3), and the accelerated particle spectrum is a cut-off power law (Eq. 2.2).
However, magnetic reconnection is more efficient than relativistic shocks at
accelerating particles, meaning that a smaller power law index ($k\simeq 2.0$)
and larger microphysical parameters than the ones characteristic of the
internal shock model are considered: $\varepsilon_{\text{d}}=0.35$,
$\varepsilon_{A}=0.5$ and $\varepsilon_{e}=0.5$ [56, 57, 71, 72]. Table 1
summarizes the benchmark jet properties adopted for this model. Moreover, the
magnetic field is [40]

$B^{\prime}=\sqrt{\frac{2\tilde{L}_{\text{iso}}}{c\,\Gamma^{2}\,R_{\text{ICMART}}^{2}}\,\frac{\sigma}{\sigma+1}}\
.$ (2.6)

##  3 Acceleration and cooling timescales

We present in this section the acceleration and cooling timescales used to
compute the maximum energy ($E^{\prime}_{A,\max}$) up to which a nucleus of
mass number $A$ and charge $Z$ inside the jet can be accelerated (see Eq.
2.2). We highlight that any equation presented in the following section is
also valid for protons and neutrons, which correspond to $(A=1,~{}Z=1)$ and
$(A=1,~{}Z=0)$ respectively.

The inverse of the acceleration timescale is

$\tau^{\prime-1}_{\text{acc}}=\frac{c\,Ze\,B^{\prime}}{E^{\prime}_{A}}\ ,$
(3.1)

with $E^{\prime}_{A}$ being the energy of the nucleus ${}^{A}_{Z}X$. We need
to compare $\tau^{\prime-1}_{\text{acc}}$ with the inverse of the total
cooling timescale given by

$\tau^{\prime-1}_{\rm{cool}}=\sum_{i}\tau^{\prime-1}_{i}\ ;$ (3.2)

$\tau_{i}$ stands for the adiabatic cooling, synchrotron, inverse Compton,
Bethe-Heitler pair creation, and collisions with thermal protons [73, 74]:

$\displaystyle\tau^{\prime-1}_{\text{ad}}$ $\displaystyle=\frac{c\Gamma}{R}\
,$ (3.3) $\displaystyle\tau^{\prime-1}_{\text{synch}}$
$\displaystyle=\frac{4Z^{4}\,\sigma_{T}m_{e}^{2}E^{\prime}_{A}B^{\prime
2}}{8\pi\times 3m_{A}^{4}c^{3}}\ ,$ (3.4)
$\displaystyle\tau^{\prime-1}_{\text{IC}}$
$\displaystyle=\frac{3(m_{e}c^{2})^{2}\,Z^{2}\sigma_{T}c}{16{\gamma^{\prime}_{A}}^{2}({\gamma^{\prime}_{A}}-1)\beta^{\prime}_{A}}\int_{E^{\prime}_{\gamma,\min}}^{E^{\prime}_{\gamma,\max}}F(E^{\prime}_{\gamma},\gamma^{\prime}_{A})\;n^{\prime}_{\gamma}(E^{\prime}_{\gamma})\frac{dE^{\prime}_{\gamma}}{{E^{\prime}_{\gamma}}^{2}}\
,$ (3.5) $\displaystyle\tau^{\prime-1}_{\text{BH}}$
$\displaystyle=\frac{\alpha
r_{e}^{2}cZ^{2}(m_{e}c^{2})^{2}}{E^{\prime}_{A}}\int_{2}^{\frac{2\gamma^{\prime}_{A}E^{\prime}_{\gamma,\max}}{m_{e}c^{2}}}n^{\prime}_{\gamma}\left(\frac{\epsilon^{\prime}_{\gamma}}{2\gamma^{\prime}_{A}}\right)\frac{\varphi(\epsilon^{\prime}_{\gamma})}{{\epsilon^{\prime}_{\gamma}}^{2}}d\epsilon^{\prime}_{\gamma}\
,$ (3.6) $\displaystyle\tau^{\prime-1}_{Ap}$
$\displaystyle=cn^{\prime}_{p}\;\sigma_{Ap}(E^{\prime}_{A})\ ,$ (3.7)

and $\tau^{\prime-1}_{A\gamma}\equiv R^{\prime}_{A\gamma}$ is the
photohadronic interaction rate:

$R^{\prime}_{A\gamma}=\frac{c}{2{\gamma^{\prime}_{A}}^{2}}\int_{\frac{E^{\mathrm{th}}_{\gamma,r}}{2\gamma^{\prime}_{A}}}^{+\infty}\frac{n^{\prime}_{\gamma}(E^{\prime}_{\gamma})}{{E^{\prime}_{\gamma}}^{2}}\int_{E^{\mathrm{th}}_{\gamma,r}}^{2\gamma^{\prime}_{A}E^{\prime}_{\gamma}}E_{\gamma,r}\;\sigma_{A\gamma}(E_{\gamma,r})\;dE_{\gamma,r}\;dE^{\prime}_{\gamma}\
.\ $ (3.8)

In the equations above
$\epsilon^{\prime}_{\gamma}=E^{\prime}_{\gamma}/(m_{e}c^{2})$, $E_{\gamma,r}$
is the energy of the photon in the reference frame of the nucleus,
$\gamma^{\prime}_{A}=E^{\prime}_{A}/(m_{A}\,c^{2})$,
$r_{e}=e^{2}/(m_{e}\,c^{2})\approx 2.82\times 10^{-13}$ cm is the classical
radius of the electron, $\alpha\approx 1/137$ is the fine structure constant,
and $n^{\prime}_{p}=\tilde{L}_{\mathrm{iso}}/(4\pi
R_{\mathrm{IS}}^{2}m_{p}c^{3}\Gamma^{2})$ is the density of thermal protons in
the jet. The functions $\varphi(\epsilon^{\prime}_{\gamma})$ and
$F(E^{\prime}_{\gamma},\gamma^{\prime}_{A})$ model the Bethe-Heitler and
inverse Compton processes and are respectively defined in Eqs. (9.36) and
(9.37) of Ref. [74] and Eqs. (13)–(15) of Ref. [75]. The hadronic cross-
section is $\sigma_{Ap}$; we follow Ref. [76] for the proton-proton cross-
section and Ref. [77] for the extension of such cross section to any nucleus.
The total inelastic photohadronic cross-section is $\sigma_{A\gamma}(E_{r})$
and is modeled as described in Appendix A. Note that we rely on the
photohadronic interaction rate instead of the photohadronic cooling timescale
since any isotope undergoing a photohadronic interaction is
photodisintegrated, leaving behind a new isotope which is then consistently
treated as another nuclear species in our Monte-Carlo algorithm (see Sec. 5).
This approach warrants the conservation of the nucleon number and energy. For
the sake of simplicity, the set of cooling timescales and photohadronic
interaction rate is dubbed “cooling timescale” hereafter.

For illustrative purposes, Fig. 1 shows a comparison among the cooling and
acceleration timescales obtained for protons and iron, assuming a jet with
characteristic properties summarized in Table 1 for the internal shock model.
The maximum energy at which $p$ or 56Fe can be accelerated is given by the
crossing between the acceleration and total cooling timescales. The energy for
iron useful to compute the neutrino emission is the energy per nucleon. In
this example, $E^{\prime}_{p,\max}=1.9\times 10^{9}$ GeV and
$E^{\prime}_{\text{Fe},\max}/A_{\text{Fe}}=4.2\times 10^{8}$ GeV, i.e.
acceleration of iron is slightly less efficient than proton acceleration. This
result is coherent with the fact that photonuclear interactions add an
important contribution to the radiative processes associated to iron, while
this is not the case for protons.

Figure 1: Acceleration (Eq. 3.1), inverse cooling (Eqs. 3.2–3.7) timescales,
and photohadronic interaction rate (Eq. 3.8) for protons (on the left) and
iron (on the right) as functions of the comoving proton energy and normalized
iron energy respectively. We consider a jet with characteristic properties as
in Table 1 for the internal shock model. The maximum comoving energy is
obtained computing the energy for which
$\tau^{\prime-1}_{\rm{cool}}=\tau_{\rm{acc}}^{\prime-1}$ (cf. thick black and
red curves). One can see that the normalized maximum energy for iron is
smaller than the proton one because of photonuclear interactions contributing
to the total cooling of iron and being irrelevant for protons.

##  4 Neutrino production from proton and neutron photohadronic interactions

Neutrinos can be produced thanks to photohadronic interactions of protons and
neutrons or nuclei. In this section, we focus on $p\gamma$ interactions,
following Refs. [78, 50]. We adopt the same approach for what concerns
$n\gamma$ interactions after applying charge symmetry on the mesons. We also
illustrate how to compute the neutrino fluence expected at Earth.

###  4.1 Photohadronic interactions

Figure 2: Flow chart illustrating the approach adopted to compute neutrino
production in jets loaded with protons only ($p$-jet) and jets loaded with
iron (Fe-jet).

The physical processes leading to the production of neutrinos, discussed in
this section, are illustrated in Fig. 2. Photohadronic interactions can lead
to the production of mesons, if the energy of the photon in the reference
frame of the proton is larger than $140$ MeV [79]. The production channels can
be classified as follows [78]:

Resonances: $\displaystyle\;p+\gamma\rightarrow(\Delta,N,\rho)\rightarrow
n+\pi^{+}$ (4.1) Direct production: $\displaystyle\;p+\gamma\rightarrow
n+\pi^{+}$ (4.2) Multi-pion production: $\displaystyle\;p+\gamma\rightarrow
n+\pi_{1}+\pi_{2}+...$ (4.3)

The production rate for the particle $\mathcal{M}$ (with
$\mathcal{M}=\pi^{+}$, $\pi^{-}$, $K^{+}$ or $K^{-}$) is defined as the number
of particles created per volume, per unit of time and per energy interval (in
units of cm-3 s-1 GeV-1) and it is given by:

$Q^{\prime}_{\mathcal{M}}(E^{\prime}_{\mathcal{M}})=\sum_{\alpha}\int
n^{\prime}_{p}(E^{\prime}_{p})\;\Gamma^{\alpha}_{p\gamma\rightarrow\mathcal{M}}\;\frac{d{n^{\prime\alpha}_{p\gamma\rightarrow\mathcal{M}}}}{dE^{\prime}_{p}}(E^{\prime}_{p},E^{\prime}_{\mathcal{M}})\;dE^{\prime}_{p}\
.$ (4.4)

The sum runs over all possible interactions $\alpha$, see Eqs. (4.1)–(4.3).
The seed proton density per unit energy is $n^{\prime}_{p}(E^{\prime}_{p})$
(in units of cm-3 GeV-1), while
$\Gamma^{\alpha}_{p\gamma\rightarrow\mathcal{M}}(E^{\prime}_{p})$ is the
production rate for a given interaction (in units of s-1), and
${d{n^{\prime\alpha}_{p\gamma\rightarrow\mathcal{M}}}}/{dE^{\prime}_{p}}$
represents the probability density for production of a meson with energy
$E^{\prime}_{\mathcal{M}}$ given that the incident proton has an energy
$E^{\prime}_{p}$.

In order to simplify the numerical solution of the integral above, we assume
that all pions are produced with a certain energy that does not depend on the
parent proton energy:

$\frac{d{n^{\prime\alpha}_{p\gamma\rightarrow\mathcal{M}}}}{dE^{\prime}_{p}}\simeq
M_{\mathcal{M}}^{\alpha}\times\delta(E^{\prime}_{\mathcal{M}}-\chi_{\mathcal{M}}^{\alpha}\;E^{\prime}_{p})\
,$ (4.5)

where $M_{\mathcal{M}}^{\alpha}$ is the multiplicity of the interaction and
$\chi_{\mathcal{M}}^{\alpha}$ describes the amount of energy transferred to
the produced meson. We also assume an isotropic photon distribution. Hence,
Eq. (4.4) becomes:

$Q^{\prime}_{\mathcal{M}}(E^{\prime}_{\mathcal{M}})=\sum_{\alpha}n^{\prime}_{p}\left(\frac{E^{\prime}_{\mathcal{M}}}{\chi_{\mathcal{M}}^{\alpha}}\right)\;\frac{m_{p}c^{2}}{E^{\prime}_{\mathcal{M}}}M_{\mathcal{M}}^{\alpha}\;\int_{E_{\gamma,r}^{\mathrm{th},\alpha}/2}^{+\infty}f^{\alpha}(y)\;n^{\prime}_{\gamma}\left(y\frac{\chi_{\mathcal{M}}^{\alpha}m_{p}c^{2}}{E^{\prime}_{\mathcal{M}}}\right)dy\
,$ (4.6)

with

$f^{\alpha}(y)=\frac{1}{2y^{2}}\;\int_{E_{\gamma,r}^{\mathrm{th},\alpha}}^{2y}\;E_{\gamma,r}\;\sigma^{\alpha}(E_{\gamma,r})\;dE_{\gamma,r}\
.$ (4.7)

The parameters characterizing each interaction $\alpha$ (namely
$\chi_{\mathcal{M}}^{\alpha}$, $M_{\mathcal{M}}^{\alpha}$, $f^{\alpha}(y)$ and
$E_{\gamma,r}^{\mathrm{th},\alpha}$) can be found in Ref. [78].

###  4.2 Neutrino production

The mesons ($\pi^{+}$, $\pi^{-}$, $K^{+}$ and $K^{-}$) produced in the
reactions illustrated above can decay into neutrinos through the following
channels:

$\displaystyle\pi^{+}$ $\displaystyle\rightarrow$
$\displaystyle\mu^{+}+\nu_{\mu}\rightarrow
e^{+}+\nu_{e}+\bar{\nu}_{\mu}+\nu_{\mu}$ (4.8) $\displaystyle K^{+}$
$\displaystyle\rightarrow$ $\displaystyle\mu^{+}+\nu_{\mu}\rightarrow
e^{+}+\nu_{e}+\bar{\nu}_{\mu}+\nu_{\mu}$ (4.9)

and similar reactions hold for the antiparticle conjugates. Other channels of
kaon decay are neglected since their branching ratio is smaller and they would
generate lower energy neutrinos.

The neutrino production for the flavor $\nu_{\beta}$ from the decay of
$\mathcal{M}$ is given by [78]:

$Q^{\prime}_{\nu_{\beta}}(E^{\prime}_{\nu_{\beta}})=\int_{E^{\prime}_{{\nu_{\beta}}}}^{+\infty}Q^{\prime}_{\mathcal{M}}(E^{\prime}_{\mathcal{M}})\left[1-\exp{\left(-\frac{\tau^{\prime}_{\text{cool}}}{\tau^{\prime}_{l}}\right)}\right]\frac{1}{E^{\prime}_{\mathcal{M}}}F_{\mathcal{M}\rightarrow\nu_{\beta}}\left(\frac{E^{\prime}_{\nu_{\beta}}}{E^{\prime}_{\mathcal{M}}}\right)dE^{\prime}_{\mathcal{M}}\
.$ (4.10)

with $\nu_{\beta}$ representing one neutrino (or antineutrino,
$\bar{\nu}_{\beta}$) species. The exponential term in the equation above takes
into account the cooling of the decaying particle, with
$\tau^{\prime}_{\text{cool}}$ and $\tau^{\prime}_{l}$ being the total cooling
timescale and the mean lifetime of $\mathcal{M}$, respectively. The function
$F_{\mathcal{M}\rightarrow\nu_{\beta}}\left({E^{\prime}_{\nu_{\beta}}}/{E^{\prime}_{\mathcal{M}}}\right)$
encapsulates the reaction kinematics and is defined as in Ref. [50]. In
particular, the muon decay into neutrinos depends on the muon helicity, as
illustrated in Ref. [50].

###  4.3 Neutrino fluence at Earth

En route to the Earth, neutrinos change their flavor, with transition and
survival probabilities defined as follows [80]:

$\displaystyle
P_{\nu_{e}\rightarrow\nu_{\mu}}=P_{\nu_{\mu}\rightarrow\nu_{e}}$
$\displaystyle=$ $\displaystyle\frac{1}{4}\;\sin^{2}{2\,\theta_{12}}\ ,$
(4.11) $\displaystyle P_{\nu_{\mu}\rightarrow\nu_{\mu}}$ $\displaystyle=$
$\displaystyle\frac{1}{8}\;(4-\sin^{2}{2\,\theta_{12}})\ ,$ (4.12)
$\displaystyle P_{\nu_{e}\rightarrow\nu_{e}}$ $\displaystyle=$ $\displaystyle
1-\frac{1}{2}\;\sin^{2}{2\,\theta_{12}}\ ,$ (4.13)

with $\theta_{12}=33.5^{\circ}$ [81]. Transition and survival probabilities
are the same assuming antineutrinos instead of neutrinos. Hence the neutrino
fluence $F_{\nu_{\beta}}$ (in units of cm-2 GeV-1) is:

$F_{\nu_{\beta}}(E_{\nu_{\beta}})=V^{\prime}_{s}\;t_{\text{dur}}\;\frac{(1+z)^{2}}{4\pi
d_{L}(z)^{2}}\sum_{\alpha}P_{\nu_{\alpha}\rightarrow\nu_{\beta}}\,Q^{\prime}_{\nu_{\alpha}}\left(\frac{E_{\nu_{\beta}}(1+z)}{\Gamma}\right)\
.$ (4.14)

Within a flat $\Lambda$CDM cosmology, the luminosity distance is given by

$d_{L}\;(z)=\frac{c\;(1+z)}{H_{0}}\int_{0}^{z}\frac{dz^{\prime}}{\sqrt{\Omega_{\Lambda}+\Omega_{M}(1+z^{\prime})^{3}}}\
,$ (4.15)

where $H_{0}=67.4$ km s-1 Mpc-1, $\Omega_{M}\simeq 0.315$, and
$\Omega_{\Lambda}\simeq 0.685$ [82].

##  5 Neutrino production from nuclear photohadronic interactions

The outflows of GRB jets may contain a certain fraction of nuclei [30]. Nuclei
can be synthetized along the jet, loaded at the base of the jet, or produced
in the stellar ejecta and entrained in the jet. However, photodisintegration
from high-energy photons and spallation from protons and neutrons may
disintegrate such nuclei. In fact, Ref. [31] concluded that the survival of
nuclei up to the region of neutrino production is unlikely in HL-GRBs in the
internal shock and photospheric models, but nuclei could easily survive in
magnetized jets or LL-GRBs due to the smaller photon density. In this section,
we introduce the method adopted to model jets containing nuclei through
photohadronic interactions.

###  5.1 Monte Carlo simulation of nuclear cascades

In order to model the neutrino emission, interactions of nuclei in the jet
should be taken into account. To this purpose, Ref. [32] numerically solved a
system of partial differential equations to compute the spectrum of any
nucleus as a function of time and energy; while Ref. [47] relied on the SOPHIA
event generator to compute energy loss and photodisintegration processes of
accelerated nuclei. In the following, we first model nuclear disintegration
and nuclear decay. Then, along the lines of Ref. [47], we outline a Monte-
Carlo algorithm that generates nuclear cascades, the Nuclear Cascade Monte
Carlo (NCMC). The steps leading to the production of neutrinos, presented in
this section, are illustrated in Fig. 2.

####  5.1.1 Nuclear photodisintegration

Nuclei undergo photon-nucleus and nucleus-nucleus collisions. However, since
the photon density in the jet is several orders of magnitude larger than the
density of nuclei, we only focus on photon-nucleus interactions and neglect
any nucleus-nucleus collision.

Photonuclear interactions are fundamentally different than $p\gamma$ and
$n\gamma$ interactions. Indeed, nuclei are made of bonded nucleons. Photons
with enough energy can break the nuclear bond, leading to the emission of one
or several nucleons. We should distinguish between two ranges for the energy
of the incident photon in the rest frame of the nucleus: The Giant Dipole
Resonance (GDR) with
$E_{\gamma,r}\in[E^{\mathrm{th}}_{\gamma,r},140\;\text{MeV}]$ and the
photomeson (PM) range for $E_{\gamma,r}>140\;\text{MeV}$ [79]. The GDR
threshold energy $E^{\mathrm{th}}_{\gamma,r}$ varies as a function of $(Z,A)$,
but it is approximately $10$ MeV [83].

Data from photonuclear reactions are relatively scarce compared to
$(p,n)\gamma$ interactions. However, a few models for the cross-section and
branching ratios for any nucleus with mass number $A\leq 56$ exist. Major
improvements on the estimation of these cross-sections have been made thanks
to the Monte-Carlo event generation software Geant4 [84]. We adopt fits from
Refs. [85] and [79] for the GDR and PM ranges of the cross section. Regarding
the branching ratios of these processes, we draw from Refs. [86] and [79] for
the giant dipole resonance and the photomeson energy ranges, as illustrated in
Appendix A (this appendix presents an extensive description of the modeling of
the cross section, the branching ratios and the adopted parameters).

####  5.1.2 Nuclear decay

Any stable nucleus that looses nucleons may become a radioactive isotope. We
have refined the NCMC by taking into account both lifetimes and radioactive
decay mechanisms of any isotope of mass number $A\leq 56$. To this purpose, we
have retrieved nuclear decay data from Ref. [87] (cf. Appendix A.2).

When a nucleus undergoes photodisintegration, we compute its new isotope
lifetime at rest $t_{l}$, randomly extracting a time from a distribution
following the exponential probability density:

$p(t_{l})=\frac{1}{\bar{t_{l}}}\,e^{-\frac{t_{l}}{\bar{t_{l}}}}\ ,$ (5.1)

with $\bar{t_{l}}$ being the mean decay time of the isotope at rest. The
lifetime in the reference frame of the jet is

$t^{\prime}_{l}=\frac{E^{\prime}_{A}}{m_{A}c^{2}}\;t_{l}\ ,$ (5.2)

where $t^{\prime}_{l}$ is the comoving lifetime of the isotope in the jet
frame. Assuming that the isotope is created at $t^{\prime}=0$, if
photodisintegration does not happen at $t^{\prime}_{l}$, the isotope decays as
illustrated in Appendix A.2.

Following the computed branching ratios, we can see that it is unlikely for a
nucleus to stray from the valley of stability, as shown in Fig. 3. Since we
are interested in the high-energy tail of the spectral energy distribution of
nuclei, the lifetime of radioactive nuclei in the comoving frame is large (see
Eq. 5.2). This implies that the nucleus is relatively stable and neutrino
emission from $\beta^{+}/\beta^{-}$ processes from radioactive decays is
negligible with respect to photohadronic interactions [32]. The only exception
is the $\beta^{-}$ decay of free neutrons, which we take into account given
that the related radioactive isotope can be abundant, see for example Fig. 3.

Figure 3: Contour plot of the fraction of the total energy (initially injected
in 56Fe) carried by the different isotopes in the plane spanned by $N$ and $Z$
for our benchmark GRB (cf. Table 1 for Fe-jet). The results for the internal
shock, photospheric and ICMART models are plotted from left to right,
respectively. In order to guide the eye, 56Fe, 28Si, 16O, and $p$ are marked
with red stars. The internal shock model does not favor efficient nuclear
cascades with respect to the other two models. The ICMART model leads to the
most efficient cascades, which in turn determine the largest neutrino emission
as shown in Fig. 4.

####  5.1.3 Monte Carlo algorithm and particle spectra

Accelerated nuclei can be injected anytime during the interaction time
interval,
$t^{\prime}_{\text{shock}}=t^{\prime}_{\text{dur}}/N_{\text{shock}}$. This
implies that we can compute a random interaction time for each Monte-Carlo
event as $t^{\prime}_{\text{eff}}=t^{\prime}_{\text{shock}}\;U(0,1)$, with
$U(0,1)$ being the uniform probability density function between $0$ and $1$.
For each Monte-Carlo event, $t^{\prime}_{\text{eff}}$ is split into $N$ equal
time steps $\Delta t^{\prime}=t^{\prime}_{\text{eff}}/N$. Hence, the
probability of a photonuclear interaction for a nucleus ${}^{A}_{Z}X$ of
energy $E^{\prime}_{A}$, during the time interval $\Delta t^{\prime}$ in the
comoving frame of the jet is:

$p_{A\gamma}(E^{\prime}_{A},A,Z)=R^{\prime}_{A\gamma}(E^{\prime}_{A},A,Z)\;\Delta
t^{\prime}\ ,$ (5.3)

with $R^{\prime}_{A\gamma}$ being the photonuclear interaction rate computed
according to Eq. (3.8). At each time step, nuclei can interact or not with
photons according to Eq. (5.3). If an interaction occurs, one of the physical
processes described in Sec. 5.1.1 is randomly selected according to the
branching ratios. In parallel, nuclei loose energy at each time step,
following the adiabatic cooling; excluding the photodisintegration cooling
which is already taken into account through the NCMC, adiabatic processes
always dominate the total cooling timescale, as shown in Fig. 1. If no
interaction takes place, when the nucleus is unstable it radioactively decays,
as described in Sec. 5.1.2.

We assume that, after a photonuclear interaction that generates
$(N_{p},N_{n})$ protons and neutrons, only nucleons and no nuclei can be
emitted. For instance, the loss of $2$ protons and $2$ neutrons generates $2$
free protons and $2$ free neutrons, and not ${}_{2}^{4}$He. The energy of the
produced nucleon $\mathcal{N}=(p,n)$ is
$E^{\prime}_{\mathcal{N}}={E^{\prime}_{A}}/{A}$, with $E^{\prime}_{A}$ being
the energy of the parent particle ${}^{A}_{Z}X$ in the comoving frame of the
jet.

For any accelerated nuclear species, the spectral energy distribution
$n^{\prime}_{A}(E^{\prime}_{A})$ (in units of cm-3 GeV-1) of the seed
population is adopted to draw the injected nuclear species in the NCMC. For
simplicity, we only inject pure isotopes. In order to do this, we randomly
select the initial energy of the seed nuclei following a uniform probability
in logarithmic scale. Then, for any given energy, each seed nucleus is
associated to a certain density computed following
$n^{\prime}_{A}(E^{\prime}_{A})$. The nuclear cascade is then generated and we
retrieve the density and final energy of each isotope eventually created. We
then add each contribution to the final spectra. For more details, we refer
the reader to Appendix A.3.

###  5.2 Nuclear photohadronic interactions

It is now possible to compute the neutrino production associated to each
isotope, while for the production of mesons from protons and neutrons, we rely
on Sec. 4. For nuclei with mass number $A\leq 56$, the meson production
$Q^{\prime}_{\mathcal{M}}(E^{\prime}_{\mathcal{M}})$ (in units of cm-3 s-1
GeV-1) is given by [78]:

$Q^{\prime}_{\mathcal{M}}(E^{\prime}_{\mathcal{M}})=\sum_{\Delta
E}n^{\prime}_{A}\left(\frac{E^{\prime}_{\mathcal{M}}}{\chi_{\mathcal{M}}^{\Delta
E}}\right)\;\frac{m_{A}c^{2}}{E^{\prime}_{\mathcal{M}}}M_{\mathcal{M}}^{\Delta
E}\;\int_{\Delta E}f^{\Delta
E}(y)\;n^{\prime}_{\gamma}\left(y\frac{\chi_{\mathcal{M}}^{\Delta
E}m_{A}c^{2}}{E^{\prime}_{\mathcal{M}}}\right)dy\ ,$ (5.4)

where $\mathcal{M}$ is any meson (i.e., pions and kaons) and the sum runs over
$\Delta E$, which is an energy interval where the multiplicity and the
fraction of the energy transferred to produced mesons are assumed to be
constant (see Appendix B).

Each interaction type is characterized by three parameters:
$\chi_{\mathcal{M}}^{\Delta E}$ represents the amount of energy transferred
from the seed nucleus to the newly created meson
($E^{\prime}_{\mathcal{M}}=\chi_{\mathcal{M}}^{\Delta E}E^{\prime}_{A}$); the
multiplicity, $M_{\mathcal{M}}^{\Delta E}$, describes how many mesons are
created on average; $f^{\Delta E}(y)$ is computed as follows:

$f^{\Delta E}(y)=\frac{1}{2y^{2}}\;\int_{E_{\gamma,r}^{\mathrm{th},\Delta
E}}^{2y}\;E_{\gamma,r}\;\sigma_{A\gamma\rightarrow\mathcal{M}}(E_{\gamma,r})\;dE_{\gamma,r}\
,$ (5.5)

with $E_{\gamma,r}^{\mathrm{th},\Delta E}$ and
$\sigma_{A\gamma\rightarrow\mathcal{M}}(E_{\gamma,r})$ being the threshold
energy and the inelastic cross-section of photohadronic interactions for meson
production, respectively.

In order to compute the cross section, previous work (e.g. Ref. [32]) relied
on the Single Particle Model, according to which
$\sigma_{A\gamma\rightarrow\mathcal{M}}(E_{\gamma,r})\simeq
A\;\sigma_{p\gamma\rightarrow\mathcal{M}}(E_{\gamma,r})$. However, since Ref.
[79] showed that the Single Particle Model can be inaccurate in the photomeson
regime, we adopt the empirical model proposed in Ref. [79], which provides
more accurate results for nuclei up to 56Fe. The function $f^{\Delta E}(y)$ is
then computed using the empirical model total cross-section in the photomeson
range (see Appendix B for details). The branching ratio for the photohadronic
production is extracted from estimations of the cross-sections from Ref. [79].
Due to the lack of data, we model $M_{\mathcal{M}}^{\Delta E}$ and
$\chi_{\mathcal{M}}^{\Delta E}$ adapting the method employed to compute
$p\gamma$ interactions in Ref. [78]. The average values of
$M_{\mathcal{M},p\gamma}^{\Delta E}$ and $\chi_{\mathcal{M},p\gamma}^{\Delta
E}$ are taken into account for each $\Delta E$, and we define
$\chi_{\mathcal{M}}^{\Delta E}=\chi_{\mathcal{M},p\gamma}^{\Delta E}\;/A$,
$M_{\mathcal{M}}^{\Delta E}=M_{\mathcal{M},p\gamma}^{\Delta E}$.

For the sake of simplicity, we assume that each nucleus contains an equal
number of protons and neutrons. This approximation is reasonable given that it
is unlikely for a nucleus to stray away from the valley of stability (which
approximately corresponds to $N=Z=A/2$ for $A\leq 56$). This symmetry enables
us to get comparable production of $\mathcal{M}$ and $\bar{\mathcal{M}}$,
reducing the computational time. For details on the inputs adopted for the
different parameters, we refer the interested reader to Appendix B. Finally,
the decay of each meson into neutrinos is computed as in Sec. 4.2 and the
neutrino fluence expected at Earth follows Sec. 4.3, respectively.

##  6 Impact of the jet composition on the neutrino emission

In this section, we present our findings on the dependence of the neutrino
signal on the jet composition. We investigate how the jet composition impacts
the neutrino emission for the internal shock, photospheric and ICMART jet
models. We then explore how the jet composition affects the neutrino emission
across the allowed jet parameter space.

###  6.1 The role of the jet composition in kinetic and Poynting flux
dominated jets

In order to facilitate a comparison, we consider the neutrino fluence obtained
for a jet loaded with protons ($p$-jet) and a jet loaded with 56Fe (Fe-jet)
for the internal shock, photospheric and ICMART jet models. Jets with mixed
composition or loaded with intermediate isotopes would have a neutrino fluence
falling in between these two extreme cases. We also note that it is highly
unlikely to have a jet mostly composed by iron, yet we investigate this case
in order to explore the largest variability of the expected neutrino signal
according to the jet composition. We otherwise assume identical jet properties
(cf. Table 1); in particular, the total energy fraction going into accelerated
protons is the same as the one going into accelerated iron. This implies that
the same power law index is chosen for both accelerated protons and iron
nuclei (see Sec. 6.2 for an extensive discussion about this hypothesis).

We note that the survival of nuclei up to the region of neutrino production
may be unlikely for our benchmark model, especially for kinetic dominated jets
[31]. Yet, we rely on this benchmark jet to assess the impact of the jet
composition for illustrative purposes.

For the internal shock model, in order to test our NCMC approach, we compare
our results with the three benchmarks cases presented in Ref. [32]. In all
cases, we find excellent agreement (results not shown here). We also find
comparable neutrino fluences using the Single Particle Model, yet we prefer to
employ the more accurate photohadronic model in the rest of this paper (see
Sec. 5.2).

Figure 4: Left: Spectral energy density of nuclei as a function of the
particle energy for our benchmark GRB (cf. Table 1) in the case of Fe-jet, for
the internal shock, photospheric, and ICMART models, from top to bottom
respectively. Right: $\nu_{\mu}+\bar{\nu}_{\mu}$ neutrino fluence as a
function of the neutrino energy for our benchmark GRB, in the case of a
$p$-jet (indigo curve) and Fe-jet (the ochre curve represents the total
neutrino emission). The relative contributions to the neutrino fluence from
iron, all other lighter nuclei and nucleons are plotted in blue, green and red
respectively (cf. spectra of the corresponding parent particles in the left
panels). If the density of iron nuclei is larger than the one of all other
nuclei, neutrinos are produced mostly through iron photohadronic interactions.
Since meson production from nuclei is less efficient than the one from protons
and neutrons, the neutrino fluence for Fe-jet is smaller than for $p$-jet,
independent of the jet model. Nuclear cascades are more efficient for the
photospheric and ICMART models than for the internal shock one; this leads to
larger neutrino fluence than in the internal shock model for Fe-jet.

Figure 4 shows the spectral energy density of nuclei (left panels) and the
neutrino fluence (right panels) for the internal shock model, the photospheric
model and the ICMART model, from top to bottom respectively. Comparing Fe-jet
(ochre line) with the $p$-jet (indigo line), we can see that the neutrino
fluence is lower for Fe-jet than for the $p$-jet up to one order of magnitude,
independent of the jet model. For the internal shock case, the density of
nuclei is larger than the one of nucleons (cf. blue/green vs. red curves in
the left panels), but neutrino production from nuclei is comparable to the one
driven by protons/neutrons. Indeed, since the meson production from nuclei is
not so efficient as the one from protons or neutrons, the neutrino fluence is
smaller than what could be expected given the associated nuclei densities. The
shape of the neutrino fluence is overall the same for jets, independent of the
composition. In particular, the maximum of
$E^{2}_{\nu}F_{\nu_{\mu}+\bar{\nu}_{\mu}}$ is achieved at about the same
energy $E^{M}_{\nu}\simeq 5\times 10^{7}$ GeV.

For the photospheric model, cascades are more efficient than in the internal
shock model, due to the larger photon luminosity which makes nuclear cascades
more efficient, leading to larger proton and neutron densities (cf. Fig. 3
which displays the fraction of the total energy carried by lighter isotopes
and nucleons for Fe-jet). The ICMART model shows efficient nuclear cascades,
as shown in Fig. 3, which lead to the largest neutrino fluence across all
three jet models. However, compared to the photospheric model, the maximum
fluence is achieved at the same energy for both $p$\- and Fe-jets, and no
sharp cutoff at high energies is observed. This is explained because the
crossing between the acceleration and cooling timescales happens at larger
energies in the ICMART model since the magnetic field is stronger than in the
photospheric model. Hence, the acceleration processes are more efficient and
the maximum energy per nucleon inside an iron nucleus is similar to the
maximum energy found for protons.

Our findings are in general agreement with the ones presented in Fig. 7 of
Ref. [32], our internal shock case would correspond to their “empty cascade”
scenario, while our photospheric and ICMART models would correspond to their
“populated cascade” scenario. However, our results differ for their
“inefficient cascade” scenario because of the empirical cross-section model
adopted in this work (vs. the Single Particle Model adopted in Ref. [32]).

Figure 5: Nuclear cascade efficiency in the plane spanned by the radius
($R_{\mathrm{IS}}$) of the interaction region and the isotropic luminosity
($\tilde{L}_{\rm{iso}}$) for the photospheric model. Left: The contours
represent the region where the related isotopes (iron in black, any other
lighter nucleus in magenta, and protons and neutrons in yellow) contribute the
most to the total density. Right: Same as the left panel, but it represents
the dominant contribution to the total neutrino fluence. Dashed orange lines
indicate isocontours of $E_{\nu}^{2}F_{\nu_{\mu}+\bar{\nu}_{\mu}}$ (in units
of GeV cm-2). For low luminosity and/or large radius, nuclear cascades are not
efficient and iron dominates both the energy density of nuclei and the
neutrino fluence. As the luminosity increases and/or the radius decreases,
nuclear cascades are more efficient, and intermediate isotopes can dominate
the nuclear energy density with a resulting larger neutrino fluence. However,
protons and neutrons are responsible for most of the total neutrino fluence
despite their lower density.

Figure 5 highlights the variation of nuclear cascade efficiency in the plane
spanned by the radius of interaction and the GRB luminosity. For illustrative
purposes, the photospheric model is adopted as a proxy since the overall
structure would be the same for the internal shock model as well, while the
ICMART model would correspond to a slice in this plot drawn for constant
radius. For low luminosity and/or large radius, nuclear cascades are not
efficient and iron dominates both the particle density and the neutrino
fluence. As discussed earlier (see Fig. 4), we expect that the total neutrino
fluence in this region is approximately reduced by one order of magnitude
compared to a $p$-jet. As the luminosity increases and/or the radius
decreases, nuclear cascades become more efficient, and we reach a region where
intermediate isotopes dominate the nuclear energy density. However, except for
a thin region (cf. the right panel of Fig. 5, magenta region), protons and
neutrons are responsible for most of the total neutrino fluence despite their
lower density. This is because nucleons are more efficient at producing
neutrinos from photohadronic interactions than nuclei.

###  6.2 Dependence of the neutrino fluence on the injection spectral index of
the parent nuclei and the jet Lorentz factor

We now extend our investigation of the dependence of the neutrino emission on
the jet composition beyond our benchmark GRB jet, considering the variation of
the neutrino fluence as a function of the bulk Lorentz factor ($\Gamma$) and
the injection spectral index of the seed accelerated particle spectra ($k$).
We choose to compute the neutrino fluence for $\Gamma=10$–$1000$. Regarding
the power law index, for relativistic shocks (internal shock and photospheric
models), theoretical work predicts $k\simeq 2.2$ [88] and particle-in-cell
simulations find $k\simeq 2.0$–$2.5$ [58, 59]; for magnetized jets, particle-
in-cell simulations suggest $k\simeq 2.0$ [60, 61]. It is crucial to highlight
that most of the existing work does not consider populations of heavier nuclei
because particle-in-cell simulations with heavier nuclei are computationally
intensive. In Ref. [59], it has been found that the power law index does not
depend strongly on the mass of the ion for ion masses ranging from
$m_{i}=25\;m_{e}$ to $m_{i}=1600\;m_{e}\simeq m_{p}$. Thus, we extrapolate
from this result, and consider that the power law index for an iron population
would be similar to the one of a proton population and explore the efficiency
of neutrino production for $k=2.0$–$3.0$ for $p$\- and Fe-jets.

Figure 6: Dependence of the neutrino fluence on the injection spectral index
of the accelerated parent nucleus/nucleon spectra. Top:
$\nu_{\mu}+\bar{\nu}_{\mu}$ neutrino fluence as a function of the neutrino
energy for $k=2.0$–$3.0$, using our benchmark Fe-jet (cf. Table 1) for the
internal shock, photospheric, and ICMART jet models; the other jet parameters
are defined as in Table 1, in particular $\Gamma=300$ is fixed. Bottom: Ratio
between the Fe-jet neutrino fluence and the $p$-jet one as a function of the
neutrino energy. A smaller $k$ allows for a larger neutrino fluence, because
it favors larger densities of seed nuclei; such trend holds independent of the
jet composition.

The top panels of Fig. 6 show the neutrino fluence for varying $k$ for the
internal shock, photospheric and ICMART models; all other jet parameters are
kept fixed as indicated in Table 1. Decreasing the power law index tends to
increase the neutrino fluence, because it allows for larger densities of seed
nuclei. We can also note that decreasing $k$ slightly increases the energy
$E_{\nu}^{M}$ corresponding to the maximum neutrino fluence. For our benchmark
jet, we also observe that $k$ mostly scales the neutrino fluence in a similar
fashion for $p$\- and Fe-jets, as evident from the bottom panels of Fig. 6.

Figure 7: Same as Fig. 6, but for investigating the dependence of the neutrino
fluence on $\Gamma$ ($k=2.2$ is fixed for both the internal shock and the
photospheric models, while $k=2.0$ is used for the ICMART case). Decreasing
$\Gamma$ increases the neutrino fluence, with the larger spread occurring for
the internal shock and photospheric jet models. In fact, for these models, the
interaction radius scales as $\Gamma^{2}$, which implies that for larger
$\Gamma$, internal shocks happen at a radius where both photon and seed nuclei
have lower densities. However, this is not the case for the ICMART model, for
which the radius of interaction is independent of the bulk Lorentz factor.
This explains the milder variation of the neutrino fluence on $\Gamma$ in the
internal shock model. A strong variation of the neutrino fluence as a function
of $\Gamma$ occurs according to the composition.

Figure 7 explores the impact of $\Gamma$ on the neutrino emission. For fixed
composition, from the top panels of Fig. 7, we can see that decreasing
$\Gamma$ tends to increase the neutrino fluence, in a way that depends on the
jet model. For the internal shock and photospheric models, changes in $\Gamma$
affect the radius of interaction
($R_{\text{IS}}=2\Gamma^{2}c\tilde{t}_{\text{var}}$); when $\Gamma$ decreases,
internal shocks happen at a radius where both photon and seed nuclei have
larger densities. However, this is not the case for the ICMART model, for
which the radius of interaction is independent of the bulk Lorentz factor
(internal shocks trigger magnetic reconnection events assumed to happen at a
nearly constant radius). This explains the milder dependence of the neutrino
fluence on $\Gamma$ in the internal shock model. The energy corresponding to
the maximum neutrino fluence ($E_{\nu}^{M}$) is strongly shifted to lower
energies as $\Gamma$ decreases; this is because higher photon densities are
achieved as $\Gamma$ decreases, hence the photodisintegration timescale
increases and crosses the acceleration timescale at lower energies. This
implies that the cutoff energy of the parent spectrum decreases, which also
forces the neutrino fluence to peak at lower energies. This effect is less
prominent in Poynting-flux jets than in kinetically-dominated jets since the
acceleration process is more efficient in the former thanks to the larger
magnetic fields.

The bottom panels of Fig. 7 also highlight a strong variation of the neutrino
fluence as a function of $\Gamma$ according to the composition. The ratio
between the neutrino fluence of Fe-jet and $p$-jet is a direct measurement of
the efficiency of nuclear cascades. For large $\Gamma$, this ratio saturates
near $10^{-1}$ because of the less efficient production of neutrinos from 56Fe
nuclei; then, there is a transition from inefficient to efficient cascades as
$\Gamma$ decreases due to higher photon densities. This transition is more or
less sharp depending on the model, following the dependence on the interaction
radius on $\Gamma$. For small $\Gamma$, the ratio saturates around $1$, which
is consistent with the fact that the Fe-jet with highly efficient cascades of
iron into nucleons is equivalent to the $p$-jet (disintegrated 56Fe leads to
$26$ protons and $30$ neutrons; photohadronic interactions are equivalent for
protons and neutrons as long as no distinction between neutrino and
antineutrino is considered).

Figure 8 quantifies the dependence of the neutrino energy ($E_{\nu}^{M}$) at
which the fluence peaks in Fe-jets and $p$-jets in the parameter space spanned
by $k$ and $\Gamma$. In order to avoid degeneracies in the computation of
$E_{\nu}^{\text{M}}$ for the cases where the neutrino fluence reaches a
plateau around the spectral peak, we compute the two energies
($E_{\nu}^{90\%,1}$ and $E_{\nu}^{90\%,2}$) where the fluence is $90\%$ of its
maximum. Then, the maximum energy is computed as
$\log_{10}(E_{\nu}^{M})={1}/{2}\,[\log_{10}(E_{\nu}^{90\%,1})+\log_{10}(E_{\nu}^{90\%,2})]$.
The contours of the ratio of $E_{\nu}^{M}$ for the Fe-jet and $p$-jet in Fig.
8 highlight that $E_{\nu}^{M}$ does not change much as a function of the jet
composition. However, $E_{\nu}^{M}$ has non-trivial dependencies on the jet
parameters according to the jet model. We note that the normalization of the
neutrino fluence does strongly depend on the jet composition, as illustrated
in Figs. 6 and 7.

Figure 8: Contour plots of the ratio between the neutrino energy at which the
neutrino fluence peaks for the Fe-jet and the p-jet in the plane spanned by
$k$ and $\Gamma$. From left to right, the contour plots for the internal
shock, photospheric, and ICMART models are illustrated, respectively. The
energy at which the neutrino fluence peaks does not depend on the jet
composition strongly for all considered jet models, however the normalization
of the spectral energy distribution does (cf. Figs. 6 and 7).

##  7 Dependence of the diffuse neutrino flux on the jet composition

In this work, we do not explore the detection prospects of neutrinos, since we
are interested in assessing the model dependence on the jet composition. On
the other hand, in order to gauge the variation of the neutrino signal as a
function of the jet composition for different jet models, we consider the
diffuse emission of neutrinos from HL- LL- and sGRBs. The diffuse neutrino
emission is defined as follows:

$F^{d}_{\nu_{\beta}}(E_{\nu_{\beta}})=\int_{z_{\text{min}}}^{z_{\text{max}}}\frac{c}{4\pi
H_{0}\Gamma}\frac{1}{\sqrt{\Omega_{\Lambda}+(1+z)^{3}\Omega_{M}}}\;R_{\text{GRB}}(z)\;\frac{4\pi\Gamma{d_{L}(z)}^{2}}{(1+z)^{3}}F_{\nu_{\beta}}(E_{\nu_{\beta}},\,z)\;dz\
,$ (7.1)

where $R_{\text{GRB}}(z)$ is the cosmological rate (in units of Gpc-3 yr-1),
which is different for HL-GRBs, LL-GRBs and sGRBs. For HL-GRBs, we model the
rate following Ref. [89, 90], for LL-GRBs we consider Refs. [91, 92], and for
sGRBs we follow Ref. [93]:

$\displaystyle R_{\rm{HL-GRB}}(z)$ $\displaystyle=\rho_{0,\rm{HL-
GRB}}\begin{cases}(1+z)^{2.1}&\rm{for}\;z<z^{\rm{HL-GRB}}_{\star}\\\
(1+z^{\rm{HL}}_{\star})^{2.1+0.7}\;(1+z)^{-0.7}&\rm{for}\;z>z_{\star}\end{cases}\
,$ (7.2) $\displaystyle R_{\rm{LL-GRB}}(z)$ $\displaystyle=\rho_{0,\rm{LL-
GRB}}\left[(1+z)^{-34.0}+\left(\frac{1+z}{5000}\right)^{3.0}+\left(\frac{1+z}{9}\right)^{35.0}\right]^{-1/10}\
,$ (7.3) $\displaystyle R_{\rm{sGRB}}(z)$
$\displaystyle=\begin{cases}10^{\alpha\,z^{3}-\beta\,z^{2}+\gamma\,z+\delta}&\rm{for}\;z<z^{\rm{sGRB}}_{\star}\\\
10^{-\eta\,z+\mu}&\rm{for}\;z>z^{\rm{sGRB}}_{\star}\end{cases}\;\text{Gpc${}^{-3}$~{}yr${}^{-1}$}\
,$ (7.4)

with $\rho_{0,\rm{HL-GRB}}=200$–$2000$ Gpc-3 yr-1, $\rho_{0,\rm{LL-
GRB}}=0.5$–$0.8$ Gpc-3 yr-1, $z^{\rm{HL-GRB}}_{\star}=3.6$, and
$z^{\rm{sGRB}}_{\star}=6.0$. We have $\alpha=[0.049,\;0.047,\;0.048$],
$\beta=[0.636,\;0.609,\;0.586$], $\gamma=[1.9,\;1.9,\;2.0$],
$\delta=[0.01,\;0.49,\;0.98$], $\eta=[0.6,\;0.53,\;0.23$], and
$\mu=[3.0,\;3.6,\;3.6$]. For each of the parameters of $R_{\mathrm{sGRB}}$, we
indicate fit values corresponding to the lower, average, and higher redshift
rate.

Table 2: Jet parameters assumed for HL-GRBs [54, 94], LL-GRBs [95, 96, 91] and sGRBs [95, 97, 98]. All other jet parameters are modeled as specified in Table 1 according to the jet model.  | HL-GRB | LL-GRB | sGRB
---|---|---|---
$\tilde{E}_{\text{iso}}$ |  $4.5\times 10^{54}$ erg |  $1.5\times 10^{52}$ erg |  $1.05\times 10^{53}$ erg
$\Gamma$ | $300$ | $10$ | $300$
$t_{\text{dur}}$ |  $30$ s |  $1000$ s |  $1$ s
$t_{\text{var}}$ |  $0.5$ s |  $100$ s |  $0.01$ s
Figure 9: Diffuse $\nu_{\mu}+\bar{\nu}_{\mu}$ emission for HL-GRBs (ochre),
LL-GRBs (purple) and sGRBs (light blue) for the internal shock, photospheric
and ICMART models, from left to right, respectively, computed assuming the
average redshift rate. For each jet model and GRB type, the band marks the
uncertainty due to the jet composition ($p$-jets are plotted as solid lines
and Fe-jets as dotted lines; any mixed composition is represented by the color
band). For each GRB type, the color band is shadowed for the cases where
nuclear contamination is expected to be disfavored. The uncertainty due to the
jet composition is smaller or at most comparable to the uncertainty due to
redshift rate. The latter is represented by solid lines for the $p$-jet.
Independent of the GRB type, the diffuse neutrino flux is expected to be lower
for the internal shock model.

As for the jet parameters entering the neutrino fluence, we differentiate
among the three GRB types as summarized in Table 2 and otherwise use the
parameters in Table 1. For each GRB type, we compute the diffuse emission
considering $p$-jet and Fe-jet.

Figure 9 shows the expected diffuse emission for HL-GRBs, LL-GRBs and sGRBs,
for the three jet models, from left to right, respectively. For fixed jet
model, for each GRB type, the band represents the uncertainty due to the jet
composition (with the lower emission corresponding to $p$-jet and the larger
one to Fe-jet). In order to compare the uncertainty due to the jet composition
with the one due to the redshift rate, we plot three solid lines for each GRB
type, representing the diffuse flux of a $p$-jet assuming the lowest, average
and highest GRB rate. Independent of the GRB type, the diffuse neutrino
emission is expected to be smaller for the internal shock model.

In order to allow for a comparison across different jet models, we use the
same jet parameters (e.g. isotropic energy and total duration of the burst)
for the internal shock, photospheric and ICMART models. This implies that, due
to the smaller microphysical parameters typical of the internal shock model,
the overall gamma-ray luminosity predicted within the internal shock model is
smaller than the one of the other two models, and about two orders of
magnitude below the peak of the luminosity function for each GRB family [94,
99, 100, 98]. While this choice might seem conservative, we are interested in
exploring the variation of the neutrino signal across jet models and
composition relying on the same benchmark jet parameters.

While we refrain from an exploration of the detection prospects, that would
require a dedicated exploration of the GRB population and survival of nuclei
according to the jet properties, it is interesting to note that the
uncertainty on the GRB rate is generally larger or at most comparable to the
uncertainty due to the composition. Magnetically dominated jets and low-
luminosity ones are the most likely to guarantee the survival of heavy nuclei
up to the interaction region [31]. For this reason, we use darker (lighter)
shaded bands to distinguish between the GRB types that are likely (unlikely)
to be loaded with heavier nuclei at the interaction region.

##  8 Discussion and outlook

It is likely that the GRB jet is not only composed by composed of protons, but
also heavier nuclei, either synthetized or entrained in the jet [30, 31]. In
this paper, we have explored the impact of the jet composition on the neutrino
emission for kinetic dominated jets (according to the internal shock and
photospheric models) as well as for Poynting flux dominated jets (considering
the ICMART model). To this purpose, we have built a Monte-Carlo algorithm that
simulates nuclear cascades for any nucleus with $A\leq 56$.

We find that the nuclear composition affects the neutrino fluence in
quantitatively different ways according to the jet model, although the overall
qualitative trend is similar across models. The neutrino fluence can be lower
up to one order of magnitude for the Fe-jet case with respect to the $p$-jet
scenario because the production of mesons from nuclei is less efficient than
the one from protons. The efficiency of neutrino production is linked to the
jet model; the ICMART model, because of its fixed particle interaction radius,
guarantees efficient nuclear cascades leading to an overall larger neutrino
fluence across all three jet models, even for the Fe-jet case.

Independent of the jet model, the survival of nuclei and inefficient nuclear
cascades lead to a lower neutrino fluence due to the lower efficiency of
nuclear photohadronic interactions. Conversely, if nuclei are disintegrated,
the neutrino fluence may be similar to the one from a jet loaded with protons.
The sketch in Fig. 10 summarizes how the jet luminosity and the radius of
neutrino production, within a one-zone jet model, affect the nuclear cascade
efficiency and therefore the neutrino production (left and right panels,
respectively); see also Fig. 5.

Figure 10: Left: Sketch of the efficiency of nuclear cascades for our one-zone
jet model in the plane spanned by the radius of neutrino production and the
jet luminosity. The nuclear species dominating the particle density are
highlighted in color (protons and neutrons in red, any nucleus lighter than
iron in green, and iron in light blue). Right: Same as the left panel, except
for the fact that the different regions are colored depending on which
isotopes contribute the most to the total neutrino fluence. The efficiency of
nuclear cascades is the largest for small interaction radii and larger jet
luminosity. Notably, when nuclear cascades are efficient, the neutrino
production is driven by proton and neutrons in a larger region of the
parameter space (cf. green band that is smaller for the right panel).

Exploring the jet parameter space, we showed that the power law index of the
seed spectra of nuclei mostly affects the normalization of the spectral energy
distribution of neutrinos, without altering its shape. On the other hand, the
bulk Lorentz factor of the jet is important to assess the efficiency of
nuclear cascades and is thus highly correlated with both the composition and
the GRB model. Our work focuses on neutrinos only, complementing dedicated
work in this direction centered on ultra-high-energy cosmic rays, e.g. Refs.
[32, 34].

We have explored two extreme cases in terms of composition: either assuming
that the jet is fully made out of protons or pure 56Fe. Any other mixed
composition, involving nuclei lighter than 56Fe is expected to fall in between
the two extreme scenarios considered here. However, it is important to keep in
mind that, according to the jet luminosity and Lorentz boost factor, heavier
nuclei may not survive in the jet [31, 101]; we do not explore the conditions
for survival of nuclei in this paper, rather focus on assessing the
uncertainty linked to the jet composition across jet models. Further work
should be dedicate to explore the neutrino and cosmic ray emission from jets
with mixed composition, guided by self-consistent simulations of the source
and their ejecta composition. Nuclear cascades from GRBs could also contribute
to the spectrum of ultra-high-energy cosmic rays [101, 102, 37, 32].

Modeling the diffuse neutrino emission from HL-GRBs, LL-GRBs and sGRB, for our
three jet models and varying jet composition, we find that also for Fe-jets,
the diffuse neutrino flux from LL-GRBs may dominates the overall GRB diffuse
emission, consistently with the findings of Refs. [95, 103] for $p$-jets. Our
results on the diffuse neutrino emission are in overall good agreement with
Ref. [38], which predicted that the neutrino production from nuclei should be
up to one order of magnitude lower in the limit of inefficient nuclear
cascades. However, the uncertainty on the jet composition is generally smaller
or at most comparable with the one associated to the GRB cosmological rate.

## Acknowledgments

We are grateful to Annika Rudolph for involvement in the very early stages of
this project and Walter Winter for helpful comments on the manuscript. This
research project has received support from the Villum Foundation (Project No.
37358), the Deutsche Forschungsgemeinschaft through Sonderforschungsbereich
SFB 1258 “Neutrinos and Dark Matter in Astro- and Particle Physics” (NDM), and
the École Normale Supérieure Paris-Saclay scholarship program. V.D.L. also
thanks the Particle Astrophysics group members at the Niels Bohr Institute for
their warm hospitality.

##  Appendix A Details on the Monte Carlo simulation of nuclear cascades

In this appendix, we present details on the modeling of the photonuclear cross
section for any nuclei with mass number $A\leq 56$, as well as the half-life
and radioactive decay channels of unstable isotopes. We also describe the
method adopted to construct the particle spectral energy distributions and the
modeling of the nuclear photohadronic interactions.

###  A.1 Modeling of the nuclei photonuclear cross section

In order to model the photonuclear cross sections we distinguish among three
different energy regimes:

  * Giant Dipole Resonance 1. We consider the energy range between $0$ and $30$ MeV as being the Giant Dipole Resonance 1 (GDR1) one; it coincides with the GDR peak. We use the fit from Ref. [85] to compute the total inelastic cross-section for any nucleus of mass number $A\leq 56$ in this energy region:

$\displaystyle\sigma_{A\gamma}(E_{\gamma,r},A)=\sum_{i=1,2,4,8}\frac{e^{i\,(\rho_{i}-z)}}{1+e^{3i\,(\tau_{i}-z)}}+f_{r}\,(r_{H}+r_{\Delta})+s_{p}\,f_{p}\,h_{p}\
,$ (A.1)

with $z=\ln{(E_{\gamma,r}/1\,\text{MeV})}$, and the fit parameters $\rho_{i}$,
$\tau_{i}$, $f_{r}$, $r_{H}$, $r_{\Delta}$, $s_{p}$, $f_{p}$, and $h_{p}$
being defined as in Ref. [85]. Figure 11 compares this fit to the GDR model
presented in Ref. [86], displaying a satisfactory agreement in the GDR1 energy
range.

In the GDR1 energy range, the emission of only one or two nucleons is
possible. The fit of Ref. [85] does not provide this information, hence the
data from Ref. [86] are used to compute the branching ratio $P_{\mathcal{N}}$
for $\mathcal{N}=Z+N$ nucleons emission . The branching ratio $P$ is a
function of the number $(Z,N)$ of protons and neutrons emitted and is computed
as follows:

$P(Z,N)=P_{\mathcal{N}}(Z+N)\,\binom{Z+N}{Z}\,\frac{1}{2^{Z+N}}\ ,$ (A.2)

where the binomial distribution is equivalent to the hypothesis that protons
and neutrons have the same probability of being emitted and all emission
channels are independent.

Figure 11: Comparison among different models of the 56Fe photonuclear
inelastic cross section. The blue, orange and green curves refer to the models
presented in Refs. [86], [85] and [79], respectively. The dashed vertical
lines highlight the different energy ranges considered in this work (GDR1,
GDR2, and PM).

  * Giant Dipole Resonance 2. The second energy range is the Giant Dipole Resonance 2 (GDR2) one, which spans from $30$ to $150$ MeV. This range coincides with the transition between the GDR peak and the photomeson range. The fit in Eq. (A.1) is still used to compute the cross-section in this energy range and is displayed in Fig. 11. A notable improvement on the modeling of the cross section in this energy range has been achieved in Ref. [85], with respect to Ref. [86]. An estimation of the branching ratios $P_{\mathcal{N}}$ can be found in Ref. [86] (see Eq. (A.2)).

  * Photomeson regime. The photomeson energy range (PM) is defined for energies above $150$ MeV, where the production of mesons is possible thanks to interactions with photons. Ref. [79] provides the cross-section and the estimation of the branching ratio in the PM range. We have performed a polynomial fit of the total inelastic cross-section (in microbarn):

$\sigma_{A\gamma}^{PM}(E_{\gamma,r},A)=A\left(\frac{A}{56}\right)^{\alpha(E_{\gamma,r})-1}\begin{cases}\sigma_{\text{low}}(E_{\gamma,r})&\text{if
}0.18\;\text{GeV}\leq E_{\gamma,r}\leq 1\;\text{GeV}\\\
\sigma_{\text{mid}}(E_{\gamma,r})&\text{if }1\;\text{GeV}\leq E_{\gamma,r}\leq
10^{4}\;\text{GeV}\ .\\\ \sigma_{\text{high}}(E_{\gamma,r})&\text{if
}E_{\gamma,r}\geq 10^{4}\;\text{GeV}\end{cases}$ (A.3)

Using $z=\log_{10}(E_{\gamma,r}/1\text{GeV})$, the expressions for the cross
section above in the three energy ranges are respectively:

$\displaystyle\sigma_{\text{low}}({E_{\gamma,r}})$ $\displaystyle=$
$\displaystyle 2720.9z^{5}-1979.5z^{4}-1139.4z^{3}+902.4z^{2}-165.1z+139.1$
(A.4) $\displaystyle\sigma_{\text{mid}}({E_{\gamma,r}})$ $\displaystyle=$
$\displaystyle-2.02z^{5}+25.32z^{4}-112.23z^{3}+212.91z^{2}-188.22z+156.26\ ;$
(A.5) $\displaystyle\sigma_{\text{high}}({E_{\gamma,r}})$ $\displaystyle=$
$\displaystyle 39.1$ (A.6)

while the function $\alpha(E_{\gamma,r})$ is defined as follows

$\alpha(E_{\gamma,r})=\begin{cases}1&\text{if }E_{\gamma,r}\leq
0.2\;\text{GeV}\\\
0.004z^{5}-0.023z^{4}+0.019z^{3}+0.018z^{2}-0.043z+0.97&\text{if
}0.2\;\text{GeV}\leq E_{\gamma,r}\leq 10^{3}\;\text{GeV}\\\ 2/3&\text{if
}E_{\gamma,r}\geq 10^{3}\;\text{GeV}\end{cases}$ (A.7)

This PM parametrization of the cross section is displayed in Fig. 11.

The cross-sections of each process averaged over the PM energy range are used
to compute the branching ratios for the five types of photonuclear reactions:
direct proton production, direct neutron production, multi neutron production,
spallation and pion production. The branching ratio of a process is estimated
as the ratio of its cross-section over the total cross-section.

###  A.2 Radioactive decay of unstable isotopes

For each isotope, we extract the mean lifetime (or half-life) and decay
channel from the International Atomic Energy Agency database [87]. Figure 12
shows the half-life (left panel) and decay type (right panel) for isotopes of
mass number $A\leq 56$. We can see that the lifetimes can vary by orders of
magnitudes. Moreover, the isotope stability decreases exponentially as the
isotope position in the plane spanned by $N$ and $Z$ strays from the valley of
stability (which corresponds to $Z=N=A/2$). An unstable isotope can undergo
five types of radioactive decays: $\beta^{+}$ decay, $\beta^{-}$ decay,
$\alpha$ decay, single proton emission and single neutron emission. Single
proton and neutron emission are unusual decay types that describe the behavior
of the most proton-rich/neutron-rich isotopes. Such isotopes are so unstable
that they emit a proton/neutron almost instantaneously. No precise measurement
of their lifetime is available; this explains why they are not highlighted in
color in the left panel of Fig. 12. Consequently, if a nucleus arrives in the
proton/neutron emission region, the emission of single nucleons can be
integrated out until the isotope reaches a $\beta^{+}$, $\beta^{-}$, $\alpha$
or stable state (as indicated in the contour scale of the right panel).

Figure 12: Contour plots of the half-life (left panel) and radioactive decay
channels (right panel) for nuclei of mass number $A\leq 56$ [87] in the plane
spanned by $N$ and $Z$. The lifetime of different isotopes can largely vary.
The isotope stability decreases exponentially as the isotope strays from the
valley of stability ($Z=N=A/2$; corresponding to the black region on the right
panel); in this case, the isotope can undergo five different types of
radioactive decays: $\beta^{+}$ decay, $\beta^{-}$ decay, $\alpha$ decay,
single proton emission and single neutron emission.

###  A.3 Spectral energy distributions

In order to generate the spectral energy distribution for each nucleon, we
first select the seed isotope energy randomly out of a logarithmic
distribution:

$\log_{10}(E^{\prime}_{0})=\log_{10}(E^{\prime}_{\text{min}})+(\log_{10}(E^{\prime}_{\text{max}})-\log_{10}(E^{\prime}_{\text{min}}))\;U(0,1)\
,$ (A.8)

with $U(0,1)$ being the uniform probability density function between $0$ and
$1$, $E^{\prime}_{\text{min}}$ and $E^{\prime}_{\text{max}}$ the minimum and
maximum energies.

The spectral density $n^{\prime}_{A}(E^{\prime}_{A})$ is a set of $N_{s}$
points located at energies $E^{\prime}_{s,i}$ with $0\leq i\leq N_{s}-1$.
Hence, to each nucleus with initial energy $E^{\prime}_{0}$, we can associate
a number density:

$d_{0}(E^{\prime}_{0})=\int_{E^{\prime}_{s,i_{0}}}^{E^{\prime}_{s,i_{0}+1}}n^{\prime}_{A}(E^{\prime})dE^{\prime}\approx
n^{\prime}_{A}(E^{\prime}_{s,i_{0}})\;(E^{\prime}_{s,i_{0}+1}-E^{\prime}_{s,i_{0}})\
,$ (A.9)

with $i_{0}$ such that $E^{\prime}_{s,i_{0}}\leq E^{\prime}_{0}\leq
E^{\prime}_{s,i_{0}+1}$. If a test nucleus undergoes a photonuclear reaction
that generates $(N_{p},N_{n})$ free protons and neutrons, this creates the
nucleon density $d_{\mathcal{N}}=N_{\mathcal{N}}\;d_{0}(E^{\prime}_{0})$, with
$\mathcal{N}=(p,n)$ representing a proton or a neutron. Note that these
particle densities are independent of the energy at which the test nucleus
interacts with a photon and the energy of the produced nucleons.

We then add these nucleons together to populate the spectral energy
distribution for each species. Nucleons have an energy $E^{\prime}_{f}$ at the
end of the Monte-Carlo run. Their density is computed as

$d_{\mathcal{N}}(E^{\prime}_{f})=\int_{E^{\prime}_{s,i_{f}}}^{E^{\prime}_{s,i_{f}+1}}n^{\prime}_{\mathcal{N}}(E^{\prime})dE^{\prime}\approx
n^{\prime}_{\mathcal{N}}(E^{\prime}_{s,i_{f}})\;(E^{\prime}_{s,i_{f}+1}-E^{\prime}_{s,i_{f}})\
,$ (A.10)

with $i_{f}$ such that $E^{\prime}_{s,i_{f}}\leq
E^{\prime}_{f}\leq^{\prime}E_{s,i_{f}+1}$. Using
$d_{\mathcal{N}}=N_{\mathcal{N}}\;d(E^{\prime}_{0})$ and Eqs. (A.9) and
(A.10), we get

$n^{\prime}_{\mathcal{N}}(E^{\prime}_{s,i_{f}})=N_{\mathcal{N}}\;n^{\prime}_{A}(E^{\prime}_{s,i_{0}})\;\frac{E^{\prime}_{s,i_{0}+1}-E^{\prime}_{s,i_{0}}}{E^{\prime}_{s,i_{f}+1}-E^{\prime}_{s,f_{0}}}\
.$ (A.11)

We note that these calculations presented in this section can be easily
generalized to compute the spectrum of the associated photodisintegrated
nuclei.

##  Appendix B Nuclear photohadronic model

In order to estimate the meson production, the function $f$ (cf. Eq. (5.4))
should be identified for each energy interval $\Delta
E=[y_{\min},\;y_{\max}]$:

$f^{\Delta E}(y)=\frac{1}{2y^{2}}\;\int_{E^{\mathrm{th},\Delta
E}_{\gamma,r}}^{2y}\;E_{\gamma,r}\;\sigma_{A\gamma\rightarrow\mathcal{M}}(E_{\gamma,r})\;dE_{\gamma,r}$
(B.1)

To this purpose, we compute the total function ($f^{\mathrm{tot}}$) using the
total cross-section fit from Eq. (A.3).

The function $f^{\mathrm{tot}}$ does not only describe the meson production,
but also all the other processes in this energy range. That is why it needs to
be scaled by the branching ratios of meson production. Using data for the
$\pi^{0}$ production cross-section [79], we compute an approximation of the
product of the multiplicity and the branching ratio and find
$M_{\pi^{0}}\;{\sigma_{A\gamma\rightarrow\pi^{0}}^{\mathrm{PM}}}/{\sigma_{A\gamma}^{\mathrm{PM}}}\approx
0.18$. Then, using data for $(p,n)$–$\gamma$ interactions from Ref. [78] and
assuming that they are still valid for nuclei, we compute the ratios
$f_{\pi^{+}}/f_{\pi^{0}}$ and $f_{\pi^{-}}/f_{\pi^{0}}$ as approximations of
the multiplicity ratio $M_{\pi^{+}}/M_{\pi^{0}}$ and
$M_{\pi^{-}}/M_{\pi^{0}}$.

Figure 13:  Ratio of the function $f_{\pi^{+}}/f_{\pi^{0}}$ and multiplicity
ratio $M_{\pi^{+}}/M_{\pi^{0}}$ in the nuclear photohadronic interaction model
as functions of $y$.

This is shown in Fig. 13. We highlight that
$M_{\pi^{+}}/M_{\pi^{0}}=M_{\pi^{-}}/M_{\pi^{0}}$ because $50\%$ of
interactions are expected to happen with a proton and the rest with a neutron
of the nucleus. Indeed, the parameters describing $p\gamma$ and $n\gamma$
photohadronic interactions are the same, except that produced mesons have
opposite charge [78]. As a consequence, we have symmetry between the
production of $\pi^{+}$ and the one of $\pi^{-}$. Following Ref. [78], we also
compute approximations of the fraction of energy transferred to the produced
meson $\chi_{\mathcal{M}}$. Both $M_{\pi}$ and $\chi_{\mathcal{M}}$ need to be
constant as functions of the energy for a given energy range. Thus, we
consider seven energy intervals $\Delta E$, each associated to a range where
the multiplicity and the energy fraction is approximately constant, as
summarized in Table 3.

Table 3: Multiplicities and fraction of the energy transferred to pions for the seven considered energy ranges. Intervals $\Delta E$ | $y_{\min}=E^{\mathrm{th}}_{\gamma,r}/2$ | $y_{\max}=E^{\max}_{\gamma,r}/2$ | $M_{\pi^{+}}/M_{\pi^{0}}=M_{\pi^{-}}/M_{\pi^{0}}$ | $\chi_{\pi}\times A$
---|---|---|---|---
1 | 0.1 GeV | 0.25 GeV | 0.25 | 0.22
2 | 0.25 GeV | 1.75 GeV | 0.6 | 0.2
3 | 1.75 GeV | 2.5 GeV | 0.86 | 0.25
4 | 2.5 GeV | 25 GeV | 0.95 | 0.18
5 | 25 GeV | 250 GeV | 0.91 | 0.25
6 | 250 GeV | 2500 GeV | 0.89 | 0.25
7 | 2500 GeV | $+\infty$ | 0.88 | 0.3

The product of the multiplicity and the function $f^{\Delta E}$, taking into
account both branching ratio and multiplicity, is:

$M_{\pi^{+}}\,f^{\Delta
E}(y)\approx\left(\frac{M_{\pi^{+}}}{M_{\pi^{0}}}\right)^{\Delta
E}\left(M_{\pi^{0}}\;\frac{\sigma_{A\gamma\rightarrow\pi^{0}}^{\mathrm{PM}}}{\sigma_{A\gamma}^{\mathrm{PM}}}\right)\left(f^{\mathrm{tot}}(y)-f^{\mathrm{tot}}\left(y^{\Delta
E}_{\mathrm{th}}\right)\left(\frac{y^{\Delta
E}_{\mathrm{th}}}{y}\right)^{2}\right).$ (B.2)

On the right hand side of the equation above, the second term is subtracted
from $f^{\mathrm{tot}}(y)$ since the latter is an integral with lower limit
equal to $E^{\mathrm{th},\mathrm{PM}}_{\gamma,r}=0.2$ GeV (and not from
$E^{\mathrm{th},\Delta E}_{\gamma,r}$):

$\displaystyle f^{\mathrm{tot}}(y)$ $\displaystyle=$
$\displaystyle\frac{1}{2y^{2}}\;\int_{E^{\mathrm{th},\mathrm{PM}}_{\gamma,r}}^{2y}\;E_{\gamma,r}\;\sigma_{A\gamma}(E_{\gamma,r})\;dE_{\gamma,r}=f^{\mathrm{tot}}\left(y^{\Delta
E}_{\mathrm{th}}\right)\left(\frac{y^{\Delta E}_{\mathrm{th}}}{y}\right)^{2}+$
$\displaystyle+$ $\displaystyle\frac{1}{2y^{2}}\;\int_{E^{\mathrm{th},\Delta
E}_{\gamma,r}}^{2y}\;E_{\gamma,r}\;\sigma_{A\gamma}(E_{\gamma,r})\;dE_{\gamma,r}\
.$

Finally, knowing the parameters $M_{\pi}$, $\chi_{\pi}$ and $f^{\Delta E}$, we
use Eq. (5.4) to compute the pion production.

## References

  * [1] D. Frederiks et al., _Properties of the Extremely Energetic GRB 221009A from Konus-WIND and SRG/ART-XC Observations_ , _Astrophys. J. Lett._ 949 (2023) L7, [2302.13383].
  * [2] D. Eichler, M. Livio, T. Piran and D. N. Schramm, _Nucleosynthesis, Neutrino Bursts and Gamma-Rays from Coalescing Neutron Stars_ , _Nature_ 340 (1989) 126–128.
  * [3] B. Paczynski, _Cosmological gamma-ray bursts_ , _Acta Astron._ 41 (1991) 257–267.
  * [4] E. Nakar, _Short-Hard Gamma-Ray Bursts_ , _Phys. Rept._ 442 (2007) 166–236, [astro-ph/0701748].
  * [5] LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL collaboration, B. P. Abbott et al., _Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A_ , _Astrophys. J. Lett._ 848 (2017) L13, [1710.05834].
  * [6] LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT collaboration, B. P. Abbott et al., _Multi-messenger Observations of a Binary Neutron Star Merger_ , _Astrophys. J. Lett._ 848 (2017) L12, [1710.05833].
  * [7] K. P. Mooley et al., _A mildly relativistic wide-angle outflow in the neutron star merger GW170817_ , _Nature_ 554 (2018) 207, [1711.11573].
  * [8] A. Goldstein et al., _An Ordinary Short Gamma-Ray Burst with Extraordinary Implications: Fermi-GBM Detection of GRB 170817A_ , _Astrophys. J. Lett._ 848 (2017) L14, [1710.05446].
  * [9] S. E. Woosley and J. S. Bloom, _The Supernova Gamma-Ray Burst Connection_ , _Ann. Rev. Astron. Astrophys._ 44 (2006) 507–556, [astro-ph/0609142].
  * [10] M. Modjaz, Y. Q. Liu, F. B. Bianco and O. Graur, _The Spectral SN-GRB Connection: Systematic Spectral Comparisons between Type Ic Supernovae, and broad-lined Type Ic Supernovae with and without Gamma-Ray Bursts_ , _Astrophys. J._ 832 (2016) 108, [1509.07124].
  * [11] Y.-H. Yang et al., _A lanthanide-rich kilonova in the aftermath of a long gamma-ray burst_ , _Nature_ 626 (Feb., 2024) 742–745, [2308.00638].
  * [12] E. Troja et al., _A nearby long gamma-ray burst from a merger of compact objects_ , _Nature_ 612 (2022) 228–231, [2209.03363].
  * [13] J. H. Gillanders et al., _Heavy element nucleosynthesis associated with a gamma-ray burst_ , 2308.00633.
  * [14] H. Sun, B. Zhang and Z. Li, _Extragalactic High-energy Transients: Event Rate Densities and Luminosity Functions_ , _Astrophys. J._ 812 (2015) 33, [1509.01592].
  * [15] F. Virgili, E. Liang and B. Zhang, _Low-Luminosity Gamma-Ray Bursts as a Distinct GRB Population:A Monte Carlo Analysis_ , _Mon. Not. Roy. Astron. Soc._ 392 (2009) 91, [0801.4751].
  * [16] B. Paczynski, _Gamma-ray bursters at cosmological distances_ , _Astrophys. J. Lett._ 308 (1986) L43–L46.
  * [17] G. Drenkhahn and H. C. Spruit, _Efficient acceleration and radiation in Poynting flux powered GRB outflows_ , _Astron. Astrophys._ 391 (2002) 1141, [astro-ph/0202387].
  * [18] M. J. Rees and P. Mészáros, _Unsteady outflow models for cosmological gamma-ray bursts_ , _Astrophys. J. Lett._ 430 (1994) L93–L96, [astro-ph/9404038].
  * [19] A. Spitkovsky, _Particle acceleration in relativistic collisionless shocks: Fermi process at last?_ , _Astrophys. J. Lett._ 682 (2008) L5, [0802.3216].
  * [20] A. M. Beloborodov and P. Mészáros, _Photospheric Emission of Gamma-Ray Bursts_ , _Space Sci. Rev._ 207 (2017) 87–110, [1701.04523].
  * [21] A. Pe’er and F. Ryde, _Photospheric Emission in Gamma-Ray Bursts_ , _Int. J. Mod. Phys. D_ 26 (2017) 1730018, [1603.05058].
  * [22] H. C. Spruit, F. Daigne and G. Drenkhahn, _Large scale magnetic fields and their dissipation in grb fireballs_ , _Astron. Astrophys._ 369 (2001) 694, [astro-ph/0004274].
  * [23] A. Rudolph, I. Tamborra and O. Gottlieb, _Subphotospheric Emission from Short Gamma-Ray Bursts: Protons Mold the Multimessenger Signals_ , _Astrophys. J. Lett._ 961 (2024) L7, [2309.08667].
  * [24] E. Waxman, _Gamma-ray bursts: The Underlying model_ , _Lect. Notes Phys._ 598 (2003) 393, [astro-ph/0303517].
  * [25] P. Mészáros, _Astrophysical Sources of High Energy Neutrinos in the IceCube Era_ , _Ann. Rev. Nucl. Part. Sci._ 67 (2017) 45–67, [1708.03577].
  * [26] T. Piran, _The physics of gamma-ray bursts_ , _Rev. Mod. Phys._ 76 (2004) 1143–1210, [astro-ph/0405503].
  * [27] E. Moore, B. Gendre, N. B. Orange and F. H. Panther, _Constraints on the ultra-high energy cosmic ray output of gamma-ray bursts_ , _Mon. Not. Roy. Astron. Soc._ 530 (2024) 555–559, [2303.13781].
  * [28] E. Guarini, I. Tamborra, R. Margutti and E. Ramirez-Ruiz, _Transients stemming from collapsing massive stars: The missing pieces to advance joint observations of photons and high-energy neutrinos_ , _Phys. Rev. D_ 108 (2023) 083035, [2308.03840].
  * [29] D. Boncioli, D. Biehl and W. Winter, _On the common origin of cosmic rays across the ankle and diffuse neutrinos at the highest energies from low-luminosity Gamma-Ray Bursts_ , _Astrophys. J._ 872 (2019) 110, [1808.07481].
  * [30] A. M. Beloborodov, _Nuclear composition of gamma-ray burst fireballs_ , _Astrophys. J._ 588 (2003) 931–944, [astro-ph/0210522].
  * [31] S. Horiuchi, K. Murase, K. Ioka and P. Mészáros, _The survival of nuclei in jets associated with core-collapse supernovae and gamma-ray bursts_ , _Astrophys. J._ 753 (2012) 69, [1203.0296].
  * [32] D. Biehl, D. Boncioli, A. Fedynitch and W. Winter, _Cosmic-Ray and Neutrino Emission from Gamma-Ray Bursts with a Nuclear Cascade_ , _Astron. Astrophys._ 611 (2018) A101, [1705.08909].
  * [33] T. Pitik, I. Tamborra and M. Petropoulou, _Neutrino signal dependence on gamma-ray burst emission mechanism_ , _JCAP_ 05 (2021) 034, [2102.02223].
  * [34] J. Heinze, D. Biehl, A. Fedynitch, D. Boncioli, A. Rudolph and W. Winter, _Systematic parameter space study for the UHECR origin from GRBs in models with multiple internal shocks_ , _Mon. Not. Roy. Astron. Soc._ 498 (2020) 5990–6004, [2006.14301].
  * [35] A. Rudolph, M. Petropoulou, Ž. Bošnjak and W. Winter, _Multicollision Internal Shock Lepto-hadronic Models for Energetic Gamma-Ray Bursts (GRBs)_ , _Astrophys. J._ 950 (2023) 28, [2212.00765].
  * [36] K. Murase, _Prompt High-Energy Neutrinos from Gamma-Ray Bursts in the Photospheric and Synchrotron Self-Compton Scenarios_ , _Phys. Rev. D_ 78 (2008) 101302, [0807.0919].
  * [37] X.-Y. Wang, S. Razzaque and P. Mészáros, _On the Origin and Survival of UHE Cosmic-Ray Nuclei in GRBs and Hypernovae_ , _Astrophys. J._ 677 (2008) 432–440, [0711.2065].
  * [38] K. Murase and J. F. Beacom, _Neutrino Background Flux from Sources of Ultrahigh-Energy Cosmic-Ray Nuclei_ , _Phys. Rev. D_ 81 (2010) 123001, [1003.4959].
  * [39] K. Toma, X.-F. Wu and P. Mészáros, _A Photosphere-Internal Shock Model of Gamma-Ray Bursts: Case Studies of Fermi/LAT Bursts_ , _Mon. Not. Roy. Astron. Soc._ 415 (2011) 1663–1680, [1002.2634].
  * [40] B. Zhang and H. Yan, _The Internal-Collision-Induced Magnetic Reconnection and Turbulence (ICMART) Model of Gamma-Ray Bursts_ , _Astrophys. J._ 726 (2011) 90, [1011.1197].
  * [41] O. Bromberg, E. Nakar, T. Piran and R. Sari, _The propagation of relativistic jets in external media_ , _Astrophys. J._ 740 (2011) 100, [1107.1326].
  * [42] S. Kobayashi, T. Piran and R. Sari, _Can internal shocks produce the variability in GRBs?_ , _Astrophys. J._ 490 (1997) 92–98, [astro-ph/9705013].
  * [43] F. Daigne and R. Mochkovitch, _Gamma-ray bursts from internal shocks in a relativistic wind: temporal and spectral properties_ , _Mon. Not. Roy. Astron. Soc._ 296 (1998) 275, [astro-ph/9801245].
  * [44] D. Guetta, M. Spada and E. Waxman, _Efficiency and spectrum of internal gamma-ray burst shocks_ , _Astrophys. J._ 557 (2001) 399, [astro-ph/0011170].
  * [45] M. Bustamante, K. Murase, W. Winter and J. Heinze, _Multi-messenger light curves from gamma-ray bursts in the internal shock model_ , _Astrophys. J._ 837 (2017) 33, [1606.02325].
  * [46] A. Rudolph, J. Heinze, A. Fedynitch and W. Winter, _Impact of the Collision Model on the Multi-messenger Emission from Gamma-Ray Burst Internal Shocks_ , _Astrophys. J._ 893 (2020) 72, [1907.10633].
  * [47] N. Globus, D. Allard, R. Mochkovitch and E. Parizot, _UHECR acceleration at GRB internal shocks_ , _Mon. Not. Roy. Astron. Soc._ 451 (2015) 751–790, [1409.1271].
  * [48] L. Sironi and A. Spitkovsky, _Particle Acceleration in Relativistic Magnetized Collisionless Electron-Ion Shocks_ , _Astrophys. J._ 726 (2011) 75, [1009.0024].
  * [49] P. Crumley, D. Caprioli, S. Markoff and A. Spitkovsky, _Kinetic simulations of mildly relativistic shocks – I. Particle acceleration in high Mach number shocks_ , _Mon. Not. Roy. Astron. Soc._ 485 (2019) 5105–5119, [1809.10809].
  * [50] P. Lipari, M. Lusignoli and D. Meloni, _Flavor Composition and Energy Spectrum of Astrophysical Neutrinos_ , _Phys. Rev. D_ 75 (2007) 123005, [0704.0718].
  * [51] X.-G. Wang et al., _How bad or Good are the External Forward Shock Afterglow Models of Gamma-ray Bursts?_ , _Astrophys. J. Suppl._ 219 (2015) 9, [1503.03193].
  * [52] IceCube collaboration, M. G. Aartsen et al., _Extending the search for muon neutrinos coincident with gamma-ray bursts in IceCube data_ , _Astrophys. J._ 843 (2017) 112, [1702.06868].
  * [53] G. Ghirlanda, F. Nappo, G. Ghisellini, A. Melandri, G. Marcarini, L. Nava et al., _Bulk Lorentz factors of Gamma-Ray Bursts_ , _Astron. Astrophys._ 609 (2018) A112, [1711.06257].
  * [54] H. Zitouni, N. Guessoum, K. M. ALQassimi and O. Alaryani, _Distributions of Pseudo-Redshifts and Durations (Observed and Intrinsic) of Fermi GRBs_ , _Astrophys. Space Sci._ 363 (2018) 223, [1810.04124].
  * [55] ANTARES collaboration, A. Albert et al., _Constraining the contribution of Gamma-Ray Bursts to the high-energy diffuse neutrino flux with 10 yr of ANTARES data_ , _Mon. Not. Roy. Astron. Soc._ 500 (2020) 5614–5628, [2008.02127].
  * [56] W. Deng, H. Li, B. Zhang and S. Li, _Relativistic MHD simulations of collision-induced magnetic dissipation in Poynting-flux-dominated jets/outflows_ , _Astrophys. J._ 805 (2015) 163, [1501.07595].
  * [57] L. Sironi, M. Petropoulou and D. Giannios, _Relativistic Jets Shine through Shocks or Magnetic Reconnection?_ , _Mon. Not. Roy. Astron. Soc._ 450 (2015) 183–191, [1502.01021].
  * [58] D. Groselj, L. Sironi and A. Spitkovsky, _Long-term Evolution of Relativistic Unmagnetized Collisionless Shocks_ , _Astrophys. J. Lett._ 963 (2024) L44, [2401.02392].
  * [59] L. Sironi, A. Spitkovsky and J. Arons, _The Maximum Energy of Accelerated Particles in Relativistic Collisionless Shocks_ , _Astrophys. J._ 771 (2013) 54, [1301.5333].
  * [60] H. Zhang, L. Sironi, D. Giannios and M. Petropoulou, _The Origin of Power-law Spectra in Relativistic Magnetic Reconnection_ , _Astrophys. J. Lett._ 956 (2023) L36, [2302.12269].
  * [61] L. Sironi and A. Spitkovsky, _Relativistic Reconnection: an Efficient Source of Non-Thermal Particles_ , _Astrophys. J. Lett._ 783 (2014) L21, [1401.5471].
  * [62] D. Band et al., _BATSE observations of gamma-ray burst spectra. 1. Spectral diversity._ , _Astrophys. J._ 413 (1993) 281–292.
  * [63] D. Gruber et al., _The Fermi GBM Gamma-Ray Burst Spectral Catalog: Four Years Of Data_ , _Astrophys. J. Suppl._ 211 (2014) 12, [1401.5069].
  * [64] L. Amati, _The E(p,i) - E(iso) correlation in grbs: updated observational status, re-analysis and main implications_ , _Mon. Not. Roy. Astron. Soc._ 372 (2006) 233–245, [astro-ph/0601553].
  * [65] J. N. Bahcall and P. Mészáros, _5-GeV to 10-GeV neutrinos from gamma-ray burst fireballs_ , _Phys. Rev. Lett._ 85 (2000) 1362–1365, [hep-ph/0004019].
  * [66] K. Murase, K. Kashiyama and P. Mészáros, _Subphotospheric Neutrinos from Gamma-Ray Bursts: The Role of Neutrons_ , _Phys. Rev. Lett._ 111 (2013) 131102, [1301.4236].
  * [67] K. Kashiyama, K. Murase and P. Mészáros, _Neutron-Proton-Converter Acceleration Mechanism at Subphotospheres of Relativistic Outflows_ , _Phys. Rev. Lett._ 111 (2013) 131103, [1304.1945].
  * [68] X.-Y. Wang and Z.-G. Dai, _Prompt TeV neutrinos from dissipative photospheres of gamma-ray bursts_ , _Astrophys. J. Lett._ 691 (2009) L67–L71, [0807.0290].
  * [69] E. Guarini, I. Tamborra and O. Gottlieb, _State-of-the-art collapsar jet simulations imply undetectable subphotospheric neutrinos_ , _Phys. Rev. D_ 107 (2023) 023001, [2210.03757].
  * [70] B. Zhang, _The Physics of Gamma-Ray Bursts_.  Cambridge University Press, 2018, 10.1017/9781139226530.
  * [71] G. R. Werner, D. A. Uzdensky, M. C. Begelman, B. Cerutti and K. Nalewajko, _Non-thermal particle acceleration in collisionless relativistic electron–proton reconnection_ , _Mon. Not. Roy. Astron. Soc._ 473 (2018) 4840–4861, [1612.04493].
  * [72] F. Guo, X. Li, H. Li, W. Daughton, B. Zhang, N. Lloyd-Ronning et al., _Efficient Production of High-energy Nonthermal Particles During Magnetic Reconnection in a Magnetically Dominated Ion–electron Plasma_ , _Astrophys. J. Lett._ 818 (2016) L9, [1511.01434].
  * [73] S. Gao, K. Asano and P. Mészáros, _High Energy Neutrinos from Dissipative Photospheric Models of Gamma Ray Bursts_ , _JCAP_ 11 (2012) 058, [1210.1186].
  * [74] C. D. Dermer and G. Menon, _High Energy Radiation from Black Holes: Gamma Rays, Cosmic Rays, and Neutrinos_.  Princeton University Press, 2009.
  * [75] F. C. Jones, _Inverse compton scattering of cosmic-ray electrons_ , _Phys. Rev._ 137 (Mar, 1965) B1306–B1311.
  * [76] E. Kafexhiu, F. Aharonian, A. M. Taylor and G. S. Vila, _Parametrization of gamma-ray production cross-sections for pp interactions in a broad proton energy range from the kinematic threshold to PeV energies_ , _Phys. Rev. D_ 90 (2014) 123014, [1406.7369].
  * [77] A. M. Lebedev, S. A. Slavatinskii and B. V. Tolkachev, _Interaction cross section and energy conserved by a nucleon colliding with complex nuclei_ , _JETP_ 46 (Jun, 1964) 2151.
  * [78] S. Hummer, M. Ruger, F. Spanier and W. Winter, _Simplified models for photohadronic interactions in cosmic accelerators_ , _Astrophys. J._ 721 (2010) 630–652, [1002.1310].
  * [79] L. Morejon, A. Fedynitch, D. Boncioli, D. Biehl and W. Winter, _Improved photomeson model for interactions of cosmic ray nuclei_ , _JCAP_ 11 (2019) 007, [1904.07999].
  * [80] L. A. Anchordoqui et al., _Cosmic Neutrino Pevatrons: A Brand New Pathway to Astronomy, Astrophysics, and Particle Physics_ , _JHEAp_ 1-2 (2014) 1–30, [1312.6587].
  * [81] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, _The fate of hints: updated global analysis of three-flavor neutrino oscillations_ , _JHEP_ 09 (2020) 178, [2007.14792].
  * [82] Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6, [1807.06209].
  * [83] F. W. Stecker and M. H. Salamon, _Photodisintegration of ultrahigh-energy cosmic rays: A New determination_ , _Astrophys. J._ 512 (1999) 521–526, [astro-ph/9808110].
  * [84] GEANT4 collaboration, S. Agostinelli et al., _GEANT4–a simulation toolkit_ , _Nucl. Instrum. Meth. A_ 506 (2003) 250–303.
  * [85] M. V. Kossov, _Approximation of photonuclear interaction cross-sections_ , _European Physical Journal A_ 14 (Jan., 2002) 377–392.
  * [86] J. L. Puget, F. W. Stecker and J. H. Bredekamp, _Photonuclear Interactions of Ultrahigh-Energy Cosmic Rays and their Astrophysical Consequences_ , _Astrophys. J._ 205 (1976) 638–654.
  * [87] IAEA, _Table of nuclides - nuclear strcuture and decay data_ , 2024.
  * [88] A. Achterberg, Y. A. Gallant, J. G. Kirk and A. W. Guthmann, _Particle acceleration by ultrarelativistic shocks: Theory and simulations_ , _Mon. Not. Roy. Astron. Soc._ 328 (2001) 393, [astro-ph/0107530].
  * [89] J.-J. Wei, X.-F. Wu, F. Melia, D.-M. Wei and L.-L. Feng, _Cosmological tests using gamma-ray bursts, the star formation rate and possible abundance evolution_ , _Mon. Not. Roy. Astron. Soc._ 439 (2014) 3329–3341, [1306.4415].
  * [90] D. Wanderman and T. Piran, _The luminosity function and the rate of Swift’s Gamma Ray Bursts_ , _Mon. Not. Roy. Astron. Soc._ 406 (2010) 1944–1958, [0912.0709].
  * [91] R.-Y. Liu, X.-Y. Wang and Z.-G. Dai, _Nearby low-luminosity GRBs as the sources of ultra-high energy cosmic rays revisited_ , _Mon. Not. Roy. Astron. Soc._ 418 (2011) 1382, [1108.1551].
  * [92] H. Yuksel, M. D. Kistler, J. F. Beacom and A. M. Hopkins, _Revealing the High-Redshift Star Formation Rate with Gamma-Ray Bursts_ , _Astrophys. J. Lett._ 683 (2008) L5–L8, [0804.4008].
  * [93] O. S. Salafia, M. E. Ravasio, G. Ghirlanda and I. Mandel, _The short gamma-ray burst population in a quasi-universal jet scenario_ , _Astron. Astrophys._ 680 (2023) A45, [2306.15488].
  * [94] G.-X. Lan, J.-J. Wei, H.-D. Zeng, Y. Li and X.-F. Wu, _Revisiting the luminosity and redshift distributions of long gamma-ray bursts_ , _Mon. Not. Roy. Astron. Soc._ 508 (2021) 52–68, [2109.00766].
  * [95] I. Tamborra and S. Ando, _Diffuse emission of high-energy neutrinos from gamma-ray burst fireballs_ , _JCAP_ 09 (2015) 036, [1504.00107].
  * [96] A. Rudolph, Ž. Bošnjak, A. Palladino, I. Sadeh and W. Winter, _Multiwavelength radiation models for low-luminosity GRBs and the implications for UHECRs_ , _Mon. Not. Roy. Astron. Soc._ 511 (2022) 5823–5842, [2107.04612].
  * [97] H. Ito, O. Just, Y. Takei and S. Nagataki, _A Global Numerical Model of the Prompt Emission in Short Gamma-ray Bursts_ , _Astrophys. J._ 918 (2021) 59, [2105.09323].
  * [98] D. Wanderman and T. Piran, _The rate, luminosity function and time delay of non-Collapsar short GRBs_ , _Mon. Not. Roy. Astron. Soc._ 448 (2015) 3026–3037, [1405.5878].
  * [99] E. J. Howell, D. M. Coward, G. Stratta, B. Gendre and H. Zhou, _Constraining the rate and luminosity function of Swift gamma-ray bursts_ , _Mon. Not. Roy. Astron. Soc._ 444 (2014) 15–28, [1407.2333].
  * [100] X. Dai, _Intensity Distribution and Luminosity Function of the Swift Gamma-Ray Bursts_ , _Astrophys. J. Lett._ 697 (2009) L68–L71, [0812.4466].
  * [101] B. T. Zhang, K. Murase, S. S. Kimura, S. Horiuchi and P. Mészáros, _Low-luminosity gamma-ray bursts as the sources of ultrahigh-energy cosmic ray nuclei_ , _Phys. Rev. D_ 97 (2018) 083010, [1712.09984].
  * [102] K. Murase, K. Ioka, S. Nagataki and T. Nakamura, _High-energy cosmic-ray nuclei from high- and low-luminosity gamma-ray bursts and implications for multimessenger astronomy_ , _Physical Review D_ 78 (July, 2008) .
  * [103] K. Murase, K. Ioka, S. Nagataki and T. Nakamura, _High Energy Neutrinos and Cosmic-Rays from Low-Luminosity Gamma-Ray Bursts?_ , _Astrophys. J. Lett._ 651 (2006) L5–L8, [astro-ph/0607104].