Search for neutrinos from decaying dark matter with IceCube

With the observation of high-energy astrophysical neutrinos by the IceCube Neutrino Observatory, interest has risen in models of PeV-mass decaying dark matter particles to explain the observed flux. We present two dedicated experimental analyses to test this hypothesis. One analysis uses 6 years of IceCube data focusing on muon neutrino ‘track’ events from the Northern Hemisphere, while the second analysis uses 2 years of ‘cascade’ events from the full sky. Known background components and the hypothetical flux from unstable dark matter are fitted to the experimental data. Since no significant excess is observed in either analysis, lower limits on the lifetime of dark matter particles are derived: we obtain the strongest constraint to date, excluding lifetimes shorter than 1028s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{28}\hbox { s}$$\end{document} at 90% CL for dark matter masses above 10TeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10~\hbox {TeV}$$\end{document}.

cles to explain the observed flux. We present two dedicated experimental analyses to test this hypothesis. One analysis uses 6 years of IceCube data focusing on muon neutrino 'track' events from the Northern Hemisphere, while the second analysis uses 2 years of 'cascade' events from the full sky. Known background components and the hypothetical flux from unstable dark matter are fitted to the experimental data. Since no significant excess is observed in either analy-sis, lower limits on the lifetime of dark matter particles are derived: we obtain the strongest constraint to date, excluding lifetimes shorter than 10 28 s at 90% CL for dark matter masses above 10 TeV.
We present two dedicated analyses to test whether the description of the observed neutrino flux can be improved by an additional component from heavy (m DM > 10 TeV) dark matter decays as an alternative to bottom-up scenarios of astrophysical acceleration [33]. Such heavy particles are receiving increased attention because the classic WIMP paradigm of weak-scale mass dark matter is disfavoured by the negative results in searches for new physics at the LHC [34], in direct DM detection experiments [35][36][37][38][39], and in searches for DM annihilation into neutrinos [40,41] or gamma-rays [42][43][44][45][46].
Our results significantly improve upon the best previous experimental bounds on decaying dark matter obtained with gamma rays [44][45][46][47], neutrinos [48], and those derived from high-energy cosmic rays and the cosmic microwave background radiation [4,5].

IceCube detector and event selections
IceCube is a cubic-kilometer ice Cherenkov detector located at the South Pole, situated between 1450 and 2450 m below the surface [49]. Charged particles produced in neutrino interactions with the Antarctic ice or the bedrock below are detected by the Cherenkov light they emit, allowing the reconstruction of the originating neutrino's direction and energy [50].
The presented analyses use two different event samples. The first analysis is based on 6 years of ν μ charged-current data collected between 2009 and 2015, i.e., track-like events from the Northern Hemisphere. More details can be found in Ref. [2]. The second analysis uses 2 years of data collected between 2010 and 2012. The event selection is based on a previous study [51], modified to select only cascade events from the full sky which are produced in NC interactions or CC interactions of ν e or ν τ . Note that in the following no distinction is made between particles and anti-particles; the labels neutrino and lepton include the respective anti-particles and the used cross-sections incorporate both particles and antiparticles.
The two analysis samples are statistically independent, and while the track sample contains a much larger number of events, the full-sky coverage and better energy resolution of the cascade sample (see Table 1) lead to comparable sensitivities.

Analysis
To test whether the observed flux of high-energy neutrinos (partly) arises from heavy decaying dark matter, a forwardfolding likelihood fit of the distribution of reconstructed energy and direction is performed on both datasets, similar to Refs. [2,51]. The total observed flux is modelled as a sum of background and signal flux components. Each of these components is described by a parametrized flux template that depends on the fitted model parameters.

Flux components
Cosmic-ray air showers produce secondary mesons which decay into charged leptons and neutrinos. These atmospheric neutrinos are the main source of background in both data samples. They can be further divided into conventional atmospheric neutrinos produced by the decay of pions and kaons and prompt neutrinos produced by the decay of charmed mesons. This latter flux is sub-dominant at high energies and has not been separately identified yet [2]. Atmospheric neutrino flux predictions are taken from Refs. [52,53] for the conventional (modified to account for the cosmic-ray knee [2]) and prompt component, respectively. From the Southern Hemisphere, cosmic-ray induced atmospheric muons can also penetrate the ice, reach the detector and mimic a neutrino signal. After application of appropriate event selections, the atmospheric muon contamination is negligible in the tracklike sample and below 10% in the cascade sample. Astrophysical neutrinos from cosmic rays interacting in or near their production sites constitute a second background flux to the targeted signal of neutrinos from decaying dark matter. Since the origin of cosmic rays is unknown, an exact modelling of this astrophysical flux is not possible. A generic parametrization of these astrophysical neutrinos as an isotropic flux with a power-law energy spectrum agrees well with present measurements [1,2] and is therefore used in the fitting. The spectral index γ and the flux normalization Φ astro are taken as free parameters in the fit.
When heavy dark matter decays into standard model particles, neutrinos are necessarily expected in the final state [54]. Observing these neutrinos would thus constitute an indirect probe of the scenario of decaying dark matter. The energy spectrum, d N ν /d E ν , of the expected neutrinos depends on the exact decay mechanism and is model dependent. In this analysis, several "hard" (e.g., dark matter decaying directly into neutrinos [8,55,56]) and "soft" (e.g., neutrinos produced in the subsequent hadronic decay-chain of standard model particles [6]) decays are used as benchmark channels. Their spectra were simulated with PYTHIA 8.1 [57] and are shown in Fig. 1.
At Earth, the neutrino flux from dark matter decays has to be subdivided into a galactic and an extragalactic component. The expected energy distribution of the galactic component Φ Halo follows the initial decay spectrum. Its angular distribution incorporates the (uncertain) distribution of dark matter in the Milky Way halo via the line-of-sight integral [58]. The Burkert halo profile [59,60] with best-fit parameters from Ref. [61] is used as a benchmark and other halo profiles are considered as systematic uncertainties. The extragalactic neutrino flux from dark matter Φ Cosm. is expected to be isotropic and to have a red-shifted decay spectrum in energy. This flux is calculated adopting the ΛCDM cosmological model with parameters from Ref. [62]. The total signal flux is computed as the sum of both fluxes assuming that a single dark matter particle constitutes the observed dark matter in the universe. Additionally, neutrino mixing is applied with parameters from Ref. [63], the effects are shown in Fig. 2.
The total flux depends on two fit parameters: the mass m DM , which determines the energy cut-off, and the lifetime τ DM of the dark matter particle, which determines the normalization. Explicitly, it is given by

Likelihood analysis
In order to find the combination of the flux components that describe the data best, a forward-folding likelihood fit is performed. Flux templates, as a function of the fit parameters, are generated from a dedicated simulation of the detector response (see Refs. [2,51] for more details) and then compared to the observed event distributions in reconstructed energy E, right-ascension α, and zenith angle θ . Given a set of observed events, N , and the predicted number of events, μ i (ξ ), the Poisson likelihood is calculated and the fit parameters ξ are optimized, namely, While a binned likelihood method is used in the analysis of the track-like events, an unbinned approach is used in the analysis of the cascade sample, which corresponds to the limit of infinitesimal bin size.
To quantify the statistical significance of the best fit result, a test statistic is defined as the ratio of the maximum likelihood values for the background-only case (atmospheric and astrophysical fluxes) and for the background-plus-signal case (i.e., including the additional flux from dark matter decay), namely: Since the signal-plus-background case has additional degrees of freedom (four vs. two physical fit parameters), the TS value will always be positive. The observed TS value of the best-fit result is then compared to pseudo-experiments of the background and different signal hypotheses to construct confidence intervals.

Systematics
The systematic uncertainties of the two analyses arise from the modelling of the dark matter halo, the detector and the background fluxes. The dominant systematic uncertainty is the poorly understood dark matter distribution in our galactic halo. To investigate the resulting effect, the Burkert halo parameters are varied within intervals of one standard deviation while keeping their correlation fixed, by selecting β 2 = −0.5 (see discussion in Ref. [61]). In addition, the impact of a different halo profile, namely the Navarro-Frenk-White [64,65] profile, with best-fit parameters from Ref. [61], on the fit results is studied. The total effect of these halo model variations on the derived lifetime limit is ± 10%. This value is consistent across all the masses and decay channels and between the two analyses. The uncertainty on the extragalactic flux component, which arises from the average extragalactic dark matter density, is on the order of a few percent [62] and is thus not considered here. Detector simulation and background flux uncertainties are treated differently between the two analyses. In the analysis of track-like events, several nuisance parameters are fitted simultaneously in order to absorb deviations from the baseline expectation (see Ref. [2] for more details). They include the normalization of the prompt atmospheric flux, cosmicray flux model uncertainties, relative contribution from pion and kaon decays to the atmospheric fluxes, optical properties of the glacial ice, and the optical efficiency of the detector.
In the analysis of the cascade-like events, prompt atmospheric flux uncertainties [51], errors in the event reconstruction due to ice model uncertainties [66], a 10% uncertainty on the optical efficiency of the detector, and the impact of the finite simulation statistics are taken into account. The data are reanalyzed under different assumptions within the systematic uncertainties and the spread of the resulting limits is taken as the overall systematic uncertainty.

Fit results
To address the question of whether the observed flux of cosmic neutrinos can be described significantly better by including a component from decaying dark matter, the hard decay channels DM → H + ν (cascades) and DM → Z + ν (tracks) are fitted to the respective data. A dark matter signal would be expected to show up in both analyses. Also note, that the observable energy distributions are smeared out due to the limited detector resolution and the cosmological red-shift. It is therefore sufficient to fit these single decay channels in order to test whether a contribution from dark matter is present and multiple tests are not necessary. The obtained best-fit results and the corresponding p-values with respect to the background only hypothesis are listed in Table  2. The fits of the background-only hypothesis agree well with the results in Refs. [2,51]. Small differences arise due to a different choice of bins (tracks) and the altered selection (cascades). The corresponding best-fit distributions in reconstructed energy are shown in Figs. 3 and 4 together with the experimental data. Note that different energy estimators are used in the sub-samples (data-taking seasons) of the track analysis [67]. It is therefore not possible to show the experimental data in one histogram.

Interpretation of the fit results
Although the best-fit result in both analyses includes a nonzero dark matter component, the results are not significant (as both p-values are above 1%). More degrees of freedom in the modelling of the astrophysical flux, e.g. adding a second Cascade analysis: best-fit energy distribution for the signal hypothesis (components stacked to illustrate the dark matter component), with the best fit parameters listed in Table 2. The fit is performed on un-binned data, but for visualization purposes a binning is applied in the figure Track analysis: best-fit energy distribution. While the lowenergy events are well described by the conventional atmospheric component, the high-energy events are modelled by a combination of a weak diffuse astrophysical flux and a component from decaying dark matter (best-fit parameters in Table 2). The figure shows data recorded between 2012 and 2014 as they are based on the same energy estimator (see [67] for more details). The remaining years are fitted simultaneously but are not shown here component, would further reduce the significance. Thus, the result is not interpreted as a signal of dark matter decay. Furthermore, a dark matter signal should be constant in time but the fit of the track-like events shows fluctuations; see cal normalization is significantly reduced compared to previous results [2]. A dark matter only scenario, where the normalization of the astrophysical flux is zero, is however disfavoured by 2ΔL L H 1 compared to the best-fit point.
As expected, the best-fit dark matter mass that induces a cut-off in the energy spectrum is found to be independent of the diffuse astrophysical normalization while the dark matter normalization is anti-correlated. Dark matter lifetime / sec Fig. 7 Dark matter lifetime limits for all considered decay channels. For the Z + ν/H + ν channel, the limit was combined (solid grey line) as described in the text. Between m DM ∼ 10 5 GeV and m DM = 1.5 × 10 7 GeV the limit is obtained from the more sensitive track analysis. The limit from the cascade analysis is shown as a dashed line and turns out to be stronger above m DM ∼ 5 × 10 7 GeV

Lifetime limits
Since no significant dark matter signal is observed, lower limits on the lifetime of the dark matter particle (corresponding to upper limits on its signal strength) are derived. In order to combine the two analyses, the lower limit on the lifetime is extracted from the respective analysis with the better sensitivity (median limit obtained from background pseudoexperiments) at each dark matter mass. The hard decay channels Z + ν (track analysis) and H + ν (cascade analysis) are treated as the same channel because the resulting neutrino spectra are indistinguishable within energy resolutions. Further, limits for the decay channels νν, τ + τ − , μ + μ − , W + W − and bb are calculated only in the cascade analysis because the energy resolution of the track analysis is not sufficient to differentiate those channels from each other. The resulting lower limits on the dark matter lifetime are shown in Fig. 7. Note that for the bb decay channel, the lower limit on the lifetime increases steeply with the dark matter mass because QCD fragmentation generates a soft tail of low-energy neutrinos (see Fig. 1) which become increasingly relevant for large dark matter mass. Furthermore, no limit on the lifetime is calculated in this channel for m DM below 10 5 GeV because the resulting decay spectrum becomes similar to the atmospheric background fluxes and the respective uncertainties would have a major effect on the obtained limit. The enhanced limits at m DM ∼ 10 7 GeV, correspond to the nonobservation of electron neutrino events from the expected Glashow resonance [68]. For the track-like sample, all nuisance parameters are fitted to their expectation values within one standard deviation, and the effect on the signal hypothe- bb μ + μ − Fig. 8 Comparison of the lower lifetime limits with results obtained from gamma-ray telescopes: HAWC (Dwarf Spheroidal Galaxies) [44], HAWC (Galactic Halo/Center) [45] and Fermi/LAT [47] sis is found to be negligible. For the cascade-like sample, the overall impact of the systematics is roughly 10-15% for dark matter masses below 5 PeV and 1% for those above it. The limits shown here include a degradation due to 1σ systematic variation.

Conclusions
Two analyses on statistically independent datasets searching for a contribution from decaying dark matter to the astrophysical neutrino flux have been presented. It has been shown that the observed high-energy neutrino flux can be described equally well by a combination of a dark matter component and a diffuse astrophysical flux with a power-law energy spectrum. However, neither analysis identified a significant dark matter excess in the data, and models in which the cosmic neutrinos flux arises entirely from dark matter decay are disfavoured.
From the non-observation of a dark matter signal, lower limits are set on the lifetime of dark matter particles with mass above 10 4 GeV. For such heavy particles these limits are presently the strongest on the dark matter lifetime (see Fig. 8).