Geometric compatibility of IceCube TeV-PeV neutrino excess and its galactic dark matter origin

We perform a geometric analysis for the sky map of the IceCube TeV-PeV neutrino excess and test its compatibility with the sky map of decaying dark matter signals in our galaxy. We have found that a galactic decaying dark matter component in general improve the goodness of the fit of our model, although the pure isotropic hypothesis has a better fit than the pure dark matter one. We also consider several representative decaying dark matter, which can provide a good fit to the observed spectrum at IceCube with a dark matter lifetime of around 12 orders of magnitude longer than the age of the universe.

The angular resolutions of those events play an important role for identifying the geometric origin of the neutrino excess. The IceCube collaboration has performed a point source analysis for the 4 years cascade events and found that there is no significant evidence of spatial clustering and the p-value for the hypothesis of a uniform event distribution is 18% [7]. Curious about the possible linkage between the TeV-PeV neutrino excess at IceCube and the mysterious DM in our universe, in this paper we analyze the IceCube data with a special attention on its geometric distributions, and study the statistical significance of its potential DM origin, which prefers to have more signal events around the galactic center because of the DM spatial profile.
One could consider DM annihilation as an explanation. However, due to the unitarity bound [14,15], we found that the annihilation rate for a DM mass around one PeV is about four orders of magnitude lower than the required one for the IceCube data. Therefore, we concentrate on the decaying DM case, which can match to the required rate for a DM lifetime of 10 28 -10 29 s. In this paper we do not provide a theoretical understanding of the DM mass scale and the decay lifetime, but we want to point out that a heavy DM with a non-thermal history has been widely predicted in many models [16][17][18][19].
In our following statistical analysis, we will perform both the Kolmogorov-Smirnov test and the likelihood ratio test using the official release of the IceCube data for the four years [7]. The KS test will answer how compatible of different model hypotheses with the experimental data. The likelihood test will concentrate on comparing different model hypotheses directly. Those two tests are complimentary to each other: the KS test could exclude all model hypothesis using the features in experimental data; the likelihood test may tell us the individual weights of a combination of models.

Geometric analysis for decaying dark matter based on the Kolmogorov-Smirnov test
Our main goal is to study the compatibility of the neutrino sky map from DM and the sky map of the observed events at IceCube. The signal distribution from DM decays depends on the DM spacial profile in our galaxy. For the Einasto profile [38], one has ρ DM (r) = ρ e  0.3 GeV cm −3 is the approximate DM density in the solar system. The neutrino signal from DM decay is calculated by the line-of-sight integral along a given direction [18] where the integral of s is along the line of sight and the relation between r and s is r 2 = s 2 + r 2 − 2s r cos l cos b, where −90 • ≤ b < 90 • and −180 • ≤ l < 180 • as the latitude and longitude angles in the galactic coordinate. τ DM is the DM lifetime and m DM is the DM mass. The normalized neutrino differential spectrum is dN/(N dE ν ). The integrated neutrino flux from DM is For the integrated time of 1347 days and 10 m 2 · sr acceptance area for the energy around 100 TeV, there could be around 20 events observed at IceCube. The geometric distribution of the IceCube events is represented in the equatorial coordinate. We, therefore, translate the DM generated event distribution from the galactic coordinate in the latitude and longitude angles (b, l) to the equatorial coordinate in the declination angle and the right ascension angle (δ, α) (see ref. [39] for details). We define the DM probability distribution using the normalized flux with the DM event sky map shown in the left panel of figure 1. For all or subsets of the observed 53 events from IceCube, we construct the data probability distribution using the solid-angular error σ i for each event by assuming a Gaussian distribution

JHEP01(2016)161
where ∆R(δ i , α i ; δ, α) is the angular distance between the points (δ i , α i ) and (δ, α) on the sphere. In the right panel of figure 1, we show the sky map of the observed N = 53 events at IceCube after implementing the angular resolution for each event. Comparing these two maps, one can see that both have a concentration of events around the galactic center direction. On the other hand, the DM sky map has very few events in the right and upper corner, while the IceCube data map has some population in this region.
To quantify the similarity of the two sky maps in figure 1, we perform a statistical test to calculate the p-value of the hypothesis of decaying DM as an explanation of IceCube neutrino excess. We first use a two-dimensional version of the Kolmogorov-Smirnov (KS) test statistics (TS) [40] to study the compatibility between the data and the DM hypothesis. We will use the maximum likelihood-ratio test as well later. The KS test statistics is defined as the largest absolute difference between cumulative probability distributions of the data and the model. It takes better account of the relation among data points than the traditional likelihood-ratio test.
To make the definition of the TS less sensitive to the integration directions, we consider a set of four possible integration regions, for a given boundary choice (δ 0 , α 0 ). The TS or the difference of the cumulative probability distributions is defined by Choosing the largest value for all possible boundary choices, we have the KS test statistics as To calculate the p-value for the decaying DM as an explanation for the data, we generate random event maps by choosing a random (according to the model profile) rightascension angle but keeping the same declination angle and resolution of the event in the data. Notice that scrambling the events only in the right-ascencion does not introduce any bias in the analysis, having one of the coordinates fixed, it only incorporates an equal shift in the χ 2 for every random realization and the data. If the data is not in a very unlikely realization the effect of using only right-ascencion in the p-value calculation in general it reduce the power for the test statistic moreover we chose to do that to avoid possible bias introduced by the systematic effects such that earth absorption in the zenith distribution. This same technique is used for the point sources analysis done by the IceCube collaboration [5,6] with the same data.
In  test statistics TS(DM) of DM against the observed 53 events at IceCube. The p-value, or the probability of having TS(DM) smaller than the TS value from a random event map, is 36.8% for the Einasto model withᾱ = 0.17. To test how good the observed 53 events agree with a isotropic geometrical distribution, we perform the same calculation by assuming a isotropic model (in the right panel of figure 2) and found that the p-value for a isotropic distribution is 49.8% for all 53 events.
Since the atmospheric backgrounds are dominated in lower energies [8,9], a bigger fraction of the observed events could be from DM signals if only relatively high energy events are selected. Therefore, we also test the geometric distributions for the 34 events with E 50 TeV. We show the p-values for all 53 events and the 34 events with E 50 TeV in table 1. One can see that the p-values are fairly insensitive to the energy cut. In the last row of table 1, we also show the p-values for only the cascade events considering the fact that the track events could have an origin from the atmospheric muon background. From table 1, one can already see that there is no dramatic difference betweenᾱ = 0.25 andᾱ = 0.17 cases. This is due to the poor angular resolution of cascade events such that the peaked center of the DM profiles can not be resolved. The increase of the p-values for the isotropic distribution from all 53 events to 39 cascade events is due to the extremely good resolution of the 14 track events.

JHEP01(2016)161 3 Geometric analysis using the likelihood ratio test
In the last section we saw that all the hypothesis are statistically compatible with the different sets of data. The next question to ask would be how much of the signal events may come from the different hypothesis and how can we exclude only Isotropic or only DM against a combination of the two. In this section we check the compatibility of the data with the DM profile using a likelihood ratio test, which was used by the IceCube collaboration in their point source analysis [41]. We first treat the isotropic distribution as the null hypothesis with an alternative isotropic plus DM hypothesis. We define the likelihood function as where B i is the isotropic background contribution and S i is the signal DM contribution. 2 The n s is the number of signal events and will be used to maximize the likelihood. We use the observed data locations and errors convoluted by the DM probability distribution to calculate the signal contribution S i as Here, x i is a vector, defined in the (δ, α) plane, from the location of the observed event and σ i is the corresponding angular error. We show the log-likelihood function as a function of signal strength n s /N in figure 3 for three cases: all 53 events, 34 events with E 50 TeV and 39 cascade events. We can see from figure 3 that the preferred values of n s are positive, which suggests that a combination of DM plus isotropic distributions fit the data better than the isotropic-only fit. Comparing the best fitted values of n s for α = 0.17 andᾱ = 0.25, one can see that a larger value ofᾱ or a flatter DM profile prefers more DM signal events.
To quantify the p-value of the data to reject the isotropic-only hypothesis against the isotropic plus DM hypothesis, we calculate the test statistic as TS = max

JHEP01(2016)161
the current data. On the other hand the data is still not enough to exclude the DM hypothesis. In this case we would expect some galactic and extragalactic DM fluxes that may be comparable [42]. This estimation is still consistent with the best fit number of signal events that we see in the current data, table 2. In the future the IceCube [43] extension may bring enough statistics to potentially exclude the DM origin of the majority neutrino events.

Neutrino spectra from dark matter decays
The energy spectrum of the IceCube neutrino excess has interesting features [5]. First, there are three events with deposited energy above 1 PeV [6]. Secondly, there is an potential energy cutoff below the Glashow resonance energy. Although a wide range of the energy spectrum can be fit by an E −2 feature [5], it is still interesting to explore potential DM produced spectra from particle physics.
To fit the observed spectrum at IceCube, one also needs to consider different detector acceptances at different energies. For different flavors of neutrinos, the acceptance areas vary a lot with the largest one for the electron neutrino. In our analysis below, we don't distinguish different flavors of neutrinos and use the averaged acceptance areas in terms of flavors and declination angles [5], which for the flavor part is equivalent to assume the 1:1:1 combination. This effective area is only slightly different from ref. [24].
Because the uncertainties on the acceptance areas and the large statistical errors, the current IceCube data is not sufficient to distinguish spectra among different particle physics models. So, we consider several representative decaying DM models and study their fit to the observed energy spectrum. We consider candidate models according to the operator dimensions of DM coupling to SM particles.
At the renormalizable level and for a fermion DM χ, we consider the operator λHL L χ for DM coupling to the Higgs field in the SM or λH LLL χ in the lepton-specific two-Higgs doublet models, which has DM decays as χ → h + l and χ → l + H L → l + τ + + τ − , respectively. Setting the fermion DM mass to 8 PeV, it is possible to fit the spectra observed by IceCube. For a scalar DM, one can have the renormalizable coupling to the SM Higgs boson as simple as µ XHH † , which simply mediates the decay of X → 2h. Because the neutrino spectrum is relatively soft, the mass of DM has to be raised to 20 TeV for a reasonable fit. Beyond the renormalizable level, one could have DM mainly couple to two leptons via m τ Xτ + τ − /Λ, so the decay channel is X → τ + τ − . Such DM with 8 TeV mass can fit the IceCube spectrum well. We show the fitted spectra in figure 5 after using PYTHIA [44] for SM particles decay and hadronization. (see [19,20] for other spectra from DM decays).

Conclusions and discussion
Our geometrical analysis has already shown that a combination of the galactic DM contribution and an isotropic spectrum provides the best fit to the data. Because there are neutrino events pointing to the direction where the DM density is expected to be quite JHEP01(2016)161  Figure 5. The fitted spectra for several DM decay channels. The black and solid line is the atmospheric backgrounds [8,9]. For the two fermion DM cases, the DM mass is 8 PeV and both lifetimes are τ χ = 2.3 × 10 28 s. For the scalar DM decays into two higgses, the DM mass is 20 PeV and the lifetimes are 4.6 × 10 28 s. The DM decaying into τ − + τ + is 8 PeV, with lifetime 2.3 × 10 28 s, respectively. low, a purely galactic DM origin for the 53 events is excluded unless a flatter DM spacial profile like the isothermal one is used. It is also important to notice that extra-galactic DM may also produce an extra homogeneous contribution to the signal. A more detailed analysis of the DM distributions in the universe is needed to address this possibility.
If the significance of the source at the center of our galaxy keeps growing as the data accumulates, decaying dark matter will remain a strong candidate for at least part of the TeV-PeV neutrino events observed by the IceCube experiment. The IceCube has more data to be analyzed and collected, so a more robust conclusion can be drawn in the coming years. Other than IceCube, another neutrino telescope, ANTARES [45], has reached a comparable sensitivity in some declination angle region. A geometric test for the compatibility between the neutrinos (excess) observed in ANTARES and a decaying DM will be interesting.
Beyond the neutrino signal from DM, one could also search for other correlated and for sure model-dependent cosmic ray signatures from the DM decays at other experiments like Fermi LAT [46,47], PAMELA [48,49], AMS-02 [50] and HESS [50]. In the few respective models considered, additional photons, positrons and antiprotons can be produced at the same time when a neutrino signal is generated. One can easily check that the neutrino yield is considerably higher than the photon, positron and antiproton yields in every bin. Furthermore, because of the long DM lifetime of 10 28 -10 29 s, the predicted photon, positron and antiproton fluxes have been checked to satisfy the current cosmic ray constraints (see ref. [35] for the gamma ray prediction at Fermi-LAT for the hadronuclear origin of the IceCube neutrino excess).
The PeV scale DM considered here is definitely beyond the scope of high energy collider searches. If additional interactions exist between DM and quarks, the direct detection experiments may see a signature [51]. If the IceCube excess is indeed due to decaying DM, a new avenue to understanding the DM properties will be opened.