Low energy IceCube data and a possible Dark Matter related excess

In this Letter we focus our attention on the IceCube events in the energy range between 60 and 100 TeV, which show an order 2-sigma excess with respect to a power-law with spectral index 2. We analyze the possible origin of such an excess by comparing the distribution of the arrival directions of IceCube events with the angular distributions of simply distributed astrophysical galactic/extragalactic sources, as well as with the expected flux coming from DM interactions (decay and annihilation) for different DM profiles. The statistical analysis performed seems to disfavor the correlation with the galactic plane, whereas rules out the DM annihilation scenario only in case of small clumpiness effect. The small statistics till now collected does not allow to scrutinize the cases of astrophysical isotropic distribution and DM decay scenarios. For this reason we perform a forecast analysis in order to stress the role of future Neutrino Telescopes.


Introduction
The observation of astrophysical neutrinos made by the Ice-Cube experiment (IC) [1] has been an important step in the field of Neutrino Astronomy, whose impact on both Particleand Astro-Physics has still to be unveiled. Recently [2] the IC Collaboration has delivered in a preliminary work the results of four years data. In Figure 1 (upper panel) we report the so-called excess in the number of events as a function of the neutrino energy, which is the number of events detected in IC once one has subtracted the background (atmospheric conventional and prompt neutrinos plus muons [2]), and an astrophysical component characterized by a E −γ ν power-law, with γ = 2 being considered as benchmark prediction [1,2,3]. The plot seems to suggest the presence of an excess of events in the energy range 60 -100 TeV that has some tension, order 2σ, with the simple E −2 ν contribution. As can be seen from Ref. [2], such an excess partially disappears if one considers a steeper astrophysical component with power-law E −2.58 ν , even though such an exponent for a neutrino flux can be hardly explained.
In general, the cosmic ray spectrum is characterized by a power-law [4], which is usually understood in terms of acceleration from shock fronts, the so-called Fermi mechanism [4,5]. An astrophysical neutrino flux is expected to be produced during the hadronic matter acceleration and interaction with radiation (pγ) or with gas (pp). Thus its dependence on energy should be, at the source, mostly related to the differential spectrum of charged cosmic rays and to the pions production efficiency. If needed, one should take into account propagation effects so that, for example, the galactic cosmic ray spectrum at the Earth becomes ∝ E −(γ+δ) (γ + δ ≈ 2.7) up to the ankle at 1 ∼ 10 PeV [6]. The quantity δ depends on galactic magnetic fields and is evaluated through cosmic rays secondary to primary ratio measurements [4,7]. On the other hand, galactic astrophysical neutrinos, which are not affected by magnetic fields, have a flux ∝ E −γ ν , with γ ≈ 2. Even considering extragalactic sources to be significant (a reasonable hypothesis given the high galactic latitudes of some events) it is not easy to justify a steep flux [8]. Models with pγ interactions produce peaked spectra [9], so one could have a steep flux depending on the position of the peak. However, in this case one expects a peak in the spectra in the region of several PeV like for Active Galactic Nuclei [10]. Another kind of pγ source, Gamma Ray Burst, provides a flux with an upper limit (given by searches for correlation with observed GRB) more than one order of magnitude below the observed flux [11]. Moreover, hard and smooth spectra are expected in pp scenarios [12], like for extragalactic Supernova Remnants [13]. More generally, theoretical models of acceleration mechanism for hadronic matter produce a flux that should be at most as soft as Independently of the detailed theoretical model, we have, most importantly, observational constraints coming from multimessenger approaches like [16], which combine IC data with Fermi γ-ray measurements, and give strong bounds to a hadronuclear origin (while leave room for hidden pγ sources [17]). In fact, a spectrum as soft as the one obtained fitting IC data implies, when assuming sources transparent to radiation with respect to two-photon annihilation, an expected γ-ray flux bigger than the measured one [16,17]. The previous considerations support the assumption of an astrophysical E −2 ν power-law used to obtain the excess reported in Figure 1. Starting from the excess in the number of events one can obtain the corresponding quantity for the whole neutrino flux (summed on all flavors) once the effective area of the detector is taken into account [1]. In Figure 1 (lower panel) we report such a flux as a function of the neutrino energy. As already discussed in Ref. [18], we assume for simplicity the equality between the deposited and neutrino energy due to low statistics at our disposal. At the energy scale O(100) TeV, this is not strictly true for neutral current interactions [19]. When a significant statistics is collected, the average ratio between the two energies, which is of the order of (97% σ CC + 23% σ NC )/(σ CC + σ NC ) ∼ 75%, could be applied.
In this Letter we assume that the above excess, mainly concentrated in the energy range 60 -100 TeV, has a genuine physical origin. Under this ansatz, it is worth pursuing, for this energy bin, a study in order to unveil the nature of such an excess. We perform our analysis assuming as null hypothesis one of the following alternatives for the source of the IC data: i) astrophysical, which can be investigated by studying, in first approximation, the correlation with the galactic plane or with an isotropic distribution for galactic or extragalactic astrophys-ical sources, respectively; ii) induced by Dark Matter via decay or annihilation, hence related to the first or second power of the particular Dark Matter (DM) density profile adopted. Moreover, even though the small number of events already detected does not allow to exclude all DM scenarios, one can perform a forecast analysis in order to determine the required statistics.
In order to compare the IC observations with possible DM predictions we consider both decaying and stable Dark Matter cases. In the first case, 60 − 100 TeV neutrinos detected at Ice-Cube would be originated directly from the decay of the DM particles, while for a stable DM particle neutrinos are only produced via annihilation. In both cases the resulting neutrino flux would be composed by both a galactic and an extragalactic DM component. Different approaches proposed in literature [18,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34] have studied the possible presence of DM hints in the PeV range, namely for the most energetic IC events. Here we take an alternative point of view, assuming that PeV events have bottom-up origin and considering for lower energy data a possible top-down origin due to DM particles with mass scale O(100) TeV. Independently of the mass scale and of the DM couplings, neutrinos originated from DM would have an angular distribution that is more peaked around the Galactic Center where a higher DM density is expected. This is true in particular when assuming an annihilating DM, because of the squared enhancement factor. Of course, this effect is dependent on the assumed DM galactic halo profiles; for example, one could take the Navarro-Frenk-White profile (NFW) [35] or different distributions like the Isothermal profile (Isoth.), which implies a more isotropic flux.

The analysis
In order to infer about the physical origin of the excess we compare the angular distribution of the observed events in the energy bin 60 -100 TeV (in the following we discuss the procedure to take into account the presence of the background and of the experimental errors) with the angular distributions of astrophysical galactic sources (galactic plane) and extragalactic ones (isotropic distribution), as well as with the expected flux coming from DM interactions (decay and annihilation). In this approach the astrophysical E −2 power-law contribution is regarded just as an additional term to the background events counting for atmospheric neutrinos and muons. For this reason hereafter we denote as background the sum of atmospheric neutrinos, muons and neutrinos coming from the astrophysical E −2 power-law. Moreover, due to the small number of events collected till now in the energy bin under study, in this analysis we take the simplicity assumption to consider just one additional component to neutrino background at a time (alternative scenarios i) or ii) of previous section) to explain the excess. This allows us to be more predictive even though more involved scenarios can be proposed where the excess in neutrino flux can be explained in terms of several components of different origin. Other analyses have already been presented in literature with different assumptions [29,36]. In particular, in Ref. [29] it has been studied the possibility that the whole neutrino spectrum has a DM origin, while Ref. [36] analyzed only the high energy neutrino spectrum (E ν > 150 TeV) considering also mixed components to the neutrino flux. Differently from previous analyses [29,36], we also take into account the angular efficiency of the IC detector for all neutrino flavors [44]. In Figure 2 it is reported the normalized IC effective area, averaged in the energy range considered (60 − 100 TeV), as function of sin δ, where δ is the declination in the equatorial coordinates system. In the region sin δ > 0 (North Pole) the effective area is mainly affected by the Earth absorption for all neutrino flavors, whereas for sin δ < 0 (South Pole) only the muon effective area is slightly dependent on the declination. In order to simplify the description, in the following we obtain the expected neutrino angular distributions for the analyzed cases without explicitly considering the IC effective area, which will be instead considered for the real analysis.
In case of astrophysical scenarios, namely galactic plane and isotropic distribution, the expected angular distributions in the arrival directions are Note that the galactic plane angular distribution depends on the Galactic latitude b only (l being the Galactic longitude), where b gal represents the angular size of the Galactic disk. It is reasonable to assume that the galactic neutrino angular distribution has the same characteristics of the galactic gamma-ray one [37].
Using the Fermi-LAT template [38], the quantity b gal has been varied in the range [2 • , 4 • ]. We assume for simplicity that the astrophysical sources are uniformly distributed in the Galactic disk. Indeed, it is worth observing that, however, such an approximation does not dramatically change the qualitative result of our analysis, which is mainly affected by the low statistics at our disposal.
In case of decaying DM scenario, the expected differential neutrino flux results to be a sum of the contribution due to the Galactic halo, and of an Extra Galactic term In the previous expressions M and τ stand for mass and lifetime of DM particles 2 , whereas ρ h (r) is the DM galactic halo profile, being r(s, cos θ) = √ s 2 + R 2 − 2sR cos θ with R 8.5 kpc (the Sun-Galactic Center distance) and cos θ ≡ cos b cos l. Moreover we take Ω DM = 0.2685, ρ c = 5.5 × 10 −6 GeVcm −3 and H(z) = H 0 Ω Λ + Ω m (1 + z) 3 with h ≡ H 0 /100 km s −1 Mpc −1 , Ω Λ = 0.6825 and Ω m = 0.3175 [39]. The quantity f (E) ≡ dN/dE is the neutrino energy spectrum produced by the decay of a DM particle, whose expression depends on the specific model considered. Since the excess is localized in the energy range 60 -100 TeV, it is reasonable to assume that the neutrino energy spectrum is almost peaked in the same range. This means to assume a negligible neutrino energy spectrum for energy larger than 100 TeV. Note that a not vanishing neutrino energy spectrum for energy smaller that 60 TeV does not provide any contribution neither to the galactic component nor to the extragalactic one for the energy bin considered. Therefore, independently of the shape of the energy distribution f (E), the redshift integral in Eq. (5) has an upper limit that is at most equal to z max = 100 TeV/60 TeV − 1. This is determined by the relation E = E (1 + z) when the conditions E ≤ 100 TeV and 60 TeV ≤ E ≤ 100 TeV are applied. Therefore, by integrating Eq. (5) in the 60 -100 TeV range one obtains dJ EG dΩ where E m = 60 TeV and E M = 100 TeV. There exist two interesting limits: i) the neutrino energy spectrum is fully contained in the energy bin considered, then f (E m ) f (E M ) 0; ii) the neutrino energy spectrum is wider than the energy bin, hence the energy distribution f (E) can be considered almost flat within it, namely f (E m ) f (E M ). This implies that the angular distribution takes the form with where α = 1 and 0 in the first and second case respectively, so providing β 1 = 0.43/H 0 and β 0 = 0.56/H 0 . In the following, to be more conservative we consider only the case corresponding to a larger extragalactic contribution (α = 0), since in this case the isotropic cosmological contribution results to be more competitive with the Galactic term. It is worth observing that including explicitly the IC effective area in Eq. (6) does not change this result, and the case with α = 0 still represents the most conservative scenario. Similarly, one obtains the following angular distribution in case of annihilating DM scenario where ∆ 2 0 is the clumpiness factor [40] (see also Ref. [41]). The quantity ∆ 2 0 ranges from 10 4 to 10 8 [42] depending on the model considered. In our analysis we consider three particular cases where ∆ 2 0 is equal to 10 4 , 10 6 and 10 8 corresponding to an extragalactic contribution that is sub-dominant, comparable and dominant with respect to the galactic one, respectively. However, recent studies like [43] state that the cumpliness factor ∆ 2 0 can be as large as few times 10 6 , considering unphysical larger values for such a quantity.
In our analysis, we consider two different DM galactic halo profiles [41]: the Navarro-Frenk-White distribution where r c 20 kpc and ρ h = 0.33 GeVcm −3 , and the Isothermal distribution ρ Isoth where r c 4.38 kpc and ρ h = 1.39 GeVcm −3 . Since in each case the distributions depend on one angle only, we can perform a one-dimensional statistical test. In particular, we use two different non-parametric statistical tests: the Kolmogorov-Smirnov test (KS) [45] and the Anderson-Darling test (AD) [46]. These statistical tests make a comparison between the cumulative distribution function (CDF) of the null hypothesis distribution function and the empirical cumulative distribution function (EDF), given by where n is the number of observed events cos θ i . Note that, in case of galactic plane angular distribution, the variable cos θ has to be changed into sin b. In the Kolmogorov-Smirnov test, the Test Statistics (TS) is the maximum distance between the previous two cumulative distribution functions and it is defined as whereas in the Anderson-Darling test the Test Statistics is given by In particular, this expression is very sensitive to the difference between the functions EDF and CDF at the two endpoints, suggesting that the Anderson-Darling test is a suitable test for our analysis (note that the Galactic Center is in correspondence of cos θ = 1).
To take into account the atmospheric background, we consider all possible different choices of 5 background events among 12, namely 12!/(5! 7!) = 792 combinations. Moreover, we include in our analysis the angular uncertainty affecting the reconstruction of the arrival direction for IC events, which for the shower-like topology is very large, namely of the order of 15 • . In particular, we treat the uncertainties on declination and right ascension as maximum errors, and propagate them on the quantity cos θ. Note that for galactic plane scenario the variable to be considered is the Galactic latitude b.
To consider in our statistical tests the above angular uncertainty, for each choice of 5 background events, we consider 100 possible extractions of the 7 remaining events from their maximum error intervals using a uniform probability. In this way, for the 100 different choices of observed events we compute the corresponding TS values, which once compared with the null hypothesis TS distribution, provide a range of p-values. Such a range is finally averaged on the 792 different background combinations. In Table 1 we report such an average range for each test. As we can see from Table 1  correlation with the galactic plane is disfavored. Note that in this case, the Anderson-Darling test is not well defined since its CDF is vanishing within the region b < b gal (see Eq. (14)). It is worth observing that varying the angular size b gal in the range [2 • , 4 • ] does not significantly change the p-value range reported in the Table. Moreover, the DM annihilation scenario is already excluded from IC data for both DM halo density profiles in case of a small clumpiness factor (∆ 2 0 = 10 4 ). On the other hand, for a larger clumpiness factor (∆ 2 0 = 10 8 ) we get a result similiar to the one of the astrophysical isotropic distribution. This is due to the fact that in this case the annihilating DM angular distribution is almost isotropic. It is worth observing that due to a certain lack of events from the Galactic Center, the NFW DM profile that is more peaked in this central region results to be more in tension with the observations than the Isothermal pro-file. This results in smaller p-values for NFW with respect to Isothermal as shown in the Table, such difference is exacerbated for annihilating DM scenario.

Forecast
It is of interest to ask about the statistics required (number of events) in order to distinguish, at a certain confidence level, a DM induced distribution from an isotropic one. To answer this question we perform a forecast analysis restricted to decaying DM scenario and annihilating DM one with ∆ 2 0 = 10 6 that are not already excluded by present data. For a given number of events, we generate 10 5 sets of data (in the 60 -100 TeV energy range) according to the isotropic distribution, and perform the two statistical tests under null hypothesis that the data samples come from a decaying DM distribution or from an annihilating DM one. For simplicity we assume that each data sample is not affected by the background. To include the background effect in the forecast analysis one can simply increase our "predictions" by a factor of ∼ 12/7 as suggested by present data.
By varying in the set of 10 5 data samples we get a distribution of p-value for which it can be defined the p-value at 68% Confidence Level (C.L.). This value represents the upper bound for p-values in 68% of cases. In Figure 3 we report the p-value at 68% C.L. as function of the number of signal events (no background) in case of decaying DM scenario. As expected, the Anderson-Darling statistical test (solid lines) is more appropriate than the Kolmogorov-Smirnov one (dashed lines). Indeed, the p-value falls down to zero very rapidly. Assuming that the p-value required to exclude a model is O(10 −3 ), we see that the decaying DM scenario will be completely excluded only when a O(200) number of signal events is collected in the energy bin 60 -100 TeV. It is worth noticing that the NFW density profile, since more spatially concentrated around the Galactic Center, requires a small number of signal events, namely O(100), to be excluded with respect to the Isothermal profile. The forecast analysis in case of annihilating DM scenario with intermediate value of clumpiness ∆ 2 0 = 10 6 is shown in Figure 4. In particular, we have obtained that in order to exclude such a scenario the required number of signal events is O(300). Such a huge statistics, even though cannot be reached in the present experimental set up, could be eventually reached in future Neutrino Telescopes [47,48].

Conclusions and outlook
In this Letter we have analyzed the 60 − 100 TeV bin of IceCube data that seems to suggest the presence of a ∼ 2σ excess once the sum of the background (atmospheric neutrinos and muons) and an astrophysical component described by a E −2 ν power-law is subtracted. In order to get information on the possible origin of such an excess, we have compared the distribution of the arrival directions of IceCube data with the angular distributions of simply distributed astrophysical galactic/extragalactic sources (galactic plane/isotropic distribution), as well as with the expected flux coming from DM interactions (decay and annihilation) for different DM profiles. The statistical analysis performed by means of Kolmogorov-Smirnov and Anderson-Darling tests of hypothesis seems to disfavor the correlation with the galactic plane, whereas excludes the annihilating DM scenario for both NFW and Isothermal density profiles in case of a small clumpiness factor (∆ 2 0 = 10 4 ). The small number of events till now collected does not allow to distinguish between the case of an isotropic distribution (astrophysical extragalactic sources and annihilating DM scenario with a large clumpiness factor) and the remaining cases (DM decay scenarios and DM annihilation ones with ∆ 2 0 = 10 6 ). In this concern we have performed a forecast analysis and we have found that O(200) (O(300)) signal events are required in order to exclude the decaying DM scenario (annihilating DM scenario with ∆ 2 0 = 10 6 ). If the lack of neutrino events towards the Galactic Center is confirmed by future data, the NFW decaying DM scenario will be more easily excluded.