Sterile neutrino dark matter from right-handed neutrino oscillations

We study a scenario where sterile neutrino (either warm or cold) dark matter (DM) is produced through (nonresonant) oscillations among right-handed neutrinos (RHNs) and can constitute the whole DM in the Universe, in contrast to the conventional sterile neutrino production through its mixing with the left-handed neutrinos. The lightest RHN can be sterile neutrino DM whose mixing with left-handed neutrinos is sufficiently small while heavier RHNs can have non-negligible mixings with left-handed neutrinos to explain the neutrino masses by the seesaw mechanism. We also demonstrate that, in our scenario, the production of sterile RHN DM from the decay of a heavier RHN is subdominant compared with the RHN oscillation production due to the X-ray and small-scale structure constraints.


I. INTRODUCTION
While it has been established that neutrinos are massive due to the discovery of neutrino oscillations [1,2], their precise properties, however, are still under active investigation. An analogous (and even more perplexing) story applies to dark matter (DM) whose nature remains unknown despite the ever-growing evidence for its existence from the astrophysical observables. An intriguing possibility regarding these mysteries would be to introduce right-handed neutrinos (RHNs), which can address not only the neutrino mass and DM but also their potential roles in the inflation and baryon asymmetry production [3][4][5][6][7][8][9].
We, in this article, seek a possibility for a sterile RHN to make up the whole DM in the Universe and, in particular, propose the new production mechanism of sterile RHN DM through the mixing among RHNs. This is in contrast to the conventional mechanisms requiring the sterile RHN DM to couple to left-handed neutrinos which suffer from the severe tension between the bounds from the X-ray observation and the small-scale structure data [10][11][12][13][14][15]. These constraints, however, heavily depend on their production mechanisms and many possibilities have been explored to produce the desired DM abundance in addition to the conventional nonresonant/resonant active-sterile neutrino conversion mechanisms [6][7][8][9][16][17][18][19][20].
Our scenario is distinguishable from such alternative scenarios in that it still uses a simple oscillation between the thermal heavy RHN and DM, and yet it demonstrates the totally different features from the Dodelson-Widrow scenario such as the occurrence of the production peak above/around the electroweak which is of great advantage in circumventing the Lyman-α bounds due to the redshifting of DM momentum. After outlining our setup in Sec. II, we illustrate our scenario in Sec. III for a simple example of two RHNs. Section IV then demonstrates the concrete realization where we introduce a RHN mass matrix whose off-diagonal term can arise from the scalar field vacuum expectation value so that we can explain the light neutrino masses by the seesaw mechanism while avoiding the tight X-ray bounds. Section V is devoted to the discussion/conclusion.

II. SETUP
The Lagrangian we study is the standard model (SM) with three Majorana RHNs, given by L = L SM + L N where L SM is the SM Lagrangian and L N reads where H, L, and ν R are, respectively, the Higgs doublet, lepton doublet and RHN. For simplicity, we concentrate on the case of three RHNs. We begin with the field basis where y ν y † ν is diagonal, denoted as y diag ν so that y diag ν y diag † ν becomes a 3 × 3 diagonal matrix. M N is, in general, a nondiagonal matrix, which we call the interaction basis. A familiar seesaw mechanism for the mass of left-handed neutrino ν L reads, in terms of its Dirac mass m diag By taking the rotation of Eq. (2), the Yukawa coupling y ν is in general a nondiagonal matrix while the neutrino masses, M ν and M N , are simultaneously diagonalized. We call this field basis the mass basis. Thus, we obtain where R is an arbitrary 3 × 3 complex orthogonal matrix satisfying R T R = 1 [21]. The mixing between ν L and N is then parametrized by Θ = θ † U * R , and The oscillations among RHNs can take place when their mass and interaction bases differ. We, in the following discussions, consider three RHNs with their masses M diag The lightest neutrino mass (m 1 for the NH case, and m 3 for the IH case) is taken as a free parameter. In our discussions below, whenever it is not necessary to distinguish the mass orderings, m 1 refers to the lightest mass for brevity.

III. DM PRODUCTION THROUGH RHN OSCILLATION
We now check if enough abundance of RHN DM ν R1 can be produced from the RHN oscillations. In our scenario, the heavy RHNs ν R2 and ν R3 explain the lefthanded neutrino masses by the seesaw mechanism and they can have sizable neutrino Yukawa couplings to be in the thermal equilibrium at a sufficiently high temperature. ν R1 , on the other hand, has a sufficiently small coupling to the SM species, so that its production is dominated by the conversion from heavier RHNs. For clarity of the following quantitative discussion, we focus on the ν R1 abundance produced only from its mixing with ν R2 because ν R3 plays the same role as ν R2 in producing ν R1 .
The relevant reactions for the ν R2 thermalization are the scatterings caused by Yukawa interaction, The Boltzmann equation for ν R1 [23] reads where C ν R1 represents the collision term integrated over the ν R1 momentum given by The ratios between the rescaled (i.e., divided by the Yukawa couplings) reaction rates and the Hubble parameter are shown (the actual reaction rates are obtained by multiplying the Yukawa couplings). The solid curves are for M2 = 1 GeV and the dashed curves are for M2 = 1 TeV.
Here P is the oscillation probability given by is the decay width, andσ is the reduced cross section for the ν R2 collisions with the kinematical cut s min of the Mandelstam variable s, and K 1 is the modified Bessel function of the first kind. 2 ν R1 is efficiently produced when the collision terms are large. 3 Figure 1 shows Γ i /H where Γ i represents the rescaled reaction rates for the process i by taking the neutrino Yukawa coupling as unity (so that the curves can be easily scaled by multiplying the Yukawa coupling 2 A factor 1/2 in P comes from averaging out the RHN oscillation because the oscillation timescale is much shorter than the collision timescale involving ν R2 . More quantitatively, this averaging is justified for T 10 6 GeV and/or where g represents a gauge coupling for a relevant gauge interaction. As we will discuss later, y 2 ν of order 10 −14 is required for GeV-scale RHN to reach the thermal equilibrium and it is automatically realized by enforcing the seesaw mechanism. The finite temperature effects on the RHN mixing angle θ N are suppressed by the neutrino Yukawa couplings in our scenario and we simply consider a constant θ N in our estimation. The cases when these approximations are not applicable are left for the future work. 3 Some of collision terms, such as ν R H → LV , possess the infrared divergences, which are regulated by the thermal mass of the propagator in our analysis for T > T C (T C is the critical temperature of the electroweak phase transition and we take T C = 160 GeV) [24][25][26][27].
of interest). For illustration purpose, we define the reaction rates Γ i = γ col ν R2 (i)/n γ , where n γ = 2T 3 /π 2 is the radiation number density and γ col ν R2 (i) are the collision terms involving the gauge bosons [γ col ν R2 (gauge)] and the top quarks [γ col ν R2 (top)]. The inverse decay rate is given by Γ ID = γ ID ν R2 /n γ . The figure shows the plots for M 2 = 1 GeV (solid) and for M 2 = 1 TeV (dashed), and we note that the inverse decay takes place only for the latter because of the kinematics, namely, the (inverse) decay is available only for M 2 M h with M h being the Higgs mass. The actual reaction rates can be obtained by multiplying these rescaled reaction rates by (y ν y † ν ) 22 . We can see, from Fig. 1, that N 2 can reach the thermal equilibrium (Γ i /H 1) when (y ν y † ν ) 22 is larger than O(10 −13 ) for M 2 = 1−10 3 GeV, which is also in the desired numerical range to explain the neutrino masses by the seesaw mechanism.
The produced ν R1 (interaction state) constitutes the DM N 1 (mass eigenstate), 4 and the current N 1 relic number density can be estimated, in terms of the yield parameter Y N1 ≡ n N1 /s (s is the entropy density), by integrating the Boltzmann equation from T RH , the reheating temperature, to the current temperature where we have taken the limits T RH → ∞, T 0 → 0, and whereỸ 0 N1 is the rescaled yield parameter, defined by factoring out the oscillation probability and the Yukawa coupling,Ỹ 0 22 ). We found the following simple fitting formula to grasp the characteristic features of the DM abundance in our scenario This behavior matches our expectation because, as emphasized in referring to Fig. 1, the most efficient production occurs when the production rate reaches maximal with respect to the Hubble expansion rate.Ỹ 0 N1 is hence 4 The produced ν R1 is composed of N 1 and N 2 which propagate with different velocities. As the ν R1 energy gets redshifted, these two mass states are eventually well separated and thus ν R1 is expected to mostly develop the N 1 component as long as M 1 M 2 , although the oscillation property may call for a careful study [28]. little dependent on M 2 when M 2 is smaller than M h , because N 2 is dominantly produced via the inverse decay in this case, and thus the temperature at which the production rate becomes maximal is at T M h . For M 2 M h , on the other hand, the SM particles possess the thermal mass and the production rate becomes maximal around T ∼ M 2 , which leads to some power dependence of the yield parameter on M 2 . This is illustrated through a concrete example in the next section.

IV. BENCHMARK MODEL
We here discuss a possible realization of our scenario. Let us begin with a simple mass matrix given by where m and M 0 are taken to be M 0 m The resultant N 1 abundance in the NH case is then given by Ω N1 h 2 0.12 m 2 0.01 eV while, in the IH case, m 2 should be replaced by m 1 . In the case of M 0 m 2 /M 2 , we can take θ N (M 1 /M 2 ) 1/2 due to M 1 m 2 /M 2 , and thus we obtain For this simplified case, Fig. 2 shows Ω N1 h 2 as a function of M 2 for various M 1 taken from 100 keV to 100 MeV in both the NH and IH cases which are depicted by solid and dashed curves, respectively. 5 The green band in the figure indicates the observed value of the DM abundance given by Ω DM h 2 = 0.1197 ± 0.0022 [29].
On the other hand, since Θ 2 11 depends on θ N and we need a relatively large θ N for our scenario to work, the N 1 is subject to the X-ray constraint given by Θ 2 11 10 −5 (keV/M 1 ) 5 [12]. One may simply expect that the X-ray bound is easily circumvented because the Yukawa coupling of ν R1 can be negligibly small. We, however, point out that the light RHN can decay into the SM particles through its oscillation to a heavier RHN. In our current setup, we obtain Θ 2 11 = M −1 1 (m 1 |R 11 | 2 + m 2 |R 12 | 2 + m 3 |R 13 | 2 ), where R ij represents the (i, j) entry of the R matrix. We can now take R 13 = 0, since there is no mixing in this component, and m 1 = 0 is experimentally allowed. However, since we have in the NH case, and m 2 is replaced by m 1 |∆m 2 32 | in the IH case. One may naively expect that this decay of light RHN through a heavier RHN is suppressed by the hierarchically large mass ratio M 1 /M 2 1. If we did not enforce the simple seesaw mechanism to obtain the desirable light neutrino masses, this would be the case and the X-ray bound could be circumvented. We, however, in our model construction stick to the seesaw mechanism to account for the observed neutrino masses, which then inevitably increase y 2 if we choose a bigger value of M 2 to result in too big an X-ray decay rate. To keep the virtue of explaining the observed neutrino masses by the simple type-I seesaw mechanism and yet not to lose the attractive feature of simple RHN oscillation production, we now discuss a time-dependent RHN mixing to evade the X-ray constraint mentioned above.
Such a time-dependent RHN mixing can be achieved by utilizing the dynamics of a real scalar filed φ. Let us here consider the two flavor case for simplicity, but the extension to the three flavor system is straightforward. In the two flavor case, we impose Z 2 symmetry under 5 It should be noted that, in Fig. 2, tosc/t col 1 is achieved for T 10 6 × M 2 even in the large M 2 region, so that a factor 1/2 in P by averaging out the RHN oscillation is justified. which ν R2 is even, while ν R1 and φ are odd. 6 Now the mass matrix M N in Eq. (1) is given by in the interaction basis. The dynamics of φ is governed by the equation of motionφ + 3Hφ + V (φ) = 0, where V (φ) is the potential that we take V (φ) 3H, φ is almost constant, namely, φ 2ρ φ /m φ with ρ φ the energy density of φ, and when H drops below m φ , φ starts to oscillate. As we will see below, m φ 3H is always satisfied when the N 1 production rate is maximal, and thus we take φ as a constant in this regime.
The mixing angle between ν R1 and ν R2 is given by sin θ N κφ/M 2 in the case that M 1 M 2 , and thus in the constant φ regime we obtain sin 2 2θ N 4κ 2 ρ φ /(m 2 φ M 2 2 ), where the relevant θ N is determined by ρ φ (T max ) with T max being the temperature at which the production rate becomes maximal, namely, T max ∼ T c for M 2 T c and otherwise T max ∼ M 2 . As mentioned above, φ is constant until it starts to oscillate, so we can take ρ φ (T max ) ρ φ (T osc ) with T osc given by m φ = 3H(T osc ). Then, we obtain with r g = g * (T osc )/g * (T 0 ), and r = ρ 0 φ /ρ DM with ρ 0 φ and ρ DM being the energy density of φ and dark matter at the present. Here we have used g * (T 0 ) 3.36.
We also require that φ never thermalizes by taking a sufficiently small κ not to affect the big bang nucleosynthesis, which results in κ 2 M 2 /M Pl . In addition, m φ should be smaller than H(T max ) in order for φ at T max to be constant, where H(T max ) 10 −5 eV for M 2 < T c and H(T max ) 10 −5 × (M 2 /T c ) 2 for M 2 > T c .
It is worth mentioning that the dynamics of φ may be tied to inflationary models. In particular, the condition of ρ 0 φ ρ DM implies that the initial amplitude of φ is bounded φ 4 × 10 11 GeV r g 30 On the other hand, φ could be largely displaced from the origin during inflation and its oscillation at a later time possibly dominates the dark matter energy density, in an analogous manner to the Polonyi/moduli problem [31][32][33]. To suppress φ in our case, we may utilize a relatively strong coupling between φ and inflaton, which renders the adiabatic suppression of the amplitude of the coherent oscillations [34]. Its actual dynamics, however, depends on the inflationary models and how φ couples to the inflaton, which we leave unspecified for the future work.
Finally let us comment on the θ N at the present, which is relevant for the decay of N 1 . Below T osc , since ρ φ drops as a matter energy density, we obtain and therefore a sufficiently small mixing to avoid the Xray constraint can be achieved.

V. DISCUSSION/CONCLUSION
Before concluding our discussions, let us briefly point out another potentially interesting production mechanism: the production of N 1 from a heavier RHN decay. We can consider the decay of N 2 (and/or N 3 ) which is thermally decoupled while it is relativistic (otherwise N 2 number density would be too small due to the Boltzmann suppression). N 1 abundance then can be estimated as Ω N1 h 2 10 −10 Θ 2 11 10 −12 where we used the branching fraction of N 2 decay for the process N 2 → N 1 +(mesons, leptons), Br( , and the ratio of g * accounts for the change in the effective degrees of freedom from the N 2 freeze-out epoch to the present time. This production contribution is hence subdominant compared with RHN oscillation production in the parameter region of our interest. Let us next mention the small-scale structure constraints applicable to our scenario. We here discuss the Lyman-α forest constraints which can give the lower limit on the DM mass from the DM free streaming scale λ F S ∼ 1 Mpc(keV/M 1 )( p/T /3.15) [35]. Too large a free streaming scale can be excluded due to the suppression of small-scale structure formation. The average momentum of N 1 produced by the nonresonant oscillation of thermalized N 2 can be estimated as p 1 ∼ 2.8T , analogous to the conventional (nonresonant) active-sterile oscillation scenario. Taking account of momentum redshifting by a factor (g * (T N2→N1 )/g * (T MeV)) −1/3 due to the change in the effective degrees of freedom, Lyman-α data leads to the RHN DM mass bound M 1 10 keV for our scenario [14] (when N 2 → N 1 occurs most efficiently before the QCD phase transition which is the case for the parameter range discussed so far). Such a DM mass range can be realized in our scenario as explicitly demonstrated through the concrete examples in the last section while being compatible with both the right relic abundance and seesaw mechanism.
Among the possible extensions of our DM scenarios, we plan to study the leptogenesis as well as the neutrino observables such as the neutrinoless double beta decay in our future work. For instance, even though we have focused on the DM production in this article, the neutrino Yukawa couplings in our model can be further constrained by seeking the production of desirable baryon asymmetry in the Universe. The realization of leptogenesis when N 2 and N 3 are heavy enough and/or are degenerate in their masses with sufficient CP violations [5,6,24] will be explored in our forthcoming paper. The CP phases in the neutrino Yukawa couplings are of great importance not only for the leptogenesis but also for the DM production in our scenario, and the presented production mechanism for the RHN DM could uncover a new connection between DM and leptogenesis to bring considerable opportunities for subsequent studies.

ACKNOWLEDGMENTS
This work was supported by IBS under the project code IBS-R018-D1. We thank A. Kamada, A. Merle and T. Asaka for useful discussions and, in particular, the anonymous referee for the constructive suggestions.