# Estimating the Neutrino Flux from Choked Gamma-Ray Bursts

Michela Fasano  Silvia Celli  Dafne Guetta  Antonio Capone  Angela Zegarelli
and Irene Di Palma

###### Abstract

The strong constraints from the Fermi-LAT data on the isotropic gamma-ray
background suggest that the neutrinos observed by IceCube might possibly come
from sources that are hidden to gamma-ray observations. A possibility recently
discussed in the literature is that neutrinos may come from jets of collapsing
massive stars which fail to break out of the stellar envelope, and for this
reason they are known as choked jets, or choked Gamma-Ray Bursts (GRBs). In
this paper, we estimate the neutrino flux and spectrum expected from these
sources, focusing on Type II SNe. We perform detailed calculations of
$p\gamma$ interactions, accounting for all the neutrino production channels
and scattering angles. We provide predictions of expected event rates for
operating neutrino telescopes, such as ANTARES and IceCube, as well as for the
future generation telescope KM3NeT. We find that for GRB energies channeled
into protons spanning between $\sim 10^{51}-10^{53}$ erg, choked GRBs may
substantially contribute to the observed astrophysical neutrino flux, if their
local rate is $\sim 80-1$ Gpc-3yr-1 respectively.

##  1 Introduction

The diffuse flux of high-energy neutrinos observed by the IceCube experiment
[1] has not shown yet the presence of a significant event clustering [2]. In
other words, the angular distribution of the events collected so far does not
highlight any specific location of the sky to search for correlated
counterparts, e.g. in the electromagnetic domain or in the cosmic-ray
anisotropy. This might indicate that the majority of neutrino sources is of
extra-galactic origin.

The possible production mechanisms of neutrinos are respectively proton-proton
(pp) collisions and photo-meson (p$\gamma$) interactions. In both cases,
neutrinos would emerge from both the direct decay of charged pions as well as
from muon decay, while gamma rays result from the neutral pion sudden decay.
As such, within hadronic sources, a similar amount of high-energy
electromagnetic radiation accompanies the production of neutrinos.

However, in the presence of intense radiation fields intrinsic to the source,
the emitted gamma rays are likely to be absorbed in the pair production
process (see e.g. [3]). In addition, during their propagation towards Earth,
multi-TeV gamma rays can initiate electromagnetic cascades on both the extra-
galactic background light (EBL) and the cosmic microwave background (CMB), and
be reprocessed towards lower energies. By comparing the isotropic diffuse
gamma-ray background (IGRB), observed by Fermi [4], with the astrophysical
neutrino flux, a strong tension among the two was remarked [5], in the
hypothesis that the pp interaction in transparent gamma-ray sources is
originating both fluxes. This constraint motivates the alternative study of
p$\gamma$ sources, opaque to high-energy gamma rays, as main neutrino sources.
Choked GRBs belong to such a population [6]. These GRBs are believed to be
originated from core collapse of massive stars resulting in a relativistic
fireball jet, which is unable to break through the stellar envelope. A choked
jet deposits all its energy in a cocoon, that eventually breaks out from the
star releasing energetic material at sub-relativistic velocities. Such fast
moving material has a unique signature that can be detected in early time
supernova spectra (see [7, 8, 9]). Neutrino emissions from choked GRBs have
hence attracted much attention [10, 11, 12, 13], because of the possibility of
explaining the IceCube diffuse flux without incurring into inconsistencies
with the IGRB.

Estimates of neutrino fluxes from astrophysical sources as a result of the
$p\gamma$ interaction have been performed in the past both with analytical
[14, 15] and numerical calculations [16, 17, 18, 19, 20, 21], but no
simulation for choked GRBs sources has been attempted so far. To study the
expected neutrino flux on Earth from such $\nu$ bright (while $\gamma$ dark)
sources, we have developed a Monte Carlo (MC) code to simulate the particles
kinematics and interactions inside the jet, allowing for a detailed
description of the micro-physics involved in the photo-meson process. Once the
particle spectra at the source are obtained, our code calculates propagation
through cosmic distances from the source to Earth, thus allowing a reliable
estimate of the number of neutrino events expected in current and future
neutrino telescopes.

The paper is organized as follows: Section 2 briefly describes the jet
dynamics inside the source in order to characterize the interacting particle
spectra. Section 3 presents an outline of the main structure of our simulation
for the p$\gamma$ interaction. Calculations leading to the estimate of
neutrino events in telescopes from a single source are contained in Section 4.
Section 5 describes the prediction for the diffuse neutrino spectrum coming
from choked sources. Section 6 summarizes and discusses the results obtained.
Moreover, Appendix A illustrates the detailed kinematic calculations adopted
in the code.

##  2 Choked jets and primary particle spectra

At the end of their lives, massive stars collapse into a compact object, and
it is widely believed that aforesaid phenomena can produce energetic jets [22,
23]. GRBs are usually associated to type Ic SNe [24], which lack of the
hydrogen envelope and have lost part of the helium envelope as well. However,
much less is known concerning choked GRBs. There have been many investigations
in the literature about their progenitor nature. The case of Type Ib/c SNe,
corresponding to the scenario of failed GRBs, would imply a duration of the
central engine less than 100 s, which is typical for long GRBs [6, 25, 20, 21,
26]. Several authors have also considered the possibility that jets may be
choked in the stellar envelope of red/blue supergiant stars, hence associated
with Type II SNe (see e.g. [11], but also [6] consider that this model is
possible). In this case the duration of the central engine could be longer
than for typical long GRBs [27, 28]. Note that this model may be problematic
as most of the ordinary SNe may not have rapid rotation (see e.g. [29]). While
the details of the progenitor model remain uncertain, preliminary calculations
suggest that a relativistic jet can be launched along the progenitor rotation
axis [6].

We consider a jet propagating inside the progenitor source and define
$t_{\rm{jet}}$ as the jet lifetime. Such a parameter strongly depends on the
physical process which powers the jet. We note that long GRBs are typically
observed with durations shorter than 100 s. However, choked GRBs might behave
significantly differently from standard ones. For instance, $10^{5}-10^{6}$ s
lifetimes are rather expected if the jet is powered by accretion onto a
central black hole or by the spin down luminosity of a magnetar [11]. Hence
$t_{\rm{jet}}$ should be considered as a free parameter of the model, possibly
constrained by e.g. neutrino observations.
While the jet is making its way through the star, it can be slowed down in a
termination shock to a Lorentz factor much lower than its original value. If
the jet crossing time $t_{\rm{cross}}$ is longer than $t_{\rm{jet}}$, then the
jet will result choked inside the stellar envelope, namely it will not be able
to break through the star surface [22]. In other words, when the central
engine powers the jet for a shorter time than the one required to travel
across the star radius, the jet is stalled inside the source. Since the jet
crossing time reads as

$t_{\rm cros}\simeq 1.1\times
10^{4}\,\rm{s}\,R^{2}_{13.5}\,L_{\mathrm{iso,50}}^{-1/2}\,\rho_{\rm
H,-7}^{1/2}\,,$ (2.1)

it is clear that extended stellar envelopes of radius $R$ are required to
choke a luminous jet. Here (and elsewhere) $A_{n}$ denotes $A/10^{n}$, e.g.
$L_{\mathrm{iso,50}}^{1/8}=(L_{\rm{iso}}/10^{50})^{1/8}$, and cgs units are
used for normalization. This might be the case of type II SN progenitors,
namely core collapse of red supergiant stars, possessing thick stellar
envelopes (see e.g. [6] for a detailed description of this source class). Such
stars are surrounded by helium and hydrogen envelopes, which have radii
respectively equal to $r_{\rm{He}}\sim 10^{11}\,\rm{cm}$ and
$r_{\rm{H}}\gtrsim 10^{13}\,\rm{cm}$. Note that the source density in the
hydrogen envelope is low enough, $\rho_{\rm{H}}\sim 10^{-7}\,\rm{g/cm^{3}}$,
as to allow for particle acceleration to relativistic velocities, possibly
followed by photomeson interactions beyond the helium radius. Alternatively,
choked GRBs might originate from low luminous jets propagating into thinner
stellar structures, as those belonging to massive stars which have lost a
substantial fraction of their external envelopes through their final
evolutionary stages. While the origin of choked jets remains so far elusive,
in this work we will consider type II SN as a reference scenario, by adopting
within our simulation the following benchmark values characterizing the jet
features: the Lorentz factor $\Gamma=100$, the lifetime
$t_{\rm{jet}}=10^{3}\,\rm{s}$ and the isotropic luminosity
$L_{\rm{iso}}=10^{50}\,\rm{erg/s}$, where the two latter values provide a GRB
isotropic energy given by $E_{\rm{iso}}=t_{\rm{jet}}\times
L_{\rm{iso}}=10^{53}\,\rm{erg}$. We plan to investigate the type Ib/c scenario
in a future work.

Internal shocks (IS) are expected to develop in the jet interior, inward of
the termination shock, at a distance from the inner engine equal to

$R_{\rm IS}\simeq 2\Gamma^{2}c\delta t\simeq 6\times 10^{12}~{}{\rm
cm}\,\Gamma_{2}^{2}\,\delta t_{-2}\,,$ (2.2)

where $\delta t$ represents the minimum variability timescale of the central
source. Depending on the initial opening angle $\theta$ of the jet and the
properties of the environment where it propagates, the jet can result in a
conical shape or in a cylindrical one [30]. The former case, also known as
collimated jet, is verified if the jet’s head Lorentz factor $\Gamma_{\rm
h}\sim 1$ is smaller than $\theta^{-1/3}$. A non relativistic jet head
satisfies the above condition for e.g. $\theta=0.2$ rad, which we will
consider in the following. In addition, in such a case, the location of the
jet head at the end of the jet lifetime can be approximated with [11]

$R_{\rm h}\simeq 9\times 10^{12}~{}{\rm
cm}\,L_{\mathrm{iso,50}}^{1/4}\,t_{\rm{jet,\,3}}^{1/2}\,\rho_{\rm
H,-7}^{-1/4}\,,$ (2.3)

In between internal shocks and the jet head position, the collimation shock
(CS) separates the shocked jet from the unshocked one, and it is found at [31,
20]

$R_{\rm CS}\simeq 8.2\times 10^{12}~{}{\rm cm}\,t_{\rm
jet,3}^{2/5}\,L_{\mathrm{iso,50}}^{3/10}\,(\theta_{\rm
j}/0.2)^{2/5}\,\rho_{\rm H,-7}^{-3/10}\,.$ (2.4)

While the collimation shock is most likely radiation dominated, this is not
the case for internal shocks, that are instead collisionless, and can hence
accelerate protons and electrons to $\gtrsim 10^{5}$ GeV [6]. In fact,
internal shocks are such that the co-moving size of the upstream flow is
smaller than the mean free path of thermal photons in the upstream flow [20].
This constraint translates into a conservative bound of

$l_{\rm u}<l_{\rm dec}\longrightarrow\tau=n_{\rm p,IS}\sigma_{\rm T}(R_{\rm
IS}/\Gamma)<1\,,$ (2.5)

which for our benchmark parameters is verified as

$\tau=0.013\Gamma_{2}^{-3}L_{\mathrm{iso,50}}^{3/4}\,t_{\rm{jet,\,3}}^{-1/2}\,\rho_{\rm
H,-7}^{1/4}<1\,.$ (2.6)

Hence efficient acceleration of protons is expected to be realized at internal
shocks, as also indicated from the comparison among the acceleration time of
protons $t_{\rm acc}$ and the characteristic dynamical time of the system,
e.g. $t_{\rm jet}$. The former reads as [11]

$t_{\rm acc}(E_{p})=\phi\frac{E_{p}}{q_{e}Bc}\simeq
0.05\,\rm{s}\,\phi_{1}\,\epsilon^{-1/2}_{b,-1}\,E_{p,15}\,L_{\mathrm{iso,50}}^{-1/4}\,\Gamma_{2}\,t_{\rm{jet,\,3}}^{1/2}\,\rho_{\rm
H,-7}^{-1/4}\,,$ (2.7)

where $\phi$ is the number of gyro-radii required to e-fold the particle
energy (here given in eV units). This estimate shows that multi-PeV protons
can actually be accelerated at internal shocks within a time much shorter than
the jet lifetime. Note that the strength of the magnetic field in the internal
shock reads as

$B\simeq 2.9\times
10^{4}\,\rm{G}\,\epsilon^{-1/2}_{b,-1}\,\Gamma_{2}^{-1}\,L_{\mathrm{iso,50}}^{1/4}\,t_{\rm{jet,\,3}}^{-1/2}\,\rho_{\rm
H,-7}^{1/4}\,.$ (2.8)

Though little is known on the acceleration process in GRBs, we assume a first
order Fermi process for primary charged particles [32], producing a
characteristic power-law differential energy distribution as:

$\frac{dN_{p}}{dE_{p}}=k_{\mathrm{MC}}\,E_{p}^{-2}\,,$ (2.9)

where $k_{\mathrm{MC}}$ is a normalization factor which will be discussed in
detail in Section 4. Electrons are expected to lose all their energy into
synchrotron radiation. However, because of the large optical thickness due to
Thomson scattering, these photons thermalize in the reverse shock region (RS)
to a temperature $T_{\gamma}$ such that [25]:

$kT_{\gamma}\simeq
313\,\rm{eV}\,\,L_{\rm{iso,\,50}}^{1/8}\,\epsilon_{e,\,-1}^{1/4}\,t_{\rm{jet,\,3}}^{-1/4}\,\rho_{\rm
H,-7}^{1/8}\,,$ (2.10)

$k$ being the Boltzmann constant, and $\epsilon_{e}$ the fraction of
$E_{\rm{iso}}$ deposited into electrons. For the fraction of energy channeled
into electrons and magnetic field we consider the equipartition condition
$\epsilon_{\rm e}=\epsilon_{\rm b}=0.1\,E_{\rm{iso}}$.

Inside the RS, the photon spectrum is well represented by a blackbody
distribution, reading as:

$\frac{dN_{\gamma}}{dE_{\gamma}}=\frac{8}{19(kT_{\gamma})^{3}}\frac{E_{\gamma}^{2}}{e^{\frac{E_{\gamma}}{kT_{\gamma}}}-1}\,,$
(2.11)

with a peak in the energy distribution given by the Wien’s displacement law

${E}^{\mathrm{\,max}}_{\gamma,\,\mathrm{RS}}\simeq 2.82\,kT_{\gamma}\sim
881\,\rm{eV}\,,$ (2.12)

and a number density following the Planck distribution

$n_{\gamma,{\rm
RS}}=\int_{0}^{\infty}\frac{dN_{\gamma}}{dE_{\gamma}}dE_{\gamma}=16\pi\xi(3)\left(\frac{kT_{\gamma}}{hc}\right)^{3}\simeq
19.2\pi\left(\frac{kT_{\gamma}}{hc}\right)^{3}\,,$ (2.13)

where $\xi(3)\simeq 1.202$ is the Riemann zeta function. A fraction of thermal
photons will then escape in the IS, since it is optically thin, where their
peak energy becomes:

${E}^{\mathrm{\,max}}_{\gamma,\,\mathrm{IS}}=\Gamma_{\rm{IR}}{E}^{\mathrm{\,max}}_{\gamma,\,\mathrm{RS}}\sim
88.1\,\mathrm{keV},$ (2.14)

where $\Gamma_{\rm{IR}}\sim\Gamma$ (if $\Gamma\gg 1$ [11]) is the Lorentz
factor of the IS with respect to the RS. The fraction $f_{\rm esc}$ of photons
getting in the IS frame by Thomson scattering from the location $R_{\rm h}$ of
the jet’s head can be computed following [20], namely

$f_{\rm esc}=\frac{1}{\tau_{\rm h}}=\frac{1}{n_{\rm h}R_{\rm h}\sigma_{\rm
T}}\,,$ (2.15)

where $n_{\rm h}$ represents the number density of electrons in the shocked
head region, which is related to the number density $n_{\rm j}$ of the
unshocked jet through the compression ratio at the shock, i.e. $n_{\rm
h}=4\Gamma n_{\rm j}$ for a strong shock and $n_{\rm j}=L_{\rm iso}/(4\pi
c\Gamma^{2}R_{\rm h}^{2}m_{\rm p}c^{2})$. Considering the jet head radius from
Eq. (2.3) we can derive

$f_{\rm esc}=1.9\times
10^{-3}\,\Gamma_{2}L_{\mathrm{iso,50}}^{-3/4}\,t_{\rm{jet,\,3}}^{1/2}\,\rho_{\rm
H,-7}^{-1/4}\,.$ (2.16)

Hence, the photon density in the internal shock frame reads as

$n_{\gamma,{\rm IS}}=\Gamma n_{\gamma,{\rm RS}}f_{\rm esc}\simeq 1.9\times
10^{20}\,\rm{cm^{-3}}\,\Gamma_{2}^{2}\,L_{\mathrm{iso,50}}^{-3/8}\,t_{\rm{jet,\,3}}^{-1/4}\,\epsilon_{e,\,-1}^{3/4}\,\rho_{\rm
H,-7}^{1/8}\,.$ (2.17)

##  3 Photomeson interaction and neutrino production

The production of charged and neutral mesons inside a jet can take place
through the resonant production of a $\Delta^{+}$ by highly relativistic
protons that interact with photons of the ambient radiation field, namely the
thermal photons from the jet head. The code simulates the photo-meson
interaction for each pairs of protons and photons which satisfy the threshold
condition. The secondary particles emerging from different interaction
channels are also simulated. Their energy losses in the form of synchrotron
radiation, inverse Compton scattering off the dense thermal radiation field,
synchrotron self Compton and adiabatic losses are accounted for. As we will
show in the following, we identify synchrotron and adiabatic losses as the
main channel for muon losses, while negligible losses affect charged pions. In
order to account for possible energy losses of protons due to photomeson
interactions, we compared the proton acceleration time with the $p\gamma$ loss
time, the latter reading as

$t_{\rm cool,p\gamma}=\frac{1}{n_{\gamma,\rm
IS}c\sigma_{p\gamma}K_{p\gamma}}\simeq
0.0029\,\rm{s}\,\epsilon^{-3/4}_{e,-1}\,\Gamma_{2}^{-2}\,L_{\mathrm{iso,50}}^{3/8}\,t_{\rm{jet,\,3}}^{1/4}\,\rho_{\rm
H,-7}^{-1/8}\,,$ (3.1)

for a cross section $\sigma_{p\gamma}=1\times 10^{-28}$ cm2 and inelasticity
$K_{p\gamma}=0.6$ [33]. From such a comparison, we obtain a proton cut-off of

$E^{\rm IS}_{\rm cut,15}\simeq
0.066\,\phi^{-1}_{1}\,\epsilon^{1/2}_{b,-1}\,\epsilon^{-3/4}_{e,-1}\,L_{\mathrm{iso,50}}^{5/8}\,\Gamma^{-3}_{2}\,t_{\rm{jet,\,3}}^{-1/4}\,\rho_{\rm
H,-7}^{1/8}\,,$ (3.2)

where again the proton energy is given in eV units, namely a cut-off at
$E_{\rm cut}=66$ TeV arises in the IS frame. This affects the differential
energy spectrum of protons, which is modified for $E_{\rm p}>E_{\mathrm{cut}}$
with respect to the simple power law of Fermi acceleration given in Eq. (2.9),
as

$\frac{dN_{\rm p}}{dE_{\rm
p}}=k_{\mathrm{MC}}\,E_{p}^{-2}\,\,e^{-\frac{(E_{\rm
p}-E_{\mathrm{cut}})}{2E_{\mathrm{cut}}}}\,.$ (3.3)

The kinematics of particles emerging from the $\Delta^{+}$ production is
evaluated in the center of mass rest frame, and then a Lorentz boost is
applied to transform them to the IS frame through the center of mass Lorentz
factor $\gamma^{*}$. All the center of mass quantities are indicated with an
asterisk $*$. Ultimately, the observed energies at the source are obtained by
boosting the particles by the jet Lorentz factor $\Gamma$. In the center of
mass frame we use natural units where $\hbar$ = c = 1. Moreover, we indicate
with $\beta_{a}$ the speed of the generic particle $a$
($\gamma_{a}=1/\sqrt{1-\beta_{a}^{2}}$), and with $E_{a}$, $p_{a}$, $m_{a}$
its energy, momentum and mass. For detailed kinematic calculations see
Appendix A.

Once a proton and a photon are extracted, in order to proceed with the
interaction a condition must be met, i.e. the production of the $\Delta^{+}$
resonance, given by the threshold requirement (see Eq. (3.4)). In the IS, the
condition for a successful interaction can be written as:

$(m_{p}+m_{\pi})^{2}=({E}^{\rm{IS}}_{p}+{E}^{\rm{IS}}_{\gamma})^{2}-({\vec{p}}^{\,\rm{IS}}_{p}+{\vec{p}}^{\,\rm{IS}}_{\gamma})^{2}$
$m_{p}^{2}+m_{\pi}^{2}+2m_{p}m_{\pi}=m_{p}^{2}+2{E}^{\rm{IS}}_{p}{E}^{\rm{IS}}_{\gamma}(1-\cos\theta)$
${E}^{\rm{IS}}_{p}\geq\frac{{m_{\pi}}^{2}+2m_{p}m_{\pi}}{2{E}^{\rm{IS}}_{\gamma}(1-\cos\theta)}\,,$
(3.4)

where $\theta$ indicates the angle between protons and photons in the IS. When
the threshold requirement in Eq. (3.4) is not satisfied, a new photon is
extracted. In fact, as a consequence of the extremely high number density of
thermal photons (see Eq. (2.17)), all generated protons will eventually
interact. In order to compute the proton path inside the thermal photon gas,
it is necessary to invert the path probability distribution:

$\rm{dP(x_{\rm{p}})}=\frac{1}{\lambda_{p\gamma}}\
e^{-x_{\rm{p}}/\lambda_{p\gamma}}\ \rm{dx_{\rm{p}}}\,,$ (3.5)

where $\lambda_{p\gamma}=1/(\sigma_{p\gamma}n_{\gamma})$ is the interaction
length. The proton path obtained through the inversion of Eq. (3.5) peaks at
$\sim 10^{7}$ cm for $\sigma_{p\gamma}=550\,\mu\rm{b}$ (the $\Delta^{+}$
resonance peak), while for the out of resonance region value of
$\sigma_{p\gamma}$, namely $\sim 120\,\mu\rm{b}$, the proton path peaks at
$\sim 10^{8}$ cm. This means that a proton will interact, on average, after
travelling $10^{7.5}$ cm inside the source, which has an external hydrogen
envelope radius of $10^{13}$ cm. Thus, we can safely assume that all protons
interact inside the source. In other words, choked jets behave as pure
calorimeters with respect to protons above the pion production thresholds,
such that all their energy goes into secondary light particles, including
gamma rays and neutrinos. This result is consistent with findings from other
authors as well (see e.g. [21]). The high-energy radiation emitted is expected
to be reprocessed at lower energies in pair production processes with existing
jet radiation field. Hence observations of neutrinos would provide a unique
constrain on the energy channeled into accelerated protons.

Inside the IS region, the following relation is also satisfied:

$\cos{\theta}=\frac{{-p}^{\rm{IS}}_{L}}{\sqrt{({p}^{\rm{IS}}_{L})^{2}+({p}^{\rm{IS}}_{T})^{2}}}\,,$
(3.6)

${p}^{\rm{IS}}_{L}$ and ${p}^{\rm{IS}}_{T}$ being, respectively, the
longitudinal and transverse momenta of the photons in the IS. Defining $\xi$
as the angle between protons and photons in the RS, we obtain the following
algebraic equations for the photon momenta:

$-{p}^{\rm{IS}}_{L}=-\Gamma_{\rm{IR}}\ E_{\gamma,\,\rm{RS}}\ (\cos\xi\ +\
\beta_{\rm{IR}})\,,$
${p}^{\rm{IS}}_{T}={p}^{\rm{RS}}_{T}=\sqrt{E_{\gamma,\,\rm{RS}}^{2}\
(1-\cos^{2}\xi)}\,.$

For all simulated interactions, the cosine value obtained from Eq. (3.6) is
$\sim 1$, meaning that the collisions are realized head-on.

Once the interaction occurs, the Monte-Carlo algorithm takes into account the
various interaction channels. Within our simulation we considered three
possible channels for the interaction, depending on the photon energy in the
proton’s rest frame, $\epsilon_{r}$:

$s^{\rm{IS}}=2{E}^{\rm{IS}}_{p}{E}^{\rm{IS}}_{\gamma}(1-\cos\theta)\sim
4{E}^{\rm{IS}}_{p}{E}^{\rm{IS}}_{\gamma}=m^{2}_{p}+2m_{p}{\epsilon}_{r}\,,$

$s^{\rm{IS}}$ being the square of the invariant energy of the system in the IS
frame, such that

$\epsilon_{r}=\frac{4{E}^{\rm{IS}}_{p}{E}^{\rm{IS}}_{\gamma}-m_{p}^{2}}{2m_{p}}\,.$
(3.7)

The p$\gamma$ cross section as a function of $\epsilon_{r}$ is characterized
by a very peaked resonance region around the threshold, mainly due to the two
reaction channels we are studying (namely the ones giving the $p+\pi^{0}$ and
$n+\pi^{+}$ final states): in our simulation, the interaction probability
$P_{\rm{int}}$ in this region is set equal to 1. Above the threshold, the
cross section is dominated by multiple pion final states, and it becomes
approximately constant for $\epsilon_{r}>$ 10 GeV. For a detailed
characterization of the p$\gamma$ cross section see e.g. [34]. Hence, multiple
values were used for $P_{\rm{int}}$ within the code, depending on the cross-
section.

The boundaries between each interaction channel in our simulation are given by
the following numerical values:

  * •

$0.2\,\mathrm{GeV}\leq\epsilon_{r}<0.5\,\mathrm{GeV}$: this energy range
corresponds to the peak of the cross section, where the interaction occurs
with $P_{\rm{int}}=1$ and the $\Delta^{+}$ resonance is produced at rest:

$p+\gamma\to\Delta^{+}\to\begin{cases}n+\pi^{+}\\\ p+\pi^{0}\end{cases};$
(3.8)

  * •

$0.5\,\mathrm{GeV}\leq\epsilon_{r}<2\,\mathrm{GeV}$: such interval matches the
secondary peak range of the cross section, when the $\Delta^{+}$ resonance is
produced along with an additional pion ($\mathcal{N}$ is a nucleon, while both
$\pi^{\rm{A}}$ and $\pi^{\rm{B}}$ can be neutral, positively or negatively
charged pions):

$\displaystyle p+\gamma\to\pi^{\mathrm{A}}+$ $\displaystyle\Delta^{+}$
$\displaystyle\to\,\,\mathcal{N}+\pi^{\mathrm{B}}\,.$ (3.9)

Here, the interaction occurs with two possible values of probability:
$P_{\rm{int}}=0.6$ when $\epsilon_{r}\leq 1\,\rm{GeV}$, otherwise the
interaction probability is set equal to 0.4;

  * •

$2\,\mathrm{GeV}\leq\epsilon_{r}\leq 100\,\mathrm{GeV}$: this energy window
corresponds to the plateau in the cross section, in which the three types of
pions are simultaneously produced. This range is often referred to as
multipion channel, where “$\ldots$” indicates other particles created in the
interaction:

$p+\gamma\to\pi^{0}+\pi^{+}+\pi^{-}+\ldots$ (3.10)

In this region $P_{\rm{int}}=0.3$.

The code selects how to progress the interaction based on which range
$\epsilon_{r}$ occupies. The interested reader is referred to Appendix A for a
more detailed discussion on branching ratios and kinematics in the three
channels.
The pion decay products include leptons and photons:

$\displaystyle{\pi}^{+}$ $\displaystyle\to\mu^{+}+\nu_{\mu}\to
e^{+}+\nu_{e}+\bar{\nu}_{\mu}+\nu_{\mu}\,,$ (3.11) $\displaystyle{\pi}^{-}$
$\displaystyle\to\mu^{-}+\bar{\nu}_{\mu}\to
e^{-}+\bar{\nu}_{e}+\nu_{\mu}+\bar{\nu}_{\mu}\,,$ (3.12)
$\displaystyle{\pi}^{0}$ $\displaystyle\to\gamma+\gamma\,.$ (3.13)

The muon decay is taken into account, and characterized through the Michel
parameters, which describe the phase space distribution of leptonic decays of
charged leptons. However, energy losses in the dense radiation fields or in
the jet magnetic field might affect the spectrum of both pions and muons,
possibly suppressing the flux of emerging neutrinos, if these losses develop
on shorter timescales than the particle decay time. Concerning the muon, its
lifetime at rest is $\tau_{\mu}\sim 2.2\times 10^{-6}$ s, while for the
charged pion its lifetime at rest is $\tau_{\pi^{\pm}}\sim 2.6\times 10^{-8}$
s. Additionally, while the jet expands, part of the kinetic energy is lost in
the form of heat, resulting into further energy losses for the particles as
well. The characteristic timescale for adiabatic losses can be taken the same
as the dynamical timescale of the system, namely

$t_{\rm dyn}=\frac{R_{\rm IS}}{\Gamma c}=2\Gamma\delta t\simeq
2\,\rm{s}\,\Gamma_{2}\,\delta t_{-2}\,,$ (3.14)

shorter than for instance the jet lifetime ($t_{\rm jet}=1000$ s). This scale
has to be compared with the time taken to emit synchrotron radiation, or to
energize the thermal or electron-emitted synchrotron photons to high-energy,
through inverse Compton processes.
For synchrotron radiation, muons loose energy in a characteristic time of

$\tau_{\mu}^{\rm{sync}}(E^{\rm
IS}_{\mu})=0.89\,\rm{s}\,\,E^{\rm{IS},\,-1}_{\mu,\,15}\epsilon_{b,\,-1}^{-1}\,\Gamma_{2}^{2}\,L_{\mathrm{iso,50}}^{-1/2}\,t_{\rm{jet,\,3}}\,\rho_{\rm
H,-7}^{-1/2}\,,$ (3.15)

while charged pions of

$\tau_{\pi^{\pm}}^{\rm{sync}}(E^{\rm
IS}_{\pi^{\pm}})=2.60\,\rm{s}\,\,E^{\rm{IS},\,-1}_{\pi^{\pm},\,15}\epsilon_{b,\,-1}^{-1}\,\Gamma_{2}^{2}\,L_{\mathrm{iso,50}}^{-1/2}\,t_{\rm{jet,\,3}}\,\rho_{\rm
H,-7}^{-1/2}.$ (3.16)

For inverse Compton we note that the IC interaction for muons and pions always
proceed in the Klein Nishina regime, as these secondary particles emerge from
the p$\gamma$, which only occurs at energies above threshold. We follow the
approach from [6], where the IC loss time of a pion in the IS frame off the
thermal field reads as

$\tau^{\rm IC}_{\rm KN}(E_{\pi}^{\rm
IS})\simeq\frac{2}{m_{e}^{2}c^{4}}\frac{\Gamma^{2}E_{\pi}^{\rm
IS}\epsilon^{\rm RS}_{\gamma}}{c\sigma_{\rm T}n^{\pi}_{\gamma}}\,,$ (3.17)

where $n^{\pi}_{\gamma}=2.9\times 10^{17}$ cm-3 with our benchmark values and
$\sigma_{p\gamma}=5\times 10^{-28}$ cm2. Plugging numbers in, we obtain for
charged pions a timescale of

$\tau^{\rm IC}_{\rm KN}(E_{\pi}^{\rm IS})\simeq 1.2\times 10^{4}~{}{\rm
s}\left(\frac{E_{\pi}^{\rm IS}}{1~{}{\rm
TeV}}\right)\left(\frac{\Gamma}{100}\right)^{2}\left(\frac{\epsilon^{\rm
RS}_{\gamma}}{0.881~{}{\rm keV}}\right)\left(\frac{n^{\pi}_{\gamma}}{2.9\times
10^{17}~{}{\rm cm}^{-3}}\right)^{-1}\,.$ (3.18)

An analogous expression holds for the IC loss time of muons. Thus we can argue
that the effect of energy losses due to inverse Compton processes do not
affect the pions and muon decays.

Finally, for the evaluation of the loss time in SSC process of pions and muons
off the synchrotron radiation field emitted by primary electrons, we simply
note that it can always be neglected with respect to the synchrotron loss time
of the same particles, because of the equipartition condition. To summarize
the results obtained with the benchmark parameters we set, charged pions do
not suffer any significant loss process within their lifetime, namely they
always decay freely. In turn, muons with energies higher than $E_{\mu}^{\rm
IS}\simeq 100$ TeV decay with a lower probability because they are affected by
adiabatic losses, and finally above $E_{\mu}^{\rm IS}\simeq 500$ TeV by
synchrotron losses. These conditions are applied to muons in all interaction
channels, and the corresponding neutrino spectra obtained through the
simulation depend on such a suppression. It should be noted, however, that the
cut applied to protons at the level of $E_{\rm cut}\simeq 66$ TeV in the IS
frame reflects into a cut on the pion and muon energies stronger than the cut
due to their own adiabatic and synchrotron losses.

The particle spectra at the source obtained through our MonteCarlo simulation
are shown in Figure 1. In the top left panel we present the energy spectrum of
the interacting protons, as well as the secondary particles emerging from the
$\Delta^{+}\longrightarrow p+\pi^{0}$ decay channel. In this study we focus on
neutrino production, hence we do not propagate the resulting gamma rays
further in the jet, namely we neglect pair production processes of the high-
energy photons with the thermal radiation field of the jet. We defer
discussion about absorption of gamma rays in the chocked jet to a future work.
In the top right panel, we show the secondary particles from all $\Delta^{+}$
decays that give rise to $\pi^{+}$, while in the bottom panel we show the
spectrum of $\pi^{-}$ obtained from $\Delta^{+}$ decay in the multi-pion decay
channel. Note that the latter is extremely similar to the $\pi^{+}$ channel,
while it contains far less events due to the fact that this channel is not
present at the resonance peak, where most of the particles are produced. These
spectra are not normalized yet: as a consequence, the y axis is in arbitrary
units.

Figure 1: Particle spectra at the source arising from all interaction
channels, as seen in the co-moving reference frame. Top left: $\pi^{0}$
channel; Top right: $\pi^{+}$ channel; Bottom: $\pi^{-}$ channel.

##  4 Number of events expected in neutrino telescopes

Transient phenomena such as GRBs are very interesting for a potential neutrino
detection, since these analyses are almost background free given the well-
defined space-time window provided by multi-messenger constraints. To compute
the neutrino flux expected on Earth from a source at a cosmological distance,
we scale the MC-obtained fluxes for three factors: the energy channeled into
protons at the source, the dilution of the neutrino flux to Earth, and
neutrino oscillations.

The first factor comes from energetic considerations for the source: as
mentioned in Section 2, we considered a choked GRB with an isotropic energy
release of $E_{\rm iso}=10^{53}$ erg. We also assume that a comparable amount
of energy $E_{\rm p,jet}$ is delivered to the protons accelerated inside the
jet [6]. So, for the accelerated protons spectrum at the source in the IS
frame we have:

$k_{\rm MC}\int_{E_{\rm min}}^{E_{\rm max}}E_{\rm p}\frac{dN_{\rm p}}{dE_{\rm
p}}dE_{\rm p}=E_{\rm p,jet}\,,$ (4.1)

where $E_{\rm{min}}=100$ GeV and $E_{\rm{max}}=10^{8}$ GeV. Solving this
integral with the acceleration spectrum of protons, namely the pure power law
distribution given in Eq. (2.9), leads us to an estimate of the scaling factor
$k_{\mathrm{MC}}$, which provides the exact fluence of secondary particles at
the source in the IS frame, including gamma rays and neutrinos, normalized as
to take into account the amount of energy channeled into accelerated protons.
Moreover, particles travel through cosmological distances before arriving on
Earth, such that their flux undergoes a dilution during the propagation. This
factor can be calculated as $(1+z)^{2}/(4\pi d_{\rm L}(z)^{2})$ [35, 36],
$d_{\rm L}(z)$ being the source luminosity distance. For the choked GRB in our
model, we assume a redshift $z$ of $\sim 1$, from which we can evaluate
$d_{\rm L}(z=1)=2.06\times 10^{26}$ m.

Neutrino telescopes adopt muon neutrino and muon antineutrino fluxes at Earth
for astronomical studies, thanks to the good angular resolution obtained in
this event sample by reconstructing the long track length of the induced muon.
Neutrino flavor on Earth differs from the one at the source due to the effect
of neutrino oscillations. Since GRBs are extra-galactic sources, we mediate
neutrino oscillations over cosmic distances: starting from a flavor ratio of
$\nu_{e}$ : $\nu_{\mu}$ : $\nu_{\tau}$ = 1 : 2 : 0 at the source (as the
result of the $p\gamma$ interaction), one should expect a flavor ratio of 1 :
1 : 1 on Earth. The oscillation probabilities $P_{ff^{\prime}}$ of neutrinos
changing flavor ($f\to f^{\prime}$) are given by [37]:

$P_{ff^{\prime}}=\sum_{j}\left|U_{fj}^{2}\right|\left|U_{f^{\prime}j}^{2}\right|,$
(4.2)

so that the oscillated spectrum for the lepton of flavor $f$ reads as:

$E_{\nu_{f}}\frac{dN_{\nu_{f}}}{dE_{\nu_{f}}}=P_{ff}\left(E_{\nu_{f}}\frac{dN_{\nu_{f}}}{dE_{\nu_{f}}}\right)^{\rm{\,source}}+P_{f^{\prime}f}\left(E_{\nu_{f^{\prime}}}\frac{dN_{\nu_{f^{\prime}}}}{dE_{\nu_{f^{\prime}}}}\right)^{\rm{\,source}}\,,$
(4.3)

where the right-hand side quantities represent the generic neutrino spectrum
at the source in the co-moving frame, obtained through the MC simulation with
the $\Gamma$-Lorentz boost from the IS frame. This formula is symmetric in the
exchange $f\leftrightarrow f^{\prime}$, and it is valid for neutrinos and
antineutrinos.

Figure 2: All-flavor time-integrated neutrino flux on Earth. Left: neutrinos
from the resonance peak (green dashed line), from the region above the peak
(orange dot-dashed line) and their sum (red line) are compared in this plot.
Right: neutrinos from choked jets with different Lorentz factors.

Defining the total neutrino energy spectrum on Earth as

$E_{\nu}\frac{dN_{\nu}}{dE_{\nu}}=E_{\nu_{e}}\frac{dN_{\nu_{e}}}{dE_{\nu_{e}}}+E_{\nu_{\mu}}\frac{dN_{\nu_{\mu}}}{dE_{\nu_{\mu}}}+E_{\nu_{\tau}}\frac{dN_{\nu_{\tau}}}{dE_{\nu_{\tau}}}\,,$
(4.4)

the all-flavor time-integrated neutrino flux on Earth can be written as:

$E_{\nu}\varphi_{\nu}^{\rm{Earth}}=k_{\mathrm{MC}}\,\frac{(1+z)^{2}}{4\pi\,d_{\rm
L}(z)^{2}}E_{\nu}\frac{dN_{\nu}}{dE_{\nu}}\,.$ (4.5)

The all-flavor time-integrated neutrino flux on Earth can be seen in Figure 2,
where we also accounted for the redshift dependence of observed energy. In the
left panel, neutrinos coming from the delta resonance peak region
($0.2\,\mathrm{GeV}\leq\epsilon_{r}<0.5\,\mathrm{GeV}$) are compared with
neutrinos coming form the energy range exceeding $0.5\,\mathrm{GeV}$. This
result is particularly interesting because it indicates that focusing merely
on the $\Delta^{+}$ resonance peak for neutrino production leaves out a non
negligible number of secondary neutrinos coming from higher interaction
channels. Furthermore, the right panel of Figure 2 compares the all-flavor
neutrino spectra coming from choked jets with different Lorentz factors: the
higher the $\Gamma$, the higher the neutrino energy peak. Conversely, the
absolute normalization of the neutrino time-integrated flux
$E_{\nu}\varphi_{\nu}$ (see Eq. (4.5)) is independent of $\Gamma$, as $E_{\rm
p,jet}$ was set in the IS frame.
We are interested in estimating the number of muon neutrino events expected
from an individual source in different neutrino telescopes, since these
constitute the best channel of events for astronomical studies due to their
excellent angular resolution. Hence, we proceed by oscillating neutrinos
towards Earth, as to obtain the flux of muon neutrino flavor. The values for
the oscillation probabilities $P_{\mu\mu}$ and $P_{e\mu}$ to be used in Eq.
(4.3) are estimated using the latest results for the mixing angles, for
neutrino mass hierarchy in the case of normal ordering [38]. The muon neutrino
time-integrated neutrino flux observed on Earth obtained within our simulation
is shown in Figure 3, by using the following set of benchmark parameters:
$z=1$, $E_{\rm{p,\,jet}}=E_{\rm{iso}}=10^{53}$ erg, $L_{\rm{iso}}=10^{50}$
erg/s, $t_{\rm jet}=10^{3}$ s, and $\Gamma=100$.

Figure 3: Muon neutrino and antineutrino time-integrated flux expected on
Earth from a single source at z=1. Oscillations of neutrinos are here
accounted for.

Using the energy distribution of muon neutrinos thus obtained, an estimate of
the number of cosmic neutrino events expected in neutrino telescopes can be
performed. Convolving the detector effective area for neutrinos,
$A_{\mathrm{eff}}^{\nu}(E_{\nu},\delta)$, coming from declination $\delta$,
with the differential number of neutrinos per surface element expected on
Earth from the source, the number of events during the GRB event
($t_{\rm{jet}}=1000\,\rm{s}$) is obtained through the following relation:

$N_{\rm events}(\delta)=\int\frac{dN_{\rm
events}}{dE_{\nu}}(E_{\nu},\delta)\,dE_{\nu}=\int
A_{\mathrm{eff}}^{\nu}(E_{\nu},\delta)\,\Phi_{\nu}^{\rm{Earth}}(E_{\nu})\,dE_{\nu}\,.$
(4.6)

This evaluation has been performed for ANTARES [39], KM3NeT-ARCA [40] and
IceCube [41]. Note that the effective area adopted in this work for KM3NeT
refers to the trigger level, being the only available information at the time
of writing. In turn, the effective areas adopted for ANTARES and IceCube refer
to the analysis level, as these result from detailed studies performed in the
search for neutrino sources by the same Collaborations. The expected number of
events per energy bin is shown as a function of the energy in Figure 4, while
the total number of events expected in each neutrino telescope during the
choked GRB event is given in Table 1. The calculations show that current
generation detectors have non negligible chances to detect neutrinos from
individual sources of this kind, in the assumption that choked GRBs can
efficiently convert their kinetic energy into accelerated protons. In
particular, if the rate of occurrence of such explosions in the Universe is
high enough, it would be possible for the population of choked GRBs to
contribute to the observed diffuse astrophysical neutrino flux, whose sources
remain so far unidentified.

Figure 4: Number of expected events in ANTARES [39], KM3NeT-ARCA [40] and IceCube [41] from a characteristic choked GRB located at different declinations. Detector | $\delta$ | $N_{\mathrm{events}}$
---|---|---
ANTARES | $0^{\circ}<\delta<45^{\circ}$ | 2$\times 10^{-3}$
| $-45^{\circ}<\delta<0^{\circ}$ | 3$\times 10^{-3}$
| $-90^{\circ}<\delta<-45^{\circ}$ | 5$\times 10^{-3}$
KM3NeT-ARCA | Mean $\delta$ | 2$\times 10^{-1}$
IceCube | $0^{\circ}<\delta<90^{\circ}$ | 1$\times 10^{-1}$
| $0^{\circ}<\delta<30^{\circ}$ | 2$\times 10^{-1}$
Table 1: Expected number of muon neutrino events in several detectors from an
individual choked GRB with $L_{\rm iso}=10^{50}$ erg/s, t${}_{\rm jet}=1000$
s, $E_{\rm{p,\,jet}}=10^{53}$ erg in the IS frame and $\Gamma=100$ located at
different declination bands. Effective areas for each detector were taken from
[39, 40, 41].

##  5 Diffuse neutrinos from choked GRBs

To estimate the diffuse neutrino spectrum, we assume that the rate of choked
jets $R(z)$ at redshift $z$ follows the star formation rate even if there are
indications that the evolution of normal GRBs may deviate from the star
formation rate $\rho(z)$ [42, 43], i.e.,

$R(z)=R_{0}\rho(z),$ (5.1)

where the star formation rate is [44]

$\rho(z)=\frac{(1+z)^{2.7}}{1+[(1+z)/2.9]^{5.6}}\,,$ (5.2)

and $R_{0}$ is the local rate of choked GRBs in $\rm{Gpc}^{-3}\,\rm{yr}^{-1}$
. In this paper we do not consider a luminosity function for choked GRBs, but
we assume a constant luminosity. A more detailed study on the redshift
distribution and luminosity function effects on the diffuse neutrino flux is
out of the aim of the present paper but it will be a subject for future works.

The diffuse neutrino flux is calculated via integrating the neutrino spectrum
over the redshift from 0 to 8, i.e. [36]:

$E_{\nu_{\mu}}^{{\rm obs}}\phi_{\nu_{\mu}}(E_{\nu_{\mu}}^{{\rm
obs}})=\frac{c}{4\pi
H_{0}}\int_{0}^{8}\,E_{\nu_{\mu}}\frac{dN_{\nu_{\mu}}}{dE_{\nu_{\mu}}}((1+z)E^{{\rm
obs}}_{\nu_{\mu}})\frac{\frac{\Omega}{4\pi}R_{0}\rho(z)dz}{(1+z)\sqrt{\Omega_{\Lambda}+\Omega_{\rm
M}(1+z)^{3}}}\,,$ (5.3)

where the cosmological parameters are adopted as $H_{0}=70~{}{\rm
km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{\rm M}=0.3$, and $\Omega_{\Lambda}=0.7$,
$dN_{\nu_{\mu}}(E_{\nu_{\mu}})/dE_{\nu_{\mu}}$ is the differential neutrino
spectrum at the source, while $\Omega$ is the solid angle of the jet,
$\Omega=2\pi(1-\cos\alpha)$, and $\alpha=0.2$ rad $=11.5^{\circ}$ is its
aperture.
In order to obtain the best value for the rate of choked GRBs, we compared the
predicted flux of neutrino on Earth with the diffuse neutrino flux measured by
IceCube. As experimental data we used the results obtained by IceCube both
with the analysis of HESE events [45] and the parametrization of the
astrophysical diffuse muon neutrino flux [46]. The flux predicted with the
simulation was evaluated for several values of $R_{0}$: the best agreement
between data and our model has been obtained for $R_{0}=(1.0\pm 0.4)$ Gpc-3
yr-1 (with a reduced $\chi^{2}\simeq 2.5$). Fig. 5 shows this comparison: the
central blue line represents the flux foreseen for the selected value of
$R_{0}$, while the blue band represents the 68$\%$ confidence level
uncertainty on $R_{0}$.
Recent IceCube results [47] indicate that transient sources characterized by
hard spectra ($dN_{\nu}/dE_{\nu}\propto E_{\nu}^{-2.13}$), with a local rate
larger than $\sim 50$ Gpc-3yr-1, might possibly originate the observed diffuse
neutrino flux. Within our model, this constraint implies that the energy
converted into proton kinetic energy should be smaller than few times
$10^{51}$ erg, though this constraint is quite sensitive to the assumed
spectral shape of neutrinos and should be hence considered with caution. To
test the case of a less energetic choked GRB population being able to explain
the observed cosmic neutrino flux, we also performed a simulation
characterized by the following values of model parameters: $L_{\rm
iso}=10^{49}$ erg/s, $t_{\rm jet}=200$ s, $\Gamma=100$, $t_{\rm var}=0.005$ s
and $E_{\rm{p,\,jet}}=2\times 10^{51}$ erg in the IS frame. Such a realization
is consistent with the model presented in Sec. 2, though it differs in the
maximum energy of accelerated protons, which now reads $E^{\rm IS}_{\rm
cut}\simeq 25$ TeV, and in the maximum energy of muons, amounting to
$E_{\mu}^{\rm IS}\simeq 50$ TeV because of synchrotron cooling. As a result,
we obtained the yellow curve shown in Fig. 5 for the diffuse neutrino flux
contribution emerging from such a source population, which is at the level of
the IceCube data provided a local rate equal to $R_{0}=(80\pm 30)$ Gpc-3 yr-1.
The latter result is also in agreement with the transient analysis multiplet
constraint for hard spectra sources [47].

Figure 5: Diffuse neutrino spectra from choked GRBs expected on Earth,
compared with the IceCube results on the HESE data (black squares) [45] and on
$\nu_{\mu}$ tracks (red line) [46]. The blue line refers to the GRB simulation
with $L_{\rm iso}=10^{50}$ erg/s, $t_{\rm jet}=10^{3}$ s and $\Gamma=100$ that
is in better agreement with IceCube data as a function of $R_{0}$. The yellow
line refers to the simulation with $L_{\rm iso}=10^{49}$ erg/s, $t_{\rm
jet}=200$ s and $\Gamma=100$, again in better agreement with the observed
cosmic neutrino flux as a function of $R_{0}$. Note that the inferred local
rate of choked GRBs is consistent with the multiplet constraint from hard
transient sources [47] only in the less energetic case. In turn, a model with
local rate of 1 Gpc-3 yr-1 is ruled out by the same constraints [12].
Additionally, the dashed grey line shows the Waxmann & Bahcall upper bound
[48].

##  6 Results and Discussion

Recent stacking analyses performed by the current neutrino telescopes lead to
the conclusion that classical GRBs are not the main sources of the diffuse
astrophysical neutrino flux detected by IceCube [49, 50]. In addition, also
targeted searches for neutrinos in spatial and timing coincidence with the
prompt emission of those GRBs reported evidence for a lack of correlation in
neutrino data [51]. Such results have increased even more the interest in an
alternative GRB class opaque to gamma rays, choked GRBs, which would also be
consistent with the IGRB observed by Fermi [4].
We have estimated the neutrino fluxes from individual choked GRBs and the
number of expected events for present and future neutrino telescopes. This
result was obtained for the following set of benchmark parameters: $z=1$,
$L_{\rm{iso}}=10^{50}$ erg/s, $\Gamma=100$, $t_{\rm jet}=1000$ s and assuming
that $E_{\rm p,jet}=10^{53}$ erg in the IS frame. We find that, for the
declination in which each detector performs better, the number of muon
neutrino events expected from an individual choked GRB source is equal to
$5\times 10^{-3}$ for ANTARES, and $2\times 10^{-1}$ for IceCube and KM3NeT-
ARCA (see Table 1). These estimates scale almost linearly with the amount of
energy channeled into protons. The number of neutrino events expected from
each of these sources is in agreement with the lack of spatial anisotropies in
the sky map distribution of IceCube neutrinos. Moreover, the lack of
association of the latter events with known electromagnetic counterparts is to
date consistent with an extra-galactic origin of such neutrinos from photon-
dark sources, like choked GRBs.
We have also estimated the diffuse neutrino flux from the population of choked
GRBs and compared it to the IceCube latest results [45, 46]. Through a best
fit procedure we found that, if the energy that goes into accelerated protons
is $E_{\rm p,jet}\sim 10^{51}-10^{53}$ erg in the IS frame, then IceCube data
are reproduced for a local rate of choked GRBs of about $80-1$ Gpc-3yr-1,
respectively. Note that, if $E_{\rm p,jet}$ is decreased, the local rate of
choked GRBs has to increase in order to reproduce the IceCube data: this is
consistent with what found in the literature by other authors [52, 53].
In order to associate the origin of the IceCube data to choked GRBs, the
“smoking gun” would be the detection of the latter at electromagnetic bands
different than gamma rays. Preliminary studies have shown that these sources
may emit in the X-ray, optical and UV band [54, 55, 56]. We plan to extend
this work and explore the possibility to perform a follow up of choked
neutrino sources in these bands.

## Acknowledgments

The authors thank the anonymous referee for the constructive comments, which
significantly improved the manuscript. AC, SC, IDP and AZ acknowledge the
support from Sapienza Università di Roma through the grant ID
RM120172AEF49A82.

##  Appendix A Appendix - Analytic calculations for photo-meson interactions
of protons with thermal photons

To develop the kinematics relations used in the code, we consider three
reference systems: the source co-moving reference system, the IS reference
system, and the center of mass system (these kinematic quantities are
respectively indicated: plain, with the superscript IS, with an asterisk *).
We use natural units where $\hbar$ = c = 1, and we indicate with $\beta_{a}$
the speed of a generic particle ($\gamma_{a}=1/\sqrt{1-\beta_{a}^{2}}$), with
$E_{a}$ its energy and with $m_{a}$ its mass. The invariant energy of the
system, $\sqrt{s}$, is conserved and defined by the energy-momentum four-
vector:

$s=(E_{\gamma}+E_{p})^{2}-(p_{\gamma}+p_{p})^{2}=$
$=E_{p}^{2}+E_{\gamma}^{2}+2E_{p}E_{\gamma}-p_{p}^{2}-E_{\gamma}^{2}-2p_{p}E_{\gamma}\cos\theta=$
$=m_{p}^{2}+2E_{p}E_{\gamma}(1-\cos\theta)\,.$ (A.1)

This relation holds in all reference frames. The minimum energy for a proton
in the IS frame to produce a $\Delta^{+}$ resonance when interacting with a
photon of energy $E^{\rm{IS}}_{\gamma}$ is shown in Eq. (3.4): this is the
threshold condition to produce the $\Delta^{+}$ resonance.

The center of mass Lorentz factor is given by:

$\gamma^{*}=\frac{E^{\rm{IS}}_{p}+E^{\rm{IS}}_{\gamma}}{\sqrt{s}}\simeq\frac{E^{\rm{IS}}_{p}}{\sqrt{s}}\,.$
(A.2)

Given the enormous energies of the protons with respect to the photons, the
angle between the direction along which the Lorentz boost of the center of
mass is considered and the propagation direction of the proton is negligible.
For simplicity, we will consider the Lorentz boost along the direction of
observation.

To get the particles spectra in the source co-moving frame, one needs to boost
the energies in the IS by the jet’s Lorentz factor, $\Gamma$. In this
reference system the following relation holds:

$E_{p}=\Gamma\,E^{\rm{IS}}_{p}.$ (A.3)

Note that these quantities have not yet been corrected for the redshift due to
cosmological distances.

###  A.1 Delta Resonance Region

$p+\gamma\to\Delta^{+}\to\begin{cases}n+\pi^{+}\\\ p+\pi^{0}\end{cases}\,.$
(A.4)

In this region the invariant energy of the system $\sqrt{s}$ is fixed equal to
the $\Delta^{+}$ resonance mass. Once formed, the resonance decays at rest in
the center of mass system, and since the Lorentz boost is enormous, the decay
products in the IS system continue along the same direction of the incident
proton. At the resonance, the simulation code, takes into account the Clebsch-
Gordon coefficients, and the two decay channels are produced with the
following ratio:

$\frac{\mathrm{BR}(p\gamma\to p\pi^{0})}{\mathrm{BR}(p\gamma\to
n\pi^{+})}=\frac{2/3}{1/3}=2\,.$ (A.5)

By solving the interaction kinematics in this energy region (which corresponds
to $0.2\,\mathrm{GeV}\leq\epsilon_{r}<0.5\,\mathrm{GeV}$), we evaluate the
energies of the secondary particles in the center of mass frame:

$E_{\pi}^{*}=\sqrt{(p_{\pi}^{*})^{2}+m_{\pi}^{2}}=\frac{s+m_{\pi}^{2}-m_{N}^{2}}{2\sqrt{s}}\,,$
(A.6) $E_{N}^{*}=\frac{s-m_{\pi}^{2}+m_{N}^{2}}{2\sqrt{s}}\,,$ (A.7)

N being the nucleon produced together with the pion (proton or neutron).

In order to calculate the pion energy in the IS reference frame we must apply
the Lorentz transformations that are given by (considering the IS moving
towards the center of mass system along the boost direction, the x axis, with
speed -$\beta^{*}$):

$\begin{cases}p^{\rm{IS}}_{x}=\gamma(p_{x}^{*}+\beta^{*}E^{*})\\\
p^{\rm{IS}}_{y}=p_{y}^{*}\\\ p^{\rm{IS}}_{z}=p_{z}^{*}\\\
E^{\rm{IS}}=\gamma(E^{*}+\beta^{*}p_{x}^{*})\end{cases}\,,$ (A.8)

where the speed and the Lorentz factor of the resonance can be written as
$\beta^{*}\simeq 1$ and $\gamma=\gamma^{*}$. Considering that the momentum of
the produced particles in a two-body decay is equal to:

$p_{N}^{*}=p_{\pi}^{*}=\frac{\sqrt{s^{2}-2s(m_{\pi}^{2}+m_{N}^{2})+(m_{\pi}^{2}-m_{N}^{2})^{2}}}{2\sqrt{s}}\,,$
(A.9)

we finally get the pion energy in the co-moving reference frame
($\beta_{\pi}=1$):

$E_{\pi}=\Gamma\gamma^{*}(E_{\pi}^{*}+p_{x}^{*})=\frac{E_{p}}{\sqrt{s}}(E_{\pi}^{*}+p^{*}\cos\phi)\,,$
(A.10)

where $p^{*}$ is the common value of $p_{N}^{*}=p_{\pi}^{*}$, and
$0^{\circ}\leq\phi\leq 180^{\circ}$ is the angle, with respect to the proton
direction, at which the pion is emitted from the $\Delta^{+}$ decay in its
center of mass reference system.

###  A.2 Charged Pion Decays

${\pi}^{+}\to\mu^{+}+\nu_{\mu}\,,$ (A.11)
${\pi}^{-}\to\mu^{-}+\bar{\nu}_{\mu}\,.$ (A.12)

We evaluate the kinematics of the $\pi^{+}$ decay, considering that the same
relations will hold in the $\pi^{-}$ decay, which will emerge only in the out
of resonance region (see Sections A.5, A.6). The muon and the neutrino
energies in the $\pi^{+}$ center of mass system (Equation A.11) are:

$E_{\mu^{+}}^{*}=\frac{m_{\pi^{+}}^{2}+m_{\mu^{+}}^{2}}{2m_{\pi^{+}}},$ (A.13)
$E_{\nu_{\mu}}^{*}=\frac{m_{\pi^{+}}^{2}-m_{\mu^{+}}^{2}}{2m_{\pi^{+}}}.$
(A.14)

For the neutrino, $p_{\nu}=E_{\nu}$, so the momentum of the decay products is:

$p_{\nu_{\mu}}^{*}=p_{\mu^{+}}^{*}=\frac{m_{\pi^{+}}^{2}-m_{\mu^{+}}^{2}}{2m_{\pi^{+}}}\,.$
(A.15)

Using Lorentz transformations, we get the energies in the co-moving reference
frame, using the charged pion Lorentz factor
$\gamma_{\pi^{+}}=E_{\pi^{+}}/m_{\pi^{+}}$:

$E_{\mu^{+}}=\gamma_{\pi^{+}}(E_{\mu^{+}}^{*}+p_{\mu^{+}}^{*}\cos\delta),$
(A.16)
$E_{\nu_{\mu}}=\gamma_{\pi^{+}}(E_{\nu_{\mu}}^{*}-p_{\nu_{\mu}}^{*}\cos\delta),$
(A.17)

where $\delta$ is the angle between the muon and the decaying pion,
$0^{\circ}\leq\delta\leq 180^{\circ}$.

###  A.3 Neutral Pion Decay

${\pi}^{0}\to\gamma+\gamma\,.$ (A.18)

In the reference frame in which the pion is at rest, one has:

$p_{\gamma}^{*}=E_{\gamma}^{*}=\frac{1}{2}m_{\pi^{0}}.$ (A.19)

So the energies of the photons resulting from this decay process are:

$E_{\gamma}=\gamma_{\pi^{0}}(E_{\gamma}^{*}\pm
p_{\gamma}^{*}\cos\zeta)=\frac{E_{\pi^{0}}}{2}(1\pm\cos\zeta),$ (A.20)

where $\gamma_{\pi^{0}}=E_{\pi^{0}}/m_{\pi^{0}}$ (and the decay products are
emitted with an angle $0^{\circ}\leq\zeta\leq 180^{\circ}$). The $\pm$ sign
refers to the opposite directions of emission of the two photons.

###  A.4 Muon Decay

$\mu^{-}\to e^{-}+\bar{\nu}_{e}+\nu_{\mu}\,,$ (A.21) $\mu^{+}\to
e^{+}+\nu_{e}+\bar{\nu}_{\mu}\,.$ (A.22)

In the following we refer with $\nu_{\mu}$ and $\nu_{e}$ to indicate muon and
electron neutrinos and antineutrinos. The high-energy neutrino flux reaching
the detectors includes both the $\nu_{\mu}$ produced by the charged pion decay
(see Eq. (A.12)), and the $\nu_{e}$ and $\nu_{\mu}$ produced by the decay of
the muon.

As a consequence of the non-conservation of parity, muons coming from pion
decays are strongly polarized: those coming from the $\pi^{\pm}$ decay have a
right-handed helicity. In the system in which the muon is at rest, the
particle distribution as a function of energy and of the emission angle
$\omega$ is given by:

$\frac{\mathrm{d}N}{\mathrm{d}x\mathrm{d}\omega}=\frac{1}{4\pi}(f_{0}(x)-f_{1}(x)\cos\omega)\,,$
(A.23)

where $x=2E_{l}^{*}/m_{\mu}$, the Bjorken x, is the fraction of the available
energy ($m_{\mu}$) carried by the lepton l (l =
$e^{-},\,\bar{\nu}_{e},\,\nu_{\mu}$), thus: $0\leq x\leq 1$. The variable
$\omega$ defines the angle between the daughter particle and the muon spin.
The functions appearing in the above equation can be obtained through the
Michel parameters, which describe the phase space distribution of leptonic
decays of charged leptons. With these parameters, one can obtain:

  * •

for the muon neutrino

$\begin{cases}f_{0}(x)=2x^{2}(3-2x)\\\ f_{1}(x)=2x^{2}(1-2x)\end{cases}\,;$
(A.24)

  * •

for the electron neutrino

$\begin{cases}f_{0}(x)=12x^{2}(1-x)\\\ f_{1}(x)=12x^{2}(1-2x)\end{cases}\,.$
(A.25)

In the co-moving reference system, Equation (A.23) becomes:

$\frac{dN}{dy}=\frac{1}{\beta_{\mu}}(g_{0}(y,\beta_{\mu})-P_{\mu}g_{1}(y,\beta_{\mu}))\,,$
(A.26)

where

$y\simeq\frac{E_{l}}{E_{\mu}}\,,$ (A.27)

is the Bjorken y, namely the fraction of energy carried by the lepton l in the
co-moving frame, and $P_{\mu}$ is the muon spin projection along the direction
of motion in the co-moving frame. $P_{\mu}$ is defined as ($m_{\pi}$ is the
charged pion mass):

$P_{\mu}=\frac{1}{\beta_{\mu}}\left(\frac{2E_{\pi}(m_{\mu}/m_{\pi})^{2}}{E_{\mu}(1-(m_{\mu}/m_{\pi})^{2})}-\frac{1+(m_{\mu}/m_{\pi})^{2})}{1-(m_{\mu}/m_{\pi})^{2})}\right).$
(A.28)

Considering that in our case $\beta_{\mu}\simeq 1$, we have:

  * •

for the muon neutrino

$\begin{cases}g_{0}(y)=\frac{5}{3}-3y^{2}+\frac{4}{3}y^{3}\\\
g_{1}(y)=\frac{1}{3}-3y^{2}+\frac{8}{3}y^{3}\end{cases}\,;$ (A.29)

  * •

for the electron neutrino

$\begin{cases}g_{0}(y)=2-6y^{2}+4y^{3}\\\
g_{1}(y)=-2+12y-18y^{2}+8y^{3}\end{cases}\,.$ (A.30)

The energy of each neutrino in the co-moving reference frame has been obtained
by multiplying the muon energy defined by Eq. (A.16) with a value y extracted
according to Equation (A.26).

This procedure is valid regardless of the decaying muon’s charge, and it was
used to obtain the neutrino energies coming from the decays shown in Eq.
(A.21) and Eq. (A.22).

###  A.5 Delta Resonance Secondary Peak Region

In the photon energy range
$0.5\,\mathrm{GeV}\leq\epsilon_{r}<2\,\mathrm{GeV}$, the $\Delta^{+}$ is
accompanied by a pion:

$p+\gamma\to\Delta^{+}+\pi^{\mathrm{A}}\,.$ (A.31)

The invariant energy of the system is no longer equal to the mass of the
$\Delta^{+}$ resonance, but it’s given by Eq. (A.1). In the center of mass
frame, the following kinematic relations hold:

$E^{*}_{\Delta^{+}}=\frac{s+m^{2}_{\Delta^{+}}-m^{2}_{\pi^{\rm{A}}}}{2\sqrt{s}}\,,$
(A.32)
$E^{*}_{\pi^{\rm{A}}}=\frac{s+m^{2}_{\pi^{\rm{A}}}-m^{2}_{\Delta^{+}}}{2\sqrt{s}}\,.$
(A.33)

Moreover, $p^{*}_{\pi^{\rm{A}}}$ is given by Eq. A.9 with $\pi=\pi^{\rm{A}}$
and $N=\Delta$. The $\pi^{\rm{A}}$ energy in the co-moving frame can be
written as ($0^{\circ}\leq\eta\leq 180^{\circ}$):

$E_{\pi^{\rm{A}}}=\frac{E_{p}}{\sqrt{s}}(E^{*}_{\pi^{\rm{A}}}+p^{*}_{\pi^{\rm{A}}}\cos\eta)\,.$
(A.34)

The charged and neutral pions are produced with the following branching
ratios:

  * •

$1/15$ of the $A$ pions are neutral, $\pi^{0}$, which decay in two photons
(Eq. (A.18)) with energies given by Eq. (A.20);

  * •

$8/15$ of the $A$ pions are positive charged, $\pi^{+}$, which decay in a muon
neutrino and a positive muon (Eq. (A.11), which in turn decays in a positron,
an electron neutrino and a muon antineutrino, Eq. (A.22)). Energies for these
particles are given by Eqs. (A.16), (A.17), and (A.27);

  * •

$2/5$ of the $A$ pions are negative charged, $\pi^{-}$, which decay in a muon
antineutrino and a muon (Eq. (A.12), which in turn decays in an electron, an
electron antineutrino and a muon neutrino, Eq. (A.21)). Energies for these
particles are also given by Eqs. (A.16), (A.17), and (A.27).

In order to consider all the processes that can originate neutrinos and
photons, we have to take into account the $\Delta^{+}$ decay:

$\Delta^{+}\to\mathcal{N}+\pi^{\mathrm{B}}\,.$ (A.35)

In the $\Delta^{+}$ resonance rest frame, we have:

$|p^{*}_{\pi^{\mathrm{B}}}|=|p^{*}_{\mathcal{N}}|=\frac{\sqrt{[(m^{2}_{\Delta^{+}}-(m_{\mathcal{N}}+m_{\pi^{\mathrm{B}}})^{2})(m^{2}_{\Delta^{+}}-(m_{\mathcal{N}}-m_{\pi^{\mathrm{B}}})^{2})]}}{2m_{\Delta^{+}}}\,,$
(A.36)
$E^{*}_{\pi^{\mathrm{B}}}=\frac{s+m^{2}_{\pi^{\rm{B}}}-m^{2}_{\mathcal{N}}}{2m_{\Delta^{+}}}\,.$
(A.37)

Finally, the $B$ pion energy in the co-moving frame is given by:

$E_{\pi^{\mathrm{B}}}=\frac{E_{\Delta^{+}}}{m_{\Delta^{+}}}(E_{\pi^{\mathrm{B}}}^{*}+p^{*}_{\pi^{\mathrm{B}}}\cos\rho)\,,$
(A.38)

where $0^{\circ}\leq\rho\leq 180^{\circ}$ is the emitting angle of the pion,
and $E_{\Delta^{+}}\sim E_{p}-E_{\pi^{\mathrm{A}}}$. The code extracts
randomly the branching ratios for $\pi^{\mathrm{B}}$, which are:

  * •

$18/45$ of $B$ pions are neutral pions;

  * •

$19/45$ of $B$ pions are positive charged pions;

  * •

$8/45$ of $B$ pions are negative charged pions.

The energies of the particles emerging from the pion decay are obtained as
already described in the above sections.

###  A.6 Multipion Production Region

When the photon energy exceeds 2 GeV, most of the energy lost by the proton
($\sim 0.6\,E_{p}$) is split equally among three pions, and the neutral and
charged pions are approximately produced in equal numbers (this is a
simplified treatment, as found in e.g. [14, 15, 19]). So the single pion
energy is equal to $E_{\pi}=0.2\,E_{p}$, and the particles resulting from the
pion decay are simulated according with the procedures described in the
previous sections.

## References

  * [1] M. G. Aartsen et al. [IceCube Collaboration], Evidence for high-energy extraterrestrial neutrinos at the icecube detector, Science 342 (2013) 1242856 [arXiv:1311.5238].
  * [2] M. G. Aartsen et al. [IceCube Collaboration], Time-integrated Neutrino Source Searches with 10 years of IceCube Data, Phys. Rev. Lett. 124 (2020) 051103 [arXiv:1910.08488].
  * [3] S. Celli, A. Palladino and F. Vissani, Neutrinos and $\gamma$-rays from the Galactic Center Region After H.E.S.S. Multi-TeV Measurements, EPJC 77 (2017) 66 [arXiv:1604.08791].
  * [4] M. Ackermann et al. [Fermi-LAT Collaboration], The spectrum of isotropic diffuse gamma-ray emission between 100 MeV and 820 GeV, ApJ 799 (2015) 86 [arXiv:1410.3696].
  * [5] K. Murase, D. Guetta and M. Ahlers, Hidden Cosmic-Ray Accelerators as an Origin of TeV-PeV Cosmic Neutrinos, Phys. Rev. Lett. 116 (2016) 071101 [arXiv:1509.00805].
  * [6] P. Meszaros and E. Waxman, TeV neutrinos from successful and choked Gamma-Ray Bursts, Phys. Rev. Lett. 87 (2001) 171102 [arXiv:astro-ph/0103275v3].
  * [7] T. Piran, E. Nakar, P. Mazzali and E. Pian, Relativistic Jets in Core Collapse Supernovae, ApJL 871 (2019) L5 [arXiv:1704.08298].
  * [8] L. Izzo et al., Signatures of a jet cocoon in early spectra of a supernova associated with a $\gamma$-ray burst, Nature 565 (2019) 324 [arXiv:1901.05500].
  * [9] Y. Q. Xue et al., A magnetar-powered X-ray transient as the aftermath of a binary neutron-star merger, Nature 568 (2019) 198 [arXiv:1904.05368].
  * [10] P. B. Denton and I. Tamborra, Exploring the properties of choked Gamma-ray Bursts with IceCube’s high-energy neutrinos, ApJ 855 (2018) 37 [arXiv:1711.00470].
  * [11] H. He, A. Kusenko, S. Nagataki, Y. Fan and D. Wei, Neutrinos from Choked Jets Accompanied by Type-II Supernovae, ApJ 856 (2018) 119 [arXiv:1803.07478].
  * [12] A. Esmaili and K. Murase, Constraining high-energy neutrinos from choked-jet supernovae with IceCube high-energy starting events, JCAP 12 (2018) 008 [arXiv:1809.09610].
  * [13] D. Guetta, R. Rahin, I. Bartos and M. Della Valle, Constraining the fraction of core-collapse supernovae harboring choked jets with high-energy neutrinos, MNRAS 492 (2020) 843 [arXiv:1906.07399].
  * [14] A. M. Atoyan and C. D. Dermer, Neutral beams from blazar jets, ApJ 586 (2003) 79 [arXiv:astro-ph/0209231v2].
  * [15] S. R. Kelner and F. A. Aharonian, Energy spectra of gamma-rays, electrons and neutrinos produced at interactions of relativistic protons with low energy radiation, Phys. Rev. D 78 (2008) 034013 [arXiv:0803.0688].
  * [16] F. Lucarelli, Master Thesis at Sapienza University of Rome, Rivelazione di neutrini di alta energia da sorgenti puntiformi extragalattiche con l’apparato Cherenkov sottomarino NESTOR (1998) https://www.roma1.infn.it/people/capone/AHEN/Thesis/FabrizioLucarelli/tesilaureaflucarelli.pdf
  * [17] K. Murase and S. Nagataki, High energy neutrino emission and neutrino background from gamma-ray bursts in the internal shock model, PRD, 73 (2006) 063002.
  * [18] K. Murase, Prompt high-energy neutrinos from gamma-ray bursts in photospheric and synchrotron self-Compton scenarios, PRD 78 (2008) 101302.
  * [19] S. Hümmer, M. Rüger, F. Spanier and W. Winter, Simplified models for photohadronic interactions in cosmic accelerators, ApJ 721 (2010) 630 [arXiv:1002.1310v3].
  * [20] K. Murase and K. Ioka, TeV-PeV Neutrinos from Low-Power Gamma-Ray Burst Jets inside Stars, PRL 111 (2013) 121102.
  * [21] N. Senno, K. Murase and P. Meszaros, Choked jets and low-luminosity gamma-ray bursts as hidden neutrino sources, PRD 93 (2016) 083003.
  * [22] A. I. MacFadyen, S. E. Woosley and A. Heger, Supernovae, Jets, and Collapsars, ApJ 550 (2001) 410 [arXiv:astro-ph/9910034].
  * [23] M. A. Aloy et al., Relativistic Jets from Collapsars, ApJ 531 (2000) 119.
  * [24] S. E. Woosley, R.G. Eastman and B. P. Schmidt, Gamma-Ray Bursts and Type IC Supernova SN 1998BW, ApJ 516 (1999) 2.
  * [25] S. Razzaque, P. Meszaros and E. Waxman, TeV neutrinos from core collapse supernovae and hypernovae, Phys. Rev. Lett. 93 (2004) 181101 [arXiv:astro-ph/0407064v4].
  * [26] E. Sobacchi and Y. E. Lyubarsky, Instability induced by recollimation in highly magnetized outflows, MNRAS 480 (2017) 4948S.
  * [27] D. Xiao and Z. G. Dai, Neutrino Emission in the Jet Propagation Process, APJ 790 (2014) 59X.
  * [28] P. Kumar, R. Narayan and J. L. Johnson, Properties of Gamma-Ray Burst Progenitor Stars, Science 321 (2008) 376K.
  * [29] H. T. Janka, T. Melson and A. Summa, Physics of Core-Collapse Supernovae in Three Dimensions: A Sneak Preview, Annual Review of Nuclear and Particle Science 66 (2016) 341J.
  * [30] O. Bromberg, E. Nakar, T. Piran and R. Sari, The Propagation of Relativistic Jets in External Media, ApJ 740 (2011) 100B.
  * [31] A. Mizuta and K. Ioka Opening Angles of Collapsar Jets, ApJ 777 (2013) 162.
  * [32] E. Fermi, On the Origin of the Cosmic Radiation, Phys. Rev. 75 (1949) 1169.
  * [33] A. Atoyan and C. D. Dermer, High-Energy Neutrinos from Photomeson Processes in Blazars, PRL 87 (2001) 1102A.
  * [34] L. Morejon, A. Fedynitch, D. Boncioli, D. Biehl and W. Winter, Improved photomeson model for interactions of cosmic ray nuclei, JCAP 11 (2019) 007 [arXiv:1904.07999].
  * [35] P. Baerwald, S. Hümmer, W. Winter, Systematics in the interpretation of aggregated neutrino flux limits and flavor ratios from gamma-ray bursts, Astroparticle Physics 35 (2012) 508-529 [arXiv:1107.5583].
  * [36] K. Murase, High energy neutrino early afterglows from gamma-ray bursts revisited, Phys. Rev. D 76 (2007) 123001 [arXiv:0707.1140v3].
  * [37] F. L. Villante and F. Vissani, How precisely neutrino emission from supernova remnants can be constrained by gamma ray observations?, Phys. Rev. D 78 (2008) 103007 [arXiv:0807.4151].
  * [38] C. Mascaretti and F. Vissani, On the relevance of prompt neutrinos for the interpretation of the IceCube signals, JCAP 08 (2019) 004M [arXiv:1904.11938].
  * [39] S. Adrian-Martìnez et al. [ANTARES Collaboration], Search for cosmic neutrino point sources with four years of data of data from ANTARES telescope, ApJ 760 (2012) 53.
  * [40] S. Adrian-Martìnez et al. [KM3NeT Collaboration], KM3NeT 2.0 Letter of Intent for ARCA and ORCA, J. Phys. G: Nucl. Part. Phys. 43 (2016) 084001.
  * [41] M. G. Aartsen et al. [IceCube Collaboration], Searches for Extended and Point-like Neutrino Sources with Four Years of IceCube Data, ApJ 796 (2014) 2.
  * [42] N. M. Lloyd-Ronning, A. Aykutalp and J. L. Johnson, On the cosmological evolution of long gamma-ray burst properties, MNRAS 488 (2019) 4.
  * [43] V. Petrosian, E. Kitanidis and D. Kocevski, Cosmological Evolution of long Gamma-Ray Bursts and the Star Formation Rate, ApJ 806 (2015) 44.
  * [44] P. Madau and M. Dickinson Cosmic Star-Formation History, ARAA 52 (2014) 415 [arXiv:1403.0007].
  * [45] R. Abbasi et al. [IceCube Collaboration], The IceCube high-energy starting event sample: Description and flux characterization with 7.5 years of data, [arXiv:2011.03545v1].
  * [46] J. Stettner, Measurement of the Diffuse Astrophysical Muon-Neutrino Spectrum with Ten Years of IceCube Data, Pos 1017 (2019).
  * [47] M. G. Aartsen et al. [IceCube Collaboration], Constraints on Minute-Scale Transient Astrophysical Neutrino Sources, PRL 122 (2019) 051102.
  * [48] J. Bahcall and E. Waxman, High Energy Astrophysical Neutrinos: the Upper Bound is Robust, Phys. Rev. D 64 (2001) 023002 [arXiv:hep-ph/9902383].
  * [49] M. G. Aartsen et al. [IceCube Collaboration], An all sky search for three flavors of neutrinos from gamma-ray bursts with the IceCube neutrino observatory, The Astrophysical Journal 824 (2016).
  * [50] A. Albert et al. [ANTARES Collaboration], Constraining the contribution of Gamma-Ray Bursts to the diffuse neutrino flux with the ANTARES dataset (2007-2017), MNRAS 500 (2021) 5614.
  * [51] S. Adrian-Martìnez et al. (ANTARES Collaboration), Search for high-energy neutrinos from bright GRBs with ANTARES, MNRAS 469 (2017) 1.
  * [52] K. Murase and E. Waxman, Constraining high-energy cosmic neutrino sources: Implications and prospects, PRD 94 (2016) 103006.
  * [53] N. Senno, K. Murase and P. Meszaros, Constraining high-energy neutrino emission from choked jets in stripped-envelope supernovae, JCAP 01 (2018) 025.
  * [54] C. Irwin, E. Nakar and T. Piran, The propagation of choked jet outflows in power-law external media, MNRAS 489 (2019) 2844.
  * [55] R. Perna, D. Lazzati and M. Caniello, Electromagnetic Signatures of Relativistic Explosions in AGN Disks, ApJL 906 (2021) 7P.
  * [56] J. P. Zhu et al., Neutron Star Mergers in AGN Accretion Disks: Cocoon and Ejecta Shock Breakouts, ApJL 906L (2021) 11Z.