Decaying Leptophilic Dark Matter at IceCube

We present a novel interpretation of IceCube high energy neutrino events (with energy larger than 60 TeV) in terms of an extraterrestrial flux due to two different contributions: a flux originated by known astrophysical sources and dominating IceCube observations up to few hundreds TeV, and a new flux component where the most energetic neutrinos come from the leptophilic three-body decays of dark matter particles with a mass of few PeV. Differently from other approaches, we provide two examples of elementary particle models that do not require extremely tiny coupling constants. We find the compatibility of the theoretical predictions with the IceCube results when the astrophysical flux has a cutoff of the order of 100 TeV (broken power law). In this case the most energetic part of the spectrum (PeV neutrinos) is due to an extra component such as the decay of a very massive dark matter component. Due to the low statistics at our disposal we have considered for simplicity the equivalence between deposited and neutrino energy, however such approximation does not affect dramatically the qualitative results. Of course, a purely astrophysical origin of the neutrino flux (no cutoff in energy below the PeV scale - unbroken power law) is still allowed. If future data will confirm the presence of a sharp cutoff above few PeV this would be in favor of a dark matter interpretation.

energy larger than 60 TeV) in terms of an extraterrestrial flux due to two different contributions: a flux originated by known astrophysical sources and dominating IceCube observations up to few hundreds TeV, and a new flux component where the most energetic neutrinos come from the leptophilic three-body decays of dark matter particles with a mass of few PeV. Differently from other approaches, we provide two examples of elementary particle models that do not require extremely tiny coupling constants. We find the compatibility of the theoretical predictions with the IceCube results when the astrophysical flux has a cutoff of the order of 100 TeV (broken power law). In this case the most energetic part of the spectrum (PeV neutrinos) is due to an extra component such as the decay of a very massive dark matter component. Due to the low statistics at our disposal we have considered for simplicity the equivalence between deposited and neutrino energy, however such approximation does not affect dramatically the qualitative results. Of course, a purely astrophysical origin of the neutrino flux (no cutoff in energy below the PeV scale -unbroken power law) is still allowed. If future data will confirm the presence of a sharp cutoff above few PeV this would be in favor of a dark matter interpretation.

Introduction
After more than 80 years from its first evidence in the Coma galaxy cluster by Fritz Zwicky, the nature of Dark Matter (DM) still remains an open question. Elementary particle physics can provide interesting schemes where to allocate viable DM candidates, one of the most attractive and simple scenarios being the Weakly Interacting Massive Particle (WIMP) scenario with a DM mass in the range O(1) GeV -O(100) TeV [1] 1 and interaction rates of the order of weak interactions. Even though such schemes naturally emerge in the SUSY extension of the electroweak Standard Model (SM), up to now almost all indirect or direct searches have not provided any clear evidence [3], and DM observations remain linked to their indirect gravitational footprint only.
Indeed, a large amount of different theoretical frameworks has been proposed in literature, where the mass of DM candidates is spread over many order of magnitude, from about 10 −32 GeV up to 10 18 GeV, e.g. axions, KeV sterile neutrinos, majorons, or the heaviest wimpzilla (∼ 10 12 GeV). Lacking any direct DM detection at LHC experiments, it is likely that the only viable way to look for very massive DM candidates would exploit indirect searches in astrophysical observations. From this point of view, neutrino telescopes like IceCube (IC) provide a chance to observe high energy cosmic ray phenomena induced by such massive DM particles, where energetic neutrinos are produced.
The IceCube Neutrino Observatory experiment [4,5] is an excellent example of high energy neutrino astronomy. This rapidly developing branch of physics is important for different reasons [6][7][8]. In fact, by looking at neutrinos, which are neutral and weakly interacting particles, one can trace back the sources even in presence of an intergalactic background and magnetic fields. Moreover, the observation of astrophysical neutrinos is of paramount importance, both because their presence is a proof that acceleration of hadronic matter is involved, and because the features in their energy spectrum may indicate their production in a top-down framework from interactions involving heavy non standard particles.
During three years (2010-2013) IC Collaboration has observed several neutrino events in the TeV-PeV range. While the lower energy ones can be well explained in terms of atmospheric neutrinos, events in the range between 100 TeV and 2 PeV seem to be due to some extraterrestrial process. There are many astrophysical possible sources for such neutrinos [9], both from pp [10][11][12] and pγ [13] interactions; the bottom-up scenarios range from extragalactic Supernova Remnants (SNR) [14] to Active Galactic Nuclei (AGN) [15,16] and Gamma Ray Bursts (GRB) [17], all of these with specific emission spectra due to the different production environment. For example, pγ spectrum is peaked [13] while pp spectrum is flatter [12]. According to the IC analysis [5], galactic sources alone can't explain the excess and extragalactic sources are necessary. However, even assuming an extragalactic origin, it is not straightforward to fit all the data; for example, pγ AGN spectrum gives a good description of high energy events but does not satisfactory fit lower energy data [16]. On the other side, GRB predicted spectra agree with data modulo the normalization of flux, which, unfortunately, has an upper limit (given by searches for correlation with observed GRB) more than one order of magnitude below the observed value [18].
More recently it has been proposed that high energy IC events could be related to the decay of DM particles [19][20][21][22][23][24][25][26][27][28][29][30]. Note that the presence of PeV decaying DM component can also alleviate some tension among cosmological parameters estimates [31]. Such an interpretation would be supported by a lack of events above 2 PeV compatible with the decay of a massive particle with a mass of this order of magnitude, and by a distribution of arrival directions apparently uncorrelated with the galactic disk. If this interpretation is correct, IceCube will provide information about the DM mass range and cross-section with ordinary particles and thus, important hints for future DM experiments.
For an elementary particle with a mass of O(1) PeV, the maximal annihilation crosssection obtained by saturating the unitarity bound yields to a completely negligible signal at IceCube. One is then compelled to consider decaying DM instead [19].
Alternative connections between DM and IC events have been proposed in literature. In particular the boosted DM mechanism [32][33][34][35] is based on the idea that a highly energetic (boosted) population of DM particle is originated by the decay of more massive and longlived non-thermal relic, dominating the DM distribution. Such boosted particles then interact with nucleons of detector via neutral current interactions. Finally, different schemes have been proposed, where the bump in the neutrino flux at PeV observed by IceCube is explained as the s-channel enhancement of neutrino-quark scattering by a leptoquark with a mass of 0.6 TeV that couples to the τ -flavour and light quarks [36].
In the present paper, we propose an interpretation of IceCube PeV events in terms of leptophilic three-body decays of a DM particle. The operator responsible for the decay results to have quite a natural coupling and is selected out by global flavour symmetries. Such symmetries forbid at the same time the two body decays and other dangerous three body-decay operators.
The paper is organized as follow. In Section 2 we give a brief review of IC data, and in Section 3 we outline the theoretical framework for DM particle. In Section 4 we discuss the neutrino flux produced, and in Section 5 we describe our results. Section 6 contains our conclusions.

IC data
The IceCube Neutrino Observatory has been searching during three years (2010-2013) [5] for astrophysical neutrinos in the energy range from 30 TeV to 100 PeV. The observed data have been analyzed in terms of their energy spectrum, arrival direction, and flavour. The expected background arises from muons and neutrinos coming from the decays of π and K produced by cosmic ray interactions in the atmosphere. Actually, when the energy of the parent meson increases, its lifetime becomes longer and correspondingly, the interaction probability dominates on the decay, giving a suppression in the atmospheric muon and neutrino flux at high energy. No suppression is instead foreseen for the atmospheric neutrino background coming from the decay of charmed mesons, the so-called prompt component, because of the short lifetime of these mesons. Measuring the muon detection rate in a separate region of the telescope, the IC collaboration gives the following estimation of the muon and neutrino background where N all ν+ν stands for the number of all flavour neutrinos and antineutrinos and its asymmetric error is due to the prompt component.
IceCube collected 37 events in 988 days, with deposited energies ranging from 30 TeV to 2 PeV. In particular, three events were detected, with deposited energy of the order of PeV, which are the most energetic neutrino events ever detected. Among all the events, two of them, events 28 and 32, having sub-threshold signals in IceTop, seem to be part of the expected muon background (in particular, event 32 cannot be reconstructed with a single direction and energy), while three ambiguous downgoing tracks seem to be of an atmospheric origin.
It has been argued in literature [37] that, due to the present uncertainty on the overall rate and starting energy of the prompt neutrino component, by assuming a different cross section for charmed meson production and/or a slightly different cosmic ray primary flux one could modify the expected atmospheric neutrino background in such a way to reduce the necessity of an extra neutrino flux in the high energy range. In this case since about 50% of downgoing prompt neutrinos arrive together with muons, which should trigger the muon veto, upgoing events should be more than downgoing ones at these energies ("southern hemisphere" suppression). IceCube observes exactly the contrary, so it presently seems unlikely that the atmospheric background and the observed neutrino flux could be reconciled thanks to prompt neutrinos [38,39].
Another indication comes from the topology of detected events, which belong to two classes: track events, associated with the propagation of a high energy muon, and shower events, which correspond to the production of a large aggregate of secondary particles (with similar energy resolution in the two cases but, of course, better angular resolution in the first one). By assuming that data are due to a purely conventional atmospheric flux, one should count more tracks than showers (since atmospheric neutrinos are mainly muon neutrinos). On the other side, a (1/3 : 1/3 : 1/3) flavour proportion in the neutrino flux, as expected for astrophysical neutrinos, would result in only 20% ν µ CC interactions, a closer result to the IceCube finding of less tracks (24%) than showers (76%).
On the basis of this analysis the IC collaboration concludes [5] that a purely atmospheric origin of the high energy events, requiring a 3.6 times higher charm normalization, is rejected at 5.7σ. Moreover, the data are well described in terms of a global fit including background atmospheric muons and neutrinos, prompt neutrinos and an isotropic astrophysical flux, which for each flavour takes the form (quoted errors are 1-σ uncertainties) This result satisfies the Waxman-Bahcall bound for optically thin sources [40], obtained supposing that all the charged particles created by cosmic accelerators give their energy to kaons and pions.
While the unbroken E −2 hypothesis gives a fairly good description of the data energy spectrum, it predicts 3.1 more events at 2 PeV, seeming to require a softer spectrum or a high energy cutoff. A more general model of the astrophysical component by a piecewise function of the energy gives the best-fit This corresponds to the lower boundary of the total statistical and systematic uncertainty on the energy power law (with a zero charm contribution), which on the other side reaches, at 90% C.L., the previous E −2 behaviour. In the present work, by following the background analysis of ref. [5], we assume that neutrino flux is mainly dominated by the atmospheric component up to 60 TeV, whereas we consider a bottom-up neutrino contribution from known astrophysical sources (like for example extragalactic SNR) in the [60 TeV, 300 TeV] range (hereafter denoted as astrophysical neutrino flux), and, at the same time, a top-down additional component at higher energy (hereafter denoted DM neutrino flux). This would naturally produce a sharp cutoff, observed in the data around few PeV, and if confirmed by new data, represents an intriguing feature of the IC observations.

The Theoretical Framework
For a heavy fermion singlet DM candidate, χ, which is directly coupled to neutrinos, the lowest dimension coupling is a Yukawa interaction whereφ ≡ iσ 2 φ * (φ denoting the SM Higgs doublet), L α=e,µ,τ stand for the lepton doublets, and σ 2 is a Pauli matrix. As shown in [23], this new interaction term is cosmologically safe (sufficiently long lived DM), and relevant for IceCube (a DM particle specie abundant enough to fit the observed flux) for a fine tuned tiny coupling, y ∼ O(10 −30 ). Note that the interaction term in Eq. (3.1) is reminiscent of the right-handed neutrino in see-saw models, but its contribution to the lightest neutrino mass would be negligible due to the smallness of the y coupling. The operator in Eq. (3.1) would yield a sharp peak in energy. The presence of the Higgs field in Eq. (3.1) allows for an abundant production of secondary neutrinos via the decay of the Higgs particles to heavy quarks, giving an almost flat neutrino flux at energies lower than PeV [23]. Such a result could turn out to be problematic if known astrophysical contributions were included in the analysis. For instance in Ref. [14], low energy (up to 100 TeV) IC events are shown to be well fitted by means of standard extragalactic SNR.
In the present paper we try both to improve the need of an unnatural coupling and at the same time, to reproduce the IC data including an astrophysical neutrino component. This suggests to consider an interaction term that does not directly involve quarks, the Higgs field or gauge bosons (leptophilic DM), and which is a higher dimension operator. In this case, the feeble coupling can be understood in terms of a large mass scale. In particular we assume that: 1) the DM field χ is coupled to the SM particles via a leptophilic coupling; 2) the lifetime of χ is suppressed by powers of the scale of new physics entering in the non-renormalizable coupling; Dimensions DM decay operators Table 1. Gauge-invariant operators up to dimension-6 inducing fermion singlet DM decay as taken from Table 1 of Ref. [41]. The adopted notation is explained in the text.
3) there is a direct coupling to neutrinos allowing for a primary neutrino flux with energy of the order of χ mass. The multi-body final state may contain lower energy neutrinos as well, so that neutrino flux also spreads to lower energies; 4) mass and couplings of χ are determined in order to produce a neutrino flux just dominating the energy region around PeV. This implies, on the contrary, that the flux in the energy region [60 TeV, 300 TeV] is due to astrophysical sources as will be discussed in the following. We also assume for simplicity that χ represents the dominating contribution of Cold DM.
As in [41], one can list the gauge-invariant operators up to dimension-6, shown in Table 1. L and Q are the lepton and quark weak doublets and u and d the corresponding right-handed quarks. The field stands for the right-handed lepton singlet, and finally, W µν and B µν are the field strength tensors of SU (2) L and U (1) Y gauge bosons. To simplify notation, we have omitted the family index for matter fields.
Remarkably, the only operator in this list which satisfies requirements 1)-3) is the nonrenormalizable lepton portal: where {α, β, γ} are flavour indices thus labeling 27 operators, and the round brackets indicate the Lorentz contractions. The mass scale here is chosen as the Planck mass M Pl . Indeed, for this choice the order of magnitude of couplings y will not be unnaturally small, see Section 5. If this operator is the only source of DM decay, one has to invoke some selection rule which forbids the dimension-4 operator, which, as we mentioned, is compatible with IC results for an extremely small coupling only. This can be done by using global flavour symmetries, both Abelian like U f (1) and non-Abelian groups like A 4 , ∆(27), etc. We will focus in the following on two benchmark schemes (see Appendix for more details about these models) 2 : where the brackets show the flavour assignments corresponding to non vanishing couplings y αβγ . Note that expanding the SU (2) contractions, the operator of Eq. (3.2) always yields a coupling of DM with two charged leptons and one neutrino. Depending on the flavour index the charged leptons can then possibly decay producing secondary neutrinos.

DM neutrino flux
Following the method outlined in Ref. [23] we consider both the contributions to top-down neutrino flux coming from the galactic and extragalactic distributions of DM particles χ (the sum of χ andχ contributions is implicitly assumed) where the solid angle integration is on the longitude and latitude in the galactic coordinate system, l and b. In the following we will always sum the flux of neutrinos and antineutrinos of different flavours, unless explicitly stated. The galactic component due to the Milky Way halo can be written as where M χ and τ χ denote the mass and lifetime of DM particle. The quantity ρ χ (r) denotes the density profile of DM particles in our Galaxy as a function of distance from the Galactic center, r, and dN α ν+ν /dE ν stands for the energy spectrum of neutrinos and antineutrinos of flavour α produced in the decay. Note that the parameter s is related to r via the expression r(s, l, b) = s 2 + R 2 − 2sR cos b cos l , where R 8.5 kpc is the distance of the Sun from the Galactic center. In the following we assume a Navarro-Frenk-White density profile for the halo with r c = 20 kpc the critical radius and ρ χ = 0.33 GeV cm −3 . We have checked the dependence of our predictions on the choice for the density profile (Einasto, Isothermal, etc.). The results are very similar, with a variation of the best fit values of the order of few percents. The quantity dN α ν+ν /dE ν has been evaluated by means of a MonteCarlo procedure. In our code we take into account primary neutrinos (antineutrinos) and also the neutrinos (antineutrinos) produced by the decay of µ and τ leptons. The τ leptons have different decay channels that involve pions as well [3]. The MonteCarlo takes into account ∼ 90% of the τ decay width, including the two leptonic decays into muons and electrons, and semileptonic decay channels up to three pions in the final state, whose charged states eventually decay mainly producing µ and corresponding neutrinos (antineutrinos). It is worth observing that the electroweak radiative corrections produce a large number of secondary neutrinos that are typically placed at energies almost two orders of magnitude smaller than the mass of DM particle [42]. Hence, for a mass of few PeV, such low energy neutrinos do not significantly affect our analysis that is restricted to neutrino energies larger than 60 TeV. In our study, we also consider the total energy injected by DM decay in the electromagnetic (e.m.) sector, and compare such a value with the bound provided by Fermi-LAT [43].
The isotropic extragalactic component of the differential flux is given by where H(z) = H 0 Ω Λ + Ω m (1 + z) 3 is the Hubble expansion rate as a function of redshift z and ρ cr = 5.5 × 10 −6 GeV cm −3 is the critical density of the Universe. We assume a ΛCDM cosmology with parameters Ω Λ = 0.6825, Ω m = 0.3175, Ω χ = 0.2685 and h ≡ H 0 /100 km s −1 Mpc −1 = 0.6711, according to Planck experiment results [44]. We found that the galactic and extragalactic components of the DM differential flux are of the same order of magnitude. The χ-lifetime, τ χ , depends on the values assumed for y αβγ . For models 1) and 2) we get where y + and y − denote the two independent symmetric and antisymmetric couplings in A 4 [41]. In the previous expressions we neglected all final state masses, due to the large value of M χ . In the following, for simplicity we assume real couplings and the relations model 1) |y µeτ − y τ eµ | = |y eµe | ≡ y , (4.8) model 2) |y + | = |y − | ≡ y . (4.9)

DM relic abundance
As well known, the standard thermal freeze-out mechanism is not a viable option for particles with masses larger than O(100) TeV because of the unitarity bound on the cross-section [1]. Hence, considering a DM candidate of PeV mass raises the question of the production mechanism, which is crucial to determine its relic abundance. There exist several scenarios, beyond the WIMP paradigm, which make such a heavy particle a perfectly viable DM candidate (see, e.g., Ref. [45] for a review). One possibility is to consider a non-thermal production during reheating of the Universe, after the inflationary stage, in a cosmological scenario with a low reheating temperature [46,47]. After inflation, the Universe reaches its maximal temperature, T max . This temperature can be much higher than the reheating temperature, T RH , defined as where g RH ≡ g (T RH ) is the effective number of relativistic degrees of freedom at T RH , and Γ ρ the inflaton field decay rate. If the χ particle production takes place between T max and T RH , the DM abundance scales as [46]: where σv is the thermal average of the DM annihilation cross-section times the Møller flux factor. Notice that the χ particles are never in chemical equilibrium (even if they are in kinetic equilibrium), so their abundance is described by a power law instead of being exponentially suppressed as exp(−M χ /T RH ). Assuming σv ∼ M −2 χ , we obtain for g RH = 200 and M χ = 1 PeV, if χ provides the whole contribution to DM today. The previous equations are valid under the condition T max > M χ 20 T RH [48], corresponding to the fact that DM never reached local thermal equilibrium in the early Universe (which would otherwise spoil the result of eq. (4.11)). From eq. (4.12) we see that the lower bound is readily satisfied.
Another viable possibility is the production through inelastic scattering interactions among the thermal plasma and energetic particles originating from inflaton decay and taking place at temperatures between T max and T RH [49]. In this case, assuming that DM couplings are of same order as gauge couplings and that M ρ > M 2 χ /2T RH , with M ρ the mass of the inflaton, in order to account for the observed DM abundance the following condition should be satisfied [49]: (4.13) This mechanism would require quite low, thought viable, values for the reheating temperature to account for M χ ∼ 1 PeV. Finally, PeV dark matter could also be produced from inflaton decays (directly or through cascades). However, this mechanism is quite model-dependent because it strongly depends on the couplings of the inflaton to the particles of the model [50,51]. For the sake of illustration, the direct production of DM from inflaton decay gives the abundance (assuming all decays happen at T RH ): where B(ρ → DM ) is the average number of DM particles produced per inflaton decay.

Analysis and results
In order to reproduce the sharp cutoff observed in the IC neutrino data, in the present analysis we assume that the neutrino flux is given by the combination of an atmospheric component up to 60 TeV [5] and two different neutrino contributions, a bottom-up one from known astrophysical sources, like extragalactic SNR, in the [60 TeV, 300 TeV] range (hereafter denoted as J Ast ), and a top-down component at higher energies (that is the flux J χ , defined in Eq. (4.1)). Hence, for E ν ≥ 60 TeV the total neutrino flux is given by Then, the number of neutrinos in a given energy bin [E i , E i+1 ] is equal to where ∆t = 988 days is the exposure time after 3-years of IC experiment, and A α (E) is the neutrino effective area [4] for different neutrino flavour α. In order to compare the theoretical predictions of Eq. (5.2) with the observations we need to recast such expression in terms of deposited energy. In general, to statistically estimate the ratio between the deposited and neutrino energies a MonteCarlo simulation of the apparatus is required [52]. When for a bin of the deposited energy a significant statistics is collected one could apply an average ratio that results to be of the order of (σ CC 97% + σ N C 23%)/(σ CC + σ N C ) ∼ 75% (see Table 1 of [52]).
Remarkably, such a number appears to be quite stable as a function of the neutrino energy. Unfortunately, due to low statistics collected till now, this procedure would be characterized by a large uncertainty in the energy determination. For this reason we prefer in this paper to assume the simplicity ansatz that the two energies coincide. Notice that in any case an expected shift in the energy of the order of 25% is not going to change dramatically our conclusions.
Known astrophysical sources able to produce a neutrino flux similar to the one observed by IceCube are high energy accelerators [9] such as extragalactic SNR [14], Hypernova Remnants (HNR) [14], AGN [15,16], and GRB [17]. It is worth observing that, due to the large uncertainties in the parameters related to the physics of the accelerator mechanisms, there is enough room to reproduce IC neutrino flux. In particular, an extragalactic SNR can produce a neutrino flux with a cutoff O(100) TeV coming from the proton-proton hadronic interactions up to few PeV. On the other hand the HNR are sources of 100 PeV protons that can explain the PeV neutrino events [14]. While classical GRB are ruled out, low-luminosity GRB and choked jets, which cannot be triggered by GRB satellites, are totally viable as the origin of PeV neutrinos [53][54][55]. AGN can be the sources of PeV neutrinos [16,56], but if the flux is normalized to match the observed IC PeV events, the lower energy part is inconsistent with data. All these astrophysical sources have also an associated γ-ray flux that would give a significant contribution to the diffuse γ-ray background, strongly constrained by Fermi-LAT [43].
In order to parametrize the astrophysical flux one can use either an Unbroken Power Law (UPL), with a power law behavior in the whole IC region, or a Broken Power Law (BPL) where an exponential cutoff is assumed at some energy scale E 0 . Considering both options essentially covers the wide range of accelerator mechanisms related to the different astrophysical sources. Hence, using the notation adopted by the IC Collaboration, we will consider: i) Unbroken Power Law (UPL): ii) Broken Power Law (BPL):  Table 2. The marginalized 95% C.L. for parameters y and J 0 (expressed in unit of GeV cm −2 s −1 sr −1 ) corresponding to models 1) and 2), and UPL and BPL parameterizations, respectively. The last column reports the reduced χ 2 .
where γ +2 is the spectral index and J 0 is the flux normalization. In the present analysis we fix the value of E 0 to be equal to 125 TeV in agreement with the extragalactic SNR results [14]. Furthermore, we restricted the spectral index to the physical range γ ∈ [0, 1], as suggested by cosmic accelerator mechanisms. The fit has been done by means of a multi-Poisson likelihood analysis [58] in which the χ 2 takes the expression where N i is the expected number of neutrinos for energy bin provided by Eq. (5.2), while n i is the observed one. Once the atmospheric background has been subtracted, we fit the IC data by using the parametrization of the astrophysical flux in both cases of UPL and BPL, and considering model 1) and model 2) of Section 3 for the DM neutrino flux. The mass of DM particle has been varied in a range [1 PeV,10 PeV] able to produce a drop in the flux. The best fit value is found for M χ = 5.0 PeV independently of the model adopted. On the other side, after scanning the γ range, we find a best fit value γ = 1.0 for UPL (spectral index = 3.0), and γ = 0.0 for BPL (spectral index = 2.0). In Table 2, we give the marginalized 95% C.L. for parameters y and J 0 , for models 1) and 2), and UPL and BPL parameterizations, respectively. From the values of reduced χ 2 shown in the last column, we see that the experimental data slightly prefer the BPL scheme with respect to the UPL one. In each case, the analysis shows that a non-vanishing DM contribution at 2-σ level (95% C.L. for y not compatible with zero) is required. This is mainly due to the presence of the sharp cutoff in the data at high energy. By comparing the results obtained for U f (1) and A 4 , one cannot appreciate a significative difference. Indeed, the two models essentially provide similar features at the level of produced neutrino flux. For the cases reported in Table 2 we have checked that the total energy injected by DM decay in the e.m. sector is smaller than the bound provided by Fermi-LAT [43]. The neutrino events as function of E ν for model 1) are shown in Fig. 1 (first row) and compared with IC data for the two models, DM + UPL (column A) and DM+BPL (column B). In the second row we report the 68% C.L. (dashed) and 95% C.L. (solid) contours for the two parameters, y (defined in Eq. (4.8)) and J 0 (defined in Eq.s (5.3) (5.4)), in the cases DM + UPL (column -10 -  Figure 1. Results of the analysis for model 1). First row shows the neutrino events as a function of the neutrino energy E ν for the DM + UPL (column A) and DM+BPL (column B) models. The red (long-dashed) line is the best fit (background + astrophysical + DM components), and its band represents the 68% C.L. resulting from the fit. The purple (dashed) and green (solid) lines are the astrophysical and DM contributions, respectively. The black points are the IC data, and the blue region shows the upper limit for the sum of all backgrounds (see Ref. [5]). In the second row we report the 68% C.L. (dashed) and 95% C.L. (solid) contours for the two parameters y and J 0 corresponding to DM + UPL (column A) and DM+BPL (column B). The crosses are the best-fit points.
A) and DM+BPL (column B). The crosses stand for the best-fit points. The same quantities, but corresponding to the A 4 model, are reported in Fig. 2. The plots in Figs. 1 and 2 show that there is no significative difference between the two models: both two models predict the observation of neutrinos in the energy range [0.3 PeV, 1.0 PeV] and a sharp cutoff at the energy of few PeV. Notice that since in both U f (1) and A 4 symmetry cases the galactic and extragalactic components of the DM neutrino flux are of the same order of magnitude, we expect at the PeV energies an almost isotropic neutrino flux with a significant level of anisotropy near the galactic center. This is in a good agreement with the IC observations. In order to understand the reason for the strong similarity between the predictions of model 1) and 2), we have also considered two situations where the operator of Eq. (3.2) is characterized by a single term only, {α, β, γ} ≡ {e, e, e} and {α, β, γ} ≡ {τ, τ, τ }, though they cannot be realized by using the flavour symmetries considered here, as shown in the Appendix. From a phenomenological point of view it is interesting to consider such cases because they correspond to the extreme situations with respect to the number of τ leptons produced and thus, of secondary neutrinos from their decays. In the left panel of Fig. 3   report the flavour compositions at Earth of the DM neutrino flux for model 1) and 2) and the two fully diagonal cases (e, e, e) and (τ, τ, τ ), as well. The black dot represents the IC flavour analysis of Ref. [57]. Models 1) and 2) and the (τ, τ, τ ) diagonal case provide the same flavour composition at Earth, namely (f e : f µ : f τ ) ⊕ ≈ (1/3 : 1/3 : 1/3), which are the flavour ratios provided by the standard astrophysical sources. A different situation is given by the diagonal case (e, e, e) where we have only electron (anti)neutrinos as DM decay products, that leads at Earth (f e : f µ : f τ ) ⊕ ≈ (0.55 : 0.19 : 0.26). It is worth reminding that a study of shower/track composition of the PeV events can shed light on the flavour ratios and hence on the flavour structure of the coupling. Such analysis cannot be performed at the moment due to the low statistics available. In the right panel of Fig. 3 one can see that changing the flavour structure of the DM-SM coupling does not appreciably affect the DM neutrino flux.

Conclusions
During the last three years of observation the IceCube Collaboration has observed neutrino events in the TeV-PeV range. Even though lower energy events can be well explained in terms of atmospheric neutrinos, the most energetic processes detected, in the range between 100 TeV and 2 PeV, seem to originate from some extraterrestrial source. In the present analysis, by following the results on the background presented in Ref. [5], we assumed an extraterrestrial neutrino flux dominating for energies larger than 60 TeV. In particular, in order to recover the drop of the flux observed around few PeV, this high energy flux is obtained as the sum of two different components: a bottom-up neutrino contribution coming from known astrophysical sources (like for example extragalactic SNR) in the energy region [60 TeV, 300 TeV], and a top-down additional component, dominating at higher energies, originated from the decay of DM particles with mass of few PeV. The astrophysical flux was parametrized using either an unbroken power law or a broken power law with an exponential cutoff. We have considered both the options UPL and BPL in order to cover the wide range of possible astrophysical sources. The top-down term in the neutrino flux is produced by the decay of a PeV DM particle, χ, which we assume for simplicity to be the dominant cold DM term. The decay of χ has been calculated by means of a MonteCarlo procedure. A similar approach was already investigated in literature [23], but assuming the presence of a trilinear coupling for χ like the one reported in Eq. (3.1). In this case, the χ-decay induces a shower whose hadronic content yields secondary neutrinos that provide an almost flat behavior of neutrino spectrum at low energy. The main problem of such an approach is that it could be in contrast with the standard astrophysical sources since the trilinear interaction term alone totally explains the IC extraterrestrial neutrino flux. Moreover, in this scenario an unnatural tiny coupling for the interaction term of Eq. (3.1), O(10 −30 ), is required. In the present analysis we use the same approach of [23], but assume the presence of special flavour symmetries (Abelian and non-Abelian), like U f (1) and A 4 for example. Such symmetries forbid a trilinear interaction term of the type in Eq. (3.1), and allow quadrilinear leptophilic interaction terms only, like the one in Eq. (3.2). At the same time they reproduce the values of neutrino masses and mixing. The leptophilic characteristic is crucial in forbidding the production of a huge hadronic component in the decay shower that would flatten again the neutrino spectrum. Such a model allows one to improve the problem of naturalness of the coupling, since the non-renormalizable terms is weighted by the square of the new physics mass scale, that can be as large as the Planck mass. In this scenario the coupling required to recover the IC data results to be 25 orders of magnitude larger than the one for the trilinear term.
The comparison of the prediction for the neutrino events in UPL + DM and BPL + DM with the experimental observations allows the determination via a multi-Poisson likelihood approach of the free parameters of the model. In Table 2 we show the resulting parameters for the U f (1) and A 4 symmetry models. Both schemes provide a fair description of the experimental situation within the statistical uncertainty, with the BPL parameterizations providing a slightly smaller value of the reduced χ 2 . The contour plots seem to suggest at 2-σ the presence of a hard component (top-down term) in the spectrum (y not compatible with zero), whereas the astrophysical contribution could be even vanishing (J 0 compatible with zero). This is almost expected because the astrophysical contribution partially overlap with the atmospheric background (affected by a large uncertainty) up to few hundreds TeV. Nevertheless, we cannot conclude at this level of statistics that the presence of DM is really demanding. A confirmation by future data of the cutoff above few PeV is therefore particularly important. Moreover, since the galactic and extragalactic components of the DM neutrino flux are found to be of the same order of magnitude, another evidence in favor of a DM component would be the presence of a significant level of anisotropy near the Galactic Center.
We also analyzed the flavour ratios of neutrinos at Earth predicted by the two models, and found that for both models, (f e : f µ : f τ ) ⊕ ≈ (1/3 : 1/3 : 1/3). Only a coupling of χ with first lepton generation only would produce clear features in the flavour ratios, which in principle could be measurable if a large statistics of events were available.
Finally, we would like to highlight that our analysis regards only the three years IC data. The recent new IC data provide a track event with a deposited energy of about 2.6 PeV, that could be explained by increasing the mass of DM particle in our model. However, since such an event is not fully contained, the energy of the primary neutrino is unknown and could also be very high (O(10) PeV). Note that if the energy of the primary neutrino was of the order of few PeV this would not change our conclusions, differently a neutrino energy of the order of 10 PeV would be in tension with our analysis. L e , e L µ , µ L τ , τ φ χ q ψ 2 1 4 0 3 L φ χ A 4 3 3 1 1 Table 3. Charges, q ψ , for the U f (1) symmetry (upper Table). In the lower Table we report the Irreducible Representations allocating SM fields and the dark matter χ in the non-Abelian A 4 symmetric model (quarks are all singlets under A 4 ).
fact, in this case both the operatorsL α αLα χ andL α φ α should be invariant under such a symmetry, but this would imply the invariance of the term L αφ χ as well. Indeed, in terms of abelian charges we would have −2q Lα + q α + q χ = 0 and −q Lα + q φ + q α = 0, which implies −q Lα − q φ + q χ = 0. 3 However, if we mix different lepton flavours then we have more freedom to consistently define the charges of a U f (1) flavour symmetry. This can be done in such a way that only a leptophilic dimension 6 operator of Table 1 results to be invariant. A possible realization of such a scheme is shown in Table 3. In this case the only invariant DM decay operators are Pl y µeτ L µ e L τ χ + y τ eµ L τ e L µ χ + y eµe L e µ L e χ + h.c. . (6.1) The charge assignment in Table 3 is by no means unique; different choices would select different operators. Moreover with the charge assignment given in Table 3, the dimension five Weinberg operator L α L c βφφ has a structure very similar to the so-called B 4 two-zeros texture given in [59] that can fit the lepton mixing parameters, see for instance [60,61].
In principle the form of the operator is dictated by experiment. A hypothetical study of the flavour composition of PeV neutrinos in IC data would give useful hints in the definition of the possible flavour charges.

C Non-Abelian case
Another possible realization of our scheme invokes non-Abelian discrete symmetries. A simple model would employ for instance ∆(27) by assigning L, φ and to 3 Irreducible Representation while χ remains blind to the symmetry. However, for the sake of definiteness we will consider in the rest of this article a framework based on A 4 flavour symmetry with L and transforming as 3 while all the remaining fields are singlets (for more details about such a model we refer to [41]). Our conclusions are not affected by this choice. The Lagrangian of the model contains the following relevant terms: where L dim=6 is the lowest order non-renormalizable operator allowed by the symmetries of the model. It represents a Fermi-like decay interactions give by where the notation (..) 3(3 ) reflects the fact that in A 4 there are two possible contractions of two triplets into one triplet representation. We see that Eq. (6.3) not only reproduces Eq. (3.2) but it also significantly simplifies its structure. Note that in this case we have two independent couplings, namely y and y .