Ultrahigh-energy cosmic-ray nuclei and neutrinos from engine-driven supernovae
B. Theodore Zhang Kohta Murase


APS/123-QED






^1Department of Astronomy, School of Physics, Peking University, Beijing 100871, China 
^2Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China
^3Department of Physics, The Pennsylvania State University, University Park, Pennsylvania 16802, USA
^4Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, Pennsylvania 16802, USA
^5Center for Particle and Gravitational Astrophysics, The Pennsylvania State University, University Park, Pennsylvania 16802, USA
^6Yukawa Institute for Theoretical Physics, Kyoto, Kyoto 606-8502 Japan


             
             

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 ∼1% of the core-collapse supernova rate, the observed UHECR spectrum and composition can be explained with the total cosmic-ray energy ℰ_ cr∼10^51 erg. The maximum energy of UHECR nuclei can reach ∼ 10^20-10^21 eV. 
The diffuse flux of source neutrinos is predicted to be ∼ 10^-11 - 10^-10  GeV  cm^-2  s^-1  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.


Valid PACS appear here
                             
Ultrahigh-energy cosmic-ray nuclei and neutrinos from engine-driven supernovae
    Kohta Murase^3,4,5,6
    December 30, 2023
==============================================================================





§ 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 ∼3×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) <cit.>, rarer but faster SNe have been suggested by various authors <cit.> (see Ref. <cit.> for a list of transient sources for UHECRs).
Indeed, observations have revealed that some SNe have a transrelativistic component with shock velocity β_s∼0.1-1. They include engine-driven SNe associated with low-luminosity (LL) gamma-ray bursts (GRBs), such as GRB 980425 with SN 1998bw <cit.>, GRB 060218 with SN 2006aj <cit.>, GRB 100316D with SN 2010bh <cit.>, GRB 161228B with iPTF17cw <cit.>, and GRB 171205A with SN 2017iuk <cit.>. 
Such engine-driven SNe with transrelativistic ejecta have been found through radio observations; examples without GRB counterparts include SN 2009bb <cit.> (but see also Ref. <cit.>) and SN 2012ap <cit.>. 
Another possible example with low-mass subrelativistic ejecta is AT2018cow, although its origin is under debate <cit.>.  
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, ∼100 - 1000  Gpc^-3 yr^-1 <cit.>. It is also comparable to the hypernova rate <cit.>, 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 <cit.>, and the relativistic velocity component can be caused by energy injections via outflows from the engines <cit.>. 

The diffusive shock acceleration theory predicts that the acceleration time scale is ∼ E/(ZeBcβ_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 (Γ≫ 1), but it is known that ultrarelativistic, superluminal shocks are unlikely to be efficient cosmic-ray accelerators <cit.>. 
In this sense, engine-driven SNe, whose shocks are subrelativistic or mildly relativistic, Γ≲ 1-10, would be more favorable for the acceleration of UHECR nuclei. 
The magnetic field strength inferred from radio observations <cit.> also supports that the acceleration of UHECR nuclei is allowed by the Hillas criterion <cit.>.
It has been shown that the composition of UHECR nuclei becomes heavier beyond “ankle” (∼ 4 × 10^18 eV) <cit.> which requires their sources to have super-solar abundance <cit.>. 
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 <cit.>. 

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 <cit.>. 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 <cit.>.
Other source classes including galaxy clusters <cit.>, activate galactic nuclei <cit.>, classical GRBs <cit.>, tidal disruption events (TDEs) <cit.>, and new born pulsars <cit.>.
Neutrinos that directly originate from UHECR nuclei have energy with E_ν≈0.05(E/A)≃ 0.5  EeV (E/10^20  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) <cit.>, Antarctic Ross Ice Shelf Antenna Neutrino Array (ARIANNA) <cit.>, Giant Radio Array for Neutrino Detection (GRAND) <cit.>, Probe Of Extreme Multi-Messenger Astrophysics (POEMMA) <cit.>, and Trinity <cit.>. 
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 <cit.>, although there is no simple picture for the 2017 and 2014-2015 flares <cit.>. LL GRBs and engine-driven SNe will be interesting targets for such follow-up observations, as proposed by Ref. <cit.>. 

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 ∼ 1 - 2 orders of magnitude uncertainty depending on the composition, maximum acceleration energy and source redshift evolution <cit.>. 

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 Γ∼2-10, and the other is the slower transrelativistic ejecta with Γβ∼ 1. 
We revisit the “nucleus-survival problem” that was earlier studied by Ref. <cit.>, 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 <cit.> or transrelativistic ejecta <cit.> 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. <ref>.
In Sec. <ref>, we discuss the acceleration of UHECR nuclei in the RS scenario and compared to the observation results of UHECR nuclei.
In Sec. <ref>, 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. <ref> and give a summary in Sec. <ref>.

Throughout of the paper, we adopt the cgs unit and have notations Q_x ≡ Q /10^x. The cosmological parameters we use are H_0 = 67.3   km  s^-1 Mpc^-1, Ω_m = 0.315, Ω_Λ = 0.685 <cit.>. 



§ 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 <cit.>, becoming transrelativistic SNe. The material collides into the surrounding circumburst medium (CBM), leading to the formation of external reverse-forward shocks <cit.>. 

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_×. If the engine duration T_ eng is longer than the deceleration time t_ dec, corresponding to the thick shell case, the shell crossing time is approximately given by ∼ T_ eng. In more general, it is given by t_×∼ max[T_ eng,t_ 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. <ref>. 

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. <ref>.  Heavy nuclei can be extracted from the inner stellar core <cit.> or synthesized during the expansion of the outflow <cit.>. 
In this work, we adopt the nuclear composition model, Si-R 1, proposed in Ref. <cit.> as a fiducial model of the jet composition, while the transrelativistic ejecta has a composition similar to the hypernova model in Ref. <cit.>. 





 §.§ 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 <cit.>. The RS emission is important to constrain the initial Lorentz factor as well as the baryonic component in the jets <cit.>. 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 <cit.>, 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 <cit.>. 
For our purpose, we only consider the shock physics at radius R_×, where RS finishes crossing the GRB ejecta shell <cit.>. After the shock crossing, the RS light curve is expected to decline rapidly with time ∝ t^-2 as a result of the rapid cooling of electrons <cit.>. 
In the following, our calculations are based on the work in Ref. <cit.>.  

We assume that the GRB ejecta have an initial Lorentz factor Γ_0 and isotropic-equivalent ejecta energy ℰ_ k. The CBM density is given by ϱ_ cbm, and ϱ_ ej is the number density of the GRB ejecta, which can be estimated by ϱ_ ej = ℰ_ k / 4 π m_p c^2 Γ_0 (Γ_0 Δ) R^2 with Δ as the ejecta shell thickness in the engine frame. 
For a homogeneous CBM with ϱ_ cbm =  const, the RS crossing radius R_× is estimated to be R_×≃ 5.6 × 10^16ℰ_ k, 51.5^1/4ϱ_ cbm, 1^-1/4 T_4^1/4 cm, where we adopt the “thick ejecta shell” case considering Δ = c T > R_× / 2 Γ_0^2, and T = 10^4  s is the engine frame duration of the GRB ejecta <cit.>. 
This is justified when the central engine is active for a sufficiently long time.
Note that if R_× / 2 Γ_0^2 > cT, we should consider the “thin ejecta shell” Δ = R_× / 2 Γ_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 Γ_×≃ 6.3 ℰ_ k, 51.5^1/8ϱ_ cbm, 1^-1/8 T_4^-3/8, where we adopt the condition ϱ_ ej / ϱ_ cbm≪ 4 Γ_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,  Γ_×^'≈ (1/2) (Γ_× / Γ_0 + Γ_0 / Γ_×) ≃ 1.1. The magnetic field strength of the shocked GRB ejecta can be estimated assuming a fraction ϵ_B of the post-shock energy density is converted into the magnetic energy, B_×≃ 1.6 ϵ_B, -1.3^1/2ℰ_ k, 51.5^1/8ϱ_ cbm, 1^3/8 T_4^-3/8 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 ν_i = 3 e γ_i^2 B_×Γ_× / 4 π m_e c
with some characteristic Lorentz factor of electrons, γ_i.  Here ν_i represents ν_m (injection frequency), ν_a (self-absorption frequency), and ν_c (cooling frequency), respectively. The injection synchrotron frequency in the engine frame is

    ν_m    ≃    1.4 × 10^13 [(Γ_×^' - 1) / 0.1]^2  
       ×   ϵ_e, -1^2 f_e, -2^-2ϵ_B,-1.3^1/2ℰ_ k, 51.5^1/4ϱ_ cbm, 1^1/4 T_4^-3/4 Hz,

with ϵ_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. <cit.>, and the chosen value s = 2.4 is already used in previous works in order to reproduce the external reverse-forward shock emission <cit.>.
The electron cooling Lorentz factor depends on the ratio between electron radiation time scale and dynamical time scale γ_c = 6π m_e c^2 Γ_× / σ_T (Y+1)R_× B_×^2, where Y is the Compton Y parameter. 
The typical cooling frequency in the slow cooling regime is 

    ν_c ≃ 4.1 × 10^13ϵ_B,-1.3^-3/2ℰ_ k, 51.5^-1/2ϱ_ cbm, 1^-1 T_4^-1/2 Hz,
 
and the self-absorption frequency is

    ν_a    ≃    3.8 × 10^9ϵ_B,-1.3^1/5ϵ_e, -1^-1 f_e, -2^8/5ℰ_ k, 51.5^19/40
       ×   ϱ_ cbm, 1^13/40 T_4^-33/40 [(Γ_×^' - 1) / 0.1]^-1 Hz.

The latter is estimated by setting the self-absorption optical depth τ (ν_a) to unity <cit.>. 

The synchrotron emission from RS can be described as broken power law <cit.> (ν_a < ν_m < ν_c)

    dn/dε   =    n_ε,  max{[              (ε_a / ε_m)^-2/3 (ε / ε_a)                          ε_ min < ε≤ε_a;                          (ε / ε_m)^-2/3                             ε_a < ε≤ε_m;                      (ε / ε_m)^-(s+1)/2                             ε_m < ε≤ε_c; (ε_c / ε_m)^-(s+1)/2 (ε / ε_c)^-(s+2)/2                          ε_c < ε≤ε_ max ].

where n_ε,  max = L_ε,  max/4π R_×^2 c ε_m
is the normalization of the differential photon number density.
The comoving frame luminosity per unit energy is

    L_ε,  max   =   1/2πħf_e N_e √(3) e^3 B_×/m_e c^2
       =    6.9 × 10^55 f_e, -2ϵ_B,-1.3^1/2ℰ_ k, 51.5^9/8ϱ_ cbm, 1^3/8 T_4^-3/8 s^-1,

where N_e = ℰ_ k / Γ_0 m_p c^2. 
We show the comoving frame differential photon number density (blue lines) in Fig. <ref>, which are calculated from following different parameter sets:


  * Jet-A: ℰ_ k = 3 × 10^51 erg, T = 10^4  s, Γ_0 = 10, ϱ_ cbm = 10  cm^-3, ϵ_e = 0.1, f_e = 0.01, ϵ_B = 0.01, and s = 2.4.

  * Jet-B: ℰ_ k = 3 × 10^51 erg, T = 10^4  s, Γ_0 = 10, ϱ_ cbm = 1  cm^-3, ϵ_e = 0.1, f_e = 0.01, ϵ_B = 0.01, and s = 2.4.

  * Jet-C: ℰ_ k = 4 × 10^51 erg, T = 10^4  s, Γ_0 = 5, ϱ_ cbm = 10000  cm^-3, ϵ_e = 0.001, f_e = 0.01, ϵ_B = 0.001, and s = 2.1.



The detailed value of ϵ_e, f_e, and ϵ_B depends on the microphysics of the collisionless shocks, and which are still unclear from first principles <cit.>. 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 <cit.>, as shown in Fig. <ref>. We also note that our model does not overshoot either the optical or x-ray data observed for GRB 060218.















 §.§ 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, ℰ_ ej(>Γ_ ejβ_ ej) ∝ (Γ_ ejβ_ ej)^-α, where β_ ej = v_ ej/c is the normalized velocity and Γ_ ej = 1/ (1 - β_ ej^2)^1/2 is the Lorentz factor. 
The profile is usually steep α∼ 5 for ordinary SNe <cit.>, whereas GRBs usually have a much flatter profile with α≃ 0.4 <cit.>. 
SNe associated with GRBs, including LL GRBs, have intermediate properties with α∼ 2-3. This implies that most of the kinetic energy is taken by low-velocity ejecta of these SNe <cit.>.
For example, let us consider hypernovae with a total kinetic energy of ℰ_ ej∼3×10^52 erg and a typical velocity of β_ ej∼0.1 (corresponding to an ejecta mass of M_ ej∼3 M_ sun). Taking into account weights of the ejecta velocity distribution, the fast ejecta component has a typical velocity of β_ fej∼0.3 with its characteristic energy of ℰ_ fej∼3×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, ϱ_ cbm = A r^-2, where r is the radial distance, A ≡ 3 × 10^35 A_⋆ cm^-1, and A_⋆ is the mass-loss rate to wind speed ratio normalized to 10^-5 M_⊙ yr^-1 / 10^3  km s^-1 <cit.>.
Then, the deceleration radius for the fast velocity ejecta is estimated to be R_ dec≃ 2.1 × 10^18ℰ_ fej, 51.5 A_⋆, 0.7^-1β_ fej, -0.5^-2 cm. 
Here, we give a rough approximation for the magnetic field strength as B_ rs≈ (ϵ_ B 4 πϱ_ cbmβ_ fs^2 c^2)^1/2, where we have used the energy density balance, i.e., ϱ_ fejβ_ rs^2=ϱ_ cbmβ_ 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,

    ν_m    ≃    1.5 × 10^10ℰ_ fej, 51.5^-1ϵ_e, -1^2 f_e, -2^-2ϵ_ B, -0.3^1/2 A_⋆, 0.7^3/2
       ×   β_ rs, -0.5^4 β_ fej, -0.5^3  Hz,

where we assume the accelerated electrons following a power-law distribution with spectral index s = 2.4 and γ_m = ((s-2)/(s-1))(ϵ_e / f_e) (m_p / m_e) β_ rs^2.
Similarly, the electron cooling frequency is

    ν_c    ≃    1.1 × 10^14ℰ_ fej, 51.5ϵ_ B, -0.3^-3/2 A_⋆, 0.7^-5/2
       ×   β_ fej, -0.5^-5 Hz.

and the synchrotron self-absorption frequency is 

    ν_a    ≃    2.5 × 10^7 ℰ_ fej, 51.5^-1ε_e,-1^7/5 f_e, -2^-7/5ϵ_B,-0.3^1/5
       ×    A_⋆, 0.7^9/5β_ rs, -0.5^17/5β_ fej, -0.5^18/5 Hz.

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

    L_ε,  max≃ 1.2 × 10^56 f_e, -2ϵ_ B, -0.3^1/2 A_⋆, 0.7^3/2β_ fej, -0.5 s^-1.

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

TRSN-A: ℰ_ fej = 3 × 10^51 erg, β_ fej = 0.8, A_⋆ = 1, ϵ_e = 0.1, f_e = 0.01, ϵ_B = 0.1 and s = 2.4.

TRSN-B: ℰ_ fej = 3 × 10^51 erg, β_ fej = 0.3, A_⋆ = 5, ϵ_e = 0.1, f_e = 0.01, ϵ_B = 0.5 and s = 2.4.



§ UHECR NUCLEI FROM ENGINE-DRIVEN SNE



 §.§ Acceleration and survival of UHECR nuclei

Charged nuclei can be accelerated to higher energies through the first-order Fermi acceleration mechanism <cit.>.
The maximum acceleration energy is determined by the condition that the acceleration time scale should be smaller than dynamical time scale t_ acc≤ t_ dyn, as well as various energy cooling time scales, t_ acc≤ t_ cool≡ 1 / (t_ ad^-1 + t_Aγ^-1 + t_ syn^-1 + t_ IC^-1).
Here t_ acc is the acceleration time, t_ dyn is the dynamical time, t_ cool is the total energy cooling time, t_ ad∼ t_ dyn is the adiabatic cooling time, t_ syn is the synchrotron cooling time, t_ IC is the inverse Compton cooling time, and t_Aγ is the photohadronic cooling time <cit.>.
In the case of the age-limited acceleration (t_ cool≥ t_ dyn), the engine frame maximum energy of nuclei accelerated in jets can be estimated as,

    E_A,  max^ jet   ≃   η^-1 Z e B_× (R_× / Γ_×) Γ_×
       ≃    7 × 10^20η^-1 (Z/26) ϵ_B,-1.3^1/2
       ×   ℰ_ k, 51.5^1/4ϱ_ cbm, 1^1/4 T_4^1/4 eV,

where η≥ 1 is the acceleration efficiency <cit.>, Z is the nucleus charge number, R_× / Γ_× is the comoving width of the shocked ejecta. While the maximum energy of nuclei accelerated in transrelativistic ejecta is estimated to be,

    E_A,  max^ fej   ≃   η^-1 Z e B_ rs R_ decβ_ fej
       ≃    8.3 × 10^19η^-1 (Z/26) ϵ_ B, -0.3^1/2
       ×    A_⋆, 0.7^1/2β_ fej, -0.5^2  eV.

We see that the acceleration energy of heavy nuclei can exceed ∼ 10^20 eV in both jets and transrelativistic ejecta. Note that the maximum energy is sensitive to the ejecta velocity, E_A,  max^ fej∝β_ 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 <cit.>. Considering the diversity of engine-driven SNe and the uncertainty of the ejecta profile <cit.>, we simply assume that high-velocity ejecta (β_ ej≥ 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 ℰ_ CRacc = 6 × 10^50ξ_ CRacc, -0.7ℰ_ k, 51.5 erg where ℰ_ k = 3 × 10^51 erg is the kinetic energy of the jet or transrelativistic ejecta and ξ_ CRacc = 1/5 is the energy fraction of accelerated CRs. 
Then, we can derive the energy injection rate density, Q_ UHECR = ℰ_ CRescρ̇≃ 2 × 10^44ℰ_ CResc, -1ρ̇_3.5 erg Mpc^-3 yr^-1, where ρ̇ is the local event rate of engine-drivne SNe. We found the energy injection rate ensity Q_ UHECR is comparable to the required energy budget of UHECR nculei <cit.>.

In Fig. <ref>, 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_ esc≤ t_ dyn. We can see that cosmic rays around maximum energy E_A,  max^ (esc)≃ x_ esc E_A,  max^ (acc)≃ 10^18.3 x_ esc, -0.2 eV can escape efficiently, where x_ esc is the ratio of the escape boundary to the width of the shocked region <cit.>. The maximum escape energy adopted in this work can be achieved assuming the acceleration efficiency is η = 8 for the jet and η=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 τ_Aγ≡  t_ dyn / t_Aγ- int can be estimated analytically <cit.>

    τ_Aγ   ≈   2 (R_× / Γ_×) n_0 ε_0/ (1 + α)σ_ GDRΔε̅_ GDR/ε̅_ GDR(E_A/E_A0)^α - 1
       ≃    4.6 ℰ_ k, 51.5^1/4ϱ_ cbm, 1^-1/4 T_4^1/4
       ×   (ε_0 / 10^-2 eV ) (E_A / 10^20 eV )^-1/3,

where we adopt model Jet-A and Oxygen nuclei with σ_ GDR = 3.1 × 10^-27 cm^2, Δε̅_ GDR = 6  MeV, and ε̅_ GDR = 22.35 MeV <cit.>. 
The target photon number density is expressed as dn/dε = n_0 (ε / ε_0)^-α, where α = 2 / 3 is the spectral index. The effective optical depth is defined as f_Aγ≡  t_ dyn / t_Aγ, where t_Aγ is the energy loss time scale of UHECR nuclei. The relation between effective optical depth f_Aγ and optical depth τ_Aγ can be expressed as <cit.> 

    f_Aγ   =   τ_Aγκ̅_Aγ
       ≃    0.9 ℰ_ k, 51.5^1/4ϱ_ cbm, 1^-1/4 T_4^1/4(ε_0 / 10^-2 eV ) 
       ×   (E_A / 10^20 eV )^-1/3( κ̅_Aγ / 0.2 ),

where κ̅_Aγ is the inelasticity which represents average energy losses of the photodisintegration process <cit.>.
We can see that UHECR nuclei are in the partial survival regime for model Jet-A which satisfy the condition τ_Aγ >1 and f_Aγ < 1 <cit.> (see Fig. <ref>). 
The complete survival of UHECR nuclei is possible when τ_Aγ < 1 for model Jet-B (see Fig. <ref>). 

We also show the survival of UHECR nuclei for model TRSN-A (partial survival regime) in Fig. <ref> and TRSN-B (complete survival regime) in Fig. <ref>, 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_⋆ =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.













 §.§ Comparison with Auger and TA data

UHECRs have been observed on Earth by various experiments <cit.>, such as the Pierre Auger Observatory (Auger) <cit.> and Telescope Array (TA) collaboration <cit.>, 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” ∼ 4 × 10^18 eV, with a suppression at ∼ 5 × 10^19 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 <cit.>.

In this section, we aim to fit the energy spectrum and composition of UHECR nuclei detected by Auger <cit.> and TA <cit.>.
For TA data, we decrease the flux of UHECR nuclei by a fraction of ∼ 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_ max of Auger <cit.> and TA <cit.> directly, which may be affected by the particular treatment of biases and detector efficiencies <cit.>, 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 <cit.>.
For simplicity, we consider 1 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. <cit.>. This model is based on the hierarchical structure formation scenario, taking into account galaxy formation and evolution. 
We performed χ^2 analysis as in Ref. <cit.> to show the quality of our fitting. The spectrum energy fitting range is from 10^18.45  eV to 10^21.15  eV, where the uncertainties on the flux are taken from the measurements of Auger <cit.> and TA <cit.>. As in Ref. <cit.>, we take into account the effect of systematic uncertainty on the energy scale, E = E^ Auger (1 + δ_E). 
We derived the first two moments of the X_ max distributions using the EPOS-LHC hadronic interaction model <cit.>, and compared to the measurements from Auger <cit.> and TA <cit.>, respectively.

The results are shown in Fig. <ref>-<ref>, where we consider the complete survive regime, Jet-B and TRSN-B.
The change to a lighter composition in the X_ max distributions is caused by the photodisintegration process of nuclei from nearest sources, as explained in Ref. <cit.>.
We find model Jet-B can fit the Auger data very well, χ^2 /  d.o.f. = 1.5, where the best-fitting parameters is E_p,  max^ esc = 18.3 and δ_E = 0.13, as shown in Fig. <ref>.
On the other hand, model TRSN-B gives the best fit to TA data, χ^2 /  d.o.f. = 0.9, where the best-fitting parameters are E_p,  max^ esc = 18.3 and δ_E = 0.01, see Fig. <ref>. 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 ∼3σ 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 <cit.>. 
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. <cit.>. 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. <ref> and Sec. <ref>.







§ NEUTRINOS FROM ENGINE-DRIVEN SNE

The total cross section of the photomeson production process can be decomposed into three parts: Δ-resonance region at peak energy ε̅_Δ≃ 340  MeV (with a contribution from direct production), higher resonances with two or three peaks in the energy range 500  MeV < ε̅ < 1.5  GeV, multipion production region at much higher energies (> 1.5  GeV) (see, e.g., Ref. <cit.> for details). 
The box approximation is commonly used, in which σ_pγ≈σ_Δδ (ε̅ - ε̅_Δ) Δε̅_Δ. Here σ_Δ∼ 4.4 × 10^28  cm^2, ε̅_Δ∼ 0.34  GeV, and Δε̅_Δ∼ 0.2  GeV <cit.>.
Note that the branching ratio between neutron pion π^0 production channel and charged pion π^± production channel can effectively approach to 1:1:1 owing to multipion production <cit.>. The final neutrino energy spectrum from charged pion and muon decay can be calculated numerically through the code SOPHIA <cit.> or GEANT4 <cit.>. 

For nuclei, we assume that the cross section on nuclei scales linearly with their mass, σ_ Aγ (ε̅) = A σ_ pγ (ε̅) and the inelasticity is assumed to have κ_ Aγ (ε̅) = (1/A)  κ_ pγ (ε̅). 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 <cit.>. 

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) <cit.>.

In the first case, the all-flavor neutrino energy spectrum can be simply estimated using the following formula <cit.>,

    ε_ν^2 dN_ν/dε_ν≈3/8 f_ sup f_ mes (ε_A) (1 - f_Aγ(ε_A))ε_A^2 dN_A/dε_A,

where f_ sup is the suppression factor due to meson and muon cooling, f_ mes(ε_A) is the effective optical depth for the photomeson production, and (1- f_Aγ) represents the fraction of survival nuclei. 
The factor 3/8 is due to the fact that half of the energy ∼ 1/2 goes into charged pions and the pion decay convert ∼ 3/4 of their energy into neutrinos. 
Assuming each nucleon has energy ∼ 1/A of the parent nucleus with ε_A, then the effective optical depth to the photomeson production on nuclei can be approximated by f_ mes (ε_A) ≈ f_pγ (ε_A / A). The typical neutrino energy is ε_ν≈ 0.05 ε_A / A.
While in the second case, the all-flavor neutrino energy spectrum can be calculated as <cit.>

    ε_ν^2 dN_ν/dε_ν≈3/8f_ sup f_pγ(ε_p) f_Aγ(ε_A) ε_A^2 dN_A/dε_A,

where ε_p ≈ε_A / A and f_Aγ represents the fraction of photodisintegrated nuclei. In this case, the typical neutrino energy depends on the secondary proton energy ε_ν≈ 0.05 ε_p. Note that the secondary protons/neutrons take a fraction ∼ f_Aγ (ε_A) of the parent nuclei total energy.
The two cases can be combined into one formula

    ε_ν^2 dN_ν/dε_ν≈3/8 f_ sup f_pγ (ε_A / A) ε_A^2 dN_A/dε_A.

under the approximation f_ mes (ε_A) ≃ f_pγ (ε_A / A). 
The suppression factor accounts for the energy cooling on the primary pions and secondary muons, f_ sup≃ f_ sup, π f_ sup, μ <cit.>. The pion cooling suppression factor is calculated as f_ sup, π = 1 -  exp (-t_π,  cool / t_π,  decay), where t_π,  decay is the decay time scale of pions. The mean life time for a pion to decay is τ_π = 2.6 × 10^-8 s in the pion rest frame, then the pion decay time scale in the shock comoving frame is t_π,  decay = (ε_π / m_π) τ_π.  The dominant cooling process of pions are synchrotron cooling, t_π,  cool^-1∼ t_π,  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 <cit.>, where the photomeson production process is calculated using the SOPHIA package <cit.>. In CRPropa,  the mean free path of nuclei have the following scaling relations <cit.>,

    λ_A, Z^-1 (ε_A) ≃ A λ_p/n^-1(ε_A/A),

where λ_p, λ_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 ε_A-1 = (A-1) / A ε_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 <cit.>.



 §.§ Neutrino fluences

The observed neutrino flucences can be estimated using the following formula <cit.>,

    E_ν^2 ϕ_ν   =   1+z/4π d_L^2Γε^2_νdN_ν/dε_ν,

where d_L is the source luminosity distance, Γ is the source Lorentz factor and ε^2_ν dN_ν / dε_ν is the comoving frame neutrino energy spectrum inside the source. 

In Fig. <ref>, 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 ℰ_ CRacc = 6 × 10^50ξ_ CRacc, -0.7ℰ_ k/fej, 51.5 erg, as shown in Fig. <ref>.
 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_π,  cool = t_π,  decay, ε_ν, π = 2.7 × 10^18 B_1^-1 eV and ε_ν, μ = 1.5 × 10^17 B_1^-1 eV, where t_π,  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 <cit.>.

Similar to Fig. <ref>, we show the all-flavor neutrino fluences estimated assuming complete survival regime corresponding to model Jet-B (or TRSN-B) in Fig. <ref>. 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 ∼ 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. <ref> and Fig. <ref>.
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 <cit.>. (The similar trend was found by Ref. <cit.>.)
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. <ref>, 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 <cit.>.
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. 





 §.§ 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,

    E_ν^2 Φ_ν   =   1/4π∫_z_ min^z_ max dz  dV_c/dzF(z) ρ̇_̇0̇/1 + z E_ν^2 ϕ_ν,

where E_ν^2 ϕ_ν is the neutrino fluence derived in the previous section. The event rate of engine-driven SNe in the local universe is ρ̇_̇0̇ and F(z) is the redshift distribution parameter estimated from long GRBs <cit.>. The redshift dependence of the comoving volume is

    dV_c/dz = c/H_04 π d_L^2/(1+z)^2√(Ω_M (1+z)^3 + Ω_Λ),

where d_L is the luminosity distance to the source.

We show the energy spectrum of diffuse source neutrinos in Fig. <ref>.  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 (ξ_z ≃ 3 in Ref. <cit.>) and photodisintegration bound (olive line) for pure Silicon composition <cit.>. Note that in the non-evolution case, the diffuse neutrino fluxes can be ∼ 5 times lower than in the fast-evolution case.





 §.§ Cosmogenic neutrinos

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

    Φ_ν(E_ν)    ≡   dN_ν/dE_ν dA dΩ dt
       =   ∑_A^'c/4π H_0∫_z_ min^z_ maxdz F(z) ρ̇_0/(1 + z) √(Ω_M (1+z)^3 + Ω_Λ)
       ×   ∫_E^'_ min^E^'_ max dE^'dN_A^'/dE^'dη_A^'ν(E^', E_ν, z)/dE_ν,

where η_A^'ν is the neutrino yield function which reflects the fraction of neutrinos with energy E_ν originated from nuclei with mass number A^' and energy E^' at redshift z.
The value of η_A^'ν can be estimated numerically using CRPropa 3 <cit.>, where UHECR nuclei are propagated through the universe as in Sec. <ref>.
We show the results in Fig. <ref> where the flux of cosmogenic neutrinos (green line) can reach a level of a few × 10^-10 GeV cm^-2 s^-1 sr^-1 and have peak energy around 0.1 - 1  EeV. Note the small bump appeared in the lower energy range ∼ 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 ∼ 10 yr observations <cit.>.
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) <cit.> and tidal disruption events (TDEs) where a white dwarf is disrupted by an intermediate mass black hole (blue dotted line) <cit.>. The flux of cosmogenic neutrinos is sensitive to the source redshift evolution because most of the detected neutrinos are produced at redshift z ∼ 1, while the observed UHECR nuclei mainly originated from sources within the local universe ≤ 100 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 <cit.>.

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




§  DISCUSSION




 §.§ Common origin of IceCube neutrinos and UHECRs?

The flux of the observed diffuse neutrinos within energy range from ∼ 100 TeV to a few PeV by IceCube <cit.> is comparable to the Waxman-Bahcall (WB) bound for a spectral index of s=2.0 <cit.> and the nucleus-survival bound for a spectral index of s∼2.2-2.3 <cit.>. The fact that the diffuse energy fluxes of all three messengers are similar has led to the development of the unification picture <cit.>.

The flux of diffuse source neutrinos predicted in this work is ∼ 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 <cit.>, 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 <cit.>. 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 <cit.>. Transrelativistic SNe naturally originate from choked jets <cit.>, and high-energy neutrino emission from low-power choked jets can explain the IceCube data even in the 10-100 TeV range <cit.>. 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. <cit.>, the successful jet scenario <cit.> 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. <cit.> attempted to simultaneously explain the PeV neutrino flux and UHECR data. 





 §.§ Maximum acceleration energy?

We showed that the maximum energy of UHECR nuclei reaches ∼ a few×(Z / 26) 10^20 eV in the framework of RS formed by both jets and transrelativistic ejecta. 
We confirmed the findings of the original papers <cit.>, which found that the maximum energy can be as high as ∼ a few×10^20 eV for LL GRBs like GRB 06218 without violating observational data (see also Refs. <cit.>).   
A recent paper by Ref. <cit.> claimed that it may be difficult for protons or iron nuclei to reach the maximum acceleration energy ∼ 10^20 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_γ∼10^47 erg s^-1 and Γ∼1-10, the maximum energy can exceed ∼ 10^20 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. 






§ SUMMARY

In this work, we investigated UHECR nuclei and neutrinos originating from engine-driven SNe. 
In Sec. <ref>, 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 1 in Ref. <cit.> 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. <ref>, 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 ∼ 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 ∼ 10^-11 - 10^-10 GeV cm^-2 s^-1 sr^-1. This cumulative flux can be detected by GRAND after ∼ 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. <cit.>, 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. 



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. <cit.>. Then, while we are working on this manuscript, related papers came out <cit.>.  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.








































