## Density matrix calculation of the dark matter abundance in the Higgs induced right-handed neutrino mixing model

P. Di Bari , K. Farrag a a,b , R. Samanta a and Y.L. Zhou a a School of Physics and Astronomy , University of Southampton, Southampton, SO17 1BJ, U.K.

- b School of Physics and Astronomy , Queen Mary, University of London London, E1 4NS, U.K.

## August 28, 2020

## Abstract

We present new results on the calculation of the dark matter relic abundance within the Higgs induced right-handed (RH) neutrino mixing model, solving the associated density matrix equation. For a benchmark value of the dark matter mass M DM = 220TeV, we show the evolution of the abundance and how this depends on reheat temperature, dark matter lifetime and source RH neutrino mass M S , with the assumption M &lt; M S DM . We compare the results with those obtained within the Landau-Zener approximation, showing that the latter largely overestimates the final abundance giving some analytical insight. However, we also notice that since in the density matrix formalism the production is non-resonant, this allows source RH neutrino masses below the W boson mass, making dark matter more stable at large mass values. This opens an allowed region for initial vanishing source RH neutrino abundance. For example, for M S glyph[greaterorsimilar] 1 GeV, we find M DM glyph[greaterorsimilar] 20 PeV. Otherwise, for M &gt; M S W ∼ 100 GeV, one has to assume a thermalisation of the source RH neutrinos prior to the freeze-in of the dark matter abundance. This results into a large allowed range for M DM , depending on M S . For example, imposing M S glyph[greaterorsimilar] 300 GeV, allowing also successful leptogenesis, we find 0 5 . glyph[lessorsimilar] M DM / TeV glyph[lessorsimilar] 500. We also discuss in detail leptogenesis with two quasi-degenerate RH neutrinos, showing a case when observed dark matter abundance and matter-antimatter asymmetry are simultaneously reproduced. Finally, we comment on how an initial thermal source RH neutrino abundance can be justified and on how our results suggest that also the interesting case where M DM &lt; M S , embeddable in usual high scale two RH neutrino seesaw models, might be viable.

## 1 Introduction

There are different proposals for extending the Standard Model in a way to explain neutrino masses and mixing and at the same time addressing two of the most compelling cosmological puzzles: dark matter (DM) and matter-antimatter asymmetry of the universe. Such extensions are traditionally based on new physics at energy scales inaccessible with ground laboratories and, therefore, usually untestable. Moreover, one of course would like to consider models that are as minimal as possible. An attractive extension that fulfils both conditions, testability and minimality, and that provides a unified picture of neutrino masses and mixing, dark matter and the matter-antimatter asymmetry of the universe with leptogenesis is the scenario of (cold) dark matter from Higgs induced right-handed neutrino mixing (Higgs induced RHiNo DM) [1, 2]. This is based, in addition to a traditional type-I seesaw Lagrangian extension of the SM with right-handed (RH) neutrinos, on the introduction of a non-renormalizable 5-dim operator,

$$\mathcal { O } _ { A } = \frac { \lambda _ { I J } } { \Lambda } \, \Phi ^ { \dagger } \, \Phi \, \overline { N _ { I } ^ { c } } \, N _ { J } \,,$$

coupling the standard Higgs doublet to RH neutrinos and inducing a RH-RH neutrino mixing, therefore, differently from the usual RH-LH neutrino mixing already rising from the type-I seesaw Lagrangian. We will refer to this operator as the Anisimov operator [3, 4, 1]. It can be regarded as a special case of Higgs portal operator, though not strictly falling within the categories considered in [5].

The interesting feature of the Anisimov operator is that in addition to allow the production of a decoupled RH neutrino playing the role of DM particle, it also predicts a contribution from RH neutrino DM decays to the flux of very high energy neutrinos detectable at neutrino telescopes [1]. Therefore, the recent IceCube neutrino telescope discovery of a very high energy neutrino component in excess of the well known atmospheric contribution [6, 7, 8, 9, 10], prompts the question whether, in addition to an expected, though yet largely undetermined, astrophysical component, the IceCube signal might also receive a contribution of cosmological origin from DM decays. Initial analyses mainly focused on a scenario where heavy DM decays can explain the whole signal and in particular an excess of PeV neutrinos in early data [11, 12, 13]. This possibility seems now disfavoured by the latest data [14], though not completely excluded [15]. However, current IceCube data on the energy spectrum favour the presence of an extra contribution in addition to a traditional astrophysical component described by a power law with spectral index γ glyph[similarequal] 2, as predicted by the Fermi mechanism [16]. Different analyses have shown that in particular the addition of a contribution from DM decays can help in explaining the

IceCube data [17, 18, 19, 20, 15]. In particular, Higgs induced RHiNo DM also provides a good fit to the data for DM masses in the range ∼ 100 TeV-1 PeV [2]. After the IceCube discovery of very high energy neutrinos, various models have been presented that could potentially produce an excess with respect to an astrophysical component [21, 22, 23]. However, the Higgs induced RHiNo DM has the attractive features of minimality and also predictivity, since the same interactions, described by the Anisimov operator (1), can be responsible both for DM production and for its decays. Evidence of this predictive power is that already in the original proposal, prior to IceCube discovery, it was pointed out how the allowed range of DM masses could be probed by neutrino telescopes and in particular by IceCube [1]. At the same time the model explains neutrino masses and mixing within a traditional type-I seesaw mechanism and the matter-antimatter asymmetry of the universe with leptogenesis. In this way a unified picture of neutrino masses, dark matter and leptogenesis, satisfying all experimental constraints, is possible within a certain region in the DM mass-lifetime plane [2].

This intriguing phenomenological picture provides a strong motivation for a more solid calculation of the Higgs induced RHiNo DM relic abundance, a key ingredient for the determination of the allowed mass range. In [1, 2] different approximations and simplified assumptions were adopted. In particular, a simplistic Landau-Zener (LZ) approximation was used to calculate the fraction of source RH neutrino abundance that gets non-adiabatically converted into a DM RH neutrino abundance. In this paper we present results on the calculation of the relic DM abundance using the density matrix formalism.

There are analogies with the calculation of a (light) sterile neutrino abundance from active-sterile neutrino mixing [24] that can lead to a warm DM solution for keV sterile neutrinos [25]. However, the great difference and complication is that in the case of Higgs induced RHiNo DM the vacuum mixing angle vanishes and its role is replaced by a misalignment between the Yukawa interactions and the Higgs induced interactions. This depends on temperature, making the evolution of the system more complicated. In the calculation we still employ a monochromatic approximation and we leave a full momentum dependent calculation for a future investigation. In the final discussion we briefly comment on an extension of the results taking into account momentum dependence.

We also consider the dependence of the relic DM abundance on the initial conditions, in particular the dependence on the initial source RH neutrino abundance and on the reheat temperature.

The paper is structured as follows. In Section 2 we review the model but also improve different results such as the lifetime for two body decays and generalise others. For

example, we notice that the four body decays upper bound on M DM does not apply when the source RH neutrinos are lighter than the W boson. In Section 3 we introduce the density matrix formalism and the equations we solve. In Section 4 we show the evolution of the DM abundance and its dependence on T RH , τ DM and M DM /M S for a benchmark case M DM = 220TeV. We show that the LZ approximation overestimates the relic abundance by many orders of magnitude and provide some analytic insight that explains this result. For this benchmark case we also show an example of how observed DM and matter-antimatter asymmetry of the universe can be simultaneously reproduced. In Section 5 we show the bounds on M DM within different assumptions. In Section 6 we conclude, briefly discussing how the presented results can be extended in different ways. Appendix A contains a derivation of the two body decay rate, Appendix B discusses different equivalent ways to write the density matrix equation, Appendix C contains an original discussion of leptogenesis with two quasi-degenerate RH neutrinos.

## 2 Higgs induced RHiNo dark matter and the LZ approximation

Let us briefly review the Higgs induced RHiNo DM model and how the relic DM abundance is calculated and DM mass constraints derived within the LZ approximation. At the same time we improve and extend some results on the lifetime of the DM RH neutrino. The effective Lagrangian is given by the traditional type-I seesaw Lagrangian [26] with three RHneutrinos with the addition of the Anisimov operator. Before electroweak spontaneous symmetry breaking one has ( α = e, µ, τ and I, J = 1 2 3), , ,

$$- \mathcal { L } _ { M + \Lambda } = \overline { L _ { \alpha } } \, h _ { \alpha J } \, N _ { R J } \, \widetilde { \Phi } + \frac { 1 } { 2 } \, \overline { N _ { R J } ^ { c } } \, D _ { M I J } \, N _ { R J } + \frac { \lambda _ { I J } } { \Lambda } \, \Phi ^ { \dagger } \, \Phi \, \overline { N _ { R J } ^ { c } } \, N _ { R J } + \text{h.c.} \,, \quad ( 2 )$$

where L T α ≡ ( ν Lα , α L ) are the leptonic doublets, Φ is the Higgs doublet, with Φ ˜ ≡ i σ 2 Φ , glyph[star] the h αJ 's are the neutrino Yukawa couplings in the flavour basis where both charged lepton and Majorana mass matrices are diagonal, and we defined D M ≡ diag( M ,M ,M 1 2 3 ), where M 1 ≤ M 2 ≤ M 3 are the three heavy neutrino masses.

After spontaneous symmetry breaking the Higgs vev generates a neutrino Dirac mass matrix m D = v h . One of the three RH neutrinos is assumed to have vanishing Yukawa couplings so that the entries of one of the three columns in h and m D vanish. This assumption can be justified imposing, for example, a Z 2 symmetry. For this reason the seesaw formula,

$$D _ { m } \equiv \text{diag} ( m _ { 1 } = 0, m _ { 2 }, m _ { 3 } ) = U ^ { \dagger } \, m _ { D } \, \frac { 1 } { D _ { M } } \, m _ { D } ^ { T } \, U ^ { * }$$

where U is the leptonic mixing matrix, reduces to the two RH neutrino case with vanishing lightest neutrino mass m 1 = 0, so that the model strictly predicts hierarchical light neutrino masses. In the Yukawa basis the Yukawa matrix is by definition diagonal and given by D h ≡ diag( h A , h B , h C ), with h A = 0 &lt; h B &lt; h C , and the transformation between the two bases is described by a bi-unitary transformation

$$m _ { D } = V _ { L } ^ { \dagger } \, D _ { m _ { D } } \, U _ { R } \,,$$

glyph[negationslash]

where U R acts on RH neutrino fields and can be regarded as the RH neutrino mixing matrix in the absence of Higgs induced interactions described by the Anisimov operator for λ IJ = 0. In this case, the DM RH neutrino Majorana mass eigenstate, that we indicate with N DM , coincides with the Yukawa eigenstate N A and is rigorously stable but also fully decoupled, so that there would be no way to produce it. 1 When Higgs induced interactions are switched on, for λ IJ = 0, they trigger a mixing between N DM and the two coupled RH neutrinos.

At finite temperatures, the Yukawa and the Higgs induced interactions contribute to the RH neutrino self-energies, producing effective potentials that in general are not diagonal in the same basis. The misalignment between the two bases is responsible for RH neutrino mixing. For simplicity, but also because this minimises the constraints from decays, as shown in [2], it is convenient to assume that the RH neutrino mixing is dominantly between the DM RH neutrino N DM and only one of the other two RH neutrinos with non-vanishing Yukawa couplings, that we refer to as the source RH neutrino and we indicate with N S . 2 In this case we can consider a simple two-neutrino mixing formalism and we indicate the coupling between the DM and the source RH neutrino with λ DM S -. We stress again that for λ DM S -= 0 there would be no mixing, because, as we said, N DM would be completely decoupled. Notice also that within this two neutrino mixing formalism, Yukawa basis and flavour basis coincide, i.e., U R = I (see footnote 1).

1 Notice that it does not have to coincide exactly since it could be long-lived with a life-time τ DM glyph[greaterorsimilar] 10 28 s to avoid constraints from IceCube that we discuss later on. In this case, instead of a m D column with three texture zeros one would have some very small entries. However, in this case the mixing angles in U R between the dark and the other two RH neutrinos would have to be so minuscule, and correspondingly the entries in the DM RH neutrino column, that for simplicity they can be set to zero for all purposes. Therefore, the requirement that one mass eigenstate lifetime is sufficiently long-lived to play the role of DM evading experimental constraints, implies, within tiny corrections, both vanishing DM RH neutrino Yukawa couplings and also vanishing mixing with the other two RH neutrinos, encoded in U R (this was discussed in detail in [1]).

glyph[negationslash]

2 N DM and N S have to be regarded as the two RH neutrino Majorana mass eigenstates with M DM = M I and M S = M J for I = J = 1 2 3. They coincide with the energy eigenstates only if , , λ IJ = 0.

The Yukawa interactions clearly produce a diagonal contribution to the RH neutrino Hamiltonian in the Yukawa basis given by [27]

$$V _ { I J } ^ { Y } = \frac { T ^ { 2 } } { 8 \, E _ { J } } \, h _ { J } ^ { 2 } \, \delta _ { I J } \ \ ( I, J = \text{DM}, \text{S} ),$$

where E J is the energy of the Majorana mass eigenstate N J and h S ≡ √ ( h h † ) SS is the Yukawa coupling of N S to the thermal bath (while of course N DM has no Yukawa interactions, since we are assuming h DM = h A = 0). On the other hand, the Higgs induced interactions, described by the Anisimov operator, are in general non-diagonal in the Yukawa basis and they produce an effective potential

$$V _ { I J } ^ { \Lambda } \simeq \frac { T ^ { 2 } } { 1 2 \, \Lambda } \, \lambda _ { I J } \ \ ( I, J = \text{DM}, \text{S} ).$$

In the basis of Majorana mass eigenstates one has also to consider the usual diagonal kinetic contribution so that the Hamiltonian can be written as

$$\mathcal { H } _ { I J } = \left ( \begin{array} { c c } E _ { \text{DM} } & \frac { T ^ { 2 } } { 1 2 \, \Lambda } \\ \frac { T ^ { 2 } } { 1 2 \, \Lambda } & E _ { \text{S} } + \frac { T ^ { 2 } } { 8 \, E _ { \text{S} } } \, h _ { \text{S} } ^ { 2 } \end{array} \right ) \,,$$

where Λ ˜ ≡ Λ /λ DM S -and we assumed that the diagonal terms of the Higgs induced interactions can be neglected. Subtracting a contribution to H proportional to the identity, not affecting the mixing, the effective mixing Hamiltonian is then given by

$$\Delta \mathcal { H } _ { I J } \simeq \left ( \begin{array} { c c } - \frac { \Delta M ^ { 2 } } { 4 p } - \frac { T ^ { 2 } } { 1 6 p } \, h _ { S } ^ { 2 } & \frac { T ^ { 2 } } { 1 2 \widetilde { \Lambda } } \\ \frac { T ^ { 2 } } { 1 2 \widetilde { \Lambda } } & \frac { \Delta M ^ { 2 } } { 4 p } + \frac { T ^ { 2 } } { 1 6 p } \, h _ { S } ^ { 2 } \end{array} \right ) \,,$$

where we defined ∆ M 2 ≡ M 2 S -M 2 DM . If we adopt a monochromatic approximation, so that the momentum is replaced by its average value p glyph[similarequal] 3 T , the effective mixing Hamiltonian in the flavour basis becomes

$$\Delta \mathcal { H } _ { I J } \simeq \frac { \Delta M ^ { 2 } } { 1 2 T } \, \left ( \begin{array} { c c } - 1 - v _ { S } ^ { Y } & \sin 2 \theta _ { \Lambda } \\ \sin 2 \theta _ { \Lambda } & 1 + v _ { S } ^ { Y } \end{array} \right ) \,.$$

Here we have also introduced 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 ), that, as we said, is produced by the misalignment between the Yukawa and the Higgs induced interactions. 3

3 In the standard neutrino mixing case, among left handed neutrinos, this effective mixing angle would correspond to the vacuum mixing angle.

If M DM &gt; M S , implying ∆ M &lt; 2 0, there is a resonance for v Y S = -1, corresponding to a specific value of the temperature, the resonance temperature , given by

$$T _ { \text{res} } \equiv \frac { 2 \, \sqrt { | \Delta M ^ { 2 } | } } { h _ { \text{S} } } = \frac { 2 \, \sqrt { M _ { \text{DM} } ^ { 2 } - M _ { \text{S} } ^ { 2 } } } { h _ { \text{S} } } \,.$$

Since the process is highly non-adiabatic, 4 just a tiny fraction of N S 's, produced by Yukawa interactions, is converted into N DM 's at the resonance. However, since we are considering heavy DM RH neutrinos, even a tiny amount can be sufficient to reproduce the observed DM abundance. Indeed, the relic DM abundance can be expressed in terms of the DM conversion fraction ( N N DM /N N S ) res at the resonance simply as 5

$$\Omega _ { \text{DM} } \, h ^ { 2 } \simeq 1. 4 5 \times 1 0 ^ { 6 } \, \left ( \frac { N _ { N _ { S } } } { N _ { \gamma } } \right ) _ { \text{res} } \, \left ( \frac { N _ { N _ { \text{DM} } } } { N _ { N _ { S } } } \right ) _ { \text{res} } \, \left ( \frac { M _ { \text{DM} } } { \text{GeV} } \right ) \,,$$

where ( N N S /N γ ) res is the source RH neutrino-to-photon ratio at the resonance. This is an important parameter determined by the initial conditions. The latest 2018 Planck satellite results find for the DM abundance at the present time (combining temperature

4 This means that the neutrino states, that are initially produced as source RH neutrinos by Yukawa interactions, do not track the matter eigenstates N m I ( T ) and N m II ( T ) given explicitly by

$$N _ { I } ^ { m } ( T ) \ & = \ N _ { \text{DM} } \cos \theta _ { \text{A} } ^ { m } ( T ) - N _ { \text{S} } \sin \theta _ { \text{A} } ^ { m } ( T ) & ( 1 1 ) \\ N _ { I I } ^ { m } ( T ) \ & = \ N _ { \text{DM} } \sin \theta _ { \text{A} } ^ { m } ( T ) + N _ { \text{S} } \, \cos \theta _ { \text{A} } ^ { m } ( T ) \,,$$

where θ m Λ is the mixing angle in matter describing the transformation from mass-Yukawa eigenstates to matter eigenstates given by

$$\tan 2 \theta _ { \Lambda } ^ { m } = \frac { \sin 2 \theta _ { \Lambda } } { 1 + v _ { S } ^ { Y } } \,,$$

where, since θ Λ ≪ 1, we approximated cos θ Λ glyph[similarequal] 1.

5 For a generic DM scenario one can write [29]

$$\Omega _ { \text{DM} } \, h ^ { 2 } = \frac { n _ { \text{DM} 0 } \, M _ { \text{DM} } } { \varepsilon _ { c 0 } \, h ^ { - 2 } } = \frac { n _ { \gamma 0 } \, M _ { \text{DM} } } { \varepsilon _ { c 0 } \, h ^ { - 2 } \, f ( t _ { f }, t _ { 0 } ) } \, \left ( \frac { N _ { \text{DM} } } { N _ { \gamma } } \right ) _ { f } \simeq 1. 4 5 \times 1 0 ^ { 6 } \, \left ( \frac { N _ { \text{DM} } } { N _ { \gamma } } \right ) _ { f } \, \left ( \frac { M _ { \text{DM} } } { \text{GeV} } \right ), \quad ( 1 3 )$$

where ε c0 glyph[similarequal] 10 54 . h 2 GeVm -3 , n DM0 and n γ 0 glyph[similarequal] 410 7 . × 10 -6 m -3 are respectively the critical energy density, the DM number density and the relic photon number density at the present time. With the subscript 'f' we are indicating the DM abundance freezing time. We are also indicating with f ( t , t f 0 ) = N γ 0 /N f γ = g S f /g S 0 = 106 75 . × 11 43 / glyph[similarequal] 27 3 the dilution factor between the freezing time and the present . time calculated within a standard cosmological model (entropy production is negligible [30]). Possible further dilution due to the same degrees of freedom that play a role in DM genesis, in our case the heavy RH neutrinos, is included in the calculation of ( N N DM /N γ ) . f This expression gets specialised in our case, within a LZ approach, assuming t f = t res and that all the DM abundance is instantaneously produced at t res via non-adiabatic conversions. We will see in the next section how this changes when the relic abundance is calculated within a more realistic density matrix formalism.

and polarization anisotropies and gravitational lensing) [28]

$$\Omega _ { \text{DM} } \, h ^ { 2 } = 0. 1 1 9 3 3 \pm 0. 0 0 0 9 1 \,.$$

Therefore, one can see that one can reproduce the measured value for ( N N DM /N N S ) res ∼ 10 -10 (TeV /M DM ) / N ( N S /N γ ) res , indeed a tiny amount if M DM glyph[greaterorsimilar] 1 TeV, as in our case. This is a basic observation on which the mechanism relies.

A simple way to calculate ( N N DM /N N S ) res , [1, 2], is given by the LZ formula,

$$\frac { N _ { N _ { \text{DM} } } } { N _ { N _ { \text{S} } } } \Big | _ { \text{res} } \simeq \frac { \pi } { 2 } \, \gamma _ { \text{res} } \,,$$

where the adiabaticity parameter at the resonance , γ res , is defined as

$$\gamma _ { \text{res} } & \equiv \frac { | E _ { \text{DM} } ^ { \text{m} } - E _ { \text{S} } ^ { \text{m} } | } { 2 \, | \dot { \theta } _ { m } | } \Big | _ { \text{res} } \,,$$

and in our specific case one finds

$$\gamma _ { r e s } = \sin ^ { 2 } 2 \theta _ { \Lambda } ( T _ { r e s } ) \, \frac { | \Delta M ^ { 2 } | } { 1 2 \, T _ { r e s } \, H _ { r e s } } \simeq 0. 4 \, \frac { M _ { \text{Pl} } \, \sqrt { | \Delta M ^ { 2 } | } } { \widetilde { \Lambda } ^ { 2 } \, \sqrt { g _ { * } ^ { r e s } \, h _ { S } ^ { 3 } } } \,,$$

where g res glyph[star] = g SM glyph[star] + g N S glyph[star] = 106 75+ . g N S glyph[star] glyph[similarequal] 106 75 is the number of ultra-relativistic degrees . of freedom at the resonance given basically by the SM value and in the second numerical equation we used for the expansion rate at the resonance H res glyph[similarequal] 1 66 . √ g res glyph[star] T 2 res /M Pl .

Let us now take into account the constraints on neutrino masses from the seesaw formula and neutrino mixing experimental results. To this extent, it is useful to define the effective neutrino mass associated to the source RH neutrino, ˜ m S ≡ v 2 h /M 2 S S . This provides an easy way to normalise the Yukawa couplings taking automatically into account the seesaw formula and the information on neutrino masses from neutrino mixing experiments. Indeed if we define α S ≡ ˜ m /m S sol , where m sol is the solar neutrino mass scale, then necessarily, from the seesaw formula, one has α S ≥ 1. Notice that h S can be conveniently expressed in terms of α S as

$$h _ { S } ^ { 2 } = \alpha _ { S } \, \frac { m _ { \text{sol} } \, M _ { S } } { v ^ { 2 } } \,.$$

Introducing the variable z ≡ M DM /T and its value at the resonance, z res ≡ M DM /T res , one can express ∆ M 2 in terms of z res finding

$$\sqrt { | \Delta M ^ { 2 } | } = \frac { h _ { S } \, M _ { \text{DM} } } { 2 \, z _ { \text{res} } } \.$$

Using this relation in Eq. (18) and using the definition of α S , one can then conveniently express the adiabaticity parameter at the resonance as [2]

$$\gamma _ { \text{res} } \simeq \frac { 8 } { \alpha _ { \text{S} } \, z _ { \text{res} } } \, \left ( \frac { M _ { \text{DM} } } { M _ { \text{S} } } \right ) \, \left ( \frac { 1 0 ^ { 1 6 } \, \text{GeV} } { \widetilde { \Lambda } } \right ) ^ { 2 } \,.$$

Plugging this expression into the LZ formula Eq. (16), one obtains first ( N N DM /N N S ) res and then, from Eq. (14), for the DM abundance [2]:

$$\Omega _ { \text{DM} } \, h ^ { 2 } \simeq \frac { 0. 1 8 2 2 } { \alpha _ { \text{S} } \, z _ { \text{res} } } \, \left ( \frac { N _ { N _ { \text{S} } } } { N _ { \gamma } } \right ) _ { \text{res} } \left ( \frac { M _ { \text{DM} } } { \text{GeV} } \right ) \, \left ( \frac { M _ { \text{DM} } } { M _ { \text{S} } } \right ) \, \left ( \frac { 1 0 ^ { 2 0 } \, \text{GeV} } { \widetilde { \Lambda } } \right ) ^ { 2 } \,.$$

Imposing that this expression reproduces the measured value Eq. (15), one obtains the value of Λ that reproduces the observed DM abundance ˜

$$\widetilde { \Lambda } _ { \text{DM} } \simeq 1 0 ^ { 2 0 } \, \text{GeV} \, \sqrt { \frac { 1. 5 3 } { \alpha _ { \text{S} } \, z _ { \text{res} } } \, \left ( \frac { N _ { N _ { \text{S} } } } { N _ { \gamma } } \right ) _ { \text{res} } \frac { M _ { \text{DM} } } { M _ { \text{S} } } \, \frac { M _ { \text{DM} } } { \text{GeV} } } \.$$

In this expression one can see that there are five parameters: z res , α S , ( N N S /N γ ) res , M DM and M DM /M S or, alternatively, M S . However, from the relation (20), z res can actually be expressed in terms of α S , M DM and M DM /M S as

$$z _ { \text{res} } = \frac { h _ { \text{S} } \, M _ { \text{DM} } } { 2 \, \sqrt { M _ { \text{DM} } ^ { 2 } - M _ { \text{S} } ^ { 2 } } } \simeq 0. 8 5 \times 1 0 ^ { - 8 } \, \sqrt { \alpha _ { \text{S} } \, \frac { M _ { \text{S} } } { M _ { \text{DM} } } \, \left ( \frac { M _ { \text{DM} } } { \text{GeV} } \right ) } \frac { M _ { \text{DM} } / M _ { \text{S} } } { \sqrt { M _ { \text{DM} } ^ { 2 } / M _ { \text{S} } ^ { 2 } - 1 } } \,,$$

showing that there are actually only four independent parameters. This expression for z res also shows that z res glyph[lessmuch] 1, or equivalently T res glyph[greatermuch] M DM , implying that the reheat temperature T RH &gt; T res cannot be too low within this description. 6 Since there is an upper bound T RH glyph[lessorsimilar] 10 15 GeV, this implies some constraints on the allowed region of parameters.

The dependence on the initial conditions is encoded in the value ( N N S /N γ ) res . If one assumes that some mechanism is able to thermalise the source RH neutrinos prior to the resonant conversion, then ( N N S /N γ ) res = 3 4, and in this case one obtains /

$$\widetilde { \Lambda } _ { \text{DM} } \simeq 1 0 ^ { 2 0 } \, \text{GeV} \, \sqrt { \frac { 1. 1 5 } { \alpha _ { \text{S} } \, z _ { \text{res} } } \, \frac { M _ { \text{DM} } } { M _ { \text{S} } } \, \frac { M _ { \text{DM} } } { \text{GeV} } } \.$$

A more interesting case, since it does not rely on any external mechanism, is to assume that, after inflation, the N S -abundance vanishes and is then produced by the thermal

6 The condition that T RH needs to be higher than T res comes from the fact that at the resonance the Higgs must thermalise in order to produce medium effects via the effective potential (6) generated by the Anisimov operator.

bath through the Yukawa interactions. The production is described by the simple rate equation

$$\frac { d N _ { N _ { S } } } { d z _ { S } } = - ( D _ { S } + S _ { S } ) \left ( N _ { N _ { S } } - N _ { N _ { S } } ^ { \text{eq} } \right ),$$

where we defined z S ≡ M /T S = z M /M S DM , D S ≡ Γ S D / H z ( S ), S S ≡ Γ S S / H z ( S ) and indicated with Γ S D and Γ S S the source RH neutrino total decay and ∆ L = 1 scattering rates respectively and with H the expansion rate. Moreover, we are normalising the abundances in way that the thermal equilibrium N S -abundance is given by

$$N _ { N _ { S } } ^ { \text{eq} } ( z _ { S } ) = \frac { 1 } { 2 } \, \int _ { z _ { S } } ^ { \infty } d x \ x \, \sqrt { x ^ { 2 } - z _ { S } ^ { 2 } } \, e ^ { - x } \,.$$

In particular, in the ultra-relativistic equilibrium, one has N eq N S ( z S glyph[lessmuch] 1) = 1. Since we are now assuming initial vanishing N S -abundance, until decays are negligible compared to inverse decays and N N S glyph[lessmuch] N eq N S , the rate equation is approximated by

$$\frac { d N _ { N _ { S } } } { d z _ { S } } \simeq \left ( D _ { S } + S _ { S } \right ) N _ { N _ { S } } ^ { \text{eq} } \,.$$

Since ( D S + S S )( z S glyph[lessmuch] 1) glyph[similarequal] K / S 5 = 1 /z eq S , where z eq S glyph[similarequal] 0 5 . α -1 S [31, 2] is the value of z S at the time when the N S -abundance thermalises, one obtains the simple solution

$$N _ { N _ { S } } ( z < z _ { S } ^ { \text{eq} } ) \simeq \frac { z _ { S } } { z _ { S } ^ { \text{eq} } }$$

and at the resonance one has [2]

$$\left ( \frac { N _ { N _ { S } } } { N _ { \gamma } } \right ) _ { \text{res} } \simeq \frac { 3 } { 4 } \, \frac { z _ { \text{res} } } { z _ { S } ^ { \text{eq} } } \, \frac { M _ { S } } { M _ { \text{DM} } } \,.$$

Consequently, one obtains for the scale of new physics reproducing the observed DM abundance

$$\widetilde { \Lambda } _ { \text{DM} } \simeq 1 0 ^ { 2 0 } \, \text{GeV} \, \sqrt { \frac { 1. 1 5 } { \alpha _ { \text{S} } \, z _ { \text{S} } ^ { \text{eq} } } \, \frac { M _ { \text{DM} } } { \text{GeV} } }$$

showing that it is indeed convenient having expressed Λ ˜ DM in terms of z res in the general relation (23), since in this way it cancelled out.

These results show that Higgs induced interactions in Eq. (1) are potentially able to reproduce the correct DM abundance for a proper choice of parameters. However, the same Higgs induced interactions are also responsible for the N DM 's to decay at the present time, something that implies both constraints to be imposed but also an opportunity to

test the scenario, in particular by studying the very high energy neutrino flux discovered at IceCube. 7

There are two decay channels to be taken into account. The first one is the two body decay process N DM → A + glyph[lscript] S , where A is a gauge boson and glyph[lscript] S is either a charged lepton or a neutrino with a flavour composition determined by the N S Yukawa couplings [2]. This occurs because even at zero temperature, after electroweak spontaneous symmetry breaking, the Anisimov operator still generates a small (vacuum) mixing angle between N DM and N S given by 8

$$\lambda _ { 0 } = \frac { 2 \, v ^ { 2 } / \widetilde { \Lambda } } { M _ { \text{DM} } \left ( 1 - M _ { \text{S} } / M _ { \text{DM} } \right ) } \,.$$

This mixing results in the two body decay of N DM with rate 9

$$\Gamma _ { \text{DM} \to A + \ell _ { \text{S} } } = \theta _ { \text{A0} } ^ { 2 } \, \frac { h _ { \text{S} } ^ { 2 } } { 4 \, \pi } \ M _ { \text{DM} } \,.$$

Inserting the expression for θ Λ0 , one then obtains for the inverse decay rate

$$\Gamma _ { \text{DM} \rightarrow A + \ell _ { S } } ^ { - 1 } = \frac { \pi } { h _ { S } ^ { 2 } } \, \left ( \frac { \widetilde { \Lambda } } { v ^ { 2 } } \right ) ^ { 2 } \, M _ { \text{DM} } \, \left ( 1 - \frac { M _ { S } } { M _ { \text{DM} } } \right ) ^ { 2 } \,.$$

Using Eq. (19) and imposing Λ = Λ ˜ ˜ DM , with Λ ˜ DM given by Eq. (23), one then finds 10

$$\Gamma _ { \text{DM} \to A + \ell _ { S } } ^ { - 1 } \simeq 1. 2 5 \times \frac { 1 0 ^ { 2 3 } \, \text{s} } { \alpha _ { S } ^ { 2 } \, z _ { \text{res} } } \, \left ( \frac { N _ { N _ { S } } } { N _ { \gamma } } \right ) _ { \text{res} } \frac { M _ { \text{DM} } } { \text{GeV} } \, \left ( \frac { M _ { \text{DM} } } { M _ { S } } \right ) ^ { 2 } \left ( 1 - \frac { M _ { S } } { M _ { \text{DM} } } \right ) ^ { 2 }.$$

IceCube data constrains the lifetime to be longer than τ min DM ∼ 10 28 s, since otherwise an associated high energy neutrino flux would have been observed. Therefore, imposing Γ -1 DM → A + glyph[lscript] S ≥ τ min DM , one obtains a lower bound on M DM .

For initial thermal N S -abundance one has ( N N S /N γ ) res = 3 4 and using the expression / (24) for z res one obtains a lower bound that is much below the Higgs mass and that is, therefore, meaningless since we are assuming N DM to be heavier than the Higgs boson.

7 Plus of course in the next section we need to test how these results change solving the density matrix equation.

8 This expression contains a factor 2 that was missed in [2] and that of course goes in the direction to make bounds more stringent. See Appendix A for details on the derivation.

9 In the absence of the Higgs induced interactions, N DM and N S would coincide with the energy eigenstates and N DM would be stable. However, when Higgs induced interactions are turned on, they generate a small non-diagonal Majorana mass term that breaks the symmetry responsible for the vanishing of the DM RH neutrino Yukawa couplings and its stability (see Appendix A). Indeed, if one considers a Z 2 symmetry, the Anisimov operator is not invariant under this symmetry. Notice that this expression is proportional to M DM , correcting the one given in [2] (proportional to M S ).

10

Notice that GeV

-

1

glyph[similarequal]

6 7

.

×

10

-

25

s.

For initial vanishing N S -abundance one can use Eq. (30) for ( N N S /N γ ) res and in this case one obtains, in the hierarchical case M DM glyph[greatermuch] M S , the lower bound

$$M _ { \text{DM} } \geq M _ { \text{DM} } ^ { \min } \simeq 5 4 \, \text{TeV} \, \alpha _ { \text{S} } \, \tau _ { 2 8 } \, \left ( \frac { M _ { \text{S} } } { M _ { \text{DM} } } \right ) \,,$$

where we defined τ 28 ≡ τ min DM / (10 28 s).

Another important decay channel for N DM at the present time is the four body decay N DM → 3 A + glyph[lscript] S . In the narrow width approximation the decay rate is given by [2]

$$\Gamma _ { \text{DM} \rightarrow 3 A + \ell _ { S } } = \frac { \Gamma _ { S } } { 1 5 \cdot 2 ^ { 1 1 } \cdot \pi ^ { 4 } } \, \frac { M _ { \text{DM} } } { M _ { S } } \, \left ( \frac { M _ { \text{DM} } } { \widetilde { \Lambda } } \right ) ^ { 2 },$$

where Γ S = h 2 S M / S (4 π ). It is important to notice that this expression is valid for M &gt; S M W ∼ 100 GeV. For lower masses the source RH neutrino decays can occur via three body decays, corresponding to five body decays for N DM , and the decay rate is greatly suppressed and does not produce significant constraints. 11

Using again Eq. (19) to express h 2 S in terms of α S and imposing Λ = Λ ˜ ˜ DM (see Eq. (23)) one finds for the inverse decay rate

$$\Gamma _ { \text{DM} \rightarrow 3 A + \ell _ { S } } ^ { - 1 } \simeq \frac { 0. 1 5 3 \times 1 0 ^ { 4 0 } \text{s} } { \alpha _ { S } \, z _ { \text{res} } } \, \left ( \frac { N _ { N _ { S } } } { N _ { \gamma } } \right ) _ { \text{res} } \, \left ( \frac { M _ { \text{DM} } } { M _ { S } } \right ) ^ { 2 } \, \left ( \frac { \text{GeV} } { M _ { \text{DM} } } \right ) ^ { 3 } \,.$$

Imposing again that the lifetime is sufficiently long to escape IceCube constraints implies Γ -1 DM → 3 A + glyph[lscript] S ≥ 10 28 s τ 28 that this time leads to an upper bound on the DM RH neutrino mass given by

$$M _ { \text{DM} } \lesssim 5. 3 \, \text{TeV} \, \alpha _ { \text{S} } ^ { - \frac { 2 } { 3 } } \, z _ { \text{res} } ^ { - \frac { 1 } { 3 } } \, \tau _ { 2 8 } ^ { - \frac { 1 } { 3 } } \, \left ( \frac { N _ { N _ { S } } } { N _ { \gamma } } \right ) _ { \text{res} } ^ { \frac { 1 } { 3 } } \, \left ( \frac { M _ { \text{DM} } } { M _ { \text{S} } } \right ) ^ { \frac { 2 } { 3 } } \,.$$

We can again specialise this upper bound first to the case of initial thermal N S -abundance , for ( N N S /N γ ) res = 3 4, finding /

$$M _ { \text{DM} } \lesssim 4. 8 \, \text{TeV} \, \alpha _ { \text{S} } ^ { - \frac { 2 } { 3 } } \, z _ { \text{res} } ^ { - \frac { 1 } { 3 } } \, \tau _ { 2 8 } ^ { - \frac { 1 } { 3 } } \, \left ( \frac { M _ { \text{DM} } } { M _ { \text{S} } } \right ) ^ { \frac { 2 } { 3 } } \,,$$

and from this, using Eq. (24) for z res , one finds

$$M _ { \text{DM} } \lesssim 0. 3 \, \text{PeV} \, \alpha _ { \text{S} } ^ { - \frac { 5 } { 7 } } \, \tau _ { 2 8 } ^ { - \frac { 2 } { 7 } } \, \left ( \frac { M _ { \text{DM} } } { M _ { \text{S} } } \right ) ^ { \frac { 5 } { 7 } } \, \left [ \frac { M _ { \text{DM} } / M _ { \text{S} } } { \sqrt { M _ { \text{DM} } ^ { 2 } / M _ { \text{S} } ^ { 2 } - 1 } } \right ] ^ { - \frac { 2 } { 7 } } \,.$$

11 Indeed for three body decays, the cross section is phase space suppressed by the fifth power of the mass of the decaying particle. Notice that this case, for M &lt; M S W , has not been considered in [2].

One can notice again that the most conservative bound is obtained for α S = 1, since this minimises the N S Yukawa coupling making N DM more stable. Moreover again higher values of M DM /M S tend to relax also this upper bound.

However, there is an additional constraint coming from the requirement T res ≤ T RH glyph[lessorsimilar] 10 15 GeV translating into an upper bound on M DM /M S that can be derived combining Eq. (24) for z res with Eq. (41) for M max DM , obtaining

$$\frac { M _ { \text{DM} } } { M _ { \text{S} } } \lesssim 8 \times 1 0 ^ { 4 } \, \alpha _ { \text{S} } \, \tau _ { 2 8 } ^ { \frac { 1 } { 6 } } \,,$$

corresponding to an absolute upper bound on the DM mass 12

$$M _ { \text{DM} } \lesssim 1. 0 \times 1 0 ^ { 9 } \, \text{GeV} \, \tau _ { 2 8 } ^ { - \frac { 1 } { 6 } } \,.$$

Considering the case of initial vanishing N S -abundance , substituting Eq. (30) for ( N N S /N γ ) res into Eq. (39), one finds [2]

$$M _ { \text{DM} } \leq M _ { \text{DM} } ^ { \max } \simeq 6 \, \text{TeV} \, \alpha _ { \text{S} } ^ { - \frac { 1 } { 3 } } \, \tau _ { 2 8 } ^ { - \frac { 1 } { 3 } } \, \left ( \frac { M _ { \text{DM} } } { M _ { \text{S} } } \right ) ^ { \frac { 1 } { 3 } } \,.$$

This upper bound combined with the lower bound (36) identifies an allowed window for the value of the DM mass that, however, because of the lower bound on the lifetime τ ≥ τ min DM , opens up only in the hierarchical case, for sufficiently large M DM /M S . Imposing

$$\tau \simeq ( \Gamma _ { D M \to A + \ell _ { S } } + \Gamma _ { D M \to 3 A + \ell _ { S } } ) ^ { - 1 } > \tau _ { D M } ^ { \min } \simeq 1 0 ^ { 2 8 } \, s \,,$$

one finds 13 M DM /M S glyph[greaterorsimilar] 10 α τ S 28 , with the allowed region opening up when the lower bound on M DM /M S saturates at a value M glyph[star] DM glyph[similarequal] 8 TeV.

Notice that since the upper bound (44) applies only for M &gt; M S W ∼ 100 GeV, then one has to impose M DM /M S glyph[lessorsimilar] 10 3 M DM / TeV. This implies that for masses M &gt; M S W the upper bound can be relaxed only up to M DM glyph[lessorsimilar] 150 TeV α -1 2 S τ -1 2 28 .

In Fig. 1 we show in purple, for the most conservative case α S = 1, the allowed range on M DM for M &gt; M S W ∼ 100 GeV, calculated, more accurately, using Eq. (45) that also accounts for the lower bound (36) from two body decays (this, however, holds also for

12 Notice that the bound is saturated for M S = 10 4 GeV α -1 S τ -1 3 28 , so that the assumption M &gt; M S W for the four body decays decay constraints holds.

13 This result is clearly more stringent than the result M DM /M S glyph[greaterorsimilar] 2 3 . α S τ 28 found in [2] because of the more stringent lower bound on M DM as an effect of the two corrections we found to the rate Γ DM → A + glyph[lscript] S and also because we are more accurately taking the inverse of the sum of the rates to calculate the life time in the regime where the two rates are comparable.

M &lt; M S W ). In the case of initial vanishing N S -abundance, the constraint T res &lt; T RH glyph[lessorsimilar] 10 15 GeV is automatically satisfied in the region M &gt; M S W ∼ 100 GeV. For M &lt; M S W , as discussed, the upper bound on M DM from four body decays does not apply and one is left only with the lower bound from two body decays Eq. (36). However, there is still an upper bound on the reheat temperature T res &lt; T RH &lt; 10 15 GeV, from Eq. (24), that implies

$$M _ { \text{DM} } \lesssim 0. 8 5 \times 1 0 ^ { 7 } \, \text{GeV} \, \sqrt { \alpha _ { \text{s} } \, \frac { M _ { \text{s} } } { \text{GeV} } } \,.$$

In Fig. 1 we also show in orange the allowed range on M DM for M &gt; S 1 GeV that is obtained combining the two bounds.

Finally, one can impose constraints from leptogenesis [33]. As we have seen, within the scenario we discussed with M S &lt; M DM , there is quite a stringent upper bound M DM glyph[lessorsimilar] 10 9 GeV (see Eq. (43)). Moreover, the matter-antimatter asymmetry has to be necessarily generated from the decays of the source RH neutrinos, interfering with the third RH neutrino species in order to have non-vanishing CP asymmetries [2]. Since M S glyph[lessmuch] M DM glyph[lessorsimilar] 10 9 GeV, below the lower bound for successful leptogenesis in the two RH neutrino hierarchical case, M lep glyph[lessorsimilar] 10 10 GeV [34], the source and the interfering RH neutrinos have to be necessarily quasi-degenerate in order to have sizeable CP asymmetries resonantly enhanced [35]. Moreover, in order to have successful leptogenesis, the scale of generation of the asymmetry has to be necessarily above the temperature at which sphaleron processes, converting part of the lepton asymmetry into a baryon asymmetry, go out-of-equilibrium, with T off sph glyph[similarequal] 132 GeV [36]. Since in leptogenesis from decays the asymmetry is generated at a temperature that is at most half of the decaying RH neutrino mass, this requirement implies a lower bound [2] M S glyph[greaterorsimilar] 300 GeV, that can be also recast as a lower bound M DM /M S ≤ 3 3 . × 10 -3 M DM / GeV. This upper bound on M DM /M S can be easily combined with the bound on the DM lifetime Eq. (45).

In the case of initial thermal N S -abundance, it is easy to see that this lower bound on M S , combined with the upper bound (41), leads to an upper bound on M DM that is much less stringent than the one (see Eq. (43)) coming from the upper bound on T RH .

On the other hand, in the case of initial vanishing N S -abundance, one finds the approximate allowed region

$$4 \, T e V \, \alpha _ { S } ^ { \frac { 1 } { 2 } } \, \tau _ { 2 8 } ^ { \frac { 1 } { 2 } } \lesssim M _ { \text{DM} } \lesssim 2 7 \, T e V \, \alpha _ { S } ^ { - \frac { 1 } { 2 } } \, \tau _ { 2 8 } ^ { - \frac { 1 } { 2 } } \,.$$

This is clearly more stringent both than the upper bound we derived for M &gt; M S W and for M &lt; M S W , in this second case from the upper bound on the reheat temperature. This allowed region for successful leptogenesis, calculated more precisely from Eq. (45),

is shown in green in Fig. 1 for the most conservative case α S = 1 and one can notice that it is quite restricted. 14 In particular, one can notice that in this case there is quite a stringent upper bound on the DM lifetime τ DM glyph[lessorsimilar] 4 × 10 28 s. The existence of this upper bound shows that the possibility to combine DM with leptogenesis within this model will be certainly tested in the next future at neutrino telescopes. However, such a marginal allowed region legitimately questions whether a calculation of the DM abundance within the simple LZ approximation gives the correct results, thus motivating a calculation within a density matrix formalism.

## 3 Density matrix formalism

In this section we go beyond the LZ approximation and calculate the DM relic abundance within the density matrix formalism [38]. The use of density matrix in neutrino physics in the early universe has a long history. The most traditional application is the study of active-sterile neutrino mixing in the early universe [39]. In that case a comparison between the LZ approximation and the density matrix formalism was made in [40] finding quite a good agreement. The use of a density matrix formalism plays also a crucial role in the study of RH neutrino mixing in leptogenesis from neutrino oscillations [37, 41]. The density matrix formalism also proves to be very important in the description of flavour effects in leptogenesis [42].

In the absence of Higgs induced interactions, the only interactions able to produce the source RH neutrinos would be the Yukawa interactions so that the N DM 's would be completely decoupled. Therefore, Yukawa interactions would produce only source RH neutrinos (barring the third RH neutrino species for the time being). This production can be described by a density matrix normalised in terms of the source RH neutrino abundance that in the Yukawa basis would be diagonal and simply given by ( I, J = DM S) ,

$$N _ { I J } ( z ) = N _ { N _ { S } } ( z ) \left ( \begin{array} { c c } 0 & 0 \\ 0 & 1 \end{array} \right ).$$

Here again we notice that we describe the system within a monochromatic approximation where momentum dependence is integrated away. As we have seen, the abundance of source RH neutrinos, N N S , is described by the simple rate equation (26). However, when

14 It is interesting to notice that if one would consider leptogenesis from RH neutrino oscillations, the so-called ARS scenario [37], then since the asymmetry is produced when the source RH neutrino is ultra-relativistic, the source RH neutrino mass can be much lighter and this would certainly highly relax the constraint. A dedicated analysis would be certainly interesting in this respect.

Figure 1: Summary of the allowed regions in the M DM -τ DM plane obtained within the LZ approximation for α S = 1 imposing different requirements. The region below the dotted black horizontal line, for τ DM &lt; 10 28 s, is currently excluded by IceCube. The light blue region delimited by the long dashed line is the region satisfying T RH &lt; 10 15 GeV for initial thermal N S -abundance. The orange region is for M &gt; S 1 GeV and vanishing N in N S . The vertical line corresponds to the upper bound on M DM Eq. (46) from T RH &lt; 10 15 GeV. The green region satisfies the lower bound M &gt; S 300 GeV allowing also for successful leptogenesis. The red star, the black dotted line and the gray dotted line are respectively the best fit, 68% and 95% contour lines recently found in [19] analysing latest IceCube data including a contribution from DM (neutrinophilic) decays in addition to an astrophysical component with Fermi spectrum.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00521v2-figure-96d3aaff-9e40-4621-816b-98054232ec27.jpg)

the Higgs induced interactions are turned on, they develop off-diagonal terms that have to be taken into account together in principle with decoherence effects. This evolution is then described by a density matrix equation of the form [39]

$$\frac { d N _ { I J } } { d t } = - i \left [ \mathcal { H }, N \right ] _ { I J } - \begin{pmatrix} 0 & \frac { 1 } { 2 } ( \Gamma _ { D } + \Gamma _ { S } ) \, N _ { \text{DM-S} } \\ \frac { 1 } { 2 } ( \Gamma _ { D } + \Gamma _ { S } ) \, N _ { \text{S-DM} } & ( \Gamma _ { D } + \Gamma _ { S } ) \, ( N _ { N _ { S } } - N _ { N _ { S } } ^ { \text{eq} } ) \end{pmatrix} \,, \quad ( 4 9 )$$

where the first term is the Liouville-von Neumann term and the second term is the combination of the decoherence term, damping off-diagonal terms, and the repopulation (diagonal) term, describing the production of source RH neutrinos.

Clearly, without off-diagonal terms in the Hamiltonian, the density matrix equation would simply reduce to Eq. (26). Moreover again a diagonal term in H cancels out and we can replace H → ∆ H , with ∆ H given by Eq. (8).

As often done, we can express the matrices in the Pauli matrix basis using a vectorial notation. The effective Hamiltonian in Eq. (8) can then be recast as

$$\Delta \mathcal { H } = \frac { 1 } { 2 } \, \vec { V } \cdot \vec { \sigma } \,,$$

where the effective potential vector glyph[vector] V is defined as

$$\vec { V } \equiv \frac { \Delta M ^ { 2 } } { 6 \, T } \, \left ( \sin 2 \theta _ { \Lambda }, 0, - 1 - v _ { S } ^ { Y } \right ) \,.$$

The abundance density matrix is analogously recast, introducing the quantity P 0 and the polarisation vector glyph[vector] P , as [43]

$$N = \frac { 1 } { 2 } \, P _ { 0 } \, \left ( 1 + \vec { P } \cdot \vec { \sigma } \right ) \,,$$

in a way that

$$N _ { N _ { \text{DM} } } \ & = \ \frac { 1 } { 2 } \, P _ { 0 } \, \left ( 1 + P _ { z } \right ) \,, \\ N _ { N _ { \text{S} } } \ & = \ \frac { 1 } { 2 } \, P _ { 0 } \, \left ( 1 - P _ { z } \right ) \,, \\ N _ { N _ { \text{DM} } } + N _ { N _ { \text{S} } } \ & = \ \text{$P_{0}$} \,. \\ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot$$

Inserting Eqs. (50) and (52) into the density matrix equation (49), one obtains a set of equations for P 0 and glyph[vector] P 15

$$\frac { d \vec { P } } { d t } \ \ = \ \vec { V } \times \vec { P } - \left [ \frac { 1 } { 2 } ( \Gamma _ { D } + \Gamma _ { S } ) + \frac { d \ln P _ { 0 } } { d t } \right ] \, \vec { P } _ { T } - ( 1 + P _ { z } ) \, \frac { d \ln P _ { 0 } } { d t } \, \hat { z } \,,$$

$$\frac { d P _ { 0 } } { d t } \ = \ - ( \Gamma _ { D } + \Gamma _ { S } ) \, ( N _ { N _ { S } } - N _ { N _ { S } } ^ { \text{eq} } ) \,,$$

15 In Appendix B we give some details on the derivation and we also show a third alternative equivalent way to write the density matrix equation often used in the literature.

where we defined glyph[vector] P T ≡ P x x ˆ + P y y ˆ. If we explicitly unpack the first vectorial equation in terms of its components, we obtain the following set of four differential equations

$$\text{of its components, we obtain the following set of four differential equations} \\ \frac { d P _ { x } } { d t } \ & = \ - V _ { z } P _ { y } - \frac { 1 } { 2 } ( \Gamma _ { D } + \Gamma _ { S } ) \, P _ { x } - \frac { P _ { x } } { P _ { 0 } } \, \frac { d P _ { 0 } } { d t } \,, \\ \frac { d P _ { y } } { d t } \ & = \ V _ { z } \, P _ { x } - V _ { x } \, P _ { z } - \frac { 1 } { 2 } \, ( \Gamma _ { D } + \Gamma _ { S } ) \, P _ { y } - \frac { P _ { y } } { P _ { 0 } } \, \frac { d P _ { 0 } } { d t } \,, \\ \frac { d P _ { z } } { d t } \ & = \ V _ { x } \, P _ { y } - \frac { 1 + P _ { z } } { P _ { 0 } } \, \frac { d P _ { 0 } } { d t } \,, \\ \frac { d P _ { 0 } } { d t } \ & = \ - ( \Gamma _ { D } + \Gamma _ { S } ) \, ( N _ { N _ { S } } - N _ { N _ { S } } ^ { \text{eq} } ) \,. \\ \text{ing the independent variable, from } t \, \text{ to } z, \, \text{one obtains}$$

$$d t$$

Changing the independent variable, from t to z , one obtains

$$\text{anging the independent variable, from } t \, \text{to } z, \, \text{one obtains} \\ \frac { d P _ { x } } { d z } \ & = \ - \overline { V } _ { z } \, P _ { y } - \frac { 1 } { 2 } \, \frac { M _ { S } } { M _ { \text{DM} } } \, ( D + S ) \, P _ { x } - \frac { P _ { x } } { P _ { 0 } } \, \frac { d P _ { 0 } } { d z } \,, \\ \frac { d P _ { y } } { d z } \ & = \ \overline { V } _ { z } \, P _ { x } - \overline { V } _ { x } \, P _ { z } - \frac { 1 } { 2 } \, \frac { M _ { S } } { M _ { \text{DM} } } \, ( D + S ) \, P _ { y } - \frac { P _ { y } } { P _ { 0 } } \, \frac { d P _ { 0 } } { d z } \,, \\ \frac { d P _ { z } } { d z } \ & = \ \overline { V } _ { x } \, P _ { y } - \frac { 1 + P _ { z } } { P _ { 0 } } \, \frac { d P _ { 0 } } { d z } \,, \\ \frac { d P _ { 0 } } { d z } \ & = \ - \frac { M _ { S } } { M _ { \text{DM} } } \, ( D + S ) \, ( N _ { N _ { S } } - N _ { N _ { S } } ^ { \text{eq} } ) \,, \\ \text{ere we have already defined } D \text{ and } S \, \text{after} \, \text{Eq.} \, (26) \text{ and we have now also introduced}$$

where we have already defined D and S after Eq. (26) and we have now also introduced glyph[vector] V ≡ glyph[vector] V / H z ( ).

In the next section we show the evolution of the DM abundance obtained solving numerically this set of density matrix equations for a benchmark value M DM = 220TeV and for different values of T RH , τ DM and M DM /M S .

## 4 Evolution of the DM abundance from the density matrix equation

In this section we fix the DM mass to a benchmark value M DM = 220TeV and we show the evolution of the DM abundance, N N DM , solving the density matrix equation presented in the previous section. We choose this particular benchmark value for M DM since it is the best fit value of DM mass found in [19], where the authors analysed IceCube data on the high energy neutrino flux energy spectrum within a model where, in addition to an astrophysical component with a power-law spectrum with spectral index γ = 2 2, there is . an additional contribution from neutrinophilic DM decays. 16

16 The analysis does not straightforwardly translates to our model but it provides a good indication and motivation for the use of such value of the DM mass as benchmark value.

Though we fix M DM , we show how the evolution of the DM abundance depends on the other three parameters of the model: the reheat temperature T RH , the lifetime τ DM and finally the ratio M DM /M S (or equivalently M S considering that M DM is fixed).

In all plots we also show the relic value of the DM abundance, indicated with N f , obs N DM , that corresponds to the observed value for Ω DM h 2 given in Eq. (15). This can be easily derived from Eq. (13), finding 17

$$N _ { N _ { D M } } ^ { f, o b s } = ( 0. 1 0 9 7 \pm 0. 0 0 0 8 ) \times 1 0 ^ { - 6 } \, \left ( \frac { \text{GeV} } { M _ { D M } } \right ) \,.$$

In particular, for our benchmark value M DM = 220 TeV, one finds N f , obs N DM glyph[similarequal] 5 × 10 -13 , the value indicated in Figs. 2-4 that we now briefly discuss.

## 4.1 Dependence on the reheat temperature

In Fig. 2 we show the dependence of N N DM ( z ) on the reheat temperature both for an initial thermal N S -abundance (upper panel) and for an initial vanishing N S -abundance (lower panel). In particular, we show N N DM ( z ) for different values of T RH as indicated. Notice that the value τ DM = 3 46 . × 10 28 s, corresponding to Λ ˜ = 1 13 . × 10 24 GeV, is just the value that reproduces the observed DM abundance for T RH = 10 15 GeV and M S = 300 GeV in the case of initial thermal N S -abundance. In the case of initial vanishing N S -abundance this value for τ DM is too high, i.e., the coupling too small, to get the correct relic abundance even for maximum allowed T RH . It should be immediately noticed that the LZ approximation overestimates by many orders of magnitude the relic DM abundance. It should also be noticed how in the case of initial vanishing N S -abundance the freeze-in temperature T f ∼ 10 9 GeV is much below the resonant temperature T res ∼ 10 12 GeV. 18 Another interesting thing to highlight is that for initial vanishing N S -abundance, the relic value is basically independent of T RH , except for the lowest value T RH = 10 10 GeV when the production occurs close to the freezing and the relic value is not fully saturated. Therefore, T RH ∼ 10 10 GeV should be regarded as a border line value such that below this value the production is strongly suppressed since there is no time for the asymmetry to be produced. These are all features that should be addressed by an analytical description.

17 Notice that with the normalisation we choose, one has for the photon abundance N f γ = 4 3, while / for the photon abundance at the present time one has N γ 0 = N f t , t f γ ( f 0 ) glyph[similarequal] 36 4. .

18 Using Eq. (24) one finds z res glyph[similarequal] 1 5 . × 10 7 and since the decay rate in this case is dominated by four body decays, from Eq. (38) one finds (Λ ˜ DM / ˜ Λ) 2 ∼ 10 , translating into a final DM abundance in the case 7 of LZ approximation that is seven orders of magnitude higher than the observed one, as it can be noticed in the plot.

Figure 2: Evolution of the DM abundance N N DM for different values of T RH , as indicated, and for fixed values of M DM , M S and τ DM . The upper (lower) panel is for initial thermal (vanishing) N S -abundance, i.e., N in N S = 1 in the upper panel. The step function corresponds to the instantaneous LZ approximation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00521v2-figure-46e38f79-b645-4402-bf12-c1fef7bf27da.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00521v2-figure-44d25156-2b72-4bac-a436-efd1c4d47001.jpg)

## 4.2 Dependence on the lifetime

In Fig. 3, we fix the reheat temperature to the highest possible value, T RH = 10 15 GeV, and show how the evolution of the DM abundance depends on τ DM . Note how for increasing values of τ DM , corresponding to larger values of Λ, the relic DM abundance decreases and ˜ vice-versa. It can be noticed how in the case of initial vanishing N S -abundance, even for a very low (and excluded by experimental data) value τ DM = 10 24 s, corresponding to Λ = 6 4 ˜ . × 10 21 GeV, the relic DM abundance is a few orders of magnitude below the measured value. Notice also how, though the final relic DM abundance is clearly strongly depending on τ DM , the freeze-in temperature is not. It should be noticed that also in this case, for all values of τ DM , the LZ approximation overestimates the abundance by about seven orders of magnitude.

## 4.3 Dependence on M S

Finally, in Fig. 4, we fix the value of T RH = 10 15 GeV and τ DM = 3 46 . × 10 28 s while we show different evolutions of the N DM -abundance for different values of M DM /M S , or equivalently M S , explicitly: M DM /M S = 2 2 . × 10 5 , 10 3 , 10 2 , 10 1 1 1 01 , . , . corresponding respectively to M S = 1GeV 0 22TeV 2 2TeV 22TeV 200TeV 218TeV. , . , . , , , Notice again how, in the case of initial vanishing N S -abundance, there is no value of M DM /M S for which the relic N DM -abundance can reproduce the observed DM abundance, contrary to the case of initial thermal N S -abundance. In the case of initial vanishing N S -abundance and for M &gt; M S W , the relic N DM abundance does not depend on M DM /M S as far as M DM /M S glyph[greatermuch] 1, similarly to the independence of T RH shown in Fig. 3. However, in the quasi-degenerate limit, for M DM glyph[similarequal] M S , there is an increase of about one order of magnitude until full saturation. This is clearly in stark contrast with the LZ approximation, where increasing the value of M DM /M S corresponds to an increased value of T res and of the mixing angle. This translates into an increase of the relic abundance despite the fact that for higher temperature the value of the N S -abundance at the resonance decreases. In the case of density matrix equation solutions, the freeze-in temperature and the evolution N DM ( z ) is approximately independent of M DM /M S despite the fact that the N N S is not, something that suggests that there is a compensation between higher mixing angle but smaller N S -abundance for higher values of M DM /M S . Of course this compensation is absent assuming initial thermal N S -abundance since this stays constant for z glyph[lessmuch] 1 and in this way the relic abundance increases for increasing M DM /M S .

Notice also how in the quasi-degenerate limit, for M DM /M S → 1, the case as originally proposed in [1], the result from the LZ approximation tends toward the relic N DM abundance

Figure 3: Evolution of the DM abundance N N DM for different values of τ DM , as indicated, and for fixed values of M DM , M S and T RH . As in the previous figure, the upper (lower) panel is for initial thermal (vanishing) N S -abundance and the step functions describe the instantaneous LZ approximation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00521v2-figure-ebe3d198-59a4-4813-8772-a204c7ab7e72.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00521v2-figure-35cb4d9e-a706-4414-9c0a-0f224697a2ff.jpg)

from the solution of density matrix equations though it is still two orders of magnitude higher, both in the case of initial thermal and vanishing N S -abundance. Finally, let us discuss the interesting case M S = 1GeV (red lines). In the hierarchical limit, for M DM /M S glyph[greatermuch] 1, and in the LZ approximation the resonant temperature grows to very large values. In the case of initial vanishing N S -abundance, the upper bound on the reheat temperature translates into the upper bound on M DM Eq. (46). However, one can see that from the numerical solutions of the density matrix equation there is no resonant temperature and actually most of the asymmetry is produced prior to the freeze-in temperature that is independent of M DM /M S . In this way the upper bound Eq. (46) does not actually hold, and one can both lower M S and increase M DM in a way to suppress the four body decay rate for M &lt; M S W . When this happens, the same lifetime is obtained for a much lower value of ˜ Λ (or equivalently higher value of the coupling λ DM S -) and this is why for M S = 1GeV one can see that the relic abundance greatly increases. For M DM = 220 TeV this is still not enough to reproduce the observed DM abundance in the case of initial vanishing N S -abundance. However, as we will see, an allowed region at high values M DM glyph[greaterorsimilar] 20 PeV, for M S = 1GeV, opens up also for vanishing initial N S -abundance.

## 4.4 Oscillations of the RH neutrino DM abundance prior to the freeze-in

In Fig. 5, we show a log-linear plot of N DM for initial thermal N S -abundance, highlighting the oscillations of the DM abundance prior to the freeze-in. In this case, we show the evolution N N DM ( z ) for fixed M S = 300 GeV but for three different choices of T RH and τ DM , in a way that the observed DM abundance is reproduced in all three cases. As one can see, we still choose the benchmark values T RH = 10 15 GeV and τ DM = 3 46 . × 10 28 s. The vertical line indicates the resonant temperature within the LZ approximation, and it can be noticed again how this is much higher than the freeze-in temperature. For this reason, in the density matrix equation solution, though the production is much less efficient than in the LZ approximation, this occurs at much lower temperatures and it allows to increase the scale of DM and lower the scale of the source RH neutrino partly compensating the reduced efficiency.

Figure 4: Evolution of the DM abundance N N DM for different values of M DM /M S , as indicated, and for fixed values of T RH and τ DM . As in previous figures, the upper (lower) panel is for initial thermal (vanishing) N S abundance and the step functions correspond to the the instantaneous LZ approximation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00521v2-figure-227b5c74-8ffe-43fb-8f9e-776e09ca98dd.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00521v2-figure-11faebb7-81dc-45c9-bae6-dcdc770b4877.jpg)

z

Figure 5: Evolution of the DM abundance N N DM for three different choices of T RH and τ DM as indicated in a linear plot for the abundance. For the three cases A, B and C one has respectively Λ = 1 13 ˜ . × 10 24 GeV Λ = 7 4 , ˜ . × 10 23 GeV and Λ = 3 55 ˜ . × 10 23 GeV.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00521v2-figure-5abf4d00-5af7-4d08-b54d-f130dd9b4dbc.jpg)

## 4.5 Analytical insight

Let us finally provide some analytical insight on the results we obtained, aiming especially at explaining why the LZ approximation fails, overestimating by many orders of magnitude the relic dark matter abundance when M DM glyph[greatermuch] M S . 19 As one can clearly see from Fig. 5, the problem with the LZ approximation is that it requires that neutrino oscillations have already developed at the resonance, while this is clearly not true. The DM production can then be explained in terms of neutrino oscillations occurring below the resonance, when the thermal mass from Yukawa interactions can be neglected in the diagonal terms in the Hamiltonian, but still with a time-dependent mixing angle. 20 An exhaustive description will be given in a forthcoming paper [45]. Here we just notice that the LZ approximation can be used reliably only imposing the condition that at the resonance neutrino oscillations have already developed , implying z res glyph[greatermuch] z osc , where z osc is the value of z corresponding to an age of the Universe t osc = 12 πT osc / ∆ M 2 , where t osc is defined as the time when the first dip occurs. Using t osc = 1 (2 / H z ( osc )), one finds z osc glyph[similarequal] 11 M DM / M ( Pl ∆ M 2 ) 1 3 / . For the example of Fig. 5, one indeed correctly finds z osc glyph[similarequal] 3 × 10 -4 . Imposing z res glyph[greatermuch] z osc , one finds for the validity of the LZ approximation a very stringent constraint M DM -M S glyph[lessmuch] 10 -12 GeV. This is consistent with the results of Fig. 4 where the LZ limit is approached but never recovered, since for M DM = 220 TeV one would need M DM /M S -1 glyph[lessmuch] 10 17 .

This simple analytical insight allows some interesting considerations:

- · The DM production is well explained by simple two neutrino mixing below the resonance with a rapidly decreasing mixing angle in a way that an abundance of source RH neutrino oscillates into a DM neutrino abundance but this does not oscillate back. This simple mechanism is quite a novel simple production mechanism.
- · As in the case of traditional non-resonant active-sterile neutrino oscillations in the early Universe, the monocromatic approximation we adopted is expected to work quite well and full momentum dependent description is expected to produce just a small correction.
- · An interesting point is that since the production is non-resonant anyway, the case M DM glyph[lessmuch] M S should work similarly but this time with even less constrained value of

19 Some of these considerations were anticipated in [44].

20 If we again make an analogy with left handed neutrino mixing and in particular with solar neutrino oscillations, below the resonance one would recover vacuum neutrino oscillations. In our case thermal effects from the Anisimov operator are still important in generating the time-dependent mixing angle and so it would not be correct to say that one recovers 'vacuum oscillations' but the analogy can be helpful to understand what happens.

∆ M 2 and, interestingly, traditional two RH neutrino seesaw models and leptogenesis with hierarchical masses might potentially work in explaining neutrino masses and mixing and matter-antimatter asymmetry. This scenario is very interesting and will be explored in a forthcoming paper [45].

## 4.6 Unifying dark matter and leptogenesis

Within the Higgs induced RHiNo DM model, the explanation of the DM abundance can be combined with an explanation of the matter-antimatter asymmetry within leptogenesis, obtaining a unified picture of neutrino masses, dark matter and leptogenesis [2]. In this case, the source RH neutrino should interfere with a third RH neutrino species and they should be quasi-degenerate in order for the CP asymmetry to be resonantly enhanced and have successful leptogenesis much below the lower bound of 10 10 GeV holding in the hierarchical case. The observed baryon-to-photon ratio is given by [28]

$$\eta _ { B 0 } ^ { o b s } = ( 6. 1 2 \pm 0. 0 4 ) \times 1 0 ^ { - 1 0 } \,.$$

Using the same normalisation as for the RH neutrino abundances, the final B -L asymmetry is related to the baryon-to-photon ratio predicted by leptogenesis simply by η lep B 0 glyph[similarequal] 0 01 . N f B -L , so that N f , obs B -L glyph[similarequal] 6 1 . × 10 -8 is the final N B -L value needed to reproduce the observed value of η B 0 .

The evolution of the B -L asymmetry with temperature can be calculated as the sum of six contributions both on the two heavy neutrino flavours, the source and the interfering RH neutrinos, and on the three charged lepton flavours considering that the asymmetry will be generated in the three flavoured regime. We can then write [47, 48, 49, 50]

$$N _ { B - L } ( z ) = \sum _ { \alpha = e, \mu, \tau } \left ( N _ { \Delta _ { \alpha } } ^ { ( 1 ) } ( z ) + N _ { \Delta _ { \alpha } } ^ { ( 2 ) } ( z ) \right ) \,,$$

where N (1) ∆ α ( z ) and N (2) ∆ α ( z ) are the abundances of the flavoured asymmetries ∆ α ≡ B/ 3 -L α generated by the lightest and next-to-lightest RH neutrino (the source RH neutrino can be either one or the other). The flavoured asymmetries can be calculated as ( I = 1 2) ,

$$N _ { \Delta _ { \alpha } } ^ { ( I ) } ( z ) = \varepsilon _ { I \alpha } \, \kappa _ { I \alpha } ( z, K _ { I }, K _ { 1 \alpha } + K _ { 2 \alpha } ) \,,$$

where ε Iα and K Iα are respectively the CP flavoured asymmetry and the flavoured decay parameter associated to the RH neutrino N I , while κ Iα ( z, K 1 α + K 2 α ) is the efficiency factor at temperature T = M /z S S and an analytical solution of the Boltzmann equation

gives

$$\circ \cdot \cdot \cdot \cdot \cdot \\ \kappa _ { I \alpha } ( z, K _ { I }, K _ { 1 \alpha } + K _ { 2 \alpha } ) \ = \ & \ \int _ { \stackrel { z _ { S } } { z _ { S } ^ { \ln } } } ^ { z _ { S } } d z _ { I } ^ { \prime } \, ( D _ { I } + S _ { I } ) \, \left [ N _ { N _ { I } } ( z _ { S } ^ { \prime } ) - N _ { N _ { I } } ^ { \text{eq} } ( z _ { S } ^ { \prime } ) \right ] \\ & \quad \times \exp \left [ - \int _ { z _ { S } ^ { \ln } } ^ { z _ { S } ^ { \prime } } d z _ { S } ^ { \prime \prime } \, W _ { \alpha } ( z _ { S } ^ { \prime \prime } ) \right ] \,, \\ \cdot \quad \ \ \. \ \. \ \. \ \. \ \. \ \. \ \. \ \. \ \. \ \. \ \. \ \.$$

where remember that z S = z M DM /M S and where W z α ( S ) is the wash-out term acting on flavour α including inverse decays and ∆ L = 1 scatterings. Notice that in D I and S I we are implying a dependence on the total decay parameter K I = ∑ α K Iα . We refer the reader to Appendix C for more details on the calculation of the asymmetry and expressions of flavoured CP asymmetries and flavoured decay parameters. Here we just notice that the asymmetry N B -L ( z ) depends on the low energy neutrino parameters, including the low energy phases, the degeneracy δ lep ≡ | M 1 -M /M 2 | 1 and one complex angle in the orthogonal matrix that parameterises the Dirac neutrino mass matrix.

For a specific choice of these parameters that satisfy successful leptogenesis and such that the asymmetry is dominantly produced by the decays of the source RH neutrinos (see Appendix C for details), we plot in Fig. 6 the evolution of the B -L asymmetry N B -L ( z ). We also plot N N DM ( z ) for a choice of values of the parameters in the Higgs induced RHiNo DM scenario that also reproduce the correct observed DM abundance (the same values as in case A in Fig 5).

## 5 Bounds on the DM mass

In Fig. 7 we summarise the results we found for different choices of the parameters plotting the allowed regions in the M DM -τ DM plane. In the higher panel we imposed the most conservative upper bound T RH &lt; 10 15 GeV, in the central panel we imposed T RH &lt; 10 12 GeV and finally in the bottom panel we set more stringently T RH &lt; 10 10 GeV. In the upper panel, for T RH &lt; 10 15 GeV, one can see how the only way to have an allowed region for vanishing initial N S -abundance is for M &lt; M S W and in particular we show the allowed region for M &gt; S 1 GeV. As we discussed, this has the effect to suppress the four body decay rate, nullifying the upper bound on M DM . However, one can see how in this case there is a lower bound M DM glyph[greaterorsimilar] 20 PeV.

On the other hand, for initial thermal N S -abundance, values for M &gt; M S W and even values M &gt; S 300 GeV compatible with a traditional scenario of leptogenesis from decays, allowed regions exist. In particular for M &gt; S 300 GeV, one has 0 5 TeV . glyph[lessorsimilar] M DM glyph[lessorsimilar] 0 5 PeV . and lifetimes as large as 10 31 s are allowed. One should appreciate how improvement in

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00521v2-figure-07152470-a689-4a74-9522-841a8e4fea2f.jpg)

z

Figure 6: Evolution of the B -L asymmetry and DM abundance N N DM for a choice of parameters such that the final values simultaneously reproduce the observed values of matter-antimatter asymmetry (see Appendix C) and DM abundance. The N S -abundance is also shown and one can see that the plot is for initial thermal N S -abundance.

the lower bound on τ DM from neutrino telescope experiments will progressively test the scenario placing more and more stringent constraints.

In the central panel, for T RH &lt; 10 12 GeV, there is no significant reduction of the allowed regions and this is in line with what we noticed in Fig. 3: most of the DM abundance is produced prior to the freezing at T f ∼ 10 9 GeV and therefore only when T RH gets closer to 10 9 GeV one has a noticeable reduction of the relic DM abundance.

And indeed one can see that in the lower panel, for T RH &lt; 10 10 GeV, all allowed regions shrink considerably and, in particular, there is no allowed region for initial N S -abundance even for M &gt; S 1 GeV. This stringent upper bound on the reheat temperature might be motivated for example by a supersymmetric version of the scenario, requiring an avoidance of the gravitino problem [51].

## 6 Final discussion

We studied the production of the DM abundance within the Higgs induced RHiNo DM model, solving numerically density matrix equation. The results show that the LZ approximation overestimates the DM abundance by many orders of magnitude. In

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00521v2-figure-cf27b3e9-32be-4ff9-9e1d-375aa861670f.jpg)

M

DM

/

GeV

Figure 7: Allowed regions in the plane τ DM versus M DM for different conditions as indicated. The three panels correspond to three different choices of the upper bound on T RH . The red star and the contour lines around are the same as in Fig. 1.

the quasi-degenerate limit the mismatch is minimum but still the DM abundance is overestimated by two orders of magnitude. It is then clear that a calculation employing density matrix equation is crucial. Moreover the DM production occurs at temperatures much below the resonant temperature and this allows to open solutions for low values of the source RH neutrino mass. We have seen that in this way solutions for initial vanishing N S -abundance are still possible but only for M S below the W boson mass and with a stringent lower bound on M DM . In particular, imposing M &gt; S 1 GeV, we obtained M DM glyph[greaterorsimilar] 20 PeV. In this case, one cannot reproduce the matter-antimatter asymmetry within traditional leptogenesis from decays, but it poses the question whether this can be achieved considering leptogenesis from RH neutrino mixing that works indeed for GeV RH neutrino masses [37].

If one wants M DM ∼ 100 TeV, as IceCube data seem to favour, these results then motivate the possibility to consider processes able to thermalise the source RH neutrino prior to the DM abundance freezing. The existence of such processes is certainly plausible if one thinks that the non-renormalisable interactions in any case requires a UV-completion. One can for example think that the RH neutrinos at high temperatures might have extra-gauge interactions and can get produced by very heavy Z ′ bosons, a well known possibility [52]. However, even more interestingly, one can think that the Higgs induced interactions for the source RH neutrinos are actually much stronger than for the DM RH neutrino and able to thermalise the source RH neutrinos prior to the DM production. This possbility is quite attractive since it would not require additional interactions.

There is also another intriguing possibility emerging from our study. Within the LZ approximation it was necessary to impose M DM &gt; M S in order to have a resonance. However, the numerical solution of the density matrix equations show that the DM production is actually non-resonant. One can then wonder whether solutions with M DM &lt; M S might open up. These would be quite interesting since in this case the DM RH neutrino would be the lightest RH neutrino and one could embed the mechanism within traditional two RH neutrino high energy scale seesaw models.

Of course it would be also desirable to have an analytic understanding of our numerical results. 21 In particular, it would be quite useful to have an analytic expression for the final relic DM abundance and for the freeze-in temperature.

Our results also should be generalised taking into account the momentum distribution. However, since RH neutrinos do not contribute to the effective potentials, complicated back-reaction effects are excluded and including the momentum dependence should produce only corrections.

21 These have been cross-checked solving independent codes.

In conclusion, a density matrix calculation of the Higgs induced RHiNo DM relic abundance is certainly necessary and confirms that the mechanism can reproduce the observed DM abundance and simultaneously the matter-antimatter asymmetry within certain allowed regions in the space of parameters. However, it also paves the way for new interesting scenarios, motivating further investigation in different directions.

## Acknowledgments

We wish to thank Teppei Katori for useful discussions. PDB and YLZ acknowledge financial support from the STFC Consolidated Grant L000296/1. RS is supported by a Newton International Fellowship (NF 171202) from Royal Society (UK) and SERB (India). KF acknowledges financial support from the NExT/SEPnet Institute. This project has received funding/support from the European Union Horizon 2020 research and innovation programme under the Marie Skglyph[suppress] lodowska-Curie grant agreements number 690575 and 674896.

## Appendix A: Two body decay rate

In this Appendix we first derive the mixing angle induced by the Anisimov operator between the DM and the source RH neutrino at zero temperature Eq. (32) and responsible for the two body decay channel and then the resulting DM life time when two body decays dominate. After electroweak spontaneous symmetry breaking the Lagrangian Eq. (2) becomes

$$- \mathcal { L } _ { m } ^ { \ell + \nu } = \overline { \alpha _ { L } } \, m _ { \alpha } \, \alpha _ { R } + \overline { \nu _ { L \alpha } } \, m _ { D \alpha J } \, N _ { R J } + \frac { 1 } { 2 } \, \overline { N _ { R J } ^ { c } } \, M _ { J } \, N _ { R J } + \frac { 1 } { 2 } \, \overline { N _ { R J } ^ { c } } \, \delta M _ { I J } ^ { \lambda } \, N _ { R J } + \text{h.c.} \, \ \text{(A.1)}$$

where δM Λ = 2 λ IJ v / 2 Λ is the effective Majorana mass term correction generated by the Anisimov operator at zero temperature introducing off-diagonal terms in the total Majorana mass term

$$M _ { I J } ^ { \Lambda } & = M _ { I I } \, \delta _ { I J } + \delta M _ { \Lambda } \,.$$

This can be (Takagi) diagonalised by a unitary matrix U Λ T R , in a way that

$$M _ { I J } ^ { \Lambda } = U _ { R } ^ { \Lambda T } \, D _ { M ^ { \Lambda } } \, U _ { R } ^ { \Lambda } \,,$$

where D M Λ ≡ diag( M ,M ′ 1 ′ 2 ), and N RI → N ′ RJ ′ = U † RJ I ′ N RI . Since we are in a two-neutrino mixing case this matrix can be made real and parameterised in terms of just one mixing angle, explicitly

$$U _ { R } ^ { \Lambda } ( \theta _ { \Lambda } ^ { 0 } ) = \left ( \begin{array} { c c } \cos \theta _ { \Lambda } ^ { 0 } & \sin \theta _ { \Lambda } ^ { 0 } \\ - \sin \theta _ { \Lambda } ^ { 0 } & \cos \theta _ { \Lambda } ^ { 0 } \end{array} \right ) \,.$$

The mixing angle can be easily calculated from (A.3). The correction to the masses is negligible and one can approximate M ′ I glyph[similarequal] M I . The important point is that in the new primed basis of mass eigenstates the neutrino Dirac mass matrix becomes m ′ DαI ′ = m DαJ U RJI ′ and so using the parameterisation (A.5) one finds

$$m _ { D \alpha I ^ { \prime } } ^ { \prime } \simeq \left ( \begin{array} { c c c } - \theta _ { \Lambda } ^ { 0 } \, m _ { D e S } & \, m _ { D e S } \\ - \theta _ { \Lambda } ^ { 0 } \, m _ { D \mu S } & \, m _ { D \mu S } \\ - \theta _ { \Lambda } ^ { 0 } \, m _ { D \tau S } & \, m _ { D \tau S } \end{array} \right ) \,,$$

showing that the Higgs induced interactions generate small effective Yukawa couplings in the DM mass eigenstate that induce eventually its decays with a decay rate Γ DM → A + ν S = θ 2 Λ0 h 2 S M DM / (4 π ), that of course vanishes in the limit Λ ˜ →∞ .

## Appendix B: Equivalent forms for the density matrix equation

Let us show how from Eq. (49) for the density matrix equation, one can easily arrive to the equivalent form Eq. (54). First, using Eq. (50), we can express the anti-commutator as

$$- i \left [ \mathcal { H }, N \right ] _ { I J } = \frac { P _ { 0 } } { 2 } \left ( \vec { V } \times \vec { P } \right ) \cdot \vec { \sigma } \,.$$

Defining then R ≡ (Γ D + Γ )( S N N S -N eq N S ), the second term in the right-hand side of Eq. (49) can be recast as

$$\begin{pmatrix} 0 & \frac { 1 } { 2 } ( \Gamma _ { D } + \Gamma _ { S } ) \, N _ { \text{DM-S} } \\ \frac { 1 } { 2 } ( \Gamma _ { D } + \Gamma _ { S } ) \, N _ { \text{S-DM} } & ( \Gamma _ { D } + \Gamma _ { S } ) \, ( N _ { N _ { S } } - N _ { N _ { S } } ^ { \text{eq} } ) \end{pmatrix} = \frac { P _ { 0 } } { 4 } ( \Gamma _ { D } + \Gamma _ { S } ) \, \vec { P } _ { T } \cdot \vec { \sigma } + \begin{pmatrix} 0 & 0 \\ 0 & R \end{pmatrix} \, \text{ (B.2)}$$

and also one can write

$$\begin{pmatrix} 0 & 0 \\ 0 & R \end{pmatrix} = \frac { R } { 2 } \, I - \frac { R } { 2 } \, \sigma _ { z } \,.$$

With some straightforward steps one then arrive to

$$\frac { d \vec { P } } { d t } \cdot \vec { \sigma } = ( \vec { V } \times \vec { P } ) \cdot \vec { \sigma } - \frac { 1 } { 2 } \left ( \Gamma _ { D } + \Gamma _ { S } \right ) \vec { P } _ { T } \cdot \vec { \sigma } + \frac { 1 } { P _ { 0 } } \, \left ( R - \frac { d P _ { 0 } } { d t } \right ) - \frac { 1 } { P _ { 0 } } \, \left ( R \, \sigma _ { z } + \frac { d P _ { 0 } } { d t } \, \vec { P } \cdot \vec { \sigma } \right ) \,,$$

implying R = dP /dt 0 and Eq. (54).

Finally, it is easy to see that Eq. (49) is fully equivalent to another popular form for the density matrix equation expressed in terms of anti-commutators [53, 37]

$$\frac { d N _ { I J } } { d t } = - i \left [ \mathcal { H }, N \right ] _ { I J } - \frac { 1 } { 2 } \ \{ \Gamma, N \} _ { I J } + \frac { 1 } { 2 } \ \{ \Gamma, N ^ { \text{eq} } \} _ { I J }$$

where simply

$$\Gamma = \begin{pmatrix} 0 & 0 \\ 0 & \Gamma _ { D } + \Gamma _ { S } \end{pmatrix} \quad \text{and} \quad N ^ { \text{eq} } = \begin{pmatrix} N _ { \text{DM} } ^ { \text{eq} } & 0 \\ 0 & N _ { S } ^ { \text{eq} } \end{pmatrix} \.$$

## Appendix C: leptogenesis with two quasi-degenerate RH neutrinos

In this Appendix we give more details on the calculation of the B -L asymmetry within a two quasi-degenerate RH neutrino scenario and the specific case shown in Fig. 6. The flavoured decay parameters are defined as

$$K _ { I \alpha } \equiv \frac { \Gamma _ { I \alpha } + \overline { \Gamma } _ { I \alpha } } { H ( T = M _ { I } ) } = \frac { | m _ { D \alpha I } | ^ { 2 } } { M _ { I } m _ { * } } \,,$$

where the equilibrium neutrino mass

$$m _ { * } \equiv \frac { 1 6 \, \pi ^ { 5 / 2 } \, \sqrt { g _ { * } ^ { \text{SM} } } } { 3 \, \sqrt { 5 } } \, \frac { v ^ { 2 } } { M _ { \text{Pl} } } \simeq 1. 1 \, \text{meV} \,.$$

The flavoured CP asymmetries are defined as

$$\varepsilon _ { I \alpha } \equiv - \frac { \Gamma _ { I \alpha } - \overline { \Gamma } _ { I \alpha } } { \Gamma _ { I } + \overline { \Gamma } _ { I } } \simeq \frac { \overline { \varepsilon } ( M _ { I } ) } { K _ { I } } \, \left \{ \mathcal { I } ^ { \alpha } _ { I J } \, \xi ( M ^ { 2 } _ { J } / M ^ { 2 } _ { I } ) + \mathcal { J } ^ { \alpha } _ { I J } \, \frac { 2 } { 3 ( 1 - M ^ { 2 } _ { I } / M ^ { 2 } _ { J } ) } \right \} \,,$$

where

$$\bar { \varepsilon } ( M _ { I } ) \equiv \frac { 3 } { 1 6 \, \pi } \, \left ( \frac { M _ { I } \, m _ { \text{atm} } } { v ^ { 2 } } \right ) \simeq 1. 0 \times 1 0 ^ { - 1 3 } \, \left ( \frac { M _ { I } } { \text{TeV} } \right ) \,,$$

while I α IJ and J α IJ originate from the interference of tree level with one loop self-energy and vertex-diagrams and are given by

$$\mathcal { I } ^ { \alpha } _ { I J } \equiv \frac { \text{Im} \left [ m _ { D \alpha I } ^ { * } \, m _ { D \alpha J } \, ( m _ { D } ^ { \dagger } \, m _ { D } ) _ { I J } \right ] } { M _ { I } \, M _ { J } \, m _ { \text{atm} } \, m _ { * } } \, \quad \mathcal { J } ^ { \alpha } _ { I J } \equiv \frac { \text{Im} \left [ m _ { D \alpha I } ^ { * } \, m _ { D \alpha J } \, ( m _ { D } ^ { \dagger } \, m _ { D } ) _ { J I } \right ] } { M _ { I } \, M _ { J } \, m _ { \text{atm} } \, m _ { * } } \, \frac { M _ { I } } { M _ { J } } \,. \ \text{ (C.5)}$$

In the case of two RH neutrinos with quasi-degenerate masses, Eq. (60) gives for the final B -L asymmetry [2]

$$N _ { B - L } ^ { f } \simeq \frac { \bar { \varepsilon } ( M _ { 1 } ) } { 3 \, \delta _ { \text{lerp} } } \left ( \frac { 1 } { K _ { 1 } } + \frac { 1 } { K _ { 2 } } \right ) \, \sum _ { \alpha } \, \kappa ^ { f } ( K _ { 1 \alpha } + K _ { 2 \alpha } ) \, \left [ \mathcal { I } _ { 1 2 } ^ { \alpha } + \mathcal { J } _ { 1 2 } ^ { \alpha } \right ] \,,$$

where κ f is the final value of the efficiency factor that can be obtained from Eq. (62) in the limit z S → ∞ (see discussion below for an analytical expression). In order to take into account the low energy neutrino experimental data, it is convenient to introduce the orthogonal parameterisation for the neutrino Dirac mass matrix [54],

$$m _ { D } = U \, \sqrt { D _ { m } } \, \Omega \, \sqrt { D _ { M } } \,.$$

Since in our case m 1 = 0, one has m 2 = m sol glyph[similarequal] 8 6 meV and . m 3 = m atm glyph[similarequal] 50 meV, where m sol and m atm are the solar and atmospheric neutrino mass scales respectively and

we consider only normal hierarchy since this is now quite strongly favoured by the global analyses compared to inverted hierarchy [56].

The leptonic mixing matrix can be parameterised in terms of three mixing angles and three CP violating phases as

$$U = \left ( \begin{array} { c c c } c _ { 1 3 } c _ { 1 2 } & c _ { 1 3 } s _ { 1 2 } & s _ { 1 3 } e ^ { - i \delta } \\ - c _ { 2 3 } s _ { 1 2 } - s _ { 2 3 } s _ { 1 3 } c _ { 1 2 } e ^ { i \delta } & c _ { 2 3 } c _ { 1 2 } - s _ { 2 3 } s _ { 1 3 } s _ { 1 2 } e ^ { i \delta } & s _ { 2 3 } c _ { 1 3 } \\ s _ { 2 3 } s _ { 1 2 } - c _ { 2 3 } s _ { 1 3 } c _ { 1 2 } e ^ { i \delta } & - s _ { 2 3 } c _ { 1 2 } - c _ { 2 3 } s _ { 1 3 } s _ { 1 2 } e ^ { i \delta } & c _ { 2 3 } c _ { 1 3 } \end{array} \right ) \, \text{diag} ( e ^ { i \rho }, \, e ^ { i \sigma }, \, 1 ) \, \, \text{ (C.S)}$$

where c ij = cos θ ij , s ij = sin θ ij , δ is the Dirac phase and ρ and σ are the two Majorana phases. The best fit values and the 1 σ (3 σ ) confidence level (C.L.) ranges of the reactor, solar, atmospheric mixing angles and Dirac phase for normal hierarchy are given by [56]

$$\begin{array} { r c l c r c l } \psi & \psi & \psi & \psi & \psi & \psi & \psi \\ \theta _ { 1 3 } & = & 8. 6 1 ^ { \circ + 0. 1 3 ^ { \circ } } _ { - 0. 1 3 ^ { \circ } } \left ( 8. 2 2 ^ { \circ } - 8. 9 9 ^ { \circ } \right ), \\ \theta _ { 1 2 } & = & 3 3. 8 2 ^ { \circ + 0. 7 8 ^ { \circ } } _ { - 0. 7 6 ^ { \circ } } \left ( 3 1. 6 1 ^ { \circ } - 3 6. 2 7 ^ { \circ } \right ), \\ \theta _ { 2 3 } & = & 4 8. 3 ^ { \circ + 1. 1 ^ { \circ } } _ { - 1. 9 ^ { \circ } } \left ( 4 0. 8 ^ { \circ } - 5 1. 3 ^ { \circ } \right ), \\ \delta & = & - 1 3 8 ^ { \circ + 3 8 ^ { \circ } } _ { - 2 8 ^ { \circ } } \left ( - 1 5 9 ^ { \circ } - + 1 0 ^ { \circ } \right ), \end{array}$$

while there is no experimental information on the Majorana phases. Coming finally to the orthogonal matrix, in the two RH neutrino case this is given by

$$\Omega = \left ( \begin{array} { c c c } 0 & 0 & 1 \\ \cos \omega & \sin \omega & 0 \\ - \zeta \sin \omega & \zeta \cos \omega & 0 \end{array} \right ) \,,$$

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 [55].

In terms of the orthogonal parameterisation, the flavour decay parameters, in the considered two quasi-degenerate RH neutrino case, can be written as

$$K _ { I \alpha } = \left | \sum _ { j } \, \sqrt { \frac { m _ { j } } { m _ { * } } } \, U _ { \alpha j } \, \Omega _ { j I } \right | ^ { 2 } = \left | \sqrt { K _ { s o l } } \, U _ { \alpha 2 } \, \Omega _ { 2 I } + \sqrt { K _ { \text{atm} } } \, U _ { \alpha 3 } \, \Omega _ { 3 I } \right | ^ { 2 } \,,$$

where we defined K sol ≡ m sol /m glyph[star] glyph[similarequal] 8 and K atm glyph[similarequal] m atm /m glyph[star] glyph[similarequal] 45. Using Eq. (C.11) we can then write them in terms of the complex angle, obtaining

$$K _ { 1 \alpha } \ = \ \left | \sqrt { K _ { s o l } } \, U _ { \alpha 2 } \, \cos \omega - \zeta \, \sqrt { K _ { \text{atm} } } \, U _ { \alpha 3 } \, \sin \omega \right | ^ { 2 } \,,$$

$$K _ { 2 \alpha } \ = \ \left | \sqrt { K _ { \text{sol} } } \, U _ { \alpha 2 } \, \sin \omega + \zeta \, \sqrt { K _ { \text{atm} } } \, U _ { \alpha 3 } \, \cos \omega \right | ^ { 2 } \,.$$

We can also re-write I α 12 and J α 12 in terms of the orthogonal parametrisation, obtaining

$$\mathcal { I } _ { 1 2 } ^ { \alpha } \ = \ \sum _ { k, l, m = 2, 3 } \frac { m _ { k } \sqrt { m _ { l } \, m _ { m } } } { m _ { \text{atm} } \, m _ { * } } \, \text{Im} \left [ U _ { \alpha m } ^ { * } \, U _ { \alpha l } \, \Omega _ { m 1 } ^ { * } \, \Omega _ { l 2 } \, \Omega _ { k 1 } ^ { * } \, \Omega _ { k 2 } \right ]$$

$$\mathcal { J } ^ { \alpha } _ { 1 2 } \ = \ \sum _ { k, l, m = 2, 3 } ^ { k, l, m = 2, 3 } \, \frac { m _ { k } \sqrt { m _ { l } \, m _ { m } } } { m _ { \text{atm} } \, m _ { * } } \, \text{Im} \left [ U _ { \alpha m } ^ { * } \, U _ { \alpha l } \, \Omega _ { m 1 } ^ { * } \, \Omega _ { l 2 } \, \Omega _ { k 2 } ^ { * } \, \Omega _ { k 1 } \right ] \.$$

Therefore, one can see that the final asymmetry will depend on two light neutrino mass scales, five parameters in U (1 Majorana phase can be removed in the two RH neutrino case), on the common heavy neutrino mass scale M S , on the relative mass splitting δ lep and on the complex angle ω (in total eleven parameters). In order to obtain a specific solution, that we show in Fig. 6, we have taken a real U , and choose Re[ ω ] = 0 and Im[ ω ] glyph[lessmuch] 1. 22 In this case one obtains

$$\mathcal { I } _ { 1 2 } ^ { \alpha } \simeq \mathcal { J } _ { 1 2 } ^ { \alpha } \simeq \text{Im} [ \omega ] \, U _ { \alpha 2 } \, U _ { \alpha 3 } \, K _ { \text{sol} } \, \left ( \sqrt { \frac { K _ { \text{atm} } } { K _ { \text{sol} } } } + \sqrt { \frac { K _ { \text{sol} } } { K _ { \text{atm} } } } \right ) \,.$$

and

$$K _ { 1 \alpha } \simeq K _ { \text{sol} } U _ { \alpha 2 } ^ { 2 } \,, \ \ K _ { 2 \alpha } \simeq K _ { \text{atm} } U _ { \alpha 3 } ^ { 2 } \,,$$

implying K 1 glyph[similarequal] K sol and K 2 glyph[similarequal] K atm . In this way the expression for the final asymmetry Eq. (C.6) gets specialised into

$$N _ { B - L } ^ { f } & \simeq \frac { 2 } { 3 } \, \frac { \text{Im} [ \omega ] \, \bar { \varepsilon } ( M _ { S } ) } { \delta _ { \text{lep} } } \left ( 1 + \frac { K _ { \text{sol} } } { K _ { \text{atm} } } \right ) \left ( \sqrt { \frac { K _ { \text{atm} } } { K _ { \text{sol} } } } + \sqrt { \frac { K _ { \text{sol} } } { K _ { \text{atm} } } } \right ) \sum _ { \alpha } \kappa ^ { f } ( K _ { 1 \alpha } + K _ { 2 \alpha } ) \, U _ { \alpha 2 } \, U _ { \alpha 3 } \,. \\ \kappa ^ { \quad \text{-} } \, \lambda ^ { \quad \text{-} } \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdots$$

In Fig. 6 we have taken the best fit values for the three mixing angles M S = 1TeV, Im[ ω ] = 0 01 and . δ lep glyph[similarequal] 4 × 10 -10 , the correct values to reproduce N f , obs B -L glyph[similarequal] 6 1 . × 10 -8 . Notice that of course one could also have chosen a real ω and in that case one would have obtained a solution where all asymmetry stems from low energy phases, but in that case one needs Re[ ω ] glyph[greaterorsimilar] 1. However, also in the case of a purely imaginary Ω that we are considering, despite one has a vanishing total CP asymmetry, the final asymmetry does not, since the different washout in the different flavours prevents a full cancellation of the flavoured asymmetries.

22 This is actually a limit that is well motivated within models with discrete flavour symmetries: when the symmetry is conserved the orthogonal matrix coincides exactly with the permutation matrix and there is no CP violation at high energies, when the symmetry is broken, the small breaking parameter can be identified with | ω | and is responsible for CP violation [57, 58].

Finally, let us briefly discuss the calculation of the efficiency factor at temperature T = z /M S S , introduced in Eq. (62), and its final value κ f . 23

The efficiency factor is related to the B -L asymmetry simply by Eq. (61). The expression (62) comes from the solution of the set of rate equations ( I = 1 2; , α = e, µ, τ )

$$\frac { d N _ { N _ { I } } } { d z _ { \underset { \iota r } { r } } } \ = \ - ( D _ { I } + S _ { I } ) \, ( N _ { N _ { I } } - N _ { N _ { I } } ^ { \text{eq} } )$$

$$\frac { d N _ { \Delta _ { \alpha } } ^ { ( I ) } } { d z _ { S } } \ = \ \varepsilon _ { I \alpha } \left ( D _ { I } + S _ { I } \right ) \left ( N _ { N _ { I } } - N _ { N _ { I } } ^ { e q } \right ) - W _ { \alpha } \, N _ { \Delta _ { \alpha } } ^ { ( I ) } \,,$$

where the source RH neutrino is either the lightest (S = 1) or the next-to-lightest (S = 2) one and where let us recall that z S ≡ M /T S = z M /M S DM . Notice that since the two RH neutrinos are quasi-degenerate, then z S glyph[similarequal] M /T 1 glyph[similarequal] M /T 2 . The equilibrium RH neutrino abundance N eq N I is given by Eq. (27) and it is actually the same for lightest ( I = 1) and next-to-lightest ( I = 2) RH neutrino since they are quasi-degenerate. The sum of the decay and ∆ L = 1 scattering terms can be written as

$$D _ { I } ( z _ { S } ) + S _ { I } ( z _ { S } ) = D _ { I } ( z _ { s } ) \, j ( z _ { S } ) \,,$$

where

$$D _ { I } ( z _ { S } ) = K _ { I } \, z _ { S } ^ { 2 } \frac { \int _ { z _ { S } } ^ { \infty } \, d x \, \sqrt { x ^ { 2 } - z _ { S } ^ { 2 } } \, e ^ { - x } } { \int _ { z _ { S } } ^ { \infty } \, d x \, x \, \sqrt { x ^ { 2 } - z _ { S } ^ { 2 } } \, e ^ { - x } } \\ \dots \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdots$$

accounts for decays while the function j ( z S ) accounts for ∆ L = 1 scatterings and can be written as

$$j ( z _ { S } ) \simeq \left [ \frac { 0. 2 } { z _ { S } } + \frac { z _ { S } } { a } \ln \left ( 1 + \frac { a } { z _ { S } } \right ) \right ] \, \left ( 1 + \frac { 2 } { z _ { S } } \right ) \,,$$

where a glyph[similarequal] 8 π / 2 [9 ln( M /M S H )] and M H glyph[similarequal] 0 4 . T is the thermal Higgs mass. Notice that j ( z S glyph[greatermuch] 1) glyph[similarequal] 1. The wash-out terms W Iα ( z S ), the sum of the inverse decays and ∆ L = 1 scattering wash-out terms, can be also written through the function j ( z S ) as

$$W _ { \alpha } ( z _ { S } ) = j ( z _ { S } ) \, W _ { \alpha } ^ { I D } ( z _ { S } ) \,,$$

where the flavoured inverse decay wash-out rate can be expressed as

$$W _ { \alpha } ^ { I D } ( z _ { \S } ) = \frac { 1 } { 2 } \left ( K _ { 1 \alpha } + K _ { 2 \alpha } \right ) z _ { \S } ^ { 2 } \, \int _ { z _ { \S } } ^ { \infty } \, d x \sqrt { x ^ { 2 } - z _ { \S } ^ { 2 } } \, e ^ { - x } \,.$$

The integral expression Eq. (62) for the efficiency factor is simply the solution, in the form of a Laplace integral, of the rate equations (C.21) describing the evolution of the flavoured

23 This discussion extends the results found in [47, 48, 49, 50] to the case of two quasi-degenerate RH neutrinos including flavour effects.

asymmetries. Using the rate equation for the RH neutrino abundance Eq. (C.20), Eq. (62) for the efficiency factor can be recast as

$$\kappa _ { I \alpha } ( z _ { S }, K _ { I }, K _ { 1 \alpha } + K _ { 2 \alpha } ) = - \int _ { z _ { S } ^ { \ln } } ^ { z _ { S } } \, d z _ { S } ^ { \prime } \, \frac { d N _ { N _ { I } } } { d z _ { S } } \, \exp \left [ - \int _ { z _ { S } ^ { \ln } } ^ { z _ { S } ^ { \prime } } \, d z _ { S } ^ { \prime \prime } \, W _ { \alpha } ( z _ { S } ^ { \prime \prime } ) \right ] \,.$$

The wash-out is active (i.e., W α &gt; 1) in an interval z α on glyph[greaterorsimilar] z S glyph[greaterorsimilar] z α off , where z α on glyph[similarequal] 2 / √ K 1 α + K 2 α glyph[lessmuch] 1. For z S glyph[lessorsimilar] z on the wash-out can be neglected and moreover decays are not fast enough to track the equilibrium distribution so that in this regime N N I -N eq N I glyph[similarequal] z / 2 S 4. Moreover, for z S glyph[lessmuch] 1, one has D I + S I glyph[similarequal] 0 2 . K I . In this case one obtains for the efficiency factor

$$k _ { I \alpha } ( z _ { S } \lesssim z _ { \text{on} } ^ { \alpha }, K _ { I } ) \simeq \frac { 0. 2 \, K _ { I } } { 1 2 } \, z _ { S } ^ { 3 }.$$

On the other hand, for z glyph[greaterorsimilar] z α on , since we are in the strong wash-out regime with K , K 1 2 glyph[greatermuch] 1, one can approximate dN N I /dz glyph[similarequal] dN eq N I /dz and, noticing that

$$\frac { d N _ { N _ { S } } ^ { \text{eq} } } { d z _ { S } } = - \frac { 2 \, W _ { \alpha } ^ { I D } } { \left ( K _ { 1 \alpha } + K _ { 2 \alpha } \right ) z _ { S } } \,,$$

one finds a very simple analytical expression

$$\kappa _ { I \alpha } ( z _ { S }, K _ { 1 \alpha } + K _ { 2 \alpha } ) \simeq \frac { 2 } { ( K _ { 1 \alpha } + K _ { 2 \alpha } ) \, \bar { z } \, j ( \bar { z } ) } \, \left ( 1 - e ^ { - \frac { ( K _ { 1 \alpha } + K _ { 2 \alpha } ) \, \bar { z } \, j ( \bar { z } ) } { 2 } \left ( 1 - N _ { N _ { I } } ^ { \alpha q } ( z _ { S } ) \right ) } \right ), \quad ( C. 3 0 )$$

where ¯ z = min( z, z B ). At large values of z S glyph[greatermuch] z in , this simply reduces to κ Iα glyph[similarequal] 2 / [( K 1 α + K 2 α ) ¯ z j (¯)]. z

A good interpolation between (C.28) and ( C. 30), working at all values of z S , is then given by

$$\kappa _ { I \alpha } ( z _ { \mathcal { S } }, K _ { I }, K _ { 1 \alpha } + K _ { 2 \alpha } ) \simeq \frac { 0. 2 \, K _ { I } } { 1 2 } \, z _ { \mathcal { S } } ^ { 3 } \, \left ( 1 + \frac { 0. 2 \, K _ { I } \, ( K _ { 1 \alpha } + K _ { 2 \alpha } ) } { 2 4 } \, z _ { \mathcal { S } } ^ { 3 } \, \bar { z } \, j ( \bar { z } ) \right ) ^ { - 1 } \,. \quad ( \text{C.31} )$$

For the final value, using Eq. (C.30), one has

$$\kappa _ { \alpha } ^ { f } ( K _ { 1 \alpha } + K _ { 2 \alpha } ) \simeq \frac { 2 } { ( K _ { 1 \alpha } + K _ { 2 \alpha } ) \, z _ { B } \, j ( z _ { B } ) } \, \left ( 1 - e ^ { - \frac { ( K _ { 1 \alpha } + K _ { 2 \alpha } ) \, z _ { B } \, j ( z _ { B } ) } { 2 } } \right ),$$

where

$$z _ { B } ( K _ { 1 \alpha } + K _ { 2 \alpha } ) \simeq 2 + 4 \, ( K _ { 1 \alpha } + K _ { 2 \alpha } ) ^ { 0. 1 3 } \, e ^ { - 2. 5 / ( K _ { 1 \alpha } + K _ { 2 \alpha } ) } \,.$$

Note that the final value is independent of K I and, therefore, it is the same for both RH neutrinos. Since in the strong wash-out regime z B glyph[greatermuch] 1, one can approximate j ( z B ) glyph[similarequal] 1, meaning that ∆ L = 1 scatterings give a small correction that can be neglected with good

approximation. This makes the solution quite stable under the inclusion of many different subtle and complicated effects, like for example a precise calculation of the thermal Higgs mass. Also let us notice that since we are in the fully three flavoured regime with masses M I glyph[lessmuch] 10 9 GeV and since, on the other hand, M I glyph[greatermuch] 100 GeV and the asymmetry is produced from decays when RH neutrino mixing gives a negligible contribution, then Boltzmann equations are expected to work very well and a density matrix formalism with more heavy neutrino flavours [42] is not expected to be necessary to describe leptogenesis, though a dedicated analysis might be interesting to confirm this expectation.

## References

- [1] A. Anisimov and P. Di Bari, Cold dark matter from heavy right-handed neutrino mixing , Phys. Rev. D 80 (2009) 073017, [arXiv:0812.5085 [hep-ph]];
- [2] P. Di Bari, P. O. Ludl and S. Palomares-Ruiz, Unifying leptogenesis, dark matter and high-energy neutrinos with right-handed neutrino mixing via Higgs portal , JCAP 1611 (2016) no.11, 044 [arXiv:1606.06238 [hep-ph]].
- [3] A. Anisimov, Majorana Dark Matter , doi:10.1142/9789812770288 0058, hep-ph/0612024;
- [4] F. Bezrukov, D. Gorbunov and M. Shaposhnikov, On initial conditions for the Hot Big Bang , JCAP 0906 (2009) 029 [arXiv:0812.3622 [hep-ph]].
- [5] B. Patt and F. Wilczek, Higgs-field portal into hidden sectors , hep-ph/0605188.
- [6] M. G. Aartsen et al. [IceCube Collaboration], First observation of PeV-energy neutrinos with IceCube , Phys. Rev. Lett. 111 (2013) 021103 [arXiv:1304.5356 [astro-ph.HE]].
- [7] M. G. Aartsen et al. [IceCube Collaboration], Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector , Science 342 (2013) 1242856 [arXiv:1311.5238 [astro-ph.HE]].
- [8] M. G. Aartsen et al. [IceCube Collaboration], Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data , Phys. Rev. Lett. 113 (2014) 101101 [arXiv:1405.5303 [astro-ph.HE]].

- [9] M. G. Aartsen et al. [IceCube Collaboration], Search for a diffuse flux of astrophysical muon neutrinos with the IceCube 59-string configuration , Phys. Rev. D 89 (2014) 6, 062007 [arXiv:1311.7048 [astro-ph.HE]].
- [10] M. G. Aartsen et al. [IceCube Collaboration], Evidence for Astrophysical Muon Neutrinos from the Northern Sky with IceCube , Phys. Rev. Lett. 115 (2015) 8, 081102 [arXiv:1507.04005 [astro-ph.HE]].
- [11] B. Feldstein, A. Kusenko, S. Matsumoto and T. T. Yanagida, Neutrinos at IceCube from Heavy Decaying Dark Matter , Phys. Rev. D 88 (2013) no.1, 015004 [arXiv:1303.7320 [hep-ph]];
- [12] A. Esmaili and P. D. Serpico, Are IceCube neutrinos unveiling PeV-scale decaying dark matter? , JCAP 1311 (2013) 054 [arXiv:1308.1105 [hep-ph]].
- [13] K. Murase, R. Laha, S. Ando and M. Ahlers, Testing the Dark Matter Scenario for PeV Neutrinos Observed in IceCube , Phys. Rev. Lett. 115 (2015) no.7, 071301 [arXiv:1503.04663 [hep-ph]].
- [14] M. G. Aartsen et al. [IceCube Collaboration], Search for neutrinos from decaying dark matter with IceCube , Eur. Phys. J. C 78 (2018) no.10, 831 [arXiv:1804.03848 [astro-ph.HE]].
- [15] A. Bhattacharya, A. Esmaili, S. Palomares-Ruiz and I. Sarcevic, Update on decaying and annihilating heavy dark matter with the 6-year IceCube HESE data , arXiv:1903.12623 [hep-ph].
- [16] C. Y. Chen, P. S. Bhupal Dev and A. Soni, Two-component flux explanation for the high energy neutrino events at IceCube , Phys. Rev. D 92 (2015) no.7, 073001 [arXiv:1411.5658 [hep-ph]].
- [17] M. Chianese, G. Miele, S. Morisi and E. Vitagliano, Low energy IceCube data and a possible Dark Matter related excess , Phys. Lett. B 757 (2016) 251 [arXiv:1601.02934 [hep-ph]].
- [18] M. Chianese, G. Miele and S. Morisi, Dark Matter interpretation of low energy IceCube MESE excess , JCAP 1701 (2017) no.01, 007 [arXiv:1610.04612 [hep-ph]].
- [19] M. Chianese, G. Miele, S. Morisi and E. Peinado, Neutrinophilic Dark Matter in the epoch of IceCube and Fermi-LAT , JCAP 1812 (2018) no.12, 016 [arXiv:1808.02486 [hep-ph]].

- [20] Y. Sui and P. S. Bhupal Dev, A Combined Astrophysical and Dark Matter Interpretation of the IceCube HESE and Throughgoing Muon Events , JCAP 1807 (2018) no.07, 020 [arXiv:1804.04919 [hep-ph]].
- [21] S. M. Boucenna, M. Chianese, G. Mangano, G. Miele, S. Morisi, O. Pisanti and E. Vitagliano, Decaying Leptophilic Dark Matter at IceCube , JCAP 1512 (2015) no.12, 055 [arXiv:1507.01000 [hep-ph]].
- [22] M. Re Fiorentin, V. Niro and N. Fornengo, A consistent model for leptogenesis, dark matter and the IceCube signal , JHEP 1611 (2016) 022 [arXiv:1606.04445 [hep-ph]].
- [23] P. S. B. Dev, D. Kazanas, R. N. Mohapatra, V. L. Teplitz and Y. Zhang, Heavy right-handed neutrino dark matter and PeV neutrinos at IceCube , JCAP 1608 (2016) no.08, 034 [arXiv:1606.04517 [hep-ph]].
- [24] P. Di Bari and R. Foot, Active sterile neutrino oscillations in the early universe: Asymmetry generation at low | δm 2 | and the Landau-Zener approximation , Phys. Rev. D 65 (2002) 045003 [hep-ph/0103192].
- [25] S. Dodelson and L. M. Widrow, Sterile-neutrinos as dark matter , Phys. Rev. Lett. 72 (1994) 17 [hep-ph/9303287].
- [26] P. Minkowski, Phys. Lett. B 67 (1977) 421; T. Yanagida, in Proceedings of the Workshop on Unified Theory and Baryon Number of the universe, eds. O. Sawada and A. Sugamoto (KEK, 1979) p.95; P. Ramond, Invited talk given at Conference: C79-02-25 (Feb 1979) p.265-280, CALT-68-709, hep-ph/9809459; M. Gell-Mann, P. Ramond and R. Slansky, in Supergravity, eds. P. van Niewwenhuizen and D. Freedman (North Holland, Amsterdam, 1979) Conf.Proc. C790927 p.315, PRINT-80-0576; R. Barbieri, D. V. Nanopoulos, G. Morchio and F. Strocchi, Phys. Lett. B 90 (1980) 91; R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44 (1980) 912.
- [27] H. A. Weldon, Phys. Rev. D 26 (1982) 1394. H. A. Weldon, Phys. Rev. D 26 (1982) 2789. H. A. Weldon, Phys. Rev. D 28 (1983) 2007.
- [28] Y. Akrami et al. [Planck Collaboration], Planck 2018 results. I. Overview and the cosmological legacy of Planck , arXiv:1807.06205 [astro-ph.CO].
- [29] See for example P. Di Bari, Cosmology and the early universe , CRC Press, May 2018 (ISBN 9781498761703).

- [30] P. Di Bari, NATO Sci. Ser. C 503 (1997) 289 [hep-ph/9708246].
- [31] B. Garbrecht, F. Glowna and P. Schwaller, Scattering Rates For Leptogenesis: Damping of Lepton Flavour Coherence and Production of Singlet Neutrinos , Nucl. Phys. B 877 (2013) 1 [arXiv:1303.5498 [hep-ph]]. See also: D. Besak and D. Bodeker, Thermal production of ultrarelativistic right-handed neutrinos: Complete leading-order results , JCAP 03 (2012), 029 [arXiv:1202.1288 [hep-ph]]; M. Laine, Thermal right-handed neutrino production rate in the relativistic regime , JHEP 08 (2013), 138 [arXiv:1307.4909 [hep-ph]].
- [32] A. Atre, T. Han, S. Pascoli and B. Zhang, The Search for Heavy Majorana Neutrinos , JHEP 0905 (2009) 030 [arXiv:0901.3589 [hep-ph]].
- [33] M. Fukugita and T. Yanagida, Baryogenesis Without Grand Unification , Phys. Lett. B 174 (1986) 45.
- [34] S. Davidson and A. Ibarra, A Lower bound on the right-handed neutrino mass from leptogenesis , Phys. Lett. B 535 (2002) 25 [hep-ph/0202239]; W. Buchmuller, P. Di Bari and M. Plumacher, Cosmic microwave background, matter - antimatter asymmetry and neutrino masses , Nucl. Phys. B 643 (2002) 367 Erratum: [Nucl. Phys. B 793 (2008) 362] [hep-ph/0205349]. P. H. Chankowski and K. Turzynski, Phys. Lett. B 570 (2003) 198 doi:10.1016/j.physletb.2003.08.004 [hep-ph/0306059]; A. Abada, S. Davidson, A. Ibarra, F.-X. Josse-Michaux, M. Losada and A. Riotto, Flavour Matters in Leptogenesis , JHEP 0609 (2006) 010 [hep-ph/0605281]. S. Blanchet and P. Di Bari, Flavor effects on leptogenesis predictions , JCAP 0703 (2007) 018 [hep-ph/0607330]; S. Blanchet and P. Di Bari, New aspects of leptogenesis bounds , Nucl. Phys. B 807 (2009) 155 [arXiv:0807.0743 [hep-ph]]; S. Antusch, P. Di Bari, D. A. Jones and S. F. King, Leptogenesis in the Two Right-Handed Neutrino Model Revisited , Phys. Rev. D 86 (2012) 023516 [arXiv:1107.6002 [hep-ph]].
- [35] L. Covi, E. Roulet and F. Vissani, CP violating decays in leptogenesis scenarios , Phys. Lett. B 384 (1996) 169 [hep-ph/9605319].
- [36] M. D'Onofrio, K. Rummukainen and A. Tranberg, Sphaleron Rate in the Minimal Standard Model , Phys. Rev. Lett. 113 (2014) no.14, 141602 [arXiv:1404.3565 [hep-ph]].
- [37] E. K. Akhmedov, V. A. Rubakov and A. Y. Smirnov, Baryogenesis via neutrino oscillations , Phys. Rev. Lett. 81 (1998) 1359 [hep-ph/9803255].

- [38] For a pedagogical introduction to density matrix see for example R. Feynman, Statistical Mechanics: A set of Lectures , W. A. Benjamin, Inc., Reading, MA., 1972.
- [39] A. D. Dolgov, Sov. J. Nucl. Phys. 33 (1981) 700 [Yad. Fiz. 33 (1981) 1309]; R. A. Harris and L. Stodolsky, Two State Systems in Media and 'Turing's Paradox' , Phys. Lett. 116B (1982) 464; L. Stodolsky, On the Treatment of Neutrino Oscillations in a Thermal Environment , Phys. Rev. D 36 (1987) 2273. R. Barbieri and A. Dolgov, Neutrino oscillations in the early universe , Nucl. Phys. B 349 (1991) 743. G. Raffelt, G. Sigl and L. Stodolsky, Non Abelian Boltzmann equation for mixing and decoherence , Phys. Rev. Lett. 70 (1993) 2363 Erratum: [Phys. Rev. Lett. 98 (2007) 069902] [hep-ph/9209276].
- [40] P. Di Bari and R. Foot, Active sterile neutrino oscillations in the early universe: Asymmetry generation at low | δm 2 | and the Landau-Zener approximation , Phys. Rev. D 65 (2002) 045003 [hep-ph/0103192].
- [41] T. Asaka and M. Shaposhnikov, The nuMSM, dark matter and baryon asymmetry of the universe , Phys. Lett. B 620 (2005) 17 [hep-ph/0505013].
- [42] R. Barbieri, P. Creminelli, A. Strumia and N. Tetradis, Baryogenesis through leptogenesis , Nucl. Phys. B 575 (2000) 61 [hep-ph/9911315]; A. Abada, S. Davidson, F. X. Josse-Michaux, M. Losada and A. Riotto, Flavor issues in leptogenesis , JCAP 0604 (2006) 004 [hep-ph/0601083]; A. De Simone and A. Riotto, JCAP 0702 (2007) 005 doi:10.1088/1475-7516/2007/02/005 [hep-ph/0611357]; S. Blanchet, P. Di Bari, D. A. Jones and L. Marzola, Leptogenesis with heavy neutrino flavours: from density matrix to Boltzmann equations , JCAP 1301 (2013) 041 [arXiv:1112.4528 [hep-ph]]; K. Moffat, S. Pascoli, S. T. Petcov, H. Schulz and J. Turner, Three-flavored nonresonant leptogenesis at intermediate scales , Phys. Rev. D 98 (2018) no.1, 015036 [arXiv:1804.05066 [hep-ph]].
- [43] We adopt the same notation of N. F. Bell, R. R. Volkas and Y. Y. Y. Wong, Relic neutrino asymmetry evolution from first principles , Phys. Rev. D 59 (1999) 113001 [hep-ph/9809363].
- [44] P. Di Bari, D. Marfatia and Y. L. Zhou, Gravitational waves from neutrino mass and dark matter genesis , [arXiv:2001.07637 [hep-ph]].
- [45] P. Di Bari, R. Samanta and Y. L. Zhou, in preparation.

- [46] T. K. Kuo and J. T. Pantaleone, Neutrino Oscillations in Matter , Rev. Mod. Phys. 61 (1989), 937.
- [47] W. Buchmuller, P. Di Bari and M. Plumacher, Leptogenesis for pedestrians , Annals Phys. 315 (2005) 305 [hep-ph/0401240].
- [48] S. Blanchet and P. Di Bari, Leptogenesis beyond the limit of hierarchical heavy neutrino masses , JCAP 0606 (2006) 023 [hep-ph/0603107].
- [49] S. Blanchet and P. Di Bari, Flavor effects on leptogenesis predictions , JCAP 0703 (2007) 018 [hep-ph/0607330].
- [50] A. Anisimov, S. Blanchet and P. Di Bari, Viability of Dirac phase leptogenesis , JCAP 0804 (2008) 033 [arXiv:0707.3024 [hep-ph]].
- [51] M. Y. Khlopov and A. D. Linde, Phys. Lett. B 138 (1984) 265; J. R. Ellis, J. E. Kim and D. V. Nanopoulos, Phys. Lett. B 145 (1984) 181; K. Kohri, T. Moroi and A. Yotsuyanagi, Phys. Rev. D 73 (2006) 123511; M. Kawasaki, K. Kohri, T. Moroi and A. Yotsuyanagi, Phys. Rev. D 78 (2008) 065011 [arXiv:0804.3745 [hep-ph]].
- [52] M. Plumacher, Baryogenesis and lepton number violation , Z. Phys. C 74 (1997) 549 [hep-ph/9604229].
- [53] G. Sigl and G. Raffelt, General kinetic description of relativistic mixed neutrinos , Nucl. Phys. B 406 (1993) 423.
- [54] J. A. Casas and A. Ibarra, Oscillating neutrinos and µ → e γ , Nucl. Phys. B 618 (2001) 171 [hep-ph/0103065].
- [55] S. Antusch, P. Di Bari, D. A. Jones and S. F. King, Leptogenesis in the Two Right-Handed Neutrino Model Revisited , Phys. Rev. D 86 (2012) 023516 [arXiv:1107.6002 [hep-ph]].
- [56] I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, Global analysis of three-flavour neutrino oscillations: synergies and tensions in the determination of θ 2 3 , δ C P , and the mass ordering , JHEP 1901 (2019) 106 [arXiv:1811.05487 [hep-ph]]. NuFIT 4.1 (2019), www.nu-fit.org.
- [57] E. E. Jenkins and A. V. Manohar, Tribimaximal Mixing, Leptogenesis, and theta(13) , Phys. Lett. B 668 (2008) 210 [arXiv:0807.4176 [hep-ph]].

- [58] E. Bertuzzo, P. Di Bari, F. Feruglio and E. Nardi, Flavor symmetries, leptogenesis and the absolute neutrino mass scale , JHEP 0911 (2009) 036 [arXiv:0908.0161 [hep-ph]].
P. Di Bari K. Farrag R. Samanta Y. L. Zhou
Density matrix calculation of the dark matter abundance in the Higgs induced right-handed neutrino mixing model
We present new results on the calculation of the dark matter relic abundance
within the Higgs induced right-handed (RH) neutrino mixing model, solving the
associated density matrix equation. For a benchmark value of the dark matter
mass $M_{\rm DM} = 220\,{\rm TeV}$, we show the evolution of the abundance and
how this depends on reheat temperature, dark matter lifetime and source RH
neutrino mass $M_{\rm S}$, with the assumption $M_{\rm S} < M_{\rm DM}$. We
compare the results with those obtained within the Landau-Zener approximation,
showing that the latter largely overestimates the final abundance giving some
analytical insight. However, we also notice that since in the density matrix
formalism the production is non-resonant, this allows source RH neutrino masses
below the W boson mass, making dark matter more stable at large mass values.
This opens an allowed region for initial vanishing source RH neutrino
abundance. For example, for $M_{\rm S} \gtrsim 1\,{\rm GeV}$, we find $M_{\rm
DM}\gtrsim 20\,{\rm PeV}$. Otherwise, for $M_{\rm S} > M_W\sim 100\,{\rm GeV}$,
one has to assume a thermalisation of the source RH neutrinos prior to the
freeze-in of the dark matter abundance. This results into a large allowed range
for $M_{\rm DM}$, depending on $M_{\rm S}$. For example, imposing $M_{\rm S}
\gtrsim 300\,{\rm GeV}$, allowing also successful leptogenesis, we find $0.5
\lesssim M_{\rm DM}/{\rm TeV} \lesssim 500$. We also discuss in detail
leptogenesis with two quasi-degenerate RH neutrinos, showing a case when
observed dark matter abundance and matter-antimatter asymmetry are
simultaneously reproduced. Finally, we comment on how an initial thermal source
RH neutrino abundance can be justified and on how our results suggest that also
the interesting case where $M_{\rm DM} < M_{\rm S}$, embeddable in usual high
scale two RH neutrino seesaw models, might be viable.