FIMP dark matter candidate(s) in a $B-L$ model with inverse seesaw mechanism

TL;DR

This paper explores a $B-L$ extended standard model with inverse seesaw mechanism, proposing a fermionic FIMP dark matter candidate produced via freeze-in, and examines both single and two-component dark matter scenarios.

Contribution

It introduces a novel $B-L$ model with inverse seesaw featuring a fermionic FIMP dark matter candidate and analyzes its relic density and production mechanisms.

Findings

KeV FIMP DM can account for relic density via $Z'$ decay.

Annihilation processes can significantly contribute to DM production.

Two-component DM scenario broadens viable parameter space.

Abstract

The non-thermal dark matter (DM) production via the so-called freeze-in mechanism provides a simple alternative to the standard thermal WIMP scenario. In this work, we consider a popular $U(1)_{B-L}$ extension of the standard model (SM) in the context of inverse seesaw mechanism which has at least one (fermionic) FIMP DM candidate. Due to the added $\mathbb{Z}_{2}$ symmetry, a SM gauge singlet fermion, with mass of order keV, is stable and can be a warm DM candidate. Also, the same $\mathbb{Z}_{2}$ symmetry helps the lightest right-handed neutrino, with mass of order GeV, to be a stable or long-lived particle by making a corresponding Yukawa coupling very small. This provides a possibility of a two component DM scenario as well. Firstly, in the absence of a GeV DM component (i.e., without tuning its corresponding Yukawa coupling to be very small), we consider only a keV DM as a single component DM, which is produced by the freeze-in mechanism via the decay of the extra $Z'$ gauge boson associated to $U(1)_{B-L}$ and can consistently explain the DM relic density measurements. In contrast with most of the existing literature, we have found a reasonable DM production from the annihilation processes. After numerically studying the DM production, we show the dependence of the DM relic density as a function of its relevant free parameters. We use these results to obtain the parameter space regions that are compatible with the DM relic density bound. Secondly, we study a two component DM scenario and emphasize that the current DM relic density bound can be satisfied for a wide range of parameter space.

Tables

Table 1. Table 1: Complete particle spectrum and their corresponding charges under various gauge groups.

Gauge Group $SU{\left(2\right)}_L$ $U{\left(1\right)}_Y$ $U{\left(1\right)}_{B-L}$

Baryon Fields $Q_L^i={(u_L^i,d_L^i)}^T$ $u_R^i$ $d_R^i$ $2$ $1$ $1$ $1/6$ $2/3$ $-1/3$ $1/3$ $1/3$ $1/3$

Lepton Fields $L_L^i={({\nu }_L^i,e_L^i)}^T$ $e_R^i$ $N_R^i$ $S_1^i$ $S_2^i$ $2$ $1$ $1$ $1$ $1$ $-1/2$ $-1$ $0$ $0$ $0$ $-1$ $-1$ $-1$ $-2$ $2$

Scalar Fields    ${\varphi }_h$ ${\varphi }_H$ $2$ $1$ $1/2$ $0$ $0$ $-1$

Equations

$$\Omega h^2=0.1199\pm 0.0027\mathrm{at}68\%\mathrm{CL},$$

$$\mathcal{L}$$

$$D_{\mu }={\partial }_{\mu }-i{g}_c{T}^{\alpha }{G}_{\mu }^{\alpha }-ig{\tau }^a{W}_{\mu }^a-i{g}_{Y}Y{B}_{\mu }-i{g}^{\prime }{Y}_{BL}{B}_{\mu }^{\prime },$$

$$\mathcal{V}({\varphi }_h,{\varphi }_H)={\mu }_h^2\text{leavevmode}{\varphi }_h^{†}{\varphi }_h+{\mu }_H^2\text{leavevmode}{\varphi }_H^{†}{\varphi }_H+{\lambda }_h({\varphi }_h^{†}{\varphi }_h{)}^2+{\lambda }_H({\varphi }_H^{†}{\varphi }_H{)}^2+{\lambda }_{hH}({\varphi }_h^{†}{\varphi }_h)({\varphi }_H^{†}{\varphi }_H),$$

$${\mu }_h^2<0,{\mu }_H^2<0,{\lambda }_h\ge 0,{\lambda }_H\ge 0\mathrm{and}{\lambda }_{hH}\ge -2\sqrt{{\lambda }_h{\lambda }_H}.$$

$${\varphi }_h=\left(\begin{array}{c}0\\ \frac{v+h}{\sqrt{2}}\end{array}\right),{\varphi }_H=\frac{v^{\prime }+H}{\sqrt{2}},$$

$${\mathcal{L}}_m^{\nu }={\mu }_S({\bar{S}}_1^c{S}_1+{\bar{S}}_2^c{S}_2)+(M_D{\bar{\nu }}_L{N}_R+M_N{\bar{N}}_R^c{S}_2+h.c.),$$

$$\displaystyle\mathcal{M}_{\nu}=\left(\begin{array}[]{cccc}0\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ M_{D}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\\ M_{D}^{T}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ M_{N}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\\ 0\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ M^{T}_{N}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ \mu_{S}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\\ 0\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ \mu_{S}\end{array}\right).$$

$$M_{S_1}={\mu }_S.$$

$$M_{{\nu }_l}$$

$$M_{{\nu }_{H,H^{\prime }}}^2$$

$$\displaystyle\mathcal{M}^{2}_{\rm scalar}=\left(\begin{array}[]{cc}2\lambda_{h}v^{2}\nobreak\leavevmode&\leavevmode\nobreak\ \lambda_{hH}\,v^{\prime}\,v\\ \lambda_{hH}\,v^{\prime}\,v\nobreak\leavevmode&\leavevmode\nobreak\ 2\lambda_{H}v^{\prime 2}\end{array}\right).$$

$$h_1$$

$$h_2$$

$$\tan 2\alpha =\frac{{\lambda }_{hH}{v}^{\prime }v}{{\lambda }_H{v}^{\prime 2}-{\lambda }_h{v}^2}.$$

$$M_{h_{1,2}}^2$$

$${\lambda }_h$$

$${\lambda }_H$$

$${\lambda }_{hH}$$

$$\frac{{\Gamma }_{Z^{\prime }}}{H(T=M_{Z^{\prime }})}<1,$$

$$\hat{L}[f]=\mathcal{C}[f],$$

$$\hat{L}=\frac{\partial }{\partial t}-Hp\frac{\partial }{\partial p},$$

$$z=\frac{M_{\mathrm{sc}}}{T},{\xi }_p={\left[\frac{g_s(T_0)}{g_s(T)}\right]}^{1/3}\frac{p}{T},$$

$$\frac{dT}{dt}=-HT{\left(1+\frac{T{g}_s^{\prime }(T)}{3g_s(T)}\right)}^{-1},$$

$$\hat{L}=zH{\left(1+\frac{T{g}_s^{\prime }}{3g_s}\right)}^{-1}\frac{\partial }{\partial z},$$

$$\hat{L}{f}_{Z^{\prime }}=\sum _{i=1,2}{\mathcal{C}}^{h_i\to Z^{\prime }{Z}^{\prime }}+{\mathcal{C}}^{Z^{\prime }\to \mathrm{all}},$$

$$Y_{Z^{\prime }}=\frac{n_{Z^{\prime }}(z)}{s}=\frac{45g}{4{\pi }^4{g}_s(M_{\mathrm{sc}}/z_0)}{\int }_0^{\infty }d{\xi }_p{\xi }_p^2{f}_{Z^{\prime }}({\xi }_p,z),$$

$$\frac{d{Y}_{S_1^1}}{dz}$$

$$⟨{\Gamma }_{Z^{\prime }\to S_1^1{S}_1^1}{⟩}_{\mathrm{NTH}}=M_{Z^{\prime }}{\Gamma }_{Z^{\prime }\to S_1^1{S}_1^1}\frac{{\text{bigints}}_0^{\infty }d{\xi }_p{\xi }_p{f}_{Z^{\prime }}({\xi }_p,z){\left({\left[\mathcal{B}(z)\frac{M_{\mathrm{sc}}}{z}\right]}^2+{\left[\frac{M_{Z^{\prime }}}{{\xi }_p}\right]}^2\right)}^{-1/2}}{{\text{bigints}}_0^{\infty }d{\xi }_p{\xi }_p^2{f}_{Z^{\prime }}({\xi }_p,z)},$$

$$\mathcal{B}(z)={\left[\frac{g_s(M_{\mathrm{sc}}/z)}{g_s(M_{\mathrm{sc}}/z_0)}\right]}^{1/3}.$$

$$⟨{\sigma v}_{f\bar{f}\to S_1^1{S}_1^1}⟩=\frac{z}{8M_f^4{M}_{\mathrm{sc}}{K_2(z{M}_f/M_{\mathrm{sc}})}^2}{\int }_{4M_f^2}^{\infty }{\sigma }_{f\bar{f}\to S_1^1{S}_1^1}\text{leavevmode}\left(s-4M_f^2\right)\sqrt{s}\text{leavevmode}{K}_1\left(\frac{z\sqrt{s}}{M_{\mathrm{sc}}}\right)ds,$$

$$Y_f^{\mathrm{eq}}=\frac{45\text{leavevmode}{g}_f}{4{\pi }^4}{\left(\frac{z{M}_f}{M_{\mathrm{sc}}}\right)}^2\frac{K_2\left(z{M}_f/M_{\mathrm{sc}}\right)}{g_s(M_{\mathrm{sc}}/z)},$$

Full text

1 Introduction

The standard model (SM) is a very successful theory in describing nature. The discovery of the last missing piece of the SM, viz., the Higgs boson, further increases its concreteness. In spite of its tremendous success, the SM can not explain a number of phenomena - two of the most important ones being the presence of dark matter (DM) and non-zero neutrino mass. Presence of DM in the universe is a very well established fact. The first indication of DM came from the observation of Galactic velocities within the Coma cluster by Fritz Zwicky in 1933 Zwicky:1933gu , followed by the observation of galaxy rotation curves by Vera Rubin in 1970 Rubin:1970zza . Subsequently, the observation of bullet cluster Clowe:2006eq firmly confirmed the presence of DM. Currently the best measurement of the amount of DM present in the universe comes from the Planck data Ade:2015xua ,

[TABLE]

where $h$ is the reduced Hubble parameter and of order unity. Unfortunately, the SM does not have any fundamental particle which can be a viable DM candidate. Therefore, to address the issue of DM from particle physics point of view, we need to extend the SM particle content and/or its gauge group. One of the most promising scenarios is to consider the DM candidate as a Weakly Interacting Massive Particle (WIMP) Gondolo:1990dk ; Srednicki:1988ce , which is produced in the early universe through the thermal freeze-out mechanism Gondolo:1990dk ; Srednicki:1988ce . However, WIMP type DM attracts stringent bounds from direct and indirect detection experiments Aprile:2018dbl ; Tan:2016zwf ; Akerib:2016vxi ; Amole:2017dex ; Aprile:2019dbj ; Ahnen:2016qkx ; Abdallah:2016ygi ; Giesen:2015ufa . In particular, a large portion of the parameter space in the spin independent/dependent WIMP-nucleon cross section and DM mass plane is ruled out by the direct detection (DD) bounds. Moreover, in near future with increasing sensitivity of the DD experiments Aprile:2018dbl ; Tan:2016zwf ; Akerib:2016vxi ; Amole:2017dex ; Aprile:2019dbj , these bounds might touch the so-called neutrino floor Drukier:1983gj ; Cushman:2013zza . In this work, we follow a non-thermal way of DM production, viz., via the freeze-in mechanism Hall:2009bx . In this scenario, the DM is very feebly interacting with the other particles, and as a result never achieves thermal equilibrium in the early universe with the cosmic soup. Hence it is named Feebly Interacting Massive Particles (FIMPs). Due to their very feeble interactions, FIMPs easily escape the above mentioned DD bounds while satisfying the measured value for the DM relic density Hall:2009bx ; Konig:2016dzg ; Biswas:2016bfo ; Biswas:2016iyh ; Biswas:2016yjr ; Biswas:2017tce ; Biswas:2017ait ; Kaneta:2016vkq ; Caputo:2018zky .

On the other hand, results of the neutrino oscillation experiments Cowan:1992xc ; Fukuda:1998mi ; Ahmad:2002jz ; Eguchi:2002dm ; An:2015nua ; RENO:2015ksa ; Abe:2014bwa ; Abe:2015awa ; Salzgeber:2015gua ; Adamson:2016tbq ; Adamson:2016xxw have confirmed oscillations between neutrino flavours. Since neutrino flavour oscillations are a clear proof of the neutrinos being massive and mixed, the neutrino oscillation experiments contradict the SM which postulates that the neutrinos are massless. Consequently, in order to explain tiny neutrino masses, one has to extend the SM by adding new particles and/or additional gauge groups.

In the present work we explain the above two puzzles by extending the SM gauge group by a $U(1)_{B-L}$ gauge symmetry as a simple (minimal) and well motivated extension of the SM, where $B$ is the baryon number and $L$ is the lepton number. In addition to the extra neutral gauge boson $Z^{\prime}$ associated with the $U(1)_{B-L}$, an extra SM singlet scalar $\phi_{H}$ (charged under $U(1)_{B-L}$ to break $B-L$ gauge symmetry spontaneously) is added in this simple extension, which leads to interesting signatures at the LHC Khalil:2006yi ; Emam:2007dy ; Huitu:2008gf ; Blanchet:2009bu ; Basso:2008iv . Moreover, nine additional SM singlet fermions ($N_{R}^{i}$ and $S^{i}_{1,2},\leavevmode\nobreak\ i=1,2,3$) are needed to explain the naturally111Here, “naturally” means the Dirac neutrino masses, $M_{D}$, have the same size as the Dirac masses of the SM fermions and, in contrary to the usual type-I seesaw mechanism Mohapatra:1980qe ; Marshak:1979fm ; Wetterich:1981bx ; Masiero:1982fi ; Mohapatra:1982xz ; Buchmuller:1991ce ; Abbas:2007ag , large Dirac neutrino Yukawa couplings, $\lambda_{d}\sim{\cal O}(0.1)$, with right-handed neutrino ($N_{R}^{i}$) masses are of order TeV. small neutrino masses through the inverse seesaw mechanism MohapatraIS1 ; Mohapatra:1986bd ; Khalil:2010iu ; GonzalezGarcia:1988rw . These additional fermions are not only required to generate the tiny neutrino masses via the inverse seesaw mechanism but are also needed for the gauge anomaly cancellation. In such a framework, three of these SM singlet fermions, $S_{1}^{i}$, are completely decoupled due to the introduction of $\mathbb{Z}_{2}$ symmetry and have naturally small mass (of order keV) according to ’t Hooft’s naturalness criterion tHooft:1979rat . Therefore, the lightest one, $S^{1}_{1}$, will be a stable particle and hence a warm DM (WDM) candidate Pagels:1981ke ; Peebles:1982ib ; Bond:1982uy ; Olive:1981ak , as discussed in Ref. El-Zant:2013nta . Moreover, since these keV mass singlet fermions are odd under the $\mathbb{Z}_{2}$ symmetry, they have no mixing with the active neutrinos and consequently are safe from the bound imposed by the x-ray observations Boyarsky:2009ix . In Ref. El-Zant:2013nta , an extra moduli field was introduced to produce this keV WDM non-thermally to achieve the correct ballpark value of relic density consistent with the WMAP and Planck observations. In the current work, without introducing any extra field contrary to Ref. El-Zant:2013nta , we successfully produce the keV WDM by the freeze-in mechanism through the decay and annihilation channels of $Z^{\prime}$. After explaining the keV FIMP WDM as a successful single component FIMP DM scenario to satisfy the correct value of the DM relic density, we study a two component FIMP DM as another possible scenario in the present model, where in addition to the FIMP WDM $S_{1}^{1}$, the lightest heavy right-handed neutrino $\nu_{H}^{1}$ can be a FIMP DM (with mass of order GeV) by tuning its corresponding Yukawa coupling to be very small Fiorentin:2016avj ; DiBari:2016guw . The GeV scale FIMP DM can be produced through the decay and annihilation processes of both the extra neutral gauge boson $Z^{\prime}$ as well as the extra $B-L$ Higgs $h^{\prime}$, while the keV FIMP WDM is produced only through the decay and annihilation processes of $Z^{\prime}$.

The rest of the paper is organized as follows. In section 2 we discuss the $B-L$ model with inverse seesaw mechanism and how the light neutrinos acquire their tiny masses. In section 3 we show that a keV sterile neutrino can be a WDM and produce the observed DM relic density as a single component FIMP DM. Section 4 is dedicated for studying two component FIMP type DM. Finally, our conclusions are given in section 5.

2 $B-L$ model with inverse seesaw scenario and neutrino masses

The gauged $B-L$ extension of the SM (BLSM) is based on the gauge group $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times U(1)_{B-L}$. By imposing $U(1)_{B-L}$, the gauge sector of the SM is extended to include a new neutral gauge boson $Z^{\prime}$ associated with the $B-L$ gauge symmetry. In addition, it has three SM singlet fermions $N_{R}^{i}$ (three right-handed neutrinos) with $B-L$ charge $=-1$ that arise as a result of anomaly cancellation conditions. Included also is an extra SM singlet scalar $\phi_{H}$ with $B-L$ charge $=-1$, while $\phi_{h}$ is the usual electroweak (EW) Higgs doublet. In order to satisfy the experimental measurements for the non-vanishing light neutrino masses with TeV scale right-handed (RH) neutrino using type-I seesaw mechanism, a very small Dirac neutrino Yukawa couplings, $\lambda_{d}\lesssim{\cal O}(10^{-6})$ must be assumed Mohapatra:1980qe ; Marshak:1979fm ; Wetterich:1981bx ; Masiero:1982fi ; Mohapatra:1982xz ; Buchmuller:1991ce ; Abbas:2007ag . Therefore, the mixing angle between the left- and right-handed neutrinos is quite suppressed, as it is proportional to $\lambda_{d}\lesssim{\cal O}(10^{-6})$. As a consequence of such small mixing angle, the interactions between the RH neutrinos and the SM particles are very suppressed, making it difficult to observe them at the LHC Khalil:2006yi ; Emam:2007dy ; Huitu:2008gf ; Blanchet:2009bu ; Basso:2008iv . Thus, we generate neutrino masses using the so-called inverse seesaw mechanism MohapatraIS1 ; Mohapatra:1986bd ; Khalil:2010iu ; GonzalezGarcia:1988rw that can naturally accommodate light neutrino masses with TeV scale RH neutrinos and large Yukawa couplings. In addition to the particle content as mentioned above, the BLSM with Inverse Seesaw (BLSMIS) has three extra pairs of SM singlet fermions $(S_{1,2}^{i},\leavevmode\nobreak\ i=1,2,3)$ with $B-L$ charge $=\mp 2$, respectively. In Table 1, we show the complete particle spectrum for the BLSMIS model with their associated charges for different gauge groups. An additional discrete symmetry has been introduced, viz., $\mathbb{Z}_{2}$. All BLSMIS particles are even under this symmetry except $S_{1}$ which is odd. Due to this symmetry, terms like $N^{c}_{R}\phi^{\dagger}_{H}S_{1}$ and $S_{1}S_{2}$, that could spoil the usual inverse seesaw mechanism, are forbidden MohapatraIS1 ; Mohapatra:1986bd ; Khalil:2010iu ; GonzalezGarcia:1988rw . The complete Lagrangian for this model is given by

[TABLE]

where $F^{\prime}_{\mu\nu}=\partial_{\mu}B^{\prime}_{\nu}-\partial_{\nu}B^{\prime}_{\mu}$ is the $U(1)_{B-L}$ field strength, $D_{\mu}$ is the covariant derivative, $\tilde{\phi_{h}}=i\sigma_{2}\phi_{h}$ and the flavor indices are omitted for simplicity. The general structure of the covariant derivative $D_{\mu}$ in the present model takes the following form

[TABLE]

where $(g_{c},T^{\alpha},G^{\alpha}_{\mu})$ are the $SU(3)_{C}$ gauge coupling, generator and the gauge field, respectively. Similarly, $(g,\tau^{a},W^{a}_{\mu})$, $(g_{Y},Y,B_{\mu})$ and $(g^{\prime},Y_{BL},B^{\prime}_{\mu})$ are the corresponding quantities for $SU(2)_{L}$, $U(1)_{Y}$ and $U(1)_{B-L}$, respectively. It is worth mentioning that a kinetic mixing term $F^{\prime}_{\mu\nu}F^{\mu\nu}$ is allowed and it leads to a non-vanishing $Z$-$Z^{\prime}$ mixing angle, $\theta^{\prime}$ Holdom:1985ag ; Chankowski:2006jk ; Krauss:2012ku . However, due to the stringent constraint from LEP experiments on the $Z$-$Z^{\prime}$ mixing angle ($|\theta^{\prime}|\lesssim 10^{-3}$) Abreu:1994ria ; Alcaraz:2006mx ; Erler:2009jh , one may neglect this term. Finally, the potential $\mathcal{V}(\phi_{h},\phi_{H})$ is given by Emam:2007dy ; Khalil:2006yi

[TABLE]

where the potential $\mathcal{V}(\phi_{h},\phi_{H})$ will be bounded from below when the following inequalities are satisfied simultaneously

[TABLE]

Here, both the scalars $\phi_{H}$ and $\phi_{h}$ acquire their non-zero vacuum expectation values (VEVs), therefore, the $B-L$ and the EW symmetries are broken spontaneously and the SM Higgs doublet $\phi_{h}$ and the $B-L$ singlet $\phi_{H}$ take the following form:

[TABLE]

where $v\simeq$ 246 GeV is the EW symmetry breaking scale and $v^{\prime}$ is the scale of $B-L$ symmetry breaking which is, in general, unknown and ranging from TeV to much higher scales. After breaking the $B-L$ and the EW symmetries spontaneously, the extra neutral gauge boson $Z^{\prime}$ acquires its mass $(M_{Z^{\prime}}=g^{\prime}v^{\prime})$ Emam:2007dy ; Khalil:2006yi 222The experimental search for $Z^{\prime}$, by LEP II Cacciapaglia:2006pk ; Carena:2004xs , leads to another constraint: $M_{Z^{\prime}}/g^{\prime}\gtrsim 7$ TeV. This constraint will easily be satisfied due to a smallness of $g^{\prime}$ which is required by the freeze-in scenario Hall:2009bx ., and the neutrino Yukawa interaction terms in Eq. (2) and in addition a very small Majorana mass $\mu_{S}$ for $S_{1,2}$ lead to the following neutrino mass terms333$\mu_{S}$ is naturally small due to ’t Hooft’s naturalness criterion tHooft:1979rat , for simplicity we assume $S_{1}$ and $S_{2}$ have the same small Majorana mass ($\mu_{S}$), and the generation of such small $\mu_{S}$ from non-renormalizable terms has been discussed in Abdallah:2011ew and radiatively in Ma:2009gu .

[TABLE]

where $M_{D}=\lambda_{d}v/\sqrt{2}$ and $M_{N}=\lambda_{s}v^{\prime}/\sqrt{2}$. Therefore, the neutrino mass matrix in the basis $(\nu^{c}_{L},N_{R},S_{2},S_{1})$ can be written as

[TABLE]

It is clearly seen that $S_{1}$ is completely decoupled and has no mixing with active neutrinos. It only interacts with the neutral gauge boson $Z^{\prime}$ with a coupling $g_{Z^{\prime}S_{1}S_{1}}=g^{\prime}$. Therefore, $S_{1}$ is free from cosmological and astrophysical constraints coming from active-sterile mixing Boyarsky:2009ix . Thus its mass is given as,

[TABLE]

After diagonalising the upper left $3\times 3$ submatrix of the neutrino mass matrix $\mathcal{M}_{\nu}$, the light and heavy neutrino masses, respectively, are given by

[TABLE]

where $\mu_{S}\ll M_{D},M_{N}$ is assumed. One can naturally obtain eV scale light neutrino masses with $\mu_{S}$ of order keV and $M_{N}$ of order TeV, keeping Yukawa coupling $\lambda_{d}$ of order one. Such large couplings between the heavy RH neutrinos and the SM particles leads to interesting implications and enhances the accessibility of TeV scale $B-L$ at the LHC Abdelalim:2014cxa ; Bandyopadhyay:2012px ; Abdallah:2012nm .

Recall that due to the added $\mathbb{Z}_{2}$ symmetry, $S_{1}$ is completely decoupled. Hence the lightest fermionic singlet, $S_{1}^{1}$, is a stable particle and hence a DM candidate. Since its mass ($=\mu_{S}$) is of order keV, hence $S_{1}^{1}$ is a warm DM (WDM) candidate El-Zant:2013nta 444The contribution of the new light degrees of freedom ($S_{1}^{i}$) to the number of effective neutrino species, $N_{\rm eff}$, has been checked using Eq. (5) in Ref. Hooper:2011aj to calculate extra effective neutrino species, $\Delta N_{\rm eff}$, and found it to be negligible.. Moreover, one can easily make the lightest heavy RH neutrino, $\nu_{H}^{1}$, stable or long-lived by taking the corresponding Yukawa coupling to be very small $\lesssim 3\times 10^{-26}({\rm GeV}/M_{N})^{1/2}$ Fiorentin:2016avj ; DiBari:2016guw . Thus, from here onwards we focus on the two component DM scenario, where, one of them is GeV scale DM, $\nu_{H}^{1}$, and the other is keV scale WDM, $S_{1}^{1}$.

It is important to note that due to the mixing term in the potential $\mathcal{V}(\phi_{h},\phi_{H})$, the squared mass matrix of the neutral Higgs bosons in the basis $(h,H)$ is non-diagonal and takes the following form:

[TABLE]

Rotating this matrix into the basis $(h_{1},h_{2})$ which is defined as follows

[TABLE]

where the mixing angle $\alpha$ takes the following form:

[TABLE]

Therefore, the masses of these two physical Higgs scalars $(h_{1},h_{2})$ are given by555Hereafter, the physical state $h_{1}$ refers to the SM-like Higgs boson and its mass $M_{h_{1}}$ is fixed at $125.5$ GeV to agree with the LHC measurements Aad:2012tfa ; Chatrchyan:2012xdj . Also, according to the measured values of Higgs boson signal strengths for its various decay modes, the mixing angle $\alpha$ should be very small, thus we have fixed it at $0.01$ rad.

[TABLE]

The quartic couplings $\lambda$’s can be written in terms of the physical masses $M_{h_{1,2}}$ as follows Basso:2010jm

[TABLE]

We have used SARAH Staub:2013tta ; Staub:2009bi ; Staub:2008uz to implement the BLSMIS and the relevant masses, couplings and decay widths have been calculated using SPheno Porod:2011nf .

3 Warm DM as FIMP

As mentioned earlier, $S^{1}_{1}$ is a WDM candidate with mass in the few keV range Adhikari:2016bei ; Bernal:2017kxu ; Abada:2014zra . We next study in detail the production of this keV DM via the freeze-in mechanism. Here $S_{1}^{1}$ is produced solely from its coupling with the extra $U(1)_{B-L}$ gauge boson $Z^{\prime}$, as mentioned in the previous section. The corresponding gauge coupling $g^{\prime}$ is taken to be extremely feeble $\sim\mathcal{O}(10^{-10})$ with the result that $S_{1}^{1}$ is never in thermal equilibrium with the cosmic soup. Due to this small $B-L$ gauge coupling, the corresponding gauge boson $Z^{\prime}$ also interacts very feebly with the cosmic soup and never attains thermal equilibrium Arcadi:2013aba ,

[TABLE]

where $\Gamma_{Z^{\prime}}$ is the total decay width of $Z^{\prime}$ and $H$ is the Hubble parameter. Therefore, we first determine the distribution function for $Z^{\prime}$ 666As $Z^{\prime}$ is not in thermal equilibrium (due to very small value of $g^{\prime}$), one can not assume a Maxwell-Boltzmann distribution function for $Z^{\prime}$. Therefore, the $Z^{\prime}$ distribution can be found by solving Eq. (24).. The general formalism to determine the distribution function of any particle (say $f$) is to solve the following Boltzmann equation:

[TABLE]

where $\hat{L}$ is the Lioville’s operator and $\mathcal{C}[f]$ is known as the collision term of $f$. If we consider an isotropic and homogeneous universe, then, using the Friedman-Robertson-Walker metric, the Lioville’s operator takes the following form:

[TABLE]

where $p$ is the absolute value of the particle’s three momentum, $|\vec{p}|$. Following Konig:2016dzg , we perform a transformation of variables, $(p,\,t)\to(\xi_{p},\,z)$, in the following way:

[TABLE]

where $g_{s}(T)$ is the effective entropy degrees of freedom (d.o.f) at temperature $T$, $M_{\rm sc}$ is an arbitrary mass scale and hereafter we take it equal to the SM-like Higgs mass ($M_{\rm sc}=M_{h_{1}}=125.5$ GeV) and $T_{0}$ is the initial temperature at which the DM relic density is taken to be zero. Therefore, using the following time-temperature relation,

[TABLE]

the Lioville’s operator defined in Eq. (25) can be simply written as

[TABLE]

where $g^{\prime}_{s}(T)$ is the derivative of $g_{s}(T)$ with respect to the temperature $T$.

Taking only the decay term for the $Z^{\prime}$ production777In principle, the collision term for annihilation diagrams should also be considered but in this class of models those annihilation diagrams have subleading contribution Biswas:2016bfo , hence we have not taken into account those effects and for simplicity we consider only the decay of $h_{1,2}$ as the $Z^{\prime}$ production mechanism. Moreover, $h_{1,2}$ are in thermal equilibrium, and consequently the usual equilibrium Boltzmann distribution function has been assumed for them Konig:2016dzg ., the Boltzmann equation of the distribution function of $Z^{\prime}$ is given by

[TABLE]

where $f_{Z^{\prime}}$ is the distribution function of $Z^{\prime}$, $\mathcal{C}^{h_{i}\rightarrow Z^{\prime}Z^{\prime}}$ is the collision term of $Z^{\prime}$ production from the decays of scalars $h_{1,2}$ and $\mathcal{C}^{Z^{\prime}\rightarrow\,{\rm all}}$ is $Z^{\prime}$ decay collision term due to all its possible decay channels. The expression of these collision terms are given in the Appendix A. Once we get the distribution function of $Z^{\prime}$ by solving Eq. (29), we then can determine its co-moving number density by using the following relation:

[TABLE]

where $n_{Z^{\prime}}$ is the $Z^{\prime}$ number density, $g$ is the internal d.o.f of $Z^{\prime}$ and the universe entropy density $s$ is given by $s=\leavevmode\nobreak\ (2\pi^{2}/45)T^{3}g_{s}(T)$ Kolb:1990vq .

From Eq. (30), one can note that the co-moving number density of $Z^{\prime}$ is directly proportional to the integrated $\xi_{p}^{2}f_{Z^{\prime}}$, i.e., larger the area under a $\xi_{p}^{2}f_{Z^{\prime}}$ curve, larger is the $Z^{\prime}$ abundance. In Fig. 1, we show the variation of $\xi_{p}^{2}f_{Z^{\prime}}$ with respect to the dimensionless parameter $\xi_{p}$ for different values of $z$ ($=M_{\rm sc}/T$). As shown in the figure, areas under the curves corresponding to $z=0.02$ and $20.0$ are different because for higher $z=20.0$ (i.e., lower temperature $T$ of the universe), $Z^{\prime}$ gets more time to be produced and it then subsequently decays into WDM and the SM fermions. But as $z$ is increased further (presently $z=100.0$), $Z^{\prime}$ starts decaying significantly and its abundance gets depleted and the area under the curve for $z=100$ is smaller than for $z=20.0$, as seen in Fig. 1. For still higher values of $z$ ($z=500.0$), $Z^{\prime}$ abundance decreases further due to decay. Thus, as $z\rightarrow\infty$, $Z^{\prime}$ will gradually decay to DM and its abundance eventually goes to zero.

Once the distribution function of $Z^{\prime}$ is computed, we can describe the production of the keV DM $S^{1}_{1}$. In the present scenario, the keV DM $S_{1}^{1}$ can be produced from the decay of $Z^{\prime}$, $Z^{\prime}\to S_{1}^{1}S_{1}^{1}$ (decay contribution), and from the annihilation processes, $f\bar{f}\to S_{1}^{1}S_{1}^{1}$ mediated by $Z^{\prime}$, where $f=l,q,\nu_{l}$. The annihilation contribution has been calculated by using micrOMEGAs Belanger:2018mqt . To determine the co-moving number density of the WDM $S_{1}^{1}$, we solve the following Boltzmann equation,

[TABLE]

where $M_{\rm Pl}=1.22\times 10^{19}$ GeV is the Planck mass, $\sqrt{g_{\star}}=\frac{g_{s}(z)}{\sqrt{g_{\rho}(z)}}\left(1-\frac{1}{3}\frac{{\rm d}\ln g_{s}(z)}{{\rm d}\ln z}\right)$, where $g_{\rho}$ is the effective energy degrees of freedom Gondolo:1990dk ; Edsjo:1997bg , and the non-thermal average of $Z^{\prime}$ decay width is defined by

[TABLE]

where

[TABLE]

The expressions of a thermal average annihilation cross section $\langle{\sigma v}_{f\bar{f}\rightarrow S_{1}^{1}S_{1}^{1}}\rangle$ and an equilibrium co-moving number density of $f$ ($Y^{\rm eq}_{f}$), appearing in Eq. (31), are given respectively by Gondolo:1990dk ; Edsjo:1997bg

[TABLE]

[TABLE]

where $g_{f}$ is the internal d.o.f of $f$, $K_{1}(z)$ and $K_{2}(z)$ are the Bessel function for first and second kind, respectively, and ${\sigma}_{f\bar{f}\rightarrow S_{1}^{1}S_{1}^{1}}$ is given in Ref. Iso:2010mv . Solving the Boltzmann equation given by Eq. (31) gives us the co-moving number density $Y_{S_{1}^{1}}$. The corresponding relic density of the WDM $S_{1}^{1}$ can be calculated by using the following formula Edsjo:1997bg ,

[TABLE]

In order to understand the relative contribution of the decay and annihilation channels we will first consider them one at a time and solve the Boltzmann equation to get the relic density. We start with taking only the decay contribution and show in the left and right panels of Fig. 2 the variation of DM relic density as a function of $z$, for different values of the initial temperature $T_{0}$ ($=M_{\rm sc}/z_{0}$) and different values of the WDM mass $M_{S_{1}^{1}}$, respectively. The horizontal magenta dashed line refers to the DM relic density measurement ($\Omega h^{2}\simeq 0.12$) Ade:2015xua . In the left panel, as long as $T_{0}$ is greater than the mass of the mother particles ($h_{1,2}$) in $h_{1,2}\to Z^{\prime}Z^{\prime}$ decay channels, the final DM relic density is insensitive to $T_{0}$, as seen for the ${{z_{0}=0.001,0.005,0.01}}$ cases. However, once $T_{0}$ drops below the mass of the mother particles (presently ${{z_{0}=2}}$ case), $Z^{\prime}$ production gets Boltzmann suppressed and consequently DM relic density is suppressed. In the right panel we show the dependence of the DM relic density on its mass ($M_{S_{1}^{1}}$) for $z_{0}=0.01$. It is clear that the relic density increases with the DM mass, as expected from Eq. (36).

For the annihilation contribution, $f\bar{f}\to S_{1}^{1}S_{1}^{1}$, there are two relevant regimes are as follows. $(i)$ The on-shell regime, where $2M_{S^{1}_{1}}<M_{Z^{\prime}}<T_{0}$, in which $Y_{S_{1}^{1}}(\propto g^{\prime 4}/\Gamma_{Z^{\prime}})$ does not depend on $T_{0}$ and $(ii)$ The EFT regime, where $2M_{S^{1}_{1}}<T_{0}<M_{Z^{\prime}}$, in which $Y_{S_{1}^{1}}(\propto g^{\prime 4}T^{3}_{0}/M_{Z^{\prime}}^{4})$ depends on $T_{0}$. In the left panel of Fig. 3 we see that as long as $T_{0}$ is greater than $M_{Z^{\prime}}$, the final DM relic density is insensitive to $T_{0}$ (on-shell regime), as seen for the $z_{0}=0.001,0.005,0.01,2.0$ cases. Once $T_{0}$ drops below $M_{Z^{\prime}}$ (presently ${{z_{0}=100.0}}$ case), then $S_{1}^{1}$ production gets the suppressed by a factor $T_{0}^{3}/M_{Z^{\prime}}^{4}$ (EFT regime). In the right panel of Fig. 3, we show the dependence of the DM relic density on its mass ($M_{S_{1}^{1}}$) for $z_{0}=0.01$ (on-shell regime). It is clear that the relic density increases with the DM mass, as expected from Eq. (36).

For the decay contribution ($Z^{\prime}\to S_{1}^{1}S_{1}^{1}$), we show in Fig. 4 the variation of the co-moving number density of $Z^{\prime}$ (left panel) and the co-moving number density of $S_{1}^{1}$ with $z$, for different values of $B-L$ gauge coupling $g^{\prime}$ (left panel) and $M_{Z^{\prime}}$ (right panel). Since $Z^{\prime}$ production is proportional to the $B-L$ gauge coupling $g^{\prime}$, larger $g^{\prime}$ results in larger $Z^{\prime}$ production and consequently a larger production of DM, as seen in the left panel of this figure. Note also that $Z^{\prime}\to S_{1}^{1}S_{1}^{1}$ decay rate is directly proportional to $g^{\prime}$, and hence increasing $g^{\prime}$ increases the decay rate of $Z^{\prime}$ and hence the abundance of $S_{1}^{1}$. Therefore, it is clear that for higher values of $g^{\prime}$ the $Z^{\prime}$ co-moving number density plateau starts bending at smaller values of $z$. On the other hand, in the right panel of Fig. 4, we see that by increasing $M_{Z^{\prime}}$ exactly opposite behavior appears for $Z^{\prime}$ production while similar behavior happens for its decay. As mentioned, $Z^{\prime}$ production mainly happens through decay of the Higgs bosons ($h_{1,2}$) and those decay modes ($\Gamma_{h_{i}\to Z^{\prime}Z^{\prime}}$) are proportional to $M^{3}_{h_{i}}/M_{Z^{\prime}}^{2}$ [see Eq. (47)]. Therefore, increasing $M_{Z^{\prime}}$ reduces the production rate of $Z^{\prime}$ as $1/M_{Z^{\prime}}^{2}$. However, its decay width is simply proportional to its mass $M_{Z^{\prime}}$ [see Eq. (• ‣ A)] and so increasing $M_{Z^{\prime}}$ results in faster decay of $Z^{\prime}$. In Fig. 5 we show similar plots for the annihilation contribution. This figure shows features similar to Fig. 4.

In Fig. 6, we show the total relic density (blue dashed-dotted line) as well as the relative contributions of the two different types of WDM production processes, decay (red dashed line) and annihilation (green dotted line). Here, for a suitably selected set of model parameters ($M_{S_{1}^{1}}=10$ keV, $M_{Z^{\prime}}=10$ GeV, $g^{\prime}=1.0\times 10^{-9}$, $M_{h_{2}}$ = 1 TeV, $z_{0}=0.01$, and $\alpha=0.01$ rad), the total WDM relic density equals the observed relic density ($\Omega h^{2}\simeq 0.12$) at the present epoch, where decay contributes $\sim 62\%$ of the WDM relic density while the rest comes from annihilation. It is worth mentioning that initially for $z\lesssim 100$, WDM is dominantly produced from the annihilation processes this is because of all ingoing particles are already in the cosmic soup, while for $z\gtrsim 100$, the decay process starts dominating, as seen in Fig. 4.

Variation of total WDM relic density ($\Omega h^{2}$) as a function of the gauge coupling $g^{\prime}$ can be seen in Fig. 7, where the BLSMIS points have been generated over the following ranges of its fundamental parameters: $1\leq M_{S_{1}^{1}}\leq 10$ keV, $1\leq M_{Z^{\prime}}\leq 100$ GeV, $10^{-12}\leq g^{\prime}\leq 10^{-8}$, $M_{h_{2}}=1$ TeV, $z_{0}=0.01$, and $\alpha=0.01$ rad. From the left panel it is clear that $\Omega h^{2}$ is inversely proportional to $M_{Z^{\prime}}$ (which is represented by the color bar). More explicitly, for a fixed $g^{\prime}$ value, larger $\Omega h^{2}$ values correspond to smaller $M_{Z^{\prime}}$ values (red points) and vice versa for the blue points. On the other hand as illustrated in the right panel, $\Omega h^{2}$ is directly proportional to the WDM mass $M_{S^{1}_{1}}$ (which is represented by the color bar). This is consistent with $\Omega h^{2}$ expression given in Eq. (36).

In Fig. 8 we show the allowed points in the $(M_{Z^{\prime}},g^{\prime})$ and $(M_{S_{1}^{1}},g^{\prime})$ planes in the left and right panels, respectively, which give the relic density consistent with a relic density upper bound of the Planck measurement ($\Omega h^{2}\leq 0.12$) Ade:2015xua . All other parameter values are allowed to vary in the range mentioned in the previous paragraph. From the figure color bars (mapped to the total WDM relic density $\Omega h^{2}$), it is clearly seen that many points ($\sim 84\%$ of the scanned points) have a small DM relic density ($\Omega h^{2}\leq 10^{-2}$). Therefore, in the next section we discuss a two component FIMP DM possibility as a well-motivated scenario to get an extra relic density contribution from the lightest heavy RH neutrino, $\nu_{H}^{1}$, as a GeV scale DM.

4 Two component FIMP dark matter

In the previous section we have studied the WDM FIMP $S_{1}^{1}$, as a single component DM. As mentioned in section 2, the lightest heavy RH neutrino, $\nu_{H}^{1}$, can be long-lived particle by making the corresponding Yukawa couplings very small $\lesssim 3\times 10^{-26}({\rm GeV}/M_{N})^{1/2}$ Fiorentin:2016avj ; DiBari:2016guw . Therefore, it can be an additional DM component of mass of order GeV. Note that any interaction between $S_{1}^{1}$ and $\nu_{H}^{1}$ is completely forbidden. Thus in the present section, we consider a two component DM scenario with two DM candidates: the WDM FIMP $S_{1}^{1}$ and the lightest heavy RH neutrino $\nu_{H}^{1}$. The dominant annihilation channels of $\nu^{1}_{H}$ pair to SM particles are mediated by $Z^{\prime}$ and $h_{i}$ ($i=1,2$).888Due to the smallness of the corresponding Yukawa coupling of $\nu_{H}^{1}$ (as assumed to be a stable DM candidate), the contribution of the channels mediated by $Z$ and $W^{\pm}$ bosons is negligible. Also, the annihilation channels mediated by the SM-like Higgs $h_{1}$ are suppressed as compared to the $h_{2}$ ones, because the coupling $\lambda_{\nu^{1}_{H}\nu^{1}_{H}h_{1}}$ is very small since it is proportional to $\sin\alpha$ which is constrained to be very small by LHC Khachatryan:2016vau . The coupling strength of $\nu^{1}_{H}$ pair with $Z^{\prime}$ is given by $g^{\prime}/2$, while with $h_{i}$ ($i=1,2$) is given as

[TABLE]

where $O_{1}=\sin\alpha$ and $O_{2}=\cos\alpha$. Therefore, $\nu_{H}^{1}$ pair annihilation is proportional to the gauge coupling $g^{\prime}$ which is taken very small in the present model. Due to this feeble gauge coupling $g^{\prime}$, $\nu_{H}^{1}$ will never reach thermal equilibrium and is produced by the freeze-in mechanism. The Boltzmann equation associated with $\nu_{H}^{1}$ production is as follows

[TABLE]

where $\langle\Gamma_{Z^{\prime}\rightarrow\nu^{1}_{H}\nu^{1}_{H}}\rangle_{\rm NTH}$ and $\langle{\sigma v}_{f\bar{f}\rightarrow\nu^{1}_{H}\nu^{1}_{H}}\rangle$ are defined as in Eqs. (32),(34), respectively, by replacing $S_{1}^{1}$ with $\nu_{H}^{1}$, while $Y^{\rm eq}_{h_{i}}$ is defined as in Eq. (35) by replacing $f$ with $h_{i}$. Thermal average of the decay width of $h_{i}$ ($i=1,2$) is defined as Gondolo:1990dk

[TABLE]

where $\Gamma_{h_{i}}$ is the total decay width of $h_{i}$. After solving the Boltzmann equation of $\nu_{H}^{1}$ production, Eq. (38), the corresponding relic density of $\nu^{1}_{H}$ can be determined by using the following relation,

[TABLE]

Finally, the total relic density of this two component DM scenario is given by

[TABLE]

where $\Omega_{S_{1}^{1}}h^{2}$ is the relic density of $S_{1}^{1}$ which is defined in Eq. (36).

It is clear that the DM production depends crucially on the DM mass and the mass of the mother particles ($M_{Z^{\prime}},M_{h_{2}}$). Assuming $M_{h_{2}}>2M_{Z^{\prime}}>4M_{S_{1}^{1}}$, we divide the $\nu_{H}^{1}$ DM spectrum into two regions according to the dominant production modes of $\nu_{H}^{1}$ DM - Region I, where $M_{Z^{\prime}}>2M_{\nu_{H}^{1}}$ and $\nu_{H}^{1}$ production is $Z^{\prime}$ dominated, and Region II, where $M_{Z^{\prime}}<2M_{\nu_{H}^{1}}$ and $\nu_{H}^{1}$ production is $h_{2}$ dominated.

4.1 Region I: $M_{Z^{\prime}}>2M_{\nu_{H}^{1}}$

For our chosen set of BLSMIS parameters ($M_{Z^{\prime}}=1$ TeV, $M_{\nu_{H}^{1}}=70$ GeV, $M_{S_{1}^{1}}=10$ keV, $g^{\prime}=9.0\times 10^{-12}$, $M_{h_{2}}$ = 5 TeV, $\alpha=0.01$ rad, and $z_{0}=0.01$), we show in Fig. 9 the variation of the total DM relic density (blue solid line) and its relative contributions for a two component DM scenario. In the figure, red solid and green dashed lines correspond to the $\nu_{H}^{1}$ relic density contributions from the decay of $h_{2}$ and $Z^{\prime}$, respectively999Due to a smallness of the mixing angle $\alpha$, DM production of the SM-like Higgs $h_{1}$ is negligible., while cyan dashed line corresponds to the annihilation contribution ($\Omega^{\rm ann}_{\nu_{H}^{1}}h^{2}$). In addition, the $S_{1}^{1}$ relic density contribution from the decay of $Z^{\prime}$ ($\Omega^{Z^{\prime}}_{S_{1}^{1}}h^{2}$) and annihilation ($\Omega^{\rm ann}_{S_{1}^{1}}h^{2}$) are presented by green dashed and cyan solid lines, respectively. Note that in region I, the relative contribution of $Z^{\prime}$ decay to $\nu_{H}^{1}$ production ($\Omega^{Z^{\prime}}_{\nu_{H}^{1}}h^{2}$) is larger than the $h_{2}$ decay contribution ($\Omega^{h_{2}}_{\nu_{H}^{1}}h^{2}$) because the latter is suppressed by a factor of their partial decays ratio ($\Gamma_{h_{2}\to\nu_{H}^{1}\nu_{H}^{1}}/\Gamma_{Z^{\prime}\to\nu_{H}^{1}\nu_{H}^{1}}$ $\simeq 12M^{2}_{\nu_{H}^{1}}M_{h_{2}}/M^{3}_{Z^{\prime}}\simeq{\cal O}(0.1)$). It is also worth noting that the relic density contribution of the keV DM ($S_{1}^{1}$) is negligible compared to the GeV DM ($\Omega_{\nu_{H}^{1}}h^{2}$) even though they have the same gauge coupling strength ($g^{\prime}$) and their mediator masses ($M_{h_{2}}$ and $M_{Z^{\prime}}$) are of the same order ($\sim$ TeV). This is simply because the relic density of a DM candidate is directly proportional to its mass [see Eqs. (36),(40)]. Therefore, the contribution of the keV mass $S_{1}^{1}$ to the DM total relic density is suppressed by a factor $\simeq M_{S_{1}^{1}}/M_{\nu_{H}^{1}}\simeq{\cal O}(10^{-7})$ as compared to the GeV mass $\nu_{H}^{1}$.

In the left and right panels of Fig. 10, we show the variation of relic density contributions of the two component DM scenario for different values of $M_{Z^{\prime}}$ and $g^{\prime}$, respectively. The top panels stand for the decay contribution while the bottom ones stand for the annihilation contribution. Again from these figures, one can easily conclude that FIMP relic density contributions are inversely proportional to the mediator mass, as illustrated in left panels, and directly proportional to coupling strength as shown in right panels. We have discussed these features before in section 3.

4.2 Region II: $M_{Z^{\prime}}<2M_{\nu_{H}^{1}}$

In discussed above, for $M_{Z^{\prime}}>2M_{\nu_{H}^{1}}$ (region I), the total relic density is dominated by $Z^{\prime}$ mediated diagrams. Now we turn to the region II, where $M_{Z^{\prime}}<2M_{\nu_{H}^{1}}$ and $Z^{\prime}$ decays to $\nu_{H}^{1}$ pair is kinematically forbidden, and consequently, $\nu_{H}^{1}$ production is $h_{2}$ dominated. Therefore in region II, a major portion of our two DM candidates, $\nu_{H}^{1}$ and $S_{1}^{1}$, is produced almost independently from the $h_{2}$ and $Z^{\prime}$ mediated processes, respectively. In other words, by fixing $M_{h_{2}}$, $g^{\prime}$ and $M_{\nu^{1}_{H}}/M_{Z^{\prime}}$ at certain values to get a significant contribution from $\Omega_{\nu_{H}^{1}}h^{2}$, one can obtain a relevant $\Omega_{S_{1}^{1}}h^{2}$ contribution independently by changing $M_{S^{1}_{1}}$ within keV range. This possibility did not exist in region I because both $\nu_{H}^{1}$ and $S_{1}^{1}$ are produced dominantly via $Z^{\prime}$ and therefore have the same number density. The only way to have comparable contribution from both in region I would be to raise the mass of $S_{1}^{1}$ to the GeV range. However, this is untenable since that will spoil the inverse seesaw mechanism scenario for generating light neutrino masses MohapatraIS1 ; Mohapatra:1986bd ; Khalil:2010iu ; GonzalezGarcia:1988rw . In region II this lacuna is remedied since here $\nu_{H}^{1}$ and $S_{1}^{1}$ are produced independently - while $\nu_{H}^{1}$ are dominantly produced from $h_{2}$, $S_{1}^{1}$ are produced from $Z^{\prime}$.

In Fig. 11, for two suitably chosen sets of BLSMIS parameters ($M_{Z^{\prime}}=10$ GeV ($2.5$ GeV), $M_{\nu_{H}^{1}}=8$ GeV ($2$ GeV), $M_{S_{1}^{1}}=10$ keV ($100$ keV), $\alpha=0.01$ rad, $g^{\prime}\simeq 2.4\times 10^{-11}$, $M_{h_{2}}=5$ TeV, and $z_{0}=0.01$), we show variation of decay and annihilation contributions to the relic density of $S_{1}^{1}$ and $\nu_{H}^{1}$, as a function of $z$. In the figure, green dashed-dotted and red dashed-double-dotted lines correspond to the $S_{1}^{1}$ relic density contributions ($\Omega^{Z^{\prime}}_{S_{1}^{1}}h^{2}$ and $\Omega^{\rm ann}_{S_{1}^{1}}h^{2}$), respectively, while blue dashed line corresponds to $\nu_{H}^{1}$ relic density contribution from decay of $h_{2}$ ($\Omega^{h_{2}}_{\nu_{H}^{1}}h^{2}$). From the figure, it is clearly seen that $S_{1}^{1}$ has a relevant relic density contribution, unlike the situation in region I. Note that for a larger $S_{1}^{1}$ mass ($M_{S_{1}^{1}}=100$ keV), the $S_{1}^{1}$ contribution to the total relic density even starts to be the dominant one, as seen in the right panel of Fig. 11.

5 Conclusion

In this work we studied two beyond SM problems, viz., the non-zero neutrino masses and the existence of the DM. In studying the tiny neutrino masses, we followed the inverse seesaw mechanism within the $B-L$ extension of the SM (BLSMIS). Six SM singlet fermions were introduced for inverse seesaw mechanism to work and three more singlet fermions (with mass of order keV) were added to cancel the $U(1)_{B-L}$ gauge anomaly. The lightest of these additional fermionic states, $S_{1}^{1}$, can be a WDM, being odd under a $\mathbb{Z}_{2}$ discrete symmetry. We studied $S_{1}^{1}$ as a FIMP WDM and showed that it could be produced via the freeze-in mechanism from the decay of the extra neutral gauge boson $Z^{\prime}$ and the on-shell annihilation processes mediated by $Z^{\prime}$. We showed that the relative contributions to the DM relic density from both the decay and the on-shell annihilation processes are more or less equal. We scanned over the relevant BLSMIS parameters by imposing the Planck constraint of the DM relic density and showed that a large portion of the parameter space gives a small contribution to the DM relic density. Therefore, we studied a two component FIMP DM as a possible scenario in the BLSMIS to get an extra contribution to the DM relic density. In this scenario, the lightest heavy RH neutrino, $\nu^{1}_{H}$, can contribute to the DM relic density as an independent DM component (with mass of order GeV). For $M_{Z^{\prime}}>2M_{\nu_{H}^{1}}$, we showed that the production of $\nu_{H}^{1}$ as a DM candidate through the $Z^{\prime}$ mediator has the dominant contribution to the total DM relic density. On the other hand for $M_{Z^{\prime}}<2M_{\nu^{1}_{H}}$, the $h_{2}$ mediated processes will contribute dominantly to $\nu_{H}^{1}$ production while $Z^{\prime}$ mediated processes will contribute dominantly to $S_{1}^{1}$ production. We emphasized that in this region both FIMP candidates ($S_{1}^{1}$ and $\nu_{H}^{1}$) can contribute to the total DM relic density.

Acknowledgements

The authors would like to thank the Department of Atomic Energy Neutrino Project of Harish-Chandra Research Institute (HRI). We also acknowledge the HRI cluster computing facility (http://www.hri.res.in/cluster/). This project has received funding from the European Union’s Horizon 2020 research and innovation programme InvisiblesPlus RISE under the Marie Sklodowska-Curie grant agreement No. 690575 and No. 674896. The authors would like to thank Abhass Kumar for fruitful discussions.

Appendix

Appendix A Analytical expression of the collision terms

For any generic process $A(\tilde{p})\rightarrow B(\tilde{p_{1}})\,C(\tilde{p_{2}})$ (where $\tilde{p}=(E_{p},\bar{p})$), the collision term takes the following form Kolb:1990vq ; Gondolo:1990dk ,

[TABLE]

From Eq. (29), we can see that the Boltzmann equation which determine the distribution function of the extra gauge boson $Z^{\prime}$ contains two collision terms one is for its production and the another one is for its decay. The expression for the two collision terms are described below.

• •

$\mathcal{C}^{Z^{\prime}\rightarrow{\rm all}}$:    Collision term for the $B-L$ gauge boson decay takes the following form after using the Eq. (42),

[TABLE]

where $r_{Z^{\prime}}=M_{Z^{\prime}}/T$ and $\Gamma_{Z^{\prime}\rightarrow{\rm all}}=\Gamma_{Z^{\prime}\rightarrow f\bar{f}}$

• $\Gamma_{Z^{\prime}\rightarrow\chi\chi}$. Expression for the each decay terms are as follows,

[TABLE]

where $f$ refers to the SM fermions and $\chi=S_{1}^{1},\nu_{H}^{1}$. $n_{c}$ and $q_{f}$ are the corresponding color and electric charges, respectively. $g_{Z^{\prime}\chi\chi}=g^{\prime}$ for $S_{1}^{1}$ and $g^{\prime}/2$ for $\nu_{H}^{1}$.

• •

$\mathcal{C}^{h_{i}\rightarrow Z^{\prime}Z^{\prime}}$:   The expression for this collision term where $h_{i}$ ($i=1,2$) decays to $Z^{\prime}$ pair takes the following form,

[TABLE]

where

[TABLE]

• •

Relevant partial decay widths of the scalars $h_{i}$ ($i=1,2$):

[TABLE]

where $O_{1}=\sin\alpha$, $O_{2}=\cos\alpha$, and the coupling $\lambda_{\nu^{1}_{H}\nu^{1}_{H}h_{i}}$ is defined in Eq. (37).