Unifying leptogenesis, dark matter and high-energy neutrinos with right-handed neutrino mixing via Higgs portal

We revisit a model in which neutrino masses and mixing are described by a two right-handed (RH) neutrino seesaw scenario, implying a strictly hierarchical light neutrino spectrum. A third decoupled RH neutrino, $N_{\rm DM}$ with mass $M_{\rm DM}$, plays the role of cold dark matter (DM) and is produced by the mixing with a source RH neutrino, $N_{\rm S}$ with mass $M_{\rm S}$, induced by Higgs portal interactions. The same interactions are also responsible for $N_{\rm DM}$ decays. We discuss in detail the constraints coming from DM abundance and stability conditions, showing that in the hierarchical case ($M_{\rm DM} \gg M_{\rm S}$) there is an allowed window on $M_{\rm DM}$, which necessarily implies a contribution from DM decays to the high-energy neutrino flux recently detected by IceCube. We also show how the model can explain the matter-antimatter asymmetry of the Universe via leptogenesis in the quasi-degenerate limit. In this case, the DM mass should be within the range 300 GeV $\lesssim M_{\rm S}<M_{\rm DM} \lesssim$ 10 PeV. We discuss the specific properties of this high-energy neutrino flux and show the predicted event spectrum for two exemplary cases. Although DM decays, with a relatively hard spectrum, cannot account for all the IceCube high-energy data, we illustrate how this extra source of high-energy neutrinos could reasonably explain some potential features in the observed spectrum. In this way, this represents a unified scenario for leptogenesis and DM that could be tested during the next years with more high-energy neutrino events.


Introduction
The possibility of explaining the dark matter (DM) and baryon asymmetry of the Universe within a unified picture is an attractive idea, intensively explored during recent years [1]. This is also motivated by the simple observation that baryons and DM give a similar contribution to the cosmic energy budget.
Moreover, since both DM and baryon asymmetry of the Universe require new physics, it is conceivable that they should be ultimately explained within a common model. An extension of the standard model (SM) is also required by neutrino masses and mixing established by neutrino oscillation experiments and, therefore, it is quite natural to look at neutrino physics as a possible source for cosmological DM and matter-antimatter asymmetry. 1 The simplest way to describe neutrino masses and mixing is by adding to the SM Lagrangian right-handed (RH) neutrino Yukawa couplings and a Majorana mass term. In the seesaw limit one obtains the seesaw formula for the low-energy neutrino mass, nicely explaining why left-handed (LH) neutrinos are much lighter compared to all other massive fermions [2].
It is then natural to think whether heavy RH neutrinos can play a cosmological role. In the traditional high-energy scale leptogenesis scenario [3], RH neutrino decays are the source of the observed baryon asymmetry. In the νMSM scenario [4], the lightest RH neutrino with O(keV) mass can play the role of DM, with the correct abundance produced by active-sterile neutrino mixing. At the same time, the mixing between the two heavier RH neutrinos can also produce the observed baryon asymmetry [5]. In this model, however, the neutrino Yukawa couplings are many orders of magnitude smaller than those of all other particles and somehow the original motivation of the seesaw is not addressed. 2 As a remedy, it was proposed [8] that if one of the RH neutrinos decouples and it is stable on cosmological time scales, it could play the role of DM and the same mixing among RH neutrinos could also reproduce the correct DM abundance. This is possible if one introduces a non-renormalizable operator λ AB φ † φ N c A N B /Λ, a simple example of Higgs portal models [9], arising from new physics at some scale Λ with additional couplings λ AB . The presence of this new interaction would enhance medium effects opening a new production mechanism from RH-RH neutrino mixing (occurring much above the electroweak scale) rather than from RH-LH neutrino mixing (occurring much below the electroweak scale), as in the νMSM. At the same time it would also be responsible for the DM RH neutrino, N DM with mass M DM , eventually 1 Curiously within errors one has Ω DM,0 /Ω B,0 m atm /m sol 5, where m atm and m sol are the atmospheric and solar neutrino mass scales, respectively. 2 However, in Ref. [6] it was shown how the specific νMSM neutrino Yukawa matrix could arise as a consequence of a lepton number symmetry slightly broken by the Majorana mass terms and Yukawa coupling constants. For more details on the νMSM model we refer the reader to Ref. [7].
to decay, producing high-energy neutrinos 3 and it was pointed out 4 that this flux could be detectable at neutrino telescopes such as IceCube [8]. Intriguingly, the IceCube detector has recently detected the first ever high-energy neutrino events of extraterrestrial origin [14,15], i.e., that cannot be accounted for by the known atmospheric neutrino flux. The energies of these events are as high as O(PeV), i.e., within the natural range of masses needed by the mechanism of cold DM from RH neutrino mixing, as we will discuss in detail. Indeed, very heavy DM decays have been proposed to account for part or all of these events [16,17] and different constraints, within different models, have been presented [18][19][20][21][22][23][24][25][26][27][28][29], showing that this could be a potential explanation of the observed events (or part of them). Encouraged by this phenomenological picture, in this paper we revisit the cold DM RH neutrino mixing scenario discussing a few important aspects and showing how to test it with high-energy neutrino detectors such as IceCube and also showing explicitly how the same set up can accommodate the matter-antimatter asymmetry of the Universe via leptogenesis, as first pointed out 5 in Refs. [8,33,34].
The paper is organised as follows. In Section 2 we review the idea of the seesaw mechanism with one decoupled RH neutrino and the cold DM RH neutrino mixing scenario. In particular we discuss in detail how in the simplest scenario, in order to reproduce the correct DM abundance with a cosmologically stable candidate, M DM has to be in the range 100 GeV-1 EeV in the hierarchical case, for M DM M S . We also discuss a few ingredients, some of which represent quite plausible possibilities, that might relax the constraints especially in the quasi-degenerate limit with M S M DM . In Section 3 we show how the observed matter-antimatter asymmetry can be explained in this model via leptogenesis and in this case the range for M DM restricts to ∼ TeV-10 PeV intriguingly overlapping with the range of neutrino energies detected at IceCube. In Section 4 we discuss the features of the predicted high-energy neutrino flux and its related event spectrum in the IceCube detector, for two exemplary cases for which the DM mass is above a few 100 TeV. We compare these results with the observed event spectrum and show that, although a DM-only signal cannot explain the 4-year IceCube spectrum, it could help to explain some of its features. Finally, in Section 5 we draw our conclusions. 3 For limits on DM decays with neutrino detectors, see Refs. [10][11][12][13]. 4 In Ref. [8] it was noticed that a possible explanation of the PAMELA excess would also imply a potential signal at the IceCube detector. 5 A cosmologically stable RH neutrino playing the role of DM has also been proposed within left-right symmetric models [30], also in combination with leptogenesis [31], and within a hybrid seesaw model [32].

The cold DM RH neutrino mixing scenario
In this section, we revisit the DM RH neutrino scenario [8], combine all constraints using the most recent data and we finally discuss some plausible viable scenarios and the allowed range of values for M DM . In particular we show that relaxing the assumption of ultra-relativistic thermal N S abundance at the resonance, the hierarchical case (M DM M S ) becomes viable extending the range of allowed values for M DM .

From the minimal seesaw Lagrangian to Higgs portal interactions
We assume the usual minimal SM extension with three RH neutrinos N i with Yukawa couplings h and a Majorana mass matrix M . After spontaneous symmetry breaking, the Higgs vacuum expectation value (vev) v generates a neutrino Dirac mass so that the neutrino mass terms can be written as (α = e, µ, τ ; i = 1, 2, 3) where where U is the leptonic mixing matrix and for the light neutrino masses we adopt the convention m 1 ≤ m 2 ≤ m 3 . From neutrino oscillation experiment global analyses we know that [35,36] where c ij = cos θ ij , s ij = sinθ ij , δ is the Dirac phase and ρ and σ are the two Majorana phases. For IO, since we use the convention m 1 < m 2 < m 3 , this has to be cyclically permuted such that  [35,36]). We assume that one of the three RH neutrinos, N DM , has Yukawa couplings small enough to guarantee its stability on cosmological time scales so that it is a potential DM candidate [8]. In this case, necessarily, the neutrino Dirac mass matrix has to be written in one of the three following forms, corresponding effectively to a two-RH neutrino model in which either the lightest RH neutrino N 1 , the next-to-lightest N 2 or the heaviest N 3 is decoupled and has to be identified with N DM . 6 These three forms for m D , with three texture zeros, can be parameterised in terms of nine physical parameters, since three phases can be always reabsorbed in the LH neutrino fields. In the Yukawa basis the Dirac mass matrix is diagonal and we can express it as The transformation from the basis where the charged lepton and Majorana mass matrices are diagonal to the Yukawa basis can be described in terms of two unitary matrices, V L and U R , acting respectively on the LH and on the RH neutrinos, explicitly 6 These special m D forms can be easily justified imposing for example a Z 2 symmetry under which N DM is odd and the other two RH neutrinos are even. The same forms have also been considered in a different context in order to have resonant leptogenesis testable at colliders [37].
Our working assumption, Eq. (7), necessarily implies h A 0 and consequently, from the seesaw formula, Eq. (2), one has m 1 0: in our scenario light neutrinos are strictly hierarchical, either normal hierarchy (NH) or inverted hierarchy (IH). The matrix U R is the RH neutrino mixing matrix and connects the Yukawa eigenstates N J to the mass eigenstates N k : . It can be regarded as the analogue of the leptonic mixing matrix for the light neutrinos and it can be similarly parameterised by three mixing angles θ R ij and three phases. If h A = 0, then N DM is strictly stable and it coincides exactly with N A , implying that the two mixing angles of N DM with the other two RH neutrinos vanish.
We can also conveniently express m D in the orthogonal parameterisation [38], where Ω is the orthogonal matrix encoding information on the RH neutrino total decay widths and total CP asymmetries. The three forms for m D in Eq. (7) necessarily imply, respectively, the three following forms for Ω: where ω is a complex angle and ζ = ±1 is a discrete parameter and the two possible values correspond to two different distinct branches of Ω, with positive and negative determinant respectively [39]. Notice that the nine (real) parameters needed to parameterise the Dirac mass matrix are in this case given by five parameters in the leptonic mixing matrix U (three mixing angles, one Dirac phase, one Majorana phase), two LH neutrino masses, m 2 and m 3 , and finally two real parameters in the complex mixing angle ω.
If Eq. (7) is assumed to hold exactly, then N DM would be strictly stable. However, in this case N DM could not be produced by any interaction, except maybe via gravitational ones, for example at the end of inflation, if it has a mass close to the inflaton mass [40]. In that case, N DM would likely be identified with the heaviest RH neutrino, N 3 . At the same time, it would be questionable whether such a particle would exist at all, not having any interaction. One could think of solving both problems by perturbing the form Eq. (7) for the Dirac mass matrix, introducing some tiny Yukawa coupling h A . In this case N DM would decay with a lifetime (after EW symmetry breaking) From the latest IceCube results, as we will discuss in detail, one has to require τ DM > τ min DM 10 28 s, so that the DM Yukawa coupling would be which is too tiny to think of any DM production mechanism via Yukawa interactions. One possible solution, the so called νMSM model [4], is to have the lightest RH neutrino sufficiently light to dominantly decay into three ordinary neutrinos, in a way that τ 1 ∝ M 5 1 . At the same time, in this way the LH-RH neutrino mixing angle is enhanced and this would induce a sizeable RH neutrino production from mixing. The conditions for the cosmological stability and the correct DM abundance can then be satisfied for a mass of the lightest RH neutrino in the keV range.
However, this solution has the disadvantage of a drastic suppression of all three Yukawa couplings compared to those of all other massive fermions. An alternative possibility is to produce N DM , not necessarily the lightest, through the mixing with the other (two) thermalised RH neutrinos [8]. However, in this case it is easy to see that the mixing cannot be the minimal mixing encoded in the U R matrix, defined by Eq. (8), which describes the mismatch between the Yukawa basis and the basis where the RH neutrino mass matrix is diagonal. This is so since the mixing angles of N DM with the other two RH neutrinos N j would be θ Aj h A /h I (I = B, C), too tiny to produce a sizeable N DM abundance, as Ω DM h 2 ∝ θ 2 Aj . For all practical purposes, we can then neglect such small mixing angles and consider h A = 0.
A way out is to introduce non-standard RH neutrino interactions originating from new physics at an effective scale Λ. At lower energies, this gives rise to a non-renormalizable effective operator, an example of Higgs portal interactions [9], which in the Yukawa basis can be written as 7 (I, J = A, B, C) [8,33,41] Let us show that when introducing this new non-standard interaction, the RH neutrino mixing can lead to the correct N DM abundance.

Estimation of the N DM abundance
In general, the new interaction couplings are non-diagonal in the Yukawa basis and this can provide an efficient source for the RH neutrino mixing. Indeed they would give a contribution 7 Here, we do not refer to any specific model generating this interaction and, therefore, we will treat Λ as a free phenomenological parameter. Below, we will give an example of a simple model able to justify the large values of Λ that we will obtain.
to the effective matter potential of the RH neutrino, which in the Yukawa basis is given by [8] V On the other hand, the Yukawa interactions clearly produce a diagonal contribution to the RH neutrino Hamiltonian, which in the Yukawa basis is given by [42] V Y J = where E J is the N J energy and RH neutrinos are assumed to be in ultra-relativistic thermal equilibrium (of course, this is not true for N DM = N A , since in this case h A = 0). In the mass eigenstates basis one also has the usual kinetic contribution. Let us assume that the mass eigenstate corresponding to N DM mixes with just one of the other two thermalised (mass eigenstate) RH neutrinos. This would play the role of the source RH neutrino, that we refer to as N S , and which is in general a linear combination of N B and N C .
In this way, we have a simple two-neutrino mixing formalism, neglecting for the time being the mixing with the third RH neutrino mass eigenstate N I (the interfering RH neutrino). We will comment on this at the end of this section, showing that indeed the mixing with this RH neutrino can be neglected, since only the mixing with one RH neutrino can satisfy simultaneously all constraints. However, in Section 3 we will see that N I plays a crucial role for leptogenesis (and of course, it is in any case necessary in order to reproduce correctly the neutrino oscillation experimental data).
The source RH neutrino has a Yukawa coupling h S ≡ (h † h) ii , where i is the index corresponding to the mass eigenstate coinciding with N S . If we introduce the effective neutrino mass m S ≡ v 2 h 2 S /M S and parametrise it as m S ≡ α S m sol , then we can write 8 In this way, following the standard procedure, we can write the Hamiltonian for the mixed RH 8 For example, if we again consider the third form for m D in Eq. (7) corresponding to the third form for Ω in Eq. (10) and identify N S with N 2 , then the effective neutrino mass m S can be expressed in terms of ω as Clearly, one has m S ≥ m sol , corresponding to α S ≥ 1. There is no upper bound for α S but values α S 1 correspond to |Ω ij | 2 1, necessarily implying a fine-tuning in the seesaw formula at the level of ∼ 1/|Ω ij | 2 .
neutrinos in the mass eigenstate basis as where we defined Λ ≡ Λ/λ AS and we (reasonably) assumed that non-standard interactions are much smaller than Yukawa interactions, such that h 2 S 2 T / Λ within the relevant temperature range (see below). In the ultra-relativistic limit we can then write E DM p + M 2 DM /(2 p) and E S p + M 2 S /(2 p). As usual, subtracting a contribution to H proportional to the identity that does not contribute to the mixing, we are left with the effective mixing Hamiltonian where we defined ∆M 2 ≡ M 2 S − M 2 DM . If we now describe the neutrino spectrum by its average momentum, p 3 T , and introduce the dimensionless effective potential v Y S ≡ T 2 h 2 S /(4 ∆M 2 ) and the effective mixing angle sin 2θ Λ (T ) ≡ T 3 /( Λ ∆M 2 ), due to the presence of the nonstandard interactions, we can recast ∆H as The energy eigenstates in matter have energies E m DM (T ) and E m S (T ). While the temperature drops down, these tend to get closer to the mass eigenstates N DM and N S and The mixing angle θ m is given by Notice that also the mixing angle in vacuum θ Λ is a function of the temperature. If ∆M 2 < 0, equivalent to having M DM > M S , there is a resonance for v Y = −1. This resonance condition is verified for a specific value of the temperature: It is also useful to introduce the quantity showing that for a fixed M DM and decreasing M S , then z res decreases (i.e. T res increases). This observation will be useful when we will discuss the hierarchical case M DM /M S 2.
At T T res , the coupled RH neutrino N Λ S (T ), the interaction eigenstate, is in thermal ultrarelativistic equilibrium (below we will see that a lower value of the abundance is also possible and even preferred) and this basically coincides with the matter eigenstate with energy E m DM . Since the process is highly non-adiabatic, N Λ S (T ) does not track the matter eigenstate while the temperature drops down. In this way, at T res , just a small fraction of N Λ S (T ) is converted non-adiabatically into the DM RH neutrinos, N Λ DM (T ), the interaction eigenstate that in the absence of mixing would coincide with the mass eigenstate N DM . Let us stress that even after spontaneous symmetry breaking, at zero temperature, there is still a tiny non vanishing mixing angle, θ 0 Λ , such that N Λ,0 DM (the genuine DM state) does not exactly coincide with N DM but also has a tiny N S component, which is one of the reasons why it is not strictly stable as we will see. The fraction of N Λ S converted into N Λ DM can be calculated using the Landau-Zener formula 9 where γ res is the adiabaticity parameter at the resonance, defined as (see, e.g., Ref. [43]) Then, a straightforward calculation gives first and then, using |v Y | res = 2 H res , we finally find 10 where H res 1.66 √ g res T 2 res /M Pl is the expansion rate at the resonance and g res is the number of degrees of freedom at the resonance. This can be assumed to have approximately the SM 9 See Ref. [33] for a derivation within the density matrix formalism. 10 Here we correct the factor in the denominator given in Ref. [8], that was 6 instead of 12.
value plus the contribution from the N S , so that g res = g SM + 7/4 = 108.5. In this way we obtain For the DM abundance one obtains the simple relation where (N DM /N γ ) res is the DM-to-photon number ratio at the end of resonant conversion. If we assume that at the resonance the N S 's are fully thermalised, then (N S /N γ ) res 3/4 and, using Eq. (25), one obtains From Eq. (24), we can express ∆M 2 in terms of z res , Plugging this expression into Eq. (30) and re-expressing h S in terms of α S , we find In the quasi-degenerate limit, the ratio M DM /M S is simply M DM /M S 1.
The assumption that N S is in ultra-relativistic thermal equilibrium at the resonance imposes both a lower bound and an upper bound on z res . The condition that they are ultra-relativistic simply requires M S /T res 3, since for larger values, the N S abundance is Boltzmann suppressed.
There is also a less trivial lower bound. In our setup, the only interactions that can thermalise N S are the Yukawa couplings. Assuming that after inflation the N S abundance is negligible 11 then N S would thermalise at z = z eq , defined as that value of z such that N N S (z eq ) = 1 for initially vanishing N S abundance. Since we are assuming N S ultra-relativistic thermal equilibrium, we then have to impose z res z eq . Inverse decays would thermalise N S at z eq (6/K S ) 1/3 0.4 (100/K S ) 1/3 [44], where K S ≡ m S /m . However, when (2 ↔ 2) scatterings involving top quarks and gauge bosons are also taken into account, the N S thermalisation is more efficient and z eq 8/K S , valid for K S 10.
Further analyses have included more processes and finite temperature effects, such as thermal masses, typically enhancing RH neutrino production and thus, going in the direction of yielding smaller values of z eq [46]. The latest analysis, employing a closed path formalism [47], finds for the total production rate Γ tot , which is equivalent to (D + S)(z 1) 0.2 K S , and implies z eq 5/K S . Notice that, in principle, z eq , and consequently z res , can be made arbitrarily small by making K S arbitrarily large, but a value K S m atm /m 50, corresponding to z eq 0.1, necessarily involves some amount of finetuning in the seesaw formula. In any case, it is convenient to treat, for the time being, z res as a free parameter since it plays a crucial role. Moreover, when we will discuss the case of dynamical N S abundance at the resonance, we will show that it is possible to have arbitrarily small values of z res , although below z eq .
Plugging now the expression for γ res , Eq. (34), into the DM abundance, Eq. (32), one obtains Latest Planck satellite results find for the DM abundance (combining temperature and polarization anisotropies and gravitational lensing) [48], This implies that the correct value of Λ to reproduce the observed DM abundance is given by showing that the mechanism can reproduce the correct DM abundance for reasonable values of Λ DM . 12 In the hierarchical case, z res h S /2 0.85 × 10 −8 α S (M DM /GeV). If z res is set to z res z eq 0.1, then M DM 10 14 GeV/α S . These values imply unacceptably high values of the reheat temperature, since T res is required to be T res 10 15 GeV/α S at the end of inflation. This was the argument used in [8] to rule out the hierarchical case. However, as we said, we will show that values z res 0.1 are actually possible if the assumption of ultra-relativistic thermal N S abundance at the resonance is relaxed. Interestingly, in light of current IceCube high energy 12 One has values Λ DM ≡ Λ/λ AS neutrino data, we will see that this will rescue the possibility to open a window for M DM M S , without any additional non-minimal ingredient or resorting to theoretical uncertainties. On the other hand, in the quasi-degenerate case, one has δ DM ≡ (M DM − M S )/M S 1, so that we can write This shows that, for a fixed value of M DM and a sufficiently small value of δ DM , it is always possible to satisfy the condition for ultra-relativistic thermal equilibrium, z res 0.1.
We have now to worry whether the same new interactions responsible for the production of the abundance can spoil the DM stability on cosmological scales giving unacceptably short lifetimes with high-energy neutrino flux in disagreement with the IceCube data. Of course, at the same time, this instability also represents an opportunity to identify a potential observable contribution to the detected IceCube high-energy neutrinos, as first pointed out in Ref. [8]. This implies that the model has predictive power and can be tested.

DM decays
The N Λ,0 DM decays can proceed through two dominant decay channels [8,33,34]. The first channel is due to the mixing itself that produces the observed DM abundance. Indeed, after electroweak spontaneous symmetry breaking, though the finite temperature-induced mixing is negligible, the operator in Eq. (13) still generates a mixing (in vacuum) between N DM and N S with mixing angle 13 where Γ S ≡ h 2 S M S /(4 π) is the total N S decay width for M DM > M Higgs 125 GeV [50]. In this way, N Λ,0 DM does not exactly coincide with the stable neutrino mass (and Yukawa) eigenstate N DM , but has a tiny (fast-decaying) N S component that would decay quickly into gauge bosons and leptons. As we describe below, the flavour composition of the produced light ν S neutrinos could play an interesting role in the analysis of the high-energy neutrino flux predicted by the mechanism. The decay rate for this process of DM decay via mixing is then simply given by [33,34] 13 It is possible to derive this expression going through the usual lines already reviewed to derive θ Λ (T ).

This can be translated into an expression for the DM lifetime
where Moreover, since the first term can be neglected and we can write, using Eqs. (17) and (33) Finally, imposing the condition in Eq. (37) on Λ in order to obtain the correct DM abundance, the DM lifetime can be written as As we discuss below, IceCube data constrain τ DM τ min DM ∼ 10 28 s and, therefore, a lower bound on M DM is obtained where we defined τ 28 ≡ τ min DM /10 28 s. There is another competing decay channel: the 4 bodydecay process GeV, the decay rate for this process is approximately given by and this implies a DM lifetime 14 In case of W ± emission, ν S would be replaced by ∓ S , i.e., For the value of Λ to reproduce the correct DM abundance, Eq. (37), and for τ DM τ min DM , one finds the upper bound This upper bound is quite stringent but it can be circumvented by requiring M DM − M S 2 M A 200 GeV (which implies the quasi-degenerate limit M DM M S ), since in this way the process is kinematically forbidden. This condition translates into another upper bound Therefore, the upper bound is given by . Along with two-and four-body decays, there would be also three-body decays, N DM → N S + A → 2A + ν S . This decay channel is, however, sub-dominant with respect to four-body decays and we will not consider it. In Fig. 1 we show the diagrams for the different processes in which N DM could decay. 15 The existence of solutions satisfying all discussed physical constraints necessarily requires Here, we describe the constraints on M DM for the quasi-degenerate case under different assumptions, as for instance on the initial N S abundance. As an independent parameter, in addition to M DM , in this case it is more convenient to use z res rather than δ DM or M DM /M S or M S and it is useful to invert Eq. (38) and write Let us now discuss different possibilities highlighting some issues and indicating the remedies. 15 Annihilations N DM +N DM → 2 A should also be considered but these are subdominant, although potentially they might give a signal in dense environments in particular cases. Figure 1: Feynman diagrams of two-, three-and four-body decays of N DM .
(i) Initial N S thermal abundance. If N S has a thermal abundance for arbitrarily small values of z eq ≡ M DM /T eq , then z res can be treated as independent free parameter. 16 In this case, an allowed window of M DM values starts to open up when which is obtained for and Since res , for z res < z max res the range of allowed masses enlarges and the upper bound gets more relaxed, though notice that for z res 10 −6 , corresponding to M max DM 500 TeV, the quasi-degenerate limit does not hold any more. However, we will extend the result to the hierarchical case in Section 2.5. In any case this scenario has the clear drawback that the required small values of z res ∼ 10 −4 imply equally small values of z eq and therefore, some additional interaction able to thermalise N S (but not M DM ).
(ii) Initial N S vanishing abundance. Alternatively, one can wonder whether, starting from an initial vanishing N S abundance, the N S Yukawa interactions could produce an ultrarelativistic thermal N S abundance, without any extra-interaction. Using z eq 5/K S = 5 m / m S = 5 m /(α S m sol ) 0.5 α −1 S , the larger the value of α S , the smaller the value of z eq and hence, the smaller the value of z res could be (recall that z res z eq ). However, the , an allowed range opens up when for α S 10 6 τ 1/2 28 . For larger values of α S , although the allowed window for M DM enlarges, the upper bound becomes even more stringent and our approximations for the calculation of the DM lifetime break down for different reasons. In such a case, N DM becomes lighter than the Higgs boson and also three-body decays become dominant. Moreover, note that for such low values of M DM sphalerons are not effective and leptogenesis is not viable. In addition, such large values of α S imply a huge amount of fine tuning in the seesaw formula. Therefore, this solution is not particularly appealing. 16 As usual, such early thermalisation can be justified in terms of extra gauge interactions associated to annihilations of Z in left-right symmetric models which thermalise N S [51]. Of course, in this case N DM should be a singlet under these new gauge interactions.
(iii) Relaxing the ultra-relativistic N S thermal abundance condition for initial N S vanishing abundance. So far, we have assumed N S to be in ultra-relativistic thermal equilibrium at the resonance, so that the correct abundance of N DM is produced with the highest value of Λ which implies the longest DM lifetime. However, the strong dependence of the lifetime 17 on z res , τ DM→S→ν+φ ∝ z −4 res , actually suggests that if the N S abundance is not too suppressed for z res < z eq , the lower N S abundance can be compensated by a stronger coupling (lower Λ) without spoiling the DM stability. Let us see this quantitatively. The N S abundance can be described in terms of the kinetic equation [44] where D ≡ Γ D /(H z) and S ≡ Γ tot S /(H z), with Γ D and Γ tot S defined as the total decay and scattering rates (2 ↔ 1 and 2 ↔ 2 processes), respectively. For z < z eq 5/K S 0.5 α −1 S and assuming an initial vanishing abundance, the first term describing decays can be neglected and the N S production from the thermal bath is described simply by (using the normalisation N eq Since (D + S)(z 1) K S /5 = 1/z eq , one obtains the simple solution which shows that for z < z eq there is a linear suppression of the N S abundance. Now, using this solution to re-write (N N S /N γ ) res = 3z res /(4z eq ) in Eq. (31), the Eq. (37) for Λ DM has to be replaced (M DM M S ) by Consequently, the lower bound, Eq. (46), valid for 3 z res ≥ z eq , becomes and the upper bound, Eq. (49), 17 This derives from θ 0 Λ ∝ z 4 res , Eq. (39) and (43), so that higher resonant temperatures imply smaller mixing between N DM and N S .
where z res has been replaced by z eq 0.5/α S . The advantage is that now z res does not depend on α S any more and can be taken arbitrarily small, albeit with the condition on the reheat temperature T RH T res = M DM /z res . Imposing again M min DM ≤ M max(A) DM and using z eq 0.5/α S , this time the allowed window opens up at which is realised for and implies δ DM δ min This shows that this solution does not satisfy the quasi-degenerate limit, δ DM 1, and has to be treated more carefully within the hierarchical case, which we discuss in the next subsection, where we see how an allowed window indeed exists and gets enlarged for δ DM 1.
If one compares this solution for M DM with Eq. (52), obtained in the case of free z res and ultra-relativistic thermal N S abundance, clearly the accessible values of M DM are more constrained, but the required small values of z res are now perfectly justified. However, there are still a few options that have to be considered and that can rise the scale for M DM in the quasi-degenerate case.
(iv) Non-thermal N S abundance. We have so far assumed either an initial vanishing N S abundance or a thermal abundance. One could think of a scenario in which, at the end of inflation, part of the inflaton energy density is transferred to N S 's, so they might have an initial abundance effectively much higher than their ultra-relativistic abundance [45].
which is obtained for TeV , for This option has the drawback that if some external mechanism generates an initial N S abundance, the same mechanism might also directly create the final DM abundance. However, there are more (appealing) ways to justify the same results with values of the new parameter ξ 1, which would allow accessing values of M DM above the value ∼ 6 TeV found in Eq. (62) for ξ = 1.
(v) Non-standard expansion rate. Another modification of the minimal scenario that can lead to an increase of the efficiency of the mechanism of DM production is given by the possibility that, at resonance, the cosmological evolution is not in the standard radiationdominated regime with H(T ) ∝ T 2 , but expansion is slower and H res is smaller, which would imply the resonant conversion to be more adiabatic (cf. Eq. (29)). For example, this could happen during a phase transition. In this case, however, entropy production could dilute the N DM abundance. Therefore, the initially produced DM abundance could be larger by a factor ξ , corresponding to an increase of Λ DM by a factor ξ = √ ξ , and thus, a longer DM lifetime (weaker couplings) is possible. Numerically, the conclusions are the same as those in point (iv), simply with a different physical interpretation for ξ. Of course, such a non-standard expansion rate is not easy to motivate and should be understood as a caveat we mentioned for completeness.
(vi) Theoretical uncertainties: improved kinetic description might enhance the efficiency of the mechanism. Our description of non-adiabatic transitions, N S → N DM , is based on the Landau-Zener formula in the monochromatic approximation at zero temperature. These results should be checked within a more rigorous quantum kinetic formalism, which 18 For such a high value of initial N S non-thermal abundance, one would have an initially N S -dominated Universe, unless one has a model with a much higher number of degrees of freedom compared to the SM at very high temperatures. In this case the calculation of the effective potentials and all consequent results, including how the bounds relax with ξ should be revisited. Therefore, the results should be considered more robust for ξ 10, corresponding to M max DM 100 TeV.
accounts for different effects such as finite temperature effects, that might reasonably go into the direction of enhancing the fraction of N S converted into N DM . From this point of view, it should be noticed from Eqs. (25), (34) and (37), that, at the resonance, one has typical values (N DM /N S ) res ∼ 10 −7 (GeV/M DM ), i.e. just a tiny fraction of N S 's is converted into N DM , so one could legitimately wonder whether subtle effects might actually play an important role in the calculation of the correct N DM abundance. In this case, a value of the parameter ξ = 1, introduced in (iv), can also be regarded as a parameterisation of the theoretical uncertainties and values ξ 10 cannot be excluded.
(vii) Running of Λ at low energies. Because of radiative corrections, arising within a specific model from the presence of possible states between the high resonance energy scale and low energies, the value of Λ might vary from very high temperatures at resonance to low energies, and in particular it can increase leading to longer DM lifetimes. This effect would reconcile the DM abundance and DM stability conditions at higher values of M DM and it can again be encoded in terms of the parameter ξ introduced in point (iv) so that the same numerical arguments apply. This possible physical justification of values ξ 1 seems to us quite plausible.
(viii) N DM as a sub-dominant DM component. If N DM constitutes only a fraction ξ < 1 of the observed DM abundance, i.e. Ω N DM = ξ Ω DM , and since the neutrino flux from N DM decays is proportional to Ω N DM /τ DM , the lower limit on τ DM is correspondingly reduced by a factor ξ . Thus, we can identify (ξ ) −2 with the parameter ξ in point (iv). In this case, although the mechanism is not able to reproduce the whole DM relic density, it would still be motivated by the possibility to reproduce the matter-antimatter asymmetry via leptogenesis with testable signatures in neutrino telescopes.
(ix) Initial vanishing N S abundance and ξ = 1. As we have seen in point (iii), for z res < z eq and initial vanishing abundance (lower dynamical N S abundance at resonance), the produced DM abundance, though suppressed, can still reproduce the observed one. In this case, small values of z res are perfectly justified. We now consider the calculation in point (iii), but introducing the parameter ξ, that encodes different possible physical effects. For ξ > 1, one obtains and thus, M min DM (ξ) 2.5 × 10 12 z 1/3 eq z 4/3 res τ 1/3  TeV , which is obtained for and δ min DM gets relaxed by a factor ξ −4/3 with respect to the value found for ξ = 1, Eq. (64). For moderate values ξ 15, one has M DM 40 TeV and δ DM 0.05, so that one is correctly in the assumed quasi-degenerate limit. We will see however that, generalising the discussion to the hierarchical case, for δ DM 1, much higher values of M DM can be easily reached even for ξ = 1.
(x) Mixing with a second thermalised RH neutrino N I ? We have assumed that the new interactions, thanks to the coupling λ AS , mix the mass eigenstate N DM only with N S , but not with N I , which implies a negligible coupling λ AI . One might wonder whether turning on such a coupling might have a beneficial effect, somehow helping to relax the tension. However, at this stage, it should be clear that things can only get worse with a second coupling, since all the same constraints would also apply to this second mixing with an effective scale Λ I ≡ Λ/λ AI . Since all constraints coming from the mixing with N S are weakened by taking the minimal value of m S , it is easy to see that it would actually be impossible to have simultaneously a minimal m I . Therefore, necessarily we have to assume λ AI λ AS (or equivalently Λ I Λ), so that the second mixing is negligible.
We have seen that, assuming the quasi-degenerate limit, M S M DM , the requirement of simultaneously reproducing the correct DM relic abundance and satisfying the stability conditions is possible, either for initial thermal N S abundance (with M DM 500 TeV for ξ = 1) or for initial vanishing N S abundance but only if ξ 1 (in this case M DM 6 ξ 2/3 TeV). On the other hand, as we going to show in Section. 2.5, even for ξ = 1, when relaxing the quasi-degeneracy assumption, much higher values of M DM are possible. TeV , GeV ξ 1/3 . However, as already mentioned for the quasi-degenerate limit, assuming an initial thermal N S abundance is not a particularly attractive case.
Considering the most interesting case of initial vanishing N S abundance, one has to start from Eq. (59) for Λ DM , multiplied by ξ to take into account one (or more) of the possible effects discussed in the previous subsection, explicitly TeV .
If we now consider the upper bound from four-body decays, the upper bound found in point (ix) for the quasi-degenerate case gets relaxed by a factor (M DM /M S ) 2/3 . Explicitly, TeV .
The situation (for initial vanishing abundance) is summarised in Fig. 2 for α S = ξ = 1. In conclusion, we can say that the scenario of DM from RH neutrino mixing, implies a natural window on M DM that is quite an interesting feature of the model, since it naturally predicts high-energy neutrinos from DM decays in the energy range explored by IceCube, as first noticed in Ref. [8]. Intriguingly, as we are going to discuss, this also links neutrinos at the high energies detected by IceCube to ∼ TeV leptogenesis.

Matter-antimatter asymmetry from leptogenesis
So far, for the N DM production, we have considered only two mixed RH neutrinos with a mass splitting δ DM . Now, we also want to take into account the presence of the third interfering RH neutrino, N I , with mass M I (in any case necessary to reproduce correctly the solar and atmospheric neutrino mass scales), in order to have an interference with the source RH neutrino N S giving rise to non-vanishing CP asymmetries for the generation of a matter-antimatter asymmetry via leptogenesis [8]. As we explained in point (x) in the previous section, it is better to have negligible mixing between N DM and N I in order not to increase DM instability.
Since the matter-antimatter asymmetry of the universe is observed today in the form of a baryon asymmetry, this is related to the baryon abundance measured by the Planck satellite given by 19 Ω B,0 h 2 = 0.02226 ± 0.00016 [48]. This can be simply converted into the baryon-tophoton number ratio using where ρ c,0 = (1.05375 ± 0.00013) × 10 −5 h 2 GeV/cm 3 is the critical energy density of the Universe at the present time, h = 0.6751 ± 0.0064 [48] is the Hubble constant, H 0 , in units of 100 km s −1 Mpc −1 and m N is the nucleon mass. Let us see how we can explain this number with leptogenesis at the energy scale enforced by the DM constraints discussed in the previous section.
Since the mixing angle θ Λ (T ) is tiny, now we can completely neglect the mixing due to the non-standard interactions responsible for N DM production and focus just on the interference between N S and N I . Moreover, since N DM has to be heavier than N S , in principle, there are two cases: either M DM = M 2 or M DM = M 3 . In the first case, leptogenesis would occur through the interference between the lightest and the heaviest RH neutrino. In the second case, it would occur via the interference of the two lightest RH neutrinos. The model is effectively a two-RH neutrino model, since N DM is basically decoupled. A detailed analysis of leptogenesis in the hierarchical case of the two-RH neutrino model [39] results in a lower bound, M 1 3 × 10 10 GeV in the NH case and M 1 10 11 GeV in the IH case, and thus, the well known lower bound [52] is in this case even more stringent. Therefore, since from the DM abundance analysis M DM = M 1 and M S needs to be well below 10 10 GeV, at most at the PeV scale, then one is necessarily lead to consider the quasi-degenerate limit for N S and N I , so that the RH neutrino CP asymmetries are enhanced [53,54].
Results on leptogenesis beyond the hierarchical limit, taking into account flavour effects and assuming a two RH neutrino model, were presented in Ref. [55]. It was shown that, in the degenerate limit, when δ lep 10 −2 [56], the lower bound on M 1 is very similar to that one for the unflavoured case, just a factor two weaker, and the asymmetry is ∝ 1/δ lep . Therefore, in our case, since the lower bound has to be relaxed by about five orders of magnitude, we can anticipate that δ lep 10 −5 . As we discussed, in order to satisfy all DM constraints δ DM 10 −2 for initial vanishing N S abundance with ξ 100 and thus, necessarily N DM = N 3 , as δ lep δ DM .
Of course, there are still two possibilities: either N S is the lightest state, i.e., M S = M 1 and M I = M 2 , or the next-to-lightest, i.e., M S = M 2 and M I = M 1 . These two possible cases for the RH neutrino mass spectrum are shown in Fig. 3. Therefore, the interference between the two lightest RH neutrinos generates the matterantimatter asymmetry. The Dirac neutrino mass matrix and, correspondingly, the orthogonal matrix, are then given by the third case in Eq. (7) and Eq. (10), respectively.
Let us verify this estimate by performing a quantitative analysis. Since, as we have seen, the upper bound M max DM requires M S 1 PeV, even allowing for a large ξ, then leptogenesis necessarily occurs in the fully three-flavoured regime [57], so that the asymmetry is the sum of the three contributions from the three charged lepton flavours. At the same time, the asymmetry is the sum of the contribution from the lightest RH neutrino, N 1 , and the contribution from the next-to-lightest RH neutrino, N 2 . Therefore, we can write the final asymmetry as The six individual different contributions can be expressed as where we have introduced the flavoured decay parameters and Γ iα and Γ iα are the flavoured decay rates into leptons and anti-leptons, respectively. The equilibrium neutrino mass is given by The efficiency factors κ(K 1α + K 2α ) can be calculated using where, in our case, x = K 1α + K 2α . This simple expression is strictly valid for initial RH neutrino thermal abundance but since, in any case, the wash-out for the two RH neutrinos adds up and it is necessarily strong in each flavour, there is basically no dependence on the initial RH neutrino abundance. Indeed, notice that, since δ lep 0.01, we are in the degenerate limit, where the wash-out of the two RH neutrinos adds up [56]. The flavoured CP asymmetries are defined as (i, j = 1, 2 and i = j) where Γ i + Γ i = α (Γ iα + Γ iα ) are the total decay rates and the decay parameters are defined They can be calculated using [53] ε iα ε(M i ) where we introduced and Since, in the degenerate limit, the efficiency factor is the same for both RH neutrino contributions (in each flavour), we can rewrite Eq. (79) as Moreover, considering that I α ij = −I α ji and that in the degenerate limit J α ij −J α ji and , the two RH neutrino contributions, for each flavour α, add up (they do not cancel out) and We can now write the different quantities using the orthogonal parameterisation, since this allows us to specify clearly the dependence on the low-energy neutrino parameters. The total and flavoured decay parameters can be written as We can also write Finally, defining and considering that the baryon-to-photon number ratio at recombination is given by η B 0.01 N f B−L , one obtains where the function f (m ν , Ω) has quite a complicated dependence on the different parameters (θ ij , δ, ρ, ω). Similar analytical expressions have been given in Ref. [39] for the hierarchical case.
In Fig. 4, we show the maximal values where we have used the best-fit values θ exp ij of Ref. [35] for the mixing angles. We show density plots for NH (left panel) and IH (right panel), assuming ζ = +1 (the contour plots for ζ = −1 are obtained from those for ζ = +1 with the transformation ω → −ω). The value of f determines the value of δ lep that is needed in order to correctly reproduce the observed asymmetry. Explicitly, for NH (IH), where f max 0.005 (4 × 10 −5 ) is the maximum value of f .

Interplay between leptogenesis and dark matter results
This result clearly confirms what we anticipated, that necessarily M DM = M 3 > M S , M I . It should also be noticed that since the wash-out is necessarily strong, with K 1α + K 2α 5, then z B (K 1α + K 2α ) 4 and this implies the bulk of the asymmetry to be generated at temperature around T lep ∼ M S /z B [44]. From this, imposing conservatively 20 2 T lep T out sph 140 GeV, one 20 The factor 2 provides a conservative lower bound, M S 300 GeV. When the lower bound is saturated, although ∼ 95% of the lepton asymmetry generated within the interval [M S /(z B + 2), M S /(z B − 2)] is not  21 In this way, we have shown that the model can explain both the DM abundance and the matterantimatter asymmetry in a unified scenario. 22 In Table 1 we summarise the results for the allowed window on M DM both for the quasidegenerate and for the hierarchical case (imposing also successful leptogenesis) for the same values of M DM /M S as in Fig. 2. converted into a baryon asymmetry, the residual small fraction of the B − L asymmetry generated outside this interval, at higher temperatures T ∼ [140, 300] GeV, is still converted into a baryon asymmetry and can reproduce the observed asymmetry by further compensating with a value of δ lep lower than the one in Eq. (97). 21 From Eq. (17) for h S , that can be extended to h I by just replacing α S with an analogous quantity α I , one can see that this corresponds to Yukawas h S , h I 10 −6 , 10 −5 (10 −5 , 10 −4 ) for initial (thermal) N S -abundance. This shows that imposing leptogenesis, also in this model one still needs some reduction of the neutrino Yukawa couplings compared to the other massive fermions, although slightly less than in the νMSM. This reduction should be addressed within some full model, able also to specify the origin of the new interactions. 22 In this analysis, we have neglected the contribution to the matter-antimatter asymmetry from N S − N I neutrino mixing itself [5,37,58]. This contribution might relax the constraint on δ lep , Eq. (97), but has no impact on the obtained constraints from DM and on the properties of the high-energy neutrinos testable at IceCube.  Table 2: Allowed windows, imposing successful DM+leptogenesis, for the indicated parameters in the case of initial vanishing abundance for ξ = 1, 10.
In Table 2 we summarise the allowed windows for the different relevant quantities of the model as indicated for the case of initial vanishing abundance and for ξ = 1, 10.
From these two tables one can see how, in the case of a mild hierarchy M DM /M S 5, one has (for ξ = 1) M S 5 TeV and M DM 25 TeV. In this way, one can think of containing all new physics below O(100 TeV), so that electroweak scale stability can be obtained with a reasonable fine-tuning without the necessity of resorting to specific solutions, such as a supersymmetric extension, to the naturalness problem. On the other hand, for M DM /M S 5 one introduces a new very high energy scale, above O(100 TeV), associated to M DM and, if one wants to address naturalness, an extension such as supersymmetry would be needed. In this case the discussed constraints would get modified along similar lines extensively studied already in the case of leptogenesis [59]. The modifications could be again parameterised in terms of a contribution to the parameter ξ, that in the case of supersymmetry would likely be quite mild (ξ SUSY ∼ O(1)).

High-energy neutrinos from DM decays and IceCube data
In this section, we finally discuss the contribution to the very high-energy neutrino flux from DM decays and its properties, comparing the predictions with the most recent IceCube highenergy starting event (HESE) data [49]. As we discussed in Section 2.3, N DM dominantly decays through two-body and four-body processes (see Fig. 1). In this section we consider two-body decays, which occur via the mixing of N S with leptons, Higgs and gauge bosons, with ratio of branching ratios at the source [50], (f e : f µ : f τ ) S , BR(N S → ± W ∓ ) : BR(N S → ν α Z,ν α Z) : BR(N S → Hν α , Hν α ) S = (2 : 1 : 1) S .
We first discuss the flavour composition of the (almost monochromatic) neutrinos directly produced by N S decays, with its distinctive features, and next we derive the event energy spectrum showing two representative choices of (τ DM , M DM ) plus an astrophysical power-law flux, that altogether predict an spectrum in good agreement with the 4-year IceCube data.

Flavour composition of hard neutrinos
Whereas gauge bosons and Higgs decays generate a softer neutrino flux with (1 : 1 : 1) S flavour composition at production, neutrinos produced directly from N S decays retain information of the Yukawa structure of the model (see Eq. (7)) and thus, are of particular interest. These hard neutrinos are also generated from four-body decays, N DM → N S + 2 A → 3 A + ν S , but in this case, their relevance is further diluted compared to the softer neutrinos. They have the highest kinematically allowed energies M DM /2, so they are expected to produce the events from DM decays with the highest energies. Therefore, statistics permitting, by analysing the events close to the high-energy endpoint of the spectrum, information about the flavour composition and thus, about the Yukawa structure, might be inferred. However, let us note that depending on the type of interaction and on the flavour of the incoming neutrino, the deposited energy in the detector (which is the current observable used by the IceCube collaboration) might be quite different from the actual neutrino energy. Indeed, only for electron neutrino and antineutrino charged-current interactions off nucleons, the energy deposited in the detector is close to the neutrino energy. This makes the discrimination of the contribution from the neutrinos at the kinematical threshold a very challenging task, certainly not possible with current data. However, despite these intrinsic experimental difficulties in detecting these (almost monochromatic) neutrinos, we still think it is interesting to discuss their flavour composition, showing that, at production, it can be quite different from the standard mechanisms. Let us then first discuss their flavour composition at production and then at Earth.

Flavour composition of hard neutrinos at production
The flavour composition of monochromatic neutrinos at production is determined by the N Sflavour branching ratios (α = e, µ, τ ) where i = 1 or 2 is the index corresponding to N S . Let us now express the f α 's in terms of the low-energy neutrino parameters using the convenient orthogonal parameterisation. Taking for m D the third form in Eq. (7), recalling that N DM = N 3 , from the orthogonal parameterisation in Eq. (9), one straightforwardly finds Notice that the denominator is m S . This expression clearly implies the upper bound Using the 3σ C.L. lower bounds on the mixing matrix elements of the global fit of [35], one obtains In the left (right) panel of Fig. 5, we show the allowed regions at 1σ C.L., 2σ C.L. and 3σ C.L. (according to the χ 2 -projections provided in Refs. [35,36]) for the flavour fractions at the source, f α,S , for NH (IH). One can see that, as expected, these plots respect the analytical bounds found in Eqs. (102). In Section 4.2, as our benchmark scenario we consider the flavour composition of these hard neutrinos to be (0 : 1 : 1) S .

Flavour composition of hard neutrinos at Earth
In order to compute the final flavour composition at the detector, for all the fluxes in this work, we assume standard neutrino oscillations without exotic interactions. In this case, the neutrino flavour states produced at the source are subject to averaged oscillations in their way to the Earth. Therefore, the probability for a neutrino to arrive at the detector with flavour α, if it was produced with flavour β, is given by and correspondingly the flavour composition at Earth, in terms of the flavour composition at the source, is given by The flavour composition at Earth of the monochromatic neutrinos from DM decays is shown in Fig. 6. Looking at the plot for NH (left panel), we find the main constraint on the flavour composition to be in the electron component, which is restricted to be f e,⊕ 1/3. This feature can be qualitatively understood in the following way. Using the fact that f e,S +f µ,S +f τ,S = 1 and the normalization of P αβ , i.e., β P αβ = α P αβ = 1, one can derive the following expression for the electron flavour content at the detector: f e,⊕ = P eτ + f e,S (P ee − P eτ ) + f µ,S (P eµ − P eτ ) , where P eτ ≈ 1/5, P ee − P eτ ≈ 1/3 and |P eµ − P eτ | ≈ 0.05. Thus, the total electron component at the detector is quite insensitive to the muon component at the source, and we arrive at the Figure 7: Flavour composition at Earth of the hard neutrinos at the high-energy end of the DM decay spectrum in our model (light blue filled area) and four reference cases for normal (left panel) and inverted (right panel) hierarchical spectrum. The large black area is the maximal allowed area for arbitrary flavour composition at the source. In order to produce this plot, the 3σ C.L. ranges of the neutrino oscillation parameters have been used. The red cross is the IceCube best fit point and the areas bounded by the grey lines are the 68% and 95% confidence regions, respectively [60].
approximate relation i.e., the electron component at the detector can be between about 0.2 and 0.6. Our model predicts for NH f e,S 1/3, i.e., also f e,⊕ 1/3. For IH the restrictions are much less pronounced, as can be seen from Fig. 6 (right panel). Finally, let us stress again that the ranges of the flavour composition shown in Figs. 5-7 only correspond to the almost monochromatic flux at the high-energy end of the DM decay spectrum. Given the relatively hard neutrino spectrum produced from gauge boson decays and the IceCube detector capabilities, it will be extremely challenging to single out this contribution.
On the other hand, this could offer important information about the specific Yukawa structure, a definite feature of the model.

Event energy spectrum
Here we consider two representative cases for the DM mass and lifetime and compute the event spectra from two-body DM decays, as described above, expected after 4 years of data taking in IceCube. Since the DM signal alone does not represent a good fit to the entire data sample, we also consider an astrophysical contribution with a power-law flux.
The flux from DM decays has two contributions: galactic and extragalactic, The decays of DM particles at cosmological distances produce a nearly isotropic neutrino and antineutrino flux, which is given by where Ω DM = 0.2618 is the fraction of DM density today, H(z) = H 0 Ω Λ + Ω m (1 + z) 3 is the Hubble expansion rate as a function of redshift, with Ω Λ = 0.6879 and Ω m = 0.3121 [48] (best fit values). The neutrino energy spectrum (of each flavour) from DM decays, dN ν /dE ν , depends on the DM decay channel and on the DM mass, but in our notation we omit to make these dependences explicit. We use the tabulated results in Ref. [61], which include electroweak corrections [62] and were computed using PYTHIA 8.135 [63], and are provided for annihilations of DM particles with mass up to M DM = 100 TeV (or DM decays up to M DM = 200 TeV) and rescale them for higher masses, as done in Ref. [17]. Let us note that this procedure tends to slightly overestimate the final flux, although the precise factor depends on the decay channel and extrapolated value of the DM mass. Thus, our quoted values for the DM lifetime would need to be correspondingly scaled down by the same factor. In addition to the extragalactic signal, DM decays in the Milky Way would also produce a flux of neutrinos and antineutrinos, which would be higher in the Southern hemisphere. Unlike the extragalactic flux, the shape of the galactic flux is not distorted by the redshifting of the signal, which implies differences also in the energy spectrum. The neutrino and antineutrino flux in a direction with galactic coordinates (l, b) produced by DM decays in our own halo is given by where ρ(r) is the DM density profile of the Milky Way as a function of the distance from the galactic centre, r. For a given distance over the line-of-sight, s, the galactocentric distance depends on the galactic longitude, l and the galactic latitude, b, as where R = 8.33 kpc is the distance from the Sun to the galactic centre [64]. For the DM distribution in our galaxy we adopt a generalised Navarro-Frenk-White [65] density profile, with a scale radius r s = 20 kpc, γ = 0.75 and a local DM density ρ(R ) = 0.42 GeV/cm 3 [66]. Let us note, however, that the linear dependence on the density of the flux from DM decays, unlike the quadratic dependence of the flux from DM annihilations, implies smaller uncertainties from the poorly known shape of the DM profile and thus, the precise choice of the DM distribution in the halo is less relevant. In any case, these uncertainties would affect the normalization and the angular dependence of the flux, but not the overall shape of the energy spectrum of the galactic signal. In Fig. 8, assuming NH and the best fit values in Eqs. (5), we show the flavour-averaged neutrino flux, (ν e + ν µ + ν τ )/3, at Earth after propagation for two different DM masses and lifetimes: M DM = 300 TeV and τ DM = 10 28 s (black curves) and M DM = 8 PeV and τ DM = 3×10 28 s (red curves). In the left panel we depict the galactic (dashed curves) and extragalactic (dot-dashed curves) contributions, as well as the total flux (solid curves), whereas in the right panel we show the soft component (dashed curves), i.e., neutrinos from gauge bosons, Higgs and leptons decays and from the related electroweak corrections, and the hard component (dotdashed curves), i.e., neutrinos produced at the decay vertex (including the related electroweak corrections). We note that the galactic flux dominates over the extragalactic contribution and that the soft component of the flux dominates over the hard one, except at the highest energies, close to the kinematical threshold.
Finally, we also compute the event energy spectra for these fluxes and compare them with the 4-year IceCube HESE data. As we note below, the DM decay signal cannot explain all the observed events, so another component of the flux is required. Therefore, in addition to the events from DM decays, we add the contribution from an astrophysical flux described by a power-law spectrum, where φ is the normalization of the flux, in units of 10 −18 GeV −1 cm −2 s −1 sr −1 . For this astrophysical neutrino flux we assume the canonical flavour composition at Earth, (1 : 1 : 1) ⊕ . In this work we do not attempt to perform a fit to this combined model (DM decays plus power-law flux), but simply to show some exemplary cases (see, however, Ref. [68]). Therefore, after choosing some representative values for the DM parameters, we fix the normalization of the astrophysical flux by imposing the total number of events in the electromagnetic(EM)equivalent deposited energy range [60 TeV − 10 PeV] to be equal to the sum of the DM decay and astrophysical signals plus the expected backgrounds. For the atmospheric muon and neutrino backgrounds we scale the 3-year (988 days) IceCube expected numbers [14] to obtain the 4-year (1347 days) expectations, i.e., we consider 3.3 atmospheric neutrino events and 0.6 atmospheric muon events. In order to compute the event spectra of the signal and the background contributions, we closely follow the approach of Ref. [69], which in turn represents an update of the detailed calculations described in Ref. [70] (see also Refs. [71,72]), with some additional improvements. In this work, we use the angular and energy information of the spectra.
In Fig. 9 we see that the DM signal (for M DM = 300 TeV and τ DM = 10 28 s) represents the dominant contribution between 40 TeV and 150 TeV, whereas the hard astrophysical power-law flux (E 2 ν dΦ a /dE ν = 1.6 × 10 −8 GeV cm −2 s −1 sr −1 ) would explain the high-energy part of the observed event spectrum. The small low-energy excess of events with respect to the 3-year results can be nicely explained by neutrinos from DM decays within the scenario described in this paper. This also implies that the astrophysical neutrino flux does not have to be as soft as the result of the fit with only such a flux [49,69]. On the other hand, this hard spectrum is in agreement with the results obtained for the through-going muon sample [15,73], with a per-flavour normalization which is slightly lower, yet compatible within errors. Moreover, let us also note that the through-going muon sample is sensitive mainly to energies from a few 100 TeV to a few PeV, which are precisely the energies in which the astrophysical flux in Fig. 9 is the dominant one. In addition, this hard astrophysical spectrum would not overshoot the γ-ray cascade limit [74][75][76], or the data from air-showers arrays in galactic cases [77,78], if pp sources (where neutrinos are mainly produced from pion decays) are to explain this flux.
Finally, in Fig. 10 we show the event spectra for a heavier DM candidate in combination with a softer power-law flux. In this case, the low-energy events can be explained by the astrophysical flux (E 2 ν dΦ a /dE ν = 6.8 × 10 −8 (E ν /100 TeV) −1 GeV cm −2 s −1 sr −1 ), whereas the prediction of the hard DM decay signal (for M DM = 8 PeV and τ DM = 3×10 28 s) is in agreement with the highest-energy data. However, notice that the (almost) monochromatic flux of hard neutrinos does not translate into a bump in the total event energy spectrum. This is partly due to the particular flavor composition we chose, f e,⊕ 0.2. On the other hand, the natural kinematical cutoff in the event spectrum from DM decays (4 PeV in this case) could also explain the non-observation of events around the Glashow resonance energy (E ν ∼ 6.3 PeV) and, in this case, the through-going muon data could be explained by the hard spectrum from DM decays 23 , instead. Let us finally note that we have not shown the results for lighter M DM O(PeV), usually quoted in the literature [17,21,22,26] when considering earlier IceCube data. However, it does not represent a better agreement with data than the case shown in Fig. 10. Given the soft event spectrum resulting from the 4-year IceCube HESE data, decays from such a DM candidate, cannot explain the entire observed event spectra anymore.

Conclusions and final remarks
As we discussed in detail, in the scenario of cold DM from RH neutrino mixing, the same new interactions are responsible both for N DM production and DM decays, with much stronger predictive power compared to models one can imagine where there is one kind of interaction responsible for production and another responsible for decays (e.g., some tiny small Yukawa coupling) where one has in any case freedom to reproduce both DM abundance and a contribution to IceCube neutrinos. Therefore, finding viable solutions which are able to accommodate leptogenesis, a good DM candidate and are testable signal with neutrino telescopes is highly non-trivial. This is thanks to the possibility to generate the N DM abundance when the N S is still not fully thermalised, an observation that makes viable the hierarchical case with M DM M S . Physically, this relaxes the bounds since N S can be light with a small Yukawa coupling for higher T res , both things helping DM stability (see Eqs. (44) and (47)) and efficiency of production (see Eq. (30)). In this way, even starting from initial vanishing N S abundance, there is an allowed range for the DM mass that, depending on the ratio M DM /M S , extends from ∼ 100 GeV to about ∼ 10 9 GeV (for ξ 5). On the other hand, the higher the value of τ min DM , the narrower the allowed range of masses. For instance from the Fig. 2 one can see that for τ min DM ∼ 10 30 s, the case with M DM /M S = 5 would be ruled out. More generally the upper bound for M DM would become more and more stringent. In addition, the existence of a DM candidate nicely combines with a two-RH neutrino scenario of (resonant) leptogenesis to successfully reproduce the correct baryon asymmetry at a scale below ∼ 10 TeV (1 PeV) for initial vanishing (thermal) N S abundance. In this case, the allowed range of values for M DM narrows to ∼ 1 TeV -10 PeV 23 However, it is not possible to explain in this way the recently announced through-going muon event with deposited energy of 2.6 PeV [79], which is most likely produced by a ∼10 PeV muon neutrino (or a tau neutrino with higher energy) [80,81], unless one considers a much heavier DM candidate. Note, as well, that the kinematical cutoff would not solve the current tension between the lack of Glashow events and the observation of this very energetic through-going muon track. in order to have T lep > T out sph so that a lepton asymmetry can be reprocessed into a baryon asymmetry.
A contribution from N DM decays to the high-energy neutrino flux can help reproducing the IceCube data in addition to an astrophysical component that is, in any case, necessary. Without performing a dedicated fit to the data, we have shown the contribution from DM decays to the energy spectrum for two exemplary masses, M DM = 300 TeV and M DM = 8 PeV, that could help explaining some of the features in the current HESE data. However, we have not investigated which case is preferred (see Ref. [68]). Nevertheless, we do not find that a mass of M DM 4 PeV, discussed in the literature, is particularly favoured with the current data and definitely it cannot explain the entire event spectrum.
In principle, neutrinos produced directly from the decays of N DM via N S retain information on the Yukawa couplings and might be singled out from the rest at energies close to M DM /2. Albeit challenging, this is an interesting possibility, as the flavour composition typically differs from conventional astrophysical flavour ratios. During the next years it will be interesting to see whether more high-energy neutrino events are detected, which could support (depending on the flavor composition) the presence of a component originating from the decays of DM RH neutrinos produced by mixing with source RH neutrinos.