Estimating the Neutrino Flux from Choked Gamma-Ray Bursts
Michela Fasano, Silvia Celli, Dafne Guetta, Antonio Capone, Angela Zegarelli, Irene Di Palma
PPrepared for submission to JCAP
Estimating the Neutrino Flux fromChoked Gamma-Ray Bursts
Michela Fasano, Silvia Celli, , Dafne Guetta, Antonio Capone , Angela Zegarelli , and Irene Di Palma , Dipartimento di Fisica dell’Universit`a La Sapienza, P. le Aldo Moro 2, I-00185 Rome, Italy Istituto Nazionale di Fisica Nucleare, Sezione di Roma, P. le Aldo Moro 2, I-00185 Rome,Italy ORT-Braude College, Carmiel, IsraelE-mail: [email protected], [email protected],[email protected], [email protected],[email protected], [email protected]
Abstract.
The strong constraints from the Fermi-LAT data on the isotropic gamma-ray back-ground suggest that the neutrinos observed by IceCube might possibly come from sourcesthat are hidden to gamma-ray observations. A possibility recently discussed in the litera-ture is that neutrinos may come from jets of collapsing massive stars which fail to breakout of the stellar envelope, and for this reason they are also known as choked jets . In thispaper, we estimate the neutrino flux and spectrum expected from these sources. We performdetailed calculations of pγ interactions, accounting for all the neutrino production channelsand scattering angles. We provide predictions of expected event rates for operating neu-trino telescopes, as ANTARES and IceCube, as well as for future generation telescopes asKM3NeT. We find that if 10% − − . − yr − . This rate is consistent with the rate of GRBs detected ingamma rays. Keywords: neutrinos – stars: jets – stars: massive – supergiants – supernovae: general a r X i v : . [ a s t r o - ph . H E ] J a n Introduction
The diffuse flux of high-energy neutrinos observed by the IceCube experiment [1] hasn’t shownyet the presence of a significant event clustering [2]. In other words, the angular distributionof the events collected so far does not highlight any specific location of the sky to search forcorrelated counterparts, e.g. in the electromagnetic domain or in the cosmic-ray anisotropy.This might indicate that the majority of neutrino sources is of extra-galactic origin.The possible production mechanisms of neutrinos are respectively proton-proton (pp)collisions and photo-meson (p γ ) interactions. In both cases, neutrinos would emerge fromboth the direct decay of charged pions as well as from muon decay, while gamma rays resultfrom the neutral pion sudden decay. As such, within hadronic sources, a similar amount ofhigh-energy electromagnetic radiation accompanies the production of neutrinos.However, in the presence of intense radiation fields intrinsic to the source, the emit-ted gamma rays are likely to be absorbed in the pair production process (see e.g. [3]). Inaddition, during their propagation towards Earth, multi-TeV gamma rays can initiate elec-tromagnetic cascades on both the extra-galactic background light (EBL) and the cosmicmicrowave background (CMB), and be reprocessed towards lower energies. By comparingthe isotropic diffuse gamma-ray background (IGRB), observed by Fermi [4], with the as-trophysical neutrino flux, [5] first noted a strong tension among the two, in the hypothesisthat the pp interaction in transparent gamma-ray sources is originating both fluxes. Thisconstrain motivates the alternative study of p γ sources, opaque to high-energy gamma rays,as main neutrino sources.Choked GRBs belong to such a population [6]. These GRBs are believed to be origi-nated from core collapse of massive stars resulting in a relativistic fireball jet, which is unableto break through the stellar envelope. A choked jet deposits all its energy in a cocoon, thateventually breaks out from the star releasing energetic material at sub-relativistic velocities.Such fast moving material has a unique signature that can be detected in early time super-nova spectra (see [7–9]). Neutrino emissions from choked GRBs have hence attracted muchattention [10–13], because of the possibility of explaining the IceCube diffuse flux withoutincurring into inconsistencies with the IGRB.Estimates of neutrino fluxes from astrophysical sources following the pγ interaction havebeen performed in the past both with analytical [14, 15] and numerical calculations [16, 17],but no simulation for choked GRBs sources has been attempted so far. To study the expectedneutrino flux on Earth from such ν bright (while γ dark) sources, we have developed a MonteCarlo (MC) code to simulate the particles kinematics and interactions inside the jet, allowingfor a detailed description of the micro-physics involved in the photo-meson process. Oncethe particle spectra at the source are obtained, our code calculates propagation throughcosmic distances from the source to Earth, thus allowing a reliable estimate of the numberof neutrino events expected in current and future neutrino telescopes.The paper is organized as follows: Section 2 briefly describes the jet dynamics insidethe source in order to characterize the interacting particles spectra. Section 3 presents anoutline of the main structure of our simulation for the p γ interaction. Calculations leadingto the estimate of neutrino events in telescopes from a single source are contained in Section4. Section 5 describes the prediction for the diffuse neutrino spectrum coming from chokedsources. Moreover, Appendix A illustrates the detailed kinematic calculations adopted in thecode. – 2 – Choked jet and primary particles spectra
At the end of their lives, massive stars collapse into a compact object, and it is widely believedthat aforesaid phenomena can produce energetic jets [18, 19]. For our progenitor source, weconsider a red supergiant star possessing a thick stellar envelope (a detailed description ofthis source can be found in [6]): the star is surrounded by its helium and hydrogen envelopes,which are assumed to have radii equal to, respectively, r He ∼ cm and r H (cid:38) cm. Thecore collapse of a red supergiant star enclosed by aforementioned envelopes will result in atype II supernova, having a central iron core of radius r Fe ∼ cm. The source densityin the hydrogen envelope, ρ H ∼ − g / cm , allows for particle acceleration to relativisticvelocities, namely the photomeson interaction arises beyond the helium radius.We consider a jet propagating inside this source and define t jet as the jet lifetime. If thejet crossing time is longer than t jet , then the jet will result choked inside the stellar envelope,namely it will not be able to break through the star surface [18]. In other words, when thecentral engine powers the jet for a shorter time than it is required to travel across the starradius, the jet is stalled inside the source.Within our simulation we adopt the following benchmark values for the jet characterizingfeatures: the Lorentz factor Γ = 100, the lifetime t jet = 100 s, and the isotropic luminosity L iso = 10 erg / s. Note that the isotropic energy of the jet is given by E iso = t jet × L iso =10 erg: this quantity, inferred from gamma-ray observations, reflects the energy of electronsaccelerated to high energies by the internal shocks. These shocks are expected to accelerateprotons even more efficiently than electrons (see e.g. [20]), therefore we assume that theaccelerated proton energy is at least comparable to that carried by accelerated electrons[21, 22]. Moreover, the energy in magnetic fields present inside the source is expected to bea fraction (cid:15) b = 0 . E iso .Protons and electrons are accelerated to high energies by the diffusive shock mecha-nism in the internal shock region ( IS ). Electrons are expected to lose all their energy intosynchrotron radiation. However, because of the large optical thickness due to Thomson scat-tering, these photons thermalize in the reverse shock region ( RS ) to a temperature of [23]( k being the Boltzmann constant): kT γ (cid:39)
741 eV L / , (cid:15) / , − t − / , ρ / , − (2.1)where (cid:15) e is the fraction of E iso deposited into electrons. Here (and elsewhere) A n denotes A/ n , e.g. L / , = ( L iso / ) / . The subscripts IS , RS and IR refer to quantities measured,respectively, in the internal shock, in the reverse shock and in between the two shock regions.Following Wien’s displacement law, the peak energy of the photon spectrum in the RS isgiven by: E max γ, RS (cid:39) . kT γ ∼ . (2.2)The thermal photons will then escape in the IS, since it is optically thin, where their peakenergy becomes: E max γ, IS = Γ IR E max γ, RS ∼
210 keV . (2.3)The energy distribution of protons produced at internal shocks due to the magneticinhomogeneities of the jet follows a power law, typical of a population of particles acceleratedby the diffusive shock mechanism: dN p dE p = k MC E − p (2.4)– 3 –here k MC is a normalization factor which will be discussed in detail in Section 4.Thermal photons are emitted by accelerated electrons in the reverse shock region, whichexpands towards the central engine. Inside the RS, the photon spectrum is well representedby a Maxwellian with a peak energy given by Eq. (2.2), hence the photon distribution readsas: dN γ dE γ = 2 (cid:114) E γ π (cid:18) kT γ (cid:19) e − EγkTγ . (2.5)When thermal photons escape into the IS, where the photo-meson interaction occurs,their peak energy is boosted to be ∼
210 keV, as obtained in Eq. (2.3) for the benchmarkset of parameters here considered.
The production of charged and neutral mesons inside a jet can take place through the resonantproduction of a ∆ + by highly relativistic protons that interact with photons of the ambientradiation field, namely the thermal photons from the jet head. The code simulates theparticle dynamics in the center of mass rest frame, and then it applies a Lorentz boost tomove to the IS frame through the center of mass Lorentz factor γ ∗ . All the center of massquantities are indicated with an asterisk ∗ . Ultimately, the observed energies at the sourceare obtained by boosting the particles by the jet Lorentz factor Γ. In the center of massframe we use natural units where (cid:126) = c = 1. Moreover, we indicate with β a the speed of thegeneric particle a ( γ a = 1 / (cid:112) − β a ), and with E a , p a , m a its energy, momentum and mass.For detailed kinematic calculations see Appendix A.Once a proton and a photon are extracted, in order to proceed with the interactiona condition must be met, i.e. the production of the ∆ + resonance, given by the thresholdrequirement. In the IS, the condition for a successful interaction can be written as:( m p + m π ) = ( E IS p + E IS γ ) − ( (cid:126)p IS p + (cid:126)p IS γ ) m p + m π + 2 m p m π = m p + 2 E IS p E IS γ (1 − cos θ ) E IS p ≥ m π + 2 m p m π E IS γ (1 − cos θ ) (3.1)where θ indicates the angle between protons and photons in the IS. When the thresholdrequirement in Eq. (3.1) is not satisfied, both the extracted particles are discarded. However,as a consequence of the extremely high number density of thermal photons, all generatedprotons will eventually interact: n γ, IS ∼ . × cm − Γ L − / , t − / , (cid:15) / , − ρ / , − . (3.2)In order to compute the proton path inside the thermal photon gas, it is necessary to invertthe path probability distribution ( λ = 1 / ( σn γ ), σ being the pγ cross section):dP(x p ) = 1 λ e − x p /λ dx p . (3.3)The proton path obtained through the inversion of Eq. (3.3) peaks at ∼ cm for σ = 550 µ b(the resonance peak), while for the out of resonance region value of σ , namely ∼ µ b, the– 4 –roton path peaks at ∼ cm. This means that a proton will interact, on average, aftertravelling 10 . cm inside the source, which has an external helium envelope radius of 10 cm.Thus, we can safely assume that all protons interact inside the source.In addition, inside the IS region, the following relation is satisfied:cos θ = − p IS L (cid:113) ( p IS L ) + ( p IS T ) (3.4) p IS L and p IS T are, respectively, the longitudinal and transverse momenta of the photons in theIS. Defining ξ as the angle between protons and photons in the RS, we obtain the followingalgebraic equations for the photon momenta: − p IS L = − Γ IR E γ, RS (cos ξ + β IR ) p IS T = p RS T = (cid:113) E γ, RS (1 − cos ξ ) . For all simulated interactions, the cosine value obtained from Eq. (3.4) is ∼
1, meaning thatthe collisions are realized head-on.Once the interaction occurs, the Monte-Carlo algorithm takes into account the variousinteraction channels. Within our simulation we considered three possible channels for theinteraction, depending on the photon energy in the proton’s rest frame, (cid:15) r : s IS = 2 E IS p E IS γ (1 − cos θ ) ∼ E IS p E IS γ = m p + 2 m p (cid:15) r s IS being the square of the invariant energy of the system in the IS frame, such that (cid:15) r = 4 E IS p E IS γ − m p m p . (3.5)The p γ cross section as a function of (cid:15) r is characterized by a very peaked resonance regionaround the threshold, mainly due to the two reaction channels we are studying (namely theones giving the p + π and n + π + final states): in our simulation, the interaction probability P int in this region is set equal to 1. Above the threshold, the cross section is dominatedby multiple pion final states, and it becomes approximately constant for (cid:15) r >
10 GeV; inthese energy regions, multiple values were used for P int within the code, depending on thecross-section. For a more detailed characterization of the p γ cross section see e.g. [24].The boundaries between each interaction channel in our simulation are given by thefollowing numerical values: • . ≤ (cid:15) r < . P int = 1 and the ∆ + resonance is produced at rest: p + γ → ∆ + → (cid:40) n + π + p + π ; (3.6) • . ≤ (cid:15) r < + resonance is produced along with an additional pion ( N is anucleon, while both π A and π B can be neutral, positively or negatively charged pions): p + γ → π A +∆ + → N + π B (3.7)– 5 –ere, the interaction occurs with two possible values of probability: P int = 0 . (cid:15) r ≤ • ≤ (cid:15) r ≤
100 GeV: this energy window corresponds to the plateau in the crosssection, in which the three types of pions are simultaneously produced. This range isoften referred to as multipion channel , where “ . . . ” indicates other particles created inthe interaction: p + γ → π + π + + π − + . . . (3.8)In this region P int = 0 . (cid:15) r occupies. SeeAppendix A for a more detailed discussion on branching ratios and kinematics in the threechannels.The pion decay products include leptons and photons: π + → µ + + ν µ → e + + ν e + ¯ ν µ + ν µ (3.9) π − → µ − + ¯ ν µ → e − + ¯ ν e + ν µ + ¯ ν µ (3.10) π → γ + γ (3.11)The muon decay was taken into account, and characterized through the Michel pa-rameters , which describe the phase space distribution of leptonic decays of charged leptons.However, the synchrotron cooling of charged pions and muons will suppress the flux of neu-trinos when the lifetimes of these particles become longer than their synchrotron coolingtimescales. In other words, these particles will lose their energy in radiation before decayinginto neutrinos. Noting that the muon’s lifetime at rest ( τ µ ∼ . × − s) is much greaterthan the charged pion’s lifetime at rest ( τ π ± ∼ . × − s), we consider that the synchrotroncooling would only affect muon decay products. This is reasonable within the benchmarkvalues here assumed for the model parameters, as we prove through the following algebraiccalculations. The charged pion lifetime can be written as a function of the parent proton’senergy in the IS as: τ π ± (cid:39) .
038 s E
ISp , , (3.12)where the proton energy is here written in units of eV, and the approximation that theproduced pion energy is E π ± (cid:39) . E p was used. The charged pion synchrotron coolingtimescale instead reads as [11]: τ sync π ± = 4 .
11 s E IS , − , (cid:15) − , − Γ L − / , t jet , ρ − / , − . (3.13)Similarly for muons (and antimuons) we used the approximation that the produced muonenergy is E µ (cid:39) . E p , thus obtaining a lifetime of τ µ (cid:39) .
15 s E
ISp , , (3.14)while its synchrotron cooling timescale instead is [11] τ sync µ = 1 .
87 s E IS , − , (cid:15) − , − Γ L − / , t jet , ρ − / , − . (3.15)We can extract the characteristic value of proton’s energy corresponding to the one at whichthe muon cooling process occurs: τ sync µ ≤ τ µ −→ E IS p ≥ . × GeV . (3.16)– 6 –uons emerging from parent protons with energies exceeding Eq. (3.16) in the IS frame willdecay with a lower probability, and consequently the neutrino production will be reduced.It is hence clear that muons will undergo synchrotron cooling in such a way to affect theneutrino spectra at energies lower than pions . The corresponding cutoff energy for themuons in the IS can be written as E IS µ, cut (cid:39) . E IS p = 1 . × GeV , (3.17)which in the co-moving frame corresponds to E µ, cut = Γ E IS µ, cut = 1 . × GeV . (3.18)Therefore, noting that every lepton produced in the muon’s decay inherits ∼ / E ν, cut ∼ E µ, cut ∼ . × GeV (3.19)Therefore, the muon spectrum in our simulation is affected by an exponential cutoff, with E µ, cut given by Eq. (3.18): when the muon energy exceeds this value, the muon will havea lower probability of decaying. This cutoff is applied to muons in all interaction channels,and the corresponding neutrino spectra obtained through the simulation all depend on thisexponential suppression.The particle spectra at the source obtained through our MonteCarlo simulation areshown in Figure 1. In the top left panel we present the energy spectrum of the interactingprotons, as well as the secondary particles emerging from the neutral pion decay of theinteraction. In this study we focus on neutrino production, hence we do not propagate theresulting gamma rays further in the jet, namely we neglet pair production processes of thehigh-energy photons with the thermal radiation field of the jet. We defer discussion aboutabsorption of gamma rays in the chocked jet to a future work. In the top right panel, we showthe secondary particles from the π + channel, while in the bottom panel the π − . Note thatthe latter is extremely similar to the π + channel, while it contains far less events due to thefact that this channel is not present at the resonance peak, where most of the particles areproduced. These spectra are not normalized yet: as a consequence, the y axis is in arbitraryunits. Transient phenomena such as GRBs are very interesting for a potential neutrino detection,since these analyses are almost background free given the well-defined space-time windowprovided by multi-messenger constraints. To compute the neutrino flux expected on Earthfrom a source at a luminosity distance D ( z ), we scale the MC-obtained fluxes for threefactors: the energy channeled into protons at the source, the dilution of the neutrino flux toEarth, and neutrino oscillations. In fact, the proton’s energy corresponding to synchrotron cooling of charged pions would sit at substan-tially higher energies, namely at τ sync π ± ≤ τ π ± −→ E IS p ≥ . × GeV . – 7 – igure 1 . Particle spectra at the source arising from all interaction channels, as seen in the co-movingreference frame. Top left : π channel; Top right : π + channel; Bottom : π − channel. The first factor comes from energetic considerations for the source: as mentioned inSection 2, we assumed for a generic choked GRB an isotropic energy release of E iso = 10 erg,which we set equal to the accelerated protons energy [6]. So, for the accelerated protonsspectrum at the source in the co-moving frame we have: k MC (cid:90) E max E min E p dN p dE p dE p = E iso (4.1)where E min = 100 GeV and E max = 10 GeV. Solving this integral leads us to an estimateof the scaling factor k MC , which gives the exact fluence of neutrinos at the source in theco-moving frame, normalized as to take into account the amount of energy channeled intoaccelerated protons.Moreover, particles travel through cosmological distances before arriving on Earth, suchthat their flux undergoes a dilution during the propagation. This factor can be calculated asΩ / πD ( z ) , where Ω = 2 π (1 − cos α ) is the solid angle covered by a jet with opening angleof α = 0 . . ◦ . For the choked GRB in our model, we assume a redshift z of ∼ D ( z = 1) = 2 . × m.Neutrino telescopes are sensitive to muon neutrino and muon antineutrino fluxes atEarth, that differ from the fluence at the source due to the effect of neutrino oscillations.Since GRBs are extra-galactic sources, we should mediate neutrino oscillations over cosmicdistances: starting from a flavor ratio of ν e : ν µ : ν τ = 1 : 2 : 0 at the source (as the resultof the pγ interaction), one should expect a flavor ratio of 1 : 1 : 1 on Earth. The oscillation– 8 – igure 2 . All-flavor neutrino fluence on Earth. Left : neutrinos from the resonance peak (greendashed line), from the region above the peak (orange dot-dashed line) and their sum (red line) arecompared in this plot.
Right : neutrinos from choked jets with different Lorentz factors. probabilities P ff (cid:48) of neutrinos changing flavor ( f → f (cid:48) ) are given by [25]: P ff (cid:48) = (cid:88) j (cid:12)(cid:12) U fj (cid:12)(cid:12) (cid:12)(cid:12) U f (cid:48) j (cid:12)(cid:12) , (4.2)so that the oscillated flux for the lepton of flavor f reads as: F ν f ( E ν f ) = P ff F source ν f ( E ν f ) + P f (cid:48) f F source ν f (cid:48) ( E ν f (cid:48) ) , (4.3)where F source ν ( E ν ) = E ν dN ν /dE ν is the generic neutrino spectrum at the source in theco-moving frame, obtained through the MC simulation. This formula is symmetric in theexchange f ↔ f (cid:48) , and it is valid for neutrinos and antineutrinos.Defining F ν ( E ν ) = F ν e ( E ν e )+ F ν µ ( E ν µ )+ F ν τ ( E ν τ ) as the total neutrino energy spectrum,the all-flavor neutrino fluence on Earth can be written as:Φ Earth ν = k MC Ω4 π D ( z ) F ν ( E ν ) (4.4)The all-flavor neutrino fluence on Earth can be seen in Figure 2, where we also accountedfor the redshift dependence of observed energy. In the left panel, neutrinos coming from thedelta resonance peak region (0 . ≤ (cid:15) r < . . + resonance peak for neutrino production leaves outa non negligible number of secondary neutrinos coming from higher interaction channels.Furthermore, the right panel of Figure 2 compares the all-flavor neutrino spectra comingfrom choked jets with different Lorentz factors: the higher the Γ, the higher the neutrinoenergy peak. Conversely, the neutrino fluence φ ν (see Eq. (4.4)) is independent of Γ, asabsolute normalization E iso was set in the co-moving energy frame.We are interested in estimating the number of muon neutrino events expected from anindividual source in different neutrino telescopes, since these constitute the best channel ofevents for astronomical studies. Hence, we proceed by oscillating neutrinos towards Earth,as to obtain the flux of single neutrino flavors. In the case of muon neutrinos, their oscillated– 9 – igure 3 . Muon neutrino and antineutrino flux expected on Earth from a single source at z=1. flux can be calculated through the relation: F ν µ ( E ν µ ) = P µµ F source ν µ ( E ν µ ) + P eµ F source ν e ( E ν e ) (4.5)The values for P µµ and P eµ to be used in Eq. (4.5) are estimated using the latest results forthe mixing angles (see e.g. [26]).The muon neutrino fluence observed on Earth obtained within our simulation is shownin Figure 3, and analytically it reads as:Φ Earth ν µ = k MC Ω4 π D ( z ) F ν µ ( E ν µ ) (4.6)where F ν µ ( E ν µ ) is the oscillated muon neutrino energy spectrum given by Eq. (4.5). Note thatthis result is obtained for the following set of fundamental parameters: z = 1, E iso = 10 erg,Γ = 100, and assuming that 100% of the GRB isotropic energy release goes into acceleratedprotons.Using the muon neutrino fluence thus obtained, an estimate of the number of cosmicneutrino events expected in neutrino telescopes can be performed. Convolving the detec-tor effective area for neutrinos, A ν eff ( E ν , δ ), coming from declination δ , with the differentialnumber of neutrinos per surface element expected on Earth from the source, the number ofevents during the GRB event ( t jet = 100 s) is obtained through the following relation: N events ( δ ) = (cid:90) A ν eff ( E ν , δ ) Φ Earth ν ( E ν ) dE ν (4.7)This evaluation has been performed for ANTARES [27], KM3NeT [28] and IceCube [29].Note that the effective area adopted in this work for KM3NeT refers to the trigger level,being the only available information at the time of writing. In turn, the effective areasadopted for ANTARES and IceCube refer to the analysis level, as these result from detailedstudies performed in the search for neutrino sources by the same Collaboration. The expectednumber of events per energy bin is shown as a function of the energy in Figure 4, while thetotal number of events expected in each neutrino telescope during the choked GRB event– 10 – igure 4 . Number of expected events in ANTARES [27], KM3NeT [28] and IceCube [29] from acharacteristic choked GRB for different declinations of the source. is given in Table 1. The calculations indicate that multiple neutrino events could emergefrom GRBs that can efficiently convert their kinetic energy into primary protons. In fact,we assumed here that the entirety of E iso goes into protons, see Eq. (4.1). Note that thenumber of expected events scales linearly with the fraction of isotropic energy channeledinto non-thermal particles. However, these neutrino events would not be accompanied byhigh-energy electromagnetic radiation because of the intense radiation field of the jet, whichwould absorbe gamma rays. Detailed calculations on the reprocessing of such radiation andpossible signatures will be object of a future study.Detector δ N events ANTARES 0 ◦ < δ < ◦ − ◦ < δ < ◦ − ◦ < δ < − ◦ δ ◦ < δ < ◦ ◦ < δ < ◦ Table 1 . Expected number of muon neutrino events in several detectors from a typical choked GRB( L iso = 10 erg/s, t jet = 100 s, and Γ = 100) located at different declination bands. Effective areasfor each detector were taken from [27–29]. To estimate the diffuse neutrino spectrum, we assume that the rate of choked jets R ( z ) atredshift z follows the star formation rate even if there are evidences that the evolution ofnormal GRBs may deviate from the star formation rate ρ ( z ) [30, 31], i.e., R ( z ) = R ρ ( z ) , (5.1)– 11 –here the star formation rate is [32] ρ ( z ) = (1 + z ) . z ) / . . , (5.2)and R is the local rate of choked GRBs in Gpc − yr − . In this paper we do not consider aluminosity function for choked GRBs but we assume a constant luminosity. A more detailedstudy on the redshift distribution and luminosity function effects on the diffuse neutrino fluxis out of the aim of the present paper but it will be a subject for future works.The diffuse neutrino flux is calculated via integrating the neutrino spectrum over theredshift from 0 to 8, i.e. [33]: E obs ν µ dN ν µ dE obs ν µ ( E obs ν µ ) = c πH (cid:90) E ν µ dN ν µ dE ν µ ((1 + z ) E ν µ ) Ω4 π R ρ ( z ) dz (1 + z ) (cid:112) Ω Λ + Ω M (1 + z ) (5.3)where the cosmological parameters are adopted as H = 70 km s − Mpc − , Ω M = 0 .
3, andΩ Λ = 0 . dN ν µ ( E ν µ ) /dE ν µ is the differential neutrino spectrum at the source, and Ω is thesolid angle of the jet.In order to obtain the best value for the rate of choked GRBs, we compared the pre-dicted flux of neutrino on Earth with the diffuse neutrino flux measured by IceCube. Asexperimental data we used the results obtained by IceCube both with the analysis of HESEevents [34] and the parametrization of the astrophysical diffuse muon neutrino flux [35]. Theflux predicted with the simulation was evaluated for several values of R : the best agreementbetween data and simulation has been obtained for R = (0 . ± .
1) Gpc − yr − (with areduced χ = 1 . R . We note that the best fit value obtained is mostly driven by data on muon neutrinos,as these are characterized by smaller uncertainty bands. The value of R is in agreementwith the Waxmann & Bahcall upper bound [36], also shown in Fig. 5. Recent stacking analyses performed by the current neutrino telescopes lead to the conclu-sion that classical GRBs are not the main sources of the diffuse astrophysical neutrino fluxdetected by IceCube [37, 38]. In addition, also targeted searches for neutrinos in spatial andtiming coincidence with the prompt emission of those GRBs reported evidence for a lack ofcorrelation of neutrino data [39]. Such results have increased even more the interest in analternative GRB class opaque to gamma-rays, choked GRBs, which would also be consistentwith the IGRB observed by Fermi [4].We have estimated the neutrino fluxes from individual choked GRBs and the numberof expected events for present and future neutrino telescopes. This result was obtained forthe following set of benchmark parameters: z = 1, E iso = 10 erg, Γ = 100, t jet = 100 s andassuming that 100% of the GRB isotropic energy release goes into accelerated protons. Wefind that, for the declination in which each detector performs better, the number of muonneutrino events occurring from an individual choked source is equal to 0.03 for ANTARES,2.6 for KM3NeT, and 2.5 for IceCube (see Table 1). These estimates scale linearly with thefraction of energy in protons, i.e. if the fraction is E p ∼ . E iso then the number of expectedevents will downscale accordingly. – 12 – igure 5 . The diffuse neutrino spectrum expected on Earth is shown as the blue line: the parameterset used for this result is shown above. Whilst E iso , Γ and t jet were taken as fixed parameters, thelocal rate of choked GRBs, R , was left free to vary, its value was then defined through a chi-squaredgoodness of fit test. The IceCube HESE data (black squares) were extracted from [34], while theIceCube ν µ tracks (red line) are from [35]. The dashed grey line shows the Waxmann & Bahcallupper bound [36]. We have estimated the diffuse neutrino flux from the choked GRBs and compared itto the IceCube latest results [34, 35]. Through a best fit procedure we found that, if theenergy of GRBs that goes in protons is between 10% − E iso , then IceCube dataare respectively reproduced for a local rate of chocked GRBs of R ∼ (3 − .
3) Gpc − yr − .These values are consistent with the local rate found in literature [31, 40] for standard GRBsthat have been detected by gamma-ray telescopes. This result is extremely interesting as itimplies that choked GRBs belong to the same population of the “standard” GRBs, thoughthey exploded in a very dense environment that “choked” their gamma-ray emission out fromthe source.If the isotropic energy released by the choked GRBs is much lower than the value usedin this paper (namely if E iso (cid:28) erg), the local rate of choked GRBs has to increaseaccordingly in order to reproduce the IceCube data. This result is consistent with whatfound in the literature by other authors [41, 42].In order to associate the origin of the IceCube data to choked GRBs, the “smokinggun” would be the detection of the latter at electromagnetic bands different than gammarays. Preliminary studies have shown that these sources may emit in the X-ray, optical andUV band [43–45]. We plan to extend this work and explore the possibility to perform a followup of choked neutrino sources in these bands.– 13 – Appendix - Analytic calculations for photo-meson interactions of pro-tons with thermal photons
To develop the kinematics relations used in the code, we consider three reference systems: thesource co-moving reference system, the IS reference system, and the center of mass system(these kinematic quantities are respectively indicated: plain, with the superscript IS , with anasterisk * ). We use natural units where (cid:126) = c = 1, and we indicate with β a the speed of ageneric particle ( γ a = 1 / (cid:112) − β a ), with E a its energy and with m a its mass. The invariantenergy of the system, √ s , is conserved and defined by the energy-momentum four-vector: s = ( E γ + E p ) − ( p γ + p p ) == E p + E γ + 2 E p E γ − p p − E γ − p p E γ cos θ == m p + 2 E p E γ (1 − cos θ ) (A.1)This relation holds in all reference frames. The minimum energy for a proton in the ISframe to produce a ∆ + resonance when interacting with a photon of energy E IS γ is shown inEq. (3.1): this is the threshold condition to produce the ∆ + resonance.The center of mass Lorentz factor is given by: γ ∗ = E IS p + E IS γ √ s (cid:39) E IS p √ s (A.2)Given the enormous energies of the protons with respect to the photons, the angle betweenthe direction along which the Lorentz boost of the center of mass is considered and thepropagation direction of the proton is negligible. For simplicity, we will consider the Lorentzboost along the direction of observation.To get the particles spectra in the source co-moving frame, one needs to boost theenergies in the IS by the jet’s Lorentz factor, Γ. In this reference system the followingrelation holds: E p = Γ E IS p . (A.3)Note that these quantities have not yet been corrected for the redshift due to cosmo-logical distances. A.1 Delta Resonance Region p + γ → ∆ + → (cid:40) n + π + p + π (A.4)In this region the invariant energy of the system √ s is fixed equal to the ∆ + resonancemass. Once formed, the resonance decays at rest in the center of mass system, and sincethe Lorentz boost is enormous, the decay products in the IS system continue along the samedirection of the incident proton. At the resonance, the simulation code, takes into accountthe Clebsch-Gordon coefficients, and the two decay channels are produced with the followingratio: BR( pγ → pπ )BR( pγ → nπ + ) = 2 / / . ≤ (cid:15) r < . E ∗ π = (cid:112) ( p ∗ π ) + m π = s + m π − m N √ s (A.6) E ∗ N = s − m π + m N √ s (A.7)N being the nucleon produced together with the pion (proton or neutron).In order to calculate the pion energy in the IS reference frame we must apply theLorentz transformations that are given by (considering the IS moving towards the center ofmass system along the boost direction, the x axis, with speed - β ∗ ): p IS x = γ ( p ∗ x + β ∗ E ∗ ) p IS y = p ∗ y p IS z = p ∗ z E IS = γ ( E ∗ + β ∗ p ∗ x ) (A.8)where the speed and the Lorentz factor of the resonance can be written as β ∗ (cid:39) γ = γ ∗ .Considering that the momentum of the produced particles in a two-body decay is equal to: p ∗ N = p ∗ π = (cid:113) s − s ( m π + m N ) + ( m π − m N ) √ s (A.9)we finally get the pion energy in the co-moving reference frame ( β π = 1): E π = Γ γ ∗ ( E ∗ π + p ∗ x ) = E p √ s ( E ∗ π + p ∗ cos φ ) (A.10)where p ∗ is the common value of p ∗ N = p ∗ π , and 0 ◦ ≤ φ ≤ ◦ is the angle, with respect tothe proton direction, at which the pion is emitted from the ∆ + decay in its center of massreference system. A.2 Charged Pion Decays π + → µ + + ν µ (A.11) π − → µ − + ¯ ν µ (A.12)We evaluate the kinematics of the π + decay, considering that the same relations will hold inthe π − decay, which will emerge only in the out of resonance region (see Sections A.5, A.6).The muon and the neutrino energies in the π + center of mass system (Equation A.11) are: E ∗ µ + = m π + + m µ + m π + , (A.13) E ∗ ν µ = m π + − m µ + m π + . (A.14)– 15 –or the neutrino, p ν = E ν , so the momentum of the decay products is: p ∗ ν µ = p ∗ µ + = m π + − m µ + m π + (A.15)Using Lorentz transformations, we get the energies in the co-moving reference frame,using the charged pion Lorentz factor γ π + = E π + /m π + : E µ + = γ π + ( E ∗ µ + + p ∗ µ + cos δ ) , (A.16) E ν µ = γ π + ( E ∗ ν µ − p ∗ ν µ cos δ ) , (A.17)where δ is the angle between the muon and the decaying pion, 0 ◦ ≤ δ ≤ ◦ . A.3 Neutral Pion Decay π → γ + γ (A.18)In the reference frame in which the pion is at rest, one has: p ∗ γ = E ∗ γ = 12 m π . (A.19)So the energies of the photons resulting from this decay process are: E γ = γ π ( E ∗ γ ± p ∗ γ cos ζ ) = E π ± cos ζ ) , (A.20)where γ π = E π /m π (and the decay products are emitted with an angle 0 ◦ ≤ ζ ≤ ◦ ).The ± sign refers to the opposite directions of emission of the two photons. A.4 Muon Decay µ − → e − + ¯ ν e + ν µ (A.21) µ + → e + + ν e + ¯ ν µ (A.22)In the following we refer with ν µ and ν e to indicate muon and electron neutrinos and antineu-trinos. The high-energy neutrino flux reaching the detectors includes both the ν µ producedby the charged pion decay (see Eq. (A.12)), and the ν e and ν µ produced by the decay of themuon.As a consequence of the non-conservation of parity, muons coming from pion decaysare strongly polarized: those coming from the π ± decay have a right-handed helicity. In thesystem in which the muon is at rest, the particle distribution as a function of energy and ofthe emission angle ω is given by:d N d x d ω = 14 π ( f ( x ) − f ( x ) cos ω ) (A.23)where x = 2 E ∗ l /m µ , the Bjorken x , is the fraction of the available energy ( m µ ) carried bythe lepton l ( l = e − , ¯ ν e , ν µ ), thus: 0 ≤ x ≤
1. The variable ω defines the angle between thedaughter particle and the muon spin. The functions appearing in the above equation can beobtained through the Michel parameters , which describe the phase space distribution ofleptonic decays of charged leptons. With these parameters, one can obtain:– 16 – for the muon neutrino (cid:40) f ( x ) = 2 x (3 − x ) f ( x ) = 2 x (1 − x ) (A.24) • for the electron neutrino (cid:40) f ( x ) = 12 x (1 − x ) f ( x ) = 12 x (1 − x ) (A.25)In the co-moving reference system, Equation (A.23) becomes: dNdy = 1 β µ ( g ( y, β µ ) − P µ g ( y, β µ )) (A.26)where y (cid:39) E l E µ (A.27)is the Bjorken y , namely the fraction of energy carried by the lepton l in the co-movingframe, and P µ is the muon spin projection along the direction of motion in the co-movingframe. P µ is defined as ( m π is the charged pion mass): P µ = 1 β µ (cid:18) E π ( m µ /m π ) E µ (1 − ( m µ /m π ) ) − m µ /m π ) )1 − ( m µ /m π ) ) (cid:19) . (A.28)Considering that in our case β µ (cid:39)
1, we have: • for the muon neutrino (cid:40) g ( y ) = − y + y g ( y ) = − y + y (A.29) • for the electron neutrino (cid:40) g ( y ) = 2 − y + 4 y g ( y ) = − y − y + 8 y (A.30)The energy of each neutrino in the co-moving reference frame has been obtained by multi-plying the muon energy defined by Eq. (A.16) with a value y extracted according to Equa-tion (A.26).This procedure is valid regardless of the decaying muon’s charge, and it was used toobtain the neutrino energies coming from the decays shown in Eq. (A.21) and Eq. (A.22). A.5 Delta Resonance Secondary Peak Region
In the photon energy range 0 . ≤ (cid:15) r < + is accompanied by a pion: p + γ → ∆ + + π A (A.31)The invariant energy of the system is no longer equal to the mass of the ∆ + resonance, butit’s given by Eq. (A.1). In the center of mass frame, the following kinematic relations hold: E ∗ ∆ + = s + m + − m π A √ s (A.32)– 17 – ∗ π A = s + m π A − m + √ s (A.33)Moreover, p ∗ π A is given by Eq. A.9 with π = π A and N = ∆. The π A energy in the co-movingframe can be written as (0 ◦ ≤ η ≤ ◦ ): E π A = E p √ s ( E ∗ π A + p ∗ π A cos η ) (A.34)The charged and neutral pions are produced with the following branching ratios: • /
15 of the A pions are neutral, π , which decay in two photons (Eq. (A.18)) withenergies given by Eq. (A.20); • /
15 of the A pions are positive charged, π + , which decay in a muon neutrino and apositive muon (Eq. (A.11), which in turn decays in a positron, an electron neutrino anda muon antineutrino, Eq. (A.22)). Energies for these particles are given by Eqs. (A.16),(A.17), and (A.27); • / A pions are negative charged, π − , which decay in a muon antineutrino and amuon (Eq. (A.12), which in turn decays in an electron, an electron antineutrino and amuon neutrino, Eq. (A.21)). Energies for these particles are also given by Eqs. (A.16),(A.17), and (A.27).In order to consider all the processes that can originate neutrinos and photons, we haveto take into account the ∆ + decay: ∆ + → N + π B (A.35)In the ∆ + resonance rest frame, we have: | p ∗ π B | = | p ∗N | = (cid:113) [( m + − ( m N + m π B ) )( m + − ( m N − m π B ) )]2 m ∆ + (A.36) E ∗ π B = s + m π B − m N m ∆ + (A.37)Finally, the B pion energy in the co-moving frame is given by: E π B = E ∆ + m ∆ + ( E ∗ π B + p ∗ π B cos ρ ) (A.38)Where 0 ◦ ≤ ρ ≤ ◦ is the emitting angle of the pion, and E ∆ + ∼ E p − E π A . The codeextracts randomly the branching ratios for π B , which are: • /
45 of B pions are neutral pions; • /
45 of B pions are positive charged pions; • /
45 of B pions are negative charged pions.The energies of the particles emerging from the pion decay are obtained as already describedin the above sections. – 18 – .6 Multipion Production Region When the photon energy exceeds 2 GeV, most of the energy lost by the proton ( ∼ . E p )is split equally among three pions, and the neutral and charged pions are approximatelyproduced in equal numbers (this is a simplified treatment, as found in e.g. [14, 15, 17]). Sothe single pion energy is equal to E π = 0 . E p , and the particles resulting from the pion decayare simulated according with the procedures described in the previous sections. References [1] M. G. Aartsen et al. [IceCube Collaboration],
Evidence for high-energy extraterrestrialneutrinos at the icecube detector, Science (2013) 1242856 [arXiv:1311.5238].[2] M. G. Aartsen et al. [IceCube Collaboration],
Time-integrated Neutrino Source Searches with10 years of IceCube Data, Phys. Rev. Lett. (2020) 051103 [arXiv:1910.08488].[3] S. Celli, A. Palladino and F. Vissani,
Neutrinos and γ -rays from the Galactic Center RegionAfter H.E.S.S. Multi-TeV Measurements, EPJC (2017) 66 [arXiv:1604.08791].[4] M. Ackermann et al. [Fermi-LAT Collaboration], The spectrum of isotropic diffuse gamma-rayemission between 100 MeV and 820 GeV, ApJ (2015) 86 [arXiv:1410.3696].[5] K. Murase, D. Guetta and M. Ahlers,
Hidden Cosmic-Ray Accelerators as an Origin ofTeV-PeV Cosmic Neutrinos, Phys. Rev. Lett. (2016) 071101 [arXiv:1509.00805].[6] P. M´esz´aros and E. Waxman,
TeV neutrinos from successful and choked Gamma-Ray Bursts,Phys. Rev. Lett. (2001) 171102 [arXiv:astro-ph/0103275v3].[7] T. Piran, E. Nakar, P. Mazzali and E. Pian, Relativistic Jets in Core Collapse Supernovae (2017) [arXiv:1704.08298].[8] L. Izzo et al.,
Signatures of a jet cocoon in early spectra of a supernova associated with a γ -rayburst, Nature (2019) 324 [arXiv:1901.05500].[9] Y. Q. Xue et al., A magnetar-powered X-ray transient as the aftermath of a binaryneutron-star merger, Nature (2019) 198 [arXiv:1904.05368].[10] P. B. Denton and I. Tamborra,
Exploring the properties of choked Gamma-ray Bursts withIceCube’s high-energy neutrinos, ApJ (2018) 37 [arXiv:1711.00470].[11] H. He, A. Kusenko, S. Nagataki, Y. Fan and D. Wei,
Neutrinos from Choked JetsAccompanied by Type-II Supernovae, ApJ (2018) 119 [arXiv:1803.07478].[12] A. Esmaili and K. Murase,
Constraining high-energy neutrinos from choked-jet supernovaewith IceCube high-energy starting events, JCAP (2018) 008 [arXiv:1809.09610].[13] D. Guetta, R. Rahin, I. Bartos and M. Della Valle, Constraining the fraction of core-collapsesupernovae harboring choked jets with high-energy neutrinos, MNRAS (2020) 843[arXiv:1906.07399].[14] A. M. Atoyan and C. D. Dermer,
Neutral beams from blazar jets, ApJ (2003) 79[arXiv:astro-ph/0209231v2].[15] S. R. Kelner and F. A. Aharonian,
Energy spectra of gamma-rays, electrons and neutrinosproduced at interactions of relativistic protons with low energy radiation, Phys. Rev. D (2008) 034013 [arXiv:0803.0688].[16] F. Lucarelli, Master Thesis at Sapienza University of Rome, Rivelazione di neutrini di altaenergia da sorgenti puntiformi extragalattiche con l’apparato Cherenkov sottomarinoNESTOR (1998) – 19 –
17] S. H¨ummer, M. R¨uger, F. Spanier and W. Winter,
Simplified models for photohadronicinteractions in cosmic accelerators, ApJ (2010) 630 [arXiv:1002.1310v3].[18] A. I. MacFadyen, S. E. Woosley and A. Heger,
Supernovae, Jets, and Collapsars, ApJ (2001) 410 [arXiv:astro-ph/9910034].[19] M. A. Aloy et al.,
Relativistic Jets from Collapsars, ApJ (2000) 119.[20] R. Blandford and D. Eichler,
Particle acceleration at astrophysical shocks: a theory of cosmicray origin, Phys. Rep. (1987) 1.[21] J. Alvarez-Mu˜niz, F. Halzen and D. W. Hooper,
High energy neutrinos from gamma raybursts: event rates in neutrino telescopes, Phys. Rev. D (2000) 093015[arXiv:astro-ph/0006027].[22] E. Waxman and J. Bahcall, High Energy Neutrinos from Cosmological Gamma-Ray BurstFireballs, Phys. Rev. Lett. (1997) 2292-2295 [arXiv:astro-ph/9701231].[23] S. Razzaque, P. M´esz´aros and E. Waxman, TeV neutrinos from core collapse supernovae andhypernovae, Phys. Rev. Lett. (2005) 109903 [arXiv:astro-ph/0407064v4].[24] L. Morejon, A. Fedynitch, D. Boncioli, D. Biehl and W. Winter, Improved photomeson modelfor interactions of cosmic ray nuclei, JCAP (2019) 007 [arXiv:1904.07999].[25] F. L. Villante and F. Vissani, How precisely neutrino emission from supernova remnants canbe constrained by gamma ray observations?, Phys. Rev. D (2008) 103007 [arXiv:0807.4151].[26] G. L. Fogli, E. Lisi, A. Marrone, A. Palazzo and A. M. Rotunno, Hints of theta > fromglobal neutrino data analysis, Phys. Rev. Lett. (2008) 141801 [arXiv:0806.2649v2].[27] S. Adrian-Martinez et al. [ANTARES Collaboration], Search for cosmic neutrino pointsources with four years of data of data from ANTARES telescope, ApJ (2012) 53.[28] S. Adrian-Martinez et al. [KM3NeT Collaboration],
KM3NeT 2.0 Letter of Intent for ARCAand ORCA, J. Phys. G: Nucl. Part. Phys. (2016) 084001.[29] M. G. Aartsen et al. [IceCube Collaboration], Searches for Extended and Point-like NeutrinoSources with Four Years of IceCube Data, ApJ (2014) 2.[30] N. M. Lloyd-Ronning, A. Aykutalp and J. L. Johnson,
On the cosmological evolution of longgamma-ray burst properties, MNRAS (2019) 4.[31] V. Petrosian, E. Kitanidis and D. Kocevski,
Cosmological Evolution of long Gamma-RayBursts and the Star Formation Rate, ApJ (2015) 44.[32] P. Madau and M. Dickinson
Cosmic Star-Formation History, ARAA (2014) 415[arXiv:1403.0007].[33] K. Murase, High energy neutrino early afterglows from gamma-ray bursts revisited,Phys. Rev. D (2007) 123001 [arXiv:0707.1140v3].[34] R. Abbasi et al. [IceCube Collaboration], The IceCube high-energy starting event sample:Description and flux characterization with 7.5 years of data, [arXiv:2011.03545v1].[35] J. Stettner,
Measurement of the Diffuse Astrophysical Muon-Neutrino Spectrum with TenYears of IceCube Data, Pos (2019).[36] J. Bahcall and E. Waxman,
High Energy Astrophysical Neutrinos: the Upper Bound isRobust, Phys. Rev. D (2001) 023002 [arXiv:hep-ph/9902383].[37] M. G. Aartsen et al. [IceCube Collaboration], An all sky search for three flavors of neutrinosfrom gamma-ray bursts with the IceCube neutrino observatory, The Astrophysical Journal (2016).[38] A. Albert et al. [ANTARES Collaboration],
Constraining the contribution of Gamma-RayBursts to the diffuse neutrino flux with the ANTARES dataset (2007-2017), MNRAS . – 20 –
39] S. Adrian-Mart´ınez et al. (ANTARES Collaboration),
Search for high-energy neutrinos frombright GRBs with ANTARES, MNRAS (2017) 1.[40] D. Guetta and T. Piran,
Do long duration gamma ray bursts follow star formation?, JCAP (2007) 003.[41] K. Murase and E. Waxman, Constraining high-energy cosmic neutrino sources: Implicationsand prospects, PRD (2016) 103006.[42] N. Senno, K. Murase and P. Meszaros, Constraining high-energy neutrino emission fromchoked jets in stripped-envelope supernovae, JCAP (2018) 025.[43] C. Irwin, E. Nakar and T. Piran, The propagation of choked jet outflows in power-law externalmedia, MNRAS (2019) 2844.[44] R. Perna, D. Lazzati and M. Caniello,
Electromagnetic Signatures of Relativistic Explosions inAGN Disks , submitted to ApJ Letters (2020).[45] J. P. Zhu et al.,
Neutron Star Mergers in AGN Accretion Disks: Cocoon and Ejecta ShockBreakouts , ApJL in press (2020)., ApJL in press (2020).