Neutrino Mass, Leptogenesis and FIMP Dark Matter in a ${\rm U}(1)_{\rm B-L}$ Model

TL;DR

This paper extends the Standard Model with a U(1)_{B-L} gauge symmetry to simultaneously explain neutrino masses, dark matter, and baryon asymmetry, proposing a viable FIMP dark matter candidate produced via freeze-in.

Contribution

It introduces a new U(1)_{B-L} model with specific scalar and neutrino sectors, detailing the production of dark matter through freeze-in and linking it to neutrino mass and leptogenesis.

Findings

Dark matter can be produced via freeze-in from RH neutrino annihilation.

The model successfully explains neutrino masses, dark matter relic abundance, and baryon asymmetry.

A parameter space where DM is mainly produced from RH neutrino annihilation was identified.

Abstract

The Standard Model (SM) is inadequate to explain the origin of tiny neutrino masses, the dark matter (DM) relic abundance and also the baryon asymmetry of the Universe. In this work to address all the three puzzles, we extend the SM by a local U$(1)_{\rm B-L}$ gauge symmetry, three right-handed (RH) neutrinos for the cancellation of gauge anomalies and two complex scalars having nonzero U$(1)_{\rm B-L}$ charges. All the newly added particles become massive after the breaking of U$(1)_{\rm B-L}$ symmetry by the vacuum expectation value (VEV) of one of the scalar fields $\phi_H$. The other scalar field $\phi_{DM}$, which does not have any VEV, becomes automatically stable and can be a viable DM candidate. Neutrino masses are generated using Type-I seesaw mechanism while the required lepton asymmetry to reproduce the observed baryon asymmetry, can be attained from the CP violating out of equilibrium decays of RH neutrinos in TeV scale. More importantly within this framework, we have studied in detail the production of DM via freeze-in mechanism considering all possible annihilation and decay processes. Finally, we find a situation when DM is dominantly produced from the annihilation of RH neutrinos, which are at the same time also responsible for neutrino mass generation and leptogenesis.

Figures

Tables

Table 1. Table 1: Charges of all particles under various symmetry groups.

Gauge Group SU(2)L U(1)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$ $2$ $1$ $1$ $-1/2$ $-1$ $0$ $-1$ $-1$ $-1$

Scalar Fields ${\varphi }_h$ ${\varphi }_H$ ${\varphi }_{DM}$ $2$ $1$ $1$ $1/2$ $0$ $0$ $0$ $2$ $n_{BL}$

Table 2. Table 2: Baryon asymmetry of the Universe generated for three different values of M N 1 subscript 𝑀 subscript 𝑁 1 M_{N_{1}} and ε 1 subscript 𝜀 1 \varepsilon_{1} .

$M_{N_1}$ [GeV]	${\epsilon }_1$	$Y_B=\frac{n_B}{S}$
$1600$	$4.4\times {10}^{-4}$	$8.7121\times {10}^{-11}$
$1800$	$2.25\times {10}^{-4}$	$8.7533\times {10}^{-11}$
$2000$	$1.8\times {10}^{-4}$	$8.5969\times {10}^{-11}$

Table 3. Table 3: Couplings of FIMP ( ϕ D ​ M subscript italic-ϕ 𝐷 𝑀 \phi_{DM} ) with Z B ​ L subscript 𝑍 𝐵 𝐿 Z_{BL} , h 1 subscript ℎ 1 h_{1} and h 2 subscript ℎ 2 h_{2} .

Vertex	Vertex Factor
$abc$	$g_{abc}$
${\varphi }_{DM}{\varphi }_{DM}^{†}Z_{BL}^{}{}_{\mu }{}^{}$	$g_{BL}{n}_{BL}{\left(p_2-p_1\right)}^{\mu }$
${\varphi }_{DM}{\varphi }_{DM}^{†}{h}_1$	$-\left({\lambda }_{Dh}v\cos \alpha +{\lambda }_{DH}{v}_{BL}\sin \alpha \right)$
${\varphi }_{DM}{\varphi }_{DM}^{†}{h}_2$	$\left({\lambda }_{Dh}v\sin \alpha -{\lambda }_{DH}{v}_{BL}\cos \alpha \right)$

Equations

$$Y_B=-\frac{8N_f+4N_{{\varphi }_h}}{22N_f+13N_{{\varphi }_h}}{Y}_{B-L}$$

$$\mathcal{L}$$

$${\mathcal{L}}_{DM}$$

$$D_{\mu }\psi =({\partial }_{\mu }+i{g}_{BL}{Q}_{BL}(\psi ){Z_{BL}}_{\mu })\psi ,$$

$$V({\varphi }_h,{\varphi }_H)={\mu }_H^2{\varphi }_H^{†}{\varphi }_H+{\mu }_h^2{\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).$$

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

$$M_{N_i}$$

$$M_{Z_{BL}}$$

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

$$h_1$$

$$h_2$$

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

$$M_{h_1}^2$$

$$M_{h_2}^2$$

$$M_{DM}^2$$

$${\lambda }_H$$

$${\lambda }_h$$

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

$${\mu }_{{\varphi }_h}^2$$

$${\mu }_{{\varphi }_H}^2$$

$${\mu }_{DM}^2$$

$${\mu }_{{\varphi }_h}^2<0,{\mu }_{{\varphi }_H}^2<0,{\mu }_{DM}^2>0,$$

$${\lambda }_h\ge 0,{\lambda }_H\ge 0,{\lambda }_{DM}\ge 0,$$

$${\lambda }_{hH}\ge -2\sqrt{{\lambda }_h{\lambda }_H},$$

$${\lambda }_{Dh}\ge -2\sqrt{{\lambda }_h{\lambda }_{DM}},$$

$${\lambda }_{DH}\ge -2\sqrt{{\lambda }_H{\lambda }_{DM}},$$

$$\sqrt{{\lambda }_{hH}+2\sqrt{{\lambda }_h{\lambda }_H}}\sqrt{{\lambda }_{Dh}+2\sqrt{{\lambda }_h{\lambda }_{DM}}}\sqrt{{\lambda }_{DH}+2\sqrt{{\lambda }_H{\lambda }_{DM}}}$$

$$+2\sqrt{{\lambda }_h{\lambda }_H{\lambda }_{DM}}+{\lambda }_{hH}\sqrt{{\lambda }_{DM}}+{\lambda }_{Dh}\sqrt{{\lambda }_H}+{\lambda }_{DH}\sqrt{{\lambda }_h}\ge 0.$$

$$\displaystyle\mathcal{M}_{D}=\left(\begin{array}[]{ccc}y_{ee}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ y_{e\mu}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ y_{e\tau}\\ \nobreak\leavevmode\nobreak\leavevmode\hfil&\nobreak\leavevmode\nobreak\leavevmode\hfil\\ y_{\mu e}+i\,\tilde{y}_{\mu e}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ y_{\mu\mu}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ y_{\mu\tau}\\ \nobreak\leavevmode\nobreak\leavevmode\hfil&\nobreak\leavevmode\nobreak\leavevmode\hfil\\ y_{\tau e}+i\,\tilde{y}_{\tau e}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ y_{\tau\mu}+i\,\tilde{y}_{\tau\mu}\nobreak\leavevmode\nobreak\leavevmode&\leavevmode\nobreak\ \leavevmode\nobreak\ y_{\tau\tau}\\ \end{array}\right)\,,$$

$$\displaystyle M=\left(\begin{array}[]{cc}0&\mathcal{M_{D}}\\ \mathcal{M_{D}}^{T}&\mathcal{M_{R}}\end{array}\right)\,\,.$$

$$m_{\nu }$$

$$m_N$$

$$J_{\mathrm{CP}}=\frac{1}{8}\sin 2{\theta }_{12}\sin 2{\theta }_{23}\sin 2{\theta }_{13}\cos {\theta }_{13}\sin \delta .$$

$$J_{CP}=\frac{\mathrm{Im}(h_{13}{h}_{23}{h}_{31})}{\Delta m_{21}^2\Delta m_{32}^2\Delta m_{31}^2}$$

$$\displaystyle\begin{array}[]{cccccc}1\,\,{\rm TeV}&\leq&M_{N_{1}}&\leq&3\,\,{\rm TeV}\,,&\\ M_{N_{1}}&<&M_{N_{3}}&\leq&15\,\,{\rm TeV}\,,&\\ 1&\leq&\dfrac{\sqrt{2}\,\,y_{ij}}{v}\,\times 10^{8}&\leq&1000&(i,j=e,\,\mu,\,\tau,i=j\neq e)\,,\\ 1&\leq&\dfrac{\sqrt{2}\,\,y_{ee}}{v}\,\times 10^{10}&\leq&100\,,&\\ 1&\leq&\dfrac{\sqrt{2}\,\,\tilde{y}_{ij}}{v}\,\times 10^{8}&\leq&1000&(i=\tau,j=e,\,\mu)\,,\\ 1&\leq&\dfrac{\sqrt{2}\,\,\tilde{y}_{\mu e}}{v}\,\times 10^{9}&\leq&1000\,.&\end{array}$$

$${\epsilon }_2$$

$${\epsilon }_1$$

$$\simeq$$

$$\frac{d{Y}_{N_1}}{dz}$$

Full text

I Introduction

The presence of non-zero neutrino mass and mixing has been confirmed by observing neutrino flavour oscillations Fukuda:1998mi ; Ahmad:2002jz among its different flavours. Neutrino experiments have measured the three intergenerational mixing angles ($\theta_{12},\theta_{23},\theta_{13}$) and the two mass square differences ($\Delta m_{21}^{2}$ and $\Delta m_{32}^{2}$) 111$\Delta m_{ij}^{2}=m_{i}^{2}-m_{j}^{2}\,\,\,\,{\rm and}\,\,\,\,\Delta m_{atm}^{2}=m_{3}^{2}-\frac{m_{1}^{2}+m_{2}^{2}}{2}$ with an unprecedented accuracy Eguchi:2002dm ; An:2015nua ; RENO:2015ksa ; Abe:2014bwa ; Abe:2015awa ; Salzgeber:2015gua ; Adamson:2016tbq ; Adamson:2016xxw . Neutrinos are massless in the Standard Model (SM) of particle physics because in SM there is no right handed (RH) counterpart of the left handed (LH) neutrinos. To generate tiny neutrino masses and their intergenerational mixing angles, as suggested by different experiments, we have to think of some new interactions and/or new particles beyond the Standard Model (BSM). Moreover, there are still some unsolved problems in the neutrino sector. For example, we do not know the exact octant of the atmospheric mixing angle $\theta_{23}$ i.e. whether it lies in the lower octant ($\theta_{23}<45^{\circ}$) or in the higher octant ($\theta_{23}>45^{\circ}$), the exact sign of $\Delta m^{2}_{32}$ which is related to the mass hierarchy between $m_{2}$ and $m_{3}$ (for the normal hierarchy (NH) $\Delta m^{2}_{32}>0$ while for the inverted hierarchy (IH) $\Delta m^{2}_{32}<0$) and also about the Dirac CP phase $\delta$, responsible for the CP violation in the leptonic sector. Recently, T2K and No$\nu$A experiments have reported their preliminary result which predicts that the value of Dirac CP phase is around $\delta_{CP}\sim 270^{\circ}$ xyz . Besides these, we do not know whether the neutrinos are Dirac fermion or Majorana fermion. Observation of neutrino less double $\beta$ decay Gando:2012zm ; Agostini:2013mzu ; Asakura:2014lma ; Albert:2014awa ; xyz1 ; KamLAND-Zen:2016pfg will confirm the Majorana nature of neutrinos and might also provide important information about the Majorana phases which could be the other source of CP violation in the leptonic sector, if the SM neutrinos are Majorana fermions.

Besides these unsolved problems in the neutrino sector, another well known puzzle in recent times is the presence of dark matter (DM) in the Universe. Many indirect evidence suggests the existence of DM. Among the most compelling evidence of DM are the observed flatness of rotation curves of spiral galaxies Sofue:2000jx , gravitational lensing Bartelmann:1999yn , the observed spatial offset between DM and visible matter in the collision of two galaxy clusters (e.g. Bullet cluster Clowe:2003tk , Abell cluster Biviano:1996bg ; Kahlhoefer:2013dca ) etc. The latter also imposes an upper bound on the ratio between self interaction and mass of DM particles, which is \frac{\sigma_{DM}}{M_{DM}}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}1\,{\rm barn}/{\rm GeV} Harvey:2015hha . Moreover, satellite borne experiments like WMAP Hinshaw:2012aka and Planck Ade:2015xua have made a precise measurement of the amount of dark matter present in the Universe from the cosmic microwave background (CMB) anisotropy Ade:2015xua and the current measured value of this parameter lies in the range $0.1172\leq\Omega_{\rm DM}h^{2}\leq 0.1226$ at 67% C.L Ade:2015xua .

Despite the compelling observational evidence for DM due to its gravitational interactions, our knowledge about its particle nature is very limited.  The only thing we know about the DM is that it is very weakly interacting and electromagnetically blind.  The SM of particle physics does not have any fundamental particle which can play the role of a cold dark matter (CDM), consequently a BSM scenario containing new fundamental stable particle(s) is required. There are earth based ongoing DM direct detection experiments, namely Xenon-1T Aprile:2012zx , LUX Akerib:2015rjg , CDMS Ahmed:2010wy ; Agnese:2014aze amongst others, which have been trying to detect the Weakly Interacting Massive Particle (WIMP) Gondolo:1990dk ; Jungman:1995df ; ArkaniHamed:2008qn type DM by measuring recoil energies of the detector nuclei scattered by the WIMPs. However, no convincing DM signal has been found yet and hence the $M_{DM}-\sigma_{\rm SI}$ plane for a WIMP type DM is now getting severely constrained. Therefore, invoking particle DM models outside the WIMP paradigm seems to be pertinent at this stage Arcadi:2017kky . In the present work we study one of the possible alternatives of WIMP, namely, the Feebly Interacting Massive Particle (FIMP) Hall:2009bx ; Yaguna:2011qn ; Molinaro:2014lfa ; Biswas:2015sva ; Merle:2015oja ; Shakya:2015xnx ; Biswas:2016bfo ; Konig:2016dzg ; Biswas:2016iyh ; Biswas:2016yjr . A major difference between the WIMP and FIMP scenarios is that while in the former the DM particle is in thermal equilibrium with the plasma in the early Universe and freezes-out when the Hubble expansion rate becomes larger than its annihilation cross section, in the FIMP case the DM is never in thermal equilibrium with the cosmic soup. This is mainly ensured by its extremely weak couplings to other particles in the thermal bath. Therefore, the number density of the FIMP is negligible in the early Universe and increases when the FIMP is subsequently produced by the decays and annihilations of other particles to which it is coupled (very feebly). This process is generally known as freeze-in Hall:2009bx .

In addition to the above two unsolved problems, another long standing enigma is the presence of more baryons over anti baryons in the Universe, which is known as the baryon asymmetry or the matter-antimatter asymmetry in the Universe. The baryon asymmetry observed in the Universe is expressed by a quantity $Y_{B}=\dfrac{\eta_{B}\,n_{\gamma}}{\rm s}$, where $\eta_{B}=n_{B}-n_{\bar{B}}$ is the excess in the number density for baryon over anti-baryon while $n_{\gamma}$ and ${\rm s}$ are the photon number density and the entropy density of the Universe, respectively. At the present epoch, $\eta_{B}=\left(5.8-6.6\right)\times 10^{-10}$ at 95% C.L. Olive:2016xmw while at $T\sim 2.73$ K, the photon density $n_{\gamma}=410.7\,\,{\rm cm}^{-3}$ Olive:2016xmw and the entropy density ${\rm s}=2891.2\,\,{\rm cm}^{-3}$ Olive:2016xmw (in natural unit with Boltzmann constant $K_{B}=1$). Therefore, the observed baryon asymmetry at the present Universe is $Y_{B}=\left(8.24-9.38\right)\times 10^{-10}$, which although small, is sufficient to produce the $\sim 5\%$ energy density (visible matter) of the Universe. To generate baryon asymmetry in the Universe from a matter-antimatter symmetric state, one has to satisfy three necessary conditions, known as the Sakharov’s conditions Sakharov:1967dj . These are i) baryon number (B) violation, ii) C and CP violation and iii) departure from thermal equilibrium. Since the baryon number (B) is an accidental symmetry of the SM (i.e. all SM interactions are B conserving) and also the observed CP violation in quark sector is too small to generate the requited baryon asymmetry, hence like the previous cases, here also one has to look for some additional BSM interactions which by satisfying the Sakharov’s conditions can generate the observed baryon asymmetry in an initially matter-antimatter symmetric Universe.

In this work, we will try to address all of the three above mentioned issues. The non-observation of any BSM signal at LHC implies the concreteness of the SM. However to address all the three problems, we need to extend the particles list and/or gauge group of SM because as already mentioned, SM is unable to explain either of them. In our model, we have extended the SM gauge group ${\rm SU}(3)_{\rm c}\times{\rm SU}(2)_{\rm L}\times{\rm U}(1)_{\rm Y}$ by a local ${\rm U}(1)_{\rm B-L}$ gauge group. The ${\rm B-L}$ extension of SM Mohapatra:1980qe ; Georgi:1981pg ; Wetterich:1981bx ; Lindner:2011it has been studied earlier in the context of dark matter phenomenology Okada:2010wd ; Okada:2012sg ; Basso:2012ti ; Basak:2013cga ; Sanchez-Vega:2014rka ; Basak:2014sza ; Guo:2015lxa ; Rodejohann:2015lca ; Okada:2016gsh ; Patra:2016ofq ; Okada:2016tci and baryogenesis in the early Universe in Refs. Buchmuller:1992qc ; Buchmuller:1996pa ; Dulaney:2010dj . Since we have imposed a local U(1) symmetry, consequently an extra gauge boson ($Z_{BL}$) will arise. To cancel the anomaly due to this extra gauge boson we need to introduce three right-handed (RH) neutrinos ($N_{i},i=1$, $2$, $3$) to make the model anomaly free. Apart from the three RH neutrinos, we have also introduced two SM gauge singlet scalars namely $\phi_{H}$ and $\phi_{DM}$, both of them are charged under the proposed ${\rm U}(1)_{\rm B-L}$ gauge group. The ${\rm U}(1)_{\rm B-L}$ symmetry is spontaneously broken when the scalar field $\phi_{H}$ takes a nonzero vacuum expectation value (VEV) and thereby generates the masses for the three RH neutrinos as well as the extra neutral gauge boson $Z_{BL}$, whose mass terms are forbidden initially due to the ${\rm U}(1)_{\rm B-L}$ invariance of the Lagrangian. The other scalar $\phi_{DM}$ does not acquire any VEV and by choosing appropriate ${\rm B-L}$ charge $\phi_{DM}$ becomes naturally stable and therefore, can serve as a viable dark matter candidate. As mention above, anomaly cancellation requires the introduction of three RH neutrinos in the present model. Therefore we can easily generate the neutrino masses by the Type-I seesaw mechanism after B-L symmetry is broken. Diagonalising the light neutrino mass matrix ($m_{\nu}$, for detail see Section III.1), we determine the allowed parameter space by satisfying the $3\sigma$ bounds on the mass square differences ($\Delta m_{12}^{2}$, $\Delta m_{atm}^{2}$), the mixing angles ($\theta_{12},\theta_{13},\theta_{23}$) Capozzi:2016rtj and also the cosmological bound on the sum of three light neutrinos masses Ade:2015xua . We also determine the effective mass $m_{\beta\beta}$ which is relevant for neutrino-less double beta decay and compare it against the current bound on $m_{\beta\beta}$ from GERDA phase I experiment Agostini:2013mzu .

Next, we explain the possible origin of the baryon asymmetry at the present epoch from an initially matter-antimatter symmetric Universe via leptogenesis. We first generate the lepton asymmetry (or ${\rm B-L}$ asymmetry, $Y_{B-L}$) from the out of equilibrium, CP violating decays of RH neutrinos. The lepton asymmetry thus produced has been converted into the baryon asymmetry by the ($B+L$) violating sphaleron processes which are effective before and during electroweak phase transition Manton:1983nd ; Klinkhamer:1984di ; Kuzmin:1985mm . When the sphaleron processes are in thermal equilibrium (10^{12}\,{\rm GeV}\,\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}T\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}10^{2} GeV, $T$ being the temperature of the Universe), the conversion rate is given by Khlebnikov:1988sr

[TABLE]

where $N_{f}=3$ and $N_{\phi_{h}}=1$, are the number of fermionic generations and number of Higgs doublet in the model, respectively.

Finally, in order to address the dark matter issue, we consider the singlet scalar $\phi_{DM}$ as a DM candidate. Since the couplings of this scalar to the rest of the particles of the model are free parameters, they could take any value. Depending on the value of these couplings, we could consider $\phi_{DM}$ as a WIMP or a FIMP. Detailed study on the WIMP type scalar DM in the present ${\rm U}(1)_{\rm B-L}$ framework has been done in Refs. Guo:2015lxa ; Rodejohann:2015lca ; Biswas:2016ewm . In most of the earlier works, it has been shown that the WIMP relic density is mainly satisfied around the resonance regions of the mediator particles. Moreover, the WIMP parameter space has now become severely constrained due to non-observation of any “real” signal in various direct detection experiments. Thus, as discussed earlier, in this situation the study of scalar DM other than WIMP is worthwhile. Therefore in this work, we consider the scalar field $\phi_{DM}$ as a FIMP candidate which, depending on its mass, is dominantly produced from the decays of heavy bosonic particles such as $h_{1}$, $h_{2}$, $Z_{BL}$ and also from the annihilations of bosonic as well as fermionic degrees of freedom present in the model (e.g. $N_{i}$, $Z_{BL}$, $h_{i}$ etc.). In particular, in Ref. Biswas:2016yjr , we have also studied a SM singlet scalar as the FIMP type DM candidate in a $L_{\mu}-L_{\tau}$ gauge extension of the SM. In that work, we have considered the extra gauge boson mass in MeV range to explain the muon $(g-2)$ anomaly. Consequently, the production of a $\mathcal{O}$(GeV) DM from the decay of $Z_{\mu\tau}$ is forbidden. Additionally, in that model due to the considered $L_{\mu}-L_{\tau}$ flavour symmetry the neutrino mass matrices (both light and heavy neutrinos) have particular shape. On the other hand, in the present work, we have extensively studied the FIMP DM production mechanism from all possible decays and annihilations other particles present in the model. Moreover, we have found that depending on our DM mass, a sharp correlation exists among the three puzzles of astroparticle physics namely neutrino mass generation, leptogenesis and DM. Furthermore, earlier in Ref.  Biswas:2016bfo , one of us, along with other collaborators, has studied the freeze-in DM production mechanism in the framework of U(1)B-L extension of the SM. However, in that article they have considered an MeV range RH neutrino as the FIMP DM candidate. Thus, in the context of DM phenomenology the current work is vastly different from Ref.  Biswas:2016bfo .

In the non-thermal scenario, most of the production of the FIMP from the decay of a heavy particle occurs when $T\sim M$, where $M$ is the mass of the decaying mother particle, which is generally assumed to be in thermal equilibrium. Therefore, the non-thermality condition of the FIMP demands that \dfrac{\Gamma}{H}<1\bigg{|}_{T\sim M} Arcadi:2013aba , which in turn imposes a severe upper bound on the coupling strengths of the FIMP. Thus the non-thermality condition requires extremely small coupling of $\phi_{DM}$ with the thermal bath (\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}10^{-10}) and hence, FIMP DM can easily evade all the existing bounds from DM direct detection experiments Aprile:2012zx ; Akerib:2015rjg ; Ahmed:2010wy .

Rest of the paper has been arranged in the following manner, in Section II we discuss the model in detail. In Section III we present the main results of the paper. In particular, we discuss the neutrino phenomenology in Section III.1, baryogenesis via leptogenesis in Section III.2 and non-thermal FIMP dark matter $\phi_{DM}$ production in Section III.3. Finally in Section IV we end with our conclusions.

II Model

The gauged ${\rm U}(1)_{\rm B-L}$ extension of SM is one of the most extensively studied BSM model so far. In this model, the gauge sector of the SM is enhanced by imposing a local ${\rm U}(1)_{\rm B-L}$ symmetry to the SM Lagrangian, where B and L represent the respective baryon and lepton number of a particle. Therefore, the complete gauged group is ${\rm SU}(3)_{\rm c}\times{\rm SU}(2)_{\rm L}\times{\rm U}(1)_{\rm Y}\times{\rm U}(1)_{\rm B-L}$. Since the ${\rm U}(1)_{\rm B-L}$ extension of SM is not an anomaly free theory, hence we need to introduce some chiral fermions to cancel the anomaly. In order to achieve this, we have considered three extra right handed (RH) neutrinos to make the proposed ${\rm B-L}$ extension anomaly free. Besides the SM particles and three RH neutrinos, we have introduced two SM gauge singlet scalars $\phi_{H}$, $\phi_{DM}$ in the theory with suitable ${\rm B-L}$ charges. One of the scalar fields namely $\phi_{H}$ breaks the proposed ${\rm U}(1)_{\rm B-L}$ symmetry spontaneously by acquiring a nonzero VEV $v_{BL}$ and thereby generates masses to all the BSM particles. We have chosen the ${\rm B-L}$ charge of $\phi_{DM}$ in such a way that the Lagrangian of our model before the ${\rm U}(1)_{\rm B-L}$ symmetry breaking does not contain any interaction term involving odd powers of $\phi_{DM}$. When $\phi_{H}$ gets a nonzero VEV, this ${\rm U}(1)_{\rm B-L}$ symmetry breaks spontaneously into a remnant $\mathbb{Z}_{2}$ symmetry under which only $\phi_{DM}$ becomes odd. The $\mathbb{Z}_{2}$ invariance of the Lagrangian will be preserved as long as the parameters of the Lagrangian are such that the scalar field $\phi_{DM}$ does not get any VEV. Under this condition, the scalar field $\phi_{DM}$ becomes absolutely stable and, in principle, can serve as a viable dark matter candidate. The respective SU(2)L, U(1)Y and ${\rm U}(1)_{\rm B-L}$ charges of all the particles in the present model are listed in Table 1.

The complete Lagrangian for the model is as follows,

[TABLE]

with $\tilde{\phi_{h}}=i\sigma_{2}\phi^{*}_{h}$. The term $\mathcal{L}_{SM}$ and $\mathcal{L}_{DM}$ represent the SM and dark sector Lagrangian, respectively. The dark sector Lagrangian $\mathcal{L}_{DM}$ containing all possible gauge invariant interaction terms of the scalar field $\phi_{DM}$, has the following form

[TABLE]

where the interactions of $\phi_{DM}$ with $\phi_{h}$ and $\phi_{H}$ are proportional to the couplings $\lambda_{Dh}$ and $\lambda_{DH}$, respectively. The fourth term in Eq. (2) represents the kinetic term for the additional gauge boson $Z_{BL}^{\mu}$ in terms of field strength tensor ${F_{BL}}_{\mu\nu}$ of the ${\rm U}(1)_{\rm B-L}$ gauge group. The covariant derivatives involving in the kinetic energy terms of the BSM scalars and fermions, $\phi_{H}$, $\phi_{DM}$ and $N_{i}$ (Eq. (2)), can be expressed in a generic form

[TABLE]

where $\psi=\phi_{DM},\phi_{H}$, $N_{i}$ and $Q_{BL}(\psi)$ represents the ${\rm B-L}$ charge of the corresponding field (listed in Table 1). The quantity $V(\phi_{h},\phi_{H})$ in Eq. (2) contains the self interaction terms of $\phi_{H}$ and $\phi_{h}$ as well as the mutual interaction term between the two scalar fields. The expression of $V(\phi_{h},\phi_{H})$ is given by

[TABLE]

After the symmetry breaking, the SM Higgs doublet $\phi_{h}$ and the BSM scalar $\phi_{H}$ take the following form,

[TABLE]

where $v$ = 246 GeV is the VEV of $\phi_{h}$, which breaks the SM gauge symmetry into a residual U(1)EM symmetry. The remaining terms in Eq. (2) are the Yukawa interaction terms for the left handed and right handed neutrinos. As mentioned in the beginning of this section, when the extra scalar field $\phi_{H}$ gets a nonzero VEV $v_{BL}$, the proposed ${\rm U}(1)_{\rm B-L}$ gauge symmetry breaks spontaneously. As a results, the Majorana mass terms for the right handed neutrinos, proportional to the Yukawa couplings $y_{N_{i}}$, are generated. In general, for a three generation of right handed neutrinos, we will have a $3\times 3$ Majorana mass matrix $\mathcal{M_{R}}$ with all off diagonal terms are present. However, in the present scenario for calculational simplicity, we have chosen a basis for the $N_{i}$ fields with respect to which $\mathcal{M_{R}}$ is diagonal. The diagonal elements, representing the masses of $N_{i}$s, are given by,

[TABLE]

Like the three right handed neutrinos, the extra neutral gauge bosons also becomes massive through the Eq. (4) when $\phi_{H}$ picks up a VEV. The mass term $Z_{BL}$ is given by

[TABLE]

When both $\phi_{h}$ and $\phi_{H}$ obtain their respective VEVs, there will be a mass mixing between the states $H$ and $H_{BL}$. The mass matrix with respect to the basis $H$ and $H_{BL}$ looks like as follows

[TABLE]

Rotating the basis states $H$ and $H_{BL}$ by a suitable angle $\alpha$, we can make the above mass matrix diagonal. The new basis states ($h_{1}$ and $h_{2}$) with respect to which the mass matrix $\mathcal{M}^{2}_{scalar}$ becomes diagonal, are some linear combinations of earlier basis states $H$ and $H_{BL}$. The new basis states, now representing two the physical states, are defined as

[TABLE]

where we denote $h_{1}$ as the SM-like Higgs boson while $h_{2}$ is playing the role of a BSM scalar field. The mixing angle between $H$ and $H_{BL}$ can be expressed in terms of the parameters of the Lagrangian (cf. Eq. (2)) as,

[TABLE]

Besides the two physical scalar fields $h_{1}$ and $h_{2}$, as mentioned earlier, there is another scalar field ($\phi_{DM}$) in the present model, which can play the role of a dark matter candidate. The masses of these three physical scalar fields $h_{1}$, $h_{2}$ and $\phi_{DM}$ are give below,

[TABLE]

where $M_{x}$222Throughout the paper we have kept the mass ($M_{h_{1}}$) of SM-like Higgs boson $h_{1}$ fixed at 125.5 GeV. denotes the mass of the corresponding scalar field $x$.

In this work, we choose $M_{h_{2}}$, $M_{DM}$, $n_{BL}$, $M_{N_{i}}$, $M_{Z_{BL}}$, $g_{BL}$, $\alpha$, $\lambda_{Dh}$, $\lambda_{DH}$ and $\lambda_{DM}$ as our independent set of parameters. The other parameters in the Lagrangian namely $\lambda_{h}$, $\lambda_{H}$, $\lambda_{hH}$, $\mu_{\phi_{h}}^{2}$ and $\mu_{\phi_{H}}^{2}$ can be expressed in terms of these variables as follows Biswas:2016ewm .

[TABLE]

where $v_{BL}$ is defined in terms of $M_{Z_{BL}}$ and $g_{BL}$ in Eq. (8).

As we already know, in the present scenario two of the three scalar fields namely $\phi_{h}$ and $\phi_{H}$ obtain VEVs. On the other hand, the remaining scalar field $\phi_{DM}$ does not have any VEV, which ensures its stability by preserving its $\mathbb{Z}_{2}$ odd parity. Therefore, the ground state of the system is ($\langle\phi_{h}\rangle$, $\langle\phi_{H}\rangle$, $\langle\phi_{DM}\rangle$) = ($v$, $v_{BL}$, 0). Now, such a ground state (vacuum) will be bounded from below when the following inequalities are satisfied simultaneously Biswas:2016ewm ,

[TABLE]

Besides the lower limits of $\lambda$s as described by the above inequalities, there are also upper limits on the Yukawa and quartic couplings arising from the perturbativity condition which demands that the Yukawa and scalar quartic couplings have to be less than $\sqrt{4\,\pi}$ ($y<\sqrt{4\,\pi}$) and $4\,\pi$ ($\lambda<4\,\pi$) respectively Chakrabarty:2015yia .

III Results

III.1 Neutrino Masses and Mixing

As mentioned earlier, the cancellation of both axial vector anomaly Adler:1969gk ; Bardeen:1969md and gravitational gauge anomaly Delbourgo:1972xb ; Eguchi:1976db , in ${\rm U}(1)_{\rm B-L}$ extended SM, requires the presence of extra chiral fermions. Hence, in the present model to cancel these anomalies we have introduced three right handed (RH) neutrinos ($N_{i}$, $i$=1 to 3). The Majorana masses for the RH neutrinos are generated only after spontaneous breaking of the proposed ${\rm B-L}$ symmetry by the VEV of $\phi_{H}$. Also in the present scenario, as stated earlier, we are working in a basis where the Majorana mass matrix for the three RH neutrinos are diagonal i.e. $\mathcal{M_{R}}={\rm diag}\,(M_{N_{1}},M_{N_{2}},M_{N_{3}})$. The expression for the mass of $i$th RH neutrino ($M_{N_{i}}$) is given in Eq. (7). On the other hand, the Dirac mass terms involving both left chiral and right chiral neutrinos, are originated when the electroweak symmetry is spontaneously broken by the VEV of SM Higgs doublet $\phi_{h}$, giving rise to a $3\times 3$ complex matrix $\mathcal{M_{D}}$. In general, one can take all the elements of matrix $\mathcal{M_{D}}$ as complex but for calculational simplicity and also keeping in mind that only three physical phases (one Dirac phase and two Majorana phases) exist for three light neutrinos (Majorana type), we have considered only three complex elements in the lower triangle part of the Dirac mass matrix $\mathcal{M_{D}}$. However, the results we have presented later in this section will not change significantly if we consider all the elements of $\mathcal{M_{D}}$ are complex. The Dirac mass matrix $\mathcal{M_{D}}$ we assume has the following structure:

[TABLE]

where $y_{ij}=\dfrac{y_{ij}^{\prime}}{\sqrt{2}}v$ ($i,j=$ e, $\mu$, $\tau$) and the Yukawa coupling $y_{ij}^{\prime}$ has been defined in Eq. (2).

Now, with respect to the Majorana basis $\left(\overline{{\nu_{\alpha}}_{L}}\leavevmode\nobreak\ \leavevmode\nobreak\ \overline{({N_{\alpha}}_{R})^{c}}\right)^{T}$ and $\left(({\nu_{\alpha}}_{L})^{c}\leavevmode\nobreak\ \leavevmode\nobreak\ {N_{\alpha}}_{R}\right)^{T}$ one can write down the Majorana mass matrix for both left and right chiral neutrinos using $\mathcal{M_{D}}$ and $\mathcal{M_{R}}$ matrices in the following way,

[TABLE]

Since $M_{D}$ and $M_{R}$ are both $3\times 3$ matrices (for three generations of neutrinos), the resultant matrix $M$ will be of order $6\times 6$ and also it is a complex symmetric matrix which reflects its Majorana nature. Therefore, after diagonalisation of the matrix $M$, we get three light and three heavy neutrinos, all of which are Majorana fermions. If we use the block diagonalisation technique, we can write the light and heavy neutrino mass matrices in the leading order as,

[TABLE]

Here $M_{R}$ is a diagonal matrix and the expression of all the elements of $m_{\nu}$ in terms of the elements of $\mathcal{M_{D}}$ and $\mathcal{M_{R}}$ matrices are given in Appendix A. After diagonalising ${m_{\nu}}$ matrix we get three light neutrino masses ($m_{i}$, $i=1,\,2,\,3$), three mixing angles ($\theta_{12}$, $\theta_{13}$ and $\theta_{23}$) and one Dirac CP phase $\delta$.

We have used the Jarlskog Invariant $J_{\rm CP}$ Jarlskog:1985ht to determine the Dirac CP phase $\delta$, which is defined as,

[TABLE]

Moreover, the quantity $J_{\rm CP}$ is related to the elements of the Hermitian matrix $h=m_{\nu}m^{\dagger}_{\nu}$ in the following way,

[TABLE]

where in the numerator ${\rm Im}(X)$ represents the imaginary part of $X$ while in the denominator, $\Delta m^{2}_{ij}=m^{2}_{i}-m^{2}_{j}$. Once we determine the quantity $J_{\rm CP}$ (from Eq. (30)) and the intergenerational mixing angles of neutrinos then one can easily determine the Dirac CP phase using Eq. (29).

In the present scenario we have twelve independent parameters coming from the Dirac mass matrix. The RH neutrino mass matrix, in principle, should bring in three additional parameters. However, as we will discuss in details in Section III.2, two of the RH neutrino masses is taken to be nearly degenerate. In particular, the condition of resonant leptogenesis requires that $M_{N_{2}}-M_{N_{1}}=\Gamma_{1}/2$, where $\Gamma_{1}$ is the tree level decay width of $N_{1}$ and is seen to be $\sim 10^{-11}$ GeV. Therefore, for all practical purposes we have $M_{N_{1}}\simeq M_{N_{2}}$, and the RH neutrino mass matrix only brings in two independent parameters, $M_{N_{1}}$ and $M_{N_{3}}$. Thus, we have fourteen independent parameters which we vary in the following ranges,

[TABLE]

We try to find the allowed parameter space which satisfy the following constraints on three mixing angles ($\theta_{ij}$) and two mass square differences ($\Delta m_{ij}^{2}$), $J_{CP}$ obtained from neutrino oscillation data and the cosmological bound on the sum of three light neutrino masses. These experimental/observational results are listed below.

• •

Measured values of three mixing angles in $3\sigma$ range Capozzi:2016rtj :

$30^{\circ}<\,\theta_{12}\,<36.51^{\circ}$, $37.99^{\circ}<\,\theta_{23}\,<51.71^{\circ}$ and $7.82^{\circ}<\,\theta_{13}\,<9.02^{\circ}$.

• •

Allowed values of two mass squared differences in $3\sigma$ range Capozzi:2016rtj :

$6.93<\dfrac{\Delta m^{2}_{21}}{10^{-5}}\,{\text{eV}^{2}}<7.97$ and $2.37<\dfrac{\Delta m^{2}_{31}}{10^{-3}}\,{\text{eV}^{2}}<2.63$ in $3\sigma$ range.

• •

Above mentioned values of the neutrino oscillation parameters also put an upper bound on the absolute value of $J_{\rm CP}$ from Eq. (29), which is $|J_{\rm CP}|\leq\,0.039$.

• •

Cosmological upper bound on the sum of three light neutrino masses i.e. $\sum_{i}m_{i}<0.23$ eV at $2\sigma$ C.L. Ade:2015xua .

While it is possible to obtain both normal hierarchy (NH) ($m_{1}<m_{2}<m_{3}$) and inverted hierarchy (IH) ($m_{3}<m_{1}<m_{2}$) in this scenario, we show our results only for NH for brevity. Similar results can be obtained for IH.

In the LP of Fig. 1, we have shown the variation of $J_{\rm CP}$ parameter (as defined in Eq. (29)) with the Dirac CP phase $\delta$. From this plot one can easily notice that there are two allowed ranges of Dirac CP phase $0^{\circ}\leq\delta\leq 90^{\circ}$ and $270^{\circ}\leq\delta\leq 360^{\circ}$ respectively which can reproduce the neutrino oscillation parameters in $3\sigma$ range. Since the Jarlskog invariant $J_{\rm CP}$ is proportional to $\sin\delta$ (Eq.  29), hence we get both positive and negative values of $J_{CP}$ symmetrically placed in the first and fourth quadrants. However, the absolute values of $J_{\rm CP}$ always lie below 0.039. Also, here we want to mention that from the recent results of T2K Abe:2017vif experiment, values of $\delta$ lying in the fourth quadrant are favourable compared to those in first quadrant. In the RP of Fig.  1, we have shown the variation of neutrino less double $\beta$ decay parameter $m_{\beta\beta}$ with the mass of lightest neutrino $m_{1}$. $m_{\beta\beta}$ is an important quantity for the study of neutrino less double $\beta$ decay as the cross section of this process is proportional to $m_{\beta\beta}=\left|\sum_{i=1}^{3}({U_{\rm PMNS}})^{2}_{\,e\,i}\,\,m_{i}\right|=(m_{\nu})_{e\,e}$ (see Appendix B for details.), where $(m_{\nu})_{e\,e}$ (Eq.  56) is the (1,1) element of light neutrino mass matrix $m_{\nu}$. The nature of this plot is very to similar to the usual plot in $m_{\beta\beta}-m_{1}$ plane for the normal hierarchical scenario DellOro:2016tmg . In the same plot, we have also shown the current bound on $m_{\beta\beta}$ from KamLand-Zen experiment KamLAND-Zen:2016pfg .

III.2 Baryogenesis via Resonant Leptogenesis

As we have three RH neutrinos in the present model, in this section we have studied the lepton asymmetry generated from the CP violating out of equilibrium decays of these heavy neutrinos at the early stage of the Universe. The ${\rm B-L}$ asymmetry thus produced is converted into the baryon asymmetry through sphaleron transitions which violate ${\rm B+L}$ quantum number while conserving the ${\rm B-L}$ charge. The sphaleron processes are active between temperatures of $\sim 10^{12}$ GeV to $\sim 10^{2}$ GeV in the early Universe. At high temperatures the sphalerons are in thermal equilibrium and subsequently they freeze-out at around $T\simeq 100-200$ GeV Plumacher:1996kc ; Iso:2010mv , just before electroweak symmetry breaking (EWSB). To produce sufficient lepton asymmetry, which would eventually be converted into the observed baryon asymmetry, one requires RH neutrinos with masses \mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr>\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr>\cr\sim\cr}}}}10^{8}-10^{9} GeV Plumacher:1996kc ; Buchmuller:2002rq . This is the well know scenario of the “normal” or “canonical” leptogenesis. However, detection of these very massive RH neutrinos is beyond the reach of LHC and other future colliders. Here we consider the RH neutrinos to be in TeV mass range to allow for their detection at collider experiments. It has been shown that with RH neutrinos in the TeV mass-scale range, it is possible to generate adequate lepton asymmetry by considering the two lightest RH neutrinos $N_{1}$ and $N_{2}$ to be almost degenerate. More specifically, we demand that $M_{N_{2}}-M_{N_{1}}\simeq{\Gamma_{1}}/{2}$, where $\Gamma_{1}$ 333The typical value of $\Gamma_{1}$ is $\sim 10^{-11}$ GeV (see Fig. 3) while $M_{N_{i}}\sim\mathcal{O}$(TeV). Hence we take $M_{N_{1}}=M_{N_{2}}$ throughout the work. is the total decay width of the lightest RH neutrino $N_{1}$. This scenario is known as Resonant leptogenesis Pilaftsis:1997jf ; Pilaftsis:2003gt ; Iso:2010mv ; Heeck:2016oda .

Fig.  2 shows the tree level as well as one loop decay diagrams of the lightest RH neutrino $N_{1}$. These diagrams are applicable for all the three RH neutrinos. Here $L$ represents the SM lepton which can either be a charged lepton or a left chiral neutrino depending on the nature of the scalar field (charged 444Since these processes occurred before EWSB hence we have both charged as well as neutral scalars in the SM. or neutral) associated in the vertex while $N_{j}$ denotes the remaining two RH neutrinos, $N_{2}$ and $N_{3}$ for the case of $N_{1}$ decay. In order to produce baryon asymmetry in the Universe we need both C and CP violating interactions, which is one of the three necessary conditions (see Sakharov conditions Sakharov:1967dj given in Section I) for baryogenesis. Lepton asymmetry generated from the out of equilibrium decay of RH neutrinos is determined by the CP asymmetry parameter ($\varepsilon_{i}$), which is given by (for details see Appendix C),

[TABLE]

In the LP of Fig.  3, we show the variation of CP asymmetry parameter $\varepsilon_{1}$, generated from the decay of RH neutrino $N_{1}$, with the mass of $N_{1}$. Here we see that for the considered ranges of $M_{N_{1}}$ ($1000\,{\rm GeV}\leq M_{N_{1}}\leq 10000\,{\rm GeV}$) and other relevant Yukawa couplings (see Eq. (LABEL:para-ranges)), the CP asymmetry parameter $\varepsilon_{1}$ can be as large as $\sim 10^{-2}$, which is significantly large compared to $\varepsilon_{1}$ in the “normal” Leptogenesis case ($\varepsilon_{1}\sim 10^{-8}$ for $M_{N_{1}}\sim 10^{10}$ GeV) Plumacher:1996kc . In the RP of Fig.  3, we plot the variation of total decay width of $N_{1}$ with $M_{N_{1}}$. From this plot, one can easily notice that in the present scenario, $\Gamma_{1}$ lies between $\sim 10^{-12}$ GeV to $10^{-9}$ GeV for the entire considered range of $M_{N_{1}}$. All the points in both panels satisfy the neutrino oscillations data in the $3\sigma$ range while the black solid line in the RP provides the upper bound on $\Gamma_{1}$, obtained from the out of equilibrium conditions for $N_{1}$ i.e. $\Gamma_{1}<3\,H(M_{N_{1}})$ Plumacher:1996kc where $H$ is Hubble parameter at $T=M_{N_{1}}$.

Next, we calculate the ${\rm B-L}$ asymmetry generated from the decays as well as the pair annihilations of the RH neutrinos $N_{1}$ and $N_{2}$. In order to calculate the net ${\rm B-L}$ asymmetry produced from the interactions of $N_{1}$ and $N_{2}$ at temperature of the Universe $T\simeq 150$ GeV (freeze-out temperature of sphaleron) we have to solve a set of three coupled Boltzmann equations. The relevant Boltzmann equations Plumacher:1996kc ; Iso:2010mv for calculating $Y_{N_{i}}$ and $Y_{B-L}$ are given below,

[TABLE]

where $Y_{X}=\dfrac{n_{X}}{\rm s}$ denotes the comoving number density of $X$, with $n_{X}$ being the actual number density and $z=\dfrac{M_{N_{1}}}{T}$. Planck mass is denoted by $M_{pl}$. The quantity $g_{\star}(z)$ is a function of $g_{\rho}$ and $g_{\rm s}$, the effective degrees of freedom related to the energy and entropy densities of the Universe respectively, and it has the following expression Gondolo:1990dk ,

[TABLE]

Before EWSB, the variation of $g_{s}(z)$ with respect to $z$ is negligible compared to the first term within the brackets and hence one can use $\sqrt{g_{\star}(z)}\simeq\dfrac{g_{\rm s}(z)}{\sqrt{g_{\rho}(z)}}$. The equilibrium comoving number density of $X$ ($X=N_{i}$, $L$), obeying the Maxwell Boltzmann distribution, is given by Gondolo:1990dk

[TABLE]

where $g_{X}$ and $M_{X}$ are the internal degrees of freedom and mass of $X$ respectively while $g_{s}\left(\frac{M_{N_{1}}}{z}\right)$ is the effective degrees of freedom related to the entropy density of the Universe at temperature $T=\dfrac{M_{N_{1}}}{z}$. ${\rm K}_{2}\left(\frac{M_{X}}{M_{N_{1}}}\,z\right)$ is the modified Bessel function of order 2. The relevant Feynman diagrams including both decay and annihilation of $N_{i}$ are shown in Figs.  2 and 4. The expression of thermal averaged decay width $\langle{\Gamma_{i}}\rangle$, which is related to total decay width $\Gamma_{i}$ of $N_{i}$ is given as

[TABLE]

The thermally average annihilation cross sections ${\langle\sigma{\rm v}\rangle}_{N_{i},\,Z_{BL}}$ and ${\langle\sigma{\rm v}\rangle}_{N_{i},\,Z_{BL}}$, appearing in Boltzmann equations (Eqs. (41) and (42)) for the processes shown in Fig.  4, can be defined in a generic form,

[TABLE]

where the $\hat{\sigma}_{N_{i},\,x}$ is related to the actual annihilation cross section ${\sigma}_{N_{i},\,x}$ by the following relation

[TABLE]

where $g_{N_{i}}=2$ is the internal degrees of freedom of RH neutrino $N_{i}$. The expression of $\hat{\sigma}_{N_{i},\,Z_{BL}}$ and $\hat{\sigma}_{N_{i},\,t,\,H_{BL}}$ for the present model is given in Ref. Iso:2010mv .

To calculate the ${\rm B-L}$ asymmetry at around $T\simeq 150$ GeV, we have to numerically solve the set of three coupled Boltzmann equations (Eqs. (41)-(43)) using Eqs. (44)-(48). However, we can reduce the two flavour analysis (when both $N_{2}$ and $N_{1}$ are separately considered) into one flavour case by considering the parameters of $\mathcal{M_{D}}$ matrix in such a way so that the decay widths of $N_{1}$ and $N_{2}$ are of the same order i.e. $\Gamma_{1}\sim\Gamma_{2}$. Hence, the CP asymmetry generated from the decays of both $N_{1}$ and $N_{2}$ are almost identical ($\varepsilon_{1}\sim\varepsilon_{2}$, see Eq. (70)-(71)). In this case, the net ${\rm B-L}$ asymmetry is equal to twice of that is being generated from the CP violating interactions of the lightest RH neutrino $N_{1}$ Iso:2010mv . Hence instead of solving three coupled differential equations we now only need to solve Eqs. (41) and (43). The results we have found by numerically solving Eqs. (41) and (43) are plotted in Fig.  5. In this plot, we have shown the variation of $Y_{N_{1}}$ and $Y_{B-L}$ with $z$ for $M_{N_{1}}=2000$ GeV, $\alpha_{BL}=3\times 10^{-4}$ and $M_{Z_{BL}}=3000$ GeV 555The considered value of $M_{Z_{BL}}$ and the corresponding gauge coupling $g_{BL}$ satisfy the upper bounds obtained from LEP Carena:2004xs ; Cacciapaglia:2006pk and more recently from LHC Guo:2015lxa as well.. While solving the coupled Boltzmann equations we have considered the following initial conditions: $Y_{N_{1}}(T_{\rm in}^{B})=Y^{\rm eq}_{N_{1}}$ and $Y_{B-L}=0$ with $T_{\rm in}^{B}$ is the initial temperature which we have taken as 20 TeV. Thereafter, the evolutions of $Y_{N_{1}}$ and $Y_{B-L}$ are governed by their respective Boltzmann equations. From Fig.  5, one can notice that initially upto $z\sim 1$ ($T\sim M_{N_{1}}$), the comoving number density of $Y_{N_{1}}$ does not change much as a result of the ${\rm B-L}$ asymmetry produced from the decay, and the annihilation of $N_{1}$ is also less. However, as the temperature of the Universe drops below the mass of $M_{N_{1}}$, there is a rapid change in the number density of $N_{1}$, which changes around six orders of magnitude between $z=1$ and $z=20$. Consequently, the large change in $Y_{N_{1}}$ significantly enhances the ${\rm B-L}$ asymmetry $Y_{B-L}$ and finally $Y_{B-L}$ saturates to the desired value around $\sim 10^{-10}$, when there are practically no $N_{1}$ left to produce any further ${\rm B-L}$ asymmetry.

The produced $B-L$ asymmetry is converted to net baryon asymmetry of the Universe through the sphaleron transitions while they are in equilibrium with the thermal bath. The quantities $Y_{B-L}$ and $Y_{B}$ are related by the following equation Khlebnikov:1988sr

[TABLE]

where $T_{\rm f}\simeq 150$ GeV is the temperature of the Universe upto which the sphaleron process, converting $B-L$ asymmetry to a net $B$ asymmetry, maintains its thermal equilibrium. The extra factor of two in the above equation is due to the equal contribution to $Y_{B-L}$ arising from the CP violating interactions of $N_{2}$ as well. Finally, we calculate the net baryon asymmetry $Y_{B}$ for three different masses of RH neutrino $N_{1}$ and CP asymmetry parameter $\varepsilon_{1}$. The results are listed in Table 2. In all three cases, the final baryon asymmetry lies within the experimentally observed range for $Y_{B}$ i.e. $(8.239-9.375)\times 10^{-11}$ at 95% C.L. Olive:2016xmw .

III.3 FIMP Dark Matter

In the present section we explore the FIMP scenario for dark matter in the Universe, by considering the complex scalar field $\phi_{DM}$ as a corresponding candidate. As described in the Section II, the residual $\mathbb{Z}_{2}$ symmetry of $\phi_{DM}$ makes the scalar field absolutely stable over the cosmological time scale and hence can play the role of a dark matter candidate. Since $\phi_{DM}$ has a nonzero ${\rm B-L}$ charge $n_{BL}$, therefore DM talks to the SM as well as the BSM particles through the exchange of extra neutral gauge boson $Z_{BL}$ and two Higgs bosons present in the model, one is the SM-like Higgs $h_{1}$ while other one is the BSM Higgs $h_{2}$. The corresponding coupling strengths, in terms of gauge coupling $g_{BL}$, ${\rm B-L}$ charge $n_{BL}$, mixing angle $\alpha$ and $\lambda$s, are listed in Table  3. As the FIMP never enters into thermal equilibrium, these couplings have to be extremely feeble in order to make the corresponding interactions nonthermal. For the case of $\phi_{DM}\,\phi^{\dagger}_{DM}\,{Z_{BL}}_{\mu}$ coupling, we will make the ${\rm B-L}$ charge of $\phi_{DM}$ extremely tiny so that this interaction enters into the nonthermal regime. In principle, one can also choose the gauge coupling $g_{BL}$ to be very small, however in the present case we will keep the values of $g_{BL}$ and $M_{Z_{BL}}$ fixed at 0.07 and 3 TeV respectively as these values reproduce the observed baryon asymmetry of the Universe (see Section III.2). Also, there is another advantage of choosing tiny $n_{BL}$ as this will make only $\phi_{DM}$ out of equilibrium while keeping $Z_{BL}$ in equilibrium with the thermal bath. Moreover, due to the nonthermal nature, the initial number density of FIMP is assumed to be negligible and as the temperature of the Universe begins to fall down, they start to be produced dominantly from the decays and annihilation of other heavy particles.

In the present scenario, we have considered all the particles except $\phi_{DM}$ to be in thermal equilibrium. Before EWSB, all the SM particles are massless666Although the SM particles acquire thermal masses before EWSB, we have neglected these masses, as in this regime this approximation will not affect the DM production processes significantly.. In this regime, production of $\phi_{DM}$ occurs mainly from the decay and/or annihilation of BSM particles namely $Z_{BL}$, $H_{BL}$, and $N_{i}$. Also, before EWSB the annihilation of all four degrees of freedom of SM Higgs doublet $\phi_{h}$ can produce $\phi_{DM}$. Feynman diagrams for all the production processes of $\phi_{DM}$ before EWSB are shown in Fig. 6.

After EWSB, all the SM particles become massive and consequently besides the BSM particles, $\phi_{DM}$ can now also be produced from the decay and/or annihilation of the SM particles as well. The corresponding Feynman diagrams are shown in Fig. 7. In generating the vertex factors for different vertices to compute the Feynman diagrams as listed in Fig. 6 and Fig. 7 we have used the LanHEP Semenov:2010qt package.

In order to compute the relic density of a species at the present epoch, one needs to study the evolution of the number density of the corresponding species with respect to the temperature of the Universe. The evolution of the number density of $\phi_{DM}$ is governed by the Boltzmann equation containing all possible number changing interactions of $\phi_{DM}$. The Boltzmann equation of $\phi_{DM}$ in terms of its comoving number density $Y_{\phi_{DM}}=\dfrac{n_{\phi_{DM}}}{\rm s}$, where $n$ and ${\rm s}$ are actual number density and entropy density of the Universe is given by

[TABLE]

where $z=\dfrac{M_{h_{1}}}{T}$, while $\sqrt{g_{\star}(z)}$, $g_{\rm s}(z)$ and $M_{pl}$ are same as those in Eqs. (41)–(43) of Section III.2. In the above equation (Eq. (50)), first term represents the contribution coming from the decays of $Z_{BL}$, $h_{1}$ and $h_{2}$. The expressions of equilibrium number density $Y^{\rm eq}_{X}(z)$ ($X$ is any SM or BSM particle expect $\phi_{DM}$) and the thermal averaged decay width $\langle\Gamma_{X\rightarrow\phi_{DM}\phi_{DM}}\rangle$ can be obtained from Eqs. (45) and (46), respectively by only replacing $M_{N_{1}}$ with $M_{X}$, the mass of decaying mother particle. As mentioned above, before EWSB, the summation in the first terms is over $h_{2}$ and $Z_{BL}$ only, as there will be no contribution from the SM Higgs decay as such trilinear vertex ($h_{1}\phi_{DM}\phi_{DM}^{\dagger}$) is absent before EWSB and after EWSB there will be contributions to the relic density of $\phi_{DM}$ from all there decaying particles. The dark matter production from the pair annihilations of SM and BSM particles are described by the second term of the Boltzmann equation. Here, summation over $p$ includes all possible pair annihilation channels namely $W^{+}W^{-},\,ZZ,\,Z_{BL}Z_{BL},\,N_{i}N_{i},\,h_{i}h_{i},\,t\bar{t}$. However before EWSB, pair annihilations of the BSM particles and SM Higgs doublet $\phi_{h}$ contribute to the production processes (i.e. $p=Z_{BL},\,N_{i},\,H_{BL},\,\phi_{h}$, see Fig.  6). The third term, which is present only after EWSB, is another the production mode of $\phi_{DM}$ from the annihilation of $h_{1}$ and $h_{2}$. The expressions of all the relevant cross sections and decay widths for computing the DM number density are given in Appendix E. The most general form of thermally averaged annihilation cross section for two different annihilating particles of mass $M_{A}$ and $M_{B}$ is given by Biswas:2016yjr ,

[TABLE]

Finally, the relic density of $\phi_{DM}$ is obtained using the following relation between $\Omega h^{2}$ and $Y_{\phi_{DM}}(0)$ Edsjo:1997bg ; Biswas:2011td ,

[TABLE]

where $Y_{\phi_{DM}}(0)$ is the value of comoving number density at the present epoch, which can be obtained by solving the Boltzmann equation.

The contribution to dark matter production processes from decays as well as annihilations of various SM and BSM particles depend on the mass of $\phi_{DM}$. Accordingly, We have divided our rest of the dark matter analysis into four different regions depending on $M_{DM}$ and the dominant production modes of $\phi_{DM}$.

III.3.1 $M_{DM}<\dfrac{M_{h_{1}}}{2}$,

 $\dfrac{M_{h_{2}}}{2}$, $\dfrac{M_{Z_{BL}}}{2}$, SM and BSM particles decay dominated region.

In this case DM is dominantly produced from the decays of all three particles namely $h_{1}$, $h_{2}$ and $Z_{BL}$. Therefore, in this case ${\rm U}(1)_{\rm B-L}$ part of the present model directly enters into the dark matter production. Moreover in this mass range, $\phi_{DM}$ can also be produced from the annihilations of SM and BSM particles, however, we find that their contributions are not as significant as those from the decays of $h_{1},\,h_{2}$ and $Z_{BL}$. In the left panel (LP) and right panel (RP) of Fig. 8, we have shown the variation of DM relic density with $z$. In LP, we have shown the dependence of DM relic density with the initial temperature $T_{in}$. Initial temperature ($T_{in}$) is the temperature upto which we have assumed that the number density of DM is zero and its production processes start thereafter. We can clearly see from the figure that as long as the initial temperature is above the mass of BSM Higgs ($M_{h_{2}}\sim 500$ GeV), the final relic density does not depend on the choice of the initial temperature and reproduces the observed DM relic density of the Universe for the chosen values of model parameters as written in the caption of Fig. 8.

If we reduce the initial temperature from $500$ GeV, i.e. for $T_{in}=251$ GeV, the decay contribution of BSM Higgs $h_{2}$ becomes less since corresponding the number density of $h_{2}$ for $T_{in}<M_{h_{2}}$ is Boltzmann suppressed (exponentially suppressed), which is clearly shown by the blue dashed-dotted line. Hence, if we reduce the initial temperature ($T_{in}$) further i.e. $T_{in}<M_{h_{2}},M_{h_{1}}$ $\sim 42$ GeV then the number densities of both SM-like Higgs $h_{1}$ as well as BSM Higgs $h_{2}$ become Boltzmann suppressed and hence, less amount of DM production will take place which is evident from the LP of Fig. 8 (represented by the yellow dashed-dot line). On the other hand in the RP of Fig. 8, we have shown the contributions to DM relic density coming from decay and annihilation. Magenta dotted horizontal line represents the present day observed DM relic density of the Universe. Green dashed line represents the total decay contribution arising from the decays of both $h_{1}$, $h_{2}$ and $Z_{BL}$ whereas the net annihilation contribution coming from the annihilation of all the SM as well as BSM particles has been shown by the blue dashed-dotted line. There is a sudden rise in the annihilation contribution which occurs around the Universe temperature $T\sim 154$ GeV (i.e. EWSB temperature). After the EWSB temperature, all the SM particles become massive and hence the sudden rise in the annihilation part because of the appearance of the following annihilation channels $W^{+}\,W^{-}$, $Z\,Z$, $h_{1}\,h_{1}$, $h_{1}\,h_{2}$. The plot clearly implies that the lion share of the contribution comes from the decay of both Higgses $h_{1}$, $h_{2}$ and $Z_{BL}$, while for the considered values of model parameters the annihilation contribution is subdominant. Moreover, in this case we cannot enhance the annihilation contribution by increasing parameters $\lambda_{Dh}$, $\lambda_{DH}$ and $n_{BL}$ as these changes will result in the over production of dark matter from the decays of $h_{1}$, $h_{2}$ and $Z_{BL}$.

In the LP of Fig. 9, we have shown how the individual decay contribution from each scalar varies with $z$. Here we consider the following values of the scalar quartic couplings $\lambda_{Dh}=8.75\times 10^{-13}$ and $\lambda_{DH}=5.88\times 10^{-14}$ and the ($B-L$) charge of $\phi_{DM}$ $n_{BL}=1.33\times 10^{-10}$. From this plot we can see that before EWSB SM-like Higgs $h_{1}$ cannot decay to a pair of $\phi_{DM}$ as in this epoch it has no coupling with the latter. In this regime the decay of BSM Higgs $h_{2}$ and $Z_{BL}$ contribute, while after EWSB even the SM-like Higgs starts contributing to the DM production and hence we get an increased relic density (right side of EWSB). Its worth mentioning here that while generating the plot in the LP of Fig.  9, we have taken the scalar quartic couplings $\lambda_{Dh}$, $\lambda_{DH}$ and ${\rm B-L}$ charge of $\phi_{DM}$ $n_{BL}$ of different strengths such that the contributions of both the scalars ($h_{1}$ and $h_{2}$) and the extra gauge boson to the DM relic density are of equal order. This is because for the case of BSM Higgs $h_{2}$ decay the coupling $\lambda_{DH}$ multiplied by the ${\rm B-L}$ symmetry breaking VEV $v_{BL}$ is relevant, while for the decay of the SM-like Higgs $h_{1}$, the product of the parameter $\lambda_{Dh}$ and the EWSB VEV $v$ is relevant and the contribution from the decay of $Z_{BL}$, DM charge $n_{BL}$ is relevant. Since in the present case $v_{BL}>v$, the magnitudes of the two quartic couplings $\lambda_{Dh}$ and $\lambda_{DH}$ are of different order (see Table 3). On the other hand, in the RP of Fig.  9, we have shown the variation of the relic density with $z$ for four different values of the DM mass $M_{DM}$. From Eq. (52), one can see that the DM relic density is directly proportional to the mass $M_{DM}$ and as a result when other relevant couplings remain unchanged $\Omega h^{2}$ increases with $M_{DM}$. This feature is clearly visible in the RP for the cases with $M_{DM}=10$ GeV (black solid line), $M_{DM}=30$ GeV (red dashed line) and 50 GeV (green dashed line) respectively. However for $M_{DM}=75$ GeV (blue dashed dot line) $\Omega h^{2}$ does not rise equally because for this value of DM mass the decay of $h_{1}$ to a pair of $\phi_{DM}$ and $\phi_{DM}^{\dagger}$ becomes kinematically forbidden and hence, there is no equal increment in this case.

In LP and RP of Fig. 10, we have shown how the relic density varies with $z$ for different values of scalar quartic couplings $\lambda_{Dh}$ and $\lambda_{DH}$, respectively. In each panel, one can easily notice that there exists a kink around the EWSB region. However in the LP, the kink occurs for a higher value of $\lambda_{Dh}$ while in the RP, the situation is just opposite. We have already seen in the LP of Fig. 9 that before EWSB only $h_{2}$ decay is contributing to DM relic density and at the EWSB region SM-like Higgs $h_{1}$ also starts contributing. A kink will always appear in the relic density curve when contribution of the SM-like Higgs boson $h_{1}$ to $\Omega h^{2}$ is larger compared to that of the BSM Higgs $h_{2}$ and extra gauge boson $Z_{BL}$ i.e. $\Gamma_{h_{1}\rightarrow\phi_{DM}\phi_{DM}^{\dagger}}>\Gamma_{h_{2}\rightarrow\phi_{DM}\phi_{DM}^{\dagger}}$, $\Gamma_{Z_{BL}\rightarrow\phi_{DM}\phi_{DM}^{\dagger}}$. The values of scalar quartic couplings $\lambda_{Dh}$ and $\lambda_{DH}$ in the LP of Fig.  9 are such that $\Gamma_{h_{2}\rightarrow\phi_{DM}\phi_{DM}^{\dagger}}$ and $\Gamma_{Z_{BL}\rightarrow\phi_{DM}\phi_{DM}^{\dagger}}$ always remain large compared to $\Gamma_{h_{1}\rightarrow\phi_{DM}\phi_{DM}}$ and hence no kink is observed in the total relic density curve. However, in the present figure (in the left panel of Fig.  10) we do have kinks around the EWSB region, because in the LP with $\lambda_{DH}=8.316\times 10^{-14}$ and $n_{BL}=1.33\times 10^{-10}$, $\Gamma_{h_{1}\rightarrow\phi_{DM}\phi_{DM}^{\dagger}}>\Gamma_{h_{2}\rightarrow\phi_{DM}\phi_{DM}^{\dagger}}$, $\Gamma_{Z_{BL}\rightarrow\phi_{DM}\phi_{DM}^{\dagger}}$ condition is satisfied only for the case with larger value of $\lambda_{Dh}=1.237\times 10^{-11}$ ($\lambda_{Dh}>>\lambda_{DH}$) while in the RP with a fixed value of $\lambda_{Dh}=1.237\times 10^{-12}$, the above condition is not maintained because $Z_{BL}$ decay channel dominates.

In the LP of Fig. 11, we have shown the allowed region in the coupling plane ($\lambda_{Dh}-\lambda_{DH}$) which reproduces the observed DM relic density ($0.1172\leq\Omega h^{2}\leq 0.1226$). In this figure, we have clearly indicated the dominant DM production processes when $M_{DM}$ varies between 10 GeV to 100 GeV i.e.  DM production from the decays of $h_{1}$, $h_{2}$ or both or entirely from the annihilations of SM particles like $W^{\pm},\,\,Z,\,\,h_{1}$ etc. The parameters which are related to the $Z_{BL}$ decay ($g_{BL}$, $n_{BL}$) have been kept fixed at 0.07 and $1.33\times 10^{-10}$ respectively, so at every time an equal amount of $Z_{BL}$ decay contribution remains present. As illustrated in the figure, when the parameter $\lambda_{Dh}$ is small compared to the other parameter $\lambda_{DH}$ then among the two scalars it is the BSM Higgs $h_{2}$ which is mainly contributing to the DM production while for the opposite case, the production of $\phi_{DM}$ becomes $h_{1}$ dominated and in between both the scalars contribute equally. Apart from that, if the mass of $\phi_{DM}$ is greater than the half of the SM-like Higgs mass (i.e. $M_{DM}>\dfrac{M_{h_{1}}}{2}$) then DM production from $h_{1}$ decay becomes kinematically forbidden. In this case, however, the production from the decays of $h_{2}$ and $Z_{BL}$ are still possible. Now, the deficit in DM production can be compensated by the production from self annihilation of SM particles like $h_{1}$, $W^{\pm}$ and $Z$ and for this we need to increase the parameter $\lambda_{Dh}$. Moreover, by increasing $\lambda_{Dh}$ (decreasing $\lambda_{DH}$ simultaneously) we can arrive a situation where DM production is entirely dominated by the annihilations of SM particles and this situation has been indicated by a pink coloured arrow in the LP of Fig. 11. On the other hand, in the RP of Fig. 11 we have presented the allowed region in $M_{h_{2}}-\alpha$ plane which satisfies the relic density bound. From this figure one can see that with the increase of $M_{h_{2}}$, the allowed values of mixing angle $\alpha$ decrease. The reason behind this decrement is related to the vacuum stability conditions as given in the Eq. (17). The region satisfying both the relic density bound as well as the vacuum stability conditions is shown by the green dots while in the other part of $M_{h_{2}}-\alpha$ plane the quantity $\mu_{\phi_{h}}^{2}$ becomes positive which is undesirable in the context of the present model (see Eq. (17)).

III.3.2 $\dfrac{M_{h_{1}}}{2}<M_{\phi_{DM}}<\dfrac{M_{h_{2}}}{2},\,\dfrac{M_{Z_{BL}}}{2}$,

BSM particles decay and SM particles annihilation dominated region.

Clearly in this mass region, DM production from the decay of SM-like Higgs $h_{1}$ is kinematically forbidden and hence DM has been produced from the decays of $h_{2}$, $Z_{BL}$ only. However, unlike the previous case, here we find significant contribution to DM relic density arising from the self annihilation of the SM particles namely, $h_{1}$, $W^{\pm}$, $Z$ and $t$. On the other hand, the annihilations of BSM particles like $Z_{BL}$, $h_{2}$ and $N_{i}$ have negligible effect on DM production processes. In Fig. 12, we have shown the variation of DM relic density with $z$ for $\dfrac{M_{h_{1}}}{2}<M_{\phi_{DM}}<\dfrac{M_{h_{2}}}{2},\,\dfrac{M_{Z_{BL}}}{2}$. Since now the decay of the $h_{1}$ to $\phi_{DM}{\phi_{DM}}^{\dagger}$ is kinematically forbidden, hence we can increase the parameter $\lambda_{Dh}$ safely and this will not overproduce DM in the Universe. Due to this moderately large value of $\lambda_{Dh}$, the annihilation channels become important. From Fig. 12 it is clearly seen that in this case the annihilation channel $h_{1}h_{1}\rightarrow\phi_{DM}\phi_{DM}^{\dagger}$ (Green dashed line) contributes significantly to the DM production. Therefore in the present case, production of DM has been controlled by the decays of $h_{2}$, $Z_{BL}$ and the self annihilations of the SM particles and thus directly relates to the ${\rm U}(1)_{\rm B-L}$ sector of this model.

III.3.3 $\dfrac{M_{h_{1}}}{2},\,\dfrac{M_{h_{2}}}{2}<M_{\phi_{DM}}<\dfrac{M_{Z_{BL}}}{2}$,

BSM particles decay and annihilation dominated region.

In this regime of the DM mass, the only surviving decay mode is the decay of ${\rm B-L}$ gauge boson $Z_{BL}$ to a pair of $\phi_{DM}$. Apart from that, depending on the choice of mass of $\phi_{DM}$ a significant fraction of DM has been produced from the self annihilation of either BSM Higgs $h_{2}$. In other word, we can say that in this region the production of DM is BSM particles dominated. In LP of Fig. 13 we show the relative contribution of dominant production modes of DM to $\Omega h^{2}$ for a chosen value of $M_{DM}=450$ GeV. From this plot one can easily notice that in the case when $M_{DM}<M_{h_{2}}$, the almost entire fraction of DM is produced from the decay of $Z_{BL}$ (green dashed line) and self annihilation of BSM Higgs $h_{2}$ (solid turquoise line). This is because, as in this case the production of $\phi_{DM}$ from $h_{2}$ decay is kinematically forbidden hence one can increase the parameter $\lambda_{DH}$ so that the annihilation channel $h_{2}h_{2}\rightarrow\phi_{DM}{\phi_{DM}}^{\dagger}$, which is mainly proportional to $\lambda_{DH}^{2}$ (due to four point interaction) becomes significant.

On the other hand, in the RP we have considered a situation where almost the entire DM has been produced from the decay of ${\rm B-L}$ gauge boson. For this, we have chosen $M_{DM}>M_{h_{2}}$ and a larger value of $n_{BL}=6.2\times 10^{-11}$. Similar to the previous case (i.e. $M_{DM}<M_{h_{2}}$) here also, the production of $\phi_{DM}$ from $h_{2}$ decay still remains forbidden. However, as the sum of final state particles masses are larger than that of initial state hence, in this case $h_{2}h_{2}$ annihilation mode becomes suppressed. Moreover, to make the contribution of $h_{2}$ annihilation even more suppressed we have reduced the quartic couplings $\lambda_{Dh}$ and $\lambda_{DH}$ accordingly. As a result other annihilation channels e.g. $Z_{BL}Z_{BL}$, $N_{i}N_{i}$ also become inadequate as these channels are mediated by the exchange of $h_{1}$ and $h_{2}$. Although, RH neutrinos can annihilate to $\phi_{DM}{\phi_{DM}}^{\dagger}$ via $Z_{BL}$, we cannot increase the contribution of $Z_{BL}$ mediated diagrams because for that one has to further increase the ${\rm B-L}$ charge of $\phi_{DM}$ ($n_{BL}$), which results in an over production of DM in the Universe from $Z_{BL}$ decay. From the RP of Fig. 13, one can easily notice that in this situation $Z_{BL}$ decay is the most dominant DM production channel (red dashed line) while the total contributions from the annihilations of $h_{2}$, $Z_{BL}$ and $N_{i}$ are negligible. Therefore, for the entire mass range of $\phi_{DM}$ i.e. $\dfrac{M_{h_{1}}}{2},\,\dfrac{M_{h_{2}}}{2}<M_{\phi_{DM}}<\dfrac{M_{Z_{BL}}}{2}$, the DM production processes are always related to the ${\rm U}(1)_{\rm B-L}$ sector of the present model by receiving a sizeable contribution from $Z_{BL}$ decay.

In Fig.  14, we have shown the allowed region (green coloured points) in $M_{Z_{BL}}-g_{BL}$ plane which reproduces the observed DM relic density. While generating this plot we have varied $250\,{\rm GeV}\leq M_{DM}\leq$ 5000 GeV and $10^{-11}\leq n_{BL}\leq 10^{-8}$. In this region as mentioned above dominant contributions to DM relic density arise from $Z_{BL}$ decay and annihilation of BSM Higgs $h_{2}$. In this figure, the black solid line represents the current upper bound Chatrchyan:2012oaa ; Aad:2014cka ; Guo:2015lxa on $g_{BL}$ for a particular mass of $Z_{BL}$ from LHC 777To get the bound in $M_{Z_{BL}}-g_{BL}$ plane from LHC, ATLAS and CMS collaborations consider the Drell-Yan processes ($p\,p\rightarrow Z_{BL}\rightarrow\bar{l}\,l$, with $l$ = e or $\mu$) and by searching the dilepton resonance they put lower bound on $M_{Z_{BL}}$ for a particular value of extra gauge coupling $g_{BL}$. while the limit Carena:2004xs ; Cacciapaglia:2006pk ; Schael:2013ita from LEP 888LEP consider the processes $e^{+}\,e^{-}\rightarrow\bar{f}\,f$ ($f\neq e$) above the Z-pole mass and by measuring its cross section they put lower limit on the ratio between the gauge boson mass and guage coupling, which is $\frac{M_{Z_{BL}}}{g_{BL}}\geq 6-7$ TeV. has been indicated by the red solid line respectively. Therefore, the region below the red and black solid line is allowed by the collider experiments like LHC and LEP. The benchmark value of $g_{BL}$, $M_{Z_{BL}}$ (= 0.07, 3000 GeV) for which we have computed the baryon asymmetry in the previous section (Section III.2) is highlighted by a blue coloured star. Hence, in this regime the extra gauge boson $Z_{BL}$ immensely takes part in achieving the correct ballpark value of the DM relic density and also at the same time $Z_{BL}$ plays a significant role to obtain the observed value of the matter-antimatter asymmetry of the Universe.

III.3.4 $M_{\phi_{DM}}>\dfrac{M_{h_{1}}}{2},\,\dfrac{M_{h_{2}}}{2},\,\dfrac{M_{Z_{BL}}}{2}$,

BSM particles annihilation dominated region.

Finally, in this range of DM mass the entire production of $\phi_{DM}$ from the decays of $h_{1}$, $h_{2}$ and $Z_{Bl}$ become kinematically inaccessible. Therefore, in this case all three parameters namely $\lambda_{Dh}$, $\lambda_{DH}$ and $n_{BL}$ become free and we can make sufficient increment to these parameters so that either scalar medicated ($h_{1}$, $h_{2}$) or gauge boson mediated ($Z_{BL}$) annihilation processes of $N_{i}$, $Z_{BL}$ or both can be the dominant contributors in DM production.

Similarly, in the LP and RP of Fig. 15, we have shown two different situations where the DM production are dominated by scalar ($h_{1}$, $h_{2}$) mediated diagrams and gauge boson $Z_{BL}$ mediated diagrams respectively. In the LP, by keeping the $n_{BL}$ value low and adjusting the parameters $\lambda_{Dh}$ and $\lambda_{DH}$ one can achieve the correct value DM relic density and on the other hand, in the RP we have kept the values of $\lambda_{Dh}$ and $\lambda_{DH}$ sufficiently low and by suitably adjusting the DM charge $n_{BL}$ we have achieved the correct value of the DM relic density. Therefore, in this region, a strong correlation exists among the neutrino sector, ${\rm U}(1)_{\rm B-L}$ sector and DM sector as the entire DM is now being produced from $N_{i}N_{i}$ and $Z_{BL}Z_{BL}$ annihilations.

In Fig. 16, we have shown the allowed parameter space in $M_{DM}-M_{N_{1}}$ plane by DM relic density. In order to generate this plot we have varied DM mass in the range $1500\,{\rm GeV}\leq M_{DM}\leq 3000\,{\rm GeV}$, RH neutrino masses $1500\,{\rm GeV}\leq M_{N_{i}}\leq 10000\,{\rm GeV}$ ($i=1,\,2$), $M_{N_{1}}<M_{N_{3}}\leq M_{N_{1}}+5000$ GeV and $10^{-10}\leq n_{BL}\leq 10^{-8}$. Other relevant parameters have been kept fixed at $\lambda_{Dh}=7.212\times 10^{-13}$, $\lambda_{DH}=8.316\times 10^{-12}$, $M_{Z_{BL}}$ = 3000 GeV, $g_{BL}$ = 0.07, $M_{h_{2}}$ = 500 GeV, $\alpha=10^{-5}$ As discussed above, in this regime ($M_{DM}>\dfrac{M_{h_{1}}}{2},\,\dfrac{M_{h_{2}}}{2},\,\dfrac{M_{Z_{BL}}}{2}$) $\phi_{DM}$ is dominantly produced from the annihilations of $Z_{BL}$ and RH neutrinos. From this plot one can observe that in this high DM mass range to obtain the observed DM relic density, the mass of the lightest RH neutrino cannot be larger than $\sim 6000$ GeV. Analogous to the Fig. 14, here also we have indicated the benchmark point for which we have computed baryon asymmetry in the previous section (Section III.2) by a blue coloured star. Therefore, in this case RH neutrinos are very actively taking part in all three processes we have considered in this work namely DM production processes, tiny neutrino mass generation and also the generation of required lepton asymmetry to reproduce the observed baryon asymmetry of the Universe.

From the above four regions, which are based on the mass of our FIMP DM, it is evident that in the first region DM production mainly happens from the decay of $h_{1}$, $h_{2}$ and $Z_{BL}$ and all annihilations are subdominant. Therefore, in this region only the extra neutral gauge boson ($Z_{BL}$), BSM Higgs ($h_{2}$) and SM-like Higgs ($h_{1}$) are taking part in the DM relic density estimates and there is no significant role of the RH neutrinos. In the second region, SM-like Higgs decay does not contribute to DM production processes, hence one can safely increase the quartic coupling $\lambda_{Dh}$ and consequently $h_{1}h_{1}$ annihilation contribution increases. Similar to the previous regime, here also RH neutrinos have less importance in determining the DM relic density. In the third region, the only decay mode that involves in DM production is $Z_{BL}\rightarrow\phi_{DM}{\phi_{DM}}^{\dagger}$. Since all other decay modes correspond to $h_{1}$ and $h_{2}$ are kinematically forbidden, hence we can increase both the quartic couplings $\lambda_{Dh}$ and $\lambda_{DH}$ appropriately which eventually enhance the annihilation contribution from the BSM Higgs significantly. Moreover, due to the increment of quartic couplings in this region $Z_{BL}Z_{BL}$ and $N_{i}N_{i}$ annihilation channels start contributing in the DM production processes. Lastly in region four, due to the high value of the DM mass no decay process contributes to DM relic density and only the BSM particles annihilation contributes. Therefore, in this region by properly adjusting the extra gauge coupling $g_{BL}$, one can get a sizeable fraction of DM production from the annihilation of RH neutrinos. Since apart from the masses of the involving particles, the annihilation of RH neutrinos mediated by $Z_{BL}$ depends on the extra (B-L) gauge coupling $g_{BL}$ solely. Thus, depending on the mass range of our FIMP DM, we can say that the different model parameters and the additional BSM particles (e.g. $Z_{BL}$, $N_{i}$, $h_{2}$) are fully associated to the DM production processes in the early Universe.

III.3.5 Analytical Estimates

So far, we have solved the full Boltzmann equation (Eq. 50) for a FIMP $\phi_{DM}$ numerically. Apart from this, one can estimate the FIMP relic density (or comoving number density) by using the approximate analytical formula. Let us consider a FIMP ($\phi_{DM}$) which is produced from the decay of a particle $A$ i.e., $A\rightarrow\phi_{DM}\,{\phi_{DM}}^{\dagger}$, where $A$ in the present model can be $h_{1}$, $h_{2}$ or $Z_{BL}$. The contribution of $A$ to the FIMP relic density at the present epoch, considering the effect of both $\phi_{DM}$ and ${\phi_{DM}}^{\dagger}$, is given by Hall:2009bx ,

[TABLE]

where $M_{A}$ and $g_{A}$ are mass and internal degrees of freedom of the mother particle $A$, respectively, while $\Gamma_{A}$ is the decay width of the process $A\rightarrow\phi_{DM}\,{\phi_{DM}}^{\dagger}$. The analytic expressions for $\Gamma_{A}$ corresponding to $h_{2}$, $h_{1}$ and $Z_{BL}$ are given in Eqs. (76), (79) and (81) in the Appendices. Moreover, $g_{\rho}$ and $g_{\rm s}$, as define earlier, are the degrees of freedom related to the energy and entropy densities of the Universe, respectively. Let us now compare the analytical result with the numerical value which we obtain by solving the Boltzmann equation Eq. (50). For this, let us consider a situation when a significant fraction of our FIMP candidate ($\phi_{DM}$) is produced from the decay mode of BSM Higgs i.e., $h_{2}\rightarrow\phi_{DM}\phi_{DM}^{\dagger}$. Substituting the values of model parameters given in the caption of Fig. 9 to Eq. (53), we get the contribution of $h_{2}$ to DM relic density, which is

[TABLE]

where we consider $g_{\rho}=g_{\rm s}\approx 100$ and $g_{A}=1$. This can be compared to the contribution of $h_{2}$ obtained from exact numerical estimate shown in the LP of Fig. 9 which is,

[TABLE]

Therefore, from the above two estimates it is clearly evident that the analytical result agrees well with the full numerical result. Similarly, for the other decay modes also (i.e. $h_{1}$, $Z_{BL}$) one can match the analytical and numerical results.

IV Conclusion

In this work we considered a local ${\rm U}(1)_{\rm B-L}$ extension of the SM and to cancel the additional anomalies associated with this gauge symmetry we introduced three RH neutrinos ($N_{i}$, $i$ = 1, 2, 3). Besides the three RH neutrinos, we also introduced two SM gauge singlet scalars $\phi_{H}$ and $\phi_{DM}$. The scalar field $\phi_{H}$, being charged under ${\rm U}(1)_{\rm B-L}$, takes a nonzero VEV and breaks the proposed ${\rm B-L}$ symmetry spontaneously. Moreover, as the scalar field $\phi_{DM}$ has also a nonzero ${\rm B-L}$ charge, one can adjust this charge suitably so that after symmetry breaking the model has left with a residual $\mathbb{Z}_{2}$ symmetry and only $\phi_{DM}$ behaves as a odd particle under this leftover symmetry. This makes $\phi_{DM}$ absolutely stable over the cosmological time scale and hence acts as a dark matter candidate. After spontaneous breaking of ${\rm U}(1)_{\rm B-L}$ gauge symmetry, all RH neutrinos and extra neutral gauge boson $Z_{BL}$, acquired mass. Due to the presence of the three RH neutrinos in the model, we easily generated Majorana masses for the three light neutrinos by the Type I seesaw mechanism. This model is also able to explain baryogenesis via leptogenesis, where we generated the lepton asymmetry in the Universe from out of equilibrium, CP violating decays of two degenerate RH neutrinos and converted this lepton asymmetry to the observed baryon asymmetry through the sphaleron transitions.

In explaining the neutrino masses by Type I seesaw mechanism, we considered a complex Dirac mass matrix $\mathcal{M_{D}}$ and a diagonal Majorana mass matrix $\mathcal{M_{R}}$ for the RH neutrinos. In determining the allowed model parameter space, we used the measured values of neutrino oscillation parameters namely three mixing angles ($\theta_{12}$, $\theta_{13}$ and $\theta_{23}$) and two mass square differences ($\Delta m_{21}^{2}$, $\Delta m_{atm}^{2}$) in their current $3\sigma$ range. In particular, in the current model we could reproduce the whole allowed $3\sigma$ range of the neutrino oscillation parameters by different combinations of the relevant model parameters. The Dirac CP phase was constrained to lie within two distinct regions. One is the entire first quadrant ($0^{\circ}-90^{\circ}$) while the other one spans the entire fourth quadrant ($270^{\circ}-360^{\circ}$). However, if we considered the T2K result on Dirac CP phase then the values of $\delta$ lying in the fourth quadrant are more favourable compared to those in the first quadrant. We also computed the magnitudes of the Jarlskog invariant $J_{\rm CP}$ and found that the values of $J_{\rm CP}$, for the model parameters which satisfy neutrino oscillation data, always lie below $0.039$. Finally, we calculated the values of $m_{\beta\beta}$, the quantity relevant to neutrino less double $\beta$ decay, for the allowed model parameter space.

Since we allowed complex Yukawa couplings in the Dirac mass matrix $\mathcal{M_{D}}$, the decays of RH neutrinos were CP violating. We took the masses of the RH neutrinos in the TeV range and worked in the parameter space where the lightest two RH neutrino states were nearly degenerate, with their masses separated by their tree level decay width. This scenario led to resonant leptogenesis (or TeV scale leptogenesis) for the production of observed baryon asymmetry in the Universe from the out of equilibrium decays of RH neutrinos. We generated the observed baryon asymmetry for three different values of RH neutrino masses namely $M_{N_{1}}=1600$ GeV, 1800 GeV and 2000 GeV, respectively, where required values of CP asymmetry parameter parameter ($\varepsilon_{1}$) were $4.4\times 10^{-4}$, $2.25\times 10^{-4}$ and $1.8\times 10^{-4}$, respectively. These values of $M_{N_{1}}$ and $\epsilon_{1}$ were also seen to be allowed by the neutrino oscillation data.

Lastly, we studied the DM phenomenology by considering a FIMP type DM candidate $\phi_{DM}$. We took into account all the production modes of $\phi_{DM}$ (both before and after EWSB) arising from the annihilations and decays of SM as well as BSM particles. We found that depending on the mass of $\phi_{DM}$, the production processes of $\phi_{DM}$ can be classified into four distinct categories. These are (1) SM and BSM particles decay dominated region, (2) BSM particles decay and SM particles annihilation dominated region, (3) BSM particles annihilation and $Z_{BL}$ decay dominated region and finally (4) BSM particles annihilation dominated region. The first region is characterised by $M_{DM}<\dfrac{M_{h_{1}}}{2},\,\dfrac{M_{h_{2}}}{2},\,\dfrac{M_{Z_{BL}}}{2}$ and here DM is mainly produced from the decays of $h_{1}$, $h_{2}$ and $Z_{BL}$. In the second region, DM mass is concentrated between $\dfrac{M_{h_{1}}}{2}$ and min$\left[\dfrac{M_{h_{2}}}{2},\dfrac{M_{Z_{BL}}}{2}\right]$ i.e. $\dfrac{M_{h_{1}}}{2}<M_{DM}<\dfrac{M_{h_{2}}}{2},\,\dfrac{M_{Z_{BL}}}{2}$. In this case, $h_{2}$, $Z_{BL}$ decays and $h_{1}h_{1}$, $W^{+}W^{-}$, $ZZ$ annihilations act as the dominant production modes of FIMP DM. In the third region where $\dfrac{M_{h_{1}}}{2},\,\dfrac{M_{h_{2}}}{2}<M_{DM}<\dfrac{M_{Z_{BL}}}{2}$, DM has mainly been produced from $Z_{BL}$ decay and also from the annihilation of BSM Higgs $h_{2}$ (for $M_{DM}<M_{h_{2}}$). Finally, in the last region all three decay modes become kinematically forbidden as $M_{DM}>\dfrac{M_{h_{1}}}{2},\,\dfrac{M_{h_{2}}}{2},\,\dfrac{M_{Z_{BL}}}{2}$ and hence entire DM is produced from the self annihilations of $Z_{BL}$ and right handed neutrinos ($N_{i}$). Therefore in all four regions the ${\rm U}(1)_{\rm B-L}$ gauge boson has played a significant role in DM production while the effects of right handed neutrinos are important in the last two regions only. We also found that, since for a FIMP candidate $\phi_{DM}$ the observed DM relic density ($0.1172\leq\Omega_{h^{2}}\leq 0.1226$) is generated via the freeze-in mechanism, this puts upper bounds on the scalar and gauge portal couplings of $\phi_{DM}$ to restrict its over production i.e. \lambda_{Dh}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}10^{-11}, \lambda_{DH}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}10^{-10} and n_{BL}\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}10^{-8}. Hence, due to such extremely feeble couplings $\phi_{DM}$ can easily evade all the constrains coming from any terrestrial DM direct detection experiment.

In conclusion, our spontaneously broken local ${\rm U}(1)_{\rm B-L}$ extension of the SM with three additional RH and two additional scalars can explain the three main evidences for physics beyond the SM, viz., small neutrino masses, matter-antimatter asymmetry of the Universe and dark matter. Tiny neutrino masses and all mixing angles can be obtained via Type I seesaw mechanism where we chose a certain pattern for the real and complex Yukawa couplings. The model gave a definite prediction for the CP violating phase to be measured in the next generation long baseline experiments. The dark matter candidate is a scalar which is neutral under the SM gauge group and has a nonzero ${\rm B-L}$ charge. DM is made stable by virtue of a remnant $\mathbb{Z}_{2}$ symmetry arises after the spontaneous breaking of ${\rm U}(1)_{\rm B-L}$ gauge symmetry. This can be achieved by imposing a suitable ${\rm B-L}$ charge on $\phi_{DM}$ so that the Lagrangian does not contain any odd term of $\phi_{DM}$. This scalar DM can easily be taken as a FIMP candidate which is produced from the decays and annihilations of SM and BSM particles. Therefore, even if the WIMP type DM is ruled out in near future from direct detection experiments this present variant of ${\rm U}(1)_{\rm B-L}$ scenario with FIMP DM will still survive. Further, since $g_{BL}$ is of the order of SM gauge couplings, this model has the potential to be tested in the LHC or in other future collider experiments by detecting ${\rm B-L}$ gauge boson $Z_{BL}$ from its SM decay products. Moreover, considering the masses of RH neutrinos in TeV scale allow us to simultaneously explain the baryon asymmetry of the Universe from resonant leptogenesis, FIMP DM production via Freeze-in mechanism and also neutrino masses and mixing from TeV scale Type-I seesaw. Thus, all three phenomena addressing in this article are interconnected to each other.

V Acknowledgement

SK and AB also acknowledge the HRI cluster computing facility (http://cluster.hri.res.in). The authors would also like to thank the Department of Atomic Energy (DAE) Neutrino Project under the XII plan of Harish-Chandra Research Institute. 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. This project has received funding from the European Union’s Horizon 2020 research and innovation programme Elusives ITN under the Marie Sklodowska-Curie grant agreement No 674896.

Appendix A Expression for the Majorana mass matrix of light neutrinos

Here we have given the expression of all the elements of the light neutrino mass matrix $m_{\nu}$ in terms of the Yukawa couplings and the RH neutrino masses.

[TABLE]

Appendix B Neutrino less double $\beta$ decay parameter $m_{\beta\beta}$

Since the light neutrino mass matrix is Majorana in nature, it is a complex symmetric matrix. A complex symmetric matrix $m_{\nu}$ can be diagonalised by a Unitary matrix $U_{\rm PMNS}$ (defined in Eq. (12)) in the following way,

[TABLE]

Now equating the $(i,\,j)^{\rm th}$ element from both sides of the above equation, we get

[TABLE]

Since $m_{\rm dia}$ is a diagonal matrix, we can further simplify $({m_{\nu}})_{ij}$ by using $m_{\rm diag}=m_{k}\,\delta_{kk^{\prime}}$, where $m_{k}$ is the mass of $k^{\rm th}$ light neutrino. Therefore $(m_{\nu})_{ij}$ takes the following form

[TABLE]

Above equation expresses the elements of light neutrino mass matrix (Eq. (56)) in terms of the light neutrino masses, the intergenerational mixing angles and the phases. Taking $i=j=1$, we get the expression of the (1,1) element of $m_{\nu}$ i.e.

[TABLE]

which is related to the important parameter $m_{\beta\beta}$ of the neutrino less double $\beta$ decay DellOro:2016tmg as,

[TABLE]

Appendix C CP Asymmetric Parameter Calculation for Leptogenesis

The amount of lepton asymmetry generated in the out of equilibrium decay of the RH neutrino $N_{i}$ is parametrised by the CP asymmetry parameter ($\varepsilon_{i}$), which is defined as,

[TABLE]

If we consider only the tree level decay process of $N_{i}$ (first diagram in Fig.  2), there will not be any CP violation. The nonzero CP asymmetry is generated only by the interference between the tree level and the one loop level diagrams. The expression of CP asymmetry parameter ($\varepsilon_{i}$) is given by Flanz:1996fb ; Pilaftsis:1997jf ; Pilaftsis:2003gt ,

[TABLE]

where $V_{j}$ and $S_{j}$ are the contributions coming from the vertex correction and the self energy correction diagrams respectively (second and third diagrams in Fig. 2). The expressions of $V_{j}$ and $S_{j}$ have the following forms Flanz:1996fb ; Pilaftsis:1997jf ; Pilaftsis:2003gt ,

[TABLE]

with

[TABLE]

and $\Gamma_{j}$ denotes the tree level decay width of the RH neutrino $N_{j}$ (neglecting subdominant one loop corrections), which is given by

[TABLE]

Now, as mentioned in the beginning of this section, the enhancement in the CP asymmetry factor (Eq. (62)) occurs when two RH neutrinos are almost degenerate i.e. $M_{N_{2}}-M_{N_{1}}\simeq\dfrac{\Gamma_{1}}{2}$. This is know as the resonance condition. In the present scenario, we consider $M_{N_{3}}>M_{N_{2}}\simeq M_{N_{1}}$. Therefore, resonance condition is satisfied only for the two lightest RH neutrinos $N_{2}$ and $N_{1}$. Hence we can neglect the contribution of $N_{3}$ in the CP asymmetry parameter (Eq. (63)) by considering the summation over only $N_{1}$ and $N_{2}$ (i.e. $j=1,\,2$). Using the resonance condition in Eq. (64), one can easily notice that $S_{j}\sim\mathcal{O}\left(\frac{M_{N_{j}}}{\Gamma_{j}}\right)>>1$ ($j=1,\,2)$ and

[TABLE]

where we have neglected the quantity $V_{j}$ which is, in the present condition ($M_{N_{2}}-M_{N_{1}}\simeq\dfrac{\Gamma_{1}}{2}$), much smaller compared to $S_{j}$. The resonance condition leads to

[TABLE]

Using Eq. (68) in Eq. (64) we get,

[TABLE]

Now, substituting the expressions of $S_{1}$ and $S_{2}$ in Eq. (67) and using ${\rm Im}\left[(\mathcal{M_{D}}\mathcal{M_{D}}^{\dagger})^{2}_{12}\right]=-\,{\rm Im}\left[(\mathcal{M_{D}}\mathcal{M_{D}}^{\dagger})^{2}_{21}\right]$, one obtains,

[TABLE]

Appendix D Expressions of decay widths of $h_{2}$, $h_{1}$ and $Z_{BL}$

In the present work, we have considered the effect of electroweak symmetry breaking on dark matter production. After EWSB SM particles become massive and affect DM production, while before EWSB those particles have no effect. To take this effect into account we have defined an extra constant $C_{ASB}$. In all the equations, the value of the constant $C_{ASB}=0$ before the EWSB and this is equal to unity i.e. $C_{ASB}=1$ after the EWSB. Also before the EWSB, there is no mixing between the SM and BSM Higgs bosons, i.e. $\alpha=0$. The two vertices which are common to all Higgs mediated diagrams are as follows,

[TABLE]

Total decay width of $h_{2}$:

• •

$h_{2}$ $\rightarrow VV$** ($V=W^{\pm},Z$): **

[TABLE]

where $S_{V}=2\,(1)$ for $ZZ$ ($W^{+}W^{-}$) final state.

• •

**$h_{2}$ $\rightarrow$ $h_{1}$ $h_{1}$****: **

[TABLE]

• •

**$h_{2}$ $\rightarrow$ $\phi_{DM}^{\dagger}$ $\phi_{DM}$****: **

[TABLE]

• •

**$h_{2}$ $\rightarrow$ $f\bar{f}$****: **

[TABLE]

$n_{c}$** is the color charge, for leptons it is $1$ and for quarks it is $3$.**

Total decay width of the extra Higgs $h_{2}$ in the present case is,

[TABLE]

Total decay width of $h_{1}$****:

[TABLE]

Total decay width of SM-like Higgs boson is,

[TABLE]

where $\Gamma_{\rm SM}$ is the total decay width of SM Higgs boson.

Total decay width of $Z_{BL}$****:

[TABLE]

Total decay width of the extra neutral gauge boson $Z_{BL}$ is,

[TABLE]

Appendix E Analytical Expression of relevant Cross sections

Here we will give the expressions of cross sections for all relevant processes which take part in the FIMP DM production.

• •

$h_{1}$ $h_{1}$ $\rightarrow\phi_{DM}^{\dagger}\,\phi_{DM}$** :**

[TABLE]

• •

$h_{2}$ $h_{2}$ $\rightarrow\phi_{DM}^{\dagger}\,\phi_{DM}$** :**

[TABLE]

• •

$h_{1}$ $h_{2}$ $\rightarrow\phi_{DM}^{\dagger}\,\phi_{DM}$** :**

[TABLE]

• •

$W^{+}$ $W^{-}$ $\rightarrow\phi_{DM}^{\dagger}\,\phi_{DM}$** : **

[TABLE]

• •

$Z\,Z\rightarrow\phi_{DM}^{\dagger}\,\phi_{DM}$** : **

[TABLE]

• •

$t\bar{t}$  $\rightarrow\phi_{DM}^{\dagger}\,\phi_{DM}$** :**

[TABLE]

• •

$N_{i}$ $N_{i}$ $\rightarrow\phi_{DM}^{\dagger}\,\phi_{DM}\,\,(i=1,\,2,\,3)$** :**

[TABLE]

• •

$Z_{BL}\,Z_{BL}\rightarrow\phi_{DM}^{\dagger}\,\phi_{DM}$** : **

[TABLE]