Heavy right-handed neutrino dark matter and PeV neutrinos at IceCube

We discuss a simple non-supersymmetric model based on the electroweak gauge group $SU(2)_L\times SU(2)^\prime\times U(1)_{B-L}$ where the lightest of the right-handed neutrinos, which are part of the leptonic doublet of $SU(2)^\prime$, play the role of a long-lived unstable dark matter with mass in the multi-PeV range. We use a resonant $s$-channel annihilation to obtain the correct thermal relic density and relax the unitarity bound on dark matter mass. In this model, there exists a 3-body dark matter decay mode producing tau leptons and neutrinos, which could be the source for the PeV cascade events observed in the IceCube experiment. The model can be tested with more precise flavor information of the highest-energy neutrino events in future data.


Introduction
The IceCube neutrino telescope has reported 54 ultra-high energy (UHE) neutrino events with deposited energies ranging all the way from 20 TeV to 2 PeV in its 4-year dataset [1][2][3]. This constitutes a 6.4σ excess over the expected background of atmospheric muons and neutrinos. The events seem to be isotropically distributed in the sky, with no statistically significant evidence of point-like sources [4,5]. The specific origin, spectral shape and flavor composition of these events are currently unknown, but understanding all the observed features so far using conventional astrophysics alone appears to be a challenge [6][7][8]. This has inspired many speculations regarding possible new physics beyond the standard model (BSM) of particle physics that could (partly) explain the origin of these events [6,9], although the current statistics does not necessarily call for a BSM interpretation yet [10][11][12][13]. Nevertheless, if some peculiar features in the IceCube data, namely, an apparent energy gap just below PeV and a slight excess of events above PeV over the predictions from a single, unbroken power-law astrophysical neutrino flux, as well as the lack of any events near the Glashow resonance of 6.3 PeV [14] and no statistically significant correlation of the neutrino arrival directions with the galactic disk, persist with more data, it might lend support to a simple BSM interpretation in terms of a long lived supermassive (multi-PeV) particle dark matter (DM) decaying into neutrinos [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31], which is assumed to be the case in this paper. With the conventional direct and indirect detection, as well as collider searches for weakly interacting massive particle (WIMP) DM candidates in the O(GeV-TeV) range being unsuccessful so far, it is worthwhile to consider the observational prospects of non-WIMP scenarios, such as the PeV-scale decaying DM alluded to above, for which IceCube offers a unique opportunity. Another general motivation for the presence of a decaying DM component, though not necessarily applicable to our model, is due to the fact that it can possibly alleviate the tension between Planck data and low redshift astronomical measurements [32,33], although the extent of improvement is quite modest [34], requiring lifetimes comparable with or smaller than the age of the universe and possibly involving only a fraction of the DM. A PeV-scale RH neutrino DM can also be accommodated in leptogenesis models for explaining the matter-antimatter asymmetry [35,36].
Since this phenomenon involves neutrinos, it is plausible to surmise that it is related to neutrino mass physics. We take this approach here and present a ultraviolet (UV) complete model, where a heavy right-handed (RH) neutrino, which could be at the root of nonzero neutrino masses, is cosmologically stable enough to play the role of DM in our Universe and, being unstable, produces the energetic neutrinos when it decays. This model is however very different from the usual type-I seesaw models [37][38][39][40][41] for neutrino masses, where the heavy RH neutrinos decay rapidly to the SM Higgs boson and known leptons through their coupling with the light neutrinos, and therefore, cannot qualify as the DM of the Universe, unless the coupling constant is extremely small ∼ O(10 −30 ) [20,24,30] in which case the sub-eV light neutrino masses cannot be generated from type-I seesaw in the first place. Instead we contemplate that the RH neutrino N has its own SU (2) gauge interactions (from coupling to the corresponding gauge fields W , Z ) and remains secluded from the SM sector by a softly broken discrete Z 2 symmetry so that its decays to lepton plus Higgs are forbidden. It couples only superweakly to the SM sector via tiny scalar, fermion and gauge (W − W and Z − Z ) mixings. These mixings, whose smallness is guaranteed to be natural due to presence of the softly broken Z 2 symmetry as well as the large hierarchy between the electroweak (EW) scale and the PeV scale, provide a bridge between the secluded heavy fermion sector and the SM sector. They also explain the ultra-long lifetime for the DM N . This is the main new result of this paper. We then study some implications of this particular UV complete model for the IceCube neutrinos and also discuss the constraints from diffuse gamma ray emission.
This paper is organized as follows: in Section 2 we provide an outline of the model. In Section 3 we show how the relic density of the DM N arises in this model. Section 4 discusses its various decay modes and how its long lifetime arises. In Section 5, we illustrate how it explains the observed rate for the multi-TeV to PeV neutrinos at IceCube. In Section 6 we comment on the diffuse gamma ray spectrum associated with this decay. We then conclude with a summary of the results in Section 7.

The Model
The model is based on the gauge group SU Heavy SU (2) singlets : U L : 3, 1, 1, We assume that under the Z 2 symmetry the heavy multiplets Q R , Ψ R , U L , D L , E L are odd and the SM fermions are even. As a result, the SM and heavy fermions do not form masses of type U L u R in the Z 2 symmetry limit [42][43][44][45], which is crucial to accommodate a long-lived RH neutrino DM candidate in our model. On the other hand, both sectors share common SU (3) c and U (1) B−L symmetries, with obvious charge assignments, and the electric charge formula is given as in the left-right (LR) symmetric models by [46,47] In this sense, our model duplicates the SM fields except for the common SU (3) c and U (1) gauge interactions and at this stage similar to the model in Refs. [42][43][44][45].
The minimal Higgs sector of the model consists of SU (2) and SU (2) doublets, respectively denoted as which are both lepton-specific in the SM and heavy sector, analogous to the lepton-specific two Higgs doublet model (2HDM) [48][49][50][51]. The Higgs doublets are chosen to be even under the Z 2 symmetry. In the symmetry limit, the Yukawa couplings are given by the Lagrangian whereχ q = iσ 2 χ * q and similarly for χ q , with σ 2 being the second Pauli matrix. We give different vacuum expectation values (VEVs) to the doublets χ a and χ a , i.e.
with v 2 + v 2 q ≡ v EW 174 GeV and v 2 + v 2 q ≡ v ∼ O(10 PeV) to accommodate a superheavy DM for explaining the IceCube neutrino events. The VEVs v ,q are responsible for the SM charged fermion masses as usual, whereas v ,q make the new charged fermions superheavy. When the ratio of the VEVs tan β = χ 0 q / χ 0 1, the charged scalar χ ± decays predominantly into SM leptons. The reason is as follows: as in the most general lepton-specific 2HDMs, the physical charged Higgs boson is a linear combination of the doublets of form (v q χ ± − v χ ± q ). When the leptonic VEV v is much smaller than v q , the leptonic Yukawa couplings are enhanced by y v q /y q v . This is crucial to accommodate a leptophilic DM candidate in order to fit the IceCube data [29].
It is important to point out that this model is different from the conventional LR models with vector-like fermions [52][53][54][55][56][57][58] where the parity symmetry (also a discrete Z 2 symmetry) makes the Yukawa couplings equal in the left-and right-handed sectors. In our model instead, we have a different Z 2 symmetry which restricts the nature of Yukawa couplings but not the flavor structure of these couplings. To point out other differences between the two Z 2 's, • As mentioned above, the parity Z 2 symmetry allow couplings among the SM fermions and the heavy vector-like fermions, whereas in our case, in the symmetry limit, both sectors are distinct and do not mix with one another.
• In the parity Z 2 case, the SM fermion masses are given by a seesaw formula which depends quadratically on Yukawa couplings in the Lagrangian, whereas in our case, the Yukawa coupling hierarchy depends largely on the VEV ratio tan β as in general lepton-specific 2HDMs [48][49][50][51].
To understand the neutrino masses in the model, we add two SU (2) triplets with a Higgs potential of the form: with a = , q. In the above expression, we have omitted the quartic terms of the form (χ † a χ a ) 2 etc and shown only the terms relevant for heavy and light neutrino masses after spontaneous symmetry breaking by the triplet VEVs We choose the parameters of the model such that v L ∼ eV (corresponding to the soft mass parameter m ∼ 1 GeV), and v R ∼ 10 PeV (corresponding to the mass parameters m , M, M ∼ 10 PeV). Given the Yukawa interactions the ∆ term gives masses to the RH neutrinos of order of few PeV, whereas v L gives masses to the left-handed neutrinos via the usual type-II seesaw mechanism [59][60][61][62][63].
An important feature of this model is that for λ χχ = 0, there is no χ ± a − χ ± a , ∆ ± − ∆ ± or W − W mixing at the tree level (even at 1-loop level for the last two). As a result, if the lightest RH neutrino N has mass lower than the other fermions and bosons of the SU (2)sector, it will be stable in the Z 2 -symmetric limit. However, when we add small soft breaking terms to the model of the form a small W − W mixing can be generated at 1-loop level with a magnitude [64,65] with g and g the SM SU (2) L and SU (2) gauge couplings respectively, m t, b the masses of SM top and bottom quarks, and M T, B the masses of heavy top and bottom partner fermions. Since δ U, D are small soft breaking terms, the induced W − W mixing is small, e.g. ζ W W ∼ 10 −24 δ U δ D GeV −2 for PeV-scale SU (2) -breaking. Similarly ∆ ± −∆ ± mixing is also a loop effect and is expected to be of similar order. The dominant Higgs mixing connecting the heavy sector to the light is the χ ± a − χ ± a mixing which arises at tree level when λ χχ = 0 in Eq. (2.8). This mixing is of order In the fermion sector, there are also mixings induced by the δ term in Eq. (2.11). This mixing only connects the heavy and light charged leptons and is given by ζ eE ∼ δ /v . We will see in the subsequent section that both the scalar and fermion mixings are essential to allow a long lifetime for the lightest heavy neutrino N , as would be required for understanding the IceCube PeV neutrinos.
Since there is a clear separation of scales in our model, with only the SM spectrum (in a 2HDM extension) in the infra-red, and all the rest of the spectrum at or above the PeV scale, the two being linked by portal-like interactions, it is instructive to sketch an effective theory from the low-energy point of view and justify the small number coefficients of the phenomenologically important operators. At the effective theory level, our model has the following features added to the SM: the RH neutrino DM N , which is the lightest of the three RH neutrinos in the UV-complete theory presented above and its interactions with SM fields given by where Λ is of order of the mass of the heavy Higgs boson. We have chosen the RH neutrino DM to have only leptophilic interactions to be in agreement with data. Also we note that the dimension-5 neutrino mass operator in Eq. (2.13) arises from a type-II seesaw with a heavy triplet scalar at the PeV scale which has been integrated out in the low-energy effective Lagrangian.

DM Relic density
For generic thermal relic DM models, there is a generic upper limit on the DM mass from the unitarity limit on the annihilation cross section [66]. However, as noted already in [66] and explicitly demonstrated in this section, the unitarity bound of O(100) TeV can be relaxed in the case of a resonant annihilation, where the N N annihilation cross section can have a Breit-Wigner enhancement [67]. In order to determine the relic density of DM, we note that in the early universe, all the heavy particles were in equilibrium with the light SM sector particles due to the SU (3) c and U (1) B−L gauge interactions. As the universe cools, the particles of the heavy sector being heavier than the DM N , slowly annihilate away leaving the N 's in the primordial plasma. As the temperature falls below M N , the DM density gets Boltzmann-suppressed by e −M N /T . The primary annihilation channels to SM particles proceed via particles that connect the two sectors such as the neutral gauge bosons Z and Z , which mix at the tree level, and the scalar portal mediated by the λ(χ † a χ a )Tr ∆ † ∆ interaction term in the scalar potential. It is easy to see that since Z − Z mixing angle is highly suppressed by the VEV ratio v 2 EW /v 2 [56], its contribution to DM annihilation is very small and the Higgs portal dominates, which we consider below.
where we have neglected the χ a − χ a mixing which is suppressed by v EW /v , h is the SM Higgs, and φ ± and φ 0 respectively are the longitudinal mode of the SM W and Z bosons. In our case, at the EW scale we have only the four light states above in the scalar sector and all other components are at the PeV scale. 1 We have also neglected the tiny W − W mixing which arises at 1-loop level in presence of the soft breaking terms in Eq. (2.11).
The annihilation processes of interest are (see Figure 1) It turns out that the thermally averaged annihilation cross section times velocity is where we have used the DM mass M N = 2f v R , and is the resonance enhancement factor. If R ∼ O(1), then there would generally be the unitarity problem with the PeV scale DM, and this cross section is not large enough to reduce the relic density of N 's to the desired level. Thus we need a large enhancement factor of R which happens for M ∆ 0 2M N . For the sake of simplicity, we assume all other heavy products such as W W are kinematically forbidden and the exact evaluation of R involves only the decay width The relic density for the standard thermal relic reads [68] with the Planck mass M Pl = 1.22 × 10 19 GeV, x F = M N /T F 20 with T F being the freeze-out temperature, g * = 106.75 the relativistic degrees of freedom at T F , a and b the coefficients in the Taylor expansion σv = a + b v 2 + O(v 4 ). We consider first a simplified case, where the mediator scalar mass M ∆ 0 is very close to but slightly lighter than 2M N , then ∆ 0 could decay only into the SM Higgs and the longitudinal W and Z components.
and the relic density in Eq. (3.7) reads From Eq. (3.9), it is clear that we can always obtain the right relic density by appropriately choosing the four relevant model parameters λ, v R , M N and M ∆ 0 . As an illustration, the relic density for a benchmark point in the model parameter space is presented in Figure 2 as a function of the deviation from resonance given by M N − M ∆ 0 /2. Here we have set v = M ∆ 0 = 8 PeV and have calculated Ω N h 2 for different values of the quartic coupling λ. The horizontal dashed line shows the observed relic density, as measured by Planck [69].
As mentioned above, a fine tuning is required for the DM mass M N (or the mediator mass M ∆ 0 ), i.e. |M ∆ 0 − 2M N | < 0.5 GeV, whereas the quartic coupling λ also needs to be small 10 −3.5 in order to reproduce the correct thermal relic abundance. For larger values of λ, the relic density at the resonance increases, as is evident from Eq. (3.9). For smaller values of λ, on the other hand, the correct relic density can only be achieved very close to the resonance, which becomes narrower due to the smaller decay width [cf. Eq. (3.5)].
We also point out that in the parameter region away from resonance, where the primordial DM density is higher than desired, the correct value can be obtained by late decay of second lightest RH neutrino to relativistic SM fermions and resulting entropy generation which can cause dilution. We defer discussing the details of the mechanism to a forthcoming paper [70].

DM decay
In this model, two key pieces of information are important to understand the decay of DM: • All the new particles in the heavy sector are heavier than the RH neutrino DM N , which can be achieved by tuning properly the gauge, scalar and Yukawa parameters in the heavy sector. • In the limit of exact Z 2 symmetry, interactions between the heavy and light sectors involve two fields from the same sector and therefore in that limit, the N 's can only annihilate in pairs but not decay. This is very similar to R-parity in supersymmetry [71].
For N to decay to SM fields, we need to invoke soft breaking of Z 2 symmetry which can give rise to mixings between W − W , ∆ ± − ∆ ± and χ ± a − χ ± a . It turns out that the W − W and ∆ ± − ∆ ± mixings are forbidden at the tree level and arise only at loop levels, and are therefore suppressed compared to χ ± a − χ ± a mixing which can arise at the tree level due to the λ χχ term in the scalar potential given by Eq. (2.8). In addition, the small Z 2 breaking terms can also induce a E L − e R mixing. All these facts then provide a link between the DM N and the SM sector, so that N can decay into light neutrinos, thus giving a potential signal at IceCube. The lifetime of N is resultantly governed by the soft breaking parameters which can be adjusted to be reasonably small to make the lifetime τ N of N much longer than the age of the Universe.
As for the lepton sector, the Z 2 conserving and soft breaking terms relevant to the DM decay can be read off from Eqs. (2.5) and (2.11): After spontaneous symmetry breaking of the SU (2) and SU (2) L gauge symmetries at the PeV and EW scales respectively, we obtain the charged lepton mixing ζ eE ∼ δ /v . As for the charged scalar mixing in the doublet sector, i.e. χ ± a − χ ± a mixing, after symmetry breaking and applying the minimization conditions, we obtain the charged scalar mass terms: Regarding the DM decay, the two charged scalars χ and χ couple predominantly to the heavy and SM sector respectively and have a mixing of order λ χχ v EW /v , which gives rise to a 3-body decay of the DM into a light neutrino plus two SM charged leptons: N → − + ν , as shown in Figure 3, with the crosses denoting the (heavy-light) scalar and fermion mixings. Actually, one of the prompt charged leptons is produced via its mixing with the heavy charged leptons E and its flavor depends largely on the texture of y and δ . For a large tan β, the final states are mostly of τ -lepton flavor, as in the case of lepton-specific 2HDM [48][49][50][51].
The decay of any DM candidate injects energy into the intergalactic and interstellar medium in the form of quarks, leptons, photons or neutrinos, which has potential effects on a large number of cosmological and astrophysical observables [72,73]. For instance, it can delay recombination and/or contribute to reionization, leading to distortions in the cosmic microwave background (CMB). A recent analysis of the cosmic reionization due to DM decay using the Planck data gives an almost model-independent lower bound on the DM life-time of ∼ 10 26 sec [74], much larger than the actual age of the Universe τ U ∼ 4 × 10 17 sec. Similar model-independent limits were also obtained using the neutrino flux limits [75,76]. The soft breaking parameters in 3-body decay of N in our model help to push the DM lifetime to be much longer than the cosmological scale to avoid these constraints. A rough estimation of the lifetime reads (4.4) Thus, with a small symmetry breaking parameter δ MeV in Eq. (2.11), we can satisfy the cosmological constraints on our decaying DM scenario. While this is a tiny parameter, since it is a soft breaking of the discrete Z 2 symmetry, it is stable under renormalization and is therefore technically natural.

Fitting the IceCube Data
Before getting into details, we give an outline of the argument that shows how this model fits the observed event rate of the UHE neutrinos on Earth from the decay of N , as discussed above. Our arguments are very similar to the phenomenological implications of a generic unstable leptophilic DM [29]. The DM contribution to energy density in the Universe is roughly 25% of the critical density ρ c = 5.5 keV cm −3 ; so the density of a 4 PeV DM is roughly n DM ∼ 10 −12 cm −3 . If we assume that its lifetime is τ DM ∼ 10 27 sec, then the probability for each DM to decay is τ U /τ DM ∼ 10 −10 . Multiplying it by n DM , we get the number density of neutrinos from DM decay to be about n ν ∼ 10 −22 cm −3 . To get the flux of neutrinos per steradian on earth, we multiply n ν by the velocity of neutrinos v ν ∼ c (where c = 3 × 10 10 cm sec −1 is the speed of light) giving which agrees roughly with the flux required to fit the IceCube events at E ∼ 1 PeV [3]. Note that for an astrophysical E −2 flux as predicted by the Fermi shock accelaration mechanism, the required flux is at the edge of the Waxman-Bahcall bound for optically thin sources [79]. So invoking the decaying DM scenario for PeV events mitigates the situation to some extent. In order to do a detailed fitting of the IceCube data, we need the energy distribution of neutrinos in the 3-body decay of the RH neutrino N → τ + τ − ν τ . The neutrino energy distribution is similar to the case of electron energy in the muon decay (assuming massless final states) and is given by In Figure 4, we show the energy spectrum as a function of the true neutrino energy for various representative values of M N . From this, we can infer that the PeV neutrino events at IceCube can be explained by our DM decay scenario with the DM mass of about 4 PeV (see Figure 5 below). In practice, the τ -leptons also decay giving rise to secondary neutrinos, as well as electrons, muons and hadrons, which will subsequently lead to more secondary neutrinos in the IceCube detector, thereby raising the tail-end of the spectrum shown in Figure 4. In this respect, the exact neutrino spectrum seen at IceCube will crucially depend on the final-state flavors in the 3-body DM decay, and this feature can in principle be used to probe the flavor structure of the model in future data with more statistics. It should be noted here that for a larger W − W mixing at 1-loop level induced by the δ U,D terms in Eq. (2.11), the 2-body decay modes N → W ± ∓ could also be significant, with the W boson decaying further into SM hadrons or leptons. The large hadronic branching ratio of W gives rise to abundant secondary neutrino emissions from the quarks, which will lead to an almost flat neutrino flux at energies below PeV [16]. However, it turns out to be problematic, if a known astrophysical contribution is also added, which by itself provides a good fit to the lower-energy data, and any additional contribution in the lower energy bins seems to be disfavored by the 4-year IceCube data [3]. Therefore, although in the 2-body DM decay scenario the scalar sector can  be made much simpler, we will not consider this case in our analysis and mainly focus on the 3-body decay mentioned above.
To calculate the neutrino flux in our DM decay scenario, we follow the general method outlined in Refs. [24,25,27,29] and consider both galactic (G) and extragalactic (EG) DM components: where Ω is the solid angle and l, b are the longitude and latitude in the galactic coordinate system, respectively. The galactic component can be explicitly written as where the distance parameter s is related to the radial distance from the galactic center r by r(s, l, b) = s 2 + R 2 − 2sR cos b cos l, R 8.5 kpc is the distance of the Sun from the galactic center, and ρ N (r) is the DM density profile in Milky Way, for which we assume a Navarro-Frenk-White profile [77]: ρ N (r) = ρ 0 (r 0 /r)/(1 + r/r 0 ) 2 , with r 0 = 20 kpc and ρ 0 = 0.33 GeV cm −3 . 2 The quantity dN/dE ν has been computed using the numerical methods outlined in Refs. [72,80], which includes the primary neutrinos and antineutrinos from the DM decay, as well as the secondary ones produced by the τ -lepton decays (including the secondary pion decays). We have implicitly assumed the sum over all neutrino flavors, whereever applicable.
For the isotropic extragalactic component of the differential flux, we have where H(z) = H 0 Ω Λ + Ω m (1 + z) 3 is the Hubble rate as a function of the redshift z, H 0 = 67 km sec −1 Mpc −1 , ρ EG N = Ω DM ρ c , and we assume a ΛCDM cosmology with Ω Λ = 0.68, Ω m = 0.32, Ω DM = 0.27 from the Planck data [69]. From Eqs. (5.4) and (5.5), we note that the neutrino flux is inversely proportional to the product of the DM particle mass and lifetime. Thus for a fixed lifetime, the flux is inversely proportional to the DM mass due to the lower number density of DM particles.
We also include the standard pion-decay contribution to the flux of astrophysical neutrinos, which could come from known sources like active galactic nuclei [81,82] or supernova remnants [83]. Since an astrophysical component almost certainly exists and fits the lowerenergy part of the IceCube UHE event distribution quite well, it should not be outrightly discarded in favor of an entirely new physics interpretation. As an illustration, we assume a single unbroken power-law astrophysical flux: where Φ 0 = 2.2 GeV cm −2 sec −1 sr −1 and γ = 0.58 corresponding to the central value of the IceCube best-fit [3], assuming (1 : 1 : 1) flavor composition on Earth. 3 On the other hand, the DM decay in our model produces mostly τ neutrinos which, after oscillations over astronomical distances, average out to give a flavor ratio of roughly (4 : 7 : 7) on Earth [84]. The total neutrino flux is given by the sum of the DM decay and astrophysical contributions: Using this flux and following the analysis method outlined in Refs. [11,12,85], we compute the number of neutrino events in a given deposited energy bin at IceCube: where T is the exposure time, E dep (E ν ) is the electromagnetic (EM)-equivalent deposited energy for a given incoming neutrino energy E ν in the laboratory frame, A is the neutrino effective area (for a given flavor), and we have summed over all the neutrino flavors and integrated over the whole sky. Our results are shown in Figure 5 for a typical benchmark point with M N = 4 PeV and τ N = 10 28 sec. The background from atmospheric muons and neutrinos, the IceCube data points and the SM best-fit solution (including both charged and neutral current events) are taken from the 4-year IceCube analysis of Ref. [3]. 4 Figure 5 illustrates the fact that while the low-energy events can be readily explained by an astrophysical component of the neutrino flux, the apparent excess just above PeV energy and the subsequent sharp cut-off can be better understood by invoking a decaying PeV-scale DM hypothesis. Note that our decaying DM scenario with ν τ final states will mostly produce hadronic showers near the high-energy cut-off, as required to fit the current data. Our fit gives a slight excess in the bin just below PeV, which we believe is still consistent with the IceCube observations, since there are a bunch of lower-energy throughgoing track events, whose true energy could easily be large enough to fill the gap below PeV, given the fact that IceCube can only put a lower limit on the throughgoing muons. From Figure 5, we conclude that although not statistical significant yet, the spectral features in the future IceCube data could be exploited to test the existence of a superheavy DM in our Universe. In particular, the flavor information could be useful to distinguish our scenario (which predicts mostly τ -flavor final states) from other decaying DM models. A characterisitc feature of the tau-events is the 'double-bang' signature; however, for the current IceCube string separation of 120 m, this feature can be seen only for events with more than 5 PeV energy, which is not kinematically possible for our decaying DM scenario considered in Figure 5. Nevertheless, our scenario predicts 17% more track-events than the SM expectation for a (1 : 1 : 1) flavor composition on Earth, and therefore, might be used as a distinguishing feature in future when more data with more accurate flavor information is available. Finally, it should be remarked that the present IceCube is not large enough to test the decay lifetime of 10 28 sec, but there exist other critical multi-messenger tests that are feasible with current and near-future γ-ray detectors [28].
Apart from the excess at PeV-scale, the 4-year IceCube data as shown in Figure 5 also suggests an apparent excess around 100 TeV. If this becomes statistically more significant with respect to the excess at PeV scale and the apparent energy gap just below the PeV bin disappears, this can in principle be accommodated in our decaying DM scenario by simply choosing a lower value for the dark matter mass in the 300-400 TeV range (along with an appropriate mass spectrum for the other heavy states in the model) and a flatter astrophysical background (e.g. E −2 ). This has already been noted in the literature [29], and we do not repeat this analysis here. If more than one excess surface at very different energy scales, it will be difficult to accommodate all of them in a single DM decay scenario, and one might need to invoke a multi-component decaying DM model.

Gamma ray background from DM decay
One of the features of our model is that the neutrinos will be accompanied by energetic gamma rays from the associated secondary positron production and their subsequent annihilation as they propagate through the intergalactic space, as well as from the final state π 0 's produced in τ decay. Naively one would expect that the TeV gamma rays arising in this process will have a flux which is comparable to the neutrino flux i.e. E γ Φ γ ∼ 10 −7 GeV cm −2 sec −1 sr −1 . (6.1) Note that this corresponds to an absolute flux of Φ γ ∼ 10 −13 cm −2 sec −1 sr −1 for a PeV gamma ray. This seems to be just within the current bound on diffuse gamma ray flux at this energy [87][88][89][90][91][92]. However, one would also expect lower energy gamma rays as the DM decay electrons cool and lose energy before reaching the Earth. We argue below that those lower energy gamma rays are also well below the current upper limits provided by Fermi-LAT [93], HESS [94], HAWC [95] and VERITAS [96]. First we note that the electrons will produce gamma rays of energy E γ m DM /3 at an energy flux twice that of neutrinos. One then has to compute the spectrum of resulting photons from the spectrum of injected electrons. If the electrons did not lose energy sufficiently fast compared to their injection time (the age of the Universe) all electrons would be clustered around E e m DM /3 with their normalization being N e ∼ Q×τ U , with Q being the injection rate; then to compute the photon emission rate we would simply multiply N e σ T c, where σ T is the Thomson cross section.
However the electrons of PeV energy lose energy fast, mainly on the CMB, barely at the Thomson limit (PeV × E CMB m e c 2 ). Therefore one has to compute the electron distribution function by dividing the injection rate Q(γ e ) by their energy loss time scale. The rate of energy loss of electrons of Lorentz factor γ on photons of energy density ρ CMB is given byγ with corresponding time scale τ e γ e /γ e 10 20 /γ sec 3, 000 years. So the highest energy electrons lose their energy to photons of similar energy on time scale of 3,000 years. Therefore the amount of PeV gamma rays produced at the source will be the same as the rate of DM decay into neutrinos and electrons. If we fix τ N 10 28 sec to provide neutrino rate in agreement with IceCube data, the photon production rate at PeV energy will be similar.
However, we would like to compute the entire spectrum of the cooling electrons to check if there are discrepancies at other lower energies. The continuous injection of electrons will create differential power law spectrum of slope −2 (becauseγ ∝ γ 2 ) Inverse Compton scattering on these electrons, upon integration of this distribution with the Compton cross section which has approximately the form dσ/dγ ∼ σ T δ(E γ − γ 2 e ) (with being the soft photon energy that gets up Comptonized), will produce a spectrum of the following form and then E 2 γ (dΦ γ /dE γ ) ∝ E 1/2 γ GeV cm −2 sec −1 with maximum energy m DM /3. The flux at maximum energy will then be twice that of the neutrinos of that energy. However, because the spectrum decreases with decreasing energy it will not have an impact on the diffuse gamma-ray background (DGRB).
Finally, these gamma rays of the highest energy will be absorbed over distances short compared to the Hubble radius and re-inject electrons to form pairs and a cascade such as that in Refs. [97,98] . One has to include these in the calculations. The final outcome of such cascades is to produce a spectrum of photons of slope close to −2 (and E 2 γ (dΦ γ /dE γ ) ∝ E 0 γ ) down to a critical photon energy E γ,c at which space (at z 0) is transparent to the photonphoton scattering process, i.e. one at which photon-photon opacity τ γγ (E γ, c ) = 1, and a spectrum ∝ E , still below the level of DGRB, where E γ,c 20 TeV [99]. Therefore, the observed DGRB at energies ∼ 1 − 50 GeV is much too bright to be affected by these photons, if the neutrinos from this process are limited not to exceed the observed limits.
Before closing, we also note that a sizable fraction of the photons is produced in the galactic halo, just like for the neutrino flux, and is affected by a more complicated (and partial) absorption/reprocessing, in the O(100-1000) TeV energy regime. This is typically probed by extensive air shower type of detectors, and could manifest in peculiar anisotropy signals. For a detailed discussion, see e.g. Ref. [92].

Comments and Conclusion
A few comments are in order regarding the model and its implications: • If the gauge symmetry SU (2) is identified as the right-handed gauge group SU (2) R of the LR symmetric models, then it will be a variation of the conventional LR model with heavy vector-like fermions [52][53][54][55][56][57][58], with the heavy sector at a higher energy scale. However to have a viable dark matter in this case, we need to increase the scale to much higher values to avoid the 3-body N decay to off-shell W R mediated by SU (2) R gauge ineraction and to adjust the RH neutrino Yukawa couplings f appropriately so that the lightest N is lighter than all other fermions of the heavy sector.
• For the model to work, two fine-tunings are needed: (i) the usual SM Higgs fine-tuning since we have new physics at a PeV scale, and (ii) tuning of the DM mass 2M N M ∆ to get the right relic density. The second fine-tuning could be avoided if a non-thermal production mechanism for the right-handed neutrino DM is considered, whereas we need some UV-completion beyond the PeV scale to explain the first one. If the second right handed neutrino decays later than the lightest one (DM) , then we can avoid this fine tuning.
• We also note that the heavy electron partner of the right-handed DM N can also decay via the Higgs coupling with a life time ∼ 10 −8 sec, so that it would not be present after the Universe's temperature of 10 GeV and hence will not have any impact on the evolution of the Universe around BBN.
• In the heavy sector, the heavy quarks will not only form baryonic bound states among themselves but also with the light SM quarks since they share the same QCD. The question then arises if the lightest baryon involving one or more heavy quarks is stable or unstable. In the limit of exact Z 2 symmetry, there will be several stable states i.e. Qqq, QQq, QQQ (where we have not included the proton which in our notation is qqq type). Since the heavy stable baryonic states will have masses in the PeV range, their relic density is likely to far exceed the closure density unless there are resonant effect like in the case of the heavy DM N . One way out of this puzzle would be to introduce soft breaking of the Z 2 symmetry by adding mass terms connecting heavy and light quark states e.g. those in Eq. (2.11) and adjust the parameters (δ U,D ) to make them decay above the QCD phase transition temperature. For this purpose, a value of δ U,D ≥ 10 −4 GeV is sufficient.
To conclude, we have presented a UV-complete model for PeV-scale decaying dark matter in terms of the lightest right-handed neutrinos which are part of an extra SU (2) doublet. The cosmological stability of this DM is guaranteed by a discrete Z 2 symmetry, which is only softly broken by explicit small mass terms. The apparent excess PeV events at IceCube over the standard expectations from a single unbroken astrophysical power-law flux can be understood as due to the 3-body decay of the PeV-scale DM into tau neutrinos ν τ . Our scenario is consistent with other observational constraints, such as the diffuse gamma ray flux and CMB constraints on decaying DM. This model can in principle be tested in future IceCube data with more statistics and more accurate information on the flavor ratio of the observed neutrino events, since it predicts an enhancement of the track events over cascades. The collider and other laboratory tests of this model seem unfeasible, since the non-SM sector lies at or above the PeV scale.