A Model of Metastable EeV Dark Matter

We propose a model where a long-lived pseudoscalar EeV particle can be produced with sufficient abundance so as to account for the cold dark matter density, despite having a Planck mass suppressed coupling to the thermal bath. Connecting this state to a hidden sterile neutrino sector through derivative couplings, induced by higher dimensional operators, allows one to account for light neutrino masses while having a lifetime that can be much larger than the age of the Universe. Moreover, the same derivative coupling accounts for the production of dark matter in the very first instant of the reheating. Given the sensitivity of the IceCube and ANITA collaborations, we study the possible signatures of such a model in the form of Ultra-High-Energy Cosmic Rays in the neutrino sector, and show that such signals could be detected in the near future.

On the other hand, the absolute stability of DM is usually justified by imposing a symmetry. Discrete symmetries are the most popular (R-parity in supersymmetry [40], a Z 2 symmetry in SO (10) unification [27,41], a Z 2 symmetry in Higgs [10,11] or Z-portal [12] models) and can arise from broken gauge symmetries which are exact at the Planck scale. This is not the case for continuous global symmetries which are generically violated at the Planck scale [42][43][44][45]. In this case, the decay of DM is rendered possible through Planck-suppressed operators, as argued in [46]. Due to its very specific signature (monochromatic final states for a 2-body decay), a metastable candidate is regularly evoked when specific detection signals are claimed. For example, a positron excess [47], photon lines [48], high energy neutrinos in IceCube [49,50] or ultra-high energy neutrinos in ANITA [51,52]. However, in each case, the interpretation of a signal as a dark matter detection has to deal with a severe issue: justifying a long lifetime (and thus extremely tiny couplings) while at the same time finding a production mechanism able to produce the dark matter in a sufficiently large amounts to account for the PLANCK determined density of dark matter [1] (implying a coupling which is not so tiny).
At first sight, it would seem that Planck-suppressed couplings of DM particles to Standard Model (SM) states could be sufficient for explaining why dark matter may be long lived on cosmological time scales. For example, one may naively expect the decay width of dark matter to be of order Γ where m denotes the DM mass 1 . However, given the current limit from indi- 1 We will use the reduced Planck mass throughout the paper, arXiv:2003.02846v2 [hep-ph] 9 Mar 2020 rect gamma [53] positron [54] or neutrino [55] detection (τ = Γ −1 10 29 seconds) one would require m 10 keV which reaches the limit from Lyman-α or structure formation constraints [56]. It is, moreover, not an easy task to produce the requisite abundance of DM particles with such feeble couplings. Even the FIMP scenario necessitates couplings of the order of 10 −11 [23], i.e., much larger than m 2 M 2 P .
A potentially more natural way to couple DM to the Standard Model (SM) bath with Planck suppressed couplings is through the neutrino sector, for which there are already strong mass constraints m ν 0.15 eV [57]. Indeed, several constructions invoke a new massive scalar [58] to justify the neutrino mass through a dynamical process similar to the Higgs mechanism applied in the right-handed neutrino sector.
In this work, we show that by combining the violation of global continuous symmetries at the Planck scale, while coupling DM to the neutrino sector, one can generate a large DM lifetime in compliance with the actual experimental constraints for ∼ 1 PeV dark matter masses 2 . The paper is organized as follows. In section II we present our model and we compute the DM lifetime and relic abundance in section III. We propose experimental signatures in section IV. In section V, we propose a top-down model which incorporates all of the needed components for our EeV DM candidate and its coupling to the observational sector. Our conclusions are given in section VI. Appendix A contains additional details on the computation of the decay rates and Appendix B gives more detail on the UV microscopic model containing additional particles and interactions which generate, in the IR, the effective model with appropriate mass scales and couplings that we analyze in the bulk of the paper.

Motivations
We begin with some motivation for the general and more detailed models we present below. The models we 2 Note that the lifetime of the Universe is ∼ 4 × 10 17 seconds, corresponding to ∼ 6.5 × 10 41 GeV −1 whereas limits from indirect detection reach τ 10 29 seconds (∼ 10 53 GeV −1 ).
are proposing rely on a derivative coupling of a DM candidate, a, to matter. Indeed, axionic couplings of the type α M P ∂ µ a appears in several ultraviolet constructions. For instance, in models with string or higher-dimensional inspired moduli fields T = t + ia (see [36] for a more detailed study), they can couple to a sterile sector through the kinetic term as with Z s = 1+ βs M P t+i αs M P γ 5 a and α s , β s real for simplicity. After an integration by parts, the Lagrangian will contain terms which are of the form we consider below.
We can also find such couplings in the Majoron model. Consider a Lagrangian of the type written using a two-component notation and where φ = χe ia M P is the Majoron. After a redefinition of phases, ν s → e − ia 2M P ν s , the kinetic term, −iν sσ µ ∂ µ ν s , produces a coupling of the type given in Eq. (3) [50,59].
Even in string constructions, where we can define the moduli superfield in term of the Grassmannian variables θ andθ by T +T = 2t + 2θσ µθ ∂ µ a, we can show that a term 1 t 2 M P ∂ µ aν s σ µ ν s appears once expanding the Kähler metric as function of matter fields. In this case, α s can be identified as 1 t 2 10 −2 − 10 −3 in KKLT-like models [60]. As one can see, several ultraviolet constructions contains couplings of the type (3) which we use below.

The Lagrangian
Our goal in this section is to build a minimal model of metastable EeV dark matter. By minimal, we mean that we introduce the fewest number of new fields beyond those in the Standard Model with neutrino masses. We assume that DM is a pseudo-scalar field. As alluded to above, the most economical and natural way to proceed is to couple the pseudo-scalar to a sterile neutrino (ν s ) and/or right-handed neutrino (ν R ) sector 3 as is the case for the pseudo-scalar part of the Majoron. We consider the following Lagrangian with In (7), we have included a Yukawa coupling, y s for the sterile neutrino to a SU (2) L Standard Model doublet, L and the Higgs doublet, H, giving rise to a Dirac mass term. We also include a Majorana mass term for ν s . In (8), in addition to coupling the pseudo-scalar to ν R , we add the standard Dirac and Majorana mass terms needed for the see-saw mechanism [62].
As it will be important later when we discuss the production of dark matter during reheating, we also introduce an inflaton, Φ, and couple it to both the sterile sector and the SM, where f corresponds to a SM fermion. Finally, α is a coupling 1 that represents the physics behind the Planck suppressed terms.
The Lagrangian terms in Eqs. (7) and (8) lead to the following neutrino mass matrix: where m s GeV is the SM Higgs vacuum expectation value. We assumed flavor-diagonal couplings in the SM neutrino sector for simplicity and suppressed flavor indices. We also assume the following mass hierarchy (that will be justified in section V) After diagonalization, we can define the 3 mass eigenstates ν 1 , ν 2 , ν 3 by 4 Using the approximation msM R (m R D ) 2 which implies that and where the last inequality in (12) assumes a SM-like neutrino mass of m 2 = 0.05 eV.
In the ν 1 , ν 2 basis, we can rewrite the Lagrangian couplings of a and Φ to the light neutrino sector as and As one can see, our framework is similar to a double seesaw mechanism, and the coupling of the dark matter to the standard model will be highly dependent on the mixing angle θ. Even if couplings of the form in Eq. (7) may seem adhoc, they can in fact be justified by highscale motivated models, an example of which is given in section V.

III. THE CONSTRAINTS
In this section, we consider several necessary constraints on the model. These include constraints on the lifetime from indirect detection searches, constraints on the DM abundance -that is we require a viable production mechanism, and cosmological constraints on the sterile sector from contributions to the effective number of neutrino degrees of freedom, N eff .

Lifetime constraints
The first constraint we apply to the model is on the lifetime of the dark matter candidate a. To be a viable DM candidate, a should at least live longer than the age of the Universe. However, as was shown in [50], when dealing with long-lived decays of particles to the neutrino sector, many body final state decays can dominate over two-body decays when a spin flip makes the amplitude proportional to the neutrino mass in the final state. This is reminiscent of three-body annihilation processes generated by internal brehmshtralung which dominate over two-body annihilation processes suppressed for light fermionic final states due to spin-momentum constraints.
In principle, there are two lifetime limits of importance. First, the DM lifetime (to any final state) must be longer than the age of the Universe (τ a > 4 × 10 17 s). Second, the lifetime must exceed 10 29 s when there is an observable neutrino in the final state [63,64]. In our case (see Appendix A for details), the dominant decay channel is indeed the three-body final state decay Γ a→ν1ν2h/Z and Γ a→ν1eW . All three of these modes have similar amplitudes. Note that we are interested in final states where a SM particle appears, especially an active neutrino, as that gives us the most stringent constraints from experiment 5 . The ν 1 ν 2 h final state is most important and we obtain (see Appendix A) implying that for m 1 m 2 0.05 eV. Note first the rather amazing result that a pseudo-scalar with mass 10 9 GeV, has a lifetime which greatly exceeds the age of the Universe. This is due primarily to the Planck suppressed coupling and the neutrino mass (squared) in the decay rate. Note also that the lifetime of a is determined by the mixing angle θ which we have normalized to 10 −5 requiring a relatively small Yukawa coupling of order 10 −18 from Eq. (12). The smallness of y s will be justified in section V. In this way, we avoid taking α excessively small 6 .
Limits from [64] gives τ a→ν2ν2 5 × 10 28 seconds whereas [63] obtained τ a→bb 10 29 seconds. To be as conservative as possible, we will consider the upper limit m 2 = 0.05 eV for the neutrino mass and τ a 10 29 seconds throughout our work.

Cosmological constraints
Another important constraint comes from the relic abundance of the dark matter. Unless one heavily finetunes the coupling of the inflaton to a, the direct production of a through (two-body) inflaton decay would greatly overproduce the density of a whose annihilation rate would be extremely small. That is, we cannot rely on any kind of thermal freeze-out scenario. It is however possible to produce a in sufficient quantities through 5 Note that the dominant 2-body decay has ν 1 ν 1 in the final state.
But for m 2 θ > 10 −5 m 1 , the three body partial width is always larger (see Eq. (47) in Appendix A). 6 Indeed, the way we wrote the Planck mass coupling α M P imposes α 1 to avoid large transplanckian BSM scales.
In this section, we present a microscopic scenario in which the mass hierarchy that we naturally generated by the spontaneous breaking of a global U(1) symmetry and the presen operators.
First, we introduce a set of heavy Weyl fermion pairs {ψ i , ψ i } i=1,4 and a complex scalar global U (1) symmetry with the charge assignement listed in Tab. I. The most general renorm Tab. I: Charge assignement.
one can write involving these different fields is At tree level, the presence of these heavy fermions leads to interactions between the scalar Φ species, as can be seen from Fig. 1. 4 have masses of the same order of magnitude M mass scale close to the Planck mass, one can integrate those fermions out and obtain the effe operators One thereafter obtains at low energy the Lagrangian the three-body decay of the inflaton coupled only to SM fermions and the sterile sector as in Eq. (6). This allows for a decay channel φ → aν 1 ν 1 as shown in Fig.1.
We assume here some rather generic features of the inflationary sector, and do not need to specify a particular model. We assume a coupling of the inflaton to the SM so that reheating is achieved (we assume instantaneous reheating and thermalization). If dominant, the decay rate is given by where N is an effective number of final state fermionic degrees of freedom and is similar to the total number of degrees of freedom of the Standard Model. If dominant, this decay leads to a reheating temperature given by In general, the relic abundance of dark matter produced in inflaton decay with a branching ratio B R can be expressed as [35] where all masses are expressed in GeV. In our specific case, the partial width for producing the DM candidate a is the three body decay width and the branching ratio for Φ →ν 1 ν 1 a is given by where we have assumed that the total rate is dominated by the two-body decay to SM fermions, or equivalently that N y 2 Implementing the expression for the branching ratio into Eq. (20) we obtain We note that the expression (23) does not depend on the parameter θ, in contrast to the lifetime of a (16). Indeed, the dominant decay channel of the inflaton to neutrinos involves only the lighter state, whereas mixing proportional to θ is compulsory for decays with ν 2 in the final state. We also note that to produce the Planck determined abundance of EeV DM, we need T RH ∼ 10 11 GeV when the Yukawa couplings, y Φ and y f are similar 7 .

Constraints on N eff
It is also important to consider the contribution of neutrino sector present in our model to the overall expansion rate of the universe. In principle, adding a new light degree of freedom, would increase the effective number of light neutrinos which is strongly constrained by the CMB and BBN. The current upper limit is [65] ∆N eff < 0.17 (95% CL), where ∆N eff = N eff − 3. However, a completely sterile neutrino which would never equilibrate with the SM bath would only contribute a small fraction of a neutrino to ∆N eff since its energy density gets greatly diluted compared to the energy density of SM neutrinos [66]. In cases where a light sterile neutrino (or right-handed ν R ) mixes with the active left-handed neutrinos ν L , a non-negligible contribution to N eff may result.
Using Eq. (24), we can derive an upper limit on the mixing angle, θ, by noting that interaction rates for ν 1 are the same as those of active neutrinos, ν 2 , suppressed by θ 2 . Therefore, ν 1 will decouple at a higher temperature, T d1 , than that of ν 2 , T d2 . Indeed we can appoximate T d1 θ 2/3 = T d2 = 2 MeV. As a result, the ratio of the temperatures of ν 1 and ν 2 at T d2 will be given by 7 It is worth mentioning that an alternative possibility would be to produce dark matter only through the graviton-portal but at the price of requiring a very large reheating temperature (T RH 10 14 GeV) as was shown in [29]) where N (T d1 ) is the number of degrees of freedom at T d1 and the number of degrees of freedom at T d2 is 43/4. Furthermore, the contribution to the number of neutrino degrees of freedom will be For example, the upper limit in Eq. (24), yields T 1 /T 2 < 0.64 and N (T d1 ) > 162/4 implying that the decoupling temperature should be greater than Λ QCD . That is decoupling should take place before the QCD transition in the early universe. Thus we arrive at an upper limit, for Λ QCD = 150 MeV.
As one can see, in the cosmologically viable region of interest, the value of y s (and θ) are too weak to be constrained by N eff (at the 2σ level). In contrast, the 1σ upper limit to N eff is 0.05 [65], and in that case, T 1 /T 2 < 0.47 and N (T d1 ) > 407/4 implying that the decoupling temperature should be as large as m t (that is greater than all SM masses). In this case, the limit on θ is significantly stronger, θ < (2MeV/m t ) 3/2 ≈ 4 × 10 −8 . Because the number of degrees of freedom varies slowly with temperature above Λ QCD the limit on θ varies quickly with ∆N eff . At the value of θ ≈ 1.5 × 10 −6 corresponding to α = 0.05, we would predict T d1 ≈ 15 GeV > m b , implying that ∆N eff ≈ 0.062, which may be probed in future CMB missions. In other words, demanding that our model satisfies cosmological constraints therefore favours the region with α ∼ 1 which is more natural from the model-building point of view.

Results
We note at this point that the combination of Eqs. (16) and (23) seem to point toward a natural region of the parameter space with α 10 −2 , θ 10 −6 , and T RH 10 11 GeV which corresponds to y φ y f 10 −5 from Eq. (19). In order to explore this region of the parameter space, we performed a scan on the set of parameters {y s , y R , m s }, while fixing M R = 10 12 GeV and requiring that m 2 = 0.05 eV.
We show in Fig. 2 a scan of the plane (θ, m 1 ≈ m s ) after diagonalization of the mass matrix of Eq. (9). For all points considered, we have fixed the DM mass, m a = 1 EeV, the DM lifetime, τ a = 10 29 s, the active neutrino mass, m 2 = 0.05 eV, and the inflaton mass, m Φ = 3 × 10 13 GeV. We consider three values of α as labelled. In the simple case where the inflaton couples equally to the sterile neutrino ν 1 and SM fermions (y φ = y f ), from Eq. (23) we can satisfy Ω a h 2 0.12 for different values of α by compensating with a different value the reheating temperature. For the values of α chosen, we require T RH = 10 13 , 10 11 and 10 9 GeV as indicated on the figure. The position of the lines of constant α is determined by setting the lifetime to the experimental limit of τ a = 10 29 s which can be read from Eq. 16 (with τ a ≈ Γ −1 a→ν1ν2h ). Note that for m 1 m 2 , the lifetime can be approximated by Eq. (17) which is independent of m 1 which explains why the lines are mostly vertical in the depicted plane. Moreover, since the lifetime is proportional to (αθ) 2 , the choice of α determines θ for constant τ a . For each point of the scan, the value of the corresponding Yukawa coupling y s is indicated by the colored bar. As was anticipated in the previous subsection, in the region of interest the correction to ∆N eff is quite small as compared to the 95% CL upper limit limit ∆N eff < 0.17, which corresponds to the dashed blue contour in Fig. 2.
We also see in the figure, that the coupling y s should be quite small (O (10 −18 ) in the region of interest). We will justify this small coupling in section V. Note that the reason the contribution of ν s to N eff is small, is precisely because the coupling, y s (and mixing angle) is small.

IceCube signals
One clear signature of the model discussed above would be a monochromatic neutrino signal that could be observed by the IceCube, or the ANITA collaborations. In Ref. [51], the case of a scalar DM particle decaying into light right-handed neutrinos was studied and it was shown that the decay of an EeV DM particle followed by the scattering of the RH neutrino within the Earth's crust could lead to visible signals both for ANITA and IceCube for a mixing angle and DM lifetime of order τ a /θ 2 10 27 s. As we have seen, the region of the parameter space that is favored in our model lies towards smaller values of the mixing angle θ 10 −5 and τ a 10 29 s, leading to a ratio τ a /θ 2 10 39 s. This would indicate that our model cannot be detected in searches for anomalous upward-propagating cosmic rays.
In contrast, searches for downward-propagating ultrahigh-energy (UHE) cosmic rays are better suited for signatures of the model discussed here. The IceCube collaboration has reported limits on the decay of darkmatter particles with masses reaching up to a few hundred PeV to active neutrinos. Furthermore, it was shown in Refs. [64,68] that the creation of electroweak showers from the decay of a heavy DM state into neutrinos at very high energy might be constrained at lower energy since the secondary products of such a shower might be visible in the form of a diffuse flux of neutrinos or photons at low energy. These studies led us in the previous sections to impose that the DM lifetime is larger than τ a 10 29 s. In this section we determine the region of parameter space that might be probed experimentally by IceCube, either in the form of direct scattering of UHE neutrinos within the detector, or from secondary electroweak showers which would arrive on Earth at lower energies.

Neutrino Scattering in the IceCube Detector
Let us estimate the number of events which could be detected by IceCube under the form of a monochromatic neutrino signal at ultra-high energies. For that purpose, we suppose that the dark-matter particles follow a Navarro-Frenk-White (NFW) profile [69] : where r s = 24 kpc and the dark-matter density distribution is normalized to equal ρ = 0.3 GeV cm −3 in the vicinity of the solar system [70]. Following Ref. [51], the dark-matter flux, averaged over solid angle, is The number of events predicted for IceCube, assuming a fiducial volume of V IC ≈ (1km) 3 and an exposure time of T exp = 3142.5 days, is given by the relation where the density of the ice is taken to be ρ ice = 0.92 gcm −3 , N A is Avogadro's constant and we estimate the deep-inelastic scattering cross-section σ νN of neutrinos scattering off nuclei using the results of Ref. [71] log 10 (32) Therefore, in the region of the parameter space which we have considered, τ a 10 29 s, the number of events that IceCube might see within the detector is expected to be of order one. Therefore, it is reasonable to suppose that increasing the exposure time by a factor of a few could lead to the detection of such signal in the relatively near future.

Secondary Electroweak Shower Detection
In Ref. [64], limits on the lifetime of a dark-matter particle decaying into active neutrinos have been derived from IceCube data by studying the secondary showers that would be produced by electroweak states at lower energies. We used the limit of Ref. [64] on the lifetime τ a as a function of the dark-matter mass m a in order to translate it into a limit on the mixing angle θ for a fixed set of parameters.

Results
Our results are summarized in Fig. 3 where we have fixed the value of α to one of our previous benchmark points, α = 5 × 10 −2 , and fixed the active neutrino mass to be m 2 = 0.05 eV m 1 . The green-shaded area indicates parameter values in the m a − θ plane for which the lifetime of dark matter would be shorter than the age of the universe. This occurs only at large values of both m a and θ and would be excluded by the lack of events at IceCube. The purple-shaded area excludes the region of the parameter space in which the lifetime of dark matter is shorter than the limit derived in Ref. [64]. Finally the blue lines indicate the value of the mixing angle θ as a function of the DM mass, m a which would lead to 1 event (solid line) or 10 events (dashed line) in the IceCube detector given the exposure time T exp = 3142.5 days. As one can see, our benchmark point (yellow star) corresponding to α = 5 × 10 −2 , m 2 = 0.05 eV m 1 with m a = 1 EeV and τ a = 10 29 s (corresponding to θ 1.5 × 10 −6 ) is flirting with the experimental limits we presented above, suggesting that the prospect for discovery or exclusion of this benchmark is quite high for IceCube, especially for dark-matter masses ranging from the PeV scales to EeV scales. excluded by astrophysical constraints obtained by neutrino detectors (shaded purple). In the green shaded region, the lifetime of the DM candidate is shorter than the age of the Universe. Along the solid (dashed) blue line, the number of events expected by Icecube in its exposure time is 1 (10). The yellow star indicates the benchmark point ma = 1 EeV and θ 1.5 × 10 −6 .
As we have seen, current IceCube data is already on the edge of discovery of EeV dark matter. In its next phase, starting next year, the IceCube collaboration will be able to probe the EeV scale with much better sensitivity for an observable signal.

V. TOWARDS A MICROSCOPIC APPROACH
In this section we develop a toy microscopic model that could justify our assumed hierarchy given in Eq. (10). In fact, such a hierarchy can be generated naturally by the spontaneous breaking of a global U (1) symmetry and the generation of non-renormalizable operators at lowenergy. We introduce a set of heavy Weyl fermion pairs ψ i , ψ i with i = 1, 4 and a complex scalar field S whose charges are given in Table I. One can integrate out these heavy fermions and obtain effective interactions between the scalar S and the different neutrino species, as can be seen from Fig. 4. his section, we present a microscopic scenario in which the mass hierarchy that we introduced in Sec. ?? is lly generated by the spontaneous breaking of a global U(1) symmetry and the presence of non-renormalizable ors. t, we introduce a set of heavy Weyl fermion pairs {ψ i , ψ i } i=1,4 and a complex scalar field Φ, charged under a U (1) symmetry with the charge assignement listed in Tab. I. The most general renormalizable Lagrangian that Tab. I: Charge assignement.
n write involving these different fields is e level, the presence of these heavy fermions leads to interactions between the scalar Φ and the different no species, as can be seen from Fig. 1 Assuming that the fermions {ψ i ,ψ i } i=1,4 have masses of the same order of magnitude M i ∼ M where M is some mass scale close to the Planck scale, after integrating out the heavy fermions, one obtains the effective low energy the Lagrangian in four-component notation. In the above expression, we have introduced the effective couplings where the microscopic couplings λ i are defined in Ap-pendix B.
We assume that the global U (1) symmetry is broken spontaneously at some high energy scale, and the scalar field acquires a vacuum expectation value S = 0. After spontaneous symmetry breaking, one obtains the low energy Lagrangian where we defined If we assume all of the couplings λ j i ∼ 0.1, with the exception of λ 3 L which we take to be ∼ 10 −4 , and a heavy mass scale M M P and a symmetry breaking scale of we naturally get the desired hierarchy of scales which approximates the favored parameter space of our model. Note that in addition to the model we have previously studied, there is an additional mixing term h SR S 3 M 2 pν c s ν R in the seesaw mass matrix, but we checked that for h SR ∼ 10 −3 this term does not perturb the seesaw mechanism or our mass hierarchy.
Due to the charge assignment, a coupling of the inflaton of the type Φν c s ν s in Eq. (6) is not allowed because of the neutrality of Φ under this U (1). It is then impossible to generate a sufficiently large quantity of dark matter through the decay process shown in Fig. 1. However, the production of dark matter via the 3-body decay of the inflaton could be made possible by considering a term like µ Φ Φ|S| 2 , with µ Φ being a dimensionful parameter in addition to the term h R Sν c R ν R with h R ∼ 0.1, as depicted in Fig. (5). In this case, we expect the effective inflaton decay coupling to be as we expect, M S being the mass of the heavy scalar state.

VI. CONCLUSION
The most commonly considered mass ranges for dark matter have been either WIMPs with masses between 100 GeV to 1 TeV, or axions with masses much less than an Tab. I: Charge assignement.
e involving these different fields is the presence of these heavy fermions leads to interactions between the scalar Φ and the different neutrino n be seen from Fig. 1. that the fermions {ψ i ,ψ i } i=1,4 have masses of the same order of magnitude M i ∼ M where M is some ose to the Planck mass, one can integrate those fermions out and obtain the effective, non-renormalizable er obtains at low energy the Lagrangian the four-component notations Here we have explored another regime of dark matter masses of an EeV.
In this paper, we have shown that we can reconcile the dark matter lifetime, which requires extremely reduced couplings, with a natural production mechanism. The long lifetime is possible when Planck-mass suppressed operators (generated, for example, by the breaking of a global symmetry) are combined with tiny neutrino masses. The induced lifetime respects the strongest indirect detection limits once the dark matter is coupled to the neutrino sector. Moreover, despite the feebleness of the coupling, we showed that inflaton decay into dark matter can be sufficient to produce a relic abundance compatible with PLANCK results. Furthermore, we showed that the next generation of neutrino telescopes will be able to probe such heavy dark matter in the near future. In this appendix, we provide some relevant details concerning the computation of the dark matter decay rate. From Eq. (14), we see that up to O(θ 2 ) there are two two-body final state decay channels to consider

Acknowledgments
where we neglected some threshold factors which are negligible in the limit m a m 2 , m 1 . The direct decay of a to two SM-like neutrinos is suppressed by θ 4 .
There are also three-body final state decays which involve a Higgs, W ± , or Z in the final state. These couple to neutrinos in the ν 1 , ν 2 basis through where we used H = (v h + h)/ √ 2 in unitary gauge, where h is the Higgs real scalar field. This leads to the threebody decay width with a Higgs in the final which is given by (44) where we used the relation between the Yukawa coupling and the mixing angle. Similarly, where we assumed m 1 ∼ m 2 m Z m a and e = gg / g 2 + g 2 is the electromagnetic coupling constant. Assuming the hierarchy m 2 ∼ m 1 m e m W m a , gives which is of the same order than the partial width Γ a→ν1ν2Z by using relations between couplings and weak boson masses.
As we have seen, the absence of a helicity flip in the case of the three-body decay compensates the higher power of the Yukawa coupling which arises in the decay width. Thus the ratio of the 3-to 2-body decay widths is and both 2-and 3-body decay modes could be relevant depending on the value of θ and the ratio of light neutrino masses. However, the 3-body decay always dominates over the 2-body decay when SM particles are in the final state Finally, when going to 4-body and higher decay processes, the major change in the decay width, besides complexifying the phase space volume, is an increase in the powers of the Yukawa coupling which naturally leads to smaller widths than the 3−body decay modes.

B. THE MICROSCOPIC MODEL
We provide here the microscopic Lagrangian which allows us to derive the U (1) invariant effective theory in Eq. (33) of Section V. Using a two-component notation, the most general renormalizable and U (1) invariant Lagrangian that one can write involving the fields of the model in Section V is Taking a common mass M i ≡ M for all the supermassive fermions {ψ i ,ψ i } i=1,4 , one can integrate them out to obtain the effective operators M 2 S 2H L L ν s + h.c. (51) One thereafter obtains at low energy the Lagrangian of Eq. (33), written using four-component notation and introducing the effective couplings of Eq. (34).
In full generality, certain interaction or mass terms could be added to the Lagrangian of Eq. (50) while preserving the symmetries of the model. However, we checked that the presence of such terms do not modify the structure of the effective theory we introduce in Sec. V but simply generate additional contributions to the relations of Eq. (34). * emilian.dudas@polytechnique.edu † heurtier@email.arizona.edu ‡ yann.mambrini@th.u-psud.fr § olive@umn.edu ¶ mathias.pierre@uam.es