1 Introduction↩︎

Among the numerous puzzling issues present in modern particle physics and cosmology are the origins of light neutrino masses [1]–[3] and the nature of dark matter (DM) [4]. Numerous ideas have already been explored to address these problems, but probably one of the most intriguing questions is whether these two issues are correlated. At the very heart of the neutrino mass problem lies the practical issue of testability. It is now well known that the explanation of light neutrino masses introduces a very high scale in the theory (\(\sim 10^{15}\) GeV), which is far beyond the reach of current experiments. Hence, from a model-building point of view, it is desirable that the neutrino mass problem be probed by a lower-scale theory. To address these problems, one must go beyond the Standard Model (SM) [5] and search for a minimal framework that could possibly link these two enigmas without introducing a high scale in the theory. In this spirit, the Type-I seesaw mechanism [6]–[10] seems to be a promising framework, as it introduces three additional heavy right-handed neutrinos (RHNs), \(\sim 10^{15}\) GeV, which are singlets under the SM gauge group. The newly introduced RHNs interact with the light neutrinos and the SM Higgs via Yukawa interactions. These Yukawa interactions not only control the scale of the light neutrino masses but also their interaction strength with the RHNs. It thus seems natural to study these Yukawa interactions and inquire whether they could also play a major role in addressing the DM problem. An elegant approach in this direction is to motivate one of the RHNs as the DM candidate. Numerous studies have been performed in this direction, for instance, in the original \(\nu\)MSM model [11], [12] the lightest RHN is the DM candidate with mass \(\sim\) keV, the Dodelson-Widrow (DW) mechanism [13] which utilizes the effective active-sterile mixing for the DM production, various U(1) gauge symmetric models [14], [15] which study the sterile neutrino production via decays of SM particles as well as newly introduced gauge bosons.

The present study aims to address the phenomenological implications of identifying the lightest RHN (\(N_1\)) as the dark matter candidate with the aid of the neutrino-philic two Higgs doublet model (\(\nu2\)HDM) [16], which provides the low-scale seesaw framework. In the \(\nu2\)HDM, an extra complex scalar doublet (\(\Phi_\nu\)) with a small vev (\(v_{\nu} \sim 10\) MeV) is introduced, which is solely responsible for neutrino mass generation, while the SM Higgs doublet interacts in the usual manner. The small vev opens up the possibility for the heavy RHN masses to attain lower values down to the TeV scale without requiring vanishingly small Yukawa couplings. However, it is not possible to arbitrarily reduce the scale of the RHNs as then the Yukawa couplings must also assume lower values to induce suppressed light neutrino masses. Thus, if \(N_1\) is to be motivated as a light dark matter candidate, the corresponding Dirac Yukawa couplings (\(Y^{\nu}_{i1}\)) must be significantly reduced to respect neutrino observables. This suggests that \(N_1\) should be produced via out-of-equilibrium freeze-in production, which requires feeble couplings of dark matter with the thermal bath, rather than the freeze-out mechanism, which presumes that the dark matter was in thermal equilibrium with the rest of the particles. This class of particles for possible dark matter candidates is collectively known as Feebly Interacting Massive Particles (FIMPs) [17]. Lots of explicit implementations of the FIMP framework have previously been studied [18]–[28], and with a singlet scalar [18], [23], [24] or a new singlet fermion [28] are the two simplest extensions of the SM that incorporate a FIMP because FIMPs are by definition singlets under the Standard Model (SM) gauge group. Under the freeze-in scenario, \(N_1\) gets almost decoupled from the seesaw sector and the remaining two RHNs play the dominant role in neutrino mass generation and could also participate in the creation of the matter-antimatter asymmetry. The dominant production of \(N_1\) proceeds via two-body decay processes of the extra scalars in the setup and also from the decays of the SM gauge bosons (\(W\) and \(Z\)). This possibility was partially introduced in [29], although the crucial gauge bosons decay modes were completely overlooked. The inclusion of these extra production pathways and the study of the resulting phenomenological consequences is the major goal of this work. It is found that the dominant contribution to the comoving number density (\(Y_{N_1}\)) of the dark matter comes from the gauge boson decays and is considerably larger than the combined contributions of all the scalars present in the model. Another interesting outcome of the analysis is that the dark matter relic abundance is independent of the dark matter mass (\(M_{N_1}\)) and depends heavily on the Yukawa couplings. As a result, the mass scale of the lightest active neutrino (\(m_{\nu_{1}}\)) gets fixed by the relic density requirement, thereby correlating the three major issues of neutrino mass, dark matter, and low-scale seesaw. An interesting implication of this correlation is that an enlarged parameter space opens up for the sterile neutrino dark matter mass ranging from sub-keV to \(0.2\) GeV. Considering the latest bounds coming from Lyman-alpha results, it is also found that the setup naturally accommodates both warm and cold dark matter scenarios and provides a natural explanation for the much debated 3.55 keV X-ray line.

Our work is organized as follows. In Section 2, we discuss our model in the context of the nature of dark matter and the freeze-in mechanism. In Section 3, we carefully study the production of FIMP dark matter from the decay of gauge bosons and obtain the corresponding viable parameter space for relic density. Section 4 is dedicated to the discussion on constraints arising from small-scale structure formation. Finally, we present our conclusions in Section 5.

2 THE MODEL↩︎

Originally, the neutrino-philic two Higgs doublet model (\(\nu\)​2HDM) was proposed in [16]. In this model, the scalar sector of the Standard Model (SM) is extended by an additional doublet (the neutrinophilic doublet) \(\Phi_\nu\) with the same quantum numbers as the SM Higgs doublet \(\Phi\). Three heavy right-handed neutrinos (\(N_i\)) are also added in accordance with the standard Type-I seesaw mechanism to realize light neutrino masses. However, the standard Type-I seesaw interaction \(L \tilde{\Phi} N\) is forbidden by imposing a global \(U(1)_L\) symmetry under which \(L_\Phi = 0\), \(L_{\Phi_\nu} = -1\), and \(L_N = 0\). This ensures that \(\Phi_\nu\) couples exclusively to \(N\), while \(\Phi\) couples to quarks and charged leptons as in the SM and is responsible for their masses.

The two scalar doublets as follows: \[\displaystyle \Phi = \begin{pmatrix} \phi^+ \\ \frac{v + \phi^0_r + i \phi^0_i}{\sqrt{2}} \end{pmatrix}, \quad \Phi_\nu = \begin{pmatrix} \phi^+_\nu \\ \frac{v_\nu + \phi^0_{r\nu} + i \phi^0_{i\nu}}{\sqrt{2}} \end{pmatrix}. \label{doublet}\tag{1}\]

The corresponding Higgs potential is then \[\begin{align} V(\Phi,\Phi_\nu) ={}& m^2_\Phi\,\Phi^\dagger\Phi + m^2_{\Phi_\nu}\,\Phi^\dagger_\nu\Phi_\nu + \frac{\lambda_1}{2} (\Phi^\dagger\Phi)^2 + \frac{\lambda_2}{2} (\Phi^\dagger_\nu\Phi_\nu)^2 \nonumber \\ & + \lambda_3 (\Phi^\dagger\Phi)(\Phi^\dagger_\nu\Phi_\nu) + \lambda_4 (\Phi^\dagger\Phi_\nu)(\Phi^\dagger_\nu\Phi) \nonumber \\ & - (\mu^2 \Phi^\dagger\Phi_\nu + \text{h.c.}). \label{pot} \end{align}\tag{2}\]

where a soft symmetry breaking term \(\mu^2\) has also been considered which breaks the \(U(1)_L\) symmetry explicitly but softly. In order to calculate the vev’s of the two Higgs doublets in terms of parameters of the scalar potential, thus the following minimization condition: \[\begin{align} v\, \left[ m^2_\Phi + \frac{\lambda_1}{2} v^2 + \frac{\lambda_3 + \lambda_4}{2} v^2_\nu \right] - \mu^2 v_\nu = 0 \\ v_\nu\, \left[ m^2_{\Phi_\nu} + \frac{\lambda_2}{2} v^2_\nu + \frac{\lambda_3 + \lambda_4}{2} v^2 \right] - \mu^2 v = 0. \end{align}\]

Considering the following choices for the parameters: \[m^2_\Phi < 0,\quad m^2_{\Phi_\nu} > 0, \quad |\mu^2| \ll m^2_{\Phi_\nu},\]

the two vev’s are derived as follows: \[\label{vevs} v \simeq \sqrt{-2 m^2_\Phi/\lambda_1}, \qquad v_\nu \simeq \frac{\mu^2 v}{m^2_{\Phi_\nu} + (\lambda_3 + \lambda_4)v^2/2}.\tag{3}\]

From the Eq. (3 ), it follows that \(v_\nu \sim 1\,\text{GeV}\) is obtained with \(\mu \sim 10\,\text{GeV}\) and \(m_{\Phi_\nu} \sim 100\,\text{GeV}\). An important point to note here is that the only source of \(U(1)_L\) breaking is the \(\mu^2\) term, so the radiative corrections to \(\mu^2\) are proportional to \(\mu^2\) itself and will only be logarithmically sensitive to the cutoff [30]. This ensures that the obtained VEV hierarchy \(v_\nu \ll v\) is stable against radiative corrections [31]–[33]. After spontaneous symmetry breaking (SSB), the physical Higgs bosons are obtained, which contain admixtures of the components of both doublets and are given by [34]:

\[\begin{align} H^{+} &=& \phi^+_\nu\,\cos\beta - \phi^+\,\sin\beta, \quad A = \phi^0_{i\nu}\,\cos\beta - \phi^0_i\,\sin\beta, \\ H &=& \phi^0_{r\nu}\,\cos\alpha - \phi^0_r\,\sin\alpha, \quad h = \phi^0_r\,\cos\alpha + \phi^0_{r\nu}\,\sin\alpha, \end{align}\] and the mixing angles \(\beta\) and \(\alpha\) are given as: \[\tan\beta = \frac{v_\nu}{v}, \qquad \tan 2\alpha \simeq \frac{2 v_\nu v\, [-\mu^2 + (\lambda_3 + \lambda_4) v v_\nu]}{-\mu^2 + \lambda_1 v v_\nu}.\]

After neglecting the terms of \(\mathcal{O}(v_\nu^2)\) and \(\mathcal{O}(\mu^2)\), masses of the bosons are given as follows: \[\begin{align} m^2_{H^+} &\simeq m^2_{\Phi_\nu} + \frac{1}{2}\lambda_3 v^2, \nonumber \\ m^2_A &\simeq m^2_H \simeq m^2_{H^+} + \frac{1}{2}\lambda_4 v^2, \nonumber \\ m^2_h &\simeq \lambda_1 v^2. \end{align}\]

It is to be noted that the mixing angles are suppressed by the small value of \(v_\nu\); hence, \(h\) is almost identically the 125 GeV SM Higgs boson [35], [36]. For simplicity, we adopt a degenerate mass spectrum of \(\Phi_\nu\) as \(m_{H^+} = m_H = m_A = m_{\Phi_\nu}\), and this choice is actually allowed by various constraints [37].

The relevant parts of the Yukawa interaction and the mass term for the heavy right-handed neutrinos are given as:

\[\label{yukint} -\mathcal{L}_Y \supset Y^{\nu}_{ij} \Bar{L_{i}}\tilde{\Phi}_\nu N_{R_{j}} + \frac{1}{2} \Bar{N_{R_{i}}^c} M_{i} N_{R_{i}} + \text{h.c.}\tag{4}\] where \(\tilde{\Phi}_\nu = i \sigma_2 \Phi_\nu^*\), \(L\) is the lepton doublet, and the second term forms the mass matrix of the heavy right-handed neutrinos, which is assumed to be diagonal for simplicity. After spontaneous symmetry breaking (SSB), the above interaction Lagrangian is written as:

\[-\mathcal{L}_Y = \frac{1}{2} \begin{pmatrix} \bar{\nu_L^c} & \bar{N_R} \end{pmatrix} \mathcal{M} \begin{pmatrix} \nu_L \\ N^c_R \end{pmatrix} + \text{H.c.}, \label{eq:Lnu}\tag{5}\] where \(\mathcal{M}\) is the \(6 \times 6\) mass matrix and is written as: \[\mathcal{M} = \begin{pmatrix} 0 & M_D \\ M_D^T & M_R \end{pmatrix}, \label{eq:massmatrix}\tag{6}\] \(M_D\), \(M_R\) are \(3 \times 3\) neutrino Dirac matrix and heavy Majorana neutrino mass matrix respectively. In order to get the eigenvalues, \(\mathcal{M}\) can be diagonalized using a unitary matrix \(\mathcal{U}\) as follows:

\[\mathcal{M}_{\text{diag}} = \mathcal{U}^T \mathcal{M} \mathcal{U}, \label{eq:diag}\tag{7}\] thereby giving the masses of the neutrinos in the physical basis. The mass eigenstates can be defined via the unitary rotation: \[\begin{pmatrix} \nu_L \\ N^c_R \end{pmatrix} = \mathcal{U} \begin{pmatrix} \nu_i \\ N_{Rj} \end{pmatrix}, \label{eq:rotation}\tag{8}\] and the matrix \(\mathcal{U}\) can be expressed as (expanding in terms of \(M_D M_R^{-1}\)) \[\mathcal{U} = \begin{pmatrix} U_{\nu\nu} & U_{\nu N} \\ U_{N\nu} & U_{NN} \end{pmatrix}. \label{eq:Ublocks}\tag{9}\] To leading order, one finds [38], [39] \[\begin{align} \nonumber U_{\nu\nu} &\simeq U_{\text{PMNS}}, \\ \nonumber U_{\nu N} &\simeq M_D^\dagger M_R^{-1}, \nonumber \\ U_{N\nu} &\simeq -M_R^{-1} M_D U_{\nu\nu}, \nonumber \\ U_{NN} &\simeq I. \end{align}\] Using Eq. (7 ) and Eq. (9 ), one can quantify the mixing between different components of the seesaw. The light neutrino mass matrix is given by the seesaw formula [6] as: \[m_\nu = -\frac{v_\nu^2}{2} y m_N^{-1} y^T = U_{\rm PMNS} \hat{m}_\nu U_{\rm PMNS}^T,\] where \(\hat{m}_\nu = \mathrm{diag}(m_1, m_2, m_3)\) is the diagonalized neutrino mass matrix, and \(U_{\rm PMNS}\) is the PMNS (Pontecorvo-Maki-Nakagawa-Sakata) matrix [40]–[42]. It is evident that due to the smallness of \(v_\nu\), a TeV scale of \(m_N\) could lead to light neutrino masses at the \(0.1\,\mathrm{eV}\) scale.

2.1 Sterile Neutrino Dark Matter↩︎

One of the heavy right-handed neutrinos (\(N_1\)) is motivated to be the dark matter candidate by keeping its mass much less than the masses of the \(W\), \(Z\), and \(H\) bosons. This allows for its production via heavy-light mixing (for \(W\) and \(Z\) bosons) and via Yukawa interaction (for all the scalars in the model). The masses of the remaining heavy neutrinos are assumed to lie far above the electroweak scale to satisfy neutrino observables.

A major requirement for the production of a FIMP dark matter [43] is that its interactions with other particles in the thermal bath are feeble. This implies that the dark matter was never in equilibrium with the thermal plasma throughout its thermal evolution, necessitating very small Yukawa couplings of \(N_1\) with the lepton doublet and the second Higgs doublet (\(\Phi_\nu\)).

As a result, \(N_1\) effectively gets completely decoupled and does not participate in neutrino mass generation, leading to only two massive light neutrinos. For calculating light neutrino masses and mixing, the Casas-Ibarra parametrization [38] can be employed, which expresses the Yukawa couplings in terms of experimental data on neutrino masses and mixing angles. Hence, for the Yukawa matrix \(Y^{\nu}\) entries, one can write:

\[m_D = -i U D_{\sqrt{m}} R^T D_{\sqrt{M}}, \label{pmns}\tag{10}\] where \(U\) is the PMNS mixing matrix, \(D_m\) is the diagonal active neutrino mass matrix, \(D_M\) is the diagonal heavy right-handed neutrino mass matrix, and \(R\) is a complex orthogonal matrix which, for the present scenario, can be chosen to be of the following form:

\[R = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos\theta_R & \sin\theta_R \\ 0 & -\sin\theta_R & \cos\theta_R \end{pmatrix}. \label{matix}\tag{11}\] Here, \(\theta_R\) is generally complex. This description guarantees that the setup remains compatible with current neutrino experimental observables of mass-squared differences and all mixing angles. With this choice, \(N_1\) is completely segregated from the rest of the particles in the setup and cannot be produced by any interactions. However, as discussed earlier, feeble couplings for \(N_1\) are desired, and this can be circumvented by perturbing (\(\sim\epsilon_{i1}\)) the matrix \(M_D\), resulting in very small but non-zero entries in the first column. This leads to very small couplings of \(N_1\) with the bath particles, which, if kinematically allowed, will enable production via decay of bath particles.

The tininess of the couplings depends on the stability of the dark matter candidate \(N_1\) and is controlled by relic density requirements. It also dictates the scale of the lightest neutrino mass in this setup. A crucial point not explored in previous studies [29] is that the Yukawa interactions of these heavy right-handed neutrinos also enforce their interactions with Standard Model gauge bosons through heavy-light mixing with active neutrinos. The ignorance of these interactions is often justified due to the large masses of heavy right-handed neutrinos, rendering the heavy-light mixing \(V_{ij} \sim M_{D_{ij}}/M_j\) extremely small to play any significant role in dark matter requirements. However, as will be shown, under the FIMP scenario, such small values of \(V_{ij}\) are valuable and provide an alternative production channel for the dark matter candidate.

Thus, we write below the interactions [44] of the heavy right-handed neutrinos with light neutrinos and the Standard Model gauge bosons by utilizing the heavy-light mixing:

\[\label{gauge} \begin{align} - \mathcal{L}_{g} \subset & \frac{g}{\sqrt{2}} W_{\mu}^{+} \sum_{i, j=1}^3 [\Bar{N}_{i}^c (V^{\dagger})_{ij} \gamma^{\mu} P_{L} l_{j} ] + \frac{g}{2 C_{\theta_{w}}} Z_{\mu} \\ & \times \sum_{i, j=1}^3 [\Bar{\nu_{i}} (U^{\dagger} V)_{ij} \gamma^{\mu} P_{L} N_{j}^c + \Bar{N_{i}}^c (V^{\dagger}V)_{ij} \gamma^{\mu} P_{L} N_{j}^c ] \end{align}\tag{12}\] where \(P_{L} = \frac{1 - \gamma^5}{2}\), \(g\) is the gauge coupling constant, \(C_{\theta_{w}} = \cos\theta_{w}\), \(\theta_{w}\) is the weak mixing angle, \(U\) is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, and \(V_{ij} = \frac{M_{D_{ij}}}{M_{N_j}}\) is the heavy-light mixing matrix element through which the heavy right-handed neutrinos interact with the \(W\) and \(Z\) bosons.

2.2 A Decaying Dark Matter↩︎

In order for \(N_1\) to be a successful dark matter candidate, apart from satisfying relic abundance requirements, it should be stable over cosmological time scales. The best approach is to make the dark matter completely stable, meaning no decay channels should be available. However, in the present scenario, \(N_1\) can decay through different channels due to its mixing with light neutrinos via heavy-light mixing. This leads to a decaying dark matter scenario, which generally introduces complications; nevertheless, a distinct advantage is that the decay products could enhance the detectability prospects of the dark matter candidate.

Hence, heavy-light mixing plays a key role in producing the dark matter, destabilizing it, and enhancing the testability of the setup. Regarding the decay products, depending on the mass of \(N_1\), the only possibilities are charged leptons (mainly electrons and muons), light neutrinos, and photons. The different decay possibilities, including three-body decays, are given as follows:

(a) via off-shell \(W/Z\): \[N_1 \to \ell^-_1 \ell^+_2 \nu_{\ell_2},\, N_1 \to \ell^- q_1 \bar{q}_2, \, N_1 \to \ell^- \ell^+ \nu_\ell, \, N_1 \to \nu_\ell \bar{\ell'} \ell', \,\] \[N_1 \to \nu_\ell q \bar{q}, \quad N_1 \to \nu_\ell \nu_{\ell'} \bar{\nu}_{\ell'}, \quad N_1 \to \nu_\ell \nu_\ell \bar{\nu}_\ell ;\] (b) via off-shell Higgs: \[N_1 \to \nu_\ell \bar{\ell} \ell;\] (c) radiative decay : \[N_1 \to \gamma \nu.\] Ensuring that that the expected lifetime of \(N_1\) must be greater than the age of the Universe these decay channels will lead to strong constraints on the heavy-light mixing. It turns out that the most stringent constraint is obtained from the radiative decay to photon and light neutrino (c), which translates into a bound on the active-sterile neutrino mixing \(V_{i1}\) as [45]–[48] \[\theta^2_{1} = \sum_{i=1}^{3} |V_{i1}|^2 \leq 2.8 \times 10^{-18} \left( \frac{\text{GeV}}{M_1} \right)^{5}. \label{eq:theta95bound}\tag{13}\] So, in the subsequent analysis of dark matter the above equation is used as a strict bound for the combination of \(N_1\) mass and the corresponding Yukawa couplings. In Fig. 1, we have shown the allowed (white region) and the excluded region (blue) based on the above equation.

2.3 Non-thermal Regime↩︎

One of the major requirement of the freeze-in mechanism (FIMP) is that the dark matter particle should not reach thermal equilibrium with the rest of the particles present in the thermal bath. In the setup under consideration, only the gauge singlet heavy right handed neutrinos could play the role of FIMPs. Thermal equilibrium in the early universe can be prevented only if their Yukawa interactions are sufficiently suppressed. In this section, we analyze the different processes that can produce these singlets and obtain the necessary conditions for the out of equilibrium requirement. We show that only one of them, which we choose to be \(N_1\), can play the role of a FIMP.

All the right handed neutrinos have the Yukawa interactions of the form \(N_{k}\,L \,\phi_{\nu}\) so, they can be produced via the two-body decay of the scalars present in the model and the decay width is approximately given by \(\Gamma(H_2\to N_1\,L)=M_{H_2}\,y_1^2/(8\pi)\) (see Eqs. (16 ) and (17 ) below). However, apart from the scalar decay modes, two crucial decay channels involving W and Z bosons also exist which utilizes the mixing between \(N_1\) and the light neutrinos. The fastest one out of all these available decay channels is that of Z boson decay whose decay width is given as: \(\Gamma_{Z \to N \nu_i + N \bar{\nu}_i} = \frac{1}{48\pi} m_Z |Y_{\nu i}|^2 f(m_N^2 / m_Z^2)\) where \(f(x) = (1 - x)^2 (1 + 2/x)\). Then the out of equilibrium condition is given as: \[\begin{align} \Gamma(Z \to N \nu_i + N \bar{\nu}_i)\lesssim H(T\sim M_{Z})\,, \end{align}\]

which gets translated to the following condition on the Yukawa couplings of \(N_1\) as: \[\label{ybound} \sum_{\alpha=e,\mu,\tau} |Y_{\alpha1}|^2 \lesssim 1 \times 10^{-16} \left( \frac{M_{N_{1}}}{10~\mathrm{GeV}} \right)^2,\tag{14}\]

Hence, if \(|Y_{\alpha1}|^2\) were larger than this value, \(N_1\) would be produced quite abundantly and will reach thermal equilibrium. These small values of the Yukawa couplings of the first column of \(M_D\) implies that, \(N_1\) plays no role in the mass generation of active neutrinos. This conclusion is in complete accordance with that of [49], [50], where the authors have argued that in low-scale Type-1 seesaw models with three heavy neutrinos if the mass of lightest active neutrino \(< 10^{-3}\) eV then, only one of the heavy sterile states never thermalize. However, as discussed in subsequent sections, the requirement of the correct relic density of dark matter will make the mass of the lightest neutrino to attain even further smaller values. It must be pointed out that if \(N_1\) was never in thermal equilibrium then the only mechanism of its production in the early universe is via freeze-in. However, \(N_1\) can also be produced via the decay of the heavier singlets (\(N_{2,3}\)) or via \(2\to 2\) scatterings of SM leptons or the scalars in the setup. However, all these processes are rather subdominant and do not modify the above equilibrium conditions.

Models of sterile neutrino are ought to have the active-sterile mixing (or, this is just the heavy-light mixing as in Type-1 seesaw) and then apart from the production via decay of heavy particles (freeze-in scenario) which we are considering in this work another mechanism also exists which also utilizes this heavy-light mixing. Sterile neutrinos are produced through non-resonant oscillations of active neutrinos. As a result, a primordial population of sterile neutrinos is obtained with a warm dark matter spectrum unlike the usual cold dark matter candidates. The contribution to the relic abundance of sterile neutrino is given as [51]:

\[\begin{align} \Omega_{\text{DW}} h^2 \approx 0.3 \times \left( \frac{\sin^2 2\theta}{10^{-10}} \right) \left( \frac{M_{N_1}}{100 \text{ keV}} \right)^2 \end{align}\]

where \(\theta\) is the active-sterile mixing. However, this mechanism of sterile neutrino production (Dodelson-Widrow [13] is strongly constrained by Lymna-alpha bounds [45] and also X-ray observations [52] and it is now well accepted [53], [54] that this production pathway cannot account for all the observable dark matter in the Universe. But, due to the presence of active-sterile mixing this contribution to the relic density always exists and must be taken into account. For our case, we have found that for the \(N_1\) mass range \(1\)keV - \(0.1\) GeV and the Yukawa couplings of the order \(\sim 10^{-14}\) the contribution is quite less and hence irrelevant in our study.

So, the conclusion is that in order to simultaneously account for light neutrino masses and asking for a TeV scale of seesaw only one FIMP candidate is allowed. In the next section, we will discuss the various production channels in detail and also show that it is rather easy to account for the observed dark matter relic density.

3 Dark Matter Production↩︎

In this section, we discuss the various available production channels for the dark matter candidate \(N_1\). Since \(N_1\) has a direct coupling to the lepton doublet and to the neutrinophilic doublet, the relevant production pathway is by the decays of the scalars (\(H^0,A^0,H^\pm\)) which are in equilibrium with the thermal bath. The SM gauge bosons \(W\) and \(Z\) will also take part in the production mechanism with the aid of heavy-light mixing. The \(N_1\) yield, \(Y_{N_{1}}(T)=n_{N_1}(T)/s(T)\), is computed by solving the following Boltzmann equation [14]

\[\begin{align} \frac{dY_{N_1}}{dz} &= \frac{2 M_{pl} z}{1.66 m_h^2} \frac{\sqrt{g_*(z)}}{g_s(z)} \Big[ Y_Z^{eq} \big\langle\Gamma_{Z\to N_1 \nu}\big\rangle + Y_H^{eq} \big\langle\Gamma_{H\to N_1 \nu}\big\rangle \nonumber \\ & + Y_k^{eq} \big\langle\Gamma_{k\to N_1 \nu}\big\rangle \nonumber + Y_{K^{0}}^{eq} \big\langle\Gamma_{K^{0}\to N_1 \nu}\big\rangle \nonumber \\ & + Y_{K^{\pm}}^{eq} \big\langle\Gamma_{K^{\pm}\to N_1 \nu}\big\rangle + Y_W^{eq} \big\langle\Gamma_{W^{\pm}\to N_1\ell^{\pm}}\big\rangle\Big],\label{1BE} \end{align}\tag{15}\] with \(z = m_h/T\) and the quantity \(\big\langle \Gamma_{A\to BC}\big\rangle\) represents the thermally averaged decay width [55] and is defined as: \[\begin{align} \big\langle \Gamma_{A\to BC}\big\rangle = \frac{K_1(z)}{K_2(z)} \Gamma_{A\to BC} \end{align}\]

where \(K_1(z)\) and \(K_1(z)\) are the modified Bessel functions of order 1 and 2 respectively. The function \(g_*(z)\) is given by:

\[\sqrt{g_*(z)} = \frac{g_s(z)}{\sqrt{g_\rho(z)}}\left(1 - \frac{1}{3} \frac{d \ln g_s(z)}{d \ln z}\right)\]

where, \(g_\rho(z)\) and \(g_s(z)\) are the effective degrees of freedom related to energy density (\(\rho\)) and the entropy density (s) of the universe respectively. The decay rates that enter into this expression are given as  [29], [56]: \[\begin{align} \label{GS1} \Gamma\left(H^{0}/A^{0} \to N_{1} \, \nu_{\alpha} \right) & = \frac{m_{H^{0}/A^{0}}\,\left| Y^{\nu}_{\alpha 1} \right|^{2}}{32\,\pi}\left(1-\frac{M_{1}^{2}}{m_{H^{0}/A^{0}}^{2}}\right)^{2}\;\nonumber,\\& \approx \frac{m_{H^{0}/A^{0}}\,\left| Y^{\nu}_{\alpha 1} \right|^{2}}{32\,\pi},\nonumber \\ \end{align}\tag{16}\]

\[\begin{align} \label{GS2} \Gamma\left(H^{+} \to N_{1} \, \overline{\ell_{\alpha}} \right) & = \frac{m_{H^{+}}\,\left| Y^{\nu}_{\alpha 1} \right|^{2}}{16\,\pi}\left(1-\frac{M_{1}^{2}}{m_{H^{+}}^{2}}\right)^{2}\;\nonumber,\\& \approx \frac{m_{H^{+}}\,\left| Y^{\nu}_{\alpha 1} \right|^{2}}{16\,\pi} \nonumber.\\ \end{align}\tag{17}\] \[\begin{align} \tag{18} \Gamma_{W^{\pm} \to N \ell^{\pm}_i} & =& \frac{1}{48\pi} m_W |Y_{\nu i}|^2 f(m_N^2 / m_W^2), \\ \Gamma_{Z \to N \nu_i + N \bar{\nu}_i} &=& \frac{1}{48\pi} m_Z |Y_{\nu i}|^2 f(m_N^2 / m_Z^2), \tag{19} \end{align}\] where \(f(x) = (1 - x)^2 (1 + 2/x).\) The approximations used in Eq. (16 ) and Eq. (17 ) are valid unless there is a strong mass degeneracy between \(N_1\) and one of the scalars. Note that the annihilations producing \(N_1\) are significantly suppressed (\(\sim \epsilon_i^4\)) compared to decays (\(\sim \epsilon_i^2\)), and hence are excluded from our analysis. Back reactions involving \(N_1\) are also not included since the \(N_1\) number density is vanishingly small to start with. For the same reason, terms proportional to \(Y_{N_1}\) are also dropped, which is a standard approximation for the freeze-in case [17], [57].

Now, in order to compute the relic abundance (\(\Omega_{N_{1}} h^2\)) of the sterile neutrino dark matter we need to find the value of its co-moving number density (\(Y_{N_{1}}\)) at the present epoch (\(T=T_{\infty}\), \(T_{\infty}=2.73\)K). This value (\(Y_{N_{1}}(T=T_{\infty})\)) is obtained by solving the Boltzmann eq for the number density of \(N_1\). The expression of (\(\Omega_{N_{1}} h^2\)) in terms of \(Y_{N_{1}}(T=T_{\infty})\) is given as [58]:

\[\label{relicbound} \Omega_{N_1}h^2 = 2.755\times 10^8 \bigg(\frac{M_1}{\text{GeV}}\bigg)Y_{N_1}(z_{\infty})\tag{20}\]

We have studied quantitatively the freeze-in production of dark matter (\(N_{1}\)) in this scenario by numerically solving the Boltzmann Eq. (15 ). We start with the standard assumption of the freeze-in scenario that at high temperatures the comoving number density (\(Y_{N_{1}}\)) of dark matter is zero. So, the Boltzmann equation is solved with the initial condition \(Y_{N_1}=0\) for \(T \gg m_H\).

In Fig. 1, we show the variation of the comoving number density (\(Y_{N_{1}}\)) of \(N_1\) versus \(x\), where \(v_{\nu}=20\) MeV, \(M_{N_{1}}=0.1\) MeV, and the corresponding Yukawa couplings \(Y^{\nu}_{1i} \sim 10^{-14}\). The total \(Y_{N_{1}}\) is obtained by solving the complete Boltzmann Eq. (15 ), considering all relevant production pathways of \(N_1\), which consists of two-body decays from all the scalars, \(W\), and \(Z\) bosons. From the nature of the curves, we can observe the standard FIMP behavior: at high temperatures, \(N_1\) has a vanishing \(Y_{N_{1}}\), but as the Universe expands and cools down (i.e., with increasing \(x\)), \(N_1\) slowly gets produced over time and \(Y_{N_{1}}\) gradually increases, with maximum production occurring around \(x \sim 2-5\). As the temperature drops further to lower values, \(Y_{N_{1}}\) gets “freezes-in" and thereby attains a constant value.

For the chosen values of the parameters, the solid maroon line corresponds to the correct relic abundance within the \(3 \sigma\) of DM relic density  [59], [60]. The contributions to \(Y_{N_{1}}\) from specific processes are also shown in Fig. 1. It is clearly visible that \(Y_{N_{1}}\) receives a dominant contribution from the decays of \(W\) and \(Z\) bosons, depicted as dashed green and dashed blue lines, respectively. The horizontal dashed black line corresponds to the contribution from the decays of all scalars combined, which is negligible in comparison to \(W\) and \(Z\) decays.

Hence, despite the presence of some extra heavy scalars in the model coming from linear combinations of the two complex scalar doublets, the major contribution proceeds via gauge boson decays. This feature was overlooked in previous studies regarding sterile neutrino production in the neutrino-philic two Higgs doublet model  [29]. This is because the production via gauge bosons utilizes the heavy-light mixing, which is generally very small as it is suppressed by heavy right-handed neutrinos. However, under the FIMP scenario, the couplings of the dark matter particle are indeed required to be quite small in order to respect the non-thermality condition, so \(W\) and \(Z\) decays should not be excluded from the analysis; in fact, as demonstrated above, they play the major role in \(N_1\) production.

In Fig. 2, we show the variation of relic abundance versus \(x\) with relative contributions coming from gauge bosons and all scalars combined for the previously chosen parameter values. It is clear from the curves that the major contribution comes from the gauge bosons, depicted by the dashed green line for the \(W\) boson and the dashed blue line for the \(Z\) boson, with the \(W\) boson contribution slightly leading over the \(Z\) boson at any given temperature. The combined scalar contribution is shown by the dashed black line, whose contribution to the relic density at a given temperature is quite low, as noted previously.

We have also depicted the total relic density by the solid maroon line, summing over all contributions, and noted that it matches the required relic abundance shown by the horizontal solid red line.

From Eq. (20 ), one notes that the relic abundance is evaluated to be proportional to the mass of the dark matter candidate (\(M_{N_{1}}\)). However, information about the relevant Yukawa couplings remains hidden in the decays which enter the Boltzmann equation. To study the intricate relationship between relic abundance, \(N_1\) mass, and the Yukawa couplings of the first column of \(M_D\) (\(Y^{\nu}_{i1}\)), we vary the Yukawa couplings and \(M_{N_1}\) simultaneously and inquire about the combinations giving rise to the correct relic abundance.

Generally, one would expect many points spanning an ample region of the parameter space. However, on the contrary, we find a straight line corresponding to those values of (\(Y^{\nu}_{i1}, M_{N_1}\)) that yield the appropriate relic satisfaction. We also note the rather interesting result that not all values of the Yukawa couplings lead to the correct relic, implying that there is not enough freedom in varying the Yukawa couplings. Nevertheless, there is substantial scope over the choices of \(M_{N_{1}}\) values, ranging from sub-keV to a few hundred MeV. This behavior suggests that the relic abundance of dark matter depends heavily on the Yukawa couplings (\(Y^{\nu}_{i1}\)) and is almost insensitive to the mass of the dark matter.

We would like to note that this can be explained by the fact that \(Y^{\nu}_{i1}\) is directly responsible for the production of \(N_1\). From the seesaw relation, it follows that \(Y^{\nu}_{i1} \propto \sqrt{m_{\nu_1} M_{N_1}}\), where \(m_{\nu_1}\) is the mass of the lightest neutrino.

Major production occurs via two-body gauge boson decays, which, based on Eqs. (18 ) and (19 ), show that \(\Gamma_{W,Z} \propto \frac{m_{\nu_1}}{M_{N_1}}\). Using analytical estimates  [17], [29] it also follows that \(Y_{N_1}(z_{\infty}) \propto \frac{m_{\nu_1}}{M_{N_1}}\). From Eq. (20 ), the relic abundance \(\Omega_{N_1} h^2 \propto M_{N_1} Y_{N_1}(z_{\infty})\), which, based on the above arguments, is directly determined by the lightest neutrino mass scale \(m_{\nu_1}\) and is independent of \(M_{N_1}\). This explains clearly the observed straight-line behavior in Fig. 3 and the limited variation of the Yukawa couplings.

In Fig. 4, we show the contour plot for the relic abundance in the \(\theta_1\)-\(M_{N_{1}}\) plane, with the solid red line. The blue region denotes the excluded part based on the constraint Eq. (14 ) for the radiative decay of the sterile neutrino. This plot depicts that not all values of the pair (\(Y^{\nu}_{i1}, M_{N_1}\)) are allowed and that they are highly constrained by the stability requirements of \(N_1\). It also shows the complete allowed parameter space for the \(N_1\) mass, which ranges from keV to \(0.2\) GeV.

Although \(N_1\) masses below the keV scale are allowed by the setup, they are excluded due to the Tremaine-Gunn bound  [61] on sterile neutrino mass. Thus, the lower limit for the \(M_{N_{1}}\) mass is taken to be \(1\) keV, while the maximum allowed limit is found to be around \(0.2\) GeV, clearly depicting an enlarged parameter space for the sterile neutrino mass.

In order to further probe the nature of this dark matter candidate and place it on firm ground, constraints arising from small-scale structure formation must be taken into consideration. This will not only further constrain the model but also provide information on whether \(N_1\) constitutes hot, warm, or cold dark matter, and is discussed in the following section.

4 Constraints from small-scale structure formation↩︎

For a truly viable dark matter analysis, obtaining the correct relic abundance is generally not sufficient. Dark matter must also allow structure formation  [62], [63], such as the formation of galaxies. For instance, dark matter should not be too hot as it would disturb cosmological perturbations and inhibit structure formation in the early Universe. However, in the standard picture, a keV-scale dark matter candidate is motivated to be a warm dark matter candidate.

Despite this, simply having a mass in the keV range and satisfying relic density requirements does not guarantee warm dark matter. This is because the capacity to form small structures depends on two major factors: (i) the mass of the dark matter and (ii) the momentum distribution function of the dark matter candidate. The distribution functions depend on the production mechanism and the details of cosmological evolution.

Thus, to determine whether the assumed candidate falls under the category of warm, hot, or cold dark matter, one needs to precisely find the corresponding distribution function and constrain it with astrophysical data  [64]. However, this is a highly non-trivial task, and various approximations are generally employed.

As an approximate measure, a parameter called the free-streaming horizon (\(\lambda_{FS}\)) is introduced  [45], through which dark matter is categorized as hot, warm, or cold. Physically speaking, \(\lambda_{FS}\) represents the average collision-free distance traveled by a dark matter particle after production, before gravitational clustering effects become significant. In order to calculate \(\lambda_{FS}\), one requires the distribution function of the dark matter.

The precise form of the distribution was calculated in  [29] for the neutrino-philic two Higgs doublet model, and the corresponding relationship between the free-streaming horizon and mass is given by:

\[r_{\text{fs}} = 0.047 \text{ Mpc} \times \left(\frac{10 \text{ keV}}{m_{N_1}}\right) \label{eq:jhepformula}\tag{21}\]

which is obtained by considering only the scalars decay and under the assumption of degenerate masses of extra scalars. However, in this study we have extra decay processes coming from W and Z bosons so, roughly speaking, the distribution function and hence the treatment of the free streaming horizon should change. However, considering new decay modes for \(N_1\) production only leads to addition of extra distribution functions corresponding to these new modes. So, the total distribution function of \(N_1\) is just the sum of the various distribution functions calculated for different decay modes and only the coefficients of the distribution functions will change while the polylogarithmic structure will remain intact. Moreover, the phase-space structures of the scalar, \(W\), and \(Z\) boson decay modes don not differ much as they produce similar daughter particles. This similarity is crucial because if the phase space structures of the various decay processes were diverse, the distribution function could have distinct maxima and it becomes difficult to define mean values, rendering the calculation of the free-streaming horizon futile. For a detailed discussion on these issues, the interested reader is referred to the discussions in [65], [66]. The classification of sterile neutrino dark matter via the free-streaming horizon is shown in Fig. 5, along with the latest bounds coming from Lyman-\(\alpha\) forest measurements [67]. The following points are noted:

• As noted from Fig. 4, the model predicts a light dark matter mass \(\sim 1\) keV, while from Eq. (21 ) one obtains \(r_{\text{fs}} = 0.4\) Mpc, corresponding to the hot dark matter regime, which is excluded by all observations. For hot dark matter, one must have \(r_{\text{fs}} > 0.1\) Mpc; hence, the region \(m_{N_1} < 4.7\) keV is excluded and is indicated by the beige colored region in Fig. 5.

• The region \(4.7 < m_{N_1} < 33.5\) keV also gives the desired relic abundance (see Fig. 4). Using Eq. (21 ), this interval corresponds to \(0.014 < r_{\text{fs}} < 0.1\) Mpc, which characterizes warm dark matter. Thus, the model accommodates a warm dark matter candidate if \(M_{N_1}\) lies within this range. The yellow region in Fig. 5 depicts the allowed interval for warm dark matter as per Eq. (21 ); however, recent Lyman-\(\alpha\) forest measurements further constrain this region, with the latest bound shown by the solid green line providing an updated lower limit on the warm dark matter mass. An interesting feature of the setup is its feasibility to explain the claimed X-ray line at 3.55 keV, since \(M_{N_1} = 7.1\) keV is allowed by the relic density and corresponds to warm dark matter. This is clearly shown in Fig. 5 by a vertical dashed orange line.

• Lastly, for the case \(r_{\text{fs}} < 0.014\) Mpc, Eq. (21 ) implies \(m_{N_1} > 33.6\) keV. This regime, up to 0.2 GeV, can explain all the observed dark matter in the Universe and corresponds to cold dark matter. This region is shown in Fig. 5 as the cyan region and is allowed by all observational data.

So, as one goes to higher mass scale of \(N_1\), it follows from Eq. (21 ) that one obtains colder dark matter candidate as a result of decreasing \(r_{FS}\) (Table-I).However, a novel feature of this setup is that it is able to accommodate both warm and cold dark matter options and an explanation [68], [69] for the much talked \(3.55\) keV X-ray line.

5 Discussion and Conclusions↩︎

Considering the neutrino-philic two Higgs doublet model, we have shown that one can simultaneously account for the masses of the active neutrinos, demand a low-scale seesaw, and explain the dark matter problem, where exactly one of the singlet fermions, \(N_1\), can be out of equilibrium in the early Universe and act as a FIMP. We discussed important phenomenological implications of this model that were overlooked in previous studies, especially the dominant role of gauge bosons in dark matter production. Due to the FIMP nature of \(N_1\), it becomes almost completely decoupled from the thermal bath; as a result, one of the light neutrinos is extremely light, with its mass scale being fixed by relic density requirements. Specifically, we studied the evolution of the comoving number density of \(N_1\) as the Universe expands and demonstrated the dominance of \(W\) and \(Z\) boson decays over scalar decays. The same feature was observed in the evolution of the relic abundance, with the \(W\) boson contribution slightly leading over the \(Z\) boson, while the scalar contribution remained negligible. The respective values of the pair (\(Y^{\nu}_{i1}\), \(M_{N_1}\)) were scrutinized by the relic requirements, and it was concluded that the relic density is almost independent of the mass of the dark matter while the Yukawa couplings are tightly constrained. The model predicts a wide range for \(N_1\) masses, spanning sub-keV to \(0.2\) GeV, owing to the small vev (\(v_\nu \sim 10\) MeV) of the new doublet, which helps widen the available parameter space and accounts for all observable dark matter present in the Universe. According to the latest Lyman-\(\alpha\) results, the model accommodates both warm and cold dark matter options and provides a plausible explanation for the debated \(3.55\) keV X-ray line. We point out that for the mass range \((1-10)\) MeV, the sterile neutrino will undergo three-body decays via off-shell \(W\) and \(Z\) bosons, producing an electron-positron pair and a light neutrino. The annihilation of the electron-positron pair could offer a feasible explanation for the observation of the \(511\) keV line seen by INTEGRAL/SPI and serves as a smoking-gun signature for the non-thermal cold dark matter option present in this model. Other pathways towards testability stem from the prediction of a very light scale for one of the active neutrinos, which could falsify the setup through ongoing (or future) experiments such as KATRIN [70] and PROJECT-8 Collaborations  [71].

Acknowledgments↩︎

We thank Prof. Rathin Adhikari for valuable guidance. S.K. also gratefully acknowledge Prof. Sushant Ghosh for providing the necessary computational facility used in this work.

References↩︎