# The role of jet-cocoon mixing, magnetization and shock breakout in neutrino
and cosmic-ray emission from short GRBs

Ore Gottlieb11affiliation: School of Physics & Astronomy, Tel Aviv University,
Tel Aviv 69978, Israel. E-mail<EMAIL_ADDRESS>22affiliation: Center for
Interdisciplinary Exploration & Research in Astrophysics (CIERA), Physics &
Astronomy, Northwestern University, Evanston, IL 60208, USA and Noémie
Globus33affiliation: ELI Beamlines, Institute of Physics, Czech Academy of
Sciences, 25241 Dolní Břežany, Czech Republic. E-mail<EMAIL_ADDRESS>44affiliation: Center for Computational Astrophysics, Flatiron Institute,
Simons Foundation, New-York, NY10003, USA 55affiliation: Astrophysical Big
Bang Laboratory (ABBL), Cluster for Pioneering Research, RIKEN, 2‐1 Hirosawa,
Wakō, Saitama 351‐0198, Japan

###### Abstract

We perform GRMHD and RMHD simulations of weakly and highly magnetized gamma-
ray burst (GRB) jets propagating in binary neutron star (BNS) merger ejecta.
Using the simulations, we first find that mixing between the jet and cocoon,
which is present in all types of jets, inhibits the formation of
subphotospheric collisionless shocks. However, we show that a mild
magnetization may lead to the formation of collisionless subshocks which allow
efficient proton acceleration. We consider shear acceleration and diffusive
shock acceleration at collimation shocks, internal shocks, shock breakout, and
external shocks, to provide the first estimate for neutrino and cosmic-ray
(CR) signals from self-consistent simulations of GRBs in BNS mergers. We find
that short GRBs do not produce detectable neutrino signals with current-day
facilities. Shock breakout yields $\sim 10$ PeV neutrinos at viewing angles
$\sim 20^{\circ}$, independent of the jet magnetization. However, a neutrino
signal from shock breakout is well below detection limits of current
detectors. Such a signal would allow a coincident neutrino-$\gamma-$ray
detection, providing a testable prediction for shock breakout as a neutrino
production site. Using the numerical modeling that fits GW170817 afterglow
emission, we find that blast waves in BNS mergers can account for 5%-10% of
the Galactic CR luminosity in the PeV-EeV energy range. Based on these
estimates, the observed level of CR anisotropy places a constraint on the
distance of the latest Galactic binary neutron star merger to $\lesssim 3$
kiloparsecs.

cosmic rays, neutrinos, gamma-ray bursts, neutron star mergers

##  1 Introduction

The year 2017 was a breakthrough in multi-messenger astronomy, with the
detection of the binary neutron star (BNS) merger GW170817 through
gravitational waves (GW) in advanced LIGO/Virgo (Abbott et al., 2017) and
electromagnetic emission (see e.g., Nakar, 2020; Margutti & Chornock, 2020,
for review). However, although $\gamma$-ray bursts (GRB) are considered to be
a promising source of neutrinos (e.g., Waxman & Bahcall, 1997), no neutrino
associated with the event has been detected by IceCube (Albert et al., 2017),
presumably due to the large viewing angle with respect to the jet axis,
$\,{\theta_{\rm obs}}={20^{\circ}}^{+8^{\circ}}_{-7^{\circ}}$ (Mooley et al.,
2018). Future neutrino detectors will enable us to connect neutrino and
electromagnetic signals. With the majority of BNS mergers with a detectable
signal of jets expected at viewing angles outside the jet opening angle
(Gottlieb et al., 2019), understanding the dependency of the neutrino flux
upon the viewing angle is of utmost importance.

The GW170817’s multiwavelength campaign showed evidence for a structured jet,
as the origin of the $\gamma$-ray emission and the first months of X-ray and
radio emission originated in a mildly relativistic ($\Gamma\sim 3$) cocoon
(e.g., Gottlieb et al., 2018). The cocoon, which is inflated during the
interplay between the jet and the merger’s ejecta, applies pressure on the jet
and could form a collimation shock at the jet’s base. The jet may successfully
break out of the ejecta or be choked inside. Regardless of the jet’s fate, the
cocoon eventually breaks out of the ejecta, radiating its internal energy upon
breakout, in what is likely to be the origin of GRB170817A (Gottlieb et al.,
2018). The question then arises if there is a specific signature in energetic
particles associated with the fate of the jet that is connected to the nature
of the central engine, which is still under debate. One possible jet
production mechanism is the magnetic extraction of the spin energy of a Kerr
black hole (BZ, Blandford & Znajek, 1977). If the inner regions of the disk
are dominated by neutrino cooling, the large amount of MeV neutrinos
annihilate in the polar region to produce a hot, pressure-driven ${\rm
e}^{+}{\rm e}^{-}$ outflow that may interfere with the magnetic energy
extraction (Globus & Levinson, 2014). Therefore, the jet could be either
weakly or highly magnetized, depending on the plasma injection processes near
the black hole (BH) horizon. As the launching mechanism is still under debate,
unmagnetized, mildly and highly magnetized jets should all be investigated for
neutrino production.

Previous studies considered neutrino emission in GRB jets at their collimation
and internal shocks (e.g., Murase & Ioka, 2013; Globus et al., 2015; Biehl et
al., 2018; Kimura et al., 2018a), and shear turbulence between the jet and the
cocoon (e.g., Rieger & Duffy, 2004; Sahayanathan, 2009; Rieger & Duffy, 2016).
Additional studies also examined neutrino emission from sources outside the
jet such as supernova shock breakout (e.g., Waxman & Loeb, 2001; Katz et al.,
2011; Kashiyama et al., 2013) and cocoons of successful or choked jets (e.g.,
Ostrowski, 2000; Senno et al., 2016; He et al., 2018). These studies estimated
analytically the structure of the shock, energy reservoir and magnetic
dissipation at the aforementioned production sites. However, the jet dynamics
involves a non-linear evolution of the complex jet-cocoon system which makes
relativistic magnetohydrodynamic simulations an essential tool for studying
the jet evolution in BNS mergers (e.g., Christie et al., 2019; Fernández et
al., 2019; Kathirgamaraju et al., 2019). One particularly important property
that emerges in numerical simulations is the mixing between the jet and the
cocoon. In unmagnetized and weakly magnetized jets, the mixing comes about
following the development of hydrodynamic instabilities on the jet-cocoon
boundary, loading the jet with baryons and reducing its Lorentz factor (e.g.,
Gottlieb et al., 2021a, b; Matsumoto et al., 2021). The mixing, which has been
ignored in previous works, loads the jet with baryons, and thus inhibits the
formation of collisionless shocks which are necessary for neutrino production.

BZ-launched GRB jets have a large degree of magnetization which has a
tremendous effect on the neutrino production. First, the magnetization
significantly alters the jet structure such that the energy is distributed
differently between the jet and the cocoon (Gottlieb et al., 2020a). Second,
and more importantly, the magnetization needs to be neither too high nor too
low for neutrino production below the photosphere. On one hand, the jet needs
to be sufficiently magnetized to support the formation of subphotospheric
collisionless shocks. Upon launching, the jet has to initially punch through
the optically-thick ejecta, forming radiation-mediated shocks. In such shocks
the scattering mean free path is much larger than the Larmor radius implying
that such shocks are inefficient for particle acceleration (see e.g., Levinson
& Bromberg, 2008). If the jet is magnetized, strong subshocks may emerge
within the radiation-mediated shock, on a Larmor radius scale (Beloborodov,
2017, hereafter B17). On the other hand, in the presence of very strong
magnetic fields, pions, which constitute the main neutrino production channel,
cool down by synchrotron much faster than their half-life time, thereby
significantly mitigating the neutrino production (see e.g., Lipari et al.,
2007). Therefore, the diffusive shock acceleration (DSA) efficiency and hence,
the neutrino outcome are strongly dependent upon the evolution of the
magnetization, which can be studied with numerical simulations.

CR production from both mildly and ultra relativistic DSA has been studied in
GRB jets following their interaction with the surrounding medium in the
afterglow phase (Vietri, 1995; Dermer & Atoyan, 2006). In the context of BNS
mergers, it has been shown that Galactic BNS (GBNS) remnants may significantly
contribute to the observed CR in the energy range between the knee and the
ankle, the so called ”shin” region (Kimura et al., 2018b). However, despite
the valuable information about the blast wave structure and energetics that
was provided by GW170817 afterglow, the contribution of the jet-cocoon blast-
wave to CR remained unexplored.

In this Letter, we investigate particle production in short GRBs (sGRBs) by
performing 3D relativistic magneto-hydrodynamic (RMHD) and general
relativistic magneto-hydrodynamic (GRMHD) simulations of structured jets
propagating in a BNS ejecta, using GW170817 observations for calibrating the
energetics of the jet-cocoon system. We consider neutrino production by DSA in
jets with different magnetizations in regions where strong collisionless
subshocks can form. That is: collimation shocks, internal shocks and breakout
of the forward shock from the ejecta. We also examine shear acceleration
between the jet and the cocoon. Finally, we consider DSA at the external
relativistic blast-wave to estimate CR production.

In §2, we show that previous criteria for the formation of collisionless
shocks in the jet are not satisfied in numerical simulations, and magnetic
fields are necessary for the formation of collisionless subshocks. In §3, we
briefly describe the setup of our simulations and where collisionless shocks
may form. In §4 and §5 we present the first neutrino and CR calculations that
are based on realistic sGRB simulations, estimated at a variety of production
sites and for different types of jets. In §6 and §7 we summarize and discuss
the implications of our results.

##  2 Criterion for efficient CR acceleration in radiation-mediated shocks

Subphotospheric shocks are leading candidates for neutrino production by
virtue of Fermi acceleration. These shocks are mediated by photons in the
optically thick ejecta which supports $p\gamma$ interactions between the
trapped photons and the accelerated protons. However, efficient Fermi
acceleration can only take place in collisionless shocks where particle
interactions are mediated by plasma processes, raising the question of whether
subphotospheric shocks can efficiently accelerate CR. Murase & Ioka (2013)
suggested that a collisionless shock can form inside the optically thick
region if the shock’s comoving upstream medium is optically thin, namely
$n\sigma_{T}R^{\prime}\lesssim 1$, where $n$ is the proper number density,
$\sigma_{T}$ is the Thomson cross section, and $R^{\prime}$ is the proper size
of the shock. This criterion can be applied to both oblique collimation shocks
and internal shocks which reside inside the optically thick medium. According
to this criterion, Murase & Ioka (2013) found that neutrino production is
disfavored in high-luminosity choked jets, leaving low-luminosity jets as
potential neutrino sources. However, due to the mixing between the jet and the
cocoon material, which is present in all types of jets, we find that the above
inequality is not satisfied in any type of sGRB jets, whether it is high or
low luminosity (Gottlieb et al., 2021b), choked or successful. Thus,
collisionless shocks and a neutrino signal are unlikely to emerge under this
criterion in sGRBs, as we discuss next.

For typical parameters, the criterion reads that an unmagnetized collimation
shock forms a collisionless shock if $L_{{\rm
iso},48}R_{s,10}\Gamma_{2}^{-3}\lesssim 1$, where $Q_{x}$ denotes the value of
the quantity $Q$ in units of $10^{x}$ times its c.g.s. units. The weakest
observed sGRB jets maintain $\gamma$-ray luminosity of approximately $L_{{\rm
iso},48}\approx 10$ (e.g., Shahmoradi & Nemiroff, 2015). That implies that in
order to obtain collisionless collimation shocks according to the above
criterion, either i) the collimation shock reaches $R_{s,10}\gg 1$, which is
unlikely for typical sGRB durations, particularly in low luminosity jets that
likely propagate at subrelativistic velocities inside the ejecta; or ii) the
Lorentz factor is of a few hundreds early on, but the mixing between the jet
and the dense ejecta (which is necessary for collimation) inhibits such
relativistic velocities at small radii. Thus, the hydrodynamic evolution of
sGRB jets seems to rule out the possibility of having neutrino emission from
collimation shocks of hydrodynamic jets.

We also find that strong internal shocks, which seem to necessitate a variable
engine, do not support the formation of collisionless shocks according to the
aforementioned criterion. The reason lies in recent results of Gottlieb et al.
(2020b) who found intense mixing from the cocoon into the jet during the low
power episodes of the jet. The mixing reduces the Lorentz factor and
substantially increases the baryon loading in the jet such that the comoving
optical depth increases dramatically, for all jets in general, and due to more
intense mixing, in low luminosity jets in particular.

If instead, the jet is magnetized, B17 showed that under certain conditions,
strong subphotospheric collisionless subshocks may form within the radiation-
mediated shock on a Larmor radius scale. In particular, the collimation and
internal shocks that occur below the photosphere may develop such
collisionless subshocks and efficiently accelerate CR in the jet. B17 found
that $\chi\equiv 2\beta_{p}$ dictates the shock dissipation efficiency as
follows ($\beta_{p}$ is the ratio between the thermal to magnetic pressure,
$p_{t}/p_{m}$). For mildly relativistic shocks, if the upstream magnetization
is $0.1\lesssim\sigma_{u}\equiv\frac{B^{\prime 2}}{4\pi\rho c^{2}}\lesssim 1$
($B^{\prime}$ being the plasma proper magnetic field, $\rho$ is the mass
density) and $\chi\lesssim 2$, the photons do not have sufficient energy to
fully decelerate the upstream to the downstream velocity, subsequently a
collisionless subshock forms to satisfy the jump conditions at the shock. In
the case of $\chi\gg 1$, the magnetization is too weak to sustain a strong
subshock, and for $\chi\ll 1$ the magnetization is too high, and the kinetic
energy dissipation efficiency is low. This criterion implies that
collisionless shocks can emerge during the jet evolution inside the optically-
thick ejecta if the jet magnetization is $\chi\sim 1$, and $\sigma_{u}\gtrsim
0.1$. The latter also promises that the collisionless subshock dominates the
energy in the mildly relativistic shock (see e.g., figure 11 in Levinson &
Nakar, 2019).

At highly magnetized, relativistic shocks, Fermi acceleration is known to be
inefficient. Using PIC simulations, Sironi et al. (2013) found a critical
magnetization of $\sigma\gtrsim 1\%$, above which the particles cannot travel
back upstream even when moving at $\sim c$. They showed that at the critical
magnetization efficient acceleration can take place only if the upstream
Lorentz factor is $\Gamma_{u}\lesssim 5$. We thus consider only mildly-
magnetized and mildly relativistic shocks (however further studies of CR
acceleration efficiency in this regime are needed, see e.g., Crumley et al.,
2019). Assuming the formation of strong subshocks in this regime, we consider
a kinetic energy dissipation efficiency of $\epsilon_{d}=0.2$ (e.g., Mimica &
Aloy, 2010; Komissarov, 2012).

##  3 Numerical simulations overview

We calculate the jet structure and the resulting CR and neutrino emission
during different phases of the jet evolution, up to the afterglow phase. We
consider a variety of jet parameters, from hydrodynamic to highly magnetized
jets, as well as jets which successfully break out from the ejecta and those
that are choked inside. In all simulations, the initial setup includes a
homologous, expanding, subrelativistic ($v<0.2$c) cold ejecta that emerges in
the aftermath of the merger. The ejecta mass is chosen to be $M\approx
0.05\,{\rm M_{\odot}}$, as inferred from GW170817 observations, and we choose
its radial mass density profile to be $\rho(r)\propto r^{-2}$. In the
hydrodynamic models the dense ejecta is embedded in a lighter mildly
relativistic homologous tail ejecta that moves at $v<0.6$c.

Our simulations suggest that the early jet evolution can be largely divided
into two types of jets: highly magnetized and unmagnetized/weakly magnetized.
The setups of both cases is described below. For the reader interested in the
technicalities of the simulations, more details are provided in Appendix A.

Figure 1: Magnetized jets in expanding ejecta. Panels (a-f): Meridional maps
of a highly magnetized successful jet (model $\,{\it HM}$) inside the ejecta
0.4 s after the jet launch (0.9 s after the merger). Shown are: Logarithmic of
the comoving magnetic field (a), magnetization (b), (non logarithmic) Lorentz
factor (c), kinetic energy density (d), mass density (e), and $\chi\equiv
2p_{t}/p_{m}$ (f). Black circles mark examples of internal shocks. Panels
(g-h): The evolution of the luminosity (g) and the maximal terminal Lorentz
factor (product of magnetization and Lorentz factor upon launching,
$\sigma_{0}\Gamma_{0}$) (h) as measured on the horizon, for the case of highly
magnetized jet. Dashed red line in the luminosity (g) panel reflects the
average luminosity. Panel (i-n): same as panels (a-f) (but Lorentz factor is
logarithmic) for a weakly magnetized jet (model $\,{\it WM}$) 1s after the jet
injection (1.6s after the merger), when the jet already broke out from the
ejecta. White ellipses mark the collimation shock. In the panels of $\chi$ we
only present the jet and the cocoon, since $p_{t}=0$ in the unshocked
surrounding medium.

(i) Highly magnetized (model $\,{\it HM}$) jet (performed with the 3D GPU-
accelerated general relativistic magneto-hydrodynamic (GRMHD) code h-amr
(Liska et al., 2019) in Kerr-Schild coordinates, using adaptive mesh
refinement (AMR) and local adaptive time-stepping). The jet launching is
affected by the accreted mass, and thus in the GRMHD simulation, where
accretion is included, we also consider a more realistic angular
distribution111We also apply this angular dependence in model $\,{\it CU}$.:
$\rho(r,\theta)\propto r^{-2}(\frac{1}{4}+{\rm sin}^{3}\theta)$ such that most
of the mass lies close to the accretion disk on the equatorial plane222Note
that our ejecta is a single component with angular and radial profiles, and
the accretion disk forms by the angular momentum of the system. This is in
contrast to previous works (e.g., Christie et al., 2019; Fernández et al.,
2019) which included the disk as part of their initial conditions.. The rest
of the initial conditions include a rapidly spinning (dimensionless spin
$a=0.8$) Kerr BH with mass $\,{M_{{\rm BH}}}=3\,{\rm M_{\odot}}$, embedded in
the aforementioned expanding ejecta. The ejecta radial density drops linearly
from 10 gravitational radii, $10r_{g}=10\,{\cal G}\,{M_{{\rm BH}}}c^{-2}$
towards $r_{g}$. The disk angular velocity is $10^{-3}~{}{\rm s^{-1}}$ for
$r/r_{g}<70$, and $5(r_{g}/r)^{2}~{}{\rm s^{-1}}$ otherwise. The initial
vector-potential profile is a dipole-like with a core of radius $\sim
100r_{g}$, and maximum magnetic field of $10^{13}$ G.

The jet is launched self-consistently following the accretion of a magnetized
ejecta onto a central BH, and maintains a luminosity of $L\approx 10^{50}{\rm
erg~{}s^{-1}}$, similar to that in GW170817. The jet is initially collimated
by the winds from the disk, but most of the collimation process is done by the
cocoon as can be seen in http://www.astro.tau.ac.il/~ore/NSM_GRMHD.html). We
follow the jet propagation up to $5\times 10^{9}$ cm, 0.4 second after its
launch. Fig. 1 depicts maps of the jet structure (panels a-f) and its
luminosity and maximal terminal Lorentz factor on the event horizon (panels
g-h). Movies of the jet evolution can be found in
http://www.astro.tau.ac.il/~ore/NSM_GRMHD.html. The simulation features a
small, highly magnetized collimation shock/nozzle that forms at the base of
the jet. The strong magnetic field at the collimation shock of
$B^{\prime}\approx 10^{12}$ G (panel b) leads to pion and kaon cooling much
faster than their decay times, and thus does not support efficient neutrino
production.

Our focus in model $\,{\it HM}$ is on neutrino production from internal shocks
(black circles in Fig. 1a,d). We find that within the typical GRB engine time
of $\lesssim 0.5$ s, internal shocks with the required upstream magnetization
$\sigma_{u}\gtrsim 0.1$ and $\chi\lesssim 2$ (Fig. 1f) primarily take place at
$r\lesssim 3\times 10^{9}$ cm where $\Gamma\lesssim 3$, thereby being mildly
relativistic, and mildly magnetized shocks with $B^{\prime}\approx 10^{11}$ G.
At later times and larger radii, if the jet engine work-time is sufficiently
long, the jet accelerates to ultra-relativistic velocities ($\Gamma\gtrsim
100$). At larger radii, the magnetization drops, and thus the shocks become
ultra relativistic and weakly-magnetized such that the formation of
collisionless shocks is suppressed, and thus we do not consider this phase in
our calculation.

(ii) Hydrodynamic (models $\,{\it SU1}$, $\,{\it SU2}$, $\,{\it SU3}$) and
weakly magnetized (model $\,{\it WM}$) jets (performed with the 3D
relativistic magneto-hydrodynamic RMHD simulations code pluto (Mignone et al.,
2007) with constrained transport scheme that keeps $\nabla\cdot B=0$). We
choose the jet parameters that are motivated by sGRB observations. The jets
are injected with maximal initial magnetization $\sigma_{0}$ at $t=t_{d}$ for
a total jet duration of $t_{j}$, into the expanding ejecta with a Lorentz
factor $\Gamma_{0}\approx 0.7\theta_{j}^{-1}$ (Mizuta & Ioka, 2013; Harrison
et al., 2018). The jet initially accelerates under its thermal pressure to
relativistic velocities within an opening angle of $\theta_{j}$. As the fluid
elements encounter the cocoon, the jet becomes collimated by forming a
collimation shock of a size $R_{s}$. Following the collimation shock crossing,
the jet elements decelerate to mildly relativistic Lorentz factor
$\Gamma\sim\theta_{j}^{-1}$ (e.g., Bromberg et al., 2011; Gottlieb et al.,
2019). The emerging weakly magnetized ($\,{\it WM}$) jet-cocoon structure is
shown in panels (i-n) in Fig. 1.

Previous numerical simulations of weakly magnetized jets in a dense medium
(Gottlieb et al., 2020a) suggested that a substantial magnetic dissipation
takes place at the collimation shocks, such that the jet magnetization drops
to values of $\sigma\lesssim 0.1$ following the shocks, disfavoring particle
acceleration at internal shocks later on333Note however that the jet
magnetization could be amplified at the collimation and internal shocks to
$\sigma\gtrsim 0.1$ (Gottlieb et al., 2021a).. When considering the neutrino
production at the collimation shock from our weakly magnetized jets (model
$\,{\it WM}$, white ellipses in Fig. 1k,m), we find that the upstream Lorentz
factor is marginally mildly relativistic, but for initial magnetization
$\sigma_{0}=10^{-2}$ (model $\,{\it WM}$) the magnetization is too low and
$\chi\gg 1$ at all potential sites such that no neutrino production is
expected prior to breakout. A mild magnetization of $\sigma_{0}\approx 0.1$
may allow $\chi\gtrsim 1$ at the collimation shock to support neutrino
production, however due to numerical limitations we could not perform such a
simulation. Instead, we assume in the next section that the collimation shock
properties are not considerably different between $\sigma_{0}=10^{-2}$ and
$\sigma_{0}=10^{-1}$, and consider neutrino production in weakly magnetized
jets from their collimation shocks based on our $\,{\it WM}$ model.

As the jet-cocoon outflow breaks out of the ejecta the forward shock becomes
collisionless, and the magnetization drops substantially such that the outflow
can be approximated as hydrodynamic444However, while the jets are likely to be
unmagnetized at late times, the initial magnetization does affect the emerging
jet-cocoon structure (Gottlieb et al., 2020a).. Thus, for the breakout
emission and the afterglow phase, we consider successful and choked
hydrodynamic jets. We utilize the fact that the structure of the jet at these
stages can be inferred from observations. For the shock breakout emission we
average the angular energy and velocity distributions of two models ($\,{\it
SU2}$ and $\,{\it SU3}$), and calibrate their total observed shock breakout
energy at $20^{\circ}$ by GRB170817A observations assuming the energy source
was shock breakout (see Appendix B). For the blast wave emission, we use the
jet-cocoon model that fits GW170817 afterglow light curve from Mooley et al.
(2018) (model $\,{\it SU1}$), and use an additional uncollimated choked jet
(model $\,{\it CU}$). The jet parameters of all models are shown in A.

##  4 Neutrino emission

Following the CR acceleration at the collimation shock, internal shocks,
forward shock and the shear layer, the CR interact with the background photon
field, and generate neutrino emission. The photon field in mildly relativistic
radiation-mediated shocks relaxes to a Wien spectrum as the rest-frame
downstream temperature is regulated by pairs to peak at $\sim 50$ keV
(consistent with the observed peak of $\sim 100$ keV), independent of the
Lorentz factor of the shock (e.g., Budnik et al., 2010; Nakar & Sari, 2012;
Levinson & Nakar, 2019). We thus assume the photon spectrum to be a Wien
spectrum at the collimation shock (Gottlieb et al., 2019) and during the shock
breakout (Gottlieb et al., 2018). Internal shocks may broaden the spectrum to
become a Band-like function, and thus for estimating the $p\gamma$ interaction
at the internal shocks, we assume a Band spectrum with $\alpha_{b}=-1.2$ and
$\beta_{b}=-2.3$ (Goldstein et al., 2013)555The cooling time of the internal
shocks is dominated by synchrotron cooling, thus the specific shape of the
spectrum does not significantly affect the results.. We normalize the photon
spectrum based on the available internal energy flux, $\sim 4p_{t}c$ at the
relevant acceleration site, as found from the simulation (or from GRB170817A
for the case of shock breakout). We consider both diffusive shock acceleration
in §4.1 and shear acceleration in §4.2.

###  4.1 Diffusive shock acceleration

The Fermi process accelerates protons to a power-law distribution
$dN_{p}/dE_{p}\propto E_{p}^{-s}$, where $s$ is the spectral index. We assume
a power-law index $s=2$, which roughly holds for mildly relativistic regimes,
as it is expected to be $2.0\lesssim s\lesssim 2.2$ (e.g., Blandford &
Eichler, 1987; Keshet & Waxman, 2005; Sironi et al., 2013). The rest-frame CR
spectrum is:

$E_{p}^{\prime
2}\frac{dN_{p}}{dE_{p}^{\prime}}=\frac{\epsilon_{CR}\epsilon_{d}\,{E^{\prime}_{k,\rm
iso}}}{{\rm log}(\,{E^{\prime}_{p,\rm max}}/{\rm
GeV})}e^{-\frac{E_{p}^{\prime}}{\,{E^{\prime}_{p,\rm max}}}}~{},$ (1)

where we assumed the minimal CR energy to be the proton rest-mass energy,
$m_{p}c^{2}\approx{\rm GeV}$. The upstream kinetic energy dissipation
efficiency is $\epsilon_{d}=0.2$ in mildly-magnetized shocks (§2) and
$\epsilon_{d}=1$ in shock breakout. We take $\epsilon_{CR}=0.1$ as the
fraction deposited in accelerating CR from the rest frame isotropic equivalent
dissipated kinetic energy reservoir $\epsilon_{d}\,{E^{\prime}_{k,\rm iso}}$.

The CR are accelerated by the Fermi process on a comoving timescale,

$t^{\prime}_{acc}(E_{p})\approx\frac{E^{\prime}_{p}}{\xi_{acc}qB^{\prime}c}\approx
10^{-6}\frac{E^{\prime}_{p,12}}{B^{\prime}_{6}}~{}{\rm s}~{}.$ (2)

In the following we use $\xi_{acc}=0.1$ (Globus et al., 2015).

The maximal energy to which the protons are accelerated $\,{E^{\prime}_{p,\rm
max}}$, is governed by the cooling time of the protons. In general, several
cooling processes such as Bethe-Heitler pair production, synchrotron
radiation, inverse-Compton radiation, adiabatic cooling, $pp$ and $p\gamma$
interactions should be considered. Close to the base of the jet, the high
photon number density in the region of the collimation shock renders the
photohadron interaction timescale to be the shortest one above the $p\gamma$
energy threshold. We integrate numerically the $p\gamma$ cooling time
(Berezinskii et al., 1990),

$\,{t_{p\gamma}}(E_{p})=\frac{2\gamma_{p}^{2}}{c}\bigg{[}\int_{\acute{E}_{\gamma,{\rm
min}}}^{\infty}\acute{E_{\gamma}}\kappa_{p\gamma}\sigma_{p\gamma,i}d\acute{E_{\gamma}}\int_{\frac{\acute{E_{\gamma}}}{2\gamma_{p}}}^{\infty}E_{\gamma}^{-2}h^{-1}dn^{\prime}_{\gamma,\nu}dE_{\gamma}\bigg{]}^{-1}~{},$
(3)

where $h$ is the Planck constant, $\gamma_{p}=E_{p}/m_{p}c^{2}$ is the Lorentz
factor of the proton, $n^{\prime}_{\gamma,\nu}$ is the downstream spectral
photon number density, $\acute{E}_{\gamma}$ is the photon energy in the proton
frame of reference, $\kappa_{p\gamma,\pi}$ and $\sigma_{p\gamma,\pi}$ are the
inelasticity and cross section of particle $i$ in the $p\gamma$ interaction,
respectively666For the pion cross section and elasticity we use the two-region
approximation (Atoyan & Dermer, 2003): $\sigma_{p\gamma,\pi}=340~{}\mu{\rm
b}$, and $\kappa_{p\gamma}=0.2$ when $\,{E_{\gamma,{\rm
min}}}<E_{\gamma}<3\,{E_{\gamma,{\rm min}}}$, and
$\sigma_{p\gamma,\pi}=120~{}\mu{\rm b}$, and $\kappa_{p\gamma}=0.6$ when
$E_{\gamma}>3\,{E_{\gamma,{\rm min}}}$. For the kaon cross section we use
Hümmer et al. (2010): $\sigma_{p\gamma,K^{+}}=2.0;~{}3.7;~{}2.7~{}\mu{\rm b}$
when: $6.7\,{E_{\gamma,{\rm min}}}<E_{\gamma}<8\,{E_{\gamma,{\rm
min}}};~{}8\,{E_{\gamma,{\rm min}}}<E_{\gamma}<11\,{E_{\gamma,{\rm
min}}};~{}11\,{E_{\gamma,{\rm min}}}<E_{\gamma}$, respectively..

At large distance from the jet launchpad, where the photon number density
drops, synchrotron cooling also plays a role. The synchrotron timescale is
given by

$t_{syn}(E_{p})=\frac{6\pi m_{p}^{4}c^{3}}{m_{e}^{2}\sigma_{T}B^{\prime
2}E^{\prime}_{p}}~{},$ (4)

where $m_{e}$ is the electron mass. Finally, as the outflow expands and breaks
out, the photon density drops and cooling by proton-proton ($pp$) interaction
may also become important, depending on the ejecta density,

$t_{pp}=\frac{m_{p}}{\sigma_{pp}\kappa_{pp}c\rho_{\rm bo}}~{},$ (5)

where $\rho_{\rm bo}$ is the mass density of the breakout layer,
$\kappa_{pp}\approx 0.5$ is the proton-proton inelasticity, and
$\sigma_{pp}\approx 5\times 10^{-26}~{}{\rm cm}^{2}$ is the proton-proton
cross section at the relevant energies. The total cooling time is $\,{t_{{\rm
cool}}}=(\,{t_{p\gamma}}^{-1}+t_{syn}^{-1}+t_{pp}^{-1})^{-1}$.

The photohadron process becomes important as the $p\gamma$ cross section
increases rapidly for photons with energies above $E_{\rm th}\approx{1.4\times
10^{17}{\rm eV^{2}}}/{E_{p}}$, measured in the proton comoving frame. Above
the resonance the interaction yields pions (kaons) which may decay into
neutrinos. The efficiency of the process depends on three suppression factors:
i) whether $p\gamma$ or $pp$ is the dominant cooling process:
$f_{p\gamma,i}=\frac{\,{t_{{\rm cool},i}}}{\,{t_{p\gamma,i}}}$, where
$i=\pi(K^{+})$ hereafter, and similarly for $pp$ process,
$f_{pp,\pi}=\frac{\,{t_{{\rm cool},\pi}}}{t_{pp,\pi}}$; ii) whether the
secondary particle decays before cooling down: $\,{f_{\rm sup,i}}=1-{\rm
exp}\left(-\frac{\,{t_{{\rm cool},i}}}{\,{t_{{\rm decay},i}}}\right)$, where
$\,{t_{{\rm decay},\pi}}=2.6\times 10^{-8}$ s and $\,{t_{{\rm
decay},K^{+}}}=1.2\times 10^{-8}$ s; iii) the comoving number density of the
photon field, $n_{\gamma}$: for each proton energy, we estimate the fraction
of protons that participate in the $p\gamma$ interaction by taking the number
of the mean free paths along the proper length of the photon target
$R^{\prime}$, roughly the width of the shock, which is on order of Thomson
optical depth (e.g., Nakar & Sari, 2012):

$f_{p,i}(E_{p})=R^{\prime}\sigma_{p\gamma,i}\int_{E_{\rm
th}}n_{\gamma,\nu}d\nu\approx
R^{\prime}f_{\gamma}(E_{p})n_{\gamma}\sigma_{p\gamma,i}~{},$ (6)

where $\sigma_{p\gamma,i}$ is the $p\gamma$ cross section for particle $i$,
and $f_{\gamma}(E_{p})$ is the fraction of the photons that satisfy
$E_{\gamma}>E_{\rm th}$, which we calculate numerically. If $f_{p}(E_{p})>1$,
we take $f_{p}(E_{p})=1$.

In the relevant energy range, $E\gtrsim 10$ GeV, the inelasticity of $p\gamma$
and $pp$ is $\kappa\sim 0.5$, namely half of the protons energy is deposited
in the produced pions/kaons. In the main $p\gamma$ ($pp$) channel,
$f_{p\rightarrow\pi^{\pm}}=\frac{1}{2}(\frac{2}{3})$ of the interactions
result in charged pions (Anchordoqui et al., 2014) that decay to neutrinos via
$\pi^{+}\rightarrow\mu^{+}+\nu_{\mu}\rightarrow
e^{+}+\nu_{e}+\nu_{\mu}+\bar{\nu}_{\mu}$ and
$\pi^{-}\rightarrow\mu^{-}+\bar{\nu}_{\mu}\rightarrow
e^{-}+\bar{\nu}_{e}+\nu_{\mu}+\bar{\nu}_{\mu}$. The energy is shared rather
equally between the four leptons such that
$f_{\pi^{\pm}\rightarrow\nu_{\mu}}\approx 0.5$ of the charged pion energy is
transferred to $\nu_{\mu}$ and $f_{\pi^{\pm}\rightarrow\nu_{e}}\approx 0.25$
of it in $\nu_{e}$. Averaging over the neutrino flavors and neglecting the
weak dependency upon the photon energy, the maximal rest-frame neutrino energy
is $\kappa
f_{p\rightarrow\pi^{\pm}}(f_{\pi^{\pm}\rightarrow,\nu_{\mu}}+f_{\pi^{\pm}\rightarrow,\nu_{e}})E_{p}\approx\frac{3}{16}E_{p}~{}(\frac{1}{4}E_{p})$
for $p\gamma$ ($pp$) interactions, respectively.

When strong magnetic fields are present, pions cool down rapidly and their
suppression factor becomes small. Then kaons, whose mass is larger and thus
cool slower, become important. The main channel ($f_{p\rightarrow
K^{+}}=0.63$) is $K^{+}\rightarrow\mu^{+}+\nu_{\mu}$, where $\nu_{\mu}$
carries $f_{K^{+}\rightarrow\nu_{\mu}}=0.48$ of the kaon energy (Lipari et
al., 2007). We neglect the energy of the $\mu^{+}$ since it cools down faster
than pions, as well as other $K^{+}$ channels as their contribution is at most
a few percent, Overall, the maximal rest-frame energy of neutrinos produced by
kaons is $\kappa f_{p\rightarrow
K^{+}}f_{K^{+}\rightarrow\nu_{\mu}}E_{p}\approx\frac{1}{6}E_{p}$.

The total observed neutrino energy flux is

$\begin{split}E_{\nu}^{2}\frac{dN_{\nu}}{dE_{\nu}}\approx&\bigg{[}f_{p,\pi}\,{f_{\rm
sup,\pi}}\left(\frac{3}{16}f_{p\gamma,\pi}+\frac{1}{4}f_{pp,\pi}\right)+\\\
&f_{p,K^{+}}\,{f_{\rm
sup,K^{+}}}\frac{\sigma_{p\gamma,K^{+}}}{\sigma_{p\gamma,\pi}}\frac{1}{6}f_{p\gamma,K^{+}}\bigg{]}\delta_{D}E_{p}^{\prime
2}\frac{dN_{p}}{dE_{p}^{\prime}}~{},\end{split}$ (7)

where $\delta_{D}\equiv[\Gamma(1-\beta{\rm cos}(\theta-\theta_{\rm
obs})]^{-1}$ is the Doppler factor of the shock, $\beta$ is the dimensionless
velocity, $\theta_{\rm obs}$ is the viewing angle, $\sigma_{p\gamma,\pi}$,
$\sigma_{p\gamma,K^{+}}$ are the $p\gamma$ cross sections for pions and kaons,
respectively, $\,{f_{\rm sup,i}}$ is the suppression factor, $f_{pp,i}$ and
$f_{p\gamma,i}$ are the $pp$ and $p\gamma$ dominant cooling factors,
respectively, and $f_{p,i}$ is given by Eq. 6.

###  4.2 Shear acceleration

In addition to DSA, another process through which neutrinos can be produced is
shear acceleration. In GRBs the most prominent shear region is at the
interface between the jet and the cocoon. The shear acceleration time is
(Rieger & Duffy, 2004),

$t_{acc}^{s}=\frac{1}{(4+\zeta)\chi^{s}\tau}=\frac{\xi_{acc}c}{5\chi^{s}r_{L}}~{},$
(8)

where $\tau$ is the mean scattering time, in its general form
$\tau=\tau_{0}p^{\zeta}$. In the R.H.S we used Bohm diffusion, $\zeta=1$ and
$\tau=r_{L}/\xi_{acc}c$. $\chi^{s}$ is the shear coefficient given for a
relativistic flow as (Rieger & Duffy, 2004):

$\chi^{s}=\frac{(\Gamma
c)^{2}}{15}\bigg{(}\frac{\partial\beta_{z}(r)}{\partial
r}\bigg{)}^{2}\Big{(}1+(\Gamma\beta_{z})^{2}\Big{)}~{},$ (9)

where $\Gamma$ and $\beta_{z}$ are the bulk Lorentz factor and the
dimensionless velocity component in the direction $z$ of the jet axis,
respectively, and $r$ is the radial coordinate of the jet cross-section. In
order for the acceleration to be efficient, $t_{acc}^{s}$ has to be shorter
than the cooling time, here $\,{t_{p\gamma}}$,

$\chi_{12}^{s}E_{p,12}\beta_{\bot}\gtrsim 6\times 10^{3}B_{6}~{},$ (10)

where $\beta_{\bot}$ is the dimensionless velocity component that is
perpendicular to the magnetic field. Since magnetized jets maintain $B\gg
10^{6}$ G, and since unmagnetized jets are prone to mixing which smears the
shear gradient thereby reducing $\chi^{s}$, the only possible shear
acceleration site is the unmagnetized/weakly magnetized collimation shock. Our
hydrodynamic simulations show that the shear layer between the shocked jet
material and the cocoon can yield a typical velocity gradient from $\Gamma\sim
100$ to $\beta\sim 0.1$ over $\Delta r\sim\theta_{j}R_{s}\sim 10^{8}$
cm777Note that in our hydrodynamic simulation the jet is injected at $\sim
10^{7}$ cm, in principle the conical acceleration can start at $r\approx
r_{g}$, thereby having $\Delta r$ that is smaller by roughly an order of
magnitude, such that $\chi_{12}^{s}\approx 10^{2}$., which leads to
$\chi_{12}^{s}\approx 1$. Consequently, we require $E_{p}\gtrsim 10^{4}$ TeV
to satisfy Eq. 10, however at these energies and for $B\approx 10^{6}$ G, the
Larmor radius becomes $\gtrsim 10^{7}$ cm, larger than the transversal
confinement size. We conclude that neutrino production by shear acceleration
is unlikely, as at low energies the shear acceleration timescale is longer
than the cooling time, and at high energies the Larmor radius becomes too
large.

###  4.3 Neutrino spectrum

Following §2 and §3, collisionless subshocks may form if the radiation-
mediated shocks are mildly magnetized, $\sigma_{u}\gtrsim 0.1$, as we find at
the collimation shock in model $\,{\it WM}$ and internal shocks in model
$\,{\it HM}$. If the jet is launched hydrodynamically, a collisionless shock
can only form after jet/cocoon breakout. For highly magnetized jets,
$B^{\prime}\gtrsim 10^{11}$ G at the collimation nozzle, which renders the
neutrino production inefficient as the synchrotron cooling timescale is much
shorter than their decay time, $\,{f_{\rm sup,\pi}}\ll 1$. We thus consider
neutrino emission from the following sites: collimation shock in weakly
magnetized jets (white ellipses in Fig. 1k,m), internal shocks from highly
magnetized jets (e.g. black circles in Fig. 1a,d), and shock breakout from
choked and successful hydrodynamic jets (Fig. 2 for a successful jet
breakout). In Appendix B, we describe the calibration of the latter by
GRB170817 observations.

Figure 2: Logarithmic density map ${\rm log}_{10}(\rho)~{}[{\rm g~{}cm^{-3}}]$
of the jet breaking out from the ejecta in the hydrodynamic model $\,{\it
SU2}$, 2.8 s after the jet was launched. Since in this simulation the jet
engine was working for $t_{j}=2$ s, the unshocked jet material is at
$z=0.8{\rm s}\times c=2.4\times 10^{10}~{}{\rm cm}$.

Figure 3:  The acceleration (black) and cooling ($p\gamma$ in blue and
synchrotron in red) times in successful jets, at the following production
sites for highly magnetized (HM) and weakly magnetized (WM) jets: collimation
shock (WM) at different heights, $3\times 10^{9}$ cm in solid lines and
$3\times 10^{10}$ cm in dashed lines (panel a); internal shocks (HM, panel b);
and shock breakout (hydrodynamic jets, panel c). The thin blue lines in panels
a,b reflect the $p\gamma$ cooling time when using the cross section for kaons,
showing that $p\gamma$ dominates for collimation and internal shocks also for
kaons (the synchrotron time in panel a is for the shock at $3\times 10^{9}$
cm, at $3\times 10^{10}$ it becomes negligibly long). In panel c we also show
the $pp$ cooling time (yellow) where we assume that the ejecta is composed
only of protons. The typical $p\gamma$ cooling time is shown both for the jet
(solid blue) and for the cocoon (dotted blue).  Figure 4:  Neutrino signals
from a BNS merger at 41 Mpc. Shown are the neutrino spectral fluences from the
collimation shock at different heights (green and turquoise), internal shocks
(light blue), which originate in the jet and can be observed within the jet
opening angle. In purple is shown the shock breakout characteristic emission
from the jet and the cocoon at viewing angle $\theta_{\rm obs}=20^{\circ}$
(successful and choked jet simulations produce similar neutrino emission from
the cocoon). Solid lines mark the total neutrino emission and dashed lines
mark the contribution from the decayed kaons. Upper limits (at 90% CL) from
various instruments on the neutrino spectral fluence from GW170817 during a
$\pm 500$ s window centered on the gravitational wave (GW) trigger time are
shown in grey and taken from Aartsen et al. (2020). All fluences are shown as
the per-flavor sum of neutrino and anti-neutrino fluence, assuming equal
fluence in all flavors, as expected for standard neutrino oscillation
parameters. At production sites where we use a Wien spectrum (collimation
shock and shock breakout), the fluence rises at $\sim 10^{2}$ GeV due to the
maximal photon energy, and falls at the maximal CR energy, as shown in panels
(a-c). For internal shocks, where we use a Band function, the neutrino fluence
is quasi-flat at low CR energies.

Fig. 3 depicts the acceleration (black) and relevant cooling times as a
function of the proton energy in the comoving fluid frame at the following
sites: collimation shocks for: $R_{s}=3\times 10^{9}$ cm (solid) and
$R_{s}=3\times 10^{10}$ cm (dashed) (panel a); internal shocks (panel b); jet
(solid) and cocoon (dotted) shock breakout (panel c). The $p\gamma$ cooling
time (charged pions in thick lines and kaons in thin lines) is dominant at all
production sites, both for pions and kaons, and thus the suppression factor is
always of order unity. We note that at later times of the shock breakout, the
photons density drops, causing the $p\gamma$ process to become subdominant.

|  | Hydrodynamic | weakly magnetized | Highly magnetized |  Choked†
---|---|---|---|---|---
Acceleration | Production site | Peak fluence |  Max. $E_{\nu}$ | Peak fluence |  Max. $E_{\nu}$ | Peak fluence |  Max. $E_{\nu}$ | Peak fluence |  Max. $E_{\nu}$
|  | $\left[{\rm\frac{GeV}{cm^{2}~{}s~{}sr}}\right]$ |  [${\rm GeV}$] | $\left[{\rm\frac{GeV}{cm^{2}~{}s~{}sr}}\right]$ |  [${\rm GeV}$] | $\left[{\rm\frac{GeV}{cm^{2}~{}s~{}sr}}\right]$ |  [${\rm GeV}$] | $\left[{\rm\frac{GeV}{cm^{2}~{}s~{}sr}}\right]$ |  [${\rm GeV}$]
Fermi | col. shock | - | $10^{-1}$ | $10^{3}$ | - | -
internal shocks | - | $10^{-6}$ | $10^{3}$ | $10^{-6}$ | $10^{3}$ | -
shock breakout | $10^{-4}$ | $10^{5}$ | $10^{-4}$ | $10^{5}$ | $10^{-4}$ | $10^{5}$ | $10^{-7}$ | $10^{5}$
Shear acceleration | -
Table 1:  Prospects of neutrino maximal fluence at $41$ Mpc, for different
types of jets at different production sites. In sites where neutrino
production is possible, we show the maximal neutrino energy and flux.
† Choked jets refer to the production of neutrino after the jet is choked.

The resulting neutrino fluence from a GW170817-like event at GW170817 distance
of $D=41$ Mpc (Hjorth et al., 2017) is shown in Fig. 4. A summary of the
neutrino peak fluence for different types of jets and at different production
sites is given in Table 1. It is shown that the neutrino signals are composed
of two components: a strong low energy component and a weaker high energy
tail, which correspond to the pion and kaon contributions, respectively. The
latter emerges thanks to the smaller $p\gamma$ cross section of kaons which
allows CR to accelerate to higher energies prior to the $p\gamma$ interaction
(Fig. 3a,b). Since kaons are heavier than pions, they have a longer
synchrotron cooling time by $\sim 2$ orders of magnitude. When the
magnetization is high, the pions cool down rapidly, and thus the kaon
contribution becomes important at all neutrino energies.

Upon launching, the high magnetization and photon density prevent the
production of high energy neutrinos. At later times, the magnetic field at the
internal shocks drops, but so is their magnetization such that their ability
to support a strong collisionless subshock is questionable. If the jet engine
is working for a long time ($\gtrsim 2$ s) and the jet initial magnetization
is mild, then a large ($\sim 10^{10}$ cm) collimation shock forms, where pions
decay before undergoing synchrotron cooling. This results in our highest
fluence of neutrinos up to $\sim 10$ TeV. However, such a scenario is unlikely
as it requires jet launching with subdominant magnetic fields, and atypically
long duration of jet launching.

The neutrinos emitted from the shocks inside the ejecta are confined to the
jet opening angles. In the shock breakout scenario, although the strongest
neutrino signal emerges from shocks inside the jet, a weaker signal also
emerges from the shock breakout of the mildly relativistic cocoon, at viewing
angles $\lesssim 20^{\circ}$. The latter is also a robust estimate for the
neutrino fluence, as it is calibrated by GRB170817A observations, and weakly
dependent on the jet magnetization and fate (successful and choked jets have
similar neutrino emission from the cocoon). Shock breakout is also the only
neutrino production site for hydrodynamic jets, or choked jets (see Appendix
C).

##  5 Cosmic-ray production

We consider the CR acceleration at the external shock wave. Once the
structured jet (or cocoon, if the jet is choked) completes the initial
acceleration phase, it is coasting freely in the interstellar medium (ISM).
Assuming the ISM has a constant number density $n_{0}$, and assuming an
independent hydrodynamic evolution of each angle $\theta$, the deceleration
spherical radius of the shock is

$r_{dec}(\theta)=\bigg{(}\frac{3M(\theta)}{4\pi
n_{0}m_{p}\Gamma(\theta)}\bigg{)}^{1/3}~{},$ (11)

where $M(\theta)$ is the isotropic equivalent swept mass of the outflow at
angle $\theta$. As the shock decelerates, it deposits a fraction of its
dissipated energy into accelerating electrons in the magnetic field generated
by the shock. The accelerated electrons emit synchrotron radiation to form the
afterglow signal. In the following, we use a numerical model that fits
GW170817A’s afterglow observations, as found by Mooley et al. (2018). We also
consider a case of an uncollimated choked jet (models $\,{\it SU1}$ and
$\,{\it CU}$ in Appendix A). Fig. 5 depicts maps of kinetic energy density
(left panels) and $\Gamma\beta$ (right panels) for these models.

The maximal energy at which CR can be accelerated at the shock is estimated by
equating the acceleration time
$t^{\prime}_{acc}=t^{\prime}_{L}/\xi_{acc}=E^{\prime}_{p}/(ZqB^{\prime}c\xi_{acc})$
with the escape time $t^{\prime}_{esc}=R^{\prime}/(\Gamma\beta c)$ (Hillas,
1984), where $R^{\prime}$ is the proper acceleration length, $q$ is the
elementary charge, and $\xi_{acc}\approx 0.1$ is the efficiency of the
acceleration in the mildly relativistic case (e.g., Globus et al., 2015). We
find

$\displaystyle\,{E^{\prime}_{p,\rm max}}(\theta)$ $\displaystyle=$
$\displaystyle\frac{\xi_{acc}R^{\prime}(\theta)ZqB^{\prime}(\theta)}{\Gamma(\theta)\beta(\theta)}$
(12) $\displaystyle\approx$ $\displaystyle 300\,\xi_{acc}ZB^{\prime}_{\rm
G}(\theta)R_{\rm
cm}(\theta)\left[\Gamma(\theta)\beta(\theta)\right]^{-1}\,{\rm eV}\,.$

Note the angular-dependent quantities since we perform a semi-analytic
calculation on structured jets. We assume the acceleration length to be the
proper shock width, $R^{\prime}\approx{r_{dec}}/{12\Gamma}$ (e.g., van Eerten,
2015), and use the general convention that the magnetic energy density behind
the shock can be characterized by the equipartition parameter $\epsilon_{B}$,
such that the magnetic field is
$B^{\prime}\approx\left[32\pi\epsilon_{B}n_{0}m_{p}\Gamma(\Gamma-1)\right]^{1/2}c$.
Plugging $B^{\prime}$ and Eq. 11 into Eq. 12, the maximal energy of the
accelerated particle is:

$\footnotesize\,{E^{\prime}_{p,\rm max}}(\theta)\approx 2\times
10^{7}~{}\frac{\xi_{acc}Z(\Gamma(\theta)-1)^{1/2}}{\beta(\theta)}\left(\frac{n_{0}M_{-6}(\theta)^{2}\epsilon_{B,-3}^{3}}{\Gamma(\theta)^{11}}\right)^{1/6}~{}{\rm
GeV}~{}\,,$ (13)

where as before, $Q_{x}$ denotes the value of the quantity $Q$ in units of
$10^{x}$ times its c.g.s. units, except for $M_{x}$ which is given in units of
solar masses.

As the physical quantities have an angular dependence (see Fig. 5), we
discretize the shock wave into bins in $\theta$ and $\Gamma$ (on the radial
dimension), on which the quantities do not vary much. We assume that the
angles can be treated independently of each other (e.g., no spreading), and
that the mass in the front decelerates to the Lorentz factor of the mass that
supports it behind. That is, the cumulative mass from radius $r$ is moving
with $\Gamma(r)$. We calculate the cumulative isotropic equivalent mass over
the radial direction to calculate the maximal CR energy from Eq. 13 for each
angle and proper-velocity, which corresponds to different times. We then use
Eq. 1888Note that here we perform a full calculation at each angle, and thus
we use the total energy rather than the isotropic equivalent energy. to
calculate the comoving CR spectrum at each angle and proper-velocity. We boost
the rest-frame spectra to the observer by accounting for the beaming and the
Doppler boost for each observer. We ignore the temporal evolution of the
emission, and integrate over time (radial direction) since the CR diffusion
time is longer than the emission times. Finally, we average the angle-
dependent flux over the beaming factor as CR do not propagate in straight
lines in the Galactic magnetic field, and obtain the CR spectrum from a single
source.

Figure 5:  Maps of successful (a) and uncollimated choked (b) jets during the
free coasting phase. Shown are the logarithmic kinetic energy density (left)
and the logarithmic proper-velocity (right).  Figure 6: Total CR luminosity in
successful (a) and choked jets (b) from external shocks from BNS mergers,
assuming that 90% of the energy dissipated into CR is deposited in light
components (protons, in blue), and 10% is in heavy (Fe, in pink) elements. We
assume that the rate of both successful and choked jets is the GBNS merger
rate ${\cal R}_{MW}=42^{+30}_{-14}$ Myr-1. We use $\xi_{acc}=0.1$,
$\epsilon_{CR}=0.1$, and two values of $\epsilon_{B}$, that of GW170817,
$\epsilon_{B,-3}=1$ (strong colors), and the more commonly used
$\epsilon_{B,-3}=100$ (faint colors). The horizontal bars denote the Galactic
luminosity in CR in different energy bins, estimated for protons and iron. In
the successful jet scenario (calibrated by GW170817), the BNS account for
5-10% of the total Galactic CR luminosity. This also leads to an upper limit
of the GBNS merger rate, ${\cal R}_{MW}=5\times 10^{-4}$ yr-1 (dashed lines).
For the choked jets, the CR luminosity is higher, but poorly constrained due
to the unknown rate and energetics of such events.

To estimate the total CR luminosity available from GBNS mergers, we use the
rate constrained by the known sample of GBNS systems, ${\cal R}_{MW}$
($42^{+30}_{-14}$ Myr-1; Pol et al. 2019), in agreement with LIGO estimates
${\mathcal{R}_{\rm BNS}}=320^{+490}_{-240}~{}{\rm Gpc}^{-3}{\rm yr}^{-1}$
(Abbott et al., 2021), which we assume to be the rate of both successful jets
and choked jets. Using $\epsilon_{CR}=0.1$, $\xi_{acc}=0.1$, and circum-merger
number density $n_{0}=1\,{\rm cm}^{-3}$ (Berger, 2014), we calculate the CR
luminosity from BNS mergers, assuming that protons constitute 90% of the
energy and Fe account for the remaining 10% (see Kimura et al., 2018b, for a
discussion about the CR composition from the kilonova ejecta). Note that we do
not aim to accurately fit the composition as this depends on many parameters
(e.g., the jet parameters, microphysics of the shock and BNS merger rate), but
rather compare the CR energy from our model with the observed CR in the Milky
Way, which is estimated as follows. The kinetic energy density of cosmic rays
in the Milky Way is $\rho={4\pi}{c}^{-1}\int\Phi(E)dE$. Up to the knee, the
Galactic CR flux is $\Phi(E)\sim 1.8\times 10^{4}(E/1{\rm GeV})^{-1.7}\rm
m^{-2}s^{-1}sr^{-1}$ and $\rho_{\rm GeV-PeV}\sim 1\rm eV\,cm^{-3}$. Above the
knee, $\Phi(E)\sim 1.8\times 10^{-6}(E/1{\rm PeV})^{-2.1}\rm
m^{-2}s^{-1}sr^{-1}$ and $\rho_{\rm PeV-EeV}\sim 6.8\,10^{-5}\rm eV\,cm^{-3}$.
The measured ratio of radioactive isotopes and the secondary-to-primary ratios
indicate a residence time of CR in the Galactic disk of volume
$V\sim{\pi(15\rm kpc)^{2}}(200\rm pc)$, of the order $\tau\sim 30\,{\rm
Myr}\left[{(E/Z)}/{\rm 5GeV}\right]^{-\delta}$ with $\delta\sim 0.33$.
Assuming protons, $Z=1$, the power required to maintain the observed Galactic
cosmic ray flux in the GeV-PeV range is $\rho V/\tau\sim 4.5\times
10^{39}~{}{\rm erg\,s}^{-1}$, a few percent of the total kinetic power of
supernovae. Above the knee, in the PeV-EeV energy range, the required power is
only $\sim 2.7\times 10^{37}~{}{\rm erg\,s}^{-1}$ but the sources need to
accelerate CR to higher energies than typical supernovae remnants ($\sim 100$
TeV).

Fig. 6 depicts the resulting CR luminosity for BNS merger rate ${\cal
R}_{MW}=42^{+30}_{-14}$ Myr-1 and $\epsilon_{B,-3}=1$ (dark shaded area) as
constrained by GW170817’s afterglow (Mooley et al., 2018), and
$\epsilon_{B,-3}=100$ (light shaded area). We find that blast-waves of sGRBs
may account for $\sim 5-10$% of the CR luminosity in the Milky Way in the PeV-
EeV energy range for GW170817-like events if the jet is present. The maximum
energy for the iron component, $\sim 2\times 10^{17}$ eV ($10^{18}$ eV) for
$\epsilon_{B,-3}=1$ ($\epsilon_{B,-3}=100$). Choked jets seem to be able to
account for all Galactic CR flux, owing to their quasi-spherical nature which
increases the angle averaged emission999We also note that in our simulation
the uncollimated choked jet has 2.5 times more energy than that of GW170817
(Appendix A).. However, we stress that this result is poorly constrained as we
assume the rate and the energy of choked jets to be comparable to that of
successful jets, but those quite possibly are lower for choked jets.

The dipole anisotropy for a single transient source has an amplitude
$\Delta=3r/2ct$ for an event occurring at a time $t$ in the past and located
at a distance $r$. The upper limit on the observed dipole anisotropy in the
PeV-EeV range is $\Delta\sim 0.02$. If the latest GBNS occurred at $1/{\cal
R}_{MW}$ in the past, with ${\cal R}_{MW}=42^{+30}_{-14}$ Myr-1, its distance
would be in the range 0.6-1.5 kpc (1-3 kpc) if the BNS mergers account for 10%
(5%) of the Galactic CR flux in the PeV-EeV energy range.

##  6 Summary of results

We performed one GRMHD simulation of a highly magnetized jet and additional
RMHD simulations of a weakly magnetized jets to estimate neutrino and CR
signals from hydrodynamic/magnetized and choked/successful jets. Our main
findings are:
$\bullet$ collisionless shocks cannot form below the photosphere due to mixing
between jet and cocoon material; collisionless subshocks can form under
certain conditions that depend on the jet magnetization;
$\bullet$ hydrodynamic jets and choked jets cannot produce neutrinos prior to
jet-breakout;
$\bullet$ if the jet is operating for atypically long time for sGRBs ($\gtrsim
2$ s) and its magnetization is mild at the collimation shock, the collimation
shock can produce $\sim 10$ TeV neutrinos;
$\bullet$ internal shocks, for which we use those that emerge in the full
GRMHD simulation, yield a signal well below detection limits;
$\bullet$ jet/cocoon shock breakout, which is considered here as a source of
neutrino in BNS for the first time, is independent of the jet magnetization
and the only production site for hydrodynamic jets. It can produce $\sim 10$
PeV neutrinos to large angles, on order of the cocoon opening angle of $\sim
20$∘; however, the emission is $\sim 100$ times lower than collimation shocks
at 1 PeV and would require events within our local group or more sensitive
neutrino detectors;
$\bullet$ the shear layer does not produce high energy neutrinos;
$\bullet$ external shocks, responsible for the afterglow emission, can
accelerate CR up to $\sim 10^{18}$ eV.

To summarize, we find that sGRBs do not produce a neutrino signal that is
detectable by current-day facilities. We find that GBNS can account for 5-10%
of the Galactic CR flux in the PeV-EeV range, depending on the GBNS merger
rate, jet energetics and CR composition.

##  7 Discussion

As long as the jet is propagating inside the ejecta, its mixing with the
cocoon plays a crucial role in the jet dynamics, and inhibits the formation of
collisionless subphotospheric shocks in the jet. If the jet maintains some
non-negligible degree of magnetization, collisionless subshocks may form,
allowing CR to be accelerated efficiently if the Larmor radius is larger than
the subshock width. We found that the strongest neutrino flux emerges from the
collimation shock in the weakly-magnetized jet, however its existence is
highly questionable due to the atypical engine work-time and the fine-tuned
magnetization. Internal shocks suffer from heavy mixing at early times, when
the jet is still sufficiently magnetized to support the formation of
collisionless shocks, and thus do not yield a promising neutrino signal
either.

Our results are in agreement with those of Biehl et al. (2018) who found that
neutrino signals from sGRBs are unlikely to be detected by current-day
facilities, and with Kimura et al. (2018a) who found that collimation shocks
of sGRBs are unlikely to produce a detectable signal. However, these studies,
which have addressed the neutrino emission analytically and thus overlooked
the effects of mixing and jet magnetization on neutrino production at
subphotospheric shocks, also found internal shocks to be the most promising
neutrino production sites. We showed that in realistic jet-cocoon outflows,
the production of neutrinos at subphotospheric shocks may only be possible if
the jet is magnetized. Our magnetized jet simulations show that the required
magnetization is only present at small radii, where the high photon density
prevents CR from efficiently accelerating to high energies. Consequently, we
find that it is unlikely that subphotospheric internal shocks provide any
detectable neutrino signal.

The jet/cocoon shock breakout is considered here as a source of neutrinos in
BNS for the first time. It is generally independent of the jet initial
magnetization, and its energy is calibrated by GRB170817A observations. The
neutrinos produced by the decay of kaons in the shock breakout can reach $\sim
10$ PeV, owing to low magnetization and photon density which allow protons to
accelerate to high energies. If the shock breakout is accompanied by a
collimated jet, such that the outflow breaks out with an angular structure,
the on-axis neutrino emission ($\theta\lesssim\theta_{j}$) is $\sim 2-3$ order
of magnitude brighter than that of the cocoon ($\lesssim 20^{\circ}$), which
arises regardless of the jet’s fate. The energy stored in the shock breakout
is smaller than in the other potential production sites, and thus such
emission may only become detectable in the future. We find that the neutrino
fluence from on-axis shock breakout is comparable with that from fallback
accretion of unbound ejecta in BNS mergers (Decoene et al., 2020). Note
however that neutrinos from fallback accretion are generated over much longer
timescales, such that the neutrino flux from the jet breakout emission is much
higher, and the flux from fallback accretion is comparable with that from the
cocoon shock breakout. Finally, unlike subphotospheric production sites where
neutrinos are anticipated to arrive before the photons, neutrinos from shock
breakout are candidates for a neutrino-photon coincident detection. At
distances of $\sim 100$ Mpc, the arrival times of TeV-PeV neutrinos and
photons to Earth is expected to be essentially the same (e.g., Wei et al.,
2016), allowing a coincident detection of the two, which can be a probe of the
neutrino production site.

Our estimate for the CR emission from the jet-cocoon blast wave, calibrated by
GRB170817A observations, have several important implications. First, the
Galactic CR luminosity sets an upper limit to the GBNS merger rate to ${\cal
R}_{MW}\sim 5\times 10^{-4}$ yr-1 for successful jets. Second, the level of
anisotropy in the PeV-EeV range constrains the distance from the last GBNS
merger. In a purely turbulent magnetic field of $\sim 1\mu G$, the typical
angular spread of CR sources at distances $\lesssim$ 1 kpc is $\lesssim 50$
degrees at 50 PeV, such that it might be possible to observe an intermediate-
scale anisotropy from CR sources emitting in the shin region from BNS at such
distances. Even if unlikely that the last BNS occured at such close distance,
it is worth mentioning that at 33 PeV, an intermediate-scale anisotropy has
been reported in the direction ($l\sim$80∘, $b\sim$15∘) in Galactic
coordinates (see figure 3 of Ahlers, 2019), and owing to their natal kicks, we
expect an anisotropy from a BNS merger to be located at higher latitudes than
other Galactic transient sources. In the more likely scenario that the latest
BNS merger occurred at a larger distance, our estimates suggest a distance of
$\lesssim 3$ kpc for the latest GBNS merger. Future combined search for CR
anisotropy and X-ray line emission from radioactive elements (e.g., Wu et al.,
2019) could pinpoint the latest BNS merger in our Galaxy.

Finally, our results also have a few interesting implications to long GRBs,
which share similar jet-cocoon mixing and jet structure with sGRBs (Gottlieb
et al., 2021a). Since long GRBs operate over longer engine work-times and thus
are also more energetic, it is possible that such systems can produce more
neutrinos, particularly in the larger collimation shock. Furthermore, these
conditions may help achieving optically thin shocks as suggested by Murase &
Ioka (2013), with the possible extended envelope is a promising production
site (Nakar, 2015). It is also worth mentioning that the dynamics of long GRB
jets propagating in their progenitor stars is not significantly different in
terms of the emerging jet-cocoon structure (Gottlieb et al., 2021b). Since the
energetics times the long GRB rate is higher than that of short GRBs, it
implies that the contribution of long GRBs to the Galactic PeV-EeV CR flux is
even higher than that of short GRBs.

## Acknowledgements

We are grateful to Amir Levinson, Anatoli Fedynitch, Kohta Murase, Andrew
MacFadyen, Roger Blandford, Denis Allard, Brian Metzger, Alexander
Tchekhovskoy and Matthew Liska for helpful discussions and useful comments. OG
acknowledges support by the Israel Science Foundation Grant No. 1114/ 17 and
ERC grant (JetNS). NG acknowledges support by the Czech Science Foundation
under the grant GACR 20-19854S titled “Particle Acceleration Studies in
Astrophysical Jets”, partial support from the National Science Foundation
through AST grant 1715661, and the Simons Foundation.

## References

  * Aartsen et al. (2020) Aartsen, M., Abbasi, R., Ackermann, M., et al. 2020, arXiv preprint arXiv:2008.04323
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, The Astrophysical Journal, 848, L13, doi: 10.3847/2041-8213/aa920c
  * Abbott et al. (2021) Abbott, R., Abbott, T., Abraham, S., et al. 2021, arXiv preprint arXiv:2101.12130
  * Ahlers (2019) Ahlers, M. 2019, The Astrophysical Journal Letters, 886, L18
  * Albert et al. (2017) Albert, A., André, M., Anghinolfi, M., et al. 2017, arXiv preprint arXiv:1710.05839
  * Anchordoqui et al. (2014) Anchordoqui, L. A., Barger, V., Cholis, I., et al. 2014, Journal of High Energy Astrophysics, 1, 1, doi: 10.1016/j.jheap.2014.01.001
  * Atoyan & Dermer (2003) Atoyan, A. M., & Dermer, C. D. 2003, ApJ, 586, 79, doi: 10.1086/346261
  * Beloborodov (2017) Beloborodov, A. M. 2017, ApJ, 838, 125, doi: 10.3847/1538-4357/aa5c8c
  * Berezinskii et al. (1990) Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., & Ptuskin, V. S. 1990, Astrophysics of cosmic rays
  * Berger (2014) Berger, E. 2014, ARA&A, 52, 43, doi: 10.1146/annurev-astro-081913-035926
  * Biehl et al. (2018) Biehl, D., Heinze, J., & Winter, W. 2018, MNRAS, 476, 1191, doi: 10.1093/mnras/sty285
  * Blandford & Eichler (1987) Blandford, R., & Eichler, D. 1987, Phys. Rep., 154, 1, doi: 10.1016/0370-1573(87)90134-7
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, Monthly Notices of the Royal Astronomical Society, 179, 433
  * Bromberg et al. (2011) Bromberg, O., Nakar, E., Piran, T., & Sari, R. 2011, The Astrophysical Journal, 740, 100, doi: 10.1088/0004-637X/740/2/100
  * Budnik et al. (2010) Budnik, R., Katz, B., Sagiv, A., & Waxman, E. 2010, ApJ, 725, 63, doi: 10.1088/0004-637X/725/1/63
  * Christie et al. (2019) Christie, I. M., Lalakos, A., Tchekhovskoy, A., et al. 2019, MNRAS, 490, 4811, doi: 10.1093/mnras/stz2552
  * Crumley et al. (2019) Crumley, P., Caprioli, D., Markoff, S., & Spitkovsky, A. 2019, MNRAS, 485, 5105, doi: 10.1093/mnras/stz232
  * Decoene et al. (2020) Decoene, V., Guépin, C., Fang, K., Kotera, K., & Metzger, B. D. 2020, J. Cosmology Astropart. Phys, 2020, 045, doi: 10.1088/1475-7516/2020/04/045
  * Dermer & Atoyan (2006) Dermer, C. D., & Atoyan, A. 2006, New Journal of Physics, 8, 122, doi: 10.1088/1367-2630/8/7/122
  * Fernández et al. (2019) Fernández, R., Tchekhovskoy, A., Quataert, E., Foucart, F., & Kasen, D. 2019, MNRAS, 482, 3373, doi: 10.1093/mnras/sty2932
  * Globus et al. (2015) Globus, N., Allard, D., Mochkovitch, R., & Parizot, E. 2015, Monthly Notices of the Royal Astronomical Society, 451, 751, doi: 10.1093/mnras/stv893
  * Globus & Levinson (2014) Globus, N., & Levinson, A. 2014, The Astrophysical Journal, 796, 26
  * Goldstein et al. (2013) Goldstein, A., Preece, R. D., Mallozzi, R. S., et al. 2013, ApJS, 208, 21, doi: 10.1088/0067-0049/208/2/21
  * Goldstein et al. (2017) Goldstein, A., Veres, P., Burns, E., et al. 2017, ApJ, 848, L14, doi: 10.3847/2041-8213/aa8f41
  * Gottlieb et al. (2021a) Gottlieb, O., Bromberg, O., Levinson, A., & Nakar, E. 2021a, MNRAS, 504, 3947, doi: 10.1093/mnras/stab1068
  * Gottlieb et al. (2020a) Gottlieb, O., Bromberg, O., Singh, C. B., & Nakar, E. 2020a, MNRAS, 498, 3320, doi: 10.1093/mnras/staa2567
  * Gottlieb et al. (2019) Gottlieb, O., Levinson, A., & Nakar, E. 2019, MNRAS, 488, 1416, doi: 10.1093/mnras/stz1828
  * Gottlieb et al. (2020b) —. 2020b, MNRAS, 495, 570, doi: 10.1093/mnras/staa1216
  * Gottlieb & Loeb (2020) Gottlieb, O., & Loeb, A. 2020, MNRAS, 493, 1753, doi: 10.1093/mnras/staa363
  * Gottlieb et al. (2021b) Gottlieb, O., Nakar, E., & Bromberg, O. 2021b, MNRAS, 500, 3511, doi: 10.1093/mnras/staa3501
  * Gottlieb et al. (2019) Gottlieb, O., Nakar, E., & Piran, T. 2019, Monthly Notices of the Royal Astronomical Society, Volume 488, Issue 2, p.2405-2411, 488, 2405, doi: 10.1093/mnras/stz1906
  * Gottlieb et al. (2018) Gottlieb, O., Nakar, E., Piran, T., & Hotokezaka, K. 2018, Monthly Notices of the Royal Astronomical Society, 479, 588, doi: 10.1093/mnras/sty1462
  * Harrison et al. (2018) Harrison, R., Gottlieb, O., & Nakar, E. 2018, Monthly Notices of the Royal Astronomical Society, Volume 477, Issue 2, p.2128-2140, 477, 2128, doi: 10.1093/mnras/sty760
  * He et al. (2018) He, H.-N., Kusenko, A., Nagataki, S., Fan, Y.-Z., & Wei, D.-M. 2018, ApJ, 856, 119, doi: 10.3847/1538-4357/aab360
  * Hillas (1984) Hillas, A. M. 1984, ARA&A, 22, 425, doi: 10.1146/annurev.aa.22.090184.002233
  * Hjorth et al. (2017) Hjorth, J., Levan, A. J., Tanvir, N. R., et al. 2017, ApJ, 848, L31, doi: 10.3847/2041-8213/aa9110
  * Hümmer et al. (2010) Hümmer, S., Rüger, M., Spanier, F., & Winter, W. 2010, ApJ, 721, 630, doi: 10.1088/0004-637X/721/1/630
  * Kashiyama et al. (2013) Kashiyama, K., Murase, K., Horiuchi, S., Gao, S., & Mészáros, P. 2013, ApJ, 769, L6, doi: 10.1088/2041-8205/769/1/L6
  * Kathirgamaraju et al. (2019) Kathirgamaraju, A., Tchekhovskoy, A., Giannios, D., & Barniol Duran, R. 2019, MNRAS, 484, L98, doi: 10.1093/mnrasl/slz012
  * Katz et al. (2011) Katz, B., Sapir, N., & Waxman, E. 2011, arXiv e-prints, arXiv:1106.1898.  https://arxiv.org/abs/1106.1898
  * Keshet & Waxman (2005) Keshet, U., & Waxman, E. 2005, Phys. Rev. Lett., 94, 111102, doi: 10.1103/PhysRevLett.94.111102
  * Kimura et al. (2018a) Kimura, S. S., Murase, K., Bartos, I., et al. 2018a, Phys. Rev. D, 98, 043020, doi: 10.1103/PhysRevD.98.043020
  * Kimura et al. (2018b) Kimura, S. S., Murase, K., & Mészáros, P. 2018b, ApJ, 866, 51, doi: 10.3847/1538-4357/aadc0a
  * Komissarov (2012) Komissarov, S. S. 2012, MNRAS, 422, 326, doi: 10.1111/j.1365-2966.2012.20609.x
  * Levinson & Bromberg (2008) Levinson, A., & Bromberg, O. 2008, Phys. Rev. Lett., 100, 131101, doi: 10.1103/PhysRevLett.100.131101
  * Levinson & Nakar (2019) Levinson, A., & Nakar, E. 2019, eprint arXiv:1909.10288.  https://arxiv.org/abs/1909.10288
  * Lipari et al. (2007) Lipari, P., Lusignoli, M., & Meloni, D. 2007, Phys. Rev. D, 75, 123005, doi: 10.1103/PhysRevD.75.123005
  * Liska et al. (2019) Liska, M., Chatterjee, K., Tchekhovskoy, A., et al. 2019, arXiv e-prints, arXiv:1912.10192.  https://arxiv.org/abs/1912.10192
  * Margutti & Chornock (2020) Margutti, R., & Chornock, R. 2020, arXiv e-prints, arXiv:2012.04810.  https://arxiv.org/abs/2012.04810
  * Matsumoto et al. (2021) Matsumoto, J., Komissarov, S. S., & Gourgouliatos, K. N. 2021, MNRAS, 503, 4918, doi: 10.1093/mnras/stab828
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, The Astrophysical Journal Supplement Series, 170, 228, doi: 10.1086/513316
  * Mimica & Aloy (2010) Mimica, P., & Aloy, M. A. 2010, MNRAS, 401, 525, doi: 10.1111/j.1365-2966.2009.15669.x
  * Mizuta & Ioka (2013) Mizuta, A., & Ioka, K. 2013, The Astrophysical Journal, Volume 777, Issue 2, article id. 162, 17 pp. (2013)., 777, doi: 10.1088/0004-637X/777/2/162
  * Mooley et al. (2018) Mooley, K. P., Deller, A. T., Gottlieb, O., et al. 2018, Nature, 561, 355, doi: 10.1038/s41586-018-0486-3
  * Murase & Ioka (2013) Murase, K., & Ioka, K. 2013, Phys. Rev. Lett., 111, 121102, doi: 10.1103/PhysRevLett.111.121102
  * Nakar (2015) Nakar, E. 2015, ApJ, 807, 172, doi: 10.1088/0004-637X/807/2/172
  * Nakar (2020) —. 2020, Phys. Rep., 886, 1, doi: 10.1016/j.physrep.2020.08.008
  * Nakar et al. (2018) Nakar, E., Gottlieb, O., Piran, T., Kasliwal, M. M., & Hallinan, G. 2018, The Astrophysical Journal, Volume 867, Issue 1, article id. 18, 9 pp. (2018)., 867, doi: 10.3847/1538-4357/aae205
  * Nakar & Sari (2012) Nakar, E., & Sari, R. 2012, ApJ, 747, 88, doi: 10.1088/0004-637X/747/2/88
  * Ostrowski (2000) Ostrowski, M. 2000, MNRAS, 312, 579, doi: 10.1046/j.1365-8711.2000.03146.x
  * Pol et al. (2019) Pol, N., McLaughlin, M., & Lorimer, D. R. 2019, The Astrophysical Journal, 870, 71
  * Rieger & Duffy (2004) Rieger, F. M., & Duffy, P. 2004, The Astrophysical Journal, 617, 155
  * Rieger & Duffy (2016) Rieger, F. M., & Duffy, P. 2016, ApJ, 833, 34, doi: 10.3847/1538-4357/833/1/34
  * Sahayanathan (2009) Sahayanathan, S. 2009, Monthly Notices of the Royal Astronomical Society: Letters, 398, L49
  * Senno et al. (2016) Senno, N., Murase, K., & Mészáros, P. 2016, Phys. Rev. D, 93, 083003, doi: 10.1103/PhysRevD.93.083003
  * Shahmoradi & Nemiroff (2015) Shahmoradi, A., & Nemiroff, R. J. 2015, MNRAS, 451, 126, doi: 10.1093/mnras/stv714
  * Sironi et al. (2013) Sironi, L., Spitkovsky, A., & Arons, J. 2013, ApJ, 771, 54, doi: 10.1088/0004-637X/771/1/54
  * van Eerten (2015) van Eerten, H. J. 2015, Journal of High Energy Astrophysics, 7, 23, doi: 10.1016/j.jheap.2015.04.004
  * Vietri (1995) Vietri, M. 1995, ApJ, 453, 883, doi: 10.1086/176448
  * Waxman & Bahcall (1997) Waxman, E., & Bahcall, J. 1997, Phys. Rev. Lett., 78, 2292, doi: 10.1103/PhysRevLett.78.2292
  * Waxman & Loeb (2001) Waxman, E., & Loeb, A. 2001, Phys. Rev. Lett., 87, 071101, doi: 10.1103/PhysRevLett.87.071101
  * Wei et al. (2016) Wei, J.-J., Wu, X.-F., Gao, H., & Mészáros, P. 2016, J. Cosmology Astropart. Phys, 2016, 031, doi: 10.1088/1475-7516/2016/08/031
  * Wu et al. (2019) Wu, M.-R., Banerjee, P., Metzger, B. D., et al. 2019, The Astrophysical Journal, 880, 23

##  Appendix A Numerical simulations setup

All simulations make use of a relativistic ideal gas equation of state,
piecewise parabolic reconstruction method and an HLL Riemann solver. In the
RMHD simulations, the jet injection height is $z_{\rm{beg}}$, and thus its
injection nozzle radius is $r_{\rm noz}=z_{\rm{beg}}\theta_{j}$, where the
origin is the center of the system of the merger. The weakly magnetized jets
are injected with a toroidal magnetic field whose peak magnitude at $r=r_{\rm
noz}/2$ is set by $\sigma_{0}=10^{-2}h_{j}$, where $h_{j}$ is the jet’s
initial specific enthalpy (see Gottlieb et al., 2020a, for the full profile).
The jet quantities of each model are listed in Table 2. We employ 3D Cartesian
grids to properly account for the formation of instabilities and avoid
numerical artifacts (see e.g., Gottlieb et al., 2021b). After the jet reaches
the homologous expansion phase, the hydrodynamic jets are remapped into
axisymmetric grid for modeling the afterglow phase (see Nakar et al., 2018,
for details). The grid setups of the successful hydrodynamical and weakly
magnetized models are provided in Mooley et al. (2018) and Gottlieb et al.
(2020a) respectively101010The hydrodynamic choked jets has an identical grid
setup to the successful jets, and the weakly magnetized choked jet has an
identical setup to the successful weakly magnetized jet with two differences:
the jet launching duration $t_{j}$, and a higher grid resolution by 25% on
$\hat{x}$\- and $\hat{y}$\- axes and by 60% in $\hat{z}$-axis..

In the GRMHD simulation the vector-potential profile is

$A=A_{\phi}(r,\theta)=10^{13}\frac{{\rm sin}\theta}{r}\cdot{\rm
max}\bigg{(}\frac{r^{2}}{r^{3}+r_{c}^{3}}-\frac{r_{o}^{2}}{r_{o}^{3}+r_{c}^{3}},0\bigg{)}~{}{\rm
Mx~{}cm^{-1}}~{},$ (A1)

where $r_{c}=5\times 10^{7}$ cm and $r_{o}=10^{10}~{}{\rm cm}$. The second
term sets $B=0$ at $r\approx 15r_{g}$ and at $r=r_{o}$. The jet is generated
self-consistently, its resulting parameters are listed in Table 2. We use a
tilted metric to avoid numerical dissipation on the jet axis at $\theta=0$.
The GRMHD simulation is carried out in spherical coordinates with logarithmic
grid on $\hat{r}$-axis from $r_{g}$ to $5\times 10^{9}\,{{\rm cm}}$, and
uniform on $\hat{\theta}$\- and $\hat{\phi}-$axes . We use two levels of AMR,
which is activated at high entropy to properly resolve the jet and the cocoon.
Our effective resolution in the maximal refinement level is $1152\times
576\times 512$ cells on $\hat{r}-\hat{\theta}-\hat{\phi}$ respectively, enough
to resolve the development of magneto-rotational instabilities in the disk.

Model | $\sigma_{0}$ | Jet fate |  $t_{d}$ [s] | $\theta_{j}$ | $h_{j}$ | $z_{\rm{beg}}$ | $L_{j}[10^{50}~{}{\rm erg~{}s^{-1}}]$ |  $t_{j}$ [s] | $\Gamma_{\infty}$
---|---|---|---|---|---|---|---|---|---
$\,{\it SU1}$ | 0 | Successful | 0.2 | $4^{\circ}$ | 200 | $4.5\times 10^{8}~{}{\rm cm}$ | 1.4 | 0.8 | 2000
$\,{\it SU2}$ | 0 | Successful | 0.4 | $15^{\circ}$ | 100 | $3\times 10^{8}~{}{\rm cm}$ | 2 | 2.0 | 300
$\,{\it SU3}$ | 0 | Successful | 0.6 | $8^{\circ}$ | 100 | $7\times 10^{8}~{}{\rm cm}$ | 0.3 | 3.6 | 500
$\,{\it CU}$ | 0 | Choked | 0.2 | $23^{\circ}$ | 20 | $4.5\times 10^{8}~{}{\rm cm}$ | 3 | 1.0 | 32
$\,{\it WM}$ | $10^{-2}h_{j}$ | Successful | 0.6 | $8^{\circ}$ | 100 | $7.5\times 10^{8}~{}{\rm cm}$ | 0.2 | 3.0 | 500
$\,{\it CW}$ | $10^{-2}h_{j}$ | Choked | 0.6 | $8^{\circ}$ | 100 | $7.5\times 10^{8}~{}{\rm cm}$ | 0.2 | 0.225 | 500
$\,{\it HM}$ | $\sim 25$ | Successful | 0.5 | $\sim 12^{\circ}$ | $\sim 1$ | $r_{h}$ | $\sim 1$ | 0.4 | $\sim 100$
Table 2: Parameters of the different models. In the hydrodynamic and weakly
magnetized jets the parameters are set as boundary conditions whereas in the
highly magnetized jets the jet is launched self-consistently and the shown
parameters are average values as found in the simulation. Shown are the
initial jet magnetization $\sigma_{0}$; the specific jet enthalpy $h_{j}$; the
time delay from the merger to the jet launching, $t_{d}$, during which the
ejecta expanded homologously; the jet opening angle upon injection,
$\theta_{j}$; the jet injection height in the simulation, $z_{\rm{beg}}$,
where $r_{h}=r_{g}(1+\sqrt{1-a^{2}})$ is the event horizon; the jet total
luminosity, $L_{j}$; the engine working time, $t_{j}$; the terminal Lorentz
factor, if no mixing takes place, $\Gamma_{\infty}$.

##  Appendix B Calibration of the shock breakout energy

Once the jet-cocoon forward shock breaks out from the ejecta into an optically
thin medium, it becomes collisionless, emit $\gamma$-ray photons, and
accelerates CR which interact with the photons in the breakout layer to
generate neutrinos. The shock breakout of the mildly relativistic cocoon is a
leading candidate for the origin of the $\gamma$-ray signal in GW170817, as it
was found to be consistent with all the $\gamma-$ray signal observables
(Gottlieb et al., 2018). Therefore, we use the properties of GRB170817A to
estimate the breakout layer characteristics, assuming that cocoon shock
breakout is indeed the physical mechanism behind the $\gamma$-ray production
and that GW170817 is a typical BNS. The $\gamma$-ray isotropic equivalent
energy of GRB170817A is $E_{\rm iso,{\rm GRB17}}\approx 5\times 10^{46}$ erg
(Goldstein et al., 2017), and the breakout layer was at $R_{\rm bo}\approx
2\times 10^{11}$ cm (Gottlieb et al., 2018). To find the energy and velocity
angular structures of the jet-cocoon system, we utilize RMHD simulations of
both successful and choked jets. For the successful jet we average the
breakout properties of two high resolution simulations: $II$ in Gottlieb &
Loeb (2020) (our model $\,{\it SU2}$), and $S_{b}$ in Gottlieb et al. (2021b)
(our model $\,{\it SU3}$), and for the uncollimated choked jet we use
simulation $\,{\it CU}$ (see Appendix A for the initial parameters of the
models).

The energy originates in the range of angles between the minimal angle in
which the Lorentz factor satisfies $\Gamma^{-1}>\theta_{\rm obs}-\theta$, and
$\sim\theta_{\rm obs}$. We calibrate the energy in this range in the
simulations with the observed energy, and then calculate the energy observed
by observers at other viewing angles. Our simulations suggest that the width
of $\tau$ of a few in the merger frame is of a few percent of the shock
radius, such that the width in the shock frame is $\sim 10^{10}\Gamma^{-1}$
cm. in agreement with analytic estimates (e.g., Nakar & Sari, 2012). At
$R_{\rm bo}$ the magnetic field becomes subdominant. By assuming adiabatic
expansion of our weakly magnetized jet simulation, we find that the mass
density is $3\times 10^{-8}~{}{\rm g~{}cm^{-3}}$, and the magnetic field is
$B(R_{\rm bo})\approx 3\times 10^{6}$ G. The latter is also obtained by
applying the equipartition parameter $\epsilon_{B}=10^{-3}$ on the shock.

Finally, we estimate the $\gamma$-ray photon spectrum to be composed of two
phases of equal amount of energy. First, energy from the breakout layer is
released without expansion of the breakout layer or thermalization of the
photons, thus we use a Wien spectrum with a rest-frame peak at 50 keV, as
appropriate for mildly relativistic radiation-mediated shocks. Over time
emission from deeper layers is leaking and there is more time for the photons
to thermalize and cool to a blackbody spectrum around 1 keV. For our
accelerated CR energies, interaction with both $\gamma$-ray and 1 keV X-ray
photons can yield neutrinos.

##  Appendix C Neutrinos from choked jets

Choked jets may produce neutrinos from subphotospheric shocks before or after
they are choked. Here we show that after the jet is choked, cocoon shock
breakout is the only neutrino production site (§4.3), and no neutrino emission
is expected while the cocoon is still inside the ejecta. The reason lies in
heavy mixing of jet and cocoon material once the jet engine is shut off
(Gottlieb et al., 2021a). The heavy mixing reduces the magnetization, thereby
disfavoring the formation of collisionless shocks. A relativistic jet can be
choked due to a variety of reasons. Generally, choked jets can be sub-
categorized to collimated and uncollimated choked jets. In case of the former,
the jet engine is terminated before the jet finds its way outside the ejecta.

When the injected jet is uncollimated, the jet cannot punch through the dense
medium, even if the central engine operates for a typical GRB duration or
longer. In that case, the angular distribution of the outflow is more
isotropic. A collimated jet that is choked just before breakout has less time
to mix with the cocoon material so it can keep its magnetization high,
compared to a jet that is choked early-on or an uncollimated jet. Thus, we
perform a numerical simulation of such an optimal case (model $\,{\it CW}$),
by fine tuning the collimated jet to be choked just upon breakout, we
demonstrate that even under such conditions the mixing in the jet is too high
to allow the formation of a collisionless subshock. Fig. 7 depicts
$\hat{x}-\hat{z}$ meridional plane maps of the choked jet upon breakout. The
jet cavity is filled with baryonic matter while the magnetized cocoon
undergoes mixing with the surrounding unmagnetized ejecta such that its
magnetization drops. With characteristic $\sigma\sim 10^{-2.5}$ in the cocoon,
it is unable to form collisionless subshocks in the outflow. We conclude that
even the most optimistic choked jet model, particle acceleration is
inefficient. If the jet is choked earlier or is uncollimated in the first
place, the resulting cocoon would maintain an even lower magnetization.

Figure 7: Logarithmic mass density (left) and magnetization (right) meridional
maps of a jet that is choked upon breakout from the ejecta.