Ultrahigh-energy cosmic-ray nuclei and neutrinos from engine-driven supernovae

TL;DR

This paper explores how engine-driven supernovae can accelerate nuclei to ultrahigh energies and produce detectable neutrinos, aligning with observed cosmic-ray spectra and offering predictions for future neutrino detection.

Contribution

It presents a detailed model of UHECR acceleration and neutrino production in engine-driven supernovae, including spectra calculations and survival conditions, advancing understanding of their role as cosmic-ray sources.

Findings

UHECR nuclei can reach energies of 10^{20}-10^{21} eV.

Predicted neutrino fluxes are within detection capabilities of future observatories.

The model explains observed cosmic-ray composition and spectrum with a supernova rate of about 1%.

Abstract

Transrelativistic supernovae (SNe), which are likely driven by central engines via jets or winds, have been among candidate sources of ultrahigh-energy cosmic rays (UHECRs). We investigate acceleration and survival of UHECR nuclei in the external reverse shock scenario. With composition models used in Zhang et al. (2018), we calculate spectra of escaping cosmic rays and secondary neutrinos. If their local rate is $\sim1$% of the core-collapse supernova rate, the observed UHECR spectrum and composition can be explained with the total cosmic-ray energy ${\mathcal E}_{\rm cr}\sim10^{51}$ erg. The maximum energy of UHECR nuclei can reach $\sim 10^{20}-{10}^{21}\rm~eV$. The diffuse flux of source neutrinos is predicted to be $\sim 10^{-11} - 10^{-10}~{\rm GeV}~{\rm cm}^{-2}~{\rm s}^{-1}~{\rm sr}^{-1}$ in the 0.1-1 EeV range, satisfying nucleus-survival bounds. The associated cosmogenic neutrino flux is calculated, and shown to be comparable or even higher than the source neutrino flux. These ultrahigh-energy neutrinos can be detected by ultimate detectors such as the Giant Radio Askaryan Neutrino Detector and Probe Of Extreme Multi-Messenger Astrophysics.

Equations

$${\nu }_m$$

$${\nu }_c\simeq 4.1\times 10^{13}{ϵ}_{B,-1.3}^{-3/2}{\mathcal{E}}_{\mathrm{k},51.5}^{-1/2}{\varrho }_{\mathrm{cbm},1}^{-1}{T}_4^{-1/2}\mathrm{Hz},$$

$${\nu }_a$$

$$\frac{dn}{d\epsilon }$$

$$L_{\epsilon ,\max }$$

$${\nu }_m$$

$${\nu }_c$$

$${\nu }_a$$

$$L_{\epsilon ,\max }\simeq 1.2\times 10^{56}{f}_{e,-2}{ϵ}_{\mathrm{B},-0.3}^{1/2}{A}_{\star ,0.7}^{3/2}{\beta }_{\mathrm{fej},-0.5}{\mathrm{s}}^{-1}.$$

$$E_{A,\max }^{\mathrm{jet}}$$

$$E_{A,\max }^{\mathrm{fej}}$$

$${\tau }_{A\gamma }$$

$$f_{A\gamma }$$

$${\epsilon }_{\nu }^2\frac{d{N}_{\nu }}{d{\epsilon }_{\nu }}\approx \frac{3}{8}{f}_{\sup }{f}_{\mathrm{mes}}({\epsilon }_A)(1-f_{A\gamma }({\epsilon }_A)){\epsilon }_A^2\frac{d{N}_A}{d{\epsilon }_A},$$

$${\epsilon }_{\nu }^2\frac{d{N}_{\nu }}{d{\epsilon }_{\nu }}\approx \frac{3}{8}{f}_{\sup }{f}_{p\gamma }({\epsilon }_p)f_{A\gamma }({\epsilon }_A){\epsilon }_A^2\frac{d{N}_A}{d{\epsilon }_A},$$

$${\epsilon }_{\nu }^2\frac{d{N}_{\nu }}{d{\epsilon }_{\nu }}\approx \frac{3}{8}{f}_{\sup }{f}_{p\gamma }({\epsilon }_A/A){\epsilon }_A^2\frac{d{N}_A}{d{\epsilon }_A}.$$

$${\lambda }_{A,Z}^{-1}({\epsilon }_A)\simeq A{\lambda }_{p/n}^{-1}({\epsilon }_A/A),$$

$$E_{\nu }^2{\varphi }_{\nu }$$

$$E_{\nu }^2{\Phi }_{\nu }$$

$$\frac{d{V}_c}{dz}=\frac{c}{H_0}\frac{4\pi d_L^2}{(1+z{)}^2\sqrt{{\Omega }_M(1+z{)}^3+{\Omega }_{\Lambda }}},$$

$${\Phi }_{\nu }(E_{\nu })$$

Full text

††preprint: APS/123-QED

I Introduction

The deaths of massive stars lead to the production of cosmic rays, which have been established by cumulative evidences from studies on supernova (SN) remnants. Ordinary SNe are accompanied by a nonrelativistic shock, at which the diffusive shock acceleration occurs. They could be responsible for observed cosmic rays up to the knee at $\sim 3\times{10}^{15}$ eV or iron knee around $10^{17}$ eV, but the maximum energy cannot reach ultrahigh energies. Among various candidate sources of ultrahigh-energy cosmic rays (UHECRs) Hillas (1984); Kotera and Olinto (2011); Anchordoqui (2018), rarer but faster SNe have been suggested by various authors Murase et al. (2006); Wang et al. (2008); Budnik et al. (2008); Wang et al. (2007); Murase et al. (2008a, 2009); Chakraborty et al. (2011); Liu et al. (2011); Liu and Wang (2012); Zhang et al. (2018) (see Ref. Murase and Takami (2009) for a list of transient sources for UHECRs). Indeed, observations have revealed that some SNe have a transrelativistic component with shock velocity $\beta_{s}\sim 0.1-1$. They include engine-driven SNe associated with low-luminosity (LL) gamma-ray bursts (GRBs), such as GRB 980425 with SN 1998bw Galama et al. (1998); Kulkarni et al. (1998), GRB 060218 with SN 2006aj Soderberg et al. (2006); Campana et al. (2006), GRB 100316D with SN 2010bh Margutti et al. (2013), GRB 161228B with iPTF17cw Corsi et al. (2017), and GRB 171205A with SN 2017iuk D’Elia et al. (2018); Wang et al. (2018); Izzo et al. (2019). Such engine-driven SNe with transrelativistic ejecta have been found through radio observations; examples without GRB counterparts include SN 2009bb Soderberg et al. (2010) (but see also Ref. Nakauchi et al. (2015)) and SN 2012ap Chakraborti et al. (2015); Margutti et al. (2014). Another possible example with low-mass subrelativistic ejecta is AT2018cow, although its origin is under debate Perley et al. (2018); Margutti et al. (2018). It has been shown that the event rate of SN 2009bb-like SNe is about 1% of the core-collapse SN rate, which is not far from the true rate of LL GRBs, $\sim 100-1000\rm~{}Gpc^{-3}~{}yr^{-1}$ Chakraborty et al. (2011); Sun et al. (2015). It is also comparable to the hypernova rate Podsiadlowski et al. (2004); Guetta and Della Valle (2007), although all hypernovae are not necessarily accompanied by transrelativistic ejecta. The diversity of these explosive phenomena may originate from characteristics of central engines and/or progenitors Thompson et al. (2004), and the relativistic velocity component can be caused by energy injections via outflows from the engines Lazzati et al. (2012); Nakar (2015); Suzuki et al. (2017).

The diffusive shock acceleration theory predicts that the acceleration time scale is $\sim E/(ZeBc\beta_{s}^{2})$, where $E$ is particle energy, $Z$ is the charge, and $B$ is magnetic field strength. We can see the acceleration is easier for faster shocks ($\Gamma\gg 1$), but it is known that ultrarelativistic, superluminal shocks are unlikely to be efficient cosmic-ray accelerators Achterberg et al. (2001). In this sense, engine-driven SNe, whose shocks are subrelativistic or mildly relativistic, $\Gamma\lesssim 1-10$, would be more favorable for the acceleration of UHECR nuclei. The magnetic field strength inferred from radio observations Chakraborty et al. (2011) also supports that the acceleration of UHECR nuclei is allowed by the Hillas criterion Hillas (1984). It has been shown that the composition of UHECR nuclei becomes heavier beyond “ankle” ($\sim 4\times 10^{18}\rm~{}eV$) Aab et al. (2017a) which requires their sources to have super-solar abundance Aab et al. (2017b). Contrary to the forward shock (FS) scenario, in the internal shock and reverse shock (RS) scenarios, an outflow or ejecta composition with intermediate mass nuclei enriched is naturally achieved in terms of progenitor models Zhang et al. (2018).

If UHECRs are dominated by intermediate or heavy nuclei, source identification with UHECR observations is more challenging because of larger magnetic deflections. The requirement of nucleus survival also restricts the production of high-energy neutrinos both inside and outside the sources Murase and Beacom (2010a). Detections of neutrinos or electromagnetic counterparts from the UHECR accelerators will provide us with a unique opportunity to reveal the sources, as studied for LL GRBs and transrelativistic SNe Murase et al. (2006); Gupta and Zhang (2007); Murase et al. (2008a); Murase (2009); Murase and Beacom (2010b); Liu et al. (2011); Murase (2012); Kashiyama et al. (2013); Murase and Ioka (2013); Senno et al. (2016); Boncioli et al. (2018); Denton and Tamborra (2018). Other source classes including galaxy clusters Murase et al. (2008b); Fang and Murase (2018), activate galactic nuclei Murase (2017), classical GRBs Waxman (1995); Murase et al. (2008a); Globus et al. (2015); Biehl et al. (2017a), tidal disruption events (TDEs) Zhang et al. (2017); Biehl et al. (2017b); Guépin et al. (2017), and new born pulsars Fang et al. (2014). Neutrinos that directly originate from UHECR nuclei have energy with $E_{\nu}\approx 0.05(E/A)\simeq 0.5~{}{\rm EeV}~{}{(E/{10}^{20}~{}\rm eV)}(10/A)$. This energy range is higher than that of IceCube neutrinos, 0.1-1 PeV. Such extremely high-energy neutrinos are the main targets for planned neutrino detectors in the near future, Askaryan Radio Array (ARA) Allison et al. (2016), Antarctic Ross Ice Shelf Antenna Neutrino Array (ARIANNA) Barwick et al. (2015), Giant Radio Array for Neutrino Detection (GRAND) Alvarez-Muniz et al. (2018), Probe Of Extreme Multi-Messenger Astrophysics (POEMMA) Olinto et al. (2018), and Trinity Otte (2018). The discovery of IceCube-170922A, coinciding with a flaring blazar, TXS 0506+056, demonstrated the feasibility of neutrino-triggered follow-up observations at different wavelengths Aartsen et al. (2018a); Keivani et al. (2018); Ansoldi et al. (2018), although there is no simple picture for the 2017 and 2014-2015 flares Murase et al. (2018). LL GRBs and engine-driven SNe will be interesting targets for such follow-up observations, as proposed by Ref. Murase et al. (2006).

In addition to the source neutrinos, cosmogenic neutrinos which are produced during the propagation of UHECR nuclei in the intergalactic space are believed to be guaranteed possibilities, even though the expected neutrino flux are subject to $\sim 1-2$ orders of magnitude uncertainty depending on the composition, maximum acceleration energy and source redshift evolution Berezinsky and Zatsepin (1969); Stecker (1979); Engel et al. (2001); Takami et al. (2009); Kotera et al. (2010); Alves Batista et al. (2018).

In this work, we provide a comprehensive study of the generation and survival of UHECR nuclei and neutrinos from engine-driven SNe. For general consideration, we discuss two kinds of outflows; one is mildly relativistic jets with $\Gamma\sim 2-10$, and the other is the slower transrelativistic ejecta with $\Gamma\beta\sim 1$. We revisit the “nucleus-survival problem” that was earlier studied by Ref. Murase et al. (2008a); Wang et al. (2008); Murase and Beacom (2010a), and then calculate spectra of escaping cosmic rays and secondary neutrinos both numerically and analytically. Our results demonstrate that engine-driven SNe accompanied by LL GRBs jets Murase et al. (2008a) or transrelativistic ejecta Liu and Wang (2012) give a viable explanation for the UHECR data even in the RS scenario.

This paper is organized as follows. We first study the physics of RS formed by both jets and transrelativistic ejecta in Sec. II. In Sec. III, we discuss the acceleration of UHECR nuclei in the RS scenario and compared to the observation results of UHECR nuclei. In Sec. IV, we estimate neutrinos that are coproduced with UHECR nuclei from jets or winds considering the photomeson production process between UHECR nuclei and ambient radiation fields. We concentrate on neutrino fluences from single source and diffuse neutrinos which takes into account the contribution from all of the events in the universe and cosmogenic neutrinos. We discuss implications of our results in Sec. V and give a summary in Sec. VI.

Throughout of the paper, we adopt the cgs unit and have notations $Q_{x}\equiv Q/10^{x}$. The cosmological parameters we use are $H_{0}=67.3\ \rm km\ s^{-1}\ Mpc^{-1}$, $\Omega_{m}=0.315$, $\Omega_{\Lambda}=0.685$ Olive et al. (2014).

II The physics of reverse shock

After the core collapse of a progenitor star, outflows from the engine, either in the form of jets or winds, will push the ejecta. We expect that such engine-driven ejecta posses an enhanced fast-velocity component Izzo et al. (2019), becoming transrelativistic SNe. The material collides into the surrounding circumburst medium (CBM), leading to the formation of external reverse-forward shocks Mészáros (2006). While the FS plows the CBM directly, the RS propagates back and decelerate the outflow. Most of the RS emission occurs around the shell crossing time $t_{\times}$. If the engine duration $T_{\rm eng}$ is longer than the deceleration time $t_{\rm dec}$, corresponding to the thick shell case, the shell crossing time is approximately given by $\sim T_{\rm eng}$. In more general, it is given by $t_{\times}\sim{\rm max}[T_{\rm eng},t_{\rm dec}]$. The shocked ejecta and shocked CBM are separated by the contact discontinuity (CD) where the pressure equilibrium has been established across the interacting surface, see Fig. 1.

One of the particular features of models for acceleration by engine-driven outflows is that they may be composed by a large fraction of heavy nuclei as we have illustrated in Fig. 1. Heavy nuclei can be extracted from the inner stellar core Horiuchi et al. (2012); Zhang et al. (2018) or synthesized during the expansion of the outflow Nakamura et al. (2001); Liu and Wang (2012); Metzger et al. (2011). In this work, we adopt the nuclear composition model, Si-R I, proposed in Ref. Zhang et al. (2018) as a fiducial model of the jet composition, while the transrelativistic ejecta has a composition similar to the hypernova model in Ref. Zhang et al. (2018).

II.1 Reverse shock by jets

An observed signature of RS emission is the early optical or radio afterglow at the end of the main burst of GRBs Akerlof et al. (1999); Laskar et al. (2013, 2018). The RS emission is important to constrain the initial Lorentz factor as well as the baryonic component in the jets Zhang et al. (2003); Nakar and Piran (2004); McMahon et al. (2006). The RS emission can be suppressed or even missed if the ejecta is dominated by the magnetic energy, but we can expect strong RS emission for baryon dominated jets Zhang and Kobayashi (2005); Kumar and Zhang (2014), which are promising for engine-driven SNe.

The dynamical properties of RS can be more complicated than the self-similar evolution of FS, especially in the presence of long-lasting energy injection Uhm and Beloborodov (2007); Genet et al. (2007). For our purpose, we only consider the shock physics at radius $R_{\times}$, where RS finishes crossing the GRB ejecta shell Panaitescu and Kumar (2004); Murase (2007). After the shock crossing, the RS light curve is expected to decline rapidly with time $\propto t^{-2}$ as a result of the rapid cooling of electrons Zhang et al. (2003); Kobayashi and Zhang (2003). In the following, our calculations are based on the work in Ref. Murase (2007).

We assume that the GRB ejecta have an initial Lorentz factor $\Gamma_{0}$ and isotropic-equivalent ejecta energy $\mathcal{E}_{\rm k}$. The CBM density is given by $\varrho_{\rm cbm}$, and $\varrho_{\rm ej}$ is the number density of the GRB ejecta, which can be estimated by $\varrho_{\rm ej}=\mathcal{E}_{\rm k}/4\pi m_{p}c^{2}\Gamma_{0}(\Gamma_{0}\Delta)R^{2}$ with $\Delta$ as the ejecta shell thickness in the engine frame. For a homogeneous CBM with $\varrho_{\rm cbm}=\rm const$, the RS crossing radius $R_{\times}$ is estimated to be $R_{\times}\simeq 5.6\times 10^{16}\mathcal{E}_{{\rm k},51.5}^{1/4}\varrho_{\rm cbm,1}^{-1/4}T_{4}^{1/4}{\rm~{}cm}$, where we adopt the “thick ejecta shell” case considering $\Delta=cT>R_{\times}/2\Gamma_{0}^{2}$, and $T=10^{4}\rm~{}s$ is the engine frame duration of the GRB ejecta Panaitescu and Kumar (2004). This is justified when the central engine is active for a sufficiently long time. Note that if $R_{\times}/2\Gamma_{0}^{2}>cT$, we should consider the “thin ejecta shell” $\Delta=R_{\times}/2\Gamma_{0}^{2}$, where the thickness of the ejecta shell are dominated by the velocity spreading.

The Lorentz factor of the shocked ejecta in the engine frame is $\Gamma_{\times}\simeq 6.3~{}\mathcal{E}_{{\rm k},51.5}^{1/8}\varrho_{\rm cbm,1}^{-1/8}T_{4}^{-3/8}$, where we adopt the condition $\varrho_{\rm ej}/\varrho_{\rm cbm}\ll 4\Gamma_{0}^{2}$ for more tenuous ejecta. The Lorentz factor of the shocked ejecta viewed from the frame of the unshocked ejecta can be calculated from the addition of velocities in special relativity, $\Gamma_{\times}^{\prime}\approx(1/2)(\Gamma_{\times}/\Gamma_{0}+\Gamma_{0}/\Gamma_{\times})\simeq 1.1$. The magnetic field strength of the shocked GRB ejecta can be estimated assuming a fraction $\epsilon_{B}$ of the post-shock energy density is converted into the magnetic energy, $B_{\times}\simeq 1.6\epsilon_{B,-1.3}^{1/2}\mathcal{E}_{{\rm k},51.5}^{1/8}\varrho_{\rm cbm,1}^{3/8}T_{4}^{-3/8}{\rm~{}G}$.

Once we know the Lorentz factor and magnetic field strength of the shocked ejecta, we can constrain the RS emission spectra. The typical break frequencies measured in the engine frame can be calculated using the formula $\nu_{i}=3e\gamma_{i}^{2}B_{\times}\Gamma_{\times}/4\pi m_{e}c$ with some characteristic Lorentz factor of electrons, $\gamma_{i}$. Here $\nu_{i}$ represents $\nu_{m}$ (injection frequency), $\nu_{a}$ (self-absorption frequency), and $\nu_{c}$ (cooling frequency), respectively. The injection synchrotron frequency in the engine frame is

[TABLE]

with $\epsilon_{e}$ is the equipartition value of the thermal energy convert to electrons, $f_{e}$ is the number fraction of electrons that are accelerated. We adopt $s=2.4$ as the default electron spectral index as in Ref. Murase (2007), and the chosen value $s=2.4$ is already used in previous works in order to reproduce the external reverse-forward shock emission Meszaros and Rees (1999); Panaitescu and Kumar (2004). The electron cooling Lorentz factor depends on the ratio between electron radiation time scale and dynamical time scale $\gamma_{c}=6\pi m_{e}c^{2}\Gamma_{\times}/\sigma_{T}(Y+1)R_{\times}B_{\times}^{2}$, where $Y$ is the Compton Y parameter. The typical cooling frequency in the slow cooling regime is

[TABLE]

and the self-absorption frequency is

[TABLE]

The latter is estimated by setting the self-absorption optical depth $\tau(\nu_{a})$ to unity Panaitescu and Kumar (2004); Murase (2007).

The synchrotron emission from RS can be described as broken power law Murase (2007) ($\nu_{a}<\nu_{m}<\nu_{c}$)

[TABLE]

where $n_{\varepsilon,\rm~{}max}=L_{\varepsilon,\rm~{}max}/4\pi R_{\times}^{2}c\varepsilon_{m}$ is the normalization of the differential photon number density. The comoving frame luminosity per unit energy is

[TABLE]

where $N_{e}=\mathcal{E}_{\rm k}/\Gamma_{0}m_{p}c^{2}$. We show the comoving frame differential photon number density (blue lines) in Fig. 2, which are calculated from following different parameter sets:

• •

Jet-A: $\mathcal{E}_{\rm k}=3\times 10^{51}\rm~{}erg$, $T=10^{4}\rm~{}s$, $\Gamma_{0}=10$, $\varrho_{\rm cbm}=10\rm~{}cm^{-3}$, $\epsilon_{e}=0.1$, $f_{e}=0.01$, $\epsilon_{B}=0.01$, and $s=2.4$.

• •

Jet-B: $\mathcal{E}_{\rm k}=3\times 10^{51}\rm~{}erg$, $T=10^{4}\rm~{}s$, $\Gamma_{0}=10$, $\varrho_{\rm cbm}=1\rm~{}cm^{-3}$, $\epsilon_{e}=0.1$, $f_{e}=0.01$, $\epsilon_{B}=0.01$, and $s=2.4$.

• •

Jet-C: $\mathcal{E}_{\rm k}=4\times 10^{51}\rm~{}erg$, $T=10^{4}\rm~{}s$, $\Gamma_{0}=5$, $\varrho_{\rm cbm}=10000\rm~{}cm^{-3}$, $\epsilon_{e}=0.001$, $f_{e}=0.01$, $\epsilon_{B}=0.001$, and $s=2.1$.

The detailed value of $\epsilon_{e}$, $f_{e}$, and $\epsilon_{B}$ depends on the microphysics of the collisionless shocks, and which are still unclear from first principles Eichler and Waxman (2005). The main difference of Jet-A and Jet-B are the density of external medium. We also consider model Jet-C, which can give a well-fit to the radio afterglow of GRB 060218 Toma et al. (2007), as shown in Fig. 7. We also note that our model does not overshoot either the optical or x-ray data observed for GRB 060218.

II.2 Reverse shock by transrelativistic ejecta

The origin of transrelativistic SNe has been under debate but a natural possibility is that they originate from outflows or jets. The distribution of kinetic energy has a velocity dependence, $\mathcal{E}_{\rm ej}(>\Gamma_{\rm ej}\beta_{\rm ej})\propto(\Gamma_{\rm ej}\beta_{\rm ej})^{-\alpha}$, where $\beta_{\rm ej}=v_{\rm ej}/c$ is the normalized velocity and $\Gamma_{\rm ej}=1/(1-\beta_{\rm ej}^{2})^{1/2}$ is the Lorentz factor. The profile is usually steep $\alpha\sim 5$ for ordinary SNe Tan et al. (2001), whereas GRBs usually have a much flatter profile with $\alpha\simeq 0.4$ Lazzati et al. (2012). SNe associated with GRBs, including LL GRBs, have intermediate properties with $\alpha\sim 2-3$. This implies that most of the kinetic energy is taken by low-velocity ejecta of these SNe Margutti et al. (2013). For example, let us consider hypernovae with a total kinetic energy of ${\mathcal{E}}_{\rm ej}\sim 3\times 10^{52}$ erg and a typical velocity of $\beta_{\rm ej}\sim 0.1$ (corresponding to an ejecta mass of $M_{\rm ej}\sim 3~{}M_{\rm sun}$). Taking into account weights of the ejecta velocity distribution, the fast ejecta component has a typical velocity of $\beta_{\rm fej}\sim 0.3$ with its characteristic energy of ${\mathcal{E}}_{\rm fej}\sim 3\times{10}^{51}$ erg. In reality, the faster ejecta are decelerated earlier by the CBM. However, for simplicity, we treat the fast velocity ejecta as a single velocity component, which is sufficient for the demonstrative purpose of this work.

The CBM is assumed to have a wind-like density profile, $\varrho_{\rm cbm}=Ar^{-2}$, where $r$ is the radial distance, $A\equiv 3\times 10^{35}A_{\star}{\rm~{}cm}^{-1}$, and $A_{\star}$ is the mass-loss rate to wind speed ratio normalized to $10^{-5}\rm~{}M_{\odot}~{}yr^{-1}/10^{3}\rm~{}km~{}s^{-1}$ Panaitescu and Kumar (2004). Then, the deceleration radius for the fast velocity ejecta is estimated to be $R_{\rm dec}\simeq 2.1\times 10^{18}\mathcal{E}_{\rm fej,51.5}A_{\star,0.7}^{-1}\beta_{\rm fej,-0.5}^{-2}\rm~{}cm$. Here, we give a rough approximation for the magnetic field strength as $B_{\rm rs}\approx(\epsilon_{\rm B}4\pi\varrho_{\rm cbm}\beta_{\rm fs}^{2}c^{2})^{1/2}$, where we have used the energy density balance, i.e., $\varrho_{\rm fej}\beta_{\rm rs}^{2}=\varrho_{\rm cbm}\beta_{\rm fs}^{2}$. Note that the reverse shock velocity is measured in the ejecta rest frame. Then we can derive the characteristic synchrotron frequency in the engine frame,

[TABLE]

where we assume the accelerated electrons following a power-law distribution with spectral index $s=2.4$ and $\gamma_{m}=((s-2)/(s-1))(\epsilon_{e}/f_{e})(m_{p}/m_{e})\beta_{\rm rs}^{2}$. Similarly, the electron cooling frequency is

[TABLE]

and the synchrotron self-absorption frequency is

[TABLE]

The radiation luminosity per unit energy at characteristic frequency can be estimated to be

[TABLE]

We show the differential photon number density in Fig. 2 which is estimated using following parameter sets,

TRSN-A: $\mathcal{E}_{\rm fej}=3\times 10^{51}\rm~{}erg$, $\beta_{\rm fej}=0.8$, $A_{\star}=1$, $\epsilon_{e}=0.1$, $f_{e}=0.01$, $\epsilon_{B}=0.1$ and $s=2.4$.

TRSN-B: $\mathcal{E}_{\rm fej}=3\times 10^{51}\rm~{}erg$, $\beta_{\rm fej}=0.3$, $A_{\star}=5$, $\epsilon_{e}=0.1$, $f_{e}=0.01$, $\epsilon_{B}=0.5$ and $s=2.4$.

III UHECR nuclei from engine-driven SNe

III.1 Acceleration and survival of UHECR nuclei

Charged nuclei can be accelerated to higher energies through the first-order Fermi acceleration mechanism Bell (1978); Blandford and Eichler (1987); Waxman and Bahcall (1998). The maximum acceleration energy is determined by the condition that the acceleration time scale should be smaller than dynamical time scale $t_{\rm acc}\leq t_{\rm dyn}$, as well as various energy cooling time scales, $t_{\rm acc}\leq t_{\rm cool}\equiv 1/(t_{\rm ad}^{-1}+t_{A\gamma}^{-1}+t_{\rm syn}^{-1}+t_{\rm IC}^{-1})$. Here $t_{\rm acc}$ is the acceleration time, $t_{\rm dyn}$ is the dynamical time, $t_{\rm cool}$ is the total energy cooling time, $t_{\rm ad}\sim t_{\rm dyn}$ is the adiabatic cooling time, $t_{\rm syn}$ is the synchrotron cooling time, $t_{\rm IC}$ is the inverse Compton cooling time, and $t_{A\gamma}$ is the photohadronic cooling time Zhang et al. (2017); Murase (2007). In the case of the age-limited acceleration ($t_{\rm cool}\geq t_{\rm dyn}$), the engine frame maximum energy of nuclei accelerated in jets can be estimated as,

[TABLE]

where $\eta\geq 1$ is the acceleration efficiency Murase and Nagataki (2006), $Z$ is the nucleus charge number, $R_{\times}/\Gamma_{\times}$ is the comoving width of the shocked ejecta. While the maximum energy of nuclei accelerated in transrelativistic ejecta is estimated to be,

[TABLE]

We see that the acceleration energy of heavy nuclei can exceed $\sim 10^{20}\rm~{}eV$ in both jets and transrelativistic ejecta. Note that the maximum energy is sensitive to the ejecta velocity, $E_{A,\rm~{}max}^{\rm fej}\propto\beta_{\rm fej}^{2}$ in the latter case, and the acceleration of UHECR nuclei to ultrahigh energies is much difficult for lower-velocity ejecta. Note that most of the ejecta kinetic energy is taken by ejecta with relatively lower velocities, which may make the situation even worse Wang et al. (2007). Considering the diversity of engine-driven SNe and the uncertainty of the ejecta profile Lazzati et al. (2012), we simply assume that high-velocity ejecta ($\beta_{\rm ej}\geq 0.3$) takes a significant fraction of the total outflow energy. The situation can also be alleviated by denser CBM and/or stronger magnetic field amplification.

Another requirement is that the energy budget of UHECR nuclei from engine-driven SNe should explain the observation results. The total energy of CR nuclei per event is $\mathcal{E}_{\rm CRacc}=6\times 10^{50}\xi_{\rm CRacc,-0.7}\mathcal{E}_{\rm k,51.5}\rm~{}erg$ where $\mathcal{E}_{\rm k}=3\times 10^{51}\rm~{}erg$ is the kinetic energy of the jet or transrelativistic ejecta and $\xi_{\rm CRacc}=1/5$ is the energy fraction of accelerated CRs. Then, we can derive the energy injection rate density, $Q_{\rm UHECR}=\mathcal{E}_{\rm CResc}\dot{\rho}\simeq 2\times 10^{44}\mathcal{E}_{\rm CResc,-1}\dot{\rho}_{3.5}\rm~{}erg~{}Mpc^{-3}~{}yr^{-1}$, where $\dot{\rho}$ is the local event rate of engine-drivne SNe. We found the energy injection rate ensity $Q_{\rm UHECR}$ is comparable to the required energy budget of UHECR nculei Murase and Takami (2009); Katz et al. (2009); Zhang et al. (2018).

In Fig. 4, we show the energy spectrum of UHECR nuclei from engine-driven SNe. The escaping cosmic rays show a very hard energy spectrum if we adopt the escape-limited acceleration scenario under the condition $t_{\rm esc}\leq t_{\rm dyn}$. We can see that cosmic rays around maximum energy $E_{A,\rm max}^{\rm(esc)}\simeq x_{\rm esc}E_{A,\rm max}^{\rm(acc)}\simeq 10^{18.3}x_{\rm esc,-0.2}\rm~{}eV$ can escape efficiently, where $x_{\rm esc}$ is the ratio of the escape boundary to the width of the shocked region Ohira et al. (2010); Zhang et al. (2018). The maximum escape energy adopted in this work can be achieved assuming the acceleration efficiency is $\eta=8$ for the jet and $\eta=1$ for the transrelativistic ejecta.

The dominant energy cooling process of UHECR nuclei is photodisintegration. Using the giant dipole resonance (GDR) approximation, the optical depth of UHECR nuclei which is defined as $\tau_{A\gamma}\equiv t_{\rm dyn}/t_{A\gamma-{\rm int}}$ can be estimated analytically Murase and Beacom (2010a)

[TABLE]

where we adopt model Jet-A and Oxygen nuclei with $\sigma_{\rm GDR}=3.1\times 10^{-27}\rm~{}cm^{2}$, $\Delta\bar{\varepsilon}_{\rm GDR}=6\rm~{}MeV$, and $\bar{\varepsilon}_{\rm GDR}=22.35\rm~{}MeV$ Alves Batista et al. (2016). The target photon number density is expressed as $dn/d\varepsilon=n_{0}(\varepsilon/\varepsilon_{0})^{-\alpha}$, where $\alpha=2/3$ is the spectral index. The effective optical depth is defined as $f_{A\gamma}\equiv t_{\rm dyn}/t_{A\gamma}$, where $t_{A\gamma}$ is the energy loss time scale of UHECR nuclei. The relation between effective optical depth $f_{A\gamma}$ and optical depth $\tau_{A\gamma}$ can be expressed as Murase and Beacom (2010a)

[TABLE]

where $\bar{\kappa}_{A\gamma}$ is the inelasticity which represents average energy losses of the photodisintegration process Zhang et al. (2017). We can see that UHECR nuclei are in the partial survival regime for model Jet-A which satisfy the condition $\tau_{A\gamma}>1$ and $f_{A\gamma}<1$ Zhang et al. (2017, 2018) (see Fig. 5). The complete survival of UHECR nuclei is possible when $\tau_{A\gamma}<1$ for model Jet-B (see Fig. 6).

We also show the survival of UHECR nuclei for model TRSN-A (partial survival regime) in Fig. 8 and TRSN-B (complete survival regime) in Fig. 9, respectively. We can see that the survival of UHECR nuclei is possible in both jets and transrelativistic ejecta for various parameter sets. One limitation of our study is that in our calculation we assume an ISM-like density profile of CBM for jets and a wind-like density profile for transrelativistic ejecta. Details of the survival are affected by the ambient density. For jets, the wind-like density profile is dangerous for the conventional value $A_{\star}=1$ because the density becomes too high. Our results are not sensitive to the specific value of electron spectral index, even though it may slightly affect the electron injection frequency.

III.2 Comparison with Auger and TA data

UHECRs have been observed on Earth by various experiments Kotera and Olinto (2011); Anchordoqui (2018), such as the Pierre Auger Observatory (Auger) Aab et al. (2015) and Telescope Array (TA) collaboration Abu-Zayyad et al. (2013), which are the two largest cosmic ray detectors at present. It has been shown that the energy spectrum of UHECR nuclei becomes flatter at energy around “ankle” $\sim 4\times 10^{18}\rm~{}eV$, with a suppression at $\sim 5\times 10^{19}\rm~{}eV$. However, there may be a significant statistical discrepancy in the energy spectra between Auger and TA around the suppression energy. The TA results show a relatively higher flux, and the reason is still unclear yet The (2018).

In this section, we aim to fit the energy spectrum and composition of UHECR nuclei detected by Auger Aab et al. (2017a) and TA Fukushima (2015); Abbasi et al. (2018). For TA data, we decrease the flux of UHECR nuclei by a fraction of $\sim 13\%$ in order to be compatible with the Auger data around “ankle” region. Even though it is difficult to make a direct comparison of the mean and standard deviation of the measured $X_{\rm max}$ of Auger Aab et al. (2017a) and TA Abbasi et al. (2018) directly, which may be affected by the particular treatment of biases and detector efficiencies The (2018), we show the results derived by Auger and TA together for demonstrative purposes.

In this work, we propagate UHECR nuclei using the state-of-the-art Monte Carlo code CRPropa 3 Alves Batista et al. (2016). For simplicity, we consider $1\rm D$ propagation, where the sources are uniformly distributed in the universe at each redshift bin and neglect the effects of large-scale magnetic fields. We adopt a semi-analytic extragalactic background light (EBL) model developed by Ref. Gilmore et al. (2012). This model is based on the hierarchical structure formation scenario, taking into account galaxy formation and evolution. We performed $\chi^{2}$ analysis as in Ref. Zhang et al. (2018); Fang and Murase (2018) to show the quality of our fitting. The spectrum energy fitting range is from $10^{18.45}~{}\rm eV$ to $10^{21.15}~{}\rm eV$, where the uncertainties on the flux are taken from the measurements of Auger Aab et al. (2017a) and TA Fukushima (2015). As in Ref. Heinze et al. (2016), we take into account the effect of systematic uncertainty on the energy scale, $E=E^{\rm Auger}(1+\delta_{E})$. We derived the first two moments of the $X_{\rm max}$ distributions using the EPOS-LHC hadronic interaction model De Domenico et al. (2013); Zhang et al. (2017), and compared to the measurements from Auger Zhang et al. (2018) and TA Abbasi et al. (2018), respectively.

The results are shown in Fig. 10-11, where we consider the complete survive regime, Jet-B and TRSN-B. The change to a lighter composition in the $X_{\rm max}$ distributions is caused by the photodisintegration process of nuclei from nearest sources, as explained in Ref. Zhang et al. (2018). We find model Jet-B can fit the Auger data very well, $\chi^{2}/\rm d.o.f.=1.5$, where the best-fitting parameters is $E_{p,\rm max}^{\rm esc}=18.3$ and $\delta_{E}=0.13$, as shown in Fig. 10. On the other hand, model TRSN-B gives the best fit to TA data, $\chi^{2}/\rm d.o.f.=0.9$, where the best-fitting parameters are $E_{p,\rm max}^{\rm esc}=18.3$ and $\delta_{E}=0.01$, see Fig. 11. The difference between Jet-B and TRSN-B is the composition, in which the latter has a significant fraction of iron-group nuclei. Compared to intermediate mass group nuclei, iron-group nuclei can be accelerated to higher energies, with longer energy loss lengths. Although the discrepancy in the energy spectra between the Auger and TA data is still at the $\sim 3\sigma$ level, our results imply that it could in principle be caused by more engine-driven SNe (with more iron nuclei) in the northern sky within the local universe Globus et al. (2017). The injection of heavier nuclei depend on the angular momentum of the progenitor core. Canonical GRB preferentially occur in low-metallicity environments, but LL GRBs and engine-driven SNe may largely occur in starburst galaxies that are often metal polluted.

Note that both Jet-B and TRSN-B models essentially correspond to propagation-only models in Ref. Boncioli et al. (2018). In this work, we also consider the partial survival regime, where only a fraction of nuclei can survive and the energy spectrum of escaping UHECR nuclei may be affected by the nuclear cascade effect. However, even for Jet-A and TRSN-A models, it turns out that the final energy spectrum of UHECR nuclei is not too much affected by the nuclear cascade process if we adopt the parameters listed in the end of Sec. II.1 and Sec. II.2.

IV Neutrinos from engine-driven SNe

The total cross section of the photomeson production process can be decomposed into three parts: $\Delta$-resonance region at peak energy $\bar{\varepsilon}_{\Delta}\simeq 340\rm~{}MeV$ (with a contribution from direct production), higher resonances with two or three peaks in the energy range $500{\rm~{}MeV}<\bar{\varepsilon}<1.5{\rm~{}GeV}$, multipion production region at much higher energies ($>1.5{\rm~{}GeV}$) (see, e.g., Ref. Rachen (1996) for details). The box approximation is commonly used, in which $\sigma_{p\gamma}\approx\sigma_{\Delta}\delta(\bar{\varepsilon}-\bar{\varepsilon}_{\Delta})\Delta\bar{\varepsilon}_{\Delta}$. Here $\sigma_{\Delta}\sim 4.4\times 10^{28}{\rm~{}cm^{2}}$, $\bar{\varepsilon}_{\Delta}\sim 0.34{\rm~{}GeV}$, and $\Delta\bar{\varepsilon}_{\Delta}\sim 0.2{\rm~{}GeV}$ Murase and Beacom (2010a). Note that the branching ratio between neutron pion $\pi^{0}$ production channel and charged pion $\pi^{\pm}$ production channel can effectively approach to $1:1:1$ owing to multipion production Rachen (1996); Murase and Nagataki (2006). The final neutrino energy spectrum from charged pion and muon decay can be calculated numerically through the code SOPHIA Mucke et al. (2000) or GEANT4 Murase and Nagataki (2006).

For nuclei, we assume that the cross section on nuclei scales linearly with their mass, $\sigma_{\rm A\gamma}(\bar{\varepsilon})=A\sigma_{\rm p\gamma}(\bar{\varepsilon})$ and the inelasticity is assumed to have $\kappa_{\rm A\gamma}(\bar{\varepsilon})=(1/A)\kappa_{\rm p\gamma}(\bar{\varepsilon})$. Note that the cross section of the photomeson production on nuclei can be different from the simple scaling relations if we consider the more complicated Fermi motion of nucleons and other inmedium modifications Rachen (1996); Kampert et al. (2013); Biehl et al. (2018).

Considering the photodisintegration of UHECR nuclei, the generated neutrinos can be divided into two parts, one is the direct neutrino contribution due to the photomeson production on nuclei, and another is the indirect neutrino contribution by secondary protons/neutrons (that are stripped off from their parent nuclei due to the photodisintegration process) Murase and Beacom (2010a).

In the first case, the all-flavor neutrino energy spectrum can be simply estimated using the following formula Murase and Beacom (2010a),

[TABLE]

where $f_{\rm sup}$ is the suppression factor due to meson and muon cooling, $f_{\rm mes}(\varepsilon_{A})$ is the effective optical depth for the photomeson production, and ($1-f_{A\gamma}$) represents the fraction of survival nuclei. The factor $3/8$ is due to the fact that half of the energy $\sim 1/2$ goes into charged pions and the pion decay convert $\sim 3/4$ of their energy into neutrinos. Assuming each nucleon has energy $\sim 1/A$ of the parent nucleus with $\varepsilon_{A}$, then the effective optical depth to the photomeson production on nuclei can be approximated by $f_{\rm mes}(\varepsilon_{A})\approx f_{p\gamma}(\varepsilon_{A}/A)$. The typical neutrino energy is $\varepsilon_{\nu}\approx 0.05\varepsilon_{A}/A$. While in the second case, the all-flavor neutrino energy spectrum can be calculated as Murase and Beacom (2010a)

[TABLE]

where $\varepsilon_{p}\approx\varepsilon_{A}/A$ and $f_{A\gamma}$ represents the fraction of photodisintegrated nuclei. In this case, the typical neutrino energy depends on the secondary proton energy $\varepsilon_{\nu}\approx 0.05\varepsilon_{p}$. Note that the secondary protons/neutrons take a fraction $\sim f_{A\gamma}(\varepsilon_{A})$ of the parent nuclei total energy. The two cases can be combined into one formula

[TABLE]

under the approximation $f_{\rm mes}(\varepsilon_{A})\simeq f_{p\gamma}(\varepsilon_{A}/A)$. The suppression factor accounts for the energy cooling on the primary pions and secondary muons, $f_{\rm sup}\simeq f_{\rm sup,\pi}f_{\rm sup,\mu}$ Kimura et al. (2017a). The pion cooling suppression factor is calculated as $f_{\rm sup,\pi}=1-{\rm exp}(-t_{\pi,\rm cool}/t_{\pi,\rm decay})$, where $t_{\pi,\rm decay}$ is the decay time scale of pions. The mean life time for a pion to decay is $\tau_{\pi}=2.6\times 10^{-8}\rm~{}s$ in the pion rest frame, then the pion decay time scale in the shock comoving frame is $t_{\pi,\rm decay}=(\varepsilon_{\pi}/m_{\pi})\tau_{\pi}$. The dominant cooling process of pions are synchrotron cooling, $t_{\pi,\rm cool}^{-1}\sim t_{\pi,\rm syn}^{-1}$. The muon cooling suppression factor can be calculated using the same method.

We also perform numerical simulations in order to take into account the nuclear cascade process inside the source with the help of the publicly available Monte-Carlo code CRPropa Kampert et al. (2013); Alves Batista et al. (2016), where the photomeson production process is calculated using the SOPHIA package Mucke et al. (2000). In CRPropa, the mean free path of nuclei have the following scaling relations Kampert et al. (2013),

[TABLE]

where $\lambda_{p}$, $\lambda_{n}$ are mean free paths for protons and neutrons, separately. The scaling relations equivalent to that there is only one of the nucleons in the parent nuclei actually interact with the target photons as we have declared before where the remained nuclei have mass $(A-1)$ and energy $\varepsilon_{A-1}=(A-1)/A\varepsilon_{A}$ assuming the conservation of the Lorentz factor. For our purpose, we inject nuclei from the original point and propagate them following one direction until nuclei and the secondaries reach to the boundary of the acceleration site. In the numerical simulations, we avoid the synchrotron cooling on the intermediate products, such as pions and muons Murase and Nagataki (2006).

IV.1 Neutrino fluences

The observed neutrino flucences can be estimated using the following formula Murase (2007),

[TABLE]

where $d_{L}$ is the source luminosity distance, $\Gamma$ is the source Lorentz factor and ${\varepsilon^{2}_{\nu}}dN_{\nu}/d\varepsilon_{\nu}$ is the comoving frame neutrino energy spectrum inside the source.

In Fig. 12, we show the results of all-flavor neutrino fluences observed at Earth (red line) emitted from one engine-driven SNe located at redshift $z=0.005$. The neutrino fluence is normalized assuming the total energy of CR nuclei per event is $\mathcal{E}_{\rm CRacc}=6\times 10^{50}\xi_{\rm CRacc,-0.7}\mathcal{E}_{\rm k/fej,51.5}\rm~{}erg$, as shown in Fig. 4. The partial survival regime corresponding to model Jet-A (or TRSN-A) is assumed. The effect of synchrotron cooling on intermediate pions and muons is shown (blue line). We can see that the synchrotron cooling is important in the partial survival regime which can suppress the neutrino fluences at highest energy range. The break energy due to pion/muon cooling can be estimated when $t_{\pi,\rm cool}=t_{\pi,\rm decay}$, $\varepsilon_{\nu,\pi}=2.7\times 10^{18}B_{1}^{-1}\rm~{}eV$ and $\varepsilon_{\nu,\mu}=1.5\times 10^{17}B_{1}^{-1}\rm~{}eV$, where $t_{\pi,\rm dec}$ is the decay time scale of pions or muons. We also show the energy spectrum of direct neutrinos (dashed line) and indirect neutrinos (dot-dashed line), respectively. We find that in the regime of partial nucleus survival, the contribution from indirect neutrinos can have a dominant contribution to the observed neutrino fluences at higher energies Murase and Beacom (2010a); Biehl et al. (2017a).

Similar to Fig. 12, we show the all-flavor neutrino fluences estimated assuming complete survival regime corresponding to model Jet-B (or TRSN-B) in Fig. 13. We can see that the contribution to the neutrino fluences from indirect neutrinos can be neglected and the predicted neutrino fluences in the complete survival regime are about $\sim 1-2$ orders lower than in the partial survival regime.

We also present the neutrino fluences estimated through the numerical simulation which are indicated as black curve in Fig. 12 and Fig. 13. The neutrino fluences derived from numerical simulations which take into account the effect of nuclear cascade process are comparable to our analytic results. The difference at the higher energy range is due to the fact that the analytical approach does not take into account the contribution from multi-pion neutrino production as well as the energy spread of secondaries Mucke et al. (2000); Alves Batista et al. (2019). (The similar trend was found by Ref. Murase (2007).) In the following, we will adopt the results calculated using our analytic formula which takes into account the effect of synchrotron cooling on intermediate pions and muons.

In Fig. 14, we show all-flavor neutrino fluences considering cooling effect derived in this work for one engine-driven SNe located at redshift $z=0.005$ and the expected detection sensitivity from GRAND project Alvarez-Muniz et al. (2018). We considered five models, Jet-A/B/C and TRSN-A/B, and we can see that the neutrino fluence depends on the composition and survival regime of UHECR nuclei. We found that the neutrino fluence from engine-driven SNe is significantly lower than the detection sensitivity of GRAND, which makes it challenging to be detected as a point source of UHE neutrinos.

IV.2 Diffuse source neutrinos

Considering the higher event rate of engine-driven SNe, it is reasonable to consider their contributions to the diffuse neutrino background, including both neutrinos produced inside sources (diffuse source neutrinos) and cosmogenic neutrinos. Note that the normalization of the diffuse neutrinos is fixed under the condition that the escaped UHECR nuclei can reconcile with the observed energy spectrum and composition of UHECR nuclei measured by Auger simultaneously. The flux of diffuse source neutrinos from engine-driven SNe can be estimated using the following formula,

[TABLE]

where $E_{\nu}^{2}\phi_{\nu}$ is the neutrino fluence derived in the previous section. The event rate of engine-driven SNe in the local universe is $\dot{\rho_{0}}$ and $F(z)$ is the redshift distribution parameter estimated from long GRBs Sun et al. (2015). The redshift dependence of the comoving volume is

[TABLE]

where $d_{L}$ is the luminosity distance to the source.

We show the energy spectrum of diffuse source neutrinos in Fig. 15. The thick (thin) blue line represents diffuse neutrinos estimated from model Jet-A (Jet-B), while the thick (thin) red line represents diffuse neutrinos estimated from model TRSN-A (TRSN-B). For comparison, we also show the well-known Waxman-Bahcall bound (black line) assuming fast-evolution scenario ($\xi_{z}\simeq 3$ in Ref. Waxman and Bahcall (1998)) and photodisintegration bound (olive line) for pure Silicon composition Murase and Beacom (2010a). Note that in the non-evolution case, the diffuse neutrino fluxes can be $\sim 5$ times lower than in the fast-evolution case.

IV.3 Cosmogenic neutrinos

The flux of cosmogenic neutrinos can be estimated using the following formula,

[TABLE]

where $\eta_{A^{\prime}\nu}$ is the neutrino yield function which reflects the fraction of neutrinos with energy $E_{\nu}$ originated from nuclei with mass number $A^{\prime}$ and energy $E^{\prime}$ at redshift $z$. The value of $\eta_{A^{\prime}\nu}$ can be estimated numerically using CRPropa 3 Alves Batista et al. (2016), where UHECR nuclei are propagated through the universe as in Sec. III.2. We show the results in Fig. 16 where the flux of cosmogenic neutrinos (green line) can reach a level of a few $\times 10^{-10}\rm~{}GeV~{}cm^{-2}~{}s^{-1}~{}sr^{-1}$ and have peak energy around $0.1-1\rm~{}EeV$. Note the small bump appeared in the lower energy range $\sim\rm~{}PeV$ is due to the effect of neutron beta decay. We can see that the planned neutrino detector GRAND can reach the required sensitivity to detect cosmogenic neutrinos predicted in this work after $\sim 10\rm~{}yr$ observations Alvarez-Muniz et al. (2018). For comparison, we also show the energy spectrum of cosmogenic neutrinos predicted from different candidate sources, including radio galaxies in the shear acceleration scenario (cyan dashed line) Kimura et al. (2017b) and tidal disruption events (TDEs) where a white dwarf is disrupted by an intermediate mass black hole (blue dotted line) Zhang et al. (2017). The flux of cosmogenic neutrinos is sensitive to the source redshift evolution because most of the detected neutrinos are produced at redshift $z\sim 1$, while the observed UHECR nuclei mainly originated from sources within the local universe $\leq 100\rm~{}Mpc$. Engine-driven SNe have a fast redshift evolution, which traces the star formation history. On the other hand, the number density of TDEs usually has a negative redshift evolution Sun et al. (2015).

The main results derived in this work are summarized in Fig. 17, where the predicted energy spectra of UHECR nuclei and neutrinos from engine-driven SNe are shown together.

V Discussion

V.1 Common origin of IceCube neutrinos and UHECRs?

The flux of the observed diffuse neutrinos within energy range from $\sim 100\rm~{}TeV$ to a few PeV by IceCube Aartsen et al. (2014) is comparable to the Waxman-Bahcall (WB) bound for a spectral index of $s=2.0$ Waxman and Bahcall (1998) and the nucleus-survival bound for a spectral index of $s\sim 2.2-2.3$ Murase and Beacom (2010a). The fact that the diffuse energy fluxes of all three messengers are similar has led to the development of the unification picture Murase and Waxman (2016); Fang and Murase (2018).

The flux of diffuse source neutrinos predicted in this work is $\sim 2-3$ orders lower than the observed flux of IceCube TeV-PeV neutrinos. However, it is quite plausible that UHECR nuclei and TeV-PeV neutrinos are produced at different regions within the same source. For example, in the hybrid “two-zone” model Zhang et al. (2018), it has been suggested that UHECR nuclei mainly come from larger radii, where the survival is easy, whereas TeV-PeV neutrinos are produced efficiently near or under the photosphere Murase et al. (2006); Kashiyama et al. (2013); Murase and Ioka (2013); Senno et al. (2016). Indeed, the mechanism of prompt emission from LL GRBs may be different from that of canonical high-luminosity GRBs. The most popular explanation is shock breakout emission of transrelativistic SNe Campana et al. (2006); Nakar (2015). Transrelativistic SNe naturally originate from choked jets Nakar (2015); Senno et al. (2016), and high-energy neutrino emission from low-power choked jets can explain the IceCube data even in the 10-100 TeV range Murase and Ioka (2013); Senno et al. (2016). Our RS model studied in this work is consistent with this shock breakout scenario, and the common origin requires the two-zone picture. Note that our model is also consistent with late time observations at optical wavelengths. As earlier calculated by Ref. Murase et al. (2006), the successful jet scenario Toma et al. (2007); Ghisellini et al. (2007) may also be viable, in which prompt emission comes from the jet that breaks out of the star. Following this specific scenario, a recent paper by Ref. Boncioli et al. (2018) attempted to simultaneously explain the PeV neutrino flux and UHECR data.

V.2 Maximum acceleration energy?

We showed that the maximum energy of UHECR nuclei reaches $\sim{\rm a~{}few}\times(Z/26)10^{20}\rm~{}eV$ in the framework of RS formed by both jets and transrelativistic ejecta. We confirmed the findings of the original papers Murase et al. (2006, 2008a), which found that the maximum energy can be as high as $\sim{\rm a~{}few}\times 10^{20}$ eV for LL GRBs like GRB 06218 without violating observational data (see also Refs. Liu et al. (2011); Boncioli et al. (2018); Zhang et al. (2018)). A recent paper by Ref. Samuelsson et al. (2018) claimed that it may be difficult for protons or iron nuclei to reach the maximum acceleration energy $\sim 10^{20}\rm~{}eV$, if one considers the synchrotron modeling of observed prompt emission for GRBs. Note that even in that work, if one chooses the fiducial parameter sets based on the observations of GRB 060218, $L_{\gamma}\sim 10^{47}\rm~{}erg~{}s^{-1}$ and $\Gamma\sim 1-10$, the maximum energy can exceed $\sim 10^{20}\rm~{}eV$. Also, the mechanism of prompt emission from LL GRBs is yet uncertain, and the constraints on synchrotron emission do not directly applied in the shock breakout scenario.

VI Summary

In this work, we investigated UHECR nuclei and neutrinos originating from engine-driven SNe. In Sec. II, we showed that the acceleration and survival of UHECR nuclei is possible by RS formed by both jets and transrelativistic ejecta. We considered five models, Jet-A/B/C and TRSN-A/B, which are different by the composition and survival regime of UHECR nuclei. The models Jet-A/B/C have a composition similar to Si-R I in Ref. Zhang et al. (2018) and models TRSN-A/B have a hypernova composition with heavier iron-group nuclei. Both of the models Jet-B and TRSN-B are in the complete survival regime, and we found Jet-B is compatible with the Auger data, including the energy spectrum and composition. On the other hand, model TRSN-B can give a better fit to the TA data. It implies that the flux discrepancy between Auger and TA data around the suppression energy could in principle be caused by the enhanced number of engine-driven SNe in the northern sky. However, we should keep in mind that the survival of nuclei is sensitive to the ambient density, as in model Jet-A/C and TRSN-A, which may affect the composition of escaped UHECR nuclei. In Sec. III, we estimated neutrinos that are coproduced with UHECR nuclei. At first, we calculated the received neutrino fluences at Earth from a single engine-driven SN, which is located at distance $z=0.005$. Then we found that the detection is challenging even for planned neutrino detectors, such as GRAND. We also showed that the neutrino fluences can be $\sim 1$ orders higher in the partial survival regime than in the complete survive regime. In the partial survival regime, neutrinos produced from secondary protons/neutrons can play a dominant role in the high-energy range. Because of the relatively high event rate of engine-driven SNe compared to classical GRBs, the diffuse neutrinos can reach a flux level of $\sim 10^{-11}-10^{-10}\rm~{}GeV~{}cm^{-2}~{}s^{-1}~{}sr^{-1}$. This cumulative flux can be detected by GRAND after $\sim 10$ years observations. We found that the flux of cosmogenic neutrinos is also comparable or even higher than the diffuse source neutrinos. This means that UHECR nuclei may lose most of their energy in the intergalactic space other than inside the source. Finally, we discuss the common origin of IceCube TeV-PeV neutrinos and UHECR nuclei. Our results confirmed Ref. Zhang et al. (2018), which suggested that the hybrid “two-zone” model may give an alternative explanation of UHECR nuclei and IceCube TeV-PeV neutrinos. This picture is consistent with the most popular interpretation that the prompt emission of LL GRBs originates from shock breakout of transrelativistic SNe in a dense stellar wind.

Acknowledgements.

The work of K. M. is supported by Alfred P. Sloan Foundation and the U.S. National Science Foundation (NSF) under grants NSF Grant No. PHY-1620777. B.T.Z. is supported by High-performance Computing Platform of Peking University. The common origin of IceCube neutrinos and UHECRs are discussed in Ref. Zhang et al. (2018). Then, while we are working on this manuscript, related papers came out Boncioli et al. (2018); Samuelsson et al. (2018). One of the differences is that this work focuses on the external RS model, in which prompt emission can be attributed to shock breakout emission of transrelativistic SNe. Our model prediction is also consistent with both of the optical and x-ray data observed for GRB 060218.