Reconciling neutrino flux from heavy dark matter decay and recent events at IceCube

The IceCube detector has recently reported the observation of 28 events at previously unexplored energies. While the statistics of the observed events are still low, these events hint at the existence of a neutrino flux over and above the atmospheric neutrino background. We investigate the possibility that a significant component of the additional neutrino flux originates due to the decay of a very heavy dark matter (VHDM) particle via several possible channels into standard model particles. We show that a combination of a power law astrophysical neutrino spectrum and the neutrino flux from the decay of a DM species of mass in the range $150-400$ TeV improves the fit to the observed neutrino events than that obtained from a best-fit astrophysical flux alone. Assuming the existence of an astrophysical background described by the IC best-fit, we also show that, for the decay of even heavier DM particles ($m_{\text{DM}} \sim 1$ PeV), the same observations impose significant constraints on the decay lifetimes. Allowing the astrophysical flux normalization to vary leads to modifications of these limits, however, there is still a range of dark matter mass and lifetime that is excluded by the IC results.


I. INTRODUCTION
The IceCube detector (IC) at the South Pole has recently observed 28 events at energies 30 TeV -1.2 PeV [1]. This is the first ever observation of neutrino events at energies beyond the tens of TeV at any detector and heralds a new chapter in understanding the physics responsible for producing and accelerating fundamental particles to such high energies. The commonly accepted means for the production of fundamental particles at these socalled ultra-high energies (UHE) is at the extremely reactive cores and jets of astrophysical sources [2] beyond our galaxy, where the already fast-moving charged particles can be accelerated further due to Fermi acceleration [3] when traversing through the shock waves associated with these sources. Standard theoretical arguments [4] deduce that the diffuse flux of neutrinos arriving at the earth from the such astrophysical sources should have a power-law energy spectrum proportional to E −2 . Although, admittedly, statistics are low, the recent observations at IC hint at incident neutrino fluxes largely in keeping with these expectations -corresponding determination of the power-law ( dΦ/ dE ∝ E −α ) spectrum best matching these events indeed gives α = 2.0. However, there are some noteworthy incompatibilities between the observed events and the corresponding best-fit power-law spectrum: a) the expected number of events in the lower energies 30-150 TeV falls consistently below the observations, and b) at two high energy bins marked by a * atrib@email.arizona.edu † mary-hall-reno@uiowa.edu ‡ ina@physics.arizona.edu distinct lack of observed events, the E −2 flux with the best-fit normalization obtained by IC gives event rates that are simply too high. While it might be entirely possible, indeed even likely, given the small sampling of data we are working with, that these differences are merely statistical fluctuations which will disappear with the accumulation of more data, the systematic nature of the difference between theory and experiment, especially at the lower energy bins, suggests that the total neutrino flux incident at the earth might comprise a component additional to that from astrophysical sources. Given the extremely high energies in question, the sources that can produce a neutrino flux that should be detectable at IC are rather few. One interesting source, apart from astrophysical ones, that is capable of producing UHE neutrinos is from the interactions or decays of a very heavy dark matter (VHDM) species.
The fact that a sizable fraction, about 27%, of the universe's energy density is made of non-luminous matter is now well known, and has been recently confirmed by the PLANCK observations [5]. The most consistent explanation of this large amount of energy not accounted for by standard baryonic matter is given by postulating the existence of one or more, usually heavy, dark matter particles which are relatively inert to interactions with the standard model particles [6]. Efforts to experimentally detect such particles, called Weakly Interacting Massive Particles (WIMP), are ongoing [7][8][9][10][11][12]. While these direct search experiments are constrained to focus on detecting WIMP-like particles with masses in the range 10 GeV -1 TeV, DM species in nature might be significantly heavier. For thermal relics, the mass is constrained by unitarity to lie below ∼ 500 TeV [13,14]; however, nothing pre-vents non-thermal DM masses to be even as high as a 10 9 TeV [14]. Indirect searches such as the detection of products of DM annihilation in the sun at large volume detectors like the IC are sensitive to heavier DM species with masses above the TeV's as well.
In this paper, without going into the specific models that predict the existence of such a candidate, we assume the existence of a DM species in the universe with its mass in the range 100 TeV -1 PeV. If it exists, both annihilation and/or decay (if it is unstable) of such a VHDM particle, will lead to the production of standard model particles eventually generating neutrinos as secondaries. In this paper, we will investigate the possibility that neutrinos from the decay of such a species might be able to assuage some of the tension between observation and theoretical predictions as outlined in the previous paragraph.
Specifically we will focus on a VHDM species which is unstable, and can, therefore, decay but with lifetimes well beyond the age of the universe -present constraints from gamma ray observations tie the DM lifetimes to 10 26 s [15,16]. In general the DM could decay into different possible SM final states, including charged leptons, neutrinos, quarks and weak interaction bosons W and Z. Limiting ourselves to two-body decays, we will do a generic model independent study by analyzing the spectrum of neutrinos produced in each of the possible decay channels individually. For each of the channels, we will evaluate the number of events expected at the IC from a total flux comprising neutrinos from DM decay and an astrophysical power-law flux. Finally we will show that, for a reasonable set of parameters for each decay mode, the event rates from the total (DM + astrophysical) diffuse flux gives a significantly better statistical match to the observed events than does the power-law spectra of the astrophysical flux alone. In doing so, we will determine the DM mass and lifetimes that lead to event rates most closely in match with the IC observations. The paper is organized as follows. In sec. II we discuss models of decaying dark matter in existing literature and we calculate the resulting neutrino flux for the various decay channels. For each decay mode in the different models, we calculate the resulting neutrino flux expected at earth after standard oscillation amongst the three flavours. In sec. III, we compute the net all-sky neutrino flux incident at IC by combining the flux from DM decay obtained in the previous section with an astrophysical power-law spectrum and use this flux to calculate the event rates expected at IC. In sec. III A we do a full parameter-space scan in the DM mass-lifetime plane to show the region of the plane that is excluded by recent observations. By doing a χ 2 analysis of the predicted event rates against observations at IC, in sec. III B we derive the values of the DM mass and lifetime that give the closest match. Finally, we draw our conclusions in sec. IV.

II. DARK MATTER DECAY MODES AND THE RESULTING NEUTRINO FLUX
DM candidates in several models are unstable and can decay to standard model particles. Decaying dark matter is a favored solution for the explanation of the e ± excess in cosmic radiation seen by PAMELA [17] and ATIC [18]. Indeed, a certain class of DM models, called leptophilic DM in which the DM particles annihilate or decay preferably to the charged leptons rather than to quarks, was first proposed [19] to explain these excesses. Specifically, in the case of decay [20][21][22], the DM particle can produce an excess of leptons via DM → + − , where is one of e, µ or τ . Decays to lighter charged leptons are generally more strongly constrained from the positron excess measurements by Fermi Large Area Telescope [23], and AMS-02 [24], however, and decays to e − e + pairs are strongly disfavored (see, e.g., [15]). DM particles might also decay to neutrinos producing a sharp resonance in the spectrum at E ν = m DM /2 [25,26]. In addition heavy DM candidates from supersymmetry such as the gravitino (ψ 3/2 ) can also decay, e.g., to γν, but more preferably to W ± ∓ and Z 0 ν.
Here, we do a largely model-independent study of dark matter decay in the context of the recent events seen at IC, restricting ourselves to studying two-body decays of a bosonic DM. We consider each of the decay channels and, in each case, compute the energy spectrum, dN ν /dE ν , of neutrinos produced due to the decay of a DM particle by using the event generator PYTHIA 8.1 [27], taking care to include electroweak corrections [28]. We do not consider decays to e ± pairs because these are much more strongly constrained than any of the four channels considered above. We also do not consider decays to qq pairs as the secondary neutrino flux produced when the DM particle decays to quarks in the relevant energy range is low, and does not give observable event rates at the IC. Neutrinos are produced as secondaries in all the four channels, in the case of c) and d) due to the fragmentation of the Z and W bosons.
The total neutrino flux from the decay of DM is composed of contributions from interactions within the galactic halo and from outside the galaxy [29], i.e., where, Φ G and Φ EG represent the galactic and extragalactic components of the total neutrino flux respectively. The differential flux from DM decay in the Milky Way halo is given by where, m DM and τ DM represent the mass and lifetime of the dark matter particle, l and b represent the galactic coordinates of the place where the interaction occurs and ρ represents the density profile of dark matter within the galaxy. The integral is over the line-of-sight parameter s, and it it related to the distance from the galactic centre r by R being the distance of the sun from the GC, R = 8.5 kpc. Using the NFW [30] profile for the dark matter distribution in the galaxy, the all-sky differential neutrino flux from DM decay in the galaxy is given by (see, e.g., [29], [31]) where, Additionally, the extra-galactic component of the differential neutrino flux is given by the expression with Here, z represents the red-shift of the source, ρ c = 5.6 × 10 −6 GeV cm −3 denotes the critical density of the universe, and we have used H(z) = Ω Λ + Ω m (1 + z) 3 , and Ω Λ = 0.6825, Ω m = 0.3175, Ω DM = 0.2685 and H 0 = 67.1 km s −1 Mpc −1 from the recent PLANCK data [5].
Having already computed the neutrino spectrum per decay, dN ν /dE ν , we can now directly use it in eqns. (1), (3) and (4b) to calculate the total neutrino flux produced due to each decay mode of the DM.

Effect of neutrino oscillation
As the neutrinos propagate to the earth, the three flavors of neutrinos oscillate among themselves, leading to an averaging out of the flux ratios. Decays to τ ± pairs or any of the gauge boson leads to the production of an equal number of neutrinos of each flavor at source. Due to the effect of oscillation, a flux ratio of ν e : ν µ : ν τ = 1 : 1 : 1 remains unchanged while propagating to the earth. Decays to µ ± pairs on the other hand leads to a neutrino flux of 1 : 1 : 0 at source, which gets modified to 0.785 : 0.607 : 0.607 due to oscillation during propagation. We take the effect of oscillation into account when calculating the neutrino flux reaching the earth from the DM decay as well as that from astrophysical sources.

III. EVENT RATES AT IC FROM DM DECAY AND ASTROPHYSICAL NEUTRINO FLUXES
We assume that the astrophysical neutrino flux arrives at the earth essentially in a 1 : 1 : 1 flavor ratio [32] and follows an unbroken power law, i.e., Here, k is a normalisation constant and α 2 is the spectral index of the spectrum. We fix the spectral index at α = 2.0 as suggested by Fermi shock acceleration in the jets and cores of high energy astrophysical object, and as determined from preliminary fits to the IC data [1] [33]. To evaluate the total event rate expected at IC, we sum the total neutrino flux for each flavor at earth from DM decay and astrophysical sources, i.e., where, Φ λ Total indicates the total flux at earth for the neutrino of flavor λ, with λ = e, µ, or τ .
The IceCube has recently carried out a survey of very high energy contained events, i.e., both cascades and muon track events with their starting vertex located within the body of the detector. Based on the 28 events detected over a run-time of 662 days, the IC finds the best-fit to the astrophysical E −2 flux to be The statistical degree of match of this best-fit flux with the observed events is given by the χ 2 : the smaller this value, the better the match. For the IC best-fit flux Φ IC , χ 2 = 10.7. Note that, while IC can only measure the energy deposited in the detector, we have conservatively assumed these energies are equal to the incident neutrino energies. This is only strictly true for ν e charged-current events; for all other events the deposited energy represents the minimum energy of the incident neutrino.
Using the fluxes obtained above, along with the effective exposure areas for contained events over a period of 662 days of IC's run-time given in [1] for each neutrino flavor, we calculate the total events expected for a given set of the parameters m DM , τ DM , and k, for each of the possible decay modes of the DM particle.

A. Accessible DM parameter space
To analyze the accessible region of the m DM − τ DM parameter space, in view of the IC events, if the total flux reaching the detector is a combination of fluxes from DM decay and astrophysical sources, we fix the astrophysical flux at the IC best-fit (and given by Eq. (7)), and carry out a scan of the DM parameter space by varying the mass and lifetimes in the range 100 TeV m DM 1 PeV and 10 26 s τ DM 10 29 s respectively. For each set of parameters we calculate the total event rates expected at the IC, and compare this to the observed event spectrum at IC by means of a χ 2 test. By doing this in turn for each DM decay channel, we plot the parameter space plane and show the region of this plane disfavored at 5σ. As previously noted, the IC fit is in significant tension with especially the lower energy events (30 -100 TeV). Physically, this analysis represents the scenario where the decay of an O(100 TeV) DM particle augments the low energy events from the astrophysical flux to provide a significantly better match to the observed events, while the high energy events closer to the PeV are explained by the latter alone. It is, therefore, obvious that fixing the astrophysical flux at the IC best-fit tends to push the accessible DM parameter space closer to the low m DM regions -observed event numbers at the high energies are either already consistent with those expected from the astrophysical flux alone or, as in the case of the two high energy bins immediately below 1 PeV, somewhat less, and added events from DM decay only serves to increase the mismatch in such a scenario. This is borne out by the allowed DM parameter space at 3σ, which lies in the region around m DM ∼ 200 TeV. For heavier DM with mass m DM 500 TeV, decay lifetimes are constrained to be more than 10 27 s at 5σ. The results are shown in fig.  1.

B. Best-fit parameters
In this section, we derive the best-fit values for the DM decay and astrophysical flux parameter space.
For each of the possible DM decay channels we keep the astrophysical flux normalization fixed at the IC best-fit value, i.e., k = E 2 dΦ IC / dE = 1. respectively [34]. Thus considering each decay channel in turn, we calculate the number of events expected due to a sum total of the astrophysical flux and that from the decay. The resulting best fits and χ 2 representing the degree of match are shown in table I. Event rates corresponding to the best-fit parameters are shown in fig. 2. It is evident from the figure that, especially, at the lower energies, i.e., 30 TeV E 100 TeV, the combined DM and astrophysical flux gives a significantly better fit to the observed data than the IC best-fit astrophysical E −2 flux.
To illustrate the dependence on the astrophysical flux normalization, we also carry out a similar analysis by fixing it at a reduced value, e.g., i.e., one-half of the IC best-fit value. The results (also shown in fig. 2) illustrate the significantly better match that this gives with the high energy events while still being consistent with events at the lower energies. It is clear that, if neutrino fluxes arrive both from DM decay and astrophysical sources, a reduced normalization for the latter is more consistent with the observation, including the high energy bins immediately below 1 PeV where no events are seen. Significantly, an astrophysical flux with its normalization at the IC best-fit seems to be in some disagreement with these "gaps" in the event spectrum. The best-fit parameters obtained with this reduced astrophysical flux are shown in table II. The value of the spectral index in these analyses is determined by the event spectrum at the higher energies ( 200 TeV), where the astrophysical flux is the sole contributor to the incident neutrino flux. As such the best match to the events in these energies comes from an incident astrophysical flux that goes as E −2 , and we have, therefore, chosen to fix the α at this value for our analyses. Softer flux spectra with α 2 would lead to event rates that fall off quickly and, for large enough α, become too small to explain the two events seen at energies above 1 PeV. If the incident neutrino flux had a softer spectral index not very different from α = 2 (so that it were still consistent with the high energy events), the lower energies, where fluxes from both DM decay and the astrophysical power-law spectrum contribute, would witness a slight increase in the relative contribution of the latter over the former, and, a similar analysis would, consequently, yield correspondingly larger best-fit DM lifetimes.
TABLE I: The best fit parameters for each of the DM decay modes and a comparison with the IC best-fit power law spectrum. The astrophysical flux normalization is fixed at the IC best-fit (see Eq. (7)).  (7)). The red shaded region is ruled out at 5σ, while the green patch shows the region of parameter space consistent with data at 3σ. To compare with existing bounds on lifetimes from gamma ray observations, we show the bound obtained in [16] when the DM decays to a W ± pair (black dotted curve) in the top panel plots (DM → Z 0 Z 0 and DM → W + W − ), and that obtained when the DM decays to µ + µ − (black dot-dashed curve) in plots in the bottom panel (DM → τ + τ − and DM → µ + µ − ), with the region below the curves excluded in both cases. The best fit parameters for each of the DM decay modes and a comparison with the IC best-fit power law spectrum. A reduced normalization for the astrophysical flux, given by Eq. (9), is assumed to illustrate the better match to the observed event rates it provides, as evident from the significantly lower χ 2 values (when compared against those in table I, and certainly when compared against that obtained from the IC best-fit).

IV. CONCLUSION
We have shown that if the events seen at the ultra-high energies at the IceCube are understood to be due to a combination of neutrino fluxes from both DM decay and astrophysical sources, the observations strongly favor a DM mass in the range 150-250 TeV, with decay lifetimes of O(10 27 s). By considering important decay pathways of the DM particle, viz., i) DM → Z 0 Z 0 , ii) DM → W + W − , iii) DM → τ + τ − , and iv) DM → µ + µ − , we have done a largely model independent scan of the m DM − τ DM parameter space to check for consistency with IC data. We see that while, as expected, the channel DM → µ + µ − is the most strongly constrained of the four, the recent events still allow for consistency with a significant region in the parameter space around m DM ∼ 200 TeV and reasonable lifetimes for all four channels.
Despite the limited statistics that the observation of only a score of high energy events entails, we see that a total neutrino flux with contributions from DM decay and astrophysical power law spectrum provides for a better match to the observation than that from an astrophysical source only flux, no matter what spectral shape one assumes for the latter. Especially at the low energies around 30 -100 TeV, where the best-fit astrophysical flux from the IC analysis predicts event rates consistently less Signal-only event rates (atmospheric backgrounds have been subtracted out) corresponding to the best fit parameters for each decay channel as listed in tables I and II with the astrophysical flux normalization set at the IC best-fit (magenta lines) and at half its value (blue lines) respectively. Event rates from the IC best-fit flux are also shown (green dotted lines), as are the actual event numbers seen at IC along with their associated errors (red dots).
than that observed, the contribution of the neutrino flux from the decay of a ∼ 200 TeV DM species is significant and serves to significantly enhance the consistency between theory and experiment.
On the other hand, a neutrino flux augmented by DM decay allows a significantly reduced astrophysical flux to still be consistent with the low energy events, while, concurrently inducing notable improvement in consistency with events at the high energy (E mean 400 TeV) bins. We show, e.g., that a reduced E −2 astrophysical flux at half the IC best-fit provides a much better match with observations, including being consistent with the lack of events in some of the high energy (E mean ∼ 500 and 800 TeV) bins, to within statistical errors. The marked improvement in consistency with observations in this case is also borne out by the significantly reduced measures for the χ 2 .
As the recent events at the IC show, the existence of an astrophysical flux of neutrinos looks increasingly true. However, notwithstanding the low statistics thus far, the IC events do seem to indicate that a simple power law spectrum for the astrophysical flux cannot, by itself, ex-plain these observations completely. Accumulation of events in the IC as it collects data in the near future will serve as an important test for the consistency of the observations with the astrophysical flux predicted from standard calculations. If, akin to that seen with the first set of UHE events, inconsistencies between those predictions and observations continue to persist, it will perhaps bring into consideration the existence of additional component(s) to the neutrino flux arriving at the earth. At these remarkably high energies, one of the (rather small set of) possible neutrino sources is from the interactions, either decay or annihilation, of a dark matter species with its mass O(100 TeV). While DM annihilation would not contribute a statistically observable number of events at the IC, decays to standard model particles with lifetimes O(10 27 s) lead to detectable neutrino fluxes. Given the limited statistics at present, we have shown that a combination of neutrino fluxes from DM decay and astrophysical sources proves to be more consistent with observations than an astrophysical power-law spectrum alone. We have also shown, in addition, that it might well be possible that the magnitude of the actual astrophysical flux seen by the IC detector might be as low as one-half the present best-fit value. While, if the only component in the neutrino flux were astrophysical in nature, such a low flux would have been too small to explain the low energy events, the contribution of an additional neutrino flux from DM decay, which becomes pertinent precisely at these low energies, would be enough to make up for the deficit. In the process, we have shown that a UHE neutrino flux composed of neutrinos from DM decay, in addition to those from an astrophysical power-law spectrum with a significantly reduced magnitude, is in extremely good agreement with present observations.
With more data and consequently improved statistics in the future, the IC should be in prime position to conclusively determine if it is indeed seeing events from DM decays in addition to those from astrophysical sources or, if, like many analyses based on low statistics made prior to this, the present conclusions may get modified. To determine either way, and to distinguish the signal from a flux component made up of neutrinos from DM decays against a pure astrophysical power-law flux, however, it would be important to note the subtle changes to the shape of the event spectrum that would correspondingly occur. α = 2.0 for the rest of the paper.
[34] Although τDM does not have an observed or theoretically motivated upper bound, the neutrino flux from DM decay falls with increasing decay lifetimes, and when as large as 10 29 s, it already leads to unobservably small event rates at IC. Here, we set the upper bound for the τDM parame-ter space scan to 10 29 s for computational purposes -for the purposes of the analysis, taking even larger values of τDM is equivalent to assuming the neutrino events seen at IC are solely due to the astrophysical power-law flux.