# High-Energy Neutrinos from Gamma-Ray-Faint Accretion-Powered Hypernebulae

Navin Sridhar Department of Astronomy, Columbia University, New York, NY
10027, USA  Theoretical High Energy Astrophysics (THEA) Group, Columbia
University, New York, New York 10027, USA Navin Sridhar
<EMAIL_ADDRESS>Brian D. Metzger Department of Physics, Columbia
University, New York, NY 10027, USA  Theoretical High Energy Astrophysics
(THEA) Group, Columbia University, New York, New York 10027, USA  Center for
Computational Astrophysics, Flatiron Institute, New York, New York 10010, USA
Ke Fang Department of Physics, Wisconsin IceCube Particle Astrophysics Center,
University of Wisconsin, Madison, WI, 53706, USA

###### Abstract

Hypernebulae are inflated by accretion-powered winds accompanying hyper-
Eddington mass transfer from an evolved post-main sequence star onto a compact
object. The ions accelerated at the termination shock of the
hypernebula—formed by the interaction of the collimated fast disk winds/jet
against the slower, wide-angled winds—can engage in neutrino generation via
hadronic ($pp$) reactions, and photomesonic ($p\gamma$) interactions with the
disk thermal and Comptonized nonthermal background photons. If hypernebulae
follow the volumetric rates of steller-merger transients such as common
envelope events, the volume-integrated neutrino emission from them—depending
on the mass-transfer rate—could explain $\sim 25\%-100\%$ of the high-energy
diffuse neutrino flux observed by the IceCube Observatory and the Baikal-GVD
telescope. Notably, the time-averaged neutrino spectra from
hypernebulae—depending on the population parameters—could also provide a
reasonable description of the observed neutrino spectral shape. Additionally,
some hypernebulae could be sources of $>$100 PeV ultra-high-energy neutrinos
that could be detected by upcoming neutrino observatories. The large optical
depth of hypernebulae to Breit-Wheeler ($\gamma\gamma$) interaction attenuates
the escape of GeV-PeV gamma-rays produced along with neutrinos, which makes
them gamma-ray-faint sources of neutrinos, thus being consistent with the
Fermi observations of the isotropic gamma-ray background. Furthermore, a
population of fast radio bursts (FRBs) could also be powered by jetted hyper-
accreting engines. Even though the neutrino emission from the ms-duration
bursts is less likely to be detected, the FRB sources associated with a
persistent radio counterpart could be significant contributors to the high
energy diffuse background neutrino flux if they are powered by hypernebulae.

Neutrino Astronomy (1100); Radio transient sources (2008); Ultraluminous X-ray
sources (2164); X-ray binary stars (1811); Shocks (2086); Plasma astrophysics
(1261); High energy astrophysics (739)

##  1 Introduction

IceCube, and recently, the Baikal-GVD neutrino telescopes measured the flux
and spectrum of the diffuse high energy (TeV-100 PeV; HE, hereafter) neutrino
background, and placed measured the all-flavor all-sky 100 TeV flux at
$5\times 10^{-8}\,{\rm GeV\,s^{-1}\,cm^{-2}\,sr^{-1}}$ (Aartsen+18;
Aartsen+20a; Abbasi+21; Abbasi+22a; Baikal_Collaboration22). The sensitivity
limits of the current neutrino facilities have not allowed for the detection
of the long sought-after ultra-high-energy ($\geq$100 PeV; UHE, hereafter).
The detected HE neutrinos are believed to be produced in baryonic shocks: one
of the consequences of magnetized shocks of electron-ion plasmas is the non-
thermal acceleration of ions to relativistic energies via diffusive shock
acceleration (Blandford&Ostriker78; and also via turbulence, for example,
Comisso&Sironi22). Such relativistic ions can generate neutrinos via hardonic
($pp$), and photomesonic ($p\gamma$) interactions (Waxman&Bahcall97;
Meszaros&Waxman01; Dermer&Atoyan03; Guetta&Granot03). The contribution from
the Galactic plane to the all-sky neutrino flux is found to be
$\lesssim$5%-10% at around 100 TeV suggesting that most of the observed
apparently isotropic diffuse neutrino flux is extragalactic (Fang&Murase21).
There has been evidence for at least one gamma-ray blazar, TXS 0506+056,
contributing to the neutrino signal (IceCube_Blazar1+18; IceCube_Blazar2+18).
Recently, the ten-year point-source searches with IceCube have indicated that
the type-II Seyfert starburst galaxy NGC 1068 is the most significant (at
$4.2\sigma$) steady source of neutrinos (Aartsen+20; IceCubeCollaboration_22).
However, in their respective energy ranges, NGC 1068 and TXS 0506+056
contribute no more than $\sim$1% above 10 TeV to the overall HE diffuse
background neutrino flux of astrophysical neutrinos.

Therefore, the sources for the majority of the extragalactic background
neutrinos still remain a mystery. Various astrophysical objects have been
studied as the sources of the HE neutrinos, such as gamma-ray bursts (GRBs;
Paczynski&Hu94; Waxman&Bahcall97; Meszaros&Waxman01; Dermer&Atoyan03;
Fasano+21; Kimura_22), active galactic nuclei (AGN; Eichler79; Stecker+91;
Jacobsen+15; Murase17; Kun+22), supernovae (SNe; Razzaque+04; Murase+11;
Murase18; Chang+22), starburst galaxies (Romero&Torres03; Loeb&Waxman06;
Thompson+06), tidal disruption events (TDEs; Wang&Liu16; Dai&Fang17; Senno+17;
Lunardini&Winter17), nonrelativistic shock-powered transients (such as e.g.,
luminous red novae, fast blue optical transients; Fang+20), fast radio bursts
(FRBs; Metzger+20), protomagnetars (Bhattacharya+22), supermassive black hole
mergers (Jaroschewski+22), and so on. However, stacking analysis with IceCube
and ANTARES has shown that none of the above sources contribute a major
fraction of the all-sky neutrino flux (Aartsen+16; Albert+21; Liodakis+22).

The HE diffuse background neutrino flux is comparable to that of UHE cosmic
rays and gamma-rays, hinting at the possibility of a common origin. This is
expected because a comparable flux of gamma-rays will also be produced with
neutrinos following both the proton-proton ($pp$) and proton-photon
($p\gamma$) interactions. But, the diffuse isotropic gamma-ray background
(IGRB; subtracting the contributions from resolved sources, foreground
emission, and unresolved blazars) observed by Fermi Large Area Telescope (LAT)
between 100 MeV and 1 TeV (Ackermann+15) has disfavored gamma-ray bright
sources being dominant sources of HE neutrinos, i.e., at lower energies, there
is more neutrino flux than the corresponding gamma-ray flux if they were
produced by the same sources through the $pp$ or $p\gamma$ channels
(IceCube21). This fact that the all-sky neutrino flux at 10 TeV is larger than
the expected gamma-ray flux (from an accelerated proton with a powerlaw index
of $-2$) suggests that the dominant source of the HE diffuse background
neutrino flux must be opaque to GeV-TeV gamma-rays (Murase+16; Capanema+20;
Capanema+21; Fang+22).

Such requirements could naturally be met by the new class of astrophysical
objects called ‘hypernebulae’ (Sridhar&Metzger22). A hypernebula is an
energetic, compact “bubble” ($\lesssim 1$ pc; much smaller and powerful than
the nebulae surrounding typical ultraluminous X-ray (ULX) sources) of plasma
inflated into the circumbinary medium through shocks by powerful accretion
disk winds and jets during the dynamically unstable (or even thermal-
timescale) transfer of matter into a compact object at rates that exceed the
Eddington rate by many orders of magnitude. Such brief phases of rapid mass
transfer are expected to occur following the Roche-lobe overflow from a donor
star, prior to a common-envelope merger event (and the transient events
accompanying it such as Luminous Red Novae). Recently Sridhar+21b showed that
the flaring collimated jets launched by the compact object at the heart of
such accreting engines could power the (periodically) repeating FRBs. The
radio-bright hypernebulae surrounding them could be the source of the
persistent radio emission seen spatially coincident with FRBs (Chatterjee+17;
Niu+22), which also imparts the rotation measure onto the pulses propagating
through them (Sridhar&Metzger22). In fact, hypernebulae could already be
present in blind wide-field radio surveys such as FIRST (Becker+95) and VLASS
(Lacy+20) as off-nuclear radio sources, and could be spatially resolved by the
next generation VLA (Selina+18).

In this paper, we investigate the properties of the ions accelerated at
different shock fronts of hypernebula, and their interactions with other ions
($pp$) and the ‘background’ disk thermal and Comptonized nonthermal UV/X-ray
photons ($p\gamma$). Such interactions could produce pions which will decay
into neutrinos and gamma-rays. We show that the HE neutrinos generated by
hypernebulae could explain—depending on the mass-transfer rate—a significant
fraction ($\sim 25-100\%$) of the observed extragalactic HE diffuse neutrino
background. The dense windy environment of hypernebulae further attenuates the
$>$TeV gamma-rays produced by pion decay, thus making them gamma-ray faint. In
§2, we discuss the nature of the outflows from the accreting engine:
properties of the shocks and their associated timescales (§2.1), and the
sources of background radiation (§2.2). The model for neutrino production is
presented in §3, with discussions on pion production channels (§3.1), pion
production efficiency (§3.2), ion heating and acceleration (§3.3), and
neutrino production from pion decay (§3.4). We then present the observational
implications of our model in §4 before concluding in §5.

##  2 Properties of the engine

As discussed in Sridhar&Metzger22, there can be at least three different types
of disk outflow interactions in a hypernebula. (1) The slow, wide-angled winds
from the accretion disk can interact with the circumstellar medium (CSM) via a
sub-relativistic forward shock (FS); the location of this shell—that sweeps up
both the wind material and the CSM—controls the overall radius of the system,
and hence, the gas and X-ray photon number density within the hypernebula. (2)
The trailing slow winds hit the accumulated shell ahead of them, forming what
we call, the ‘wind termination shock’ (WTS). (3) The faster winds/jet
channeled along the instantaneous angular momentum axis of the inner accretion
disk/jet can inflate bipolar lobes upon their interaction with the swept-up
shell, mediated by what we call, the ‘jet termination shock’ (JTS). We refer
readers to Fig. 1 of Sridhar&Metzger22 for a schematic diagram of the
hypernebula and its various interaction fronts.

Our assumed fiducial model parameters of the system [and the explored range
inside square brackets] are as follows. The central compact object is assumed
to be a black hole (although an accreting neutron star would work too) with
mass $m_{\bullet}=M_{\bullet}/M_{\odot}=10$, accreting at a rate of
$\dot{m}=\dot{M}/\dot{M}_{\rm Edd}=10^{5}~{}[10^{3},10^{7}]$ 111$\dot{M}_{\rm
Edd}=L_{\rm Edd}/(0.1c^{2})$ is the Eddington accretion rate, and $L_{\rm
Edd}\simeq 1.5\times 10^{39}(m_{\bullet}/10)\,{\rm erg\,s^{-1}}$ from a
companion post-main sequence star with mass
$m_{\star}=M_{\star}/M_{\odot}=30$. The accretion disk drives wide-angled slow
winds with speeds $v_{\rm w}/c=0.03~{}[0.01,0.1]$ and kinetic luminosity
$L_{\rm w}\sim\dot{M}_{\rm w}v_{\rm w}^{2}$, and drives fast wind/jet about
the black hole spin axis with a speed $v_{\rm j}/c=0.3~{}[0.1,0.5]$
(Kobayashi+18; Urquhart+18; Pinto&Kosec22). We take the mass loss rate in the
wide angled disk winds to be $\dot{M}_{\rm w}\sim\dot{m}$ (Hashizume+15), and
the mass loss rate in the jet to be $\dot{M}_{\rm j}=\eta\dot{M}_{\rm
w}(v_{\rm w}/v_{\rm j})^{2}$, with $\eta\sim 0.5~{}[0.1,1]$ expected from the
advection dominated models of accretion disk, due to the roughly equal
gravitational energy released per radial decade in the accretion flow
(Blandford&Begelman99). We define the magnetization parameter of the jet
$\sigma_{\rm j}=0.1~{}[0.01,1]$ as the ratio of its Poynting flux to kinetic
energy flux. The density of the ‘upstream’ CSM with which the outflows
interact is taken to be $\rho_{\rm csm}=\mu nm_{\rm p}$, where $\mu=1.38$ is
the mean atomic weight for neutral solar composition gas, and
$n=10\,[1,100]\,{\rm cm}^{-3}$, typical of the environments of massive stars
(Abolmasov+08; Pakull+10; Toala&Arthur_2011).

###  2.1 Ejecta

Figure 1: Top panel: Different timescales of the system for our fiducial model
with $\dot{m}=10^{5}$ (solid curves); dotted and dashed curves represent
$\dot{m}=10^{3}$ and $\dot{m}=10^{7}$, respectively. Pink curves denote the
hypernebula expansion timescale (Eq. 5); yellow curves denote the light-
crossing timescale of the hypernebula ($t_{\rm cross}=R/c$), and the black
curves denote the time particles take to laterally cross the X-ray emission
cone (Eq. 32). Dark and light blue curves denote the interaction timescale of
protons between thermal and nonthermal photons (Eq. 3.2), respectively. Dark
brown curves denote the hadronic ($pp$) interaction timescale (Eq. 3.2). The
overall pion creation timescale (Eq. 35), not shown here, overlaps with
$t_{\rm p\gamma}^{\rm th}$ at earlier times, and with $t_{\rm pp}$ at later
times, when applicable (e.g., for $\dot{m}=10^{3}$). Bottom panel: Neutrino
production suppression factor for different accretion rates for the fiducial
model (Eq. 37).

The mass transfer from a donor star that has evolved off the main sequence,
preceding a stellar merger/common envelope event could be dynamically
unstable. The active duration of such a rapid, unstable mass-transferring
phase can be approximated as,

$t_{\rm active}\sim\frac{M_{\star}}{\dot{M}}\approx 10^{3}\,{\rm
yr}\left(\frac{m_{\star}}{30}\right)\dot{m}_{5}^{-1},$ (1)

Roughly, this corresponds to the lifetime of a hypernebula. Throughout the
paper, we adopt the shorthand notation $Y_{\rm x}\equiv Y/10^{\rm x}$ for
quantities given in cgs units, and in in units of $\dot{M}_{\rm Edd}$ for
$\dot{m}$. The ejecta freely expands until a time when the mass released in
the accretion disk-winds is comparable to the mass accumulated from the CSM by
the ejecta. The hypernebula then enters a Sedov-Taylor (Sedov_59) deceleration
phase after a time (Sridhar&Metzger22),

$t_{\rm free}\approx 80\,{\rm yr}\,\left(\frac{L_{\rm
w,42}}{n_{1}}\right)^{1/2}v_{\rm w,9}^{-5/2}.$ (2)

Unless otherwise specified, we fiducially normalize all our timescales and
analytical calculations to $t=t_{\rm free}$ throughout the rest of the paper.
The temporal evolution of the radius of the shell—bridging the two phases
across $t_{\rm free}$—can be approximated as,

$R(t)\simeq\begin{cases}v_{\rm w}t\approx 0.75\,{\rm pc}\,v_{\rm
w,9}\left(\frac{t}{80\,{\rm yr}}\right)&(t<t_{\rm free})\\\
\alpha\left(\frac{L_{\rm w}t^{3}}{\rho_{\rm csm}}\right)^{1/5}\approx
0.75\,{\rm pc}\,\left(\frac{L_{\rm
w,42}}{n_{1}}\right)^{1/5}\left(\frac{t}{80\,{\rm yr}}\right)^{3/5}&(t>t_{\rm
free}),\\\ \end{cases}$ (3)

where $\alpha\approx 0.88$ while the forward shock is adiabatic and
$\alpha\approx 0.76$ after it becomes radiative (Weaver+77)222As we show
below, for the cases of interest $\dot{m}>10^{3}$, the FS does not become
radiative by the time the accreting engine stops.. This shell expands at a
rate set by the velocity of the forward shock,

$v_{\rm fs}=\frac{dR}{dt}\approx 10^{9}\,{\rm
cm\,s}^{-1}\times\begin{cases}v_{\rm w,9}&(t<t_{\rm free})\\\
\left(\frac{L_{\rm w,42}}{n_{1}}\right)^{1/5}\left(\frac{t}{80\,{\rm
yr}}\right)^{-2/5}&(t>t_{\rm free}).\\\ \end{cases}$ (4)

The corresponding expansion timescale can be defined as,

$t_{\rm exp}\simeq\frac{R}{v_{\rm fs}}\simeq 80\,{\rm
yr}\left(\frac{t}{80\,{\rm yr}}\right),$ (5)

over which adiabatic losses occur. As the shell expands, the mass it
accumulates within is,

$M_{\rm sh}\simeq\begin{cases}\dot{M}_{\rm w}t\approx
2\,M_{\odot}\dot{m}_{5}\left(\frac{t}{80\,{\rm yr}}\right)&(t<t_{\rm free})\\\
(4\pi/3)\rho_{\rm csm}R^{3}\approx 2\,M_{\odot}L_{\rm
w,42}^{3/5}n_{1}^{2/5}\left(\frac{t}{\rm 80\,yr}\right)^{9/5}&(t>t_{\rm
free}).\end{cases}$ (6)

At later times ($t>t_{\rm free}$), the swept-up mass from the CSM grows
rapidly as $t^{9/5}$ compared to the mass of the initial disk wind. With this,
we can calculate the baryon number density inside the shell at different
times,

$\displaystyle n_{\rm p}\simeq\frac{3M_{\rm sh}}{4\pi R^{2}m_{\rm p}\mu{\rm
d}R}\approx\begin{cases}30\,{\rm cm}^{-3}\,\dot{m}_{5}v_{\rm
w}^{-3}\left(\frac{t}{80\,{\rm yr}}\right)^{-2}&(t<t_{\rm free})\\\ 10\,{\rm
cm}^{-3}\,n_{1}&(t>t_{\rm free})\\\ \end{cases},$ (7)

where $dR\sim R$ for adiabatic shocks, as it is in our case. This number
density is substantially lower than say, a magnetar wind nebula
(Fang&Metzger17), during its active neutrino production phase (which lasts for
$\lesssim 1$ yr since its birth), and has implications to the interaction
channels through which neutrinos are produced (more in §3.2).

The kinetic power in the outflowing ejecta is distributed in the following way
between different shock fronts,

$L_{\rm s}\approx\frac{9\pi}{8}R^{2}\times\begin{cases}\rho_{\rm fs}v_{\rm
fs}^{3}&({\rm FS})\\\ \rho_{\rm wts}(v_{\rm w}-v_{\rm fs})^{3}&({\rm WTS})\\\
\rho_{\rm jts}(v_{\rm j}-v_{\rm fs})^{3}&({\rm JTS})\end{cases}$ (8)

where the densities of the upstream gas entering the above expressions are
$\rho_{\rm fs}=\rho_{\rm csm}$, $\rho_{\rm wts}=\dot{M}_{\rm w}/(4\pi v_{\rm
w}R^{2})$, and $\rho_{\rm jts}=\dot{M}_{\rm j}/(4\pi v_{\rm j}R^{2})$,
respectively. At $t>t_{\rm free}$, the gas behind the shocks gets heated to a
mean energy,

$\bar{E}\approx\frac{3}{16}\mu m_{\rm p}\times\begin{cases}v_{\rm
fs}^{2}\approx 220\,{\rm keV}\,\left(\frac{L_{\rm
w,42}}{n_{1}}\right)^{2/5}\left(\frac{t}{\rm 80\,yr}\right)^{-4/5}&({\rm
FS})\\\ (v_{\rm w}-v_{\rm fs})^{2}\approx 120\,{\rm keV}\,v_{\rm
w,9}^{2}&({\rm WTS})\\\ (v_{\rm j}-v_{\rm fs})^{2}\approx 17\,{\rm
MeV}\,\left(\frac{v_{\rm j}}{0.3\,c}\right)^{2}&({\rm JTS})\end{cases}$ (9)

The corresponding cooling timescale is (e.g., Vlasov+16),

$t_{\rm cool}\simeq\frac{3\bar{E}}{2n_{\rm s}\mu\Lambda[T_{\rm s}]}$ (10)

where $n_{\rm s}=4\rho_{\rm s}/\mu m_{\rm p}$ is the post-shock density for
adiabatic index $\hat{\gamma}=5/3$, and we have assumed a cooling function
appropriate for solar-metallicity gas that includes both free-free emission
and lines $\Lambda[T_{\rm s}]=\Lambda_{\rm ff}[T_{\rm s}]+\Lambda_{\rm
line}[T_{\rm s}]$, where (Schure+09; Draine11),

$\displaystyle\Lambda_{\rm ff}[T_{\rm s}]\simeq 3\times 10^{-27}(T_{\rm
s}/{\rm K})^{0.5}\,{\rm erg~{}s}^{-1}{\rm cm}^{3},$ (11a)
$\displaystyle\Lambda_{\rm lines}[T_{\rm s}]\simeq 2.8\times 10^{-18}(T_{\rm
s}/{\rm K})^{-0.7}\,{\rm erg~{}s}^{-1}{\rm cm}^{3}.$ (11b)

The radiative efficiency of the shocks can be quantified with a parameter
(Vlasov+16; Sridhar&Metzger22),

$f_{\rm rad}=\left(1+\frac{5}{2}\frac{t_{\rm cool}}{t_{\rm exp}}\right)^{-1},$
(12)

where $f_{\rm rad}<1$ implies an adiabatic expansion and $f_{\rm rad}>1$
implies that the shock is radiative.

Figure 2: Solid, dashed, and dotted curves show the properties of the forward
shock, wind termination shock, and jet termination shock, respectively. The
colors denote different accretion rates. The time evolution of the shock power
(Eq. 8) is shown in the top panel, and the bottom panel shows the shocks’
radiative efficiency (Eq. 12). The left-facing triangles at the bottom frame
denote the active lifetime of the accreting engine (Eq. 1), and the downward-
facing triangles at the top frame denote the free-expansion timescale of the
shock (Eq. 2).

The top panel of Fig. 2 shows the time evolution of the power in the three
shocks, for different $\dot{m}$. For the purpose of calculating the temporal
evolution of various shock properties, the engine is smoothly ‘turned-off’ at
a time $t=t_{\rm active}$ when accretion stops (Eq. 1), by multiplying the
shock properties with an exponential function of the form $e^{-t/t_{\rm
active}}$. For a given $\dot{m}$, at earlier times of the expansion ($t\ll
t_{\rm free}$), the power in the JTS dominates, which stays nearly constant
(as $v_{\rm j}\gg v_{\rm w}\gtrsim v_{\rm fs}$) till the engine shuts down.
Once the shell crosses the free-expansion phase at $t>t_{\rm free}$ (only for
$\dot{m}<10^{6}$), the higher $\rho_{\rm fs}$ ($>\rho_{\rm wts}\gg\rho_{\rm
jts}$) makes the kinetic power in the FS and WTS comparable to that of the
JTS. The bottom panel of Fig. 2 shows the radiative efficiency ($f_{\rm rad}$)
of FS and WTS. $f_{\rm rad}$ of JTS is not shown because the particles are
heated to relativistic temperatures behind it, and the corresponding extremely
long cooling timescales ($\gtrsim 10^{10}$ yr) means that they always stay
adiabatic $f_{\rm rad}\ll 10^{-8}$. Generally, the $f_{\rm rad}$ of the FS and
the WTS increases with $\dot{m}$. Unlike the FS, which turns from being
adiabatic to radiative at a critical time $t=t_{\rm rad}^{\rm fs}$, the WTS
turns from being radiative to adiabatic as soon as its onset. Also note that
the FS of only those cases with $\dot{m}\lesssim 10^{4}$ turns radiative
before $t_{\rm active}$.

For the rest of our calculations, we just concentrate on the JTS, for the
following reason. Even though the FS cumulatively carries more energy when
integrated over time and volume, the energy contained within the solid angle
relevant for photohadronic interactions (see Eq. 17) would make it (and the
WTS) sub-dominant compared to the JTS, where on the other hand, all of its
power is channeled along the relevant solid angle.

The magnetic field strength at the JTS due to the flux carried from the base
of the jet is,

$\displaystyle B_{\rm jts}$ $\displaystyle=\sqrt{16\pi\sigma_{\rm j}m_{\rm
p}n_{\rm s}v_{\rm j}^{2}}\approx 4\,{\rm mG}\,\sigma_{\rm
j,-1}^{1/2}\left(\frac{v_{\rm j}}{0.3\,c}\right)^{1/2}\eta_{-1}^{1/2}$
$\displaystyle\times\begin{cases}\dot{m}_{5}^{1/2}\left(\frac{t}{\rm
80\,yr}\right)^{-1}&(t<t_{\rm free})\\\ \dot{m}_{5}^{3/10}v_{\rm
w,9}^{3/5}n_{1}^{1/5}\left(\frac{t}{\rm 80\,yr}\right)^{-3/5}&(t>t_{\rm
free})\end{cases},$ (13)

where $n_{\rm s}=4\rho_{\rm jts}/\mu m_{\rm p}$ is the post-JTS particle
number density. In addition, the charged particles accelerated behind the JTS
play a role in the non-resonant streaming instability (Bell04) that amplifies
the magnetic field such that it saturates with an energy density,

$\frac{B_{\rm Bell}^{2}}{8\pi}\sim P_{\rm acc}\left(\frac{3v_{\rm
j}}{2c}\right),$ (14)

where the ram pressure of the accelerated particles at the JTS, $P_{\rm
acc}\simeq\epsilon_{\rm rel}\rho_{\rm jts}v_{\rm j}^{2}$ is assumed to be
$\epsilon_{\rm rel}=1\%$ of the JTS ram pressure (see also Eq. 39). $B_{\rm
jts}\gg B_{\rm Bell}$ for $\sigma_{\rm j}>10^{-3}$, and we take the total
magnetic field strength near the shock to be $B_{\rm sh}=B_{\rm Bell}+B_{\rm
jts}$.

###  2.2 Background radiation

There are at least three lasting sources of photon fields333We do not consider
thermal synchrotron radiation that peaks in the radio band (Sridhar&Metzger22)
because, the required proton threshold energy of $E_{\rm thr}\sim 10$ ZeV to
interact with $\sim$GHz photons and produce pions (Eq. 26) is much larger than
the maximum energy to which protons can be accelerated at the hypernebula JTS
(Fig. 4). with which the ions accelerated at the JTS (more in §3.3) can
interact to produce neutrinos.

(1) Disk thermal emission: For accretion rates $\dot{m}\gg 1$, the strong
outflows from the accretion disk can limit the bolometric luminosity of the
blackbody emission to (Begelman+06; Poutanen+07),

$L_{\rm bol}=L_{\rm Edd}(1+\ln\dot{m})\approx 10^{40}\,{\rm
erg\,s}^{-1}(1+\ln\dot{m}_{5}).$ (15)

The isotropic-equivalent luminosity of the disk emission along the axis of the
inner accreting funnel, however, can be much greater than this value due to
the accreting funnel’s geometric collimation. For an X-ray emitting cone with
a collimation angle of $\theta_{\rm X}$, the corresponding beaming fraction
can be defined as,

$f_{\rm b,X}\equiv 2\pi[1-\cos{(\theta_{\rm X})}]/4\pi\simeq\theta_{\rm
X}^{2}/4.$ (16)

For our considered case with $\dot{m}\gg 10$, we conservatively assume that
the beaming fraction goes as $f_{\rm b,X}\approx 10^{-3}/\dot{m}_{3}$
(Sadowski&Narayan15). This beaming effect is due to the influence of dense
accretion disk winds with large optical depth along the disk plane, which
causes the photons (Eddington-limited at each radius; Eq. 15) to diffuse
inward until they are channeled along the central disk funnel’s axis (King09).
Therefore, the effective (isotropic-equivalent) X-ray luminosity—asssuming
that the kinetic luminosity of the wind is comparable to the X-ray luminosity
(Sadowski&Narayan15)—is,

$L_{\rm X}=\frac{L_{\rm bol}}{f_{\rm b,X}}\approx 10^{45}\,{\rm
erg\,s}^{-1}(1+\ln\dot{m}_{5})\dot{m}_{5}.$ (17)

The large optical depth due to the disk winds permeating the collimated
accreting funnel would mean that the photons would emerge from a ‘photosphere’
at a radius,

$r_{\rm ph}=\frac{L_{\rm X}\sigma_{\rm T}(1-v_{\rm j}/c)}{4\pi m_{\rm
p}c^{2}v_{\rm j}}\approx 10^{12}\,{\rm
cm}\,(1+\ln\dot{m}_{5})\dot{m}_{5}\left(\frac{v_{\rm j}}{0.3\,c}\right)^{-1}$
(18)

with an effective temperature of,

$\displaystyle k_{\rm B}T_{\rm eff}$ $\displaystyle=\left(\frac{L_{\rm
X}}{4\pi r_{\rm ph}^{2}\sigma_{\rm SB}}\right)^{1/4}\approx 0.06\,{\rm keV}$
$\displaystyle\times(1+\ln\dot{m}_{5})^{-1/4}\dot{m}_{5}^{-1/4}\left(\frac{v_{\rm
j}}{0.3\,c}\right)^{1/2},$ (19)

where, $\sigma_{\rm SB}=5.67\times 10^{-5}\,{\rm g\,s^{-3}\,K^{-4}}$ is the
Stefan-Boltzmann constant. This temperature is obtained by setting the Thomson
optical depth along the fast wind/jet funnel $\tau_{\rm T,j}=\int(1-v_{\rm
j}/c)n_{\rm j}\sigma_{\rm T}dr=1$, where the radial matter number density
profile along the channel is taken to be $n_{\rm j}=L_{\rm X}/(4\pi
r^{2}m_{\rm p}c^{2}v_{\rm j})$. Interestingly, hypernebulae could manifest as
the enigmatic ultraluminous supersoft X-ray sources (Kahabka_06; Micic+22).

The number density of these thermal photons from the accretion disk
photosphere at the location of the shock is,

$\displaystyle n_{\rm X}^{\rm th}$ $\displaystyle=\frac{L_{\rm X}}{4\pi
R^{2}k_{\rm B}T_{\rm eff}c}=4\times 10^{6}\,{\rm
cm}^{-3}\,(1+\ln\dot{m}_{5})^{5/4}\dot{m}_{5}^{5/4}$
$\displaystyle\times\left(\frac{v_{\rm
j}}{0.3\,c}\right)^{-1/2}\times\begin{cases}v_{\rm
w,9}^{-2}\left(\frac{t}{80\,{\rm yr}}\right)^{-2}&({\rm t<t_{\rm free}})\\\
\left(\frac{L_{\rm w,42}}{n_{1}}\right)^{-2/5}\left(\frac{t}{80\,{\rm
yr}}\right)^{-6/5}&(t>t_{\rm free})\end{cases}$ (20)

(2) Disk/Comptonized nonthermal emission: The thermal UV/X-ray disk emission
will be up-scattered by a Comptonizing corona. Given the large Thomson depth
in the polar funnel, the Comptonized nonthermal emission would be generated in
the funnel/jet walls at large distances $>r_{\rm ph}$ (e.g., powered by
processes such as magnetic reconnection Sironi&Beloborodov20; sridhar+21c;
Sridhar+23a) along the jet boundary (Ripperda+20). The Comptonized photons
follow a powerlaw spectrum $N(\epsilon)\propto\epsilon^{-q}$ with say, a
powerlaw index $q=2$, from the thermal peak energy $\epsilon_{\rm pk,th}^{\rm
eff}=3k_{\rm B}T_{\rm eff}$ to the pair creation threshold energy $2m_{\rm
e}c^{2}$ (Svensson87). Assuming that the energy contained in the nonthermal
photons is comparable to the thermal photons, the number density of the
nonthermal photons can be estimated as,

$\displaystyle n_{\rm X}^{\rm nth}$ $\displaystyle\sim\frac{L_{\rm X}}{4\pi
R^{2}\epsilon_{\rm pk,th}^{\rm eff}c\ln{(2m_{\rm e}c^{2}/\epsilon_{\rm
pk,th}^{\rm eff})}}$ $\displaystyle\approx 4\times 10^{5}\,{\rm
cm}^{-3}\,\left(\frac{v_{\rm
j}}{0.3\,c}\right)^{-1/2}(1+\ln\dot{m}_{5})^{5/4}\dot{m}_{5}^{5/4}$
$\displaystyle\times\begin{cases}v_{\rm w,9}^{-2}\left(\frac{t}{80\,{\rm
yr}}\right)^{-2}&({\rm t<t_{\rm free}})\\\ \left(\frac{L_{\rm
w,42}}{n_{1}}\right)^{-2/5}\left(\frac{t}{80\,{\rm
yr}}\right)^{-6/5}&(t>t_{\rm free})\end{cases}.$ (21)

(3) Free-free X-ray emission: The FS becomes radiative at a time $t>t_{\rm
rad}^{\rm fs}$ when $f_{\rm rad}>1$ (see bottom panel of Fig. 1). At high
temperatures of downstream gas $T>10^{7}$ K (Eq. 9), free-free cooling
dominates over line emission. The luminosity of this free-free X-ray emission
from the FS dominates the free-free emission from both the WTS and the JTS, as
the latter stay adiabatic. The number density of the FS free-free photons at
the location of the FS is,

$n_{\rm X}^{\rm ff}=\frac{\Lambda_{\rm ff}}{\Lambda}\frac{f_{\rm rad}L_{\rm
s}}{4\pi R^{2}\bar{E}c}$ (22)

where the factor $\Lambda_{\rm ff}/\Lambda$ accounts for the fraction of the
FS’s total luminosity radiated via free-free emission.

The photon number densities from the aforementioned three emission processes
are shown in Fig. 3. Despite the fact that $n_{\rm X}^{\rm ff}$ increases with
time as $\propto t^{4/5}$ whereas $n_{\rm X}^{\rm th}$ and $n_{\rm X}^{\rm
nth}$ decrease with time as $\propto t^{-6/5}$, we find $n_{\rm X}^{\rm
th}>n_{\rm X}^{\rm nth}\gg n_{\rm X}^{\rm ff}$ at all times $t<t_{\rm active}$
for all $\dot{m}$. Therefore, for the rest of our work, we consider the
interaction of just the thermal disk emission and its nonthermal Comptonized
counterpart with the JTS-accelerated protons. The number densities of the
other sources of hard X-ray/gamma-ray photons e.g., from the jet knots (Safi-
Harb+22) are negligible compared to $n_{\rm X}^{\rm th/nth}$.

Figure 3: Background X-ray photon number density at the location of the shock.
Solid curves represent the disk blackbody photons (Eq. 2.2), dashed curves
represent the Compton upscattered nonthermal photons (Eq. 2.2), and the dotted
curves represent the free-free X-ray photons (Eq. 22). Different colors
indicate different accretion rates. The left-facing triangles in the lower
frame represent the active duration of the engine ($t_{\rm active}$; Eq. 1)
for different $\dot{m}$, and the bottom-facing triangles along the upper frame
mark the end of the free expansion of the shell ($t_{\rm free}$; Eq. 2).

##  3 Model for neutrino emission

###  3.1 Pion production

Immersed in an UV/X-ray photon field (§2.2), the protons accelerated at the
JTS can engage in the photomeson process through the $\Delta^{+}$-resonance
(Berezinskii+90),

$\displaystyle p+\gamma$ $\displaystyle\longrightarrow$
$\displaystyle\Delta^{+}\ \longrightarrow\ \left\\{\begin{array}[]{l
l}p+\pi^{0}&{\rm(2/3^{rd}~{}of~{}cases)}\\\
n+\pi^{+}&{\rm(1/3^{rd}~{}of~{}cases)}\\\ \end{array}\right.\ .$ (25)

The $\Delta^{+}$-resonance decay branching ratio is such that two-thirds of
the products follow the neutral pion ($\pi^{0}$) channel, and the remaining
one-thirds will follow the charged pion ($\pi^{+}$) channel. The pion
production through the photohadronic process (Eq. 25) with thermal photons
from the accretion disk peaks at the $\Delta^{+}$-resonance at a threshold
proton energy given by (Mannheim&Biermann89; Waxman&Bahcall97),

$E_{\rm thr}\gtrsim\frac{E_{\Delta^{+}}}{k_{\rm B}T_{\rm eff}}\frac{m_{\rm
p}c^{2}}{2}\approx 2\,{\rm
PeV}(1+\ln\dot{m}_{5})^{1/4}\dot{m}_{5}^{1/4}\left(\frac{v_{\rm
j}}{0.3\,c}\right)^{-1/2},$ (26)

where $E_{\Delta^{+}}\approx 0.3$ GeV is the energy where the cross-section
peaks. However, note that $E_{\rm thr}$ can be smaller by a factor
$(\epsilon/2m_{\rm e}c^{2})^{-2}$ for photohadronic interactions with higher-
energy nonthermal photons with energies $k_{\rm B}T_{\rm
eff}{<}\epsilon{<}2m_{\rm e}c^{2}$ (more discussion around Eq. 3.4). Since
$E_{\rm thr}\gg\bar{E}$, the photohadronic reactions with thermal photons can
proceed only with the JTS-accelerated protons from the nonthermal tail with
energies $E\gg\bar{E}$ (more in §3.3).

In addition to this channel, pions will also be produced through
$pp$-interactions that follow,

$\displaystyle p+p$ $\displaystyle\longrightarrow$ $\displaystyle\Delta^{++}\
\longrightarrow\ \left\\{\begin{array}[]{l l}p+p+\pi^{+}\\\ p+p+\pi^{0}\\\
p+p+\pi^{+}+\pi^{-}\\\ p+n+\pi^{+}+\pi^{+}+\pi^{-}\\\ \end{array}\right.\ .$
(31)

Note that the other branches of Eq. 31 involving the production of Kaons
through the decay of $\phi$-mesons are not shown. In the following section, we
calculate the particle interaction timescales that will determine the
efficiency of pion production through $p\gamma$ (Eq. 25) and $pp$ channels
(Eq. 31).

###  3.2 Particle interaction timescales

The most energetic particles will propagate inside the magenetized region of
the hypernebula on the light-crossing time, $t_{\rm cross}=R/c$, while most of
the less-energetic particles that are trapped inside the nebula by the
turbulent magnetic field would escape the system on a longer, shell expansion
timescale (Eq. 5). The particles that are the most relevant for neutrino
production are the ones that are trapped within the X-ray solid angle (Eq. 16)
and can engage in photohadronic interactions. We can correspondingly define a
timescale for laterally crossing the X-ray cone as,

$t_{\rm cross}^{\rm cone}\simeq\frac{R}{v_{\rm fs}}\theta_{\rm X}\approx
0.5\,{\rm yr}\left(\frac{t}{\rm 80\,yr}\right)\dot{m}_{5}^{-1}.$ (32)

Let’s first consider the photohadronic interactions of the protons accelerated
at the JTS with the thermal and nonthermal photons (Eqs. 2.2 and 2.2). We can
define an optical depth for this interaction, $\tau_{\rm
p\gamma}^{\rm(th,nth)}=\sigma_{\rm p\gamma}\kappa_{\rm p\gamma}Rn_{\rm
X}^{\rm(th,nth)}$ where, $\sigma_{\rm p\gamma}\sim 5\times 10^{-28}\,{\rm
cm}^{2}$ is the cross-section of the photopion production
(ParticleDataGroup:2004fcd) at the $\Delta^{+}$-resonance, and $\kappa_{\rm
p\gamma}\sim 0.15$ is the average fraction of energy lost from a proton per
elastic collision. This can be deemed the photopion production efficiency:
i.e., the fraction of the shock-accelerated protons that will be converted
into pions before streaming out of the shock, so long as they satisfy the
threshold energy requirements (Eq. 26). The corresponding $p\gamma$
interaction timescale is,

$\displaystyle t_{\rm p\gamma}^{\rm(th,nth)}=\frac{t_{\rm cross}^{\rm
cone}}{\tau_{\rm p\gamma}^{\rm(th,nth)}}\approx
8(1+\ln\dot{m}_{5})^{-5/4}\dot{m}_{5}^{-7/4}\left(\frac{v_{\rm
j}}{0.3\,c}\right)^{1/2}$
$\displaystyle\times\begin{cases}(10^{2},10^{3})\,{\rm yr}\,v_{\rm
w,9}\left(\frac{t}{80\,{\rm yr}}\right)^{2}&({\rm t<t_{\rm free}})\\\
(10^{2},10^{3})\,{\rm yr}\,\left(\frac{L_{\rm
w,42}}{n_{1}}\right)^{1/5}\left(\frac{t}{80\,{\rm yr}}\right)^{8/5}&(t>t_{\rm
free})\end{cases}.$ (33)

The other relevant channel of neutrino production is the hadronuclear ($pp$)
interaction of protons accelerated by the JTS, whose interaction timescale is,

$\displaystyle t_{\rm pp}$ $\displaystyle\simeq\frac{\theta_{\rm X}}{n_{\rm
p}\sigma_{\rm pp}\kappa_{\rm pp}v_{\rm fs}}\approx 1.5\times 10^{5}\,{\rm yr}$
$\displaystyle\times\begin{cases}\dot{m}_{5}^{-3/2}v_{\rm
w}^{2}\left(\frac{t}{80\,{\rm yr}}\right)^{2}&(t<t_{\rm free})\\\
\dot{m}_{5}^{-1/2}\left(\frac{L_{\rm
w,42}}{n_{1}}\right)^{-1/5}n_{1}^{-1}\left(\frac{t}{80\,{\rm
yr}}\right)^{2/5}&(t>t_{\rm free})\end{cases},$ (34)

where $\sigma_{\rm pp}\sim 10^{-25}\,{\rm cm}^{-2}$ (at $\sim$1 EeV) and
$\kappa_{\rm pp}\sim 0.5$ (ParticleDataGroup:2004fcd). Note that $t_{\rm pp}$
increases sharply as $\propto t^{2}$ only till $t=t_{\rm free}$, beyond which
it increases moderately as $\propto t^{2/5}$, given now the constant influx of
protons into the shell from the CSM (see Eq. 7). The total pion creation rate
due to photopion and hadronuclear interactions can be written as,

$1/t_{\pi,\rm cre}=1/t_{\rm p\gamma}^{\rm th}+1/t_{\rm p\gamma}^{\rm
nth}+1/t_{\rm pp}.$ (35)

The top panel of Fig. 1 shows the aforementioned interaction timescales for
different accretion rates ($\dot{m}=10^{3}$: dotted curves, $\dot{m}=10^{5}$:
solid curves, and $\dot{m}=10^{7}$: dashed curves). At all times, for all
choices of $\dot{m}$, $t_{\rm cross}<t_{\rm exp}$. However, $t_{\rm
cross}^{\rm cone}$ is dependent on $\dot{m}$ through the beaming fraction (Eq.
16): $t_{\rm cross}^{\rm cone}<t_{\rm cross}$ for $\dot{m}\gtrsim 10^{4}$. For
low $\dot{m}\sim 10^{3}$ systems, the timescales follow the hierarchy, $t_{\rm
exp}\ll t_{\rm p\gamma}^{\rm th}<t_{\rm p\gamma}^{\rm nth}<t_{\rm pp}$ at
earlier times $t\lesssim t_{\rm free}$, and at later times $t\gg t_{\rm
free}$, the constant shell density changes the hierarchy to $t_{\rm pp}<t_{\rm
p\gamma}^{\rm th}<t_{\rm p\gamma}^{\rm nth}$. For higher $\dot{m}\sim 10^{7}$
cases, the timescale hierarchy obeys $t_{\rm p\gamma}^{\rm th}<t_{\rm
p\gamma}^{\rm nth}<t_{\rm cross}<t_{\rm pp}<t_{\rm exp}$ during most of their
lifetimes. Overall, (although not shown in Fig. 1), $t_{\rm\pi,cre}$ overlaps
with $t_{\rm p\gamma}^{\rm th}$ at earlier times and with $t_{\rm pp}$ at
later times (if applicable, e.g., for low $\dot{m}\sim 10^{3}$ cases).

Pion production happens efficiently as long as $t_{\pi,\rm cre}$ is smaller
than both the $t_{\rm exp}$ and the synchrotron cooling timescale of the
proton given by,

$\displaystyle t_{\rm p,syn}$ $\displaystyle=\frac{6m_{\rm
p}^{3}c}{\sigma_{\rm T}m_{\rm e}^{2}\bar{\gamma}_{\rm p}B_{\rm sh}^{2}}\sim
10^{10}\,{\rm yr}\,\sigma_{\rm j,-1}^{-1}\left(\frac{v_{\rm
j}}{0.3\,c}\right)^{-1}\eta_{-1}^{-1}$
$\displaystyle\times\begin{cases}\dot{m}_{5}^{-1}\left(\frac{t}{\rm
80\,yr}\right)^{2}&(t<t_{\rm free})\\\ \dot{m}_{5}^{-3/5}v_{\rm
w,9}^{-6/5}n_{1}^{-2/5}\left(\frac{t}{\rm 80\,yr}\right)^{6/5}&(t>t_{\rm
free})\end{cases}$ (36)

where we take $\bar{\gamma}_{\rm p}\sim 10^{6}$ to be the average Lorentz
factor of the accelerated protons with energy in the relevant energy range for
neutrino production ($E_{\rm thr}<E<E_{\rm max}$; more in Eqs. 41 and 26).
Comparing the curves in the top panel of Fig. 1 demonstrates the generally
weak impact $t_{\rm p,syn}$ has on any of the relevant timescales for pion
production. Given these timescales, the overall efficiency of pion production
can be quantified by defining a ‘suppression factor’,

$f^{\rm p}_{\rm sup}={\rm min}\left[1,\frac{t_{\rm cross}}{t_{\pi,\rm
cre}},\frac{t_{\rm cross}^{\rm cone}}{t_{\pi,\rm cre}},\frac{t_{\rm
p,syn}}{t_{\pi,\rm cre}}\right].$ (37)

As the bottom panel of Fig. 1 shows, the neutrino production is not suppressed
(i.e., $f_{\rm p}^{\rm sup}\sim 1$) for high $\dot{m}\gtrsim 10^{6}$ systems
throughout their lifetimes. This is because of the large background photon
density due to a stronger beaming effect that reduces the photohadronic
interaction timescales. On the other hand, for weaker accretion rates
$\dot{m}\lesssim 10^{6}$ the pion production efficiency decreases with time as
$f_{\rm sup}^{\rm p}\propto t^{-1}$ until $t=t_{\rm free}$ owing to its
expansion and smaller background photon density. Soon after, the shell attains
a constant CSM density (leading to a modestly increasing $t_{\rm pp}$), and a
consequently decreasing $t_{\rm pp}/t_{\rm exp}$ would mean that
$t_{\rm\pi,cre}$ gets dominated by hadronic interactions. This effect is most
pronounced for the long-lived $\dot{m}=10^{3}$ case where $f_{\rm p}^{\rm
sup}$ increases again at later times. Overall, $f_{\rm p}^{\rm sup}\ll 1$ for
the majority of the lifetimes of low $\dot{m}\lesssim 10^{6}$ systems.

Furthermore, pions and muons could suffer synchrotron losses (pion/muon
‘damping’) which can inhibit neutrino production if the former proceeds at a
timescale shorter than the latter. We check its effect by calculating the
critical energy above which their synchrotron losses dominate (Metzger+20),

$\displaystyle\varepsilon_{\pi,\mu}^{\rm syn}=\left(\frac{6\pi
m_{\pi,\mu}c}{\sigma_{\rm T}B_{\rm
sh}^{2}\tau_{\pi,\mu}}\right)^{1/2}\left(\frac{m_{\pi,\mu}}{m_{\rm
e}}\right)m_{\pi,\mu}c^{2}\approx(30,2)\,{\rm ZeV}\,\sigma_{\rm j,-1}^{-1/2}$
$\displaystyle\times\left(\frac{v_{\rm
j}}{0.3\,c}\right)^{-1/2}\eta_{-1}^{-1/2}\times\begin{cases}\dot{m}_{5}^{-1/2}\left(\frac{t}{\rm
80\,yr}\right)&(t<t_{\rm free})\\\ \dot{m}_{5}^{-3/10}v_{\rm
w,9}^{-3/5}n_{1}^{-1/5}\left(\frac{t}{\rm 80\,yr}\right)^{3/5}&(t>t_{\rm
free})\end{cases},$ (38)

where $m_{\pi}$ and $m_{\mu}$ are the masses of pions and muons, and
$\tau_{\pi}=2.6\times 10^{-8}$ s and $\tau_{\mu}=2.2\times 10^{-6}$ s are
their decay times (ParticleDataGroup:2004fcd). These energies—as we show below
in Eq. 41—exceed the maximum attainable proton energies during the life of any
hypernebula. Therefore, it is fair to assume that majority of the neutrino
production from the decay of pions and muons completes before their radiative
losses take over. Furthermore, due to the small proton number density at the
shock radius (see Eq. 7), the timescales of pion and muon hadronuclear
interactions, $t_{\rm\pi,p}$ and $t_{\rm\mu,p}$ greatly exceed their lifetimes
$\tau_{\pi}$ and $\tau_{\mu}$. Therefore, the efficiency of neutrino
production is not additionally suppressed.

###  3.3 Proton energization

Figure 4: Limitations on the proton and neutrino energies. Solid curves
represent the maximum proton energy (Eq. 41); the corresponding maximum
attainable neutrino energy is denoted by dotted curves. Dashed curves denote
the minimum attainable neutrino energy through the photomeson process (Eq.
3.4). Different colors represent different accretion rates. The left-facing
triangles in the lower frame represent the active duration of the engine
($t_{\rm active}$; Eq. 1) for different $\dot{m}$, and the bottom-facing
triangles along the upper frame mark the end of the free expansion of the
shell ($t_{\rm free}$; Eq. 2).

A fraction of the JTS’s energy goes into heating protons (with an average
energy $\bar{E}\simeq 17$ MeV; Eq. 9). The fraction of the total shock kinetic
luminosity transferred to the protons is taken to be,

$f_{\rm E}=\begin{cases}1/2&(E\sim\bar{E})\\\ \epsilon_{\rm
rel}&(E\gg\bar{E}),\end{cases}$ (39)

where $\epsilon_{\rm rel}\lesssim 0.01$ is the efficiency of the acceleration
process (e.g., for mildly-relativistic shocks; Fang+20).444Note that
nonthermal particle acceleration can be inefficient in highly magnetized
shocks (Sironi&Spitkovsky11). The time evolution of the proton luminosity,
$L_{\rm p}(E)=f_{\rm E}(E)L_{\rm s}$ for different $\dot{m}$ is shown in Fig.
5. We assume that protons with $E\gg\bar{E}$ are accelerated by diffusive
shock acceleration mechanism (DSA; Blandford&Ostriker78). The maximum energy
that the protons can attain through this process is determined by the
competition between the diffusive acceleration timescale of particles with
Larmor radius $r_{\rm L}$ (e.g., Caprioli&Spitkovsky14),

$t_{\rm acc}\simeq\frac{2}{3}\left(\frac{r_{\rm
L}}{c}\right)\left(\frac{v_{\rm j}-v_{\rm
fs}}{c}\right)^{-2}=\frac{2Ec}{3eB_{\rm sh}(v_{\rm j}-v_{\rm fs})^{2}},$ (40)

and the timescales corresponding to various interaction and cooling processes
(e.g., in Fig. 1). This yields a maximum ion energy of,

$E_{\rm max}=\frac{3eB_{\rm sh}(v_{\rm j}-v_{\rm fs})^{2}}{2c}{\rm
min}\left[t_{\rm cross},t_{\rm cross}^{\rm cone},t_{\rm exp},t_{\rm
p\gamma}^{\rm(th,nth)},t_{\rm pp}\right].$ (41)

Eq. 41 provides a system-size-limited maximum energy of the particles
considering that the energized particles with $r_{\rm L}>R$ (and the X-ray
cone lateral radius) would escape the accelerator (Hillas_84) and do not
contribute to the neutrino production. Solid curves in Fig. 4 show the
evolution of $E_{\rm max}(t)$ for different $\dot{m}$—complicated by the
influence of the X-ray beaming cone’s opening angle on various particle
interaction timescales. The highest energy attained in $\dot{m}=10^{3}$ and
$10^{7}$ cases is $\sim 30$ PeV, and the highest energy attained overall is,
$\sim 300$ PeV, in the $\dot{m}=10^{5}$ case.

The distribution of protons that are accelerated by the DSA into a nonthermal
distribution—with energies ${\gg}\bar{E}$ and an exponential cut-off at
$E\gtrsim E_{\rm max}$—can be described by,

$\frac{dN_{\rm p}}{dE}\sim\left(\frac{E}{\bar{E}}\right)^{-q}e^{-E/E_{\rm
max}}.$ (42)

with $q\simeq 2$ being the powerlaw index that is typically expected from DSA
(our fiducial case; Blandford&Eichler87); steeper spectra with $q\gtrsim 2$
could also arise as the magnetic fluctuations behind the FS drift away at the
local Alfvén speed (‘postcursor’; Diesing+21). Fig. 9 shows example ($q=2$)
proton spectra for different $\dot{m}$ at times $t_{\rm free}$ and $t_{\rm
active}$.

###  3.4 Pion decay to neutrinos

The pions produced through the $p\gamma$ (Eq. 25) and $pp$ (Eq. 31) channels
will eventually decay into neutrinos and gamma-rays following,

$\displaystyle\begin{aligned} \pi^{0}\longrightarrow&\gamma\gamma\\\
\pi^{+}\longrightarrow&\mu^{+}+\nu_{\mu}\\\
&\rotatebox[origin={c}]{180.0}{$\Lsh$}e^{+}+\nu_{\rm e}+\bar{\nu}_{\mu}\\\
\pi^{-}\longrightarrow&\mu^{-}+\bar{\nu}_{\mu}\\\
&\rotatebox[origin={c}]{180.0}{$\Lsh$}e^{-}+\bar{\nu}_{\rm
e}+\nu_{\mu}\end{aligned}$ (43a)

In the absence of significant pion/muon damping (Eq. 3.2), the ratio of the
neutrino flavors resulting from pion decay, post-oscillation, is $\nu_{\rm
e}:\nu_{\mu}:\nu_{\tau}=1:1:1$ by the time they reach Earth (Becker08).

Let’s now consider the photomesonic process with thermal background photons,
where the produced pions decay into neutrinos with $\sim 5\%$ of the proton’s
energy $\varepsilon\approx E/20$ (Stecker68; Kelner&Aharonian08). This sets
the lowest attainable neutrino energy to be,

$\displaystyle\varepsilon_{\rm min}$ $\displaystyle\approx{\rm
max}[\bar{E}/20,E_{\rm thr}/20]$ $\displaystyle\approx 100\,{\rm
TeV}(1+\ln\dot{m}_{5})^{1/4}\dot{m}_{5}^{1/4}\left(\frac{v_{\rm
j}}{0.3\,c}\right)^{-1/2}.$ (44)

Note that $\varepsilon_{\rm min}$ can be smaller by a factor
$(\epsilon/2m_{\rm e}c^{2})^{-2}$ for photohadronic interactions with
nonthermal photons with energies $k_{\rm B}T_{\rm eff}{<}\epsilon{<}2m_{\rm
e}c^{2}$. This would effectively place the minimum neutrino energy cut-off at
$\varepsilon_{\rm min}[\epsilon_{\rm pk,th}^{\rm eff}/k_{\rm B}T_{\rm
eff}]^{-2}\sim 10$ TeV (see Eq. 26 and also §2.2(2)). Production of even lower
energy neutrinos mediated by interactions with nonthermal photons of energies
$\epsilon>\epsilon_{\rm pk,th}^{\rm eff}$ is highly suppressed. On the other
hand, the absolute maximum neutrino energy as set by the accelerator is
$\varepsilon_{\rm max}\sim E_{\rm max}/20$. Fig. 4 shows $\varepsilon_{\rm
min}$ as the horizontal dashed line, which is independent of the age of the
hypernebula, and cares only about $k_{\rm B}T_{\rm eff}$ through $\dot{m}$
(Eq. 2.2). $\varepsilon_{\rm max}$ is denoted by dotted curves, which trace
$E_{\rm max}$ with a scaling offset of 1/20. In general, hypernebulae generate
neutrinos throughout their lifetime. An exception to this are the high
$\dot{m}>10^{7}$ systems where their short lifetimes (and therefore, a smaller
acceleration region) restrict $E_{\rm max}<\varepsilon_{\rm min}$ and prevent
any neutrino production.

The neutrino luminosities from the JTS-accelerated $p\gamma$ and $pp$
interactions are,

$\displaystyle\begin{aligned} L_{\rm\nu}^{\rm p\gamma}&=&\frac{3}{8}f_{\rm
p\gamma}L_{\rm p}&=\frac{27\pi}{64}f_{\rm p\gamma}f_{\rm E}R^{2}\rho_{\rm
jts}(v_{\rm j}-v_{\rm fs})^{3}\\\ L_{\rm\nu}^{\rm pp}&=&\frac{1}{2}f_{\rm
pp}L_{\rm p}&=\frac{9\pi}{16}f_{\rm pp}f_{\rm E}R^{2}\rho_{\rm jts}(v_{\rm
j}-v_{\rm fs})^{3}\end{aligned},$ (45a)

where the neutrino conversion efficiency factors are,

$\displaystyle\begin{aligned} f_{\rm p\gamma}^{\rm(th,nth)}&=1-e^{-\tau_{\rm
p\gamma}}=1-e^{-t_{\rm cross}^{\rm cone}/t_{\rm p\gamma}^{\rm(th,nth)}}\\\
f_{\rm pp}&=1-e^{-\tau_{\rm pp}}=1-e^{-t_{\rm cross}^{\rm cone}/t_{\rm
pp}}\end{aligned}.$ (46a)

The factor 3/8 in Eq. 45(a) accounts for the fraction of the proton energy
that goes into neutrino products (Waxman&Bahcall97; Fang&Metzger17), i.e.,
charged pions are produced in a $p\gamma$ interaction only half the time on
average, and when a pion decays, $\sim$3/4 of its energy goes to neutrinos
(see Eqs. 25 and 43). The factor 1/2 in Eq. 45(b) accounts for the fact that
pions are produced with $\sim$2/3 probability in a $pp$ interaction (see Eq.
31), and $\sim$3/4 of their energy is carried away by neutrinos. Finally, each
hypernebula provides a total neutrino energy flux per flavor that follows the
proton energy spectrum at a given time as,

$E^{2}_{\nu}\frac{dN_{\nu}}{dE_{\nu}}=\frac{E^{2}_{\rm p}}{2}\frac{dN_{\rm
p}}{dE}\left[\frac{3}{4}f_{\rm p\gamma}^{\rm(th,nth)}+f_{\rm pp}\right].$ (47)

##  4 Observational Implications

###  4.1 Luminosity and volumetric rates

Figure 5: Solid, dashed, and dotted curves show the temporal luminosity
evolution of jet termination shock (Eq. 8), and the protons and neutrinos
energized there (Eqs. 39 and 47), respectively. Different colors represent
different accretion rates. The left-facing triangles in the lower frame
represent the active duration of the engine ($t_{\rm active}$; Eq. 1) for
different $\dot{m}$, and the bottom-facing triangles along the upper frame
mark the end of the free expansion of the shell ($t_{\rm free}$; Eq. 2)

Fig. 5 shows the combined $L_{\nu}^{\rm p\gamma}+L_{\nu}^{\rm pp}$ neutrino
‘light curves’ for different $\dot{m}$. The neutrino luminosity closely
follows the profile of the JTS and proton luminosity for $\dot{m}=10^{7}$
case, but it evolves as $\propto t^{-1}$ for $\dot{m}\lesssim 10^{6}$ cases.
For even lower $\dot{m}\lesssim 10^{4}$ cases, the neutrino luminosity
increases again before turning off: this is due to the onset of hadronic
($pp$) channels of neutrino production at later times. The peak neutrino
luminosity is larger for higher $\dot{m}$: for the fiducial cases shown in
Fig. 5, the brightest emission ($L_{\nu}=10^{41}\,{\rm erg\,s^{-1}}$) is seen
from $\dot{m}=10^{7}$ during the early moments of JTS onset. Whether neutrino
emission (and its temporal evolution) can be observed from an individual
hypernebulae, however, depends on their volumetric rate.

The volumetric rate of hypernebulae can be assumed to follow common envelope
events involving a black hole/neutron star with donor stars more massive than
$m_{\star}=10$. Schroder+20 find this frequency to be $\sim 0.6\%$ core-
collapse supernova rate, corresponding to a local-universe (redshift $z=0$)
volumetric rate of ${\cal R}_{0}\sim 10^{3}\,{\rm Gpc^{-3}\,yr^{-1}}$
(Strolger+15; Vigna-Gomez+18). This rate is also in agreement (within an order
of magnitude) with that of engine-powered fast blue optical transients (FBOTs;
${\cal R}_{0}\sim 300\,{\rm Gpc^{-3}\,yr^{-1}}$) which could be failed common
envelope events (Coppejans+20; Ho+21b), the occurrence rate of massive star
binaries with mass-transfer rates $\dot{M}\gtrsim\dot{M}_{\rm Edd}$ (${\cal
R}_{0}\sim 100\,{\rm Gpc^{-3}\,yr^{-1}}$; Pavlovskii+17), and $\sim 1-10\%$ of
all the luminous red novae merger events (Karambelkar+22). Assuming an active
timescale of $\sim$10 yr (for the neutrino-brightest $\dot{m}=10^{7}$ case),
one active hypernebula can be expected within $50$ Mpc. In the most optimistic
scenario (considering $\eta\sim 1$), the peak neutrino flux from a ‘local’
($\sim$50 Mpc) source would be $10^{-10}\,{\rm
GeV\,s^{-1}\,cm^{-2}\,sr^{-1}}$, suggesting that the prospects of detecting
neutrinos from individual hypernebulae are not encouraging.

Recently, IceCube reported a flux of $E_{\nu}^{2}\Phi\nu({\rm 1\,TeV})\simeq
5\times 10^{-8}\,{\rm GeV^{-1}\,cm^{-2}\,s^{-1}}$ from NGC 1068 (at 14.4 Mpc;
IceCubeCollaboration_22). The peak flux from a single hypernebula situated at
that distance ($\sim$10% chance) would be $E_{\nu}^{2}\Phi\nu({\rm
1\,TeV})\simeq 3\times 10^{-9}\,{\rm GeV^{-1}\,cm^{-2}\,s^{-1}}$. The observed
flux of NGC 1068 could in principle, be explained by a hypernebula if either
(1) the shocks are weakly magnetized $\sigma_{\rm j}\lesssim 10^{-3}$, with a
larger $\epsilon_{\rm rel}\sim 0.1$ (Sironi&Spitkovsky11), or (2) in a
scenario where a cluster of hypernebulae migrate to the AGN core (Secunda+20).
Although, we note that other natural channels involving AGN coronae could very
well be operating (Murase+20; Inoue+22; Eichmann+22).

As far as the HE diffuse background neutrino flux is concerned, its overall
normalization can be explained by non-relativistic shock-powered transients as
long as they obey (Fang+20),

$\left(\frac{{\cal R}_{0}}{10^{5}\,{\rm
Gpc^{-3}\,yr^{-1}}}\right)\left(\frac{E_{\rm tot}}{5\times 10^{50}\,{\rm
erg}}\right)\left(\frac{\epsilon_{\rm rel}}{1\%}\right)\sim 1,$ (48)

where the total output energy is $E_{\rm tot}$, of which, a fraction
$\epsilon_{\rm rel}$ is directed towards accelerating non-thermal ions (see
Eq. 39). While the volumetric rates of hypernebulae are smaller than the
suggested value in Eq. 48 (which is normalized to that of core-collapse
supernovae; Briel+22), their longer lifetimes, and a steady injection of
energy into their surroundings brings the total energy budget of any
hypernebula (nearly independent of $\dot{m}$) to

$E_{\rm tot}\approx\int^{t_{\rm active}}L_{\rm w}dt\sim 5\times 10^{52}\,{\rm
erg}\left(\frac{M_{\star}}{30\,M_{\odot}}\right)v_{\rm w,9}^{2}.$ (49)

This makes hypernebulae excellent candidates for explaining a significant
portion of the HE diffuse background neutrino flux.

As shown in Fig. 5, neutrinos are generated from a hypernebula throughout its
lifetime, but with neutrino spectral properties constantly evolving as the
hypernebula expands (see also Fig. 9). We, therefore, calculate a
representative neutrino flux for each hypernebula model by averaging Eq. 47
over its lifetime. The integrated observed neutrino flux contributed by all
hypernebulae from different redshifts, taking into account the energy loss due
to cosmological expansion can be written as,

$E_{\nu}^{2}\Phi_{\nu}(E_{\nu})=\frac{{\cal
R}_{0}}{4\pi}\int_{0}^{6}dz\frac{c}{(1+z)^{2}H(z)}f(z)\left({E_{\nu}^{\prime}}^{2}\frac{dN_{\nu}}{dE_{\nu}^{\prime}}\right)(z),$
(50)

where $H(z)=H_{0}\left[\Omega_{\rm m}(1+z)^{3}+\Omega_{\Lambda}\right]^{1/2}$
is the Hubble constant at redshift $z$, where we take $H_{0}\simeq 70\,{\rm
km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{\rm m}=0.315$, and $\Omega_{\Lambda}=0.7$
(for a $\Lambda$CDM flat cosmology; Planck_Collaboration20).
$E_{\nu}^{\prime}=E_{\nu}(1+z)$ is the redshifted neutrino energy, $f(z)={\cal
R}(z)/{\cal R}_{0}$ captures the possibly different volumetric rates of the
source at different redshifts. We assume ${\cal R}(z=0)\sim 10^{3}\,{\rm
Gpc^{-3}\,yr^{-1}}$, and that ${\cal R}$ follows star formation history. This
lets us obtain $f(z)$ using the following empirical fit (reliable upto $z\sim
6$) of Strolger+04,

$f(z)=\frac{a\left[t^{b}e^{-t/c}+de^{d(t-t_{0})/c}\right]}{de^{-dt_{0}/c}},$
(51)

where $a=0.182$, $b=1.26$, $c=1.865$, $d=0.071$, $t=1/H_{0}(1+z)$ is the look-
back time, and $t_{0}=1/H_{0}\sim 14\,$Gyr is the age of the universe. Since
the volumetric rate of hypernebulae scales as the star formation history
($f(z)$, which increases with redshift) the total number of hypernebulae
contributing to the neutrino flux is $\sim 10\times$ what one would get if a
constant ${\cal R}={\cal R}_{0}$ is assumed.

###  4.2 Model neutrino flux

Figure 6: All-flavor HE diffuse background neutrino spectra. The IceCube
measurements shown are from high-energy cascades (turquoise markers;
Aartsen+20), and the high-energy starting event sample (grey markers;
IceCube_HESE+21). For comparison, the Fermi extra-galactic isotropic gamma-ray
background is shown with green markers (Ackermann+15). Error bars and upper
limits represent 68% confidence intervals. The grey bands of different shades
denote the sensitivity limits of current (IceCube : light grey; Aartsen+13)
and future (Trinity : medium grey; Nepomuk_Otte_19, IceCube-Gen2 : dark grey;
IceCube-Gen2+14; Aartsen+21) neutrino detectors. Left panel: Solid curves are
the total time-averaged, volume-integrated neutrino flux obtained for models
of hypernebula with different accretion rates (color-coded) but otherwise the
same fiducial parameters ($v_{\rm w}=0.03c,v_{\rm j}=0.3c,n=10,\sigma_{\rm
j}=0.1,\eta=0.5$; see §2 for more details). Dashed and dotted curves represent
the contribution of photomesonic ($p\gamma$) interaction with thermal and non-
thermal photons, respectively, and dash-dotted curves denote the contribution
of hadronic interactions ($pp$). The maximum neutrino energies for different
$\dot{m}$ are color-coded and marked as down-facing triangles along the top
frame. Right panel: The brown curve shows the flux of a model
($\dot{m}=10^{7}$, other parameters listed above the panel) representative of
the hypernebula population that provides a reasonable description of the
observed flux and slopes of the HE diffuse background neutrino flux.

Left panel of Fig. 6 shows the time-averaged, volume-integrated (Eq. 50)
neutrino flux expected from hypernebulae with different $\dot{m}$ (denoted by
different colors). The processes contributing to the overall neutrino flux are
shown with different linestyles i.e., the dominant $p\gamma$ photomesonic
reactions with thermal and nonthermal background photons, and the highly sub-
dominant hadronic $pp$ reactions, whose contribution to the neutrino flux is
relevant only at lower-energies $<$100 GeV. Also shown for comparison are the
observed HE diffuse background neutrino flux (Aartsen+20), the Fermi/LAT IGRB
(Ackermann+15) (circular markers with error bars), and the sensitivity limit
of future high-energy and UHE neutrino observatories. The peak of the neutrino
spectra in our models is primarily determined by the $p\gamma$ interactions
with the higher density thermal photons. The high-energy cut-off is determined
by the system parameters, independent of the neutrino production channel (Eq.
41 and Fig. 4). On the other hand, the low-energy turnovers are set by the
three different neutrino production channels: the turnovers at $\sim$10 TeV
and $\sim$1 TeV are due to the $p\gamma$ interactions with the thermal and
nonthermal photons, respectively (see discussion following Eq. 3.4). The third
and the lowest-energy cut-off at $\sim$135 MeV (below the shown energy range
of Fig. 6) is set by the photohadronic component of the spectrum, at the pion
creation threshold energy, corresponding to its rest-mass energy. Note that
the location of these low-energy turnovers are only mildly dependent on the
hypernebula’s parameters, and are predictive of neutrino emission from
hypernebulae (also demonstrated in Appendix C).

Left panel of Fig. 6 demonstrates that hypernebulae with $\dot{m}\sim 10^{7}$
(brown curve; with otherwise the same fiducial parameters mentioned in §2)
alone can comfortably explain all of the observed 10 TeV to PeV neutrino flux
(the chosen parameter combination underestimates the 10 TeV neutrinos, which
shall be addressed below with different model parameters). Also, the
$\dot{m}\sim 10^{6}$ and $\dot{m}\sim 10^{5}$ models, alone, can supply just
$\sim 1\%$ and $\sim 0.01\%$ of the observed flux, and thus the $\dot{m}\ll
10^{6}$ cases can be ruled out for being insignificant contributors. Although
not visibly shown, the flux of the $\dot{m}\sim 10^{3}$ case can be as high as
$\gtrsim 10^{-12}\,{\rm GeV\,s^{-2}\,cm^{-2}\,sr^{-1}}$—values comparable to
$\dot{m}\sim 10^{4}$ case—because their long lifetimes allow $pp$ channel to
dominate the waning $p\gamma$ channel for neutrino production.

In Appendix C, we explore a range of physically reasonable parameter values
(about our fiducial choice) for hypernebulae, and their resultant model
neutrino flux. Based on this exploration, we arrive at a visual ‘best-fit’
model to explain the observed neutrino flux, with the following hypernebula
parameters: $\dot{m}=10^{7},v_{\rm w}/c=0.07,v_{\rm
j}/c=0.17,\eta=1,n=1,\sigma_{\rm j}=10^{-2}$, and $q=2$. This model is shown
in the right panel of Fig. 6. The robust low-energy turnover feature of our
model provides a natural explanation for the flat/mildly rising nature of the
observed neutrino spectrum from 10–50 TeV, followed by a broken powerlaw at
$\gtrsim$100 TeV. Sridhar&Metzger22 rule out certain range of hypernebula
parameter space based on their expected radio synchrotron luminosity and by
comparing it with the VLASS data (e.g., $v_{\rm j}/c\ngtr 0.5$); the best-fit
values obtained above are consistent with this constraint. We would like to
emphasize that this best-fit is not intended to suggest that the entire
observed flux is produced by one special kind of hypernebula, instead, is
intended to offer a representative parameter space around which hypernebulae
should occupy, in order to explain all of the observed neutrino flux. While
the shown best fit model has a high-energy cut-off at a few PeV, the parameter
exploration in Appendix C reveals models that are likely to be sources of UHE
neutrinos detectable by upcoming neutrino observatories such as IceCube-Gen2
(IceCube-Gen2+14; Aartsen+21) and Trinity (Nepomuk_Otte_19).

The models provided above assume that the entire population of hypernebulae is
powered by a single $\dot{m}$ source, which certainly need not be the case.
This simplistic assumption is due to the underlying uncertainty in the
relationship between ${\cal R}_{0}$ and $\dot{m}$. We explore this further in
Fig. 7, where the different contours in the ${\cal R}_{0}-\dot{m}$ space
indicate the model neutrino flux at different energies (as a fraction of the
observed neutrino flux). Sources with $\dot{m}<4\times 10^{5}$ cannot explain
even 1% of the observed flux for any reasonable value of ${\cal R}_{0}\lesssim
10^{3}\,{\rm Gpc^{-3}\,yr^{-1}}$ (red contours). The hypernebule models with
$\dot{m}\gtrsim 10^{6}$ can explain a significant fraction ($\gtrsim 10\%$;
yellow contours) of the observed 10 TeV–1 PeV flux even with a lower local
volumetric rate of ${\cal R}_{0}={\rm few}\times 100\,{\rm
Gpc^{-3}\,yr^{-1}}$. Only $\dot{m}\sim 10^{7}$ models are able to explain all
of the observed flux (green contours). This is relevant because the life of
each hypernebulae might culminate with a highly unstable mass transfer stage
(with $\dot{m}\gtrsim 10^{7}$), before plunging into the common envelope
phase. On the other hand, in a conservative scenario where hypernebulae
accrete at a constant rate throughout their lives and their volumetric rates
are uniformly distributed between different $10^{3}\lesssim\dot{m}\lesssim
10^{7}$, then the total contribution of hypernebulae to the observed 10 TeV–1
PeV neutrino flux would be $\sim 25\%$, and will be dominated by the
$\dot{m}\sim 10^{7}$ cases. In all our analysis, we do not consider sources
accreting at rates $\dot{m}\gg 10^{7}$ because, their short lifetimes would
preclude any proton from being accelerated to energies $\gtrsim 0.1E_{\rm
thr}$ (Eq. 26) to engage in neutrino production through photomesonic channels
(the factor of 0.1 accounts for $p\gamma$ interaction with the higher energy
nonthermal photons, which requires a smaller $E_{\rm thr}$; vertical grey
region in Fig. 7). In the following section, we consider the gamma-ray
properties of the more relevant $\dot{m}=10^{7}$ scenario.

Figure 7: Neutrino flux for a range of volumetric rates ${\cal R}_{0}$
(y-axis) from accretion-powered-hypernebulae with different $\dot{m}$ (bottom
x-axis), with other parameters following our ‘best fit’ (see right panel of
Fig. 6); the top x-axis denotes the active timescale (Eq. 1) for a given
$\dot{m}$. The solid, dashed and dotted curves represent the model fluxes in
10 TeV, 100 TeV, and 1 PeV, normalized to the observed fluxes of
$10^{-7},~{}4\times 10^{-8},~{}{\rm and}~{}2\times 10^{-8}\,{\rm
GeV\,s^{-1}\,cm^{-2}\,sr^{-1}}$, respectively. Red, yellow, and green contours
denote the parameter space that can supply 1%, 10%, and 100% of the observed
flux. The top grey region denotes the parameter space excluded based on the
volumetric rate constraint of ${\cal R}_{0}\lesssim 10^{3}\,{\rm
Gpc^{-3}\,yr^{-3}}$ set by the rate of common envelope events involving
compact objects (Schroder+20). The right grey region is excluded because
extreme $\dot{m}$ events fail to produce neutrinos as their maximum proton
energies ($E_{\rm max}$; Eq. 41) is less than the minimum required proton
energy to produce neutrinos upon interaction with even the higher-energy
nonthermal photons (Eq. 26).

###  4.3 $\gamma$-ray attenuation

An inevitable byproduct of neutrino production are gamma-rays (from neutral
pion decay; see Eq. 43a), with flux and spectral properties similar to that of
neutrinos. The Fermi IGRB observations indicate a lower gamma-ray flux than
neutrinos (Ackermann+15) suggesting that the primary sources of HE diffuse
background neutrino should be faint gamma-ray emitters. In this section, we
discuss the opacity of hypernebula to the gamma-rays generated along with the
neutrinos, and the impact of the intervening medium as the gamma-rays
propagate to Earth. The emitted gamma-rays can be attenuated and reprocessed
at least by the following three intrinsic processes, and Fig. 8 shows some of
the relevant optical depths at times $t=t_{\rm free}$ and $t=t_{\rm active}$.

(1) Compton down-scattering: The gamma-ray photons can lose energy and be
attenuated by scattering off the electrons in the ejecta. The effective
optical depth for the photons to lose most of their initial energy is
(Vurm&Metzger21),

$\displaystyle\tau_{\rm C}$ $\displaystyle\simeq\frac{\dot{x}_{\rm
C}}{c}\frac{R}{c}(1+\tau_{\rm KN})\approx\tau_{\rm
T}\frac{\ln(1+x)}{(1+3x)}\left[1+\frac{3\tau_{\rm
T}}{8x}\ln\left(1+\frac{8x}{3}\right)\right]$ $\displaystyle\approx
10^{-9}\,\ln\left(\frac{x}{10^{6}}\right)\left(\frac{x}{10^{6}}\right)^{-1}$
$\displaystyle\times\begin{cases}\dot{m}_{5}v_{\rm
w,9}^{-2}\left(\frac{t}{80\,{\rm yr}}\right)^{-2}&(t<t_{\rm free})\\\
\dot{m}_{5}^{1/5}v_{\rm w,9}^{2/5}n_{\rm 1}^{4/5}\left(\frac{t}{80\,{\rm
yr}}\right)^{3/5}&(t>t_{\rm free})\end{cases},$ (52)

where $\tau_{\rm KN}$ is the Klein-Nishina optical depth for gamma-ray photons
with dimensionless energies $x\equiv E_{\gamma}/m_{\rm e}c^{2}\gg 1$
(Klein&Nishina29); the factor $1+\tau_{\rm KN}$ accounts for the diffusion
time of the photons through the ejecta. $\tau_{\rm T}\sim n_{\rm e}\sigma_{\rm
T}R$ is the Thomson optical depth with a scattering cross section $\sigma_{\rm
T}=6.6\times 10^{-25}\,{\rm cm}^{-2}$ for an electron number density $n_{\rm
e}\sim n_{\rm p}$, and the particle energy loss rate is $\dot{x}_{\rm
C}/x\approx c\sigma_{\rm T}n_{\rm e}\ln(1+x)/(1+3x)$. Clearly, even the low-
energy ($E_{\gamma}\sim 1$ GeV) gamma-ray photons do not suffer any losses due
to Compton scattering, and the higher energy photons only experience even
lesser attenuation due to this process.

(2) Bethe-Heitler process: gamma-rays can be attenuated due to pair production
upon interaction with relativistic ions in the ejecta (Bethe&Heitler34). The
corresponding optical depth is given by (Chodorowski+92),

$\displaystyle\tau_{\rm BH}=n_{\rm p}\sigma_{\rm BH}R$ $\displaystyle\approx
10^{-5}\ln\left(\frac{x}{10^{6}}\right)$
$\displaystyle\times\begin{cases}\dot{m}_{5}v_{\rm
w,9}^{-2}\left(\frac{t}{80\,{\rm yr}}\right)^{-2}&(t<t_{\rm free})\\\
\dot{m}_{5}^{1/5}v_{\rm w,9}^{2/5}n_{\rm 1}^{4/5}\left(\frac{t}{80\,{\rm
yr}}\right)^{3/5}&(t>t_{\rm free})\end{cases}.$ (53)

The scattering cross section can be approximated as Zdziarski&Svensson89,

$\sigma_{\rm BH}\simeq\frac{3}{8\pi}\alpha\sigma_{\rm
T}Z^{2}\left[\frac{28}{9}\ln(2x)-\frac{218}{27}\right]$ (54)

where $\alpha=1/137$ is the fine structure constant, $Z=1$ is the atomic
charge of the nuclei with atomic weight $A$ 555Note that if the binary is from
a massive star that lost a significant fraction of its mass prior to
overflowing its Roche Lobe, then the material being swept up at early times
might possess a different-than-solar composition. Nonetheless, this will not
qualitatively change the implication of Bethe-Heitler process to gamma-ray
escape in our case.. We can see from Fig. 8 that $\tau_{\rm BH}$ barely
depends upon $E_{\gamma}$, and the effect of Bethe-Heitler process is
negligible in attenuating the emitted gamma-rays.

(3) $\gamma\gamma\rightarrow e^{+}e^{-}$ pair production: Gamma-rays pair
produce upon interacting with ambient low-energy photons (Breit&Wheeler34).
Assuming for simplicity, a monochromatic isotropic thermal ambient UV/X-ray
photon field (with an energy given by Eq. 2.2), the energy of the gamma-ray
photon $E_{\gamma}$ at which the pair production cross-section
$\sigma_{\gamma\gamma}$ peaks is,

$E_{\gamma,\rm pk}=\frac{4m_{\rm e}^{2}c^{4}}{k_{\rm B}T_{\rm eff}}\approx
60\,{\rm GeV}(1+\ln\dot{m}_{7})^{1/4}\dot{m}_{7}^{1/4}\left(\frac{v_{\rm
j}}{0.3\,c}\right)^{-1/2}.$ (55)

The scattering cross-section can be taken as a smoothly connected piece-wise
function such that (Svensson87),

$\sigma_{\gamma\gamma}\sim\frac{3\sigma_{\rm
T}}{16}\times\begin{cases}0&(E_{\gamma}<E_{\gamma,\rm pk}/4)\\\
1&(E_{\gamma}=E_{\gamma,\rm pk})\\\ \ln(x)/x&(E_{\gamma}>E_{\gamma,\rm
pk})\end{cases}$ (56)

With this, the optical depth to $\gamma\gamma$ pair production with thermal
photons can be estimated as,

$\displaystyle\tau_{\rm\gamma\gamma}^{\rm th}$ $\displaystyle=n_{\rm X}^{\rm
th}\sigma_{\rm\gamma\gamma}R\approx
13\ln\left(\frac{x}{10^{6}}\right)\left(\frac{x}{10^{6}}\right)^{-1}(1+\ln\dot{m}_{7})^{5/4}\dot{m}_{7}^{5/4}$
$\displaystyle\times\left(\frac{v_{\rm
j}}{0.3\,c}\right)^{-1/2}\times\begin{cases}v_{\rm
w,9}^{-1}\left(\frac{t}{80\,{\rm yr}}\right)^{-1}&(t<t_{\rm free})\\\
\dot{m}_{7}^{-1/5}v_{\rm w,9}^{-2/5}n_{\rm 1}^{1/5}\left(\frac{t}{80\,{\rm
yr}}\right)^{3/5}&(t>t_{\rm free})\end{cases}.$ (57)

Likewise, one can also define an optical depth for $\gamma\gamma$ pair
production with nonthermal photons (see §2.2(2)),

$\tau_{\gamma\gamma}^{\rm nth}=n_{\rm X}^{\rm
th}\left(\frac{\epsilon}{\epsilon_{\rm pk,th}^{\rm
eff}}\right)^{-2}\sigma_{\gamma\gamma}[\epsilon]R,$ (58)

where $\sigma_{\gamma\gamma}[\epsilon]$ is constructed by taking
$E_{\rm\gamma,pk}=4m_{\rm e}^{2}c^{4}/\epsilon$ for individual photons with
energy $\epsilon$ in the nonthermal powerlaw distribution. Fig. 8 shows that
gamma-ray photons with energies in the range 1 GeV $\sim$1 PeV will be
strongly attenuated within the hypernebula by $\gamma\gamma$ pair production:
specifically, with interactions dominated by nonthermal photons in the energy
range 1 GeV-60 GeV, and by thermal photons in the energy range 60 GeV-1 PeV.
Intrinsic attenuation of $>$ PeV gamma-rays is in principle possible by their
interaction with the radio synchrotron photons. However, their small number
density $n_{\rm GHz}\lesssim 10^{-4}\,{\rm cm}^{-3}$ (Sridhar&Metzger22) would
imply a $\gamma\gamma$ pair production optical depth of
$\tau_{\gamma\gamma}^{\rm GHz}\ll 1$, and the higher energy gamma-rays
($E_{\gamma}\gtrsim$1 PeV) would still escape freely.

The escaping $>$PeV gamma-ray photons will interact with the Extragalactic
Background Light (EBL) and the Cosmic Microwave Background (CMB) radiation on
their way to Earth. The EBL is composed of the UV-Optical-IR photons
($0.1-1000\,\mu{\rm m}$) emitted by stars and galaxies since the epoch of
reionization, which can pair produce upon interaction with gamma-rays with
$E_{\gamma}\gtrsim$10 TeV (Nikishov1961; Gould&Schreder67a; Gould&Schreder67b,
see also Eq. 55). In Fig. 8, we show the optical depth of EBL and CMB photons
(Stecker+06) to gamma-rays, and show that all the escaping gamma-rays from a
hypernebula—which happen to have energies $E_{\gamma}>$PeV—will be attenuated.
At different locations, the non-linear dynamics of the pair production
involving processes such as inverse-Compton scattering ($e\gamma\rightarrow
e\gamma$) and triple pair production process ($e\gamma\rightarrow
ee^{+}e^{-}$) can initiate electromagnetic cascades until the energy of the
later generation of particles falls below the threshold to sustain these
reactions. As a result, the injected TeV-PeV gamma-rays could eventually
appear in the MeV-GeV band, with a spectral shape independent of the injected
gamma-ray injection spectrum and the target soft photon spectrum
(Berezinskii&Smirnov75), instead primarily dependent on the minimum neutrino
energy, $\varepsilon_{\rm min}$ (Fang+22). For our fiducial choice of
parameters, as $\tau_{\rm C}\ll 1$, (see Fig. 8) implying no gamma-ray
cascades in hypernebulae.

Figure 8: Optical depths of various interactions to gamma-ray photons with
different energies. Solid and dotted curves are at times $t=t_{\rm active}$
(Eq. 1) and $t=t_{\rm free}$ (Eq. 2), respectively, for our fiducial model,
but with the more relevant $\dot{m}=10^{7}$ case. Orange and yellow curves
represent $\gamma\gamma\rightarrow e^{-}e^{+}$ pair production with the
accretion disk thermal (Eq. 4.3) and nonthermal Comptonized photons (Eq. 58),
respectively; Bethe-Heitler interaction with the ions in the ejecta is
represented by red curves (Eq. 4.3); green curves represent the Compton down-
scattering in the Klein-Nishina regime with electrons in the ejecta (Eq. 4.3),
and the interaction with extragalactic background light (EBL) and cosmic
microwave background (CMB) radiation, for source at a redshift $z=0.2$, is
represented by maroon curves: the energies where EBL and CMB interactions
dominate are noted near the top frame of the figure.

##  5 Conclusions

We have developed a model for high-energy neutrino emission from
‘hypernebulae’. A hypernebula is a new class of spatially compact
($\lesssim$pc), decades-millenia long transient ‘bubble’ inflated by hyper-
Eddington accredion disk winds that are more powerful than the typical Ultra
Luminous X-ray sources (Sridhar&Metzger22). The rapid runaway mass transfer
($\dot{m}\sim 10^{3}-10^{7}$) onto a stellar mass compact object (black hole
or a neutron star) required to inflate hypernebulae are naturally expected
from an evolved companion star (e.g., imminent giant) leading up to a common-
envelope event. For the purpose of our ‘proof-of-concept’ calculations, we
consider a fiducial binary containing a massive star
($M_{\star}=30\,M_{\odot}$) accreting onto a black hole of fixed mass
($M_{\bullet}=10\,M_{\odot}$), whose active duration is varied by varying the
accretion rate ($10^{3}<\dot{m}<10^{7}$). The intrinsic properties of the
hypernebula, and the observables of the neutrino emission (light curves,
energy spectra, volume-integrated observed flux) are continuously evaluated
across the free-expansion, deceleration, and post-active phases, as a function
of the slow, wide-angled disk wind ($\dot{M}\sim\dot{M}_{\rm w},v_{\rm w}$),
and the fast, collimated disk wind/jet ($\eta=\dot{M}_{\rm j}/\dot{M},v_{\rm
j},\sigma_{\rm j}$) properties. The take-away results of this paper are as
follows:

(1) The primary source of neutrino production in this model is the adiabatic
fast wind/jet termination shock, where the shock-accelerated protons produce
pions (that eventually decay into neutrinos) upon their photohadronic
($p\gamma$) interactions with the thermal and nonthermal X-ray photons from
the accretion disk. The hadronic ($pp$) interactions play a sub-dominant role
except during the late stages of a hypernebula’s lifetime, particularly for
low $\dot{m}$ systems.

(2) The highest energy neutrinos are emitted during the end of a hypernebula’s
lifetime (when they are the largest in size). The maximum energy is found to
be $\sim$300 PeV for $\dot{m}=10^{5}$ case. UHE neutrinos are highly sought-
after over the past half-a-century, and are likely to be discovered within the
next couple of decades by upcoming neutrino observatories (BrancoValera+22)
such as CHANT (Neronov+17), Trinity (Nepomuk_Otte_19), IceCube-Gen2 detectors
(IceCube-Gen2+14; Aartsen+21), Askaryan Radio Array (Ara_Collaboration+11),
and the ARIANNA detectors (Barwick+15). Hypernebulae are promising candidates
for diffuse extragalactic ($\lesssim$EeV) UHE neutrinos.

(3) A distinct feature of hypernebulae as neutrino sources is their property
of robust minimum neutrino energy, that is set solely by the steady background
photon energies: the thermal photons set the prominent $\sim$10 TeV low-energy
spectral turnover, the nonthermal photons introduce a turnover at $\sim$1 TeV,
and the $pp$ interactions set the lowest-energy turnover at 135 MeV. IceCube-
DeepCore (Abbasi+12), KM3NeT (Adrian-Martinez+16; Fermani20, being built in
the Mediterranean sea;) and the Gigaton Volume Detector (in lake Baikal;
Allakhverdyan+21) may extend and improve the current neutrino sensitivities to
sub-TeV regime. These facilities will be able to test the hypernebula model by
better characterizing the TeV and sub-TeV spectrum of the HE diffuse
background neutrinos.

(4) Assuming that the volumetric rates of hypernebulae follow common envelope
events involving compact objects, an engine accreting at $\dot{m}=10^{7}$
could comfortably explain $\sim$100% of the HE diffuse background neutrino
flux. In a conservative scenario, if the volumetric rate of hypernebulae were
to be uniformly distributed among different accretion rate engines
$10^{3}<\dot{m}<10^{7}$, hypernebulae could explain (depending on the mass-
transfer rate) $\gtrsim$25% of the observed neutrino flux.

(5) The gamma-rays that are inevitably produced along with neutrinos are
attenuated within the hypernebulae primarily by the Breit-Wheeler
$\gamma\gamma\rightarrow e^{+}e^{-}$ annihilation process. The gamma-ray
photons in the energy range $1<E_{\gamma}/{\rm GeV}<60$ are attenuated by
nonthermal background photons and the higher energy gamma-rays
$60<E_{\gamma}/{\rm GeV}<10^{6}$ are attenuated by thermal background photons.
The escaping $>$PeV gamma-rays will be attenuated by their interaction with
the extragalactic background light. Therefore, hypernebulae are gamma-ray
faint sources of neutrino emission, thus satisfying the isotropic gamma-ray
background constraints of Fermi/LAT.

Finally, we would like to emphasize that hypernebulae are radio-synchrotron-
bright objects that can be detected as point sources with all-sky radio
surveys such as FIRST (Becker+95) or VLASS (Lacy+20). A detailed study of the
electromagnetic counterparts to hypernebulae—including light curves,
spectra—can be found in Sridhar&Metzger22. A population of FRBs could be
powered by flares along hyper-Eddington accretion disk jets (Sridhar+21b;
Bhandari+22). Recently, Abbasi+22b presented the search results for spatially-
and temporally-coincident neutrino emission from FRBs with IceCube Cascade
data, and found no significant clustering around the short-lived and beamed
FRB event. This result is in agreement with theoretical predictions
(Metzger+20). However, in a scenario of accretion-powered FRBs, the
hypernebula surrounding the central FRB engine could serve as a persistent
radio, and neutrino source, thus also contributing to the HE diffuse
background neutrino flux.

##  6 Acknowledgement

NS acknowledges the hospitality of the Department of Astronomy at Cornell
University, and the Center for Cosmology and AstroParticle Physics at The Ohio
State University for their generous hospitality, where parts of this work were
conceived. NS would like to thank John Beacom and Todd Thompson for fruitful
discussions, and thank Stephanie Wissel for sharing the details on various
Neutrino experiments’ sensitivities. N.S. acknowledges the support from NASA
(grant number 80NSSC22K0332), NASA FINESST (grant number 80NSSC22K1597), and
Columbia University Dean’s fellowship. K.F. is supported by the Office of the
Vice Chancellor for Research and Graduate Education at the University of
Wisconsin-Madison with funding from the Wisconsin Alumni Research Foundation.
K.F. acknowledges support from NASA (80NSSC22K1584) and National Science
Foundation (PHY-2110821).

##  Appendix A Timescales

A list of all the relevant timescales of the model introduced in this paper,
and their definitions, is provided in Table 1.

Table 1: Model timescales. Variable | C.f. | Definition
---|---|---
$t_{\rm cross}$ | Fig. 1 | Light-crossing timescale of the hypernebula.
$t_{\rm active}$ | Eq. 1 | Active duration of the dynamically unstable accretion phase at a mass transfer rate of $\dot{M}$.
$t_{\rm free}$ | Eq. 2 | Free-expansion timescale beyond which the wind ejecta shell starts decelerating.
$t_{\rm exp}$ | Eq. 5 | Expansion timescale of the shocked gas.
$t_{\rm cool}$ | Eq. 10 | Cooling timescale of the shock-heated gas.
$t_{\rm rad}^{\rm fs}$ | Eq. 12 | The time after which the forward shock becomes radiative.
$t_{\rm cross}^{\rm cone}$ | Eq. 32 | The time particles spend within the solid angle of background X-ray photon field.
$t_{\rm p\gamma}^{\rm(th,nth)}$ | Eq. 3.2 | Cooling time due to photomeson interaction between protons and background thermal and nonthermal photon fields.
$t_{\rm pp}$ | Eq. 3.2 | Hadronuclear interaction timescale.
$t_{\pi,\rm cre}$ | Eq. 35 | Pion creation timescale.
$t_{\rm p,syn}$ | Eq. 3.2 | Proton synchrotron cooling timescale.
$\tau_{\pi,\mu}$ | Eq. 3.2 | Decay times of pions and muons.
$t_{\rm acc}$ | Eq. 40 | Diffusive shock acceleration timescale.

##  Appendix B Proton energy spectra

The temporal evolution of the JTS-energized proton energy spectra is shown in
Fig. 9, for different $\dot{m}$. The thermal population of protons occupies a
distinct hump in the energy spectrum at $E\sim 3\bar{E}\sim 10^{-2}$ GeV (see
Eq. 9). At $t\sim t_{\rm free}$ (dotted curves), both the spectral high-energy
cutoff, $E_{\rm max}(t)$, and the spectral normalization increase with
$\dot{m}$. However, at later times $t\sim t_{\rm active}$ (solid curves), the
protons accelerated by the higher-$\dot{m}$ engines are less energetic than
the lower-$\dot{m}$ cases given the shorter lifetimes of the former, and the
resulting size-limited energization constraints (Hillas_84). This reduces the
high-energy spectral cut-off for high-$\dot{m}$ cases while retaining a
constant $E_{\rm thr}$. This limits neutrino production only to cases with
$\dot{m}\lesssim{\rm few}\times 10^{7}$ (see also Fig. 7).

Figure 9: Energy spectra of protons energized at the JTS. Solid curves are at
$t=t_{\rm active}$ (Eq. 1), and the dotted curves are at an earlier time,
$t=t_{\rm free}$ (Eq. 2). Different colors indicate different accretion rates.
The threshold energy to produce neutrinos via photohadronic interaction with
thermal background photons ($E_{\rm thr}$; Eq. 26) is indicated along the top
frame of the figure.

##  Appendix C Hypernebula parameter exploration

Figure 10: Neutrino flux models for a range of hypernebula’s physical
parameters. The fiducial model is represented by solid brown curve in all the
panels (but with $\dot{m}=10^{7}$; see also §2). The dashed and dotted curves
in each panel show how altering one parameter (that is mentioned above each
panel) changes the observed flux about the fiducial model. See Fig. 6 for a
description of the other labels and markers in each panel.

Fig. 10 shows the influence of each relevant physical parameter of hypernebula
on the model neutrino flux, for $\dot{m}=10^{7}$. As mentioned in §4.2, Fig.
10 demonstrates the robustness of the low-energy turnovers for a wide range of
hypernebula parameter space. The neutrino flux due to hadronuclear
interactions (${\rm few}\times 10^{-10}\,{\rm GeV\,s^{-1}\,cm^{-2}\,sr^{-1}}$
at $\lesssim 10^{2}$ GeV) is not strongly affected by any parameter except for
$\eta$, which proportionately scales the normalization of the model flux at
all energies as it changes only the JTS power (panel [f]). When $v_{\rm w}\ll
v_{\rm j}$, an increase in $v_{\rm w}$ increases the shell radius and
therefore the shock luminosity, but a further increase ($v_{\rm w}\lesssim
v_{\rm j}$) would mean a slower JTS (i.e., $v_{\rm j}-v_{\rm fs}$), and a
lower $n_{\rm p}$ (Eq. 7), which will sharply decrease the (photohadronic
contribution to the) neutrino flux, and $E_{\rm max}$ (panel [a]). An increase
in $v_{\rm j}$, on the other hand, monotonically increases $E_{\rm max}$, and
the overall flux normalization—primarily by boosting the neutrino production
through thermal photohadronic interactions (panel [d]). Both the upstream
density $n$ and $\sigma_{\rm j}$ act identically on the spectra (through
$B_{\rm sh}$; Eq. 2.1)—solely affecting the high-energy cut-off (panels [b]
and [e]). An increase in the proton spectral index $q$ steepens the neutrino
spectrum by decreasing the thermal photohadronic component (hump at $\sim$1
PeV) more than the nonthermal photohadronic component (hump at $\sim$1 TeV) of
the neutrino spectrum (panel [c]). The parameter exploration in Fig. 10
suggests that hypernebulae—with $v_{\rm w}/c\sim 0.3,v_{\rm j}/c\sim
0.5,n\gtrsim 100\,{\rm cm^{-3}}$, and $\sigma_{\rm j}\gtrsim 0.1$—could be
promising candidates for $>$100 PeV UHE neutrinos that could be detected by
future facilities such as IceCube-Gen2 (IceCube-Gen2+14; Aartsen+21) and
Trinity (Nepomuk_Otte_19).