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

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


Introduction
The standard model (SM) is a very successful theory in describing nature. The discovery of the last missing piece of the SM, viz., the Higgs boson, further increases its concreteness. In spite of its tremendous success, the SM can not explain a number of phenomena -two of the most important ones being the presence of dark matter (DM) and non-zero neutrino mass. Presence of DM in the universe is a very well established fact. The first indication of DM came from the observation of Galactic velocities within the Coma cluster by Fritz Zwicky in 1933 [1], followed by the observation of galaxy rotation curves by Vera Rubin in 1970 [2]. Subsequently, the observation of bullet cluster [3] firmly confirmed the presence of DM. Currently the best measurement of the amount of DM present in the universe comes from the Planck data [4], Ω h 2 = 0.1199 ± 0.0027 at 68% CL , (1.1) where h is the reduced Hubble parameter and of order unity. Unfortunately, the SM does not have any fundamental particle which can be a viable DM candidate. Therefore, to address the issue of DM from particle physics point of view, we need to extend the SM particle content and/or its gauge group. One of the most promising scenarios is to consider the DM candidate as a Weakly Interacting Massive Particle (WIMP) [5,6], which is produced in the early universe through the thermal freeze-out mechanism [5,6]. However, WIMP type DM attracts stringent bounds from direct and indirect detection experiments [7][8][9][10][11][12][13][14].
In particular, a large portion of the parameter space in the spin independent/dependent WIMP-nucleon cross section and DM mass plane is ruled out by the direct detection (DD) bounds. Moreover, in near future with increasing sensitivity of the DD experiments [7][8][9][10][11], these bounds might touch the so-called neutrino floor [15,16]. In this work, we follow a non-thermal way of DM production, viz., via the freeze-in mechanism [17]. In this scenario,

JHEP06(2019)095
the DM is very feebly interacting with the other particles, and as a result never achieves thermal equilibrium in the early universe with the cosmic soup. Hence it is named Feebly Interacting Massive Particles (FIMPs). Due to their very feeble interactions, FIMPs easily escape the above mentioned DD bounds while satisfying the measured value for the DM relic density [17][18][19][20][21][22][23][24][25].
On the other hand, results of the neutrino oscillation experiments [26][27][28][29][30][31][32][33][34][35][36] have confirmed oscillations between neutrino flavours. Since neutrino flavour oscillations are a clear proof of the neutrinos being massive and mixed, the neutrino oscillation experiments contradict the SM which postulates that the neutrinos are massless. Consequently, in order to explain tiny neutrino masses, one has to extend the SM by adding new particles and/or additional gauge groups.
In the present work we explain the above two puzzles by extending the SM gauge group by a U(1) B−L gauge symmetry as a simple (minimal) and well motivated extension of the SM, where B is the baryon number and L is the lepton number. In addition to the extra neutral gauge boson Z associated with the U(1) B−L , an extra SM singlet scalar φ H (charged under U(1) B−L to break B − L gauge symmetry spontaneously) is added in this simple extension, which leads to interesting signatures at the LHC [37][38][39][40][41]. Moreover, nine additional SM singlet fermions (N i R and S i 1,2 , i = 1, 2, 3) are needed to explain the naturally 1 small neutrino masses through the inverse seesaw mechanism [49][50][51][52]. These additional fermions are not only required to generate the tiny neutrino masses via the inverse seesaw mechanism but are also needed for the gauge anomaly cancellation. In such a framework, three of these SM singlet fermions, S i 1 , are completely decoupled due to the introduction of Z 2 symmetry and have naturally small mass (of order keV) according to 't Hooft's naturalness criterion [53]. Therefore, the lightest one, S 1 1 , will be a stable particle and hence a warm DM (WDM) candidate [54][55][56][57], as discussed in ref. [58]. Moreover, since these keV mass singlet fermions are odd under the Z 2 symmetry, they have no mixing with the active neutrinos and consequently are safe from the bound imposed by the xray observations [59]. In ref. [58], an extra moduli field was introduced to produce this keV WDM non-thermally to achieve the correct ballpark value of relic density consistent with the WMAP and Planck observations. In the current work, without introducing any extra field contrary to ref. [58], we successfully produce the keV WDM by the freeze-in mechanism through the decay and annihilation channels of Z . After explaining the keV FIMP WDM as a successful single component FIMP DM scenario to satisfy the correct value of the DM relic density, we study a two component FIMP DM as another possible scenario in the present model, where in addition to the FIMP WDM S 1 1 , the lightest heavy right-handed neutrino ν 1 H can be a FIMP DM (with mass of order GeV) by tuning its corresponding Yukawa coupling to be very small [60,61]. The GeV scale FIMP DM can be produced through the decay and annihilation processes of both the extra neutral gauge boson Z as well as the extra B − L Higgs h , while the keV FIMP WDM is produced only through the decay and annihilation processes of Z . 1 Here, "naturally" means the Dirac neutrino masses, MD, have the same size as the Dirac masses of the SM fermions and, in contrary to the usual type-I seesaw mechanism [42][43][44][45][46][47][48], large Dirac neutrino Yukawa couplings, λ d ∼ O(0.1), with right-handed neutrino (N i R ) masses are of order TeV.

JHEP06(2019)095
Lepton Fields The rest of the paper is organized as follows. In section 2 we discuss the B − L model with inverse seesaw mechanism and how the light neutrinos acquire their tiny masses. In section 3 we show that a keV sterile neutrino can be a WDM and produce the observed DM relic density as a single component FIMP DM. Section 4 is dedicated for studying two component FIMP type DM. Finally, our conclusions are given in section 5. In addition, it has three SM singlet fermions N i R (three right-handed neutrinos) with B − L charge = −1 that arise as a result of anomaly cancellation conditions. Included also is an extra SM singlet scalar φ H with B − L charge = −1, while φ h is the usual electroweak (EW) Higgs doublet. In order to satisfy the experimental measurements for the nonvanishing light neutrino masses with TeV scale right-handed (RH) neutrino using type-I seesaw mechanism, a very small Dirac neutrino Yukawa couplings, λ d O(10 −6 ) must be assumed [42][43][44][45][46][47][48]. Therefore, the mixing angle between the left-and right-handed neutrinos is quite suppressed, as it is proportional to λ d O(10 −6 ). As a consequence of such small mixing angle, the interactions between the RH neutrinos and the SM particles are very suppressed, making it difficult to observe them at the LHC [37][38][39][40][41]. Thus, we generate neutrino masses using the so-called inverse seesaw mechanism [49][50][51][52] that can naturally accommodate light neutrino masses with TeV scale RH neutrinos and large Yukawa couplings. In addition to the particle content as mentioned above, the BLSM with Inverse Seesaw (BLSMIS) has three extra pairs of SM singlet fermions (S i 1,2 , i = 1, 2, 3) with B − L charge = ∓2, respectively. In table 1, we show the complete particle spectrum for the BLSMIS model with their associated charges for different gauge groups. An additional discrete symmetry has been introduced, viz., Z 2 . All BLSMIS particles are even under this symmetry except S 1 which is odd. Due to this symmetry, terms like N c R φ † H S 1 and S 1 S 2 , that could spoil the usual inverse seesaw mechanism, are forbidden [49][50][51][52]. The complete Lagrangian for this model is given by

JHEP06(2019)095
where F µν = ∂ µ B ν − ∂ ν B µ is the U(1) B−L field strength, D µ is the covariant derivative, φ h = iσ 2 φ h and the flavor indices are omitted for simplicity. The general structure of the covariant derivative D µ in the present model takes the following form where (g c , T α , G α µ ) are the SU(3) C gauge coupling, generator and the gauge field, respectively. Similarly, (g, τ a , W a µ ), (g Y , Y, B µ ) and (g , Y BL , B µ ) are the corresponding quantities for SU(2) L , U(1) Y and U(1) B−L , respectively. It is worth mentioning that a kinetic mixing term F µν F µν is allowed and it leads to a non-vanishing Z-Z mixing angle, θ [62][63][64]. However, due to the stringent constraint from LEP experiments on the Z-Z mixing angle (|θ | 10 −3 ) [65][66][67], one may neglect this term. Finally, the potential V(φ h , φ H ) is given by [37,38] where the potential V(φ h , φ H ) will be bounded from below when the following inequalities are satisfied simultaneously Here, both the scalars φ H and φ h acquire their non-zero vacuum expectation values (VEVs), therefore, the B − L and the EW symmetries are broken spontaneously and the SM Higgs doublet φ h and the B − L singlet φ H take the following form: where v 246 GeV is the EW symmetry breaking scale and v is the scale of B−L symmetry breaking which is, in general, unknown and ranging from TeV to much higher scales. After breaking the B − L and the EW symmetries spontaneously, the extra neutral gauge boson Z acquires its mass (M Z = g v ) [37,38], 2 and the neutrino Yukawa interaction terms in eq. (2.1) and in addition a very small Majorana mass µ S for S 1,2 lead to the following neutrino mass terms 3 Therefore, the neutrino mass matrix in the basis (ν c L , N R , S 2 , S 1 ) can be written as The experimental search for Z , by LEP II [68,69], leads to another constraint: M Z /g 7 TeV. This constraint will easily be satisfied due to a smallness of g which is required by the freeze-in scenario [17]. 3 µS is naturally small due to 't Hooft's naturalness criterion [53], for simplicity we assume S1 and S2 have the same small Majorana mass (µS), and the generation of such small µS from non-renormalizable terms has been discussed in [70] and radiatively in [71].

JHEP06(2019)095
It is clearly seen that S 1 is completely decoupled and has no mixing with active neutrinos. It only interacts with the neutral gauge boson Z with a coupling g Z S 1 S 1 = g . Therefore, S 1 is free from cosmological and astrophysical constraints coming from active-sterile mixing [59]. Thus its mass is given as, After diagonalising the upper left 3 × 3 submatrix of the neutrino mass matrix M ν , the light and heavy neutrino masses, respectively, are given by where µ S M D , M N is assumed. One can naturally obtain eV scale light neutrino masses with µ S of order keV and M N of order TeV, keeping Yukawa coupling λ d of order one. Such large couplings between the heavy RH neutrinos and the SM particles leads to interesting implications and enhances the accessibility of TeV scale B − L at the LHC [72][73][74].
Recall that due to the added Z 2 symmetry, S 1 is completely decoupled. Hence the lightest fermionic singlet, S 1 1 , is a stable particle and hence a DM candidate. Since its mass (= µ S ) is of order keV, hence S 1 1 is a warm DM (WDM) candidate [58]. 4 Moreover, one can easily make the lightest heavy RH neutrino, ν 1 H , stable or long-lived by taking the corresponding Yukawa coupling to be very small 3 × 10 −26 (GeV/M N ) 1/2 [60,61]. Thus, from here onwards we focus on the two component DM scenario, where, one of them is GeV scale DM, ν 1 H , and the other is keV scale WDM, S 1 1 . It is important to note that due to the mixing term in the potential V(φ h , φ H ), the squared mass matrix of the neutral Higgs bosons in the basis (h, H) is non-diagonal and takes the following form: (2.11) Rotating this matrix into the basis (h 1 , h 2 ) which is defined as follows where the mixing angle α takes the following form: (2.13)

JHEP06(2019)095
Therefore, the masses of these two physical Higgs scalars (h 1 , h 2 ) are given by 5 (2.14) The quartic couplings λ's can be written in terms of the physical masses M h 1,2 as follows [78] We have used SARAH [79][80][81] to implement the BLSMIS and the relevant masses, couplings and decay widths have been calculated using SPheno [82].

Warm DM as FIMP
As mentioned earlier, S 1 1 is a WDM candidate with mass in the few keV range [83][84][85]. We next study in detail the production of this keV DM via the freeze-in mechanism. Here S 1 1 is produced solely from its coupling with the extra U(1) B−L gauge boson Z , as mentioned in the previous section. The corresponding gauge coupling g is taken to be extremely feeble ∼ O(10 −10 ) with the result that S 1 1 is never in thermal equilibrium with the cosmic soup. Due to this small B − L gauge coupling, the corresponding gauge boson Z also interacts very feebly with the cosmic soup and never attains thermal equilibrium [86], where Γ Z is the total decay width of Z and H is the Hubble parameter. Therefore, we first determine the distribution function for Z . 6 The general formalism to determine the distribution function of any particle (say f ) is to solve the following Boltzmann equation: whereL is the Lioville's operator and C[f ] is known as the collision term of f . If we consider an isotropic and homogeneous universe, then, using the Friedman-Robertson-Walker metric, the Lioville's operator takes the following form: Hereafter, the physical state h1 refers to the SM-like Higgs boson and its mass M h 1 is fixed at 125.5 GeV to agree with the LHC measurements [76,77]. Also, according to the measured values of Higgs boson signal strengths for its various decay modes, the mixing angle α should be very small, thus we have fixed it at 0.01 rad. 6 As Z is not in thermal equilibrium (due to very small value of g ), one can not assume a Maxwell-Boltzmann distribution function for Z . Therefore, the Z distribution can be found by solving eq. (3.2).

JHEP06(2019)095
where p is the absolute value of the particle's three momentum, | p|. Following [18], we perform a transformation of variables, (p, t) → (ξ p , z), in the following way: where g s (T ) is the effective entropy degrees of freedom (d.o.f) at temperature T , M sc is an arbitrary mass scale and hereafter we take it equal to the SM-like Higgs mass (M sc = M h 1 = 125.5 GeV) and T 0 is the initial temperature at which the DM relic density is taken to be zero. Therefore, using the following time-temperature relation, the Lioville's operator defined in eq. (3.3) can be simply written aŝ where g s (T ) is the derivative of g s (T ) with respect to the temperature T . Taking only the decay term for the Z production, 7 the Boltzmann equation of the distribution function of Z is given bŷ where f Z is the distribution function of Z , C h i →Z Z is the collision term of Z production from the decays of scalars h 1,2 and C Z → all is Z decay collision term due to all its possible decay channels. The expression of these collision terms are given in the appendix A. Once we get the distribution function of Z by solving eq. (3.7), we then can determine its co-moving number density by using the following relation: where n Z is the Z number density, g is the internal d.o.f of Z and the universe entropy density s is given by s = (2π 2 /45)T 3 g s (T ) [87]. From eq. (3.8), one can note that the co-moving number density of Z is directly proportional to the integrated ξ 2 p f Z , i.e., larger the area under a ξ 2 p f Z curve, larger is the Z abundance. In figure 1  then subsequently decays into WDM and the SM fermions. But as z is increased further (presently z = 100.0), Z starts decaying significantly and its abundance gets depleted and the area under the curve for z = 100 is smaller than for z = 20.0, as seen in figure 1. For still higher values of z (z = 500.0), Z abundance decreases further due to decay. Thus, as z → ∞, Z will gradually decay to DM and its abundance eventually goes to zero.
Once the distribution function of Z is computed, we can describe the production of the keV DM S 1 1 . In the present scenario, the keV DM S 1 1 can be produced from the decay of Z , Z → S 1 1 S 1 1 (decay contribution), and from the annihilation processes, ff → S 1 1 S 1 1 mediated by Z , where f = l, q, ν l . The annihilation contribution has been calculated by using micrOMEGAs [88]. To determine the co-moving number density of the WDM S 1 1 , we solve the following Boltzmann equation, where M Pl = 1.22 × 10 19 GeV is the Planck mass, , where g ρ is the effective energy degrees of freedom [5,89], and the non-thermal average of Z decay width is defined by  where The expressions of a thermal average annihilation cross section σv ff →S 1 1 S 1 1 and an equilibrium co-moving number density of f (Y eq f ), appearing in eq. (3.9), are given respectively by [5,89] σv ff →S 1 where g f is the internal d.o.f of f , K 1 (z) and K 2 (z) are the Bessel function for first and second kind, respectively, and σ ff →S 1 1 S 1 1 is given in ref. [90]. Solving the Boltzmann equation given by eq. (3.9) gives us the co-moving number density Y S 1 1 . The corresponding relic density of the WDM S 1 1 can be calculated by using the following formula [89], (3.14) In order to understand the relative contribution of the decay and annihilation channels we will first consider them one at a time and solve the Boltzmann equation to get the relic density. We start with taking only the decay contribution and show in the left and right panels of figure 2   measurement (Ωh 2 0.12) [4]. In the left panel, as long as T 0 is greater than the mass of the mother particles (h 1,2 ) in h 1,2 → Z Z decay channels, the final DM relic density is insensitive to T 0 , as seen for the z 0 = 0.001, 0.005, 0.01 cases. However, once T 0 drops below the mass of the mother particles (presently z 0 = 2 case), Z production gets Boltzmann suppressed and consequently DM relic density is suppressed. In the right panel we show the dependence of the DM relic density on its mass (M S 1 1 ) for z 0 = 0.01. It is clear that the relic density increases with the DM mass, as expected from eq. (3.14).
For the annihilation contribution, ff → S 1 1 S 1 1 , there are two relevant regimes are as follows. (i) The on-shell regime, where 2M S 1 1 < M Z < T 0 , in which Y S 1 1 (∝ g 4 /Γ Z ) does not depend on T 0 and (ii) The EFT regime, where 2M S 1 In the left panel of figure 3 we see that as long as T 0 is greater than M Z , the final DM relic density is insensitive to T 0 (on-shell regime), as seen for the z 0 = 0.001, 0.005, 0.01, 2.0 cases. Once T 0 drops below M Z (presently z 0 = 100.0 case), then S 1 1 production gets the suppressed by a factor T 3 0 /M 4 Z (EFT regime). In the right panel of figure 3, we show the dependence of the DM relic density on its mass (M S 1 1 ) for z 0 = 0.01 (on-shell regime). It is clear that the relic density increases with the DM mass, as expected from eq. (3.14).
For the decay contribution (Z → S 1 1 S 1 1 ), we show in figure 4 the variation of the comoving number density of Z (left panel) and the co-moving number density of S 1 1 with z, for different values of B − L gauge coupling g (left panel) and M Z (right panel). Since Z production is proportional to the B − L gauge coupling g , larger g results in larger Z production and consequently a larger production of DM, as seen in the left panel of this figure. Note also that Z → S 1 1 S 1 1 decay rate is directly proportional to g , and hence increasing g increases the decay rate of Z and hence the abundance of S 1 1 . Therefore, it is clear that for higher values of g the Z co-moving number density plateau starts bending at smaller values of z. On the other hand, in the right panel of figure 4, we see that by increasing M Z exactly opposite behavior appears for Z production while similar behavior    dashed line) and annihilation (green dotted line). Here, for a suitably selected set of model parameters (M S 1 1 = 10 keV, M Z = 10 GeV, g = 1.0 × 10 −9 , M h 2 = 1 TeV, z 0 = 0.01, and α = 0.01 rad), the total WDM relic density equals the observed relic density (Ωh 2 0.12) at the present epoch, where decay contributes ∼ 62% of the WDM relic density while the rest comes from annihilation. It is worth mentioning that initially for z 100, WDM is dominantly produced from the annihilation processes this is because of all ingoing particles are already in the cosmic soup, while for z 100, the decay process starts dominating, as seen in figure 4.
Variation of total WDM relic density (Ωh 2 ) as a function of the gauge coupling g can be seen in figure 7, where the BLSMIS points have been generated over the following ranges of its fundamental parameters: 1 ≤ M S 1 1 ≤ 10 keV, 1 ≤ M Z ≤ 100 GeV, 10 −12 ≤ g ≤ 10 −8 , M h 2 = 1 TeV, z 0 = 0.01, and α = 0.01 rad. From the left panel it is clear that Ωh 2 is inversely proportional to M Z (which is represented by the color bar). More explicitly, for a fixed g value, larger Ωh 2 values correspond to smaller M Z values (red points) and vice versa for the blue points. On the other hand as illustrated in the right panel, Ωh 2 is directly proportional to the WDM mass M S 1 1 (which is represented by the color bar). This is consistent with Ωh 2 expression given in eq. (3.14).
In figure 8 we show the allowed points in the (M Z , g ) and (M S 1 1 , g ) planes in the left and right panels, respectively, which give the relic density consistent with a relic density upper bound of the Planck measurement (Ωh 2 ≤ 0.12) [4]. All other parameter values are allowed to vary in the range mentioned in the previous paragraph. From the figure color bars (mapped to the total WDM relic density Ωh 2 ), it is clearly seen that many points (∼ 84% of the scanned points) have a small DM relic density (Ωh 2 ≤ 10 −2 ). Therefore, in the next section we discuss a two component FIMP DM possibility as a well-motivated scenario to get an extra relic density contribution from the lightest heavy RH neutrino, ν 1 H , as a GeV scale DM.

Two component FIMP dark matter
In the previous section we have studied the WDM FIMP S 1 1 , as a single component DM. As mentioned in section 2, the lightest heavy RH neutrino, ν 1 H , can be long-lived particle by making the corresponding Yukawa couplings very small 3 × 10 −26 (GeV/M N ) 1/2 [60,61]. Therefore, it can be an additional DM component of mass of order GeV. Note that any interaction between S 1 1 and ν 1 H is completely forbidden. Thus in the present section, we consider a two component DM scenario with two DM candidates: the WDM FIMP S 1 1 and the lightest heavy RH neutrino ν 1 H . The dominant annihilation channels of ν 1 H pair to SM particles are mediated by Z and h i (i = 1, 2). 8 The coupling strength of ν 1 H pair with Z 8 Due to the smallness of the corresponding Yukawa coupling of ν 1 H (as assumed to be a stable DM candidate), the contribution of the channels mediated by Z and W ± bosons is negligible. Also, the annihilation channels mediated by the SM-like Higgs h1 are suppressed as compared to the h2 ones, because the coupling λ ν 1 H ν 1 H h 1 is very small since it is proportional to sin α which is constrained to be very small by LHC [91].

JHEP06(2019)095
is given by g /2, while with h i (i = 1, 2) is given as where O 1 = sin α and O 2 = cos α. Therefore, ν 1 H pair annihilation is proportional to the gauge coupling g which is taken very small in the present model. Due to this feeble gauge coupling g , ν 1 H will never reach thermal equilibrium and is produced by the freeze-in mechanism. The Boltzmann equation associated with ν 1 H production is as follows where Γ h i is the total decay width of h i . After solving the Boltzmann equation of ν 1 H production, eq. (4.2), the corresponding relic density of ν 1 H can be determined by using the following relation, Finally, the total relic density of this two component DM scenario is given by where Ω S 1 1 h 2 is the relic density of S 1 1 which is defined in eq. (3.14). It is clear that the DM production depends crucially on the DM mass and the mass of the mother particles (   In discussed above, for M Z > 2M ν 1 H (region I), the total relic density is dominated by Z mediated diagrams. Now we turn to the region II, where M Z < 2M ν 1 H and Z decays to ν 1 H pair is kinematically forbidden, and consequently, ν 1 H production is h 2 dominated. Therefore in region II, a major portion of our two DM candidates, ν 1 H and S 1 1 , is produced almost independently from the h 2 and Z mediated processes, respectively. In other words, by fixing M h 2 , g and M ν 1 H /M Z at certain values to get a significant contribution from Ω ν 1 H h 2 , one can obtain a relevant Ω S 1 1 h 2 contribution independently by changing M S 1 1 within keV range. This possibility did not exist in region I because both ν 1 H and S 1 1 are produced dominantly via Z and therefore have the same number density. The only way to have comparable contribution from both in region I would be to raise the mass of S 1 1 to the

JHEP06(2019)095
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.