Sterile neutrino dark matter from right-handed neutrino oscillations

TL;DR

This paper proposes a novel mechanism for sterile neutrino dark matter production via right-handed neutrino oscillations, offering an alternative to traditional mixing-based models, and discusses its implications for cosmology and particle physics.

Contribution

It introduces a new scenario where sterile neutrino dark matter is generated through right-handed neutrino oscillations, differing from conventional mixing-based production methods.

Findings

Oscillation-based production can account for all dark matter.

Decay of heavier RHNs is subdominant in sterile neutrino production.

Model complies with X-ray and small-scale structure constraints.

Abstract

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.

Equations

$${\overline{\nu }}_{R}i/\partial {\nu }_R-\left[{\nu }_R^c{}^T{y}_{\nu }LH-\frac{1}{2}{\nu }_R^c{}^T{\mathcal{M}}_N{\nu }_R^c+H.c.\right],$$

$$\displaystyle\left[\begin{array}[]{c}\nu_{L}\\ \nu_{R}^{c}\end{array}\right]$$

$$y_{\nu }^{\mathrm{diag}}{y}_{\nu }^{\mathrm{diag}}{}^{†}$$

$${\Theta }^2$$

$$\frac{d{n}_{{\nu }_{R1}}}{dt}+3H{n}_{{\nu }_{R1}}=C_{{\nu }_{R1}}$$

$$C_{{\nu }_{R1}}$$

$${\gamma }_{{\nu }_{R2}}^{\mathrm{col}}$$

$${\gamma }_{{\nu }_{R2}}^{\mathrm{ID}}$$

$$Y_{N_1}^0\equiv Y_{N_1}(T=0)={\int }_0^{\infty }dT\mathcal{P}({\nu }_{R2}\to {\nu }_{R1})\frac{{\gamma }_{{\nu }_{R2}}}{sHT},$$

$${\Omega }_{N_1}{h}^2$$

$${\log }_{10}{\tilde{Y}}_{N_1}^0$$

$${\mathcal{M}}_N$$

$$U_R^{\ast }$$

$${\Omega }_{N_1}{h}^2$$

$$\displaystyle{\cal M}_{N}(\phi)=\left[\begin{array}[]{cc}M_{1}&\kappa\phi\\ \kappa\phi&M_{2}\end{array}\right]$$

$${\sin }^{2}2{\theta }_N$$

$$\varphi \lesssim 4\times 10^{11}\mathrm{GeV}{\left(\frac{r_g}{30}\right)}^{1/2}{\left(\frac{r}{10^{-4}}\right)}^{1/2}{\left(\frac{10^{-4}\mathrm{eV}}{m_{\varphi }}\right)}^{1/4}.$$

$$\frac{{\sin }^{2}2{\theta }_N(T_0)}{{\sin }^{2}2{\theta }_N(T_{\mathrm{osc}})}\simeq 1.2\times 10^{-46}{\left[\frac{r_g}{30}\right]}^{-1/4}{\left[\frac{m_{\varphi }}{10^{-4}\mathrm{eV}}\right]}^{-3/2},$$

$${\Omega }_{N_1}{h}^2$$

Full text

††preprint: CTPU-17-04

I Introduction

While it has been established that neutrinos are massive due to the discovery of neutrino oscillations Fukuda:1998mi ; Ahmad:2002jz , 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 seesaw ; Asaka:2005an ; Fukugita:1986hr ; Asaka:2005pn ; Canetti:2012kh ; Ibe:2015nfa ; kk2 .

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 Dodelson:1993je ; Shi:1998km ; Tremaine:1979we ; Boyarsky:2005us ; Horiuchi:2013noa ; ir2017 . 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 Asaka:2006ek ; Bezrukov:2009th ; kenji ; Anisimov:2008gg ; Adhikari:2016bei ; Asaka:2005pn ; Canetti:2012kh ; Ibe:2015nfa ; kk2 .

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-$\alpha$ 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 ${\cal L}={\cal L}_{\rm SM}+{\cal L}_{N}$ where ${\cal L}_{\rm SM}$ is the SM Lagrangian and ${\cal L}_{N}$ reads

[TABLE]

where $H,L$, and $\nu_{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_{\nu}y_{\nu}^{\dagger}$ is diagonal, denoted as $y_{\nu}^{\rm diag}$ so that $y_{\nu}^{\rm diag}y_{\nu}^{{\rm diag}\dagger}$ becomes a $3\times 3$ diagonal matrix. ${\cal 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 $\nu_{L}$ reads, in terms of its Dirac mass $m_{D}^{\rm diag}=y_{\nu}^{\rm diag}v$ with $v=\langle H\rangle$, ${\cal M}_{\nu}=m_{D}^{\rm diag}{}^{T}{\cal M}_{N}^{-1}m_{D}^{\rm diag}$ which can be diagonalized as ${\cal M}_{\nu}^{\rm diag}=U_{L}^{T}{\cal M}_{\nu}U_{L}$ ($U_{L}$ is the Pontecorvo-Maki-Nakagawa-Sakata matrix.111Throughout this article, we take the charged lepton Yukawa coupling to be diagonal.). The neutrino mass eigenstates are

[TABLE]

where $\theta\equiv{\cal M}_{N}^{-1}m_{D}^{\rm diag}$ and $U_{R}$ is a unitary matrix defined to diagonalize ${\cal M}_{N}$ as ${\cal M}_{N}^{\rm diag}=U_{R}^{\dagger}{\cal M}_{N}U_{R}^{*}$. By taking the rotation of Eq. (10), the Yukawa coupling $y_{\nu}$ is in general a nondiagonal matrix while the neutrino masses, ${\cal M}_{\nu}$ and ${\cal M}_{N}$, are simultaneously diagonalized. We call this field basis the mass basis. Thus, we obtain

[TABLE]

where $R$ is an arbitrary $3\times 3$ complex orthogonal matrix satisfying $R^{T}R=1$ Casas:2001sr . The mixing between $\nu_{L}$ and $N$ is then parametrized by $\Theta=\theta^{\dagger}U_{R}^{*}$, and

[TABLE]

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 ${\cal M}_{N}^{\rm diag}={\rm diag}\{M_{1},M_{2},M_{3}\}$ and take $N_{1}$ as the lightest one so that it can play the role of DM. For the active neutrino masses, we parametrize ${\cal M}_{\nu}^{\rm diag}={\rm diag}\{m_{1},m_{2},m_{3}\}$ for the normal hierarchy (NH), where $\Delta m_{21}^{2}\equiv m_{2}^{2}-m_{1}^{2}=(7.50^{+0.19}_{-0.17})\times 10^{-5}~{}{\rm eV}^{2},\Delta m_{31}^{2}\equiv m_{3}^{2}-m_{1}^{2}=(2.457^{+0.047}_{-0.047})\times 10^{-3}~{}{\rm eV}^{2}$ Gonzalez-Garcia:2015qrr . For the inverted hierarchy (IH), we take ${\cal M}_{\nu}^{\rm diag}={\rm diag}\{m_{3},m_{1},m_{2}\}$ and $\Delta m_{32}^{2}\equiv m_{3}^{2}-m_{2}^{2}=(-2.449^{+0.048}_{-0.047})\times 10^{-3}~{}{\rm eV}^{2}$. 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 $\nu_{R1}$ can be produced from the RHN oscillations. In our scenario, the heavy RHNs $\nu_{R2}$ and $\nu_{R3}$ explain the left-handed 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. $\nu_{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 $\nu_{R1}$ abundance produced only from its mixing with $\nu_{R2}$ because $\nu_{R3}$ plays the same role as $\nu_{R2}$ in producing $\nu_{R1}$.

The relevant reactions for the $\nu_{R2}$ thermalization are the scatterings caused by Yukawa interaction, $\nu_{R2}L\leftrightarrow tQ_{3},~{}\nu_{R2}t\leftrightarrow LQ_{3},~{}\nu_{R2}Q_{3}\leftrightarrow Lt$, those involving the gauge bosons, $\nu_{R2}V\leftrightarrow HL,~{}\nu_{R2}L\leftrightarrow HV,~{}\nu_{R2}H\leftrightarrow LV$ and the decay and inverse decay $\nu_{R2}\leftrightarrow LH$ [$Q_{3}(t)$ is the left (right) handed top quark, and $V$ represents the $SU(2)_{L}$ and $U(1)_{Y}$ gauge bosons].

The Boltzmann equation for $\nu_{R1}$ Dolgov:2000ew reads

[TABLE]

where $C_{\nu_{R1}}$ represents the collision term integrated over the $\nu_{R1}$ momentum given by

[TABLE]

Here ${\cal P}$ is the oscillation probability given by ${\cal P}(\nu_{R2}\to\nu_{R1})=\frac{1}{2}\sin^{2}2\theta_{N}$ ($\theta_{N}$ is the mixing angle between $\nu_{R1}$ and $\nu_{R2}$), $\Gamma(\nu_{R2}\to LH)\simeq(y_{\nu}y_{\nu}^{\dagger})_{22}M_{2}/(8\pi)$ is the decay width, and $\hat{\sigma}$ is the reduced cross section for the $\nu_{R2}$ collisions with the kinematical cut $s_{\rm min}$ of the Mandelstam variable $s$, and $K_{1}$ is the modified Bessel function of the first kind. 222A factor 1/2 in ${\cal P}$ comes from averaging out the RHN oscillation because the oscillation timescale is much shorter than the collision timescale involving $\nu_{R2}$. More quantitatively, this averaging is justified for $T\lesssim 10^{6}$ GeV and/or $\Delta M^{2}\equiv M_{2}^{2}-M_{1}^{2}\gtrsim 1~{}{\rm GeV}^{2}$ because $t_{\rm osc}/t_{\rm col}\sim({y_{\nu}^{2}}/{10^{-14}})({g^{2}}/{10^{-2}})({\rm GeV}^{2}/{\Delta M^{2}})({T}/{10^{6}~{}{\rm GeV}})^{2}$ where $g$ represents a gauge coupling for a relevant gauge interaction. As we will discuss later, $y_{\nu}^{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 $\theta_{N}$ are suppressed by the neutrino Yukawa couplings in our scenario and we simply consider a constant $\theta_{N}$ in our estimation. The cases when these approximations are not applicable are left for the future work.

$\nu_{R1}$ is efficiently produced when the collision terms are large.333Some of collision terms, such as $\nu_{R}H\to 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) Pilaftsis:2003gt ; Besak:2012qm ; DOnofrio:2014rug ; DOnofrio:2015gop . Figure 1 shows $\Gamma_{i}/H$ where $\Gamma_{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 of interest). For illustration purpose, we define the reaction rates $\Gamma_{i}=\gamma^{\rm col}_{\nu_{R2}}(i)/n_{\gamma}$, where $n_{\gamma}=2T^{3}/\pi^{2}$ is the radiation number density and $\gamma^{\rm col}_{\nu_{R2}}(i)$ are the collision terms involving the gauge bosons [$\gamma^{\rm col}_{\nu_{R2}}({\rm gauge})$] and the top quarks [$\gamma^{\rm col}_{\nu_{R2}}({\rm top})$]. The inverse decay rate is given by $\Gamma_{\rm ID}=\gamma^{\rm ID}_{\nu_{R2}}/n_{\gamma}$. 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}\gtrsim M_{h}$ with $M_{h}$ being the Higgs mass. The actual reaction rates can be obtained by multiplying these rescaled reaction rates by $(y_{\nu}y_{\nu}^{\dagger})_{22}$. We can see, from Fig. 1, that $N_{2}$ can reach the thermal equilibrium ($\Gamma_{i}/H\gtrsim 1$) when $(y_{\nu}y_{\nu}^{\dagger})_{22}$ is larger than ${\cal 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 $\nu_{R1}$ (interaction state) constitutes the DM $N_{1}$ (mass eigenstate), 444The produced $\nu_{R1}$ is composed of $N_{1}$ and $N_{2}$ which propagate with different velocities. As the $\nu_{R1}$ energy gets redshifted, these two mass states are eventually well separated and thus $\nu_{R1}$ is expected to mostly develop the $N_{1}$ component as long as $M_{1}\ll M_{2}$, although the oscillation property may call for a careful study Akhmedov:2012uu . and the current $N_{1}$ relic number density can be estimated, in terms of the yield parameter $Y_{N_{1}}\equiv n_{N_{1}}/s$ ($s$ is the entropy density), by integrating the Boltzmann equation from $T_{\rm RH}$, the reheating temperature, to the current temperature $T=T_{0}$

[TABLE]

where we have taken the limits $T_{\rm RH}\to\infty,T_{0}\to 0$, and $\gamma_{\nu_{R2}}\equiv\gamma_{\nu_{R2}}^{\rm col}+\gamma_{\nu_{R2}}^{\rm ID}$. The corresponding DM density can then be estimated in terms of the yield parameter

[TABLE]

where $\tilde{Y}^{0}_{N_{1}}$ is the rescaled yield parameter, defined by factoring out the oscillation probability and the Yukawa coupling, $\tilde{Y}^{0}_{N_{1}}\equiv Y^{0}_{N_{1}}/({\cal P}(\nu_{R2}\to\nu_{R1})(y_{\nu}y_{\nu}^{\dagger})_{22})$. We found the following simple fitting formula to grasp the characteristic features of the DM abundance in our scenario

[TABLE]

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. $\tilde{Y}^{0}_{N_{1}}$ is hence 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\simeq M_{h}$. For $M_{2}\gtrsim M_{h}$, on the other hand, the SM particles possess the thermal mass and the production rate becomes maximal around $T\sim 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

[TABLE]

where $m$ and $M_{0}$ are taken to be $M_{0}\lesssim m\ll M_{2},M_{3}$. ${\cal M}_{N}$ is then diagonalized as ${\cal M}_{N}^{\rm diag}={\rm diag}\{M_{1},M_{2},M_{3}\}$ with $M_{1}\simeq M_{0}-m^{2}/M_{2}$ by using $U_{R}$ which reads

[TABLE]

The resultant $N_{1}$ abundance in the NH case is then given by

[TABLE]

while, in the IH case, $m_{2}$ should be replaced by $m_{1}$. In the case of $M_{0}\ll m^{2}/M_{2}$, we can take $\theta_{N}\simeq(M_{1}/M_{2})^{1/2}$ due to $M_{1}\simeq m^{2}/M_{2}$, and thus we obtain

[TABLE]

For this simplified case, Fig. 2 shows $\Omega_{N_{1}}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. 555It should be noted that, in Fig. 2, $t_{\rm osc}/t_{\rm col}\ll 1$ is achieved for $T\lesssim 10^{6}\times M_{2}$ even in the large $M_{2}$ region, so that a factor 1/2 in ${\cal P}$ by averaging out the RHN oscillation is justified. The green band in the figure indicates the observed value of the DM abundance given by $\Omega_{\rm DM}h^{2}=0.1197\pm 0.0022$ Ade:2015xua .

On the other hand, since $\Theta^{2}_{11}$ depends on $\theta_{N}$ and we need a relatively large $\theta_{N}$ for our scenario to work, the $N_{1}$ is subject to the X-ray constraint given by $\Theta^{2}_{11}\lesssim 10^{-5}({\rm keV}/M_{1})^{5}$ Boyarsky:2005us . One may simply expect that the X-ray bound is easily circumvented because the Yukawa coupling of $\nu_{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 $\Theta^{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 $|R_{12}|^{2}=1/(1+(M_{1}/M_{2})\cot^{2}\theta_{N})\sim 1/2$ in our setup with $M_{1}/M_{2}\ll 1$, large $M_{1}$ is not allowed because of the X-ray constraint $\Theta^{2}_{11}\simeq m_{2}/(2M_{1})\lesssim 10^{-5}({\rm keV}/M_{1})^{5}$, where $m_{2}\simeq\sqrt{\Delta m^{2}_{21}}$ in the NH case, and $m_{2}$ is replaced by $m_{1}\simeq\sqrt{|\Delta 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}\ll 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 $\phi$. 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 which $\nu_{R2}$ is even, while $\nu_{R1}$ and $\phi$ are odd. 666 Although our setup is similar to the idea discussed in Ref. Berlin:2016bdv , the DM production scenario is quite different, since our scenario does not rely on the oscillation between active and sterile neutrinos, and thus the temperature at which the production efficiently occurs takes rather a wide range, which can imprint an observable signature on the structure formation.

Now the mass matrix ${\cal M}_{N}$ in Eq. (1) is given by

[TABLE]

in the interaction basis. The dynamics of $\phi$ is governed by the equation of motion $\ddot{\phi}+3H\dot{\phi}+V^{\prime}(\phi)=0$, where $V(\phi)$ is the potential that we take $V(\phi)\simeq(1/2)m_{\phi}^{2}\phi^{2}$. For $m_{\phi}\ll 3H$, $\phi$ is almost constant, namely, $\phi\simeq\sqrt{2\rho_{\phi}}/m_{\phi}$ with $\rho_{\phi}$ the energy density of $\phi$, and when $H$ drops below $m_{\phi}$, $\phi$ starts to oscillate. As we will see below, $m_{\phi}\ll 3H$ is always satisfied when the $N_{1}$ production rate is maximal, and thus we take $\phi$ as a constant in this regime.

The mixing angle between $\nu_{R1}$ and $\nu_{R2}$ is given by $\sin\theta_{N}\simeq\kappa\phi/M_{2}$ in the case that $M_{1}\ll M_{2}$, and thus in the constant $\phi$ regime we obtain $\sin^{2}2\theta_{N}\simeq 4\kappa^{2}\rho_{\phi}/(m_{\phi}^{2}M_{2}^{2})$, where the relevant $\theta_{N}$ is determined by $\rho_{\phi}(T_{\rm max})$ with $T_{\max}$ being the temperature at which the production rate becomes maximal, namely, $T_{\rm max}\sim T_{c}$ for $M_{2}\lesssim T_{c}$ and otherwise $T_{\rm max}\sim M_{2}$. As mentioned above, $\phi$ is constant until it starts to oscillate, so we can take $\rho_{\phi}(T_{\rm max})\simeq\rho_{\phi}(T_{\rm osc})$ with $T_{\rm osc}$ given by $m_{\phi}=3H(T_{\rm osc})$. Then, we obtain

[TABLE]

with $r_{g}=g_{*}(T_{\rm osc})/g_{*}(T_{0})$, and $r=\rho_{\phi}^{0}/\rho_{\rm DM}$ with $\rho_{\phi}^{0}$ and $\rho_{\rm DM}$ being the energy density of $\phi$ and dark matter at the present. Here we have used $g_{*}(T_{0})\simeq 3.36$.

We also require that $\phi$ never thermalizes by taking a sufficiently small $\kappa$ not to affect the big bang nucleosynthesis, which results in $\kappa^{2}\lesssim M_{2}/M_{\rm Pl}$. In addition, $m_{\phi}$ should be smaller than $H(T_{\rm max})$ in order for $\phi$ at $T_{\rm max}$ to be constant, where $H(T_{\rm max})\simeq 10^{-5}$ eV for $M_{2}<T_{c}$ and $H(T_{\rm max})\simeq 10^{-5}\times(M_{2}/T_{c})^{2}$ for $M_{2}>T_{c}$.

It is worth mentioning that the dynamics of $\phi$ may be tied to inflationary models. In particular, the condition of $\rho_{\phi}^{0}\ll\rho_{\rm DM}$ implies that the initial amplitude of $\phi$ is bounded

[TABLE]

On the other hand, $\phi$ 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 Coughlan:1983ci ; Ellis:1986zt ; Goncharov:1984qm . To suppress $\phi$ in our case, we may utilize a relatively strong coupling between $\phi$ and inflaton, which renders the adiabatic suppression of the amplitude of the coherent oscillations Linde:1996cx . Its actual dynamics, however, depends on the inflationary models and how $\phi$ couples to the inflaton, which we leave unspecified for the future work.

Finally let us comment on the $\theta_{N}$ at the present, which is relevant for the decay of $N_{1}$. Below $T_{\rm osc}$, since $\rho_{\phi}$ drops as a matter energy density, we obtain

[TABLE]

and therefore a sufficiently small mixing to avoid the X-ray 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

[TABLE]

where we used the branching fraction of $N_{2}$ decay for the process $N_{2}\to N_{1}+$(mesons, leptons), ${\rm Br}(N_{2}\to N_{1})\simeq\Gamma(N_{2}\to N_{1})/\Gamma(N_{2}\to SM)\simeq M_{2}\Theta^{2}_{11}\Theta^{2}_{22}/M_{2}\Theta^{2}_{22}\simeq\Theta^{2}_{11}$, 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-$\alpha$ forest constraints which can give the lower limit on the DM mass from the DM free streaming scale $\lambda_{FS}\sim 1~{}\mbox{Mpc}({\rm keV}/{M_{1}})({\langle p/T\rangle}/{3.15})$ kev01 . 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 $\langle p_{1}\rangle\sim 2.8T$, analogous to the conventional (nonresonant) active-sterile oscillation scenario. Taking account of momentum redshifting by a factor $(g_{*}(T_{N_{2}\rightarrow N_{1}})/g_{*}(T\ll{\rm MeV}))^{-1/3}$ due to the change in the effective degrees of freedom, Lyman-$\alpha$ data leads to the RHN DM mass bound $M_{1}\gtrsim 10$ keV for our scenario ir2017 (when $N_{2}\rightarrow 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 Fukugita:1986hr ; Asaka:2005pn ; Pilaftsis:2003gt 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.

Acknowledgements.

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.