1 Introduction

In 2012 a diffuse flux of high energy neutrinos has been discovered by IceCube, a neutrino telescope placed in the South Pole [1]. Since then, several theoretical works have been written to interpret the flavor composition observed by IceCube [2,3,4,5,6,7,8,9]. In the first few years after the detection, it was emphasized the inconsistency between the expected background for tracks contained in the High Energy Starting Events (HESE) dataset and the number of detected tracks [3]. That paper points out a “misunderstanding of the expected background events or even more compellingly, some exotic physics which deviates from the standard scenario”. Particularly when all the HESE are considered, as in [3], the main issue is that the number of expected atmospheric tracks is larger that the number of detected tracks, favoring a null track to shower ratio for astrophysical neutrinos, with no known astrophysical interpretation so far. This issue was present in the three years dataset [10] and it lingers into the 6 years dataset [11]. In [12] the study of the flavor composition has been conducted using only events having deposited energy larger than 60 TeV (where the background is expected to be negligible), showing that the usage of this subset of data restores the compatibility between the detected flavor composition and the one expected from astrophysical production mechanisms.

In order to avoid possible problems related to the atmospheric muon background that affects HESE tracks, in [4] a new method was used, with the inclusion of the flux measured with throughgoing muons (TGM). The main reason is that the throughgoing muon flux [13] represents the cleanest way to observe the flux of astrophysical muon neutrinos, due to the high energy threshold (200 TeV) and to the absence of atmospheric muons, since they have no possibility to cross the Earth and reach the detector placed in the opposite hemisphere. On the contrary HESE tracks are likely to be largely affected by atmospheric muons, as remarked in [10, 11]. Particularly the background expected for HESE tracks is larger than the number of the detected ones. On the contrary, the flux measured using throughgoing muons is: (i) free from the contamination of atmospheric muons, since they are stopped inside the Earth; (ii) in a negligible manner contaminated by conventional atmospheric neutrinos, due to the high energy threshold of 200 TeV; (iii) marginally contaminated by prompt neutrinos, at level of 20%, as can be estimated using the signalness contained in Table 4 of [14]. However HESE are also very important, since thanks to them the diffuse flux of astrophysical neutrinos has been measured for the first time [1]. Moreover they are sensible to all neutrino flavors and the atmospheric background that affects showers is relatively low, due to the innovative veto technology.

Following the previous discussion, we formulate the main hypothesis of this work:

the expected number of astrophsyical HESE tracks is computed according to the shape and the normalization suggested by the throughgoing muon flux above 200 TeV and following the shape suggested by HESE below this energy. This is what we call “baseline model” in the rest of the work. The number of astrophysical showers, instead, is extracted from the HESE dataset, accounting for the atmospheric background (that is however small for showers).

Since this method does not rely on any assumption on the background that affects HESE tracks, it is the cleanest way to compute the observed track to shower ratio of astrophysical neutrinos and to compare it with the expected theoretical ones. We also demonstrate that the way in which we extrapolate the throughgoing muon flux below 200 TeV affects only marginally (at level of 10%) the expected number of astrophysical HESE tracks. In other words, the computation of the expected number of HESE tracks depends mostly on the muon neutrino flux already measured.

The current public available data consist of 7.5 years of High Energy Starting Events (HESE) and 8 years of throughgoing muons, therefore the analysis of the flavor composition can be much more powerful and accurate compared to the past. Motivated by this argument, in this work we re-analyze the flavor composition observed by IceCube taking into account all the most recent measurements and using a spectrum that can describe the data between \(\sim 30\) TeV and few PeV. Before the IceCube measurements, the flavor composition was already considered a powerful tool to understand the origin of high energy neutrinos [15], that remains still a mystery. Indeed, up to now, only one neutrino has been associated with an identified object [16], but all the other neutrinos remain without any confirmed counterpart.

The work is structured as follows. In Sect. 2 we discuss the method used to compute the flavor composition at Earth, how to convert it into an observable quantity and how to get information on the observed track to shower ratio. In Sect. 3 we present the result and we discuss the implications in Sect. 4. We conclude the work with Sect. 5.

2 Method

2.1 The theoretical flavor composition

The propagation of cosmic TeV–PeV neutrinos [17] includes a long oscillation phase, order of \(\sim 10^{13}\) for neutrinos of 100 TeV coming from a distance of 1 Gpc. As a consequence only the average values of the oscillation probability are astrophysical observables. In this work we compute the oscillations of astrophysical neutrinos using the “natural parametrization” introduced in [18]. It permits to compute easily the uncertainties associated to the oscillation probabilities, using three Gaussian parameters \(P_0\), \(P_1\), \(P_2\), where \(P_0 \gg P_1 \simeq P_2\). The expressions of these parameters in terms of the conventional oscillation parameters (3 mixing angles and one CP violating phase) are:

$$\begin{aligned} P_0&=\frac{1}{2} \left[ (1-\epsilon )^2 \left( 1-\frac{\sin ^2 (2\theta _{12})}{2}\right) +\epsilon ^2 -\frac{1}{3} \right] , \end{aligned}$$
(1)
$$\begin{aligned} P_1&=\frac{1-\epsilon }{2}\left( \gamma \cos 2\theta _{12}+\beta \frac{1-3\epsilon }{2} \right) , \end{aligned}$$
(2)
$$\begin{aligned} P_2&=\frac{1}{2}\left[ \gamma ^2 + \frac{3}{4}\beta ^2 (1-\epsilon )^2 \right] , \end{aligned}$$
(3)

where

$$\begin{aligned} \epsilon&= \sin ^2\theta _{13}, \quad \alpha = \sin \theta _{13}\cos \delta \sin 2\theta _{12}\sin 2\theta _{23}, \\ \beta&= \cos 2\theta _{23}, \quad \gamma = \alpha - \frac{\beta }{2}\cos 2\theta _{12}(1+\epsilon ). \end{aligned}$$

Therefore \(P_0\), \(P_1\) and \(P_2\) are not single values but distributions. Even if the distributions are not gaussian, we will use the gaussian approach that it is sufficiently accurate, as explained in [18].

Following the natural parametrization the oscillation probabilities \(P_{\ell \ell ^{'}}\) (where \(\ell \) and \(\ell ^{'}\) denotes the initial and the final neutrino flavor) are given by the elements of the following matrix:

$$\begin{aligned} \mathcal {P}=\left( \begin{array}{ccc} \frac{1}{3}+2 P_0 &{} \frac{1}{3}-P_0+P_1 &{} \frac{1}{3}-P_0-P_1 \\ &{} \frac{1}{3}+\frac{P_0}{2}-P_1+P_2 &{} \frac{1}{3}+\frac{P_0}{2}-P_2 \\ &{} &{} \frac{1}{3}+\frac{P_0}{2}+P_1+P_2 \end{array} \right) . \end{aligned}$$

The matrix acts on the vector containing the flavor composition before oscillations \(\xi ^0=(\xi _e^0,\xi _\mu ^0,\xi _\tau ^0)\) just as \(\xi =\mathcal {P}\ \xi ^0\), giving the vector of fluxes observed after oscillations, \(\xi =(\xi _e,\xi _\mu ,\xi _\tau )\). The average values and the uncertainties of the natural parameters are taken from [18] (based on the knowledge of neutrino oscillations given in [19]) and they are equal to:

$$\begin{aligned} P_0= & {} 0.109 \pm 0.005, \\ P_1= & {} 0.000 \pm 0.029, \\ P_2= & {} 0.010 \pm 0.007. \end{aligned}$$

Concerning the initial flavor composition we assume three conventional astrophysical scenarios:

  • neutrinos are produced via charged pion decay, following \(\pi ^+ \rightarrow \mu ^+ \nu _\mu \rightarrow e^+ \nu _\mu \bar{\nu }_\mu \nu _e\) or \(\pi ^- \rightarrow \mu ^- \bar{\nu }_\mu \rightarrow e^- \bar{\nu }_\mu \nu _\mu \bar{\nu }_e\). In this process the flavor composition at the source is equal to \((\xi _e^0:\xi _\mu ^0:\xi _\tau ^0)=(1:2:0)\). We do not distinguish between neutrinos and antineutrinos in this work, since the only channel to observe astrophysical antineutrinos is the Glashow resonance (see Sect. 2.2) and these events are still not observed in the present neutrino telescope. Therefore current observations are only sensitive to the sum of neutrino and antineutrino fluxes.

  • neutrinos are produced by pion decay in astrophysical environment with strong magnetic fields (\(\sim 10^5\) to \(10^6\) Gauss) [20]. Under this assumption, muons lose a significant fraction of energy before decaying, therefore high energy neutrinos are only created by the first part of the previous chain, i.e. \(\pi ^+ \rightarrow \mu ^+ \nu _\mu \) or \(\pi ^- \rightarrow \mu ^- \bar{\nu }_\mu \). In this case the initial flavor composition is equal to \((\xi _e^0:\xi _\mu ^0:\xi _\tau ^0)=(0:1:0)\);

  • neutrinos are created by the decay of neutrons, according to the process \(n \rightarrow p \ e^- \bar{\nu }_e\). In this scenario the initial flavor composition is equal to \((\xi _e^0:\xi _\mu ^0:\xi _\tau ^0)=(1:0:0)\).

Using the matrix \(\mathcal {P}\) defined above we can compute the neutrino oscillations, obtaining the flavor composition at Earth as \(\xi =\mathcal {P} \xi _0\). The flavor compositions obtained at Earth are reported in Table 1 and Fig. 1. We notice that the uncertainties on the final flavor composition are not the same for all production mechanisms; this is related to the fact that the knowledge of the natural parameters is not equally good, since \(\varDelta P_1 \gg \varDelta P_0 \simeq \varDelta P_2\).

Fig. 1
figure 1

Flavor composition of astrophysical neutrinos expected at Earth after neutrino oscillations, assuming that they are produced by pion decay (left panel), pion decay with damped muons (middle panel) and neutron decay (right panel). In the last case there is overlap between the fraction of \(\nu _\mu \) and \(\nu _\tau \) at Earth and the two distributions are not distinguishable in the figure

Full size image
Fig. 2
figure 2

On the left panel: the baseline single flavor spectrum of astrophysical neutrinos is represented with a solid green curve. It is obtained following the shape and the normalization suggested by throughgoing muons above 200 TeV [11] and shape suggested by HESE below 200 TeV [25]. The green band is related to the uncertainty on the measured normalization at 100 TeV (\(\sim 30\%\) for both HESE and TGM). The blue dotted line denotes the extrapolation to lower energies following the throughgoing muon shape. On the right panel we represent the product between flux and muon neutrino effective area, in arbitrary units

Full size image
Table 1 Flavor composition expected at Earth for the three different production mechanisms, accounting for the uncertainties on the neutrino oscillations
Full size table

2.2 The theoretical track to shower ratio expected from astrophysical scenarios

The flavor composition is not a direct observable in IceCube, since only two types of event topologies are so far identified in the modern neutrino telescopes, namely tracks and showers [21]. Neutrinos are generally detected thanks to the deep inelastic scattering [22], looking at the secondary particles produced after the interaction between neutrinos and nucleons. Tracks are produced by the interaction of \(\nu _\mu \) via charged current interaction, while showers are produced by all the other processes, i.e. charged current interactions of \(\nu _e\) and \(\nu _\tau \) and neutral current interactions of whatever neutrino. In principle there are two other processes that permit to identify \(\bar{\nu }_e\) and \(\nu _\tau \): the Glashow resonance [23] and the double cascades [24], but they are still not observed.Footnote 1

The first analysis of the flavor composition observed by IceCube has been presented for the first time in [3]. However, from the previous discussion, it follows that the observable quantity is not directly the flavor composition but the ratio between the number of tracks and the number of showers, abbreviated “track to shower ratio” in the following of this work. The analysis of the track to shower ratio has been adopted in [4] and in this work we update it, using the most recent IceCube measurements after about 8 years of exposure. In order to do that we need to include information on the incident astrophysical neutrino spectrum and on the response function of the detector.

Spectrum Up to now there are measurements of the astrophysical neutrino spectrum covering different energy ranges and sky locations. Throughgoing muons, only sensible to \(\nu _\mu \) from Northern sky above 200 TeV, suggest a hard spectrum \(\propto E^{-2.2 \pm 0.1}\) [11]. On the other hand High Energy Starting Events (HESE), that are sensitive to the all flavor flux from both hemispheres, suggest a softer spectrum between \(\propto E^{-2.5 \pm 0.1}\) [25] between 30 TeV and 3 PeV.Footnote 2 Moreover let us notice that about 90% of HESE have an energy smaller than 200 TeV while all the throughgoing muons have an energy larger than 200 TeV. Therefore it is reasonable having trust of the spectral shape suggested by throughgoing muons above 200 TeV and of the spectral shape suggested by HESE below 200 TeV. This is our baseline choice for the spectrum of astrophysical neutrinos \(\phi (E)\) and it is represented in the left panel of Fig. 2. The normalization above 200 TeV replicates the normalization of the throughgoing muon spectrum given in [11]. Let us remark that due to neutrino oscillations and standard astrophysical mechanisms, we expect the same spectral shape for all flavors at Earth. Only the normalization can change according to the production mechanism, as shown in Table 1. The idea of a two component spectrum is plausible and it has been already discussed in several theoretical works [26,27,28,29,30].

Response of the detector In order to convert the theoretical flavor composition in the observable track to shower ratio, we also need to know the response of the detector to neutrinos. This information is contained in the angle-averaged effective areas \(A_\ell ^{\text {eff}}\), provided by IceCube for each neutrino flavor \(\nu _\ell \) [1]. Using all the information discussed above, we compute the parameters \(C_i\), that denote in which way the flavor composition is modified by the detector. The parameters \(C_i\) will be used in the computation of the track-to-shower ratio and they are defined as follows:

$$\begin{aligned} C_e= & {} C_0 \int _0^\infty dE \ A_e^{\text {eff}} \phi (E), \\ C_\mu ^t= & {} C_0 \ \eta \int _0^\infty dE \ A_\mu ^{\text {eff}} \phi (E), \\ C_\mu ^s= & {} C_0 \ (1-\eta ) \int _0^\infty dE \ A_\mu ^{\text {eff}} \phi (E), \\ C_\tau= & {} C_0 \int _0^\infty dE \ A_\tau ^{\text {eff}} \phi (E), \end{aligned}$$

where \(C_0=(\sum _{\ell =e,\mu ,\tau } \int _0^\infty dE \ A_\ell ^{\text {eff}} \phi (E))^{-1}\) and \(\eta =0.8\) as in [4], denoting the fraction of the muon neutrino effective area that is connected to charged current interactions.Footnote 3

The apex t or s denotes the topology of the event, namely track or shower. The values of these parameters, obtained for our baseline spectrum \(\phi (E)\), are equal to \(C_e=0.49, C_\mu ^t=0.17, C_\mu ^s=0.04, C_\tau =0.30\). Let us notice that using an ideal detector that does not modify the flavor composition, we would obtain \(C_e = C_\tau = C_\mu ^s+C_\mu ^t\); this is not true in reality due to the different energy deposited by neutrinos having different flavors.

The way to convert the flavor composition expected at Earth in the track to shower ratio r is given by the following expression, using the previous parameters \(C_\ell \):

$$\begin{aligned} r^{\text {th}}(\xi _e,\xi _\mu )= \frac{\xi _\mu C_\mu ^t}{\xi _e C_e + \xi _\mu C_\mu ^s + \xi _\tau C_\tau }. \end{aligned}$$
(4)

Let us recall that \(\xi _e+\xi _\mu +\xi _\tau =1\), therefore there are only 2 independent variables. The theoretical track to shower ratios \(r_{\text {th}}\) obtained for the three different production mechanisms using the baseline spectrum are shown in the right panel of Fig. 4; namely pion decay (orange bars), damped muons (red bars) and neutron decay (green bars). They are equal to \(0.21 \pm 0.01\) for the pion decay scenario, to \(0.29\,\pm \,0.04\) for the damped muons scenario and \(0.11\,\pm \,0.02\) for the neutron decay scenarios. Assuming \((\xi _e:\xi _\mu :\xi _\tau )=(1:1:1)\) at Earth, we obtain a track to shower ratio equal to 0.21 using our baseline spectrum. The flavor composition is always assumed to be energy-independent in the following of the work.

Before proceeding a clarification is necessary. The spectrum suggested by HESE \(\propto E^{-2.5}\) reflects the behavior of the measurements between 30 TeV and few PeV. On the other hand in our baseline model we are only using this shape for \(E < 200 \text{ TeV }\). Limiting the HESE data to the energy between 30 and 200 TeV would result in a different spectral shape, suggesting probably a softer spectrum, since the HESE above 200 TeV are in agreement with the throughgoing muons measurements (i.e. with a hard spectrum). However the analysis of the 4 years shower dataset above 1 TeV, presented in Sect. 3 of [11], suggest an \(\propto E^{-2.48 \pm 0.08}\) spectrum at lower energies. In conclusions, there are no valid reasons to use a spectrum softer than \(E^{-2.5}\) below 200 TeV.

2.3 The track to shower ratio of astrophysical neutrinos in IceCube

The expected track to shower ratio computed in the previous section has to be compared to the detected one. In order to do that we consider the most recent HESE data, presented in [31]. This dataset consists of 113 events, including 30 tracks, 81 showers and 2 not classified events (that were already present in the previous datasets), detected after 7.5 years of exposure.

The computation to predict the observed track to shower ratio is complicated, as we need to appropriately include all sources of background. Let us notice that both in [10, 11] the expected atmospheric background for HESE tracks is larger than the observed number of tracks, when all HESE are considered. This represents an issue for the computation of the track to shower ratio, since it would indicate that no astrophysical tracks are present in the HESE sample. This information was used in [3] to claim a possible tension between neutrino oscillations and IceCube measurement. On the other hand the IceCube analysis, performed using only events above 60 TeV (where the atmospheric background is expected to be negligible), claims an opposite result compared to [3], showing that the observed flavor composition is in agreement with the damped muon scenario and compatible with the pion decay [12]. However even in this case a problem remains: a large number of tracks (more than 20 tracks in the 6 years dataset [11]) is expected between 30 TeV and 60 TeV but not detected. Both these analyses depend on the assumed background for HESE tracks in the considered energy region and they give completely different results, since they use two different energy thresholds. Let us clarify that the main source of background is represented by atmospheric muons in this case, not by atmospheric neutrinos.

In [4] it has been proposed a new method, that does not require any assumption and any knowledge of the background related to HESE tracks. This method consists in the computation of the expected number of HESE tracks, using the well measured throughgoing muon flux and the muon neutrino effective area. Although this flux is only measured above 200 TeV, we demonstrate that the extrapolations to lower energies affects only marginally the expected number of HESE tracks. In other words, the most important part of the spectrum for this kind of calculation is the one already measured.

Astrophysical tracks Using the muon effective area and our baseline spectrum, we can compute the expected number of HESE tracks as follows:

$$\begin{aligned} N_t^{\text {astro}} = 4 \pi \text{ T } \ \eta \int _0^\infty \phi (E) A_\mu ^{\text {eff}} \ dE, \end{aligned}$$

where T = 7.5 years.Footnote 4 Using the baseline spectrum represented in the left panel of Fig. 2 we obtain:

$$\begin{aligned} N_t^{\text {astro}}= 9.3^{+2.6}_{-2.3}. \end{aligned}$$
(5)

The asymmetric uncertainty is related to the asymmetric uncertainty on the normalization of the throughgoing muon flux [11]. On the other hand if we extrapolate the throughgoing muon flux at lower energies following the \(E^{-2.2}\) flux, we obtain \(N_t^{\text {astro}}=8.4^{+2.1}_{-1.9}\), i.e. a discrepancy of 10%. It means that the extrapolation of the spectrum has only a minor role in this calculation and it is confirmed by the right panel of Fig. 2, in which the differential number of expected events is represented as a function of energy. Following our baseline model, the likelihood for the astrophysical tracks is then given by a function \(\mathcal {L}_t(n_t)\) having a maximum in \(n_t=9.3\) and having the integral equal to 0.68 in the interval \(n_t=9.3^{+2.6}_{-2.3}\). We choose a function \(\mathcal {L}_t(n_t)\) consisting of two pieces of not normalized Gaussian functions \(\mathcal {G}(n_t,\mu ,\sigma )\), with \(n_t=9.3\) as splitting point. The functions are characterized by having the same mean \(\mu =9.3\) and different standard deviation \(\sigma \), namely \(\sigma = 2.6\) for \(n_t\ge 9.3\) and \(\sigma =2.3\) for \(n_t < 9.3\). Then we normalize the two pieces of Gaussian function in order to obtain a continuous function. The integral of the likelihood correctly replicate the \(1\sigma \) interval found above \(N_t^{\text {astro}}= 9.3^{+2.6}_{-2.3}\), as follows:

$$\begin{aligned} \frac{\int _{n_t^1}^{n_t^2} \mathcal {L}_t(n_t) \ dn_t}{\int _{0}^{\infty } \mathcal {L}_t(n_t) \ dn_t} = 0.68, \end{aligned}$$

where \(n_t^1=7.0\) and \(n_t^2=11.9\) are the extremes of the \(1\sigma \) region of the expected astrophysical HESE tracks. The likelihood denoting the number of astrophysical tracks is represented in the left panel of Fig. 4 using a purple curve.

Astrophysical showers Here we proceed to compute the number of astrophysical showers among the 81 showers contained in the 7.5 years HESE dataset. The showers are much less affected by atmospheric background compared to tracks; this is evident from Table 4 of [10]. In order to get the background expected in 7.5 years we scale in time the background of Table 4 of [10], that refers to an exposure of 2.7 years. We obtain that after 7.5 the expected background from conventional neutrinos plus atmospheric muonsFootnote 5 is equal to:

$$\begin{aligned} N_{s}^{\text {conv}}= 8.0^{+3.0}_{-2.5}. \end{aligned}$$

Following the same table the background associated to prompt neutrinos should be \(N_{s}^{\text {prompt}} \le 20\) at 90% C.L. in 7.5 years. On the other hand that limit was derived based on [32], in which the upper limit on prompt neutrinos was \(3.8 \times \phi _{\text {ERS}}\), where \(\phi _{\text {ERS}}\) is the theoretical flux of prompt neutrinos calculated in [33]. Recently the upper limit on prompt neutrinos has been improved, reaching the level of \(1.06 \times \phi _{\text {ERS}}\) in [14]. Therefore after 7.5 years we expect that prompt neutrinos give at the best fit a null contribution to HESE showers and they can contribute at level of:

$$\begin{aligned} N_{s}^{\text {prompt}} < 5.6 \quad \text{ at } \text{90 }\%\, \hbox {C.L.} \end{aligned}$$

The likelihoods for the conventional background \(\mathcal {L}_s^{\text {conv}}\) and for prompt neutrinos \(\mathcal {L}_s^{\text {prompt}}\) are represented in Fig. 3. For conventional atmospheric showers the likelihood is constructed in order to obtain the integral equal to 0.68 in the interval \(8.0^{+3.0}_{-2.5}\); this function is constructed using two pieces of Gaussian functions, as explained below for the likelihood of astrophysical tracks. For atmospheric prompt neutrinos, instead, we consider an exponential function (being the experimental best fit equal to 0) under the assumption that the integral of the this likelihood is equal to 0.9 between 0 and 5.6.

Fig. 3
figure 3

The expected contribution of the atmospheric background to HESE showers in 7.5 years of exposure. The solid curve denotes the contribution expected from atmospheric neutrinos and atmospheric muons [10], while the dashed curve denotes the contribution expected from prompt neutrinos [10, 14]

Full size image

Now we have all the ingredient to compute the likelihood for the number of astrophysical showers \(n_s\) that contribute to HESE showers, according to the following equations:

$$\begin{aligned} \mathcal {L}_s(n_s) \propto&\ \int _0^\infty dn_s^{\text {conv}} \int _0^\infty dn_s^{\text {prompt}}\nonumber \\&\times \left( n_s+n_s^{\text {conv}}+n_s^{\text {prompt}}\right) ^{N_s}\nonumber \\&\times \exp \left[ -\left( n_s+n_s^{\text {conv}}+n_s^{\text {prompt}}\right) \right] \nonumber \\&\times \mathcal {L}_s^{\text {conv}}\left( n_s^{\text {conv}}\right) \ \mathcal {L}_s^{\text {prompt}}\left( n_s^{\text {prompt}}\right) , \end{aligned}$$
(6)

where \(N_s=81\) denotes the number of observed showers. The resulting number of astrophysical showers is equal to:

$$\begin{aligned} N_s^{\text {astro}}= 73.0^{+9.5}_{-10.2}, \end{aligned}$$
(7)

and the likelihood function \(\mathcal {L}_s(n_s)\), denoting the number of astrophysical showers, is represented in the left panel of Fig. 4 using a yellow curve.

The observed track to shower ratio can be computed using the following expression:

$$\begin{aligned} \mathcal {L}_r^{\text {obs}}(r) \propto \int _0^\infty \mathcal {L}_t(r \ n_s) \ n_s \ \mathcal {L}_s(n_s) \ dn_s, \end{aligned}$$
(8)

after the changing of variable \(n_t= r n_s\).

Fig. 4
figure 4

On the left panel: likelihood for the number of astrophysical HESE tracks (purple curve) and astrophysical HESE showers (yellow curve) after 7.5 of exposure with IceCube. On the right panel: expected astrophysical track to shower ratio for different production mechanisms (orange bars for pion decay, red bars for damped muons, green bars for neutron decay) compared to the track to shower ratio derived by IceCube observations, using the throughgoing muon flux + HESE (blue curve) and HESE above 60 TeV only (yellow region). The shaded regions show the 1\(\sigma \) interval

Full size image

3 Results

HESE + throughgoing muons The likelihood relative to the observed astrophysical track to shower ratio [see Eq. (8)], accounting for the throughgoing muon flux, is reported in the right panel of Fig. 4 as a blue curve. The observed track to shower ratio \(r_{\text {astro}}^{\text {HESE + TGM}}\) is equal to:

$$\begin{aligned} r_{\text {astro}}^{\text {HESE + TGM}} = 0.12 \pm 0.04, \end{aligned}$$

and the shaded blue region of the right panel of Fig. 4 denotes the \(1\sigma \) interval.

In order to define the compatibility between the observations and the theoretical expectations we use a statistical treatment, defining the function

$$\begin{aligned} D(\delta )^i = \int _0^\infty \mathcal {L}^{\text {obs}}_r(r+\delta ) \mathcal {L}_r^{\text {th},i}(r) dr, \end{aligned}$$

where \(\mathcal {L}^{\text {obs}}_r\) is the observed likelihood defined above and \( \mathcal {L}_r^{\text {th},i}\) is the theoretical track to shower ratio expected from the production mechanism i. Then we calculate at how many \(\sigma \) the value \(\delta =0\), i.e. the null distance between these two distributions (i.e. \(\mathcal {L}^{\text {obs}}\) and \(\mathcal {L}_r^{\text {th},i}\)), is disfavored. In order to do that we cut the distribution D in two points, at equal height, defining:

$$\begin{aligned} {\delta _1}= 0,\quad {\delta _2} \rightarrow D({\delta _2})=D(0), \end{aligned}$$

and we compute

$$\begin{aligned} \mathcal {I}^i=\frac{\int _{\delta _1}^{\delta _2} D(\delta )^i \ d\delta }{\int _{-\infty }^{\infty } D(\delta )^i \ d\delta }. \end{aligned}$$

After checking that the distributions \(D^i\) are in good approximation normally distributed, we convert the result of the previous integral in a number of \(\sigma \), using a Gaussian approach (i.e. \(0.68 \rightarrow 1\sigma , \ 0.95 \rightarrow 2\sigma , \ 0.997 \rightarrow 3\sigma \) etc.). We find that:

  • the neutron decay scenario is the best option, resulting well compatible with the observed track to shower ratio;

  • the pion decay scenario is disfavored at 2\(\sigma \);

  • the damped muon scenario is disfavored at \(2.6\sigma \).

HESE only: we also show the result that comes out from the conventional procedure, considering HESE above 60 TeV and accounting for the background. In the 7.5 years HESE dataset we find 19 tracks and 51 showers above 60 TeV.

Table 2 Summary of the results. The expected track to shower ratio for each production mechanism is reported and compared with the astrophysical track to shower ratio obtained using HESE + throughgoing muons and HESE only. The tension between observations and expectations is quoted in terms of number of sigma
Full size table

Scaling the background reported in Table 4 of [10] with the exposure, the expected background consists of \(\sim \) 6 tracks and \(\sim \) 2 showers. Following the same procedure reported in Sect. 2.3 to subtract the background and to compute the track to shower ratio, we obtain:

$$\begin{aligned} r_{\text {astro}}^{\text {HESE only}} = 0.25^{+0.11}_{-0.08}. \end{aligned}$$

The likelihood is reported in the right panel of Fig. 4 using a yellow curve, showing also the \(1\sigma \) region as a shaded yellow region. We discuss in the next section why this result is different compared to the one obtained using the throughgoing muon flux. A summary of the results is reported in Table 2.

4 Discussion

4.1 An indication of neutron decay?

The throughgoing muon flux is based on 36 tracks detected above 200 TeV after 8 years of exposure [11]. This dataset is free from atmospheric muons and negligibly contaminated by atmospheric neutrinos. It may be slightly contaminated by prompt neutrinos but it is for sure dominated by an astrophysical signal. On the other hand the 19 HESE tracks, detected after 7.5 tracks, are expected to be contaminated at level of \(\sim 30\%\) by atmospheric muons and atmospheric neutrinos. Moreover the statistic of HESE tracks is a factor 2 smaller that the statistic of throughgoing muons. For these reasons the analysis of the astrophysical track to shower ratio, performed using the throughgoing muon flux, is plausibly more accurate compared to the one performed using HESE only.

Although the pion decay has been always considered the best mechanism for the production of high energy neutrinos, the neutron decay hypothesis to explain TeV–PeV astrophysical neutrinos is plausible and it was already discussed in literature in [34]. However this paper admits that the model should be fine tuned, in order to reproduce the observe data. This is simply related to the processes involved; indeed neutrinos produced by pion decay take about 1/20 of the primary proton’s energy, while neutrinos produced by neutron decay would take about 1 / 1000 of the neutron energy. Therefore there is a factor \(\sim 50\) of difference between the energy budget available for neutrinos from pion decay versus neutrinos from neutron decay. A mechanism able to suppress the photopion production, inside the source, would be required to suppress the neutrino flux expected from the conventional pion decay. This goal may be reached with a particular choice of the target photon field inside the source (for example choosing a peculiar temperature). Another possible criticism to the neutron decay scenario would be the over production of events due to the Glashow resonance [23] due to the fact that the flux at Earth were dominated by \(\bar{\nu }_e\) in this scenario [35, 36]. On the other hand in [37] it has been shown that the spectral index and the energy cutoff play a role more important than the production mechanism in the evaluation of the expected number of resonant events. In fact, even assuming a neutron decay scenario, an energy cutoff below 6.3 PeV would nullify the possibility to observe resonant events.

We also cross checked our procedure computing the number of expected astrophysical showers (given the neutron decay as production mechanism) and comparing it with the number of astrophysical showers resulting after the background subtraction. Since assuming the neutron decay scenario the flavor composition at Earth would be roughly \((\xi _e:\xi _\mu :\xi _\tau )=(2:1:1)\), the expected number of astrophysical HESE showers can be evaluated as follows:

$$\begin{aligned} N_s^{\text {astro}} = 4 \pi \text{ T } \ \int _0^\infty \phi (E) \left[ 2 A_e^{\text {eff}}+(1-\eta ) A_\mu ^{\text {eff}} +A_\tau ^{\text {eff}}\right] \ dE, \end{aligned}$$

obtaining \(N_s^{\text {astro}}=69.3\) after 7.5 years of exposure. This result is in very good agreement with the \(\sim 73\) astrophysical showers found using the background subtraction (see Eq. 7), that is a completely independent method.

In addition to, we notice that this track to shower ratio is also compatible with the neutrino decay scenario [38], assuming normal hierarchy. This scenario has been already investigated in the past [39,40,41].

As a last remark, we notice that the normalization of the throughgoing muon flux \(\phi ^{100}_\mu \) at 100 TeV is \(1.01^{+0.26}_{-0.23}\) (see Sect. 4 of [11]) in the usual units of \(10^{-18} \mathrm \ GeV^{-1} \ cm^{-2} \ s^{-1} sr^{-1}\), while the normalization of the all flavor HESE flux \(\phi ^{100}_{3f}\) at the same energy and in the same units is \(6.7^{+1.1}_{-1.2}\) [25]. Let us remark that the normalization of the throughgoing muon flux does not require any assumption on the flavor composition, since this analysis is only sensible to muon neutrinos. On the other hand the normalization of the HESE flux requires an assumption on the production mechanism. Therefore the analysis that follows should be taken as a check of consistency, not as a conclusive result. The ratio between \(\phi ^{100}_\mu \) and \(\phi ^{100}_{3f}\) at 100 TeV is therefore equal to:

$$\begin{aligned} \frac{\phi ^{100}_\mu }{\phi ^{100}_{3f}} = 0.15 \pm 0.05. \end{aligned}$$

This can be compared to theoretical flavor fraction of muon neutrinos expected in the case of neutron decay, obtaining:

$$\begin{aligned} \xi _\mu = P_{e\mu } = 0.225 \pm 0.03. \end{aligned}$$

Also this rough estimation, based only on the flux at 100 TeV, supports the neutron decay hypothesis.

4.2 Alternative interpretation: a complex spectral shape

Excluding pessimistic hypotheses, such as a large misidentification of tracks in showers and/or the atmospheric background not under control, the result presented above can be also interpreted in a different way, that we are going to discuss in this section. Summarizing, we have seen that the number of HESE tracks obtained from our baseline model is \(N_t^{\text {astro}}= 9.3^{+2.6}_{-2.3}\) in 7.5 years (see Eq. 5). Considering only the energy range above 60 TeV, the expectation becomes roughly 80% of the previous number (see Table 4 of [10]). On the other hand 19 HESE tracks have been detected and 6 of them are expected to be background events, resulting in 13 astrophysical tracks. Therefore this number is approximately a factor 2 larger than the number expected from the baseline model, that is based on the throughgoing muon flux.

This discrepancy may suggest a spectrum much more complex than a simple power law flux. Let us assume, as an example, that the true astrophysical flux looks like the toy spectrum represented in Fig. 5 with a dashed black curve. Under the hypothesis of pion decay, this flux would give rise to \(\sim 10\) HESE tracks above 60 TeV. Considering the background, we would obtain a total of 16 HESE tracks expected versus 19 observed, with a non significant tension accounting for the Poissonian uncertainty. The spectrum of Fig. 5 would suggest a flux above 200 TeV harder than \(E^{-2}\). In [42] we found indication for a hard throughgoing muon spectrum, characterized by \(E^{-\alpha }\) with \(\alpha =1.91 \pm 0.20\). Nowadays, with the increasing of the exposure, the data seems to prefer a softer spectrum, characterized by \(\alpha \simeq 2.2 \pm 0.1\), as reported in [11]. On the other hand the hypothesis of a more complex spectral shape is worthy of being investigated, since a power law neutrino spectrum is only expected from sources in which neutrinos are produced via proton-proton interaction [43], while it is not compatible with neutrinos produced in sources dominated by \(p\gamma \) interaction [44]. For example the toy spectrum represented in Fig. 5 may be produced by two different populations of sources dominated by \(p\gamma \) interaction.

For the sake of completeness, we need to remark that all the paper is based on the assumption that the flavor composition is energy independent between roughly 10 TeV and 10 PeV. In environments with strong magnetic fields, the flavor composition may be energy dependent going from the pion decay scenario to the damped muon scenario with the increasing of the neutrino energy. However this scenario goes in the opposite direction compared to our findings, therefore it cannot be used as a possible explanation for our results.

Fig. 5
figure 5

Example of non trivial astrophysical neutrino energy spectrum characterized by two peaks. In the same figure 6 years of HESE and 4 years of cascades are reported according to [11]

Full size image

5 Conclusion

In this work we investigate the track to shower ratio suggested by astrophysical neutrinos after 8 years of observations in IceCube. We compare it with the ones expected from three theoretical scenarios, namely the pion decay, the damped muons and the neutron decay. We use the natural parametrization to compute the oscillations of astrophysical neutrinos and we take advantage of the most recent IceCube measurements, by using a broken power law spectrum that is in agreement with all the data between \(\sim 30\) TeV and few PeV. Moreover we used the flux of throughgoing muons to evaluate the expected number of astrophysical HESE tracks and we take into account that background that can affect HESE showers. We conclude that the observed track to shower ratio is fully consistent with the neutron decay scenario while it is in tension at level of \(2\sigma \) and \(2.6\sigma \) with the standard pion decay scenario and with the damped muons one, respectively. This result differs from [3], in which a null track to shower ratio was favored using all HESE, although also in that case the neutron decay scenario was the best option among the standard astrophysical mechanisms. It is also different compared to [25], in which only events above 60 TeV are considered and the damped muon scenario is the best candidate mechanism. In addition to our use of the most updated datasets, the main difference is that our work does not rely on the background that affects astrophysical HESE tracks, that represents the biggest source of uncertainties in the computation of the track to shower ratio. To tackle this problem the number of expected astrophysical HESE tracks is computed thanks to the well measured throughgoing muon flux, showing that the extrapolation below 200 TeV plays only a minor role. In principle all these three methods should give the same results; these differences may stem from the uncertainties of the poorly known atmospheric muon background.

Another possibility is that the spectrum of astrophysical neutrinos is much more complex than a power law flux. We have shown that an energy spectrum with two peaks may alleviate the tension between HESE and throughgoing muons, partially recovering the compatibility with the pion decay scenario.

Both the previous possibilities are worthy of being investigated. If the indication for a neutron decay scenario were confirmed and improved in the future, it would have an impact on the models that aim to explain the high energy neutrino emission, given the fact that in most of the models neutrinos are expected to be produced by pion decays and not by neutron decay, although the last possibility has been already considered in the scientific literature. On the other hand if the spectrum is much more complex than a power law flux, this may also have an impact on several aspects related to the interpretation of astrophysical neutrinos and to the multi-messenger connection with the diffuse flux of \(\gamma \)-rays.