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 $\mathbb{Z}_{2}$ symmetry, a SM gauge singlet fermion, with mass of order keV, is stable and can be a warm DM candidate. Also, the same $\mathbb{Z}_{2}$ symmetry helps the lightest right-handed neutrino, with mass of order GeV, to be a stable or long-lived particle by making a corresponding Yukawa coupling very small. This provides a possibility of a two component DM scenario as well. Firstly, in the absence of a GeV DM component (i.e., without tuning its corresponding Yukawa coupling to be very small), we consider only a keV DM as a single component DM, which is produced by the freeze-in mechanism via the decay of the extra $Z'$ gauge boson associated to $U(1)_{B-L}$ and can consistently explain the DM relic density measurements. In contrast with most of the existing literature, we have found a reasonable DM production from the annihilation processes. After numerically studying the DM production, we show the dependence of the DM relic density as a function of its relevant free parameters. We use these results to obtain the parameter space regions that are compatible with the DM relic density bound. Secondly, we study a two component DM scenario and emphasize that the current DM relic density bound can be satisfied for a wide range of parameter space.


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, 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 x-ray observations [59]. In Ref. [58], an extra moduli field was introduced to produce this keV WDM non-thermally to achieve the correct ballpark 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.
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 ′ .
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.

B − L model with inverse seesaw scenario and neutrino masses
The gauged B − L extension of the SM (BLSM) is based on the gauge group SU (3) C × SU (2) L × U (1) Y × U (1) B−L . By imposing U (1) B−L , the gauge sector of the SM is extended to include a new neutral gauge boson Z ′ associated with the B − L gauge symmetry. 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 non-vanishing 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., 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 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-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 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) 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]. 4 The contribution of the new light degrees of freedom (S i 1 ) to the number of effective neutrino species, N eff , has been checked using Eq. (5) in Ref. [75] to calculate extra effective neutrino species, ∆N eff , and found it to be negligible. 5 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.
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: 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 Fig. 1, we show the variation of ξ 2 p f Z ′ with respect to the dimensionless parameter ξ p for different values of z (= M sc /T ). As shown in the figure, areas under the curves corresponding to z = 0.02 and 20.0 are different because for higher z = 20.0 (i.e., lower temperature T of the universe), Z ′ gets more time to be produced and it then subsequently decays into WDM and the SM fermions. But as z is increased further (presently z = 100.0), Z ′ starts decaying significantly and its abundance gets depleted and the area under the curve for z = 100 is smaller than for z = 20.0, as seen in Fig. 1. For still higher values of z (z = 500.0), Z ′ 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 7 In principle, the collision term for annihilation diagrams should also be considered but in this class of models those annihilation diagrams have subleading contribution [19], hence we have not taken into account those effects and for simplicity we consider only the decay of h1,2 as the Z ′ production mechanism. Moreover, h1,2 are in thermal equilibrium, and consequently the usual equilibrium Boltzmann distribution function has been assumed for them [18]. 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 , (3.10) 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 Fig. 2 the variation of DM relic density as a function of z, for different values of the initial temperature T 0 (= M sc /z 0 ) and different values of the WDM mass M S 1 1 , respectively. The horizontal magenta dashed line refers to the DM relic density measurement (Ω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).
In the left panel of Fig. 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 Fig. 3, we show the dependence of the DM relic density on its mass (M S 1 1 ) for z 0 = 0.01 (on-shell regime). It is clear that the relic density increases with the DM mass, as expected from Eq. (3.14).
For the decay contribution (Z ′ → S 1 1 S 1 1 ), we show in Fig. 4 the variation of the co-moving 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 Fig. 4, we see that by increasing M Z ′ exactly opposite behavior appears for Z ′ production while similar behavior happens for its decay. As mentioned, Z ′ production mainly happens through decay of the Higgs bosons (h 1,2 ) and those decay modes . Therefore, increasing M Z ′ reduces the production rate of Z ′ as 1/M 2 Z ′ . However, its decay width is simply proportional to its   In Fig. 6, we show the total relic density (blue dashed-dotted line) as well as the relative contributions of the two different types of WDM production processes, decay (red dashed line) and annihilation (green dotted line). Here, for a suitably selected set of model parameters (M S 1 1 = 10 keV, M Z ′ = 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 Fig. 4.
Variation of total WDM relic density (Ωh 2 ) as a function of the gauge coupling g ′ can be seen in Fig. 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 Fig. 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   h 2 ), respectively, while blue dashed line corresponds to ν 1 H relic density contribution from decay of h 2 (Ω h 2 ν 1 H h 2 ). From the figure, it is clearly seen that S 1 1 has a relevant relic density contribution, unlike the situation in region I. Note that for a larger S 1 1 mass (M S 1 1 = 100 keV), the S 1 1 contribution to the total relic density even starts to be the dominant one, as seen in the right panel of Fig. 11.

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