A Multi-Component Model for the Observed Astrophysical Neutrinos

We propose a multi-component model for the observed diffuse neutrino flux, including the residual atmospheric backgrounds, a Galactic contribution (such as from cosmic ray interactions with gas), an extra-galactic contribution from pp interactions (such as from starburst galaxies) and a hard extragalatic contribution from photo-hadronic interactions at the highest energies (such as from Tidal Disruption Events or Active Galactic Nuclei). We demonstrate that this model can address the key problems of astrophysical neutrino data, such as the different observed spectral indices in the high-energy starting and through-going muon samples, a possible anisotropy due to Galactic events, the non-observation of point sources, and the constraint from the extragalatic diffuse gamma-ray background. Furthermore, the recently observed muon track with a deposited energy of 4.5 PeV might be interpreted as evidence for the extragalactic photo-hadronic contribution. We perform the analysis based on the observed events instead of the unfolded fluxes by computing the probability distributions for the event type and reconstructed neutrino energy. As a consequence, we give the probability to belong to each of these astrophysical components on an event-to-event basis.


I. INTRODUCTION
In the last few years, the IceCube experiment has opened a new era in neutrino astronomy, providing the first evidence for high-energy astrophysical neutrinos [1]. The analysis used by IceCube to claim the discovery relies on events with the interaction vertex contained in the inner detection volume, also known as high energy starting events (HESE) [1][2][3]. For this analysis, a veto based on the outer detector layer has been implemented to reduce the impact of atmospheric down-going neutrinos and muons; for instance, it rejects (mainly muon) neutrinos accompanied by an atmospheric muon. The HESE analysis includes both muon tracks and showers (cascades) from both the Southern and Northern hemispheres, with slightly different sensitive energy ranges. An alternative used before the HESE discovery -and revived recently -is the detection of through-going muons originating from interactions of ν µ from the Northern hemisphere (which are free of atmospheric muon backgrounds) outside the detector [4,5]. Similar analysis techniques are used in ANTARES [6] and the coming KM3NeT experiment [7].
The origin of the observed astrophysical neutrinos remains a mystery. While the directional event distribution is isotropic at a first glance -pointing towards an extragalactic origin of most of the astrophysical neutrinos -a more detailed analysis reveals a ∼ 2σ excess close to the Galactic plane when only HESE events with deposited energy above 100 TeV are considered [8]; see e.g. Refs. [9][10][11][12][13] for recent discussions. Concerning the spectrum, the through-going muon dataset suggests an hard spectrum close to E −2 ν where the HESE sample is described by a much softer spectrum, close to E −3 ν [8]. Consequently, the single power law hypothesis of the astrophysical flux has been questioned [13][14][15][16][17].
Stacking searches using gamma-ray catalogues of popular candidate classes, such as Gamma-Ray Bursts (GRBs) [18] and Active Galactic Nuclei (AGN) blazars [19], have revealed that the observed flux cannot be dominated by these object classes. For example, the contribution from blazars has been found to be smaller than about ∼ 20% of the observed diffuse flux [19]; see also [20].
More generically, no point sources have been resolved so far, which indicates that most of the observed extragalactic neutrinos ought to come from abundant sources with low luminosities [21][22][23]. If a large portion of the energy range is to be described, these neutrinos have to follow a power law with a spectral index between about E −2 ν and E −2.5 ν as expected from Fermi shock acceleration for the primary cosmic rays. Such a neutrino spectrum following the primary spectrum is expected for Ap or pp interactions.
Starburst galaxies are an ideal source candidate for that [24][25][26] (we will discuss other alternatives later), which however face two other challenges: a) the neutrino spectrum must not be much softer than E −2 to avoid the constraint from the observed diffuse extragalactic gamma-ray background [27,28] and b) cosmic-rays accelerated in starburst galaxies are unlikely to produce neutrinos with energies much larger than PeV energies. However a muon track with a deposited energy of 4.5 PeV and a likely neutrino energy 10 PeV has been recently detected by IceCube [5] -which requires a 100 PeV primary.
Such high-energy events could come from rarer sources producing neutrinos by Aγ or pγ interactions, which typically have hard spectra which can easily peak at PeV energies and beyond, see e.g. Ref. [29]. One example, which has been recently drawing a lot of attention, are neutrinos from (jetted) Tidal Disruption Events (TDEs) [30][31][32][33][34][35][36], which we will use as a prototype. Compared to other sources with similar spectral properties, such as AGN blazars, we will see that the TDE hypothesis can be tested by multiplet searches in the near future. In addition, a self-consistent description with ultra-high energy cosmic rays (UHECRs) can be performed [35], i.e., these neutrinos may come from the sources of the cosmic rays at the highest energies. Note that in this case the primaries include nuclei heavier than protons, as indicated by recent results of the Pierre Auger collaboration [37].
Since it is clear that a single (astrophysical) flux component cannot address all open questions, we propose a multi-component model to draw a self-consistent picture of the observed astrophysical neutrino flux which can satisfy all of these constraints. Of course, a relevant question is how many components do we really needwhich we are going to address quantitatively. Compared to many earlier works, we start off with the information on the observed events and reconstruct the probability distribution for the reconstructed neutrino energy and interaction type; see Refs. [38,39] for other such examples.
The structure of the paper is as follows: In Section II, we present the most recent IceCube dataset, namely the six years HESE dataset and the eight years throughgoing muon dataset, and we describe our procedure to obtain the probability distribution function of the reconstructed neutrino energy; for technical details, see Appendix A. We furthermore present in Section III our multi-component model, analyzing separately the contributions of atmospheric neutrinos, Galactic neutrinos, extragalactic pp neutrinos (such as neutrinos from starburst galaxies) and extragalactic pγ neutrinos (such as neutrinos from TDEs or AGNs) from a theoretical point of view. We show the results of the fit in Section IV, where we also discuss how many components are really needed, and what we can learn in the future. An eventbased assignment to the different components (probabilities) can be found in Appendix B. Finally, we summarize in Section V.

II. METHODS
In this section, we describe the used datasets and our energy reconstruction technique.

A. The IceCube datasets
Here we introduce the latest two datasets provided by the IceCube collaboration after six years and eight years of data taking [8], which we use for our analysis: high energy starting events (HESE) dataset and the throughgoing muon dataset.

High Energy Starting Events (HESE)
The HESE are the events with neutrino interaction vertex contained in the detector volume; this analysis has been used by IceCube to claim the evidence of an extraterrestrial flux of high energy neutrinos [1]. The most recent dataset includes 2078 days (5.7 years) of operation with 82 HESE detected [8]: 22 muon tracks, 58 showers and two events produced by a coincident pair of background muons from unrelated cosmic ray air showers, that have been excluded from the analysis. These events are characterized by a deposited energy larger than 30 TeV, and the most energetic HESE has deposited an energy of 2 PeV into the detector.
The flux attributed to astrophysical neutrinos has been, to a first approximation, described by a power law spectrum and an isotropic distribution. The per-flavor flux has been given by with F = 2.5 ± 0.8 and α = 2.92 +0. 33 −0.29 [8]. 1 The HESE dataset is dominated by events coming from the Southern hemisphere because above several hundred TeV the absorption probability of neutrinos crossing Earth becomes greater than 50% [40]. In addition, re-call that the veto method of atmospheric neutrinos only works for downgoing neutrinos, i.e., for neutrinos from the Southern hemisphere.
The HESE are isotropically distributed when the full dataset is considered. On the other hand, when the analysis is limited to events above ∼100 TeV, an anisotropy appears, with an accumulation of events close to the Galactic plane which is visible by bare eye; see Fig. 1. The present statistical significance of this accumulation of events is about 2σ. 2 This slight excess of events close to the Galactic plane was already present in older IceCube datasets and has been already discussed in 1 In IceCube analyses, usually the equipartition of the flavors is assumed, which is roughly expected from high energy neutrinos produced by pion decays after flavor mixing. 2 For one (isotropic) event, the probability to end up in the Galactic region is equal to sin b 0.26. Using this information, we derive the probability that a certain number of events is contained in this region, similarly to the calculation performed in Ref. [20]. We obtain that the probability that four or more events are found in the Galactic region is 72%, whereas the probability that eight events or more are found there is only 7% (which translates into a 93% confidence level, or, indeed 1.8σ for the eight events found there).
Refs. [11][12][13]41]. In this work, we will also take into account the possible existence of a Galactic component of high energy neutrinos.

Through-going muons (TGM)
The IceCube collaboration has collected a sample of charged current events produced by up-going muon neutrinos [8] from 2009 to 2017 . The through-going muons come from muon neutrinos interacting outside the detector. Thus, in order to avoid the atmospheric muon background, the field of view must be restricted to the Northern hemisphere such that atmospheric muons are absorbed in Earth.
The high energy sample (with muon reconstructed energy above ∼ 200 TeV) contains 36 such events; a purely atmospheric origin of them is excluded at more than 5σ. The most energetic event corresponds to a reconstructed muon energy of 4.5 PeV.
The corresponding cosmic muon neutrino and antineutrino flux has been obtained with a power law fit to the data: The parameters are F µ = 1.01 +0. 26 −0.23 and α = 2.19 ± 0.10. This analysis is sensitive mostly to muon neutrinos and antineutrinos plus a small contribution from the τleptons that decay into muons.
The spectral discrepancy between the HESE and TGM datasets has a significance of about 3σ and will be investigated in this work.

B. Neutrino energy reconstruction
Here we summarize our energy reconstruction mechanism, for details, see Appendix A.
As the general ansatz we start with the list of events, sorted by topology (shower-like or track-like) and dataset (HESE or TGM); see Tabs. 10, 11, and 12. For each of these events, apart from the toplogy, the sky position (and its uncertainty) is known, as well as the energy deposited in the detector in form of secondary particles and, eventually, light. The list of events contains a residual background from atmospheric neutrinos and muons passing the veto.
The final task will be to determine for each event i the reconstructed neutrino energy distribution C i (E ν ) including the probability that a shower was created from a ν e , ν τ or neutral current interaction, or a track was created from a ν µ or ν τ induced muon. Our Monte Carlo procedure to obtain this probability distribution is described in detail in Appendix A. Note that (in terminology) we go from incident (true) neutrino energy E ν to the deposited energy E dep in the detector, which leads to a distribution for the reconstructed neutrino energy E ν -which we label the same as the incident energy for the sake of readability. This reconstructed energy distribution can be used to construct binned error bars on the flux in terms of the reconstructed neutrino energy for arbitrary datasets if one deconvolves the total event rate in each bin with the effective area. This leads to similar results compared to the IceCube procedure at the highest energies (see e.g. data points in Fig. 3); however, the residual atmospheric background will show up explicitly at low energies. As another difference, we separate tracks and showers since the track sample is more affected by atmospheric backgrounds than the shower sample, which will be relevant for the construction of our multi-component model.

III. THE MULTI-COMPONENT MODEL
In this section we introduce our model, which includes (i) the atmospheric background (conventional and prompt neutrinos, atmospheric muons), since we will include it in the reconstructed data points; ii) Galactic neutrinos generated in pp collision in the Galactic disk plus Galactic neutrinos from point sources; iii) extragalactic neutrinos from pp (or Ap) interactions, such as neutrinos from starburst galaxies; and iv) neutrinos from pγ (or Aγ) interactions, such as neutrinos produced in TDEs, AGNs, or GRBs. We summarize these different components in Tab. I.

A. Residual atmospheric backgrounds
Atmospheric neutrinos and atmospheric muons passing the veto system produce a residual background for the high energy neutrino detection. The conventional atmospheric (neutrino and muon) background follows the atmospheric muon spectrum ∝ E −3.7 (determined by the initial cosmic-ray spectrum modified by interactions of pions and kaons), whereas prompt neutrinos at higher energies are characterized by an E −2.7 spectrum (following the initial cosmic-ray spectrum because the secondaries decay faster than they interact). The flavor composition of atmospheric neutrinos is to a good approximation (ν e : ν µ : ν τ ) (0 : 1 : 0) in the first case, where pion decays dominate and kaon decays give a smaller but not negligible contribution. In pion decays, muon neutrinos and muons are produced. The secondary muons have no time to decay before reaching the surface of the earth, therefore the conventional background is essentially given by ν µ andν µ . As far as the prompt neutrinos are concerned, the flavor composition is given by (ν e : ν µ : ν τ ) (1 : 1 : 0) since charmed meson decays dominate here.
Atmospheric muons contribute to the background from the Southern hemisphere (down-going events). Note that we will not consider the actual atmospheric fluxes, but a residual background passing the veto systems with a normalization describing observations.
The HESE track sample is much more affected by the atmospheric backgrounds than the HESE shower sample because atmospheric muons dominate this background (in the Southern hemisphere) and the atmospheric electron neutrino flux, leading to showers, is suppressed compared to the muon neutrino flux (in the Northern hemisphere). In order to avoid a bias on the shower sample, we assume different backgrounds for the φ µ and φ e fluxes  Normalization" refers to the (free) normalization parameter used in the fit (unless it is fixed), the column "α" to the approximate spectral index of the incident neutrino flux ∝ E −α ν , the column "Energy range" to the approximate energy range where the spectrum is found to dominate (which depends on the sky direction and event topology), and the column "Angular" to the assumed rough angular distribution. In the last column, our baseline model for the spectral shape is marked boldface.
in our model, reflecting the different systematics of these topologies: Here F e atm and F µ atm are two free (systematics) parameters of the model, whereas we fix F prompt = 0.66 to the present upper limit of the prompt neutrino flux per flavor [8]. This choice is motivated as conservative background estimate, consistent with the total background expected for shower-like events in IceCube [8]. We will see later that the prompt contribution does not affect our results because it is about an order of magnitude smaller than the astrophysical contributions.
Let us re-call that the background expected by Ice-Cube is in tension with the observations, as they expect 33.9 +9.3 −6.4 background tracks (see Tab. 3 of [42]), whereas only 22 tracks have been observed [8]. Note that in addition the astrophysical signal contributes to the tracklike events (about 20% of the astrophysical neutrinos should produce tracks, assuming pion decay as production mechanism [43]). This implies about 40 tracks are expected instead of the 22 observed, which represents a discrepancy of about (40 − 22)/ √ 22σ 3.8σ; this tension between expectation and observation was already present in the older datasets (see e.g. Tab. 4 of [2] for the three-year dataset).
In our analysis, the two systematics parameters F µ atmo and F e atmo are to be marginalized over, which means that no particular input is assumed for these. We will verify that, at the end, we will reproduce the total background event expected by IceCube [8], namely 40.8 +13.5 −8.3 events, reasonably well. While this procedure does not affect our result at the highest energies or for the astrophysical components qualitatively, the obtained values for the systematics parameters can be interpreted with respect to a possible mis-identification of events and with respect to the possible origin of the discussed track discrepancy.

B. Galactic component
The presence of a Galactic component, that could affect especially the events coming from the Southern hemisphere, has been already proposed in [11-13, 15, 44], based on the IceCube observations. A guaranteed (diffuse) flux comes from the interactions of cosmic rays with gas; however, this contribution is typically small in models in which the local cosmic ray density is assumed to be representative for the whole Milky Way [10,45]. If, however, a possible inhomogeneous cosmic ray distribution in our galaxy is taken into account, higher contributions can be obtained. In [46] this flux has been calculated taking into account different cosmic ray distributions; the most optimistic case corresponds to the so-called KRAγ model [47], in which the spectral index of cosmic rays is a function of the distance from Galactic center, namely harder close to the Galactic center and softer far away from it. In this scenario the expected number of HESE in IceCube are about six after 5.7 years of exposure. The additional contribution from Galactic point sources is considered in [41]. The expected number of events from point sources is ∼ 3 with the present exposure. Therefore the total expected number of Galactic events must not exceed 9 HESE in about 6 years of detection. These theoretical expectations are perfectly compatible with the most recent experimental limits on the Galactic flux provided by ANTARES [48] and the IceCube [49].
Following these considerations, we propose a Galactic flux that is in agreement with the present knowledge. Considering that the observed flux of Galactic cosmic rays observed at Earth is E −2.7 , we would expect that high energy neutrinos produced by pp interactions between Galactic cosmic rays and gas will follow the power law spectrum of the primary particles -being somewhat harder from the interactions ∝ E −2.6 [50]. We assume this shape in our multi-component model. Note that in a more aggressive scenario, such as the KRAγ model, the neutrino spectrum could be harder. We will also verify that our result satisfies maximum expectation from both theory and experiment.
We use an energy cutoff in the neutrino spectrum at the corresponding knee of the cosmic rays at about 3 PeV for protons and, correspondingly, 150 TeV for neutrinos 3 , considering that the average energy of neutrinos from pp interaction is 5% of the primary proton energy [50]. This cutoff is compatible with recent multi-component models for cosmic rays at around the knee [51] (see Fig. 1 in Ref. [45] for the neutrino flux), assuming that the proton component dominates the neutrino production; it is also compatible with theoretical arguments on the maximal acceleration energy of Galactic cosmic rays [52].
Consequently, the diffuse Galactic flux of high energy neutrinos is assumed to follow in which the only free parameter is the normalization F gal . In order to satisfy the present theoretical limits, we evaluate the expected number of events as using the average HESE effective area [1]. This is a good choice since different intervals of terrestrial declination are covered in the inner Galaxy region, namely , which implies that F gal 1.1 in order not to overshoot the 9 Galactic events expected in 6 years considering the most optimistic model, as described above. We will verify if our final result is in agreement with this constraint in the following of the paper.
Finally, note that we have multiplied Eq. (5) with the full solid angle 4π, which does not mean that the Galactic flux is present in the whole sky. This choice however allows for an easy comparison with the IceCube measurement, for which usually the average flux is given per steradian. In our model we take into account that the Galactic flux is present in the region of Galactic latitude |ϕ| ≤ 5 • and of longitude |λ| ≤ 45 • (inner galaxy), whereas it is not present outside. Therefore we can write the normalization in an alternative manner for the Galactic components, as follows: where Ω = 4 sin |ϕ max | · |λ max |, with |ϕ max | = 5 • and |λ max | = 45 • . 3 We assume an exponential energy cutoff for the proton spectrum, which implies, approximately, the square root of the exponential cutoff for the neutrino spectrum [50]. 4 We have checked that using the Southern sky effective area instead of the all-sky one, the result is only marginally affected; the expected number of events is 10% larger considering an E −2.6 spectrum with the energy cutoff at 150 TeV. This is naturally expected by the fact that -for a soft spectrum like this -neutrinos around 100 TeV provide the largest contribution to the events. At this energy, the effective areas are almost independent of declination -whereas they are very different at higher energies, where the absorption of neutrinos in Earth becomes relevant, penalizing neutrinos coming from Northern hemisphere.

C. Extragalactic pp neutrinos (Xpp)
Neutrinos are created via pp (or Ap) interactions in astrophysical environments rich of gas, which serves as target for the interactions. In this case the accelerated primaries with a non-thermal spectrum collide with thermal protons, producing about equal amounts of π + , π 0 and π − leading to a neutrino spectrum adopting the spectral shape of the non-thermal primaries -which is typically a power law with a spectral index of about two.
Our prototype candidate class are Starburst Galaxies (SBGs), although this is not the only possible source class. A SBG is a galaxy that has an extremely high star formation rate of the order of 10 − 100 M /year, while normal galaxies have 1−5 M /year [53]. A SBG is full of young stars, which implies that supernova events are frequent with a rate of the order of 0.03 − 0.3 per year; moreover, it must contain a large amount of gas available to be injected into the star-forming processes. The simultaneous presence of cosmic rays, accelerated in shock waves created by the exploding core-collapse SNe, and targets, i.e., large amounts of gas, makes SBGs promising sources of cosmic neutrinos produced by pp interactions. SBGs are mainly observed in the infrared band because the bulk of their emission (UV radiation emitted by young stars) is absorbed and re-emitted in this band.
The SBGs have been proposed in several papers as source of IceCube neutrinos [24][25][26] although they cannot explain the entire IceCube signal [23,28]. There have been many other alternatives discussed in the literature, which potentially also fit this category, such as radio galaxies [54], AGN winds [55], and halo/galaxy mergers [56].
In our model, we assume an E −2 spectrum with an energy cutoff for neutrinos at 1 PeV [24]: This choice of the spectral shape for SBGs is motivated by the calorimetric limit, implying that the secondary spectrum follows the primary spectrum because the interactions are the dominant process and the particle spectrum inside the SBG is not modified by energydependent escape processes in that energy range. The cutoff energy is motivated by the knee of the cosmic ray spectrum, which may be somewhat higher than in our galaxy because of larger magnetic fields; note, however, that the accelerators themselves have to produce particles with high enough energies. More extreme scenarios are e.g. discussed in Refs. [26,57]. Our proposed spectrum is also compatible with the extragalactic diffuse gamma-ray background measured by Fermi [28], see Ref. [23], applied to the non-blazar contribution. 5

D. Extragalactic pγ neutrinos (Xpγ )
Neutrinos are produced by pγ (or Aγ) interactions in astrophysical environments with high radiation densities. The accelerated primary protons (with a non-thermal spectrum) collide with photons to produce pions. In such interactions, the spectrum of the neutrinos depends on both the spectra of the primaries and the target photons.
We choose TDEs as our baseline source, since we have a model available [35] which can not only describe the PeV neutrinos, but also the UHECRs in a self-consistent way. Other candidate classes include Active Galactic Nuclei (see e.g. Refs. [58][59][60][61][62]) and (low luminosity) Gamma-Ray Bursts (see e.g. Refs. [63][64][65][66][67][68]), which however require more study to draw a fully self-consistent picture. There are ways to discriminate between TDEs and other source classes, e.g., by neutrino point source searches, which we will discuss below. Note again that the UHECR connection requires nuclei in the sources, which means that our TDE example is actually based on Aγ interactions.
Tidal disruption is the process by which a star is torn apart by the strong gravitational force of a nearby massive or supermassive black hole. About half of the stars debris remains bound to the black hole, and is ultimately accreted. TDEs as the sources of extragalactic neutrinos [30][31][32][33][34][35][36] have been recently very actively discussed in the literature. In our model we assume that the shape of the neutrino spectrum follows the best fit spectrum of [35], which we refer to as (dφ BF TDE )/(dE ν ), with a free normalization: Note that the parameter F X−pγ can be directly related to the normalization G defined in [35] Here ξ A is the baryonic loading (defined as the energy injected as nuclei over the total X-ray energy in the Swift energy range 0.4-13.5 keV) andR(0) is the local apparent rate of jetted TDEs. The reference value chosen forR(0) is the rate of white darf-intermediate mass black hole disruptions inferred from observations, Thus, the result of our combined fit can be immediately interpreted in terms of the TDE baryonic loading × local apparent rate.

IV. RESULTS
In this section we present our results, including the best-fit model and the result of our energy reconstruction spectrum for pp interactions, the observed extragalactic diffuse gamma-ray background (EGRB) constrains the maximally allowed neutrino flux. Since most of the EGRB is expected to come from AGN blazars, the constraint on the non-blazar contribution is even somewhat stronger [28]. Our obtained SBG contribution will satisfy these limits, see Ref. [23].   technique. We also demonstrate how to determine the origin of individual events at the probabilistic level.
A. Best-fit model We perform a maximum likelihood analysis using a Poissonian χ 2 with the different components introduced in Section III and the free parameters F µ atmo , F e atmo , F gal , F X-pp , and F X-pγ . The data points for the different tested data sets are reconstructed in energy as described in Appendix A, and the theoretical rates are computed with the effective areas similar to Eq. (A5).
Let us first use the shower (HESE) events and the track (HESE + TGM) events separately. We obtain the best-fit parameters listed in Tab. II for the spectral fit in this case. This set of parameters is in agreement with the present constraints: The Galactic flux does not overshoot the theoretical and the experimental expectations, and the extragalactic flux from pp interactions is compatible with the diffuse gamma-ray background. Moreover, F µ atmo and F e atmo reproduce the total atmospheric background expected by IceCube. However, they do not reproduce the expected number of tracks and showers separately. A mis-identification of tracks as showers at the level of 30%-40% could reconcile the observations and the expectations. Note again that this problem related to the background is also present in the IceCube analysis itself, since in [8] about 40 tracks are expected (from signal, conventional atmospheric neutrino background and atmospheric muons) -instead of the 22 tracks detected.
Most importantly, we find in Tab. II an almost equal normalization for the Galactic and X pp contributions, independent of the sample (tracks or showers used). Note that this result has been obtained based upon the spectral distribution of the observed events, not to the angular distribution. The evidence for X pγ > 0 depends on the sample used: whereas the best-fit is zero for showers, it is similar to the TDE best-fit in Ref. [35] for tracks -which contains the 4.5 PeV deposited energy TGM. We use the value obtained from the tracks in the following, as we need to describe both showers and tracks. The combined fit (showers + tracks) yields, instead, X pγ 0.37, which is non-vanishing but a factor of two smaller than the one suggested by the track sample only. The χ 2 /d.o.f. are reported in Tab. II as well; a value of around one means that that the model describes data very well, while at the same time it is not "over-fitting" data.
Our main result is shown in Fig. 2, where we show showers (upper panel) and tracks (lower panel) separately, including the contributions from the individual components. We use the track best-fit here since it is in good agreement with the shower one at low energy, and we expect that it is more representative for our model at high energy since the statistics is much greater above several hundreds of TeV than for the shower sample. Note that the reconstructed data points are consistent with the IceCube results, and that our figure contains the residual atmospheric background both in data and fluxes explicitly (see discussion below). From the figure it is clear that our model describes both tracks and showers extremely well. At low energies, there are no obvious contradictions with the residual atmospheric backgrounds. At high energies, the possible evidence for X pγ > 0 emerges from the TGM dataset (pink), especially the 4.5 PeV track -which is likely to be generated by an incoming neutrino above 10 PeV. As possible alternative would be to extend X pp to higher energies, but this would throw up the question if SBG can accelerate primaries to energies significantly beyond 100 PeV.
We remark again that we show shower and track events separately since they are affected by different backgrounds. Particularly the HESE track sample is more affected by the atmospheric background, since not only atmospheric neutrinos, but also atmospheric muons contribute to the background. For this reason the HESE shower dataset is more representative of the astrophysical signal at low energies, whereas the TGM dataset is more representative of the astrophysical signal at high energies due to the larger effective area related to TGM compared to HESE.
One of the questions, we need to address, is how many components we actually need. We can do that at the statistical level by removing the components one-by-one and re-performing the fit. The atmospheric backgrounds are inevitable as they are needed to the describe the spectrum below 100 TeV; this is also statistically evident. Removing the Galactic component, the quality of the spectral fit is only marginally affected. However, we have to take into account that the Galactic compo- nent also plays a role to explain the angular distribution of HESE events above 100 TeV and the spectral difference between HESE and through-going muons. The presence of a Galactic component, mainly observed from the Southern hemisphere 6 , helps to reduce the difference between the observations from South and North. Our Galactic component contributes 25% to the total signal and 30%-35% of the astrophysical events coming from South.
The X pp component provides the largest contribution to the astrophysical signal, which is about 50%. The flat E −2 spectrum is necessary to explain the middle range of energy between 100 TeV and 1 PeV. Removing this component, the quality of the combined fit goes from χ 2 min = 14.4 to 27.0, implying a 3.5σ evidence for the X pp contribution (even in the presence of X pγ and after re-marginalization). Moreover an extragalactic pp spectrum, generated in sources such as SBGs, can explain the lack of correlations between neutrinos and known sources of γ-rays, since SBGs are rarely observed in γ-rays because the photons are reprocessed inside the sources. SBGs are also compatible with the point source bounds, as they are faint enough to avoid multiplets from the same source [21][22][23].
The evidence for the X pγ contribution is much weaker from the purely statistical perspective only. 7 If one used the fit from the track sample, the pγ component would be responsible for about 5-6 events, which are the most energetic ones. Future measurements and particularly 6 70% of the Galaxy is observed from the Southern sky, whereas only 30% of it can is observed from the Northern hemisphere (see Fig. 2 of [12]). 7 Here χ 2 min = 4.7 goes to 5.7 using the through-going muon dataset, whereas it remains almost the same using both HESE and TGM, since there is a slight tension between HESE (that disfavor the pγ component) and TGM (that favor the pγ component).
At least 2 events At least 3 events IceCube-Gen2, with a 6-8 times greater exposure [69] can clarify the situation, confirming or excluding the presence of a pγ component at PeV energies. If interpreted as jetted TDE contribution, we find about 80% of the flux normalization in Ref. [35], which can be interpreted e.g. in terms of a lower baryonic loading 430, cf. Eq. (8).
Since TDEs are potentially very luminous neutrino sources, one would eventually expect several events from the same source. We compute the probability to observe at least two or three events from at least one TDE explicitly in Fig. 4 as a function of the baryonic loading, assuming that the typical TDE flux corresponds to the one in [34] (Fig. 2, upper-left panel). 8 From this figure, it is evident that the probability to observe two events from one TDE is less than 50% in the current sample, and the probability for a triplet is at the percent level. Note that the limit for individual TDEs, such as Swift J1644+57 [70], can be stronger as the effective area depends on declination. However, this method will be powerful in the future to identify the origin of the X pγ contribution and to discriminate TDEs from other classes of sources such as LL-GRBs and AGNs. Other method include the test of a possible flavor transition at the highest energies, which is expected for TDEs [34] and GRBs (e.g. Ref. [71]) due to magnetic field effects on the secondary pions, mesons, and kaons, but not for AGNs.
In Fig. 3, we compare our model to the data points of the IceCube analysis in which the background has been already subtracted [8]. We notice a good agreement between our model and the IceCube data points in the energy range between 300 TeV and several PeV, and between our reconstructed points (see pink points of Fig. 2) and the IceCube fit of the through-going muons (yellow band of Fig. 3). On the other hand, at low energies, a certain discrepancy is present: The observation of the excess at 150 TeV by IceCube could be related to an underestimation of the background, since in their data points the background is already subtracted. The gap around 600 TeV only present in the HESE shower dataset in our analysis, whereas it is neither present in our reconstructed track sample, nor in the IceCube TGM analysis. Therefore it could be simply produced by a statistical fluctuation or it could be an hint of the change of the flavor composition at very high energy, although there are not enough data to confirm the last hypothesis.
A special remark concerns the upper limit close to 7 PeV, which is in the Glashow resonance [72] region. It is plausible that the reconstruction of this data point is strongly affected by the production mechanism. In the IceCube analysis, usually φ νe = φν e is assumed; see for example Fig. 7 of [1]. This assumption applies to pp interactions -whereas the flux of electron antineutrinos can be suppressed for the pγ production of neutrinos [73][74][75][76], especially in the damped muon regime, in which the electron antineutrino flux at detection is strongly suppressed. A comparison to the track sample (lower panel in Fig. 2) is therefore safer here. We have also checked that the number of Glashow resonant events in the most optimistic case (pp interactions) is roughly 0.5 after 5.7 years of exposure using the best-fit parameters of the model and the approximation for the effective area from Ref. [42], which means that at present, our model is compatible with the non-observation of Glashow events. A similar estimate can be performed for double pulse events [77] for E ν 0.5 PeV, i.e., the detection of two separate 8 Using the IceCube HESE effective areas (corresponding to the average TDE at random position), we find that the expected number of events coming from a single TDE N 0.24 × ξ A /540. Using our multi-component model and the present IceCube exposure, i.e., 5.7 years for HESE, we find that about six events can be attributed to TDEs. This means that about 25 TDEs contribute to the flux. From this information, Fig. 4 can be explicitly computed following Ref. [20] using a Monte Carlo method. interaction vertices for ν τ . Using the effective area in Ref. [75] and assuming that only ν τ are created at the source (unrealistic, but most conservative scenario), we expect about 0.5 double pulse events after 5.7 years of exposure for our model. Again, our model is compatible with the non-observation of double pulses at present.

B. Event-by-event analysis
Using our multi-component model, we can compute the probability distribution to each observed HESE or TGM event on its possible origin. This procedure can be sketched as follows (for details, see below): • Determine the distribution of the reconstructed energy of each event as described in Appendix A; • Determine the probability to belong to each component of the multi-component model (best-fit) according to the spectrum, performing a Monte Carlo extraction of the reconstructed energy; • Determine the probability to belong to the Galactic contribution according to direction and directional resolution.
In order to include the direction of each event in our calculation, we need a function related to the compatibility of a certain direction with the Galactic plane, since the Galactic flux is present in the inner Galaxy but absent outside -including the directional resolution of the events. Note that the Galactic flux actually dominates in the inner Galaxy region up to energies of about 2 PeV, as it is illustrated in Fig. 5, whereas it is absent outside this region. Since we have assumed that the Galactic flux is present in the region |ϕ| ≤ 5 • and |λ| ≤ 45 • , we define the directional probability density function such that it is one in that region and drops exponentially outside with the directional resolution: where G is a Gaussian function similar to Eq. (A1) (but normalized in the range −∞ to +∞) with mean value m h and standard deviation σ h , and Θ is the Heaviside step function. We use two functions for Galactic latitude and longitude, i.e., h corresponds to ϕ or λ; the product of the previous two functions give the compatibility C gal ≡ C ϕ gal (ϕ) × C λ gal (λ) of each event with a Galactic origin, based on its direction. Note that we have m ϕ = 5 • and m λ = 45 • by construction, whereas the standard deviations are different for showers and tracks, for which we use average values σ ϕ = σ λ 15 • and σ ϕ = σ λ 1 • , respectively (see e.g. Ref. [1]). The physical meaning of this function is that an event inside the inner Galaxy is fully compatible with a Galactic origin (C gal = 1), whereas the probability drops exponentially to C gal = 0 going away from the inner Galaxy. To have the overall probability we have to include the information coming from the energy spectrum.
Therefore the overall probability to belong to a certain component is given by based on the spectrum at the reconstructed energy E ν . Since the reconstructed energy is not unique, we perform a Monte Carlo simulation (integration) to determine the average probability.
Using this procedure and integrating over E ν in Eq. (10), we can derive the probabilities that an event belongs to the different component. Since the list of events includes the atmospheric backgrounds, it should be clear now why we have carried the residual atmospheric background through the whole procedure. We give the individual probabilities for each event to belong to atmospheric background, Galactic contribution, X pp and X pγ in Appendix B. We find for the HESE showers very high probabilities to belong to the atmospheric background in most cases, whereas in eight cases the probability to belong to the Galactic component is higher than 50%. In 25 cases, the probability to belong to X pp is higher than 30% (up to around 50%), whereas only two events exhibit some evidence to belong to X pγ (probabilities 34% and 40%). For the HESE tracks, most events seem to belong to the atmospheric background with five exceptions with a significant probability (larger than 30%) to belong to X pp . The TGM sample, on the other hand, seems to be dominated by X pp , whereas one event may be of Galactic origin and 12 events have a probability higher than 30% to belong to X pγ -the 4.5 PeV muon track belong to that component with a 83% probability. The TGM sample is therefore the sample least dominated by the background and maybe best suited to study the extragalactic neutrino fluxes, whereas the Galactic contribution hides in the HESE showers.
This result can be also pictorially represented, see Fig. 6, where we show three different samples of the multi-component distribution produced with our Monte Carlo method based on the probability distribution for each event. The different colors corresponds to the different components as described in the figure caption, and circles represent showers and triangles represent tracks. The inner Galaxy region including the corresponding angular uncertainty for showers is also shown (dark orange lines). We notice that the events close to the Galactic plane are likely to be associated to Galactic neutrinos in all the maps (orange), whereas the assignment to the component changes in many cases. Some of the events, however, have a high probability to be of atmospheric origin (red). The extragalactic flux is dominated by X pp (green), but the assignment to this component can only be performed on a statistical basis. The X pγ contribution (blue) tends to apply to the highest energetic events, about 2-3 events in the shown representatives.
We furthermore show in Fig. 7 the unfolded IceCube energy spectrum (points with error bars) for the three samples (HESE showers, HESE tracks and TGM) separated into the different components. Here the method described in Appendix A have been used, separating the event rate contributions of the different components. Let us remark that the data points and error bars are obtained using the probability Eq. (10), which means that they depend on our model, whereas the points reported in Fig. 2 are model independent. We can use the figure to see what datasets at what energies describe the different components. For example, one can see that the residual atmospheric backgrounds are well determined at low energies from HESE data. The HESE showers constrain the Galactic contribution well around 100 TeV, whereas the X pp contribution is equally well contributed to by HESE showers, tracks, and TGMs. The evidence for X pγ is mostly driven by the TGM sample.
Of course, our analysis is only useful on a statistical basis, but it can be easily extended to future datasets. It allows one to use the information about an individual component if a theoretical model for one source class is studied. There are some limitations of our method. First of all, we do not take into account individual directional uncertainties for the events, which could be used to make the analysis more precise. In order to obtain the uncertainties in Galactic coordinates, a Monte Carlo extraction would be required since the direction and the uncertainties of each event are given in equatorial coordinates and the transformation is not trivial. For the sake of simplicity, we use in Eq. (9) the average values of the angular resolution 1 • and 15 • for tracks and showers, respectively.
Second, the effective area is actually declinationdependent, whereas we use averaged effective areas. This describes the isotropic extragalactic component, and it is a good approximation for the Galactic component (as discussed in Section III). For the atmospheric flux, we fit the residual atmospheric background passing the vetos and analysis chain with a free normalization independent for tracks and showers. While the effective area for individual events may exhibit a non-trivial declinationdependence, this dependence enters our free normalization as an effective average as long as the residual backgrounds roughly follow the atmospheric flux shape. The deviation between track event rate and residual background at low energies in Fig. 2, lower panel, may be an artifact such an energy-dependent effect.

V. SUMMARY AND CONCLUSIONS
In this work we have proposed a multi-component model to address the following challenges in current data of astrophysical neutrinos: 1. The observed through-going muon spectrum, coming from the Northern hemisphere, is considerably harder than the HESE sample.
2. The HESE dataset exhibits an anisotropy if only events above 100 TeV are considered, indicating a correlation with the Galactic plane.
3. No correlations with known sources have been observed; no point sources have been identified.
4. The limit from the observed extragalactic gammaray background must be obeyed.
5. The recent observation of a muon track with 4.5 PeV deposited energy requires a primary energy 100 PeV.
We have started from the event topology and deposited energy, and we have computed the probability distribution in terms of event type and reconstructed neutrino energy event-by-event. Consequently, we have reproduced the reconstructed energy distributions by IceCube very well. In comparison to the IceCube publications, we have included the residual atmospheric background explicitly, because we have (at the end) listed the probability for each individual event to belong to each of these different components. Using this approach, we observe two differences compared to the IceCube reconstruction: we do not observe an excess at about 150 TeV, and the gap at about 600 TeV is only present in the shower, but not in the track samples, which points towards a statistical fluctuation (compare data points in Fig. 3 to Fig. 2). Our model consists of four populations summarized in Tab. I: (residual) atmospheric contribution (including prompt contribution) passing the vetos, Galactic contribution following the cosmic-ray flux, extragalactic flux from pp interactions (X pp ), such as starburst galaxies, and extragalactic flux from Aγ or pγ interactions (X pγ ), such as Tidal Disruption Events or AGN blazars. We have taken the spectra shape for each of these components from plausible examples from the literature and we have fit the (independent) normalizations. We have demonstrated that the above challenges are solved in this model in the following way: • The through-going muon (TGM) sample is dominated by X pp and X pγ ; as a consequence, the measured energy spectrum is close to E −2 ν in the energy range between 200 TeV and 10 PeV. On the other hand the HESE dataset is more sensitive to low energy events, that are much more affected by atmospheric backgrounds and by the Galactic component, especially for events coming from the Southern hemisphere.
• Above 100 TeV, the atmospheric background is highly suppressed in the shower HESE dataset, whereas the Galactic component is visible and dominates over the other fluxes in the inner Galaxy region, justifying the 2σ excess close to the Galactic plane.
• The dominant astrophysical X pp contribution, such as from starburst galaxies, implies abundant and low luminous sources which are in consistency with current stacking or point source searches. Because of the relatively hard spectrum, the bounds from the extragalactic diffuse gamma-ray background can be easily avoided.  FIG. 7: Best-fit model (curves, for track best-fit in Tab. II) and unfolded IceCube energy spectrum (points with error bars) for HESE showers, HESE tracks and through-going muons using our method described in Appendix A. Here the data points and error bars are shown for the individual contributions given by the different components, i.e., they depend on the event splitting as given by the fit of the multi-component model.
• The observed 4.5 PeV track-like event, that is likely to be produced by a ∼ 10 PeV muon neutrino, justifies X pγ . As a baseline model, we have used jetted Tidal Disruption Events which have recently been demonstrated to be compatible with the origin of UHECRs. Since the TDE events are rare but luminous, they can be tested by point source and multiplet limits in the future to discriminate this hypothesis from alternative scenarios such as low luminosity AGN blazars or GRBs.
We note that the statistical evidence for X pp is high (3.5σ, even if X pγ is present), while that for X pγ is not significant yet from the purely statistical perspective. The role of X pγ could be taken over by the X pp component if that was able to accelerate primaries to energies larger than about 100 PeV. While theoretical arguments seem to disfavor this for starburst galaxies, there are recent hints for correlations among UHECRs and starburst galaxies at the observational level by the Auger experiment [78]. An increase of exposure, such as by IceCube-Gen2, will help to test it: since about 6-7 HESE per year are expected from X pγ at PeV energies, the non-observation of these neutrinos would rule out X pγ at 5σ after about 3 years.
We conclude that a self-consistent description of the observed diffuse neutrino flux requires multiple contributions to the diffuse astrophysical neutrino flux, and we have drawn a self-consistent picture describing these contributions by their generic characteristics in terms of spectrum, sky distribution and neutrino production mechanism. Our model solves the key challenges in current neutrino data, including multi-messenger constraints. In the future, it may be helpful to obtain neutrino constraints on the individual contributions in addition to the total spectrum to disentangle the contributions from different source classes. Our model will also help theorists to draw a self-consistent picture if only one of the proposed generic contributions is considered.

Description of the interactions
Here we define the probability density functions the connect the incident neutrino energy with the deposited energy.
ν e charged current (CC) interactions. If electron neutrinos (antineutrinos) interact via a charge current interaction, almost the entire energy is deposited into the detector as an electromagnetic shower. Assuming that the fraction x ≡ E dep /E ν of the incident neutrino energy E ν is deposited in the detector, i.e., the deposited energy E dep = x E ν with 0 ≤ x ≤ 1, we define a probability distribution function PDF e sh (x) proportional to a Gaussian For ν e CC interactions, we use m e,CC = 0.95, implying that most of the energy is deposited in the detector, and σ e,CC = 0.06 (see Tab. 1 and Fig. 1 of [79]). Our assumption can be compared to the IceCube results comparing the left panel of Fig. 8 with the left panel of Fig. 1 in [79]. The probability that a neutrino is a ν e and interacts via CC interaction producing a shower-like event is given by evaluated at the incident neutrino energy, whereσ is the cross section. The ratioσ CC /(σ CC +σ NC ) is in good approximation constant in the energy range between 100 TeV and 10 PeV, i.e., our range of interest, and it is equal to 0.71 [80]. Moreover we assume the equipartition of flavors in order to simplify the calculation, which leads to P e sh 1/3 × 0.71 23.7%. 9 ν τ charged current (CC) interactions. If tau neutrinos (antineutrinos) interact via charged currents, a fraction of energy is carried away by neutrinos generated by τ decays which is not detectable. Considering the various channels in which a τ can decay, the average quantity of energy fraction detectable is about 80% [81]. We also use a Gaussian, see Eq. (A1), with m τ,CC = 1 and σ τ,CC = 0.25 in order to reproduce the average value of the transferred energy. Therefore the uncertainty is larger here than for ν e CC interactions, see middle panel of Fig. 8.
Moreover ν τ can produce muons with a branching ratio of B.R.(τ → µ) = 17.4%. This contribution has to be subtracted, since we only discuss shower-like events in this section. The probability that a neutrino is a ν τ and interacts via charged currents producing a shower-like event is given by: which is P τ sh 19.5% Neutral current (NC) interactions. For NC interactions, which are not sensitive to flavor, a large fraction of the incoming neutrino energy is not visible because it is carried away by outgoing neutrinos. The average fraction of deposited energy is 25%, but large fluctuation are expected in this case (see Fig. 1, right panel, in [79]). We therefore use m NC = 0 (events that deposit a small fraction of energy are favored) and standard deviation σ NC = 0.3 in Eq. (A1). This choice reproduces the average fraction of deposited energy, and it also describes in good approximation the right panel of Fig. 1 in [79] (compare to the right panel of Fig. 8).
The probability that a neutrino interacts via neutral currents is given by Track-like events: ν µ . Track-like events are generated by ν µ that interact via CC interaction and produce a muon, which is visible as track in the detector. For track-like events the information on the incoming neutrino energy is very uncertain and the deposited energy only represents a lower limit to the incoming neutrino energy (the track often originates outside the detector). As for NC interactions the deposited energy is usually a small fraction of the incoming neutrino energy; therefore we use the same probability density function used for NC interaction, i.e., PDF µ tr = PDF NC sh . We have checked that this approximation reproduces the measured throughgoing muon spectrum in good approximation, as it can be seen comparing the pink points of Fig. 2 to the yellow band in Fig. 3. The probability that a neutrino is a ν µ and produces a track-like event is given by Track-like events: ν τ . Track-like events can be generated, in smaller number, by muons created in the ν τ CC interactions. The probability that a neutrinos is a ν τ which interacts via CC interaction and produces a track-like event is given by: Re-call that here the muon is generated by tau decay, which means that the muon takes about 1/3 of the ν τ energy. We therefore use PDF τ tr (x) = PDF µ tr (3x).

From deposited to reconstructed energy
We know construct the probability density function R (E dep , E ν ) (differential in E ν ) of the reconstructed neutrino energy for a given deposited energy and event type . Note that for the sake of simplicity we use the symbol E ν simultaneously for the incident (true) and reconstructed neutrino energies.
Compared to the PDFs in the previous section, the function R will depend on the spectral shape. We assume an E −2 flux for the neutrino spectrum by default, but we will discuss the consequences of different assumptions below. Because the cross section of neutrinonucleons interaction scales ∝ √ E in the energy range between 100 TeV and 10 PeV [80], the function R is given by where the last factor is the product between spectrum and cross section energy dependencies. Here denotes one of the process previously described, namely showers given by ν e CC, ν τ CC or NC interactions and tracks given by ν µ CC and ν τ CC interactions. The factor E dep /E 2 ν comes from the change of variables from x to E ν in the normalization Eq. (A2). In order to derive this formula, see also Sec. 2 of Ref. [82]. 10 The probability that an event with a certain deposited energy has been generated by a specific process becomes We have checked that these probabilities, to a very good approximation, do not depend upon the deposited energy. The probabilities that a shower-like event has been produced by ν e CC, ν τ CC, and NC interactions are 55%, 35%, 10%, respectively. The probabilities that a track-like event has been produced by ν µ CC or ν τ CC interactions are 96% and 4%, respectively. Note thatP includes the spectral re-distribution of events, as compared to P , and that we have pulled out P from the integrals because they hardly depend on energy. For example,P NC sh is relatively (compared to the other charged current showers) smaller than P NC sh because a NC shower with a certain deposited energy tends to come from a higher energy than a CC shower, where the flux is lower.

Monte Carlo simulation and deconvolved neutrino flux
Starting from the list of events with event topology (track-like or shower-like) and E dep , we construct the probability distribution for the deposited energy C i (E ν ) for the ith event as follows: 1. We determine a process that has generated the observed event from the probabilityP in Eq. (A4).
2. We extract an energy associated to the process according to the distribution R (E dep , E ν ) in Eq. (A3).
3. We repeat this procedure 10 4 times to obtain the reconstructed energy distribution C i (E ν ) for this event.
We show in Fig. 9 the distribution of the reconstructed (neutrino) energy for a shower with deposited energy of 1 PeV (left panel) and a track with deposited energy of 4.5 PeV (right panel), i.e., the most energetic event observed by IceCube up to now. The reconstructed energy of a shower-like event is close to the deposited energy and there is a small uncertainty, whereas the reconstructed energy of a track-like event is two to three times larger than the deposited one, with a much bigger uncertainty. We have checked that both our probability distribution and our reconstructed energies are in very good agreement with the theoretical estimates in Tabs. IV and V of Ref. [39], which appeared during the completion of this work. The choice of a softer spectrum instead of the assumed E −2 would reduce the uncertainty on the reconstructed energy, privileging the events that deposit most of their energy. We have verified that using an E −2.5 spectrum, the uncertainty on the reconstructed energy is reduced by about 30% for shower-like events and by about 25% for the track-like events, deteriorating both the agreement with [39] and the TGM analysis performed by IceCube [4,8].
We can now determine the total number of events in each reconstructed energy bin j covering the energy interval [E j , E j+1 ], where we choice of the binning is arbitrary. The distribution function for all events of one topology is given by where i runs over all events from a dataset, and the num-ber of events assigned to the jth bin is We assume that the uncertainty on N j ν is N j ν when N j is large; when N is small, namely smaller than one, we set a 90% C.L. limit to the expected number of N j ν using Poissonian statistics.
Thanks to the effective areas (HESE [1] and throughgoing muons [83]) it is possible to convert the N i ν into the expected flux in each energy interval, using the following formula for shower-like HESE and the following one for track-like HESE or throughgoing muons (they differ for the effective area A µ eff ): whereĒ = 10 (log 10 (Ej )+log 10 (Ej+1))/2 is the middle point on a log-scale, and η = 0.8 represents the fraction of the ν µ effective area that is connected to track-like events (see [43]). The effective areas are taken from [1] for HESE and from [83] for TGM. The reconstructed data points are, for example, shown in Fig. 2. Here the uncertainties on the flux are proportional to the uncertainties on N j ν as described above, whereas the uncertainties on the energy axis cover the interval between E j and E j+1 .
In the previous calculation we have assumed an E −2 spectrum in each energy interval between E j and E j+1 . This approximation is not accurate when we analyze low energies (especially for track-like events), in which the atmospheric neutrino spectrum E −3.7 dominates. A possible solution is to choose small energy intervals to reduce this effect. With our assumption the background related to φ µ might be underestimated in the data points at very low energies.