Signature of Collapsars as Sources for High-energy Neutrinos and $r$-process Nuclei

Collapsars produced by the collapse of massive stars into black holes have long been considered a leading candidate for powering Type Ibc supernovae and long gamma-ray bursts (GRBs) of both the bright and low-luminosity (LLGRBs) varieties Woosley (1993); MacFadyen and Woosley (1999); Woosley and Bloom (2006); Campana et al. (2006); Hjorth and Bloom (2012); Kumar and Zhang (2014a). Shocks associated with the propagation of the collapsar jet through the stellar envelope and/or the black-hole accretion disk wind were also proposed as candidate sites for producing $\sim$ TeV–PeV high energy (HE) neutrinos Mészáros and Waxman (2001); Razzaque et al. (2003, 2004); Ando and Beacom (2005); Razzaque et al. (2005); Horiuchi and Ando (2008); Enberg et al. (2009); Bartos et al. (2012); Murase and Ioka (2013); Fraija (2014); Xiao and Dai (2014); Bhattacharya et al. (2015); Varela et al. (2015); Tamborra and Ando (2016); Senno et al. (2016, 2018); Denton and Tamborra (2018a, b); He et al. (2018); Chang et al. (2022); Guarini et al. (2023). These neutrinos may contribute significantly to the diffuse flux detected by IceCube, whose astrophysical origin remains unknown despite recent reports of potential association of several events with a blazar Aartsen et al. (2018a, b), three tidal disruption events Stein et al. (2021); Reusch et al. (2022); van Velzen et al. (2021), and an active galaxy Abbasi et al. (2022a).

The ground-breaking multi-messenger observations of GW170817 linked binary neutron star mergers (BNSMs) to both short GRBs and the production of heavy elements beyond iron via the rapid neutron-capture process ($r$-process) Abbott et al. (2017a, b); Kasen et al. (2017); Drout et al. (2017); Cowperthwaite et al. (2017); Villar et al. (2017); Shibata et al. (2017); Metzger (2020); Watson et al. (2019). The astrophysical environments of BNSMs and collapsars are rather similar in that both host accretion disks and the associated jets, thereby producing GRBs and potentially HE neutrinos. Interestingly, Ref. Siegel et al. (2019) showed that the physical conditions of disk outflows in collapsars resemble those of outflows in BNSMs, and proposed that collapsars are a likely or even the dominant site for producing $r$-process nuclei. This proposal remains hotly debated and requires further theoretical and observational efforts to clarify Siegel (2019); Macias and Ramirez-Ruiz (2019); Bartos and Marka (2019); Miller et al. (2020); Fujibayashi et al. (2020); Brauer et al. (2021); Fraser and Schönrich (2022); Barnes and Metzger (2022); Just et al. (2022); Fujibayashi et al. (2022); Anand et al. (2023).

In this work, we investigate for the first time a novel connection between HE neutrinos and $r$-process nucleosynthesis in collapsars. Unstable neutron-rich heavy nuclei are produced by rapid neutron capture during the $r$-process. Their $\beta$ decay produces $\bar{\nu}_{e}$ with energy $E_{L}\sim 4$ MeV on a timescale of $\sim 1$ s everywhere inside the outflow that expands with a typical velocity $v_{\rm ej}\sim(0.05$–$0.3)c$ (see Fig. 1). Following flavor oscillations, these low energy (LE) antineutrinos can annihilate HE neutrinos of the corresponding flavor produced by shocks at radius $R$. Efficient annihilation via the $Z$-resonance requires $s=2E_{H}E_{L}(1-\cos\theta)\sim m_{Z}^{2}$, where $\theta\sim\mathcal{O}(1)$ is the intersection angle and $m_{Z}$ is the $Z$ mass. Therefore, the HE neutrino flux at $E_{H}\gtrsim\mathcal{O}(100$–1000) TeV is expected to be affected. As we will show, such annihilation can leave clear imprints on the energy spectrum and flavor composition of the emergent HE neutrinos.

The $\beta$-decay of the newly-synthesized $r$-process nuclei not only produces a profuse $\bar{\nu}_{e}$ flux, but also provides a dominant source of energy to heat the associated ejecta Metzger et al. (2010). We estimate the $\bar{\nu}_{e}$ emission using the power generated by $\beta$-decay. In our benchmark study, we simply assume that the ejecta expands with a constant velocity $v_{\rm ej}$, forming a steady spherical “wind” with a constant mass outflow rate $\dot{M}$ (see Appendix A). As the ejecta expands, the power generated by $\beta$-decay per unit mass stays almost constant for a period of $\sim T_{r}$ when the $r$-process produces neutron-rich nuclei far from stability. After neutron capture ceases, that power approximately follows a power-law decline. We take $\dot{\epsilon}_{\bar{\nu}_{e},0}\eta(t)$ as the part of the power per unit mass that is carried away by $\bar{\nu}_{e}$, where $\dot{\epsilon}_{\bar{\nu}_{e},0}$ sets the magnitude and $\eta(t)$ characterizes the time evolution with $\eta(0)\approx 1$. From numerical calculations, $\eta(t)$ can be fitted as Rosswog et al. (2014)

Ignoring nuclear shell effects and Coulomb correction, we simply take the $\bar{\nu}_{e}$ spectrum to be $\propto E^{2}(Q-E)^{2}$ with $Q$ being the $Q$-value. The energy-differential emission rate of $\bar{\nu}_{e}$ per unit mass (i.e., emissivity) at radius $r=v_{\rm ej}t$ can be estimated as

where the term in the brackets is the approximate normalized $\bar{\nu}_{e}$ spectrum with an average energy of $\langle E_{L}\rangle=Q/2=4$ MeV. Note that both $\dot{\epsilon}_{\bar{\nu}_{e},0}$ and $T_{r}$ depend on the electron fraction $Y_{e}$ of the ejecta Rosswog et al. (2014). Assuming a typical value of $Y_{e}=0.2$ for collapsar outflows and guided by the power used in Ref. Wu et al. (2019), we take $\dot{\epsilon}_{\bar{\nu}_{e},0}=5\times 10^{18}$ erg/g/s and $T_{r}=0.4$ s.

The $\bar{\nu}_{e}$ emission is approximately isotropic. Referring to Fig. 1, we estimate the $\bar{\nu}_{e}$ intensity (in units of cm^-2 s^-1 MeV^-1 sr^-1) at radius $R$ and angle $\theta$ by integrating the emissivity along the path length $l$:

where $\rho(r)=\dot{M}/(4\pi r^{2}v_{\rm ej})$ is the mass density of the ejecta at radius $r$, and $r$ is a function of $R$, $\theta$, and $\theta^{\prime}$. Unless noted otherwise, we take $\dot{M}=0.02\,M_{\odot}/{\rm s}$ Hayakawa and Maeda (2018); Fujibayashi et al. (2022) and $v_{\rm ej}=0.05c$ MacFadyen and Woosley (1999); Miller et al. (2020). The corresponding total number density of LE $\bar{\nu}_{e}$ at $R\sim 10^{10}$ cm is $\sim 10^{23}~{}{\rm cm}^{-3}$.

Following production, the LE $\bar{\nu}_{e}$ undergo flavor oscillations prior to interacting with the HE neutrinos. The oscillation scenario depends on the matter densities at the $\bar{\nu}_{e}$ emission site and the interaction site. We consider three sets of probabilities $P_{\bar{\nu}_{\alpha}}$ for the initial $\bar{\nu}_{e}$ to be a $\bar{\nu}_{\alpha}$ at the interaction site: $(P_{\bar{\nu}_{e}},P_{\bar{\nu}_{\mu}},P_{\bar{\nu}_{\tau}})\approx(0.55,0.18,0.27)$ (pure vacuum oscillations), $(0.675,0.095,0.23)$ (NH), and $(0.022,0.545,0.433)$ (IH). The latter two sets assume adiabatic flavor evolution following emission at high densities for normal (NH) and inverted (IH) neutrino mass hierarchy, respectively (see Appendix B).

Relativistic jets in collapsars induce shocks that can accelerate electrons and protons. The accelerated electrons can produce $\gamma$-rays or X-rays through synchrotron radiation or inverse Compton scattering, thereby making bright GRBs, while the HE protons can collide with photons or stellar matter to make $\pi^{\pm}$ and $K^{\pm}$, thereby producing HE neutrinos through meson decays. Shocks emerging at different sites during jet propagation lead to different scenarios for HE neutrino production Waxman and Bahcall (1997, 2000); Li et al. (2002); Guetta et al. (2004); Murase and Nagataki (2006); Murase et al. (2006); Murase (2007); Hummer et al. (2012); Zhang and Kumar (2013); Liu and Wang (2013); Bustamante et al. (2015); Nir et al. (2016); Biehl et al. (2018); Guo et al. (2020); Kimura (2022). To demonstrate the effects of LE antineutrinos from $\beta$-decay, we focus on HE neutrinos produced during jet propagation deep inside the progenitors of collapsars. Specifically, we consider proton acceleration at internal shocks before jet collimation Mészáros and Waxman (2001); Murase and Ioka (2013).

Ref. Murase and Ioka (2013) emphasized that shocks inside stars are likely mediated by radiation and in contrast to collisionless shocks Levinson and Bromberg (2008); Katz et al. (2010), may not support efficient particle acceleration and HE neutrino production. It was shown later, however, that if the jets are mildly magnetized, a strong collisionless subshock may occur within the radiation-mediated shock (RMS) Beloborodov (2017); Lundman and Beloborodov (2019) to allow particle acceleration at, e.g., mildly relativistic internal shocks Crumley et al. (2019). Indeed, using relativistic magnetohydrodynamic simulations of magnetized jets from BNSMs, Ref. Gottlieb and Globus (2021) found that subshocks can emerge at small radii and produce HE neutrinos at internal shocks, collimation shocks, and shock breakout.

Interestingly, a mild magnetization is also essential to avoid heavy baryon loading of the jets from mixing with the cocoon material Gottlieb et al. (2020, 2021); Gottlieb and Globus (2021). With negligible jet-cocoon mixing, the jet Lorentz factor could grow almost linearly with radius due to adiabatic expansion (see, e.g., Piran (2004); Meszaros (2006); Kumar and Zhang (2014b)) and reach $\gtrsim 100$ deep inside the star Mizuta and Aloy (2009); Gottlieb et al. (2021). Internal shocks caused by collisions of slow and rapid jets can occur at $R_{\rm is}\approx\Gamma_{s}^{2}c\delta t=3\times 10^{9}\Gamma^{2}_{s,1}\delta t_{-3}$ cm, where $\Gamma_{s}$ is the typical Lorentz factor of the slow jets and $\delta t$ is the variability timescale. Here and below, $A_{x}=A/10^{x}$, where $A$, if dimensional, is in cgs units. We assume that the typical Lorentz factor $\Gamma_{r}$ of the unshocked rapid jets is twice the Lorentz factor $\Gamma$ of the shocked jet. The relative Lorentz factor between the two is $\Gamma_{\rm rel}\approx\Gamma_{r}/(2\Gamma)+\Gamma/(2\Gamma_{r})=1.25$, consistent with the mildly relativistic shock.

In our scenario, a strong collisionless subshock forms within the internal RMS to facilitate particle acceleration. This scenario requires $\sigma_{u}\equiv B_{u}^{2}/(4\pi\rho_{u}c^{2})\gtrsim 0.01$ Beloborodov (2017); Lundman and Beloborodov (2019), where $B_{u}$ and $\rho_{u}$ are the proper magnetic field and mass density, respectively, of the unshocked upstream region. Due to shock compression, the magnetic field is amplified to $B_{d}=\xi B_{u}$ in the shocked downstream region with the compression ratio $\xi\sim 5$ for $\Gamma_{\rm rel}\approx 1.25$ Beloborodov (2017). Because HE neutrino flux would be strongly suppressed by synchrotron cooling of $\pi^{\pm}$ and $K^{\pm}$ in the enhanced magnetic field $B_{d}$, we focus on HE neutrino production via $pp$ and $p\gamma$ reactions in the upstream region (i.e., the unshocked jet), which was barely discussed in previous literature. For radiation-dominated jets at launching, the initial jet temperature is $\sim 1.3L^{1/4}_{\rm iso,52}$ MeV at $r=10^{7}$ cm Kumar and Zhang (2014a). The jets expand adiabatically and the temperature decreases as $r^{-1}$. At $r=10^{9}$–$10^{10}$ cm, the jet temperature remains a few keV, corresponding to a thermal energy fraction $\epsilon_{e,u}\gtrsim 0.1$. Proper values of magnetic field, proton number density, and photon temperature in the shocked and unshocked jets are given in Table 1.

Protons can be accelerated to high energies by crossing the shock fronts repeatedly. During this process, protons have significant probabilities to escape from both the upstream and downstream regions, resulting in a power-law spectrum Bell (1978). Recent simulations of mildly relativistic shocks showed that the accelerated protons in the upstream region are much more abundant than those in the downstream region (see Figs. 7 and 8 of Crumley et al. (2019)). To obtain a conservative estimate of the maximal proton energy $E_{p,{\rm max}}$, we compare the timescales of shock acceleration and different cooling processes for protons in the shocked jet (see Appendix C). For all the parameter sets adopted in this work, the $E_{p,{\rm max}}$ computed in the rest frame of the shocked jet is always sufficient to make PeV neutrinos ($E_{\nu}\sim 0.05\Gamma E_{p}$ from $pp$ and $p\gamma$ reactions).

We take the (unnormalized) spectrum of the accelerated protons to be $\phi_{p}(E_{p})=E_{p}^{-2}\exp(-E_{p}/E_{p,\rm max})$, and use PYTHIA 8.3 Bierlich et al. (2022) and SOPHIA Mucke et al. (2000) for $pp$ and $p\gamma$ reactions, respectively, to obtain the $\pi^{\pm}$ and $K^{\pm}$ yields at different centre-of-mass energies. Because the shocked jet is only mildly relativistic relative to the rapid jet ($\Gamma_{\rm rel}=1.25$), we ignore the Lorentz transformation for $E_{p}$ and $\phi_{p}(E_{p})$ in our calculations. We include cooling of $\pi^{\pm}$ and $K^{\pm}$ due to synchrotron radiation, inverse Compton scattering, $e^{\pm}$ pair production on photons ($\pi^{\pm}/K^{\pm}+\gamma\to\pi^{\pm}/K^{\pm}+e^{+}+e^{-}$), hadronic scattering, and adiabatic expansion (see Appendix C).

In our benchmark study, we take $R_{\rm is}=3\times 10^{9}$ cm, $L_{\rm iso}=10^{53}$ erg/s, $\Gamma_{r}=2\Gamma\approx 90L^{0.18}_{\rm iso,49}\approx 472$ Yi et al. (2017), $\sigma_{d}=\xi\sigma_{u}\approx 2\epsilon_{B,d}=0.25$, $\xi=5$, $\epsilon_{e,d}=0.5$, and $\epsilon_{e,u}=0.3$ (see Table 1). The corresponding $E_{p,{\rm max}}$ is $\approx 10^{6}$ GeV, and we take $E_{p,{\rm min}}=100$ GeV without affecting the results. We find that $e^{\pm}$ production on photons and synchrotron radiation dominate the cooling of $\pi^{\pm}$ and $K^{\pm}$ (see Appendix C). The benchmark HE neutrino spectra are shown in Fig. 2. The competition between $pp$ and $p\gamma$ reactions depends on the effective threshold for meson production and the number density of targets. The $pp$ reaction dominates at low proton energies and produces almost equal neutrino and antineutrino fluxes at $E\lesssim 30$ TeV (left panel) due to equal production of $\pi^{+}\,(K^{+})$ and $\pi^{-}\,(K^{-})$. Once the proton energy exceeds the effective threshold for the $p\gamma$ reaction, this reaction dominates due to the much larger number density of thermal photons and gives rise to higher fluxes of neutrinos than antineutrinos (left panel) due to dominant production of $\pi^{+}\,(K^{+})$ over $\pi^{-}\,(K^{-})$ and stronger synchrotron cooling of $\mu^{\pm}$. The latter effect also suppresses production of $\nu_{e}$ and $\bar{\nu}_{e}$ from $\mu^{\pm}$ decay compared to that of $\nu_{\mu}$ and $\bar{\nu}_{\mu}$ from $\pi^{\pm}\,(K^{\pm})$ decay (middle panel).

To cover a wide range of GRB luminosities, we compare the all-flavor neutrino spectra in the right panel of Fig. 2 using $L_{\rm iso}=10^{50}$–$10^{53}$ erg/s with other parameters unchanged. The variations with increasing $L_{\rm iso}$ are caused by the boost of neutrino energy due to a larger Lorentz factor and by more efficient cooling of charged mesons at relatively low energies due to a larger photon density or $B_{u}$. For $L_{\rm iso}=10^{53}$  erg/s, cooling by pair production and synchrotron radiation suppresses neutrino production from $\pi^{\pm}$ decay so much that $K^{\pm}$ decay produces half of the neutrinos at 100 TeV.

Simulations show that the nonrelativistic radioactive outflows occur within a few seconds of jet launching and last for a time comparable to the jet duration Just et al. (2022); Fujibayashi et al. (2022). Therefore, most of the HE neutrinos produced at $R_{\rm is}\sim 10^{9}$–$10^{10}$ cm can meet the LE antineutrinos from the radioactive outflows. The resulting annihilation takes $\lesssim$ 1 s, which is significantly shorter than the outflow duration. We use a one-dimensional Monte Carlo code to simulate the propagation of HE neutrinos in the LE antineutrino background and the accompanying processes, which include $\nu_{\alpha}\bar{\nu}_{\alpha}$ annihilation and decays of the produced $\pi^{\pm}$ and $K^{\pm}$. Note that we also include the annihilation of the regenerated HE $\nu_{\alpha}$ from the earlier annihilation processes. The radial step $\Delta R$ to the next interaction point for an HE $\nu_{\alpha}$ with energy $E_{H}$ at radius $R$ is determined by $\tau_{\alpha}(E_{H},R,\Delta R)=-\ln\lambda$, where $\tau_{\alpha}$ is the optical depth due to $\nu_{\alpha}\bar{\nu}_{\alpha}$ annihilation and $\lambda$ is a random number between 0 and 1. The relevant cross section $\sigma_{\nu_{\alpha}\bar{\nu}_{\alpha}}(s)$ peaks at $s=m_{Z}^{2}$, and we use PYTHIA 8.3 Bierlich et al. (2022) to track the final products for all $s$ values.

Assuming vacuum oscillations of LE antineutrinos prior to annihilation, we show the spectra of HE neutrinos following annihilation in Fig. 2. The $Z$-resonance starts to occur at $E\sim 300$ TeV. Beyond this energy, the HE neutrino flux drops but the HE antineutrino flux slightly increases due to production from $Z$ decay (left panel). HE neutrinos at these energies are mostly $\nu_{\mu}$, which are annihilated by the LE $\bar{\nu}_{\mu}$ (middle panel). Strikingly, the $\nu_{e}$ and $\bar{\nu}_{e}$ regenerated from such annihilation are even more abundant than those produced without annihilation. Consequently, the flux of $\nu_{e}+\bar{\nu}_{e}$ is very similar to that of the regenerated $\nu_{\tau}+\bar{\nu}_{\tau}$ at $E\gtrsim 300$ TeV following annihilation (middle panel). The above effects produce a unique flavor composition of $(F_{\nu_{e}+\bar{\nu}_{e}}:F_{\nu_{\mu}+\bar{\nu}_{\mu}}:F_{\nu_{\tau}+\bar{\nu}_{\tau}})_{\star}\approx(1:10:1)$ for $E\gtrsim 300$ TeV at source, which corresponds to a normalized composition of $\approx(0.21:0.42:0.37)$ at Earth following vacuum oscillations. The corresponding all-flavor spectrum with the suppression due to annihilation is shown in the right panel of Fig. 2 for different values of $L_{\rm iso}$. With an optical depth of 3–8 for $\nu_{\mu}\bar{\nu}_{\mu}$ annihilation for $E_{H}\sim 0.4$–1 PeV, the suppression factor can be as high as 3–3.5 (see the inset in the right panel of Fig. 2).

Although prompt HE neutrinos associated with bright GRBs are tightly constrained by IceCube Abbasi et al. (2012); Aartsen et al. (2015, 2016, 2017); Abbasi et al. (2022b), subleading contributions of precursor neutrinos from jet propagation inside stellar matter still remain possible. The most recent analyses showed that precursor neutrinos preceding the prompt $\gamma$-rays by tens of seconds are limited to $\lesssim 10$% of the diffuse flux at $E=100$ TeV (see Fig. 7 of Abbasi et al. (2022b)). On the other hand, LLGRBs, whose rate is $\sim 10$ times higher than that of classical GRBs Soderberg et al. (2006); Liang et al. (2007); Nakar (2015), could be caused by a shock breakout driven by a relativistic jet propagating into an extended envelope Campana et al. (2006). They may have a common collapsar origin and host jets with similar luminosities inside stars to those for bright GRBs Nakar (2015); Senno et al. (2016). Note that, in the shock breakout scenario, the observed $\gamma$-ray luminosities of LLGRBs ($L_{\gamma,\rm llgrb}^{\rm obs}\sim 10^{46}$–$10^{48}$ erg/s) are much lower than those of bright GRBs because the emission for LLGRBs spans a longer duration ($\Delta T_{\gamma,\rm llgrb}\sim 10^{3}$–$10^{4}$ s) as determined by the large shock breakout radius and is over a wider opening angle ($\theta_{\rm llgrb}\sim 1$ rad). Considering the duration correction, the beam factor correction, and a baryon loading factor $f_{p}\sim 10$, the true luminosity of LLGRB jets initially launched inside stars could be $L_{\rm iso}\sim L_{\gamma,\rm llgrb}^{\rm obs}f_{p}\frac{\Delta T_{\gamma,\rm llgrb}}{\Delta T_{\rm eng}}\frac{\theta^{2}_{\rm llgrb}}{\theta^{2}_{\rm j}}\sim 10^{51}$–$10^{53}$ erg/s, where $\Delta T_{\rm eng}\sim 10$–100 s is the central engine duration and $\theta_{\rm j}\sim 0.1$ rad is the initial jet opening angle at launch. As HE neutrinos are produced deep inside stars in our work, LLGRBs and bright GRBs are similar sources for HE neutrinos, except that they have different occurrence rates. Further, choked collapsar jets could be even more common than bright GRBs and LLGRBs if the central engine is not active for too long Sobacchi et al. (2017). To incorporate contributions from all similar sources, we introduce the parameter $f_{\star}$ as the ratio of the total rate of collapsars with mildly-magnetized jets to that of bright GRBs.

The unnormalized HE neutrino spectrum from a single source is shown in Fig. 2. To obtain the diffuse neutrino flux from all similar sources over cosmic history, we multiply the unnormalized spectrum by the following factor Murase and Ioka (2013)

where $f_{z}\sim 3$ accounts for evolution of source population and neutrino energy with redshift $z$ Waxman and Bahcall (1999), $\dot{n}_{\rm GRB}\sim 1$ Gpc^-3 yr^-1 is the rate of bright GRBs per unit volume, $E_{\rm iso}$ is the total isotropic energy of jets, and $f_{\rm cr}$ is the fraction of jet energy converted into HE protons. Taking $f_{\star}=10$, $f_{\rm cr}E_{\rm iso,53}\sim 1$, and assuming that contributions from sources at $z\sim 1$ dominate, we compare the all-flavor $E^{2}F_{\nu}(E)$ from our model to the IceCube data Aartsen et al. (2020) in Fig. 3. It can be seen that LLGRBs and choked jets with collapsar origin may explain the observed TeV–PeV flux at IceCube. More interestingly, the spectrum with suppression at $E\gtrsim 100$–200 TeV due to annihilation (and redshift) leads to better agreement with data than that without annihilation. Results using $\dot{M}=0.01$ and $0.05~{}M_{\odot}$/s are also shown for comparison. Note that spherical winds and vacuum oscillations of LE antineutrinos are assumed for our benchmark study. The annihilation effect could be similar or even enhanced for asymmetric winds while it will be reduced for the NH oscillation scenario with a smaller $P_{\bar{\nu}_{\mu}}$ (see Appendices A and B). We find that annihilation is significant under a range of assumptions (see Fig. 3). Our model, however, cannot account for the PeV events detected at IceCube despite the production of a small bump above $\sim 1$ PeV associated with a minimum in the timescale for cooling of $K^{\pm}$ by pair production (see Fig. 6 in Appendix C). Another source is required for these events Chen et al. (2015); Palladino et al. (2017); Chianese et al. (2017); Denton and Tamborra (2018b); Riabtsev and Troitsky (2023).

Considering collapsars as sources for both HE neutrinos and $r$-process nuclei, we have investigated a novel effect on HE neutrinos caused by their annihilation with LE antineutrinos from the $\beta$-decay of $r$-process nuclei. Assuming typical jet parameters, we have also discussed production of HE neutrinos at internal subshocks induced by mildly-magnetized jets deep inside collapsars. Bright GRBs, LLGRBs, and choked GRBs may have common origin in collapsars, and contributions to the diffuse HE neutrino flux are expected to be dominated by the latter two. Our collapsar model suggests that these sources can well account for the observed flux at IceCube, especially when annihilation with LE antineutrinos from $\beta$-decay of $r$-process nuclei is included. Specifically, an excess at 10–100 TeV and potentially, a deficit at 0.1–1 PeV can be explained consistently for the IceCube events without overproducing the diffuse $\gamma$-ray background (see e.g., Murase et al. (2016, 2020); Capanema et al. (2020); Fang et al. (2022)).

A critical test of our HE neutrino production model is detection of precursor HE neutrinos (with a lead time of $\sim 10$–$10^{3}$ s) from nearby LLGRBs or bright GRBs. For an LLGRB at 100 Mpc with canonical parameters, IceCube-Gen2 can detect $\sim 10$ cascade events within 10–100 s if the jets are mildly magnetized. This short emission would also favor the speculation that LLGRBs host very similar jets to those in bright GRBs. Our model predicts a flavor composition at Earth close to that in the so-called muon-damped scenario at 10 TeV to 1 PeV, which can be distinguished from those in the standard pion and neutron decay scenarios Ackermann (2017); Aartsen et al. (2021); Song et al. (2021) by precise measurements. The expected spectrum decays as $E^{-\gamma}$ at energies below $\sim 100$ TeV with $\gamma\sim 0.5$–1, followed by a steepening above $100$–200 TeV due to the annihilation effect. Anticipating 10 years of data at IceCube-Gen2, we find that a spectral steepening with $\Delta\gamma>1$ for our scenario can be observed at 95% CL (see Appendix D). In contrast, no clear break is expected without annihilation. The above neutrino signatures together with the $r$-process imprints in collapsar lightcurves can provide important support of collapsars as sources for both HE neutrinos and $r$-process nuclei.

We have made several simplifications in this initial study. Future improvements include detailed modeling of LE antineutrino emission, spatial and temporal variations of $r$-process production in collapsars, and self-consistent collapsar simulations for HE neutrino production. Finally, the same annihilation effect should occur in BNSMs that host short GRBs and the $r$-process, and can in principle be tested with a future nearby event.