Freeze-in production of sterile neutrino dark matter in U(1)B−L model

With the advent of new and more sensitive direct detection experiments, scope for a thermal WIMP explanation of dark matter (DM) has become extremely constricted. The non-observation of thermal WIMP in these experiments has put a strong upper bound on WIMP-nucleon scattering cross section and within a few years it is likely to overlap with the coherent neutrino-nucleon cross section. Hence in all probability, DM may have some non-thermal origin. In this work we explore in detail this possibility of a non-thermal sterile neutrino DM within the framework of U(1)B−L model. The U(1)B−L model on the other hand is a well-motivated and minimal way of extending the standard model so that it can explain the neutrino masses via Type-I see-saw mechanism. We have shown, besides explaining the neutrino mass, it can also accommodate a non-thermal sterile neutrino DM with correct relic density. In contrast with the existing literature, we have found that W± decay can also be a dominant production mode of the sterile neutrino DM . To obtain the comoving number density of dark matter, we have solved here a coupled set of Boltzmann equations considering all possible decay as well as annihilation production modes of the sterile neutrino dark matter. The framework developed here though has been done for a U(1)B−L model, can be applied quite generally for any models with an extra neutral gauge boson and a fermionic non-thermal dark matter.


Introduction
The existence of Dark Matter (DM) in the Universe is now an acceptable reality.There are various satellite borne experiments, namely WMAP [1] and Planck [2] who have already measured the current mass density (relic density) of DM in the Universe with an extremely good accuracy.Moreover, there are also some indirect evidences about the existence of dark matter such as flatness of galactic rotation curve [3], gravitational lensing of distant object [4], bullet cluster [5] etc. However the composition of DM is still unknown to us.The Standard Model (SM) of electroweak interaction does not have any fundamental particle which can play the role of DM.Hence in order to accommodate a viable dark matter candidate we need to formulate a theory beyond Standard Model (BSM) of electroweak interaction.Among the various possible BSM theories available in literature the Weakly Interacting Massive Particle (WIMP) is the most favourable class of dark matter candidates and until now neutralino in the Supersymmetric Standard Model is one of the most studied WIMPs [6].The presence of DM is also being investigated in various direct detection experiments, namely LUX [7], XENON 100 [8] etc, and no "real signal" due to a dark matter particle has been observed yet.With the increasing sensitivity of the direct detection experiments ("ton-scale") [9,10,11], the WIMP-nucleon cross section is soon to merge with the elastic neutrino-nucleon cross section [12,13].The floor mostly comprises of 8 B and 7 Be solar neutrinos [14].Hence in future our only probe to distinguish a dark matter signal (assuming that DM is a thermal WIMP) from the neutrino background will be through directional searches [15].But if we wish to move beyond this thermal WIMP scenario, there is another class of dark matter candidates which are produced through nonthermal processes at an early stage of the Universe.Such possibilities include axino [16,17,18], gravitino [19,20], very heavy dark matter candidates like WIMPzillas [21] among many others [22].Their interaction strengths with other particles (in the thermal plasma) are so feeble that they never attain thermal equilibrium.These types of dark matter candidates are known as Feebly Interacting Massive Particle or FIMP [23].In contrast with the commonly discussed WIMP scenario, the relic density of FIMP type dark matter is attained by the so called Freezein mechanism [23].Unlike the thermal Freeze-out mechanism where relic density depends on the final abundance of dark matter, in Freeze-in, DM relic density is sensitive to its initial production history (for a nice review see [24]).In literature two types of Freeze-in mechanisms are usually discussed, IR (infra-red) Freeze-in [25,26,27] and UV (ultra-violet) Freeze-in [28,29,30].Unlike the former, the DM relic density in UV Freeze-in depends explicitly on the reheat temperature (T R ).Production of the non-thermal DM candidate usually occurs via a decay of a heavy mother particle (e.g. from Inflaton decay and decay of heavy Moduli fields [31,32]).
In this work we will study a FIMP type dark matter candidate in the U(1) B−L extension of the Standard Model of particle physics.U(1) B−L extension of SM is a very well motivated BSM theory as it provides the explanation of nonzero neutrino mass through Type-I sea-saw mechanism.In this model besides the usual SM gauge (SU(3) c × SU(2) L × U(1) Y ) symmetry, an additional local U(1) B−L symmetry invariance is also imposed on the Lagrangian where B and L respectively represent the baryon and lepton number of a particle.In order to obtain an "anomaly free gauge theory", three additional right handed neutrinos (N i , i = 1, 3) are required to be added to the particle spectrum of SM.Moreover, we also require a complex scalar (Ψ) which is a singlet under the SM gauge group but possesses a suitable nonzero U(1) B−L charge.Majorana masses for the three right handed neutrinos are generated through the spontaneous breaking of the local B − L symmetry by the vacuum expectation value (VEV) of complex scalar singlet Ψ.The lightest one (N 1 ) among the three right handed neutrinos can be a viable dark matter candidate.The dark matter candidate N 1 in U(1) B−L model can be produced through both thermal as well as non-thermal processes.In the former case, the interaction strengths of DM particles with others in the early Universe are such that they are able to maintain their thermal as well as chemical equilibrium.The decoupling of the DM particles occur when their interaction rates fall short of the expansion rate of the Universe.If n eq and σv are the equilibrium number density and the thermally averaged annihilation cross section of N 1 then the decoupling condition requires neq σv H < 1 with H being the Hubble parameter.Being out of equilibrium, the relic density of N 1 freezes to a particular value which depends upon the interaction strength as well as the temperature of the Universe at which the decoupling occurred (freeze-out temperature).The thermally produced N 1 as a dark matter candidate, in the U(1) B−L extension of SM, has been studied in Refs.[33,34,35,36].In these works most of the authors have shown that the relic abundance of dark matter particle satisfied the WMAP or Planck limit only when the mass of DM is nearly half the masses of mediating scalar particles (at or near resonances).This requires significant fine tuning as there is no symmetry, in the Lagrangian, which can relate the masses of dark matter and the scalar sector particles in the above mentioned way.Hence, with respect to the above discussions, it is natural to think about a dark matter particle, in this U(1) B−L model, which is produced through some non-thermal interactions at the early stage of the Universe.Non-thermal sterile neutrino production from the oscillation of active neutrinos was first proposed by Dodelson-Widrow [37], but this idea is now in conflict with the X-ray observations [38].Other mechanisms of sterile neutrino production like Shi-Fuller mechanism [39] can alleviate some of these problems producing a colder dark matter spectrum.Several other models have also successfully discussed non-thermal sterile neutrino dark matter.They include some Supersymmetric models [40], models using warped extra-dimensions [41] and decay from charged [42] and neutral scalars [43,44] or from extra gauge bosons [45,46].Most of the studies involving production of sterile neutrino from extra gauge boson assume the gauge boson to be in thermal equilibrium with the other SM particles.However, in this work we have moved way from this assumption (details later).Several non-thermal models of sterile neutrino dark matter under the assumption of low reheating temperature have also been studied in [47,48,49] which is also not the case we are considering here.
Additionally, unlike what is usually done in building a dark matter model, we do not impose any extra symmetry to stabilise our dark matter candidate.For an O (MeV) sterile neutrino dark matter we have a dominant decay mode to e ± and ν with a very large life time (larger than the age of the Universe for the parameters we consider here) which in turn helps us to propose a possible indirect detection signal of the 511 keV line observed by INTEGRAL/SPI [50] of ESA.
The rest of the paper is organised as follows: In Section 2 we briefly describe the U(1) B−L model.In Section 3 we describe the production mechanism of non-thermal sterile neutrino dark matter in detail.Section 4 describes the Boltzmann equation(s) needed to compute the comoving number densities of both Z BL and N 1 .Calculation of relic density of sterile neutrino dark matter is given in Section 5. Section 6 deals with a possible indirect detection mode of our dark matter particle N 1 .Finally our conclusion is given in Section 7. All analytic expressions of decay widths and annihilation cross sections used in this work are listed in the Appendix.

The U(1) B−L extension of Standard Model
In the present work we have considered a minimal U(1) B−L extension of the Standard Model where the SM gauge sector is enhanced by an additional local U(1) B−L gauge symmetry with B and L are known as the baryon and lepton number of a particle.Therefore, under the U(1) B−L gauge group all SM leptons (including neutrinos) and quarks have charges −1 and 1 3 respectively.Besides the SM fields, this model requires the presence of three right handed neutrinos (N i , i = 1, to 3) with U(1) B−L charge −1 for anomaly cancellation.On the other hand, as the SM Higgs doublet (Φ) does not possess any B − L charge, hence in order to spontaneously break the local B − L symmetry one needs to introduce a scalar field which transforms nontrivially under the U(1) B−L symmetry group.As a result, the scalar sector of the present model is composed of a usual Higgs doublet (doublet under SU(2) L ) Φ and a complex scalar singlet Ψ.To generate Majorana mass terms in a gauge invariant manner for the three right handed neutrinos one needs the B − L charge of Ψ is +2.B − L symmetry is spontaneously broken when Ψ acquires VEV v BL while the remnant electroweak symmetry (SU(2) L × U(1) Y ) of the Lagrangian breaks spontaneously through the usual Higgs mechanism.In unitary gauge, the expressions of Φ and Ψ, after getting VEVs v and v BL respectively, are The gauge invariant and renormalisable Lagrangian of the scalar sector is thus given by with where are the covariant derivatives of the scalar doublet Φ and complex scalar singlet Ψ respectively while where θ is the mixing angle between the neutral scalars h and H.The expressions of mixing angle (θ) and masses (M h , M H ) of h and H are given by We have considered the physical scalar h as the SM-like Higgs boson which was discovered recently by ATLAS [51] and CMS [52] collaborations and consequently we have fixed the value of M h at 125.5 GeV.Also according to the measured values of Higgs boson signal strengths (for its various decay modes) the mixing angle θ between the SM-like Higgs boson h and extra scalar boson H should be very small.As this mixing angle does not play any significant role in the present context, we have kept θ fixed at 0.1 rad [53], throughout this work, such that it satisfies all results from both ATLAS and CMS collaborations.Besides this, in order to obtained a stable vacuum the quatic couplings of the Lagrangian (Eq. 3) must satisfy the following inequalities, Moreover, as both Φ and Ψ have nonzero VEVs, this requires µ 2 i < 0 (i = 1, 2).The gauge sector Lagrangian of the present model is given as Here, L SM gauge is the Lagrangian of the SM gauge sector while the second term represents the kinetic term for the B − L gauge bosons Z BL and in terms of Z BL the field strength tensor Z ′ µν for an abelian gauge field is defined as The gauge invariant Lagrangian for the three right handed neutrinos can be written as: where Φ = −iτ 2 Φ * and D / N = γ µ D µ N with is the covariant derivative for the right handed neutrino N i .After U(1) B−L symmetry breaking the masses of right handed neutrinos and Z BL are given by 1 In general we may also have a kinetic mixing term given by κZ µν Z ′ µν .The value of κ is however severely constrained by electroweak precision measurements (κ < ∼ 10 −4 [54]).So for calculational simplicity we have restricted ourselves to a parameter space where κ < g BL , hence neglecting its contribution.These type of scenarios where kinetic mixing term is neglected has been previously studied under the name of "Minimal/Pure"U(1) B−L model [33,34,55].
Using above two equations one can write the coupling λ R i in terms of g BL , M N i and M Z BL which is From Eq. ( 10) it is possible to generate neutrino masses via Type-I see-saw mechanism.In our analysis we want to focus on the viability of lightest sterile neutrino (N 1 ) as a dark matter candidate.So for simplicity we have neglected intergenerational mixing between the active and sterile neutrinos.The mass of the other two sterile neutrinos are also not constrained by our analysis in this work, and in principle can be very heavy aiding neutrino mass generation by the see-saw mechanism.From Eq. ( 10) and Eq. ( 13), one can find the expression of active-sterile mixing angle α i per generation as tan 2 For simplicity, throughout this work we have denoted the first generation active-sterile mixing angle α 1 by only α.
The non-observation of the extra neutral gauge boson in the LEP experiment [56,57] imposes following constraint2 on the ratio of M Z BL and g BL : In our analysis independent parameters are: Mass of the extra singlet Higgs M H , Masses of all three RH neutrinos M N i , mass of extra neutral gauge boson M Z BL , scalar mixing angle θ, the new gauge coupling g BL and active-sterile mixing angle α.In terms of our chosen independent set of model parameters, the other parameters appearing in Eq. ( 3) can be written as 3 Exploring the Non-thermal Regime Non-thermal production mechanism of dark matter has been studied for quite a long time.Their characteristic behaviour comes from the very low cross section with the Standard Model particles in the early Universe.Due to this very low cross section (lower than that of WIMPs), the nonthermal dark matter particles can never reach in thermal equilibrium with the Standard Model particles.Hence their evolution in the early Universe is studied differently than the thermal scenario.In the thermal scenario, the abundance of a relic particle (called WIMP) remains nonzero in the present epoch due to the "Freeze-out" mechanism [59], whereas in the case of a non-thermal production of DM (called FIMP), a different mechanism known as "Freeze-in" [23] is responsible for their relic abundance.In the non-thermal case, due to very low interaction cross section, the initial abundance of the dark matter is taken to be zero.As the Universe cools, they are dominantly produced by the decay of other SM/BSM particles.They can also be produced by the scattering of SM/BSM particles, but with a sub-dominant contribution.Once the non-thermal dark matter is produced, due to extremely low interaction strength, they do not thermalise with the rest of the thermal soup.Since most of the production of DM particles in the non-thermal regime occur from the decays of heavier particles, non-thermality condition will be satisfied when the rate of production from the decaying mother particle (decay width) is less than the expansion rate of the Universe at around a temperature T ∼ M, where M is the mass of the decaying particle [60].Mathematically this can be written as where, Γ is the relevant decay width and H is the Hubble parameter.However in some cases, if the production of DM particles may occur mainly from the annihilation of other particles in the thermal bath (production from decay can be forbidden due to kinematical condition or by some symmetry in the Lagrangian).Γ will then be replaced by: where, σv is the thermally averaged annihilation cross section of the particles in the thermal bath and n eq is their equilibrium number density.In this U(1) B−L model, to calculate the relic density of a non-thermal sterile neutrino dark matter (N 1 ), the principal ingredient is its production from various decay and annihilation channels.This gives the required comoving number density of N 1 upon solving the relevant Boltzmann equation.The main production channels of the sterile neutrino (via decay) are : The corresponding decay widths are given in the Appendix A.1.As discussed earlier, non-thermal dark matter particles can also be produced from the scattering of the SM/BSM particles in the thermal soup.The rate of the back reactions are negligible, since the number density of N 1 is extremely small in the early Universe.The annihilation channels along with corresponding cross sections aiding the production of N 1 are also given in the Appendix A.2.As we will see later, in the present case W ± and Z BL decays are main production channels of N 1 .Using the non-thermality condition given in Eq. ( 22) we find that the extra gauge coupling g BL and the active-sterile mixing angle α must be less than 10 −9 and 10 −7 (rad) respectively for an O(MeV) sterile neutrino with the mass of Z BL lying in 1 GeV to 100 GeV range.Although this is a simple way to estimate the order of magnitude of g BL and α required for the dark matter candidate (N 1 ) to be non-thermal, this sets a very first upper limit on these quantities.However, more stringent upper bound on α (α < ∼ 10 −9 rad) arises from the stability of DM over the cosmological time scale.
Moreover, an upper bound on the active-sterile mixing angle α is also obtained from the invisible decay of the Standard Model Z boson.Following Ref. [61] we find: In the limit when active-sterile mixing angle is small and M Z ≫ M N 1 , from the above equation we get sin 4 α < 0.007.As we will see later that for us, this condition is indeed being satisfied.
In the present scenario since M h < 2 M H , SM Higgs boson can decay invisibly only into a pair of lightest sterile neutrino N 1 .From the expression of the decay width given in Eq. ( 37) we find that it is suppressed by g 2 BL and hence very small.Thus this decay width easily satisfies the bound on invisible decay of SM Higgs boson from LHC [62].Furthermore, due to sufficiently small interaction strength with the SM particles, non-thermally produced N 1 always satisfies all the existing bounds on spin independent as well as spin dependent scattering cross sections from dark matter direct detection experiments [7].
We have mentioned earlier that for the non-thermal production of the sterile neutrinos, the coupling constant g BL should be very small (<10 −9 ).As is usually done, while considering the production of dark matter from a decay of any SM/BSM particle, the latter is implicitly assumed to be in thermal equilibrium.Hence we usually do not need to solve a system of coupled Boltzmann equations, since the equilibrium number density is assumed for the decaying mother particle.But, here due to very low interaction strength of Z BL (due to small g BL ), it will not be in thermal equilibrium with the rest of the particles.Also, the decay of Z BL is a mode of production of the our sterile neutrino dark matter N 1 .So, first, we find the comoving number density of Z BL by solving its Boltzmann equation.Then we use this to find the relic density of our sterile neutrino dark matter.Thus, in our case we have to solve a set of two coupled Boltzmann equations, one for the sterile neutrino dark matter, and another for the Z BL .
In any model with a sterile neutrino we will have an active-sterile mixing in general.Hence in such model production of the sterile neutrino via W ± decay is a very generic feature.But, it is usually not taken into account since it is suppressed by the square of the small active-sterile mixing angle.However, in this work, in our favoured parameter space, we find that a sizeable contribution (to the relic-density of N 1 ) even from the W ± decay is present (see Section 4.1).Another important feature which will be present for a generic model having an nonzero active sterile mixing is the production of sterile neutrino through the Dodelson-Widrow (DW) mechanism.Here the production of sterile neutrino occurs via the oscillations of active neutrinos to the sterile ones.But this mechanism suffers serious drawbacks from the Lyman-α bounds [63] as well as X-ray observations [38].It is now known [64,65] that sterile neutrino produced by this mechanism cannot comprise the whole of the dark matter of the Universe.The contribution arising to the relic abundance of a sterile neutrino from the DW mechanism is given by [66] where, α is the active sterile mixing angle and M N 1 is the mass of the sterile neutrino.In our case we find (see Section 5 for details) that in order to satisfy relic density, α should be less than 10 −10 rad for sterile neutrino mass lying between 1 MeV and 10 MeV.From Eq. ( 25) we see that the corresponding DW contribution to the relic density is < ∼ 1.2 × 10 −6 and hence negligible.

Boltzmann Equation
In this section, we write the two coupled Boltzmann equations that dictates the final relic abundance of the sterile neutrino dark matter N 1 .The Boltzmann equation for the evolution of Z BL which, as already discussed is very weakly interacting is given by3 : Here, Y Z BL ≡ n Z BL s is the comoving number density of the extra gauge boson with n Z BL and s being the number density of Z BL and the entropy density of the Universe respectively.Also where Λ is a mass scale and T is the temperature of the Universe.For simplicity we have taken Λ ∼ M h , the mass of SM Higgs boson while M pl is the usual Planck mass.The function g ⋆ (z) is given by: where, g ρ (z) and g s (z) are the effective degrees of freedom related to the energy density ρ and the entropy density s of the Universe respectively.The quantity Γ A→BB denotes the thermally averaged decay width for the process A → BB and its expression, in terms of decay width Γ A→BB , is given by [67] 4 : Here, K 1 (z) and K 2 (z) are the modified Bessel functions of order 1 and 2 respectively.The expressions for the relevant decay widths are given in Appendix A.3.
The SM particles acquire their masses after the process of EWSB whereas the BSM particles like U(1) B−L gauge boson Z BL and the additional Higgs boson (H) gain their masses after the breaking of U(1) B−L symmetry.Therefore, in the early Universe the main production channel of the new gauge boson is mainly through the decay of H, while the latter is in thermal equilibrium with the plasma.The first term in Eq. ( 26) denotes this contribution to the production of Z BL (i.e. increase in number density of Z BL ) and hence comes with a positive sign.Since in our case, both the masses of Z BL and H are free parameters, we have adopted the values of M H in a range such that it always satisfy the kinematical condition M H ≥ 2 M Z BL .Although in the early stage of the Universe, the decay of H is the main production channel of Z BL , in principle it can also be produced from the annihilation processes, involving both SM as well as BSM particles, like hh However, contribution of these annihilation processes is subleading to that of decay.The number density of extra gauge boson Z BL is also depleted mainly through its decay modes to N 1 N1 (other two sterile neutrinos are assumed to be heavy for simplicity), ν x νx and f f .It is denoted by the second term in the Boltzmann equation (Eq.( 26)), and as expected it comes with a negative sign, since it signifies the depletion of Z BL number density.These two competing processes (production vs. depletion) decide the final comoving number density of Z BL .Numerically solving Eq. ( 26), we graphically show the evolution of comoving number density of Z BL with z = M h T in Fig. 1.From the above plot it is seen that, the comoving number density of Z BL first rises due to the production term (first term in the R.H.S. of Eq. ( 26)) and then after a certain time it falls when the depletion term (i.e. the second term in the R.H.S. of Eq. ( 26)) begins to dominate.This situation arises because, at that time the temperature of the Universe becomes much smaller than M H (T ≪ M H ) and hence, being a non relativistic species, the equilibrium number density of H is exponentially suppressed.Therefore, the production of Z BL ceases.The middle "plateau" like portion occurs when both the production and the depletion terms are comparable and hence compensating each other.The plot is generated for the following chosen set of relevant parameters: M H = 500 GeV, M Z BL = 10 GeV and g BL = 10 −10 .Now, we proceed to write the Boltzmann equation of the lightest sterile neutrino N 1 .This will govern the number density of the dark matter candidate (N 1 ) and consequently its relic abundance at the present epoch.Similar to Eq. ( 26), the Boltzmann equation for N 1 is given by : As discussed in Eq. ( 26), since the initial abundance of the sterile neutrino dark matter N 1 is very small, the Y N 1 term in the above equation may be neglected [23,43].Here σv xx→N 1 N 1 is the thermally averaged cross section for the production of N 1 from the annihilation of x particle.The expression of σv xx→N 1 N 1 is given by [59] 5 The expressions for the relevant decay widths and annihilation cross sections are given in the Appendix A.1 and A.2 respectively.In order to get the comoving number density (Y N 1 ) of N 1 at the present epoch, we have to solve the coupled set of Boltzmann equations given in Eq. (26,28).To be more precise, the value of Y Z BL (z) for each z obtained by solving Eq. ( 26) is to be fed into Eq.(28).For understanding the physics behind this coupled set of Boltzmann equations better, let us first assume that the lightest sterile neutrinos are only produced from the decay of Z BL , whereas Z BL is produced and depleted according to Eq. (26).We plot the result in Fig. 2. In Fig. 2, initially at the early stage of the Universe there are no Z BL particles and hence no N 1 , since for simplicity we have switched off all other production channels of N 1 , except Z BL .Then, when Z BL is produced from the decay of H, we also find an increase in the number density of N 1 from the decay of Z BL .Finally, the number density of Z BL begins to fall due to its dominating decay modes (production of Z BL terminates as the the number density of H becomes negligibly small), and consequently the number density of N 1 also saturates since now there are no Z BL left to aid the production of N 1 .This plot is also generated for the following chosen set of relevant parameters: M H = 500 GeV, M Z BL = 10 GeV, M N 1 = 1 MeV and g BL = 10 −10 .In Fig. 2 we have taken the initial temperature (T i ) to be 1 TeV.The final abundances of Z BL and N 1 will not depend on this as long as T i > ∼ M H .The maximum production of Z BL from H decay occurs around a temperature of ∼ M H .However, if T i becomes less than M H then, since H is in thermal equilibrium, its own abundance will be exponentially suppressed (as it becomes non-relativistic) and thereby reducing Y Z BL and Y N 1 .All these are shown in Fig. 3 where we see as discussed above, for T i > ∼ M H , there is no change in the final values of Y Z BL and Y N 1 (red, green, blue, cyan solid lines).While for T i < ∼ M H both the final abundances are reduced (black solid line) from their previous values.
We now show the variation of Fig. 2 with different sets of chosen model parameters.In Fig. 4(a) we show the variation of Y Z BL and Y N 1 with z for two different values of U(1) B−L gauge coupling g BL .In this plot, we find that with increase in the value of g BL the number density of Z BL also increases initially.This is understandable, because the decay width Γ H→Z BL Z BL increases with g BL , and hence resulting in an increased number density of the extra gauge boson.But the depletion rate of Z BL (proportional to its total decay width) also increases with g BL , and hence will result in a faster fall of its comoving number density.This is also evident from the figure where the green line starts to fall earlier than the red one.On the other hand, the production rate of sterile neutrino dark matter (N 1 ) is proportional to Γ Z BL →N 1 N 1 and Y Z BL (see Eq. ( 28)) and both of these quantities increase with g BL .Therefore, the comoving number density of N 1 increases as the value of g BL changes from 1 × 10 −10 to 5 × 10 −10 .In Fig. 4(b), we show the variation of Y Z BL and Y N 1 with z for two different values of M H . Now, with a decrease in M H , we expect a corresponding decrease in decay width Γ H→Z BL Z BL (Z BL production rate), and hence the initial number density of Z BL will be smaller.This feature is seen in plot (b) where initially (z < ∼ 10 3 ) Y Z BL for M H = 500 GeV (green line) is larger than that for M H = 50 GeV (red line).However, the total decay width of Z BL (and consequently its depletion rate) does not depend on the mass of H. Hence both the red and green lines will start to fall off around the same time.Another noticeable change due to the variation of M H is that the width of the "plateau" becomes narrower with the decrease in mass difference between M H and M Z BL .Further, as the Y Z BL increases with an increase in M H , which in turn produces more N 1 (from the decay of Z BL ) in the final state and hence the comoving number density Y N 1 also increases with M H .In Fig. 4(c) we have shown the variation of Y Z BL and Y N 1 for three different values of sterile neutrino dark matter mass.From this plot we find that there is not much variation in Y Z BL with changing M N 1 .This is because the decay width Γ Z BL →N 1 N 1 is subdominant with respect to the other decays modes of Z BL .Increasing the mass of N 1 will lead to a further decrease of Γ Z BL →N 1 N 1 and hence will not affect the depletion rate which is dominantly controlled by the other decay channels of Z BL .But since in this case (when other production channels of N 1 are switched off) Γ Z BL →N 1 N 1 solely controls the production rate of N 1 , Y N 1 changes with M N 1 .It is seen from Fig. 4(c) that if we increase the sterile neutrino mass from 1 MeV to 4 GeV (M N 1 tends to M Z BL /2),

Solution of the complete Boltzmann equation(s) with all production and decay channels
In the previous section, we demonstrated the validity of the coupled set of Boltzmann equations needed to solve for the relic abundance of the sterile neutrino dark matter N 1 .For simplicity, we assumed that the only production channel of N 1 is the decay from Z BL .However in general, all the possible production modes of N 1 including decays as well as annihilations of SM and BSM particles, as given in Eq. ( 28), have to be taken into account.Therefore, the active-sterile mixing angle is now nonzero.The noticeable feature when the active sterile mixing is nonzero is the production of N 1 from the decay of W ± bosons (W ± → e ± N 1 ).It may a priori seem that due to small value of active-sterile mixing angle the contribution from the decay of W ± will be negligible, but we have to remember that in this non-thermal scenario, the extra gauge coupling (g BL ) is also required to be very small (∼ 10 −10 ), and hence the production of N 1 from the decay of Z BL may also compete with the former.We will show this quantitatively later.Also note the decay of W ± is solely governed by the active-sterile mixing angle α and does not depend on U(1) B−L gauge coupling g BL , hence if g BL is made very low, the only dominant production channel of N 1 will be from W ± decay.In is the existence of a "double plateau".The reason behind this is the presence of another the production mode of N 1 from W ± decay which was neglected in previous section for simplicity.The onset of N 1 production as seen in these plots is a little early than those seen in Fig. 4. The initial onset here is due to the presence of W ± decay and is independent of g BL , M Z BL and M H (Eq. ( 33)), and only depends on α which is the sole parameter that controls the W ± → e ± N 1 decay.The first plateau occurs when the number density of W boson begins to fall and hence there is a decreased rate of production of N 1 .However it again begins to rise sharply when the production from Z BL starts to dominate.Then as before, we can see by comparing with the accompanying Y Z BL lines that the second plateau results when the Z BL number density starts to deplete.It is to be noted that this "two plateau" feature will be visible only when the production from W and Z BL are comparable to each other at some point of z(= M h T ).If either one of them remains dominant for all z, then it will result in single plateau like feature.For example, the green solid line in plot (c) of Fig. 5 has only a single plateau.This is because, due to a very high value of α the decay channels of W ± remain the most dominant production mode of N 1 for all z.The corresponding variation of Z BL number density is also shown by red solid line and since α has no effect on the production and decay channels of Z BL we find no variation of Y Z BL with α.However the variation of Y N 1 is different for different values of α, as the active-sterile mixing angle α controls the production mode of N 1 from W ± decay.

Relic Density of Sterile Neutrino Dark Matter (N 1 )
In order to compute the relic abundance (Ω N 1 h 2 ) of the lightest sterile neutrino (N 1 ) we need to find the value of its comoving number density (Y N 1 ) at the present epoch (T = T 0 , T 0 ∼ 2.73 K).The value of Y N 1 (T 0 ) can be obtained by solving the two coupled Boltzmann equations (Eqs.(26,28), which we have discussed elaborately in Section 4. The expression of Ω N 1 h 2 in terms of Y N 1 (T 0 ) is given by [68], In this work, we take all decay and annihilation channels of both SM as well as BSM particles for the production of N 1 .In Fig. 6 we show the relative contributions to Ω N 1 h 2 from W ± (red solid line) and Z BL decay (green solid line) for some chosen sets of model parameters.The total relic abundance of N 1 is also shown by the blue solid line.For some combinations of model parameters we find W ± decay can be the leading production channel of N 1 (plot (a)) while for some others it can be the subleading one (plot (c)).However in all three plots (a-c) of Fig. 6, the total relic density of N 1 has the saturation value ∼ 0.12 which is in conformity with the value of dark matter relic density measured by the satellite borne experiment Planck.In plot (b) of Fig. 6 we show a situation when the relative contributions to Ω N 1 h 2 from both W ± and Z BL decays are equal.In this case we have adopted the following values of relevant model parameters: From Table 1 it is seen that for the small values of g BL ∼ 10 −11 and α ∼ 10 −10 (which are required for the non-thermality of N 1 ) the contributions of other production channels of N 1 through the decays of Z, H and h are negligible.Similarly, Table 2 shows that within this adopted ranges of model parameters the annihilation processes of SM as well as BSM particles do not contribute significantly to the production of sterile neutrino dark matter N 1 and hence we can safely consider the decays of W ± and Z BL as the two most efficient production mechanisms of N 1 .
In Fig. 7, we plot the allowed values of B − L gauge coupling g BL and active sterile mixing angle α which satisfy the relic density criteria (0.1172 ≤ Ω N 1 h 2 ≤ 0.1226) [2].During the computation of Fig. 7 we have varied the relevant parameters in the following range.
and we have kept the mass of the extra Higgs boson H fixed at 500 GeV.From Fig. 7, we see that for very small values of the extra gauge coupling g BL (∼ 10 −12 to 10 −15 ) the relic density condition of N 1 is always satisfied for a active-sterile mixing angle α ∼ 10 −10 rad.This is expected since for very small values of g BL the production of N 1 from the decay of Z BL is highly suppressed and in this situation decay of W ± becomes the principle production channel, since the latter is not suppressed by the extra gauge coupling (see Eqs. ( 33), ( 35)).Earlier works about the nonthermal production of sterile neutrino have not touched upon this point in detail (production of N 1 from W ± decay), since most of the previous authors have ignored the N 1 production mode from W ± decay in their works by assuming extremely small values of active-sterile mixing angle α.However, such an assumption needs careful attention when other couplings in the theory can also be very small.On the other hand for the higher values of g BL (∼ 10 −9 to 10 −11 ) the decay of Z BL becomes the dominant contributor and hence in this case small values of α are required to suppress the contribution of W ± decay to Ω N 1 h 2 such that the total relic density of N 1 lies within the range prescribed by the Planck experiment.Since the production of sterile neutrino from W ± decay depends only on the mixing angle α, in Fig. 7 we get only a narrow band of α which satisfies the relic density of N 1 (for small values of g BL ).But when Z BL is the main production channel of N 1 then for a given g BL and α we can make N 1 to satisfy the relic density by adjusting the Z BL mass.Hence we get a relatively wider band of allowed values of g BL for a fixed α.For the chosen mass range of sterile neutrino (i.e.O(MeV)), from Fig. 7 we find that the maximum allowed value of α is ∼ 10 −10 rad.Such a sterile neutrino is free from all the constraints arising from X-ray and BBN as seen from Fig. 1 of Ref. [24].The allowed region in M Z BL − g BL plane is shown in Fig. 8. Like the previous plot in Fig. 7, here also all the points in M Z BL − g BL plane produce the correct relic density of N 1 (0.1172 ≤ Ω N 1 h 2 ≤ 0.1226).In this case we have also varied the other relevant parameters (M N 1 , α) in the range given in Eq. (31).From this figure it is seen that the allowed values of g BL increases with M Z BL .This nature of M Z BL − g BL plane can be explained in the following way.We know that the contribution of Z BL to Y N 1 depends on both Γ Z BL →N 1 N 1 and Y Z BL (see Eq. ( 28)) where the latter quantity increases with Γ H→Z BL Z BL (Eq.( 26)) as the extra gauge bosons Z BL are produced  mainly from the decay of H. However the decay width Γ H→Z BL Z BL is suppressed by M −2 Z BL and thereby reducing the comoving number density of Z BL with its mass (see Eq. ( 55) and Fig. 5(d)).In order to keep the contribution to Ω N 1 h 2 arising from Z BL decay unaltered, this decrement in Y Z BL must be compensated by a corresponding increment in Γ Z BL →N 1 N 1 which is proportional to g 2 BL (Eq.( 35)).Hence with an increase in M Z BL , g BL should also increase to satisfy the relic density constraint.
Simulations using the standard ΛCDM cosmology requires that most of the dark matter candidates should be cold to satisfy constraints from the structure formation [69,70].In our case, to get an idea about the coldness of the sterile neutrino dark matter we try to compute its free-steaming length defined by [63]: where t in is the initial time, t 0 is the present time, v(t) is the mean velocity of the dark matter, and a(t) is the scale factor of the Universe.Following Ref. [43], the hot, cold and warm dark matters are classified as: For the case where both W ± and Z BL contribute equally to the final dark matter relic density (Fig. 6b), we have calculated the free-streaming length for dark matter produced from W ± as well as Z BL decay separately.In both the cases we find that λ W ± fs ∼ 0.003 Mpc and λ Z BL fs ∼ 0.007 Mpc.Hence following the above classification of hot, cold and warm dark matter, we find that all of our sterile neutrino is cold and thus satisfying the structure formation constraints.

A possible way of detecting the sterile neutrino Dark Matter
In 2003 INTEGRAL/SPI [50] of ESA observed an emission line at an energy of 511 keV mostly from the galactic bulge.Recently, it has been reported that the measured flux from the galactic bulge by INTEGRAL/SPI is Φ exp 511 = (0.96 ± 0.07) × 10 −3 ph cm −2 s −1 at 56σ significance [71].A possible source of this line is assumed to be the annihilation of electron and positron in the galactic core.Inspite of some astrophysical processes explaining the origin of the line [72], the sources of the galactic positrons are not clear yet.Hence a series of possible explanations have been reported in last ten years involving positrons originating from a decaying [73,74] or annihilating [75,76] dark matter.For a brief review of earlier works trying to explain 511 keV line see [77].Recently the authors of Ref. [78] have shown that the explanation of this anomalous emission line is not possible from the annihilation of thermal dark matter (WIMP) due to conflict with the latest cosmological data and they have preferred a non-thermal origin of dark matter for explaining this long standing puzzle.Earlier people have tried to explain this INTEGRAL anomaly from the decay of light sterile neutrino dark matter [79,49].Here, in the case of sterile neutrino dark matter in a non-thermal setting, we have also found that such an explanation is indeed possible.The decaying dark matter scenarios however require a more cuspy density profile than the annihilation models [80].The seed mechanism behind this 511 keV emission line is the decay of sterile neutrino (N 1 ) into a e ± pair and an active neutrino.The e ± pair thus produced, get slowed down to non relativistic velocities due to several energy loss mechanisms within the galactic bulge [81] and thereby producing 511 keV gamma-line from their pair annihilation.The mass of the sterile neutrino favourable to explain this signal is ∼ 1-10 MeV [81].In the present U(1) B−L model, there are six possible Feynman diagrams contributing to this three body decay of which those which are mediated by Z BL , h and H are sub dominant due to the suppression by the very low value of U(1) B−L gauge coupling g BL .Therefore, we have used the remaining three diagrams i.e. those mediated by Z and W ± bosons to calculate the three body decay width.The expression of matrix amplitude squared and corresponding decay width Γ N 1 →e ± ν is given in Appendix A.4.A more simpler analytical expression (using some approximation) for this three body decay width can be found in Ref. [79].The expression for the gamma ray flux obtained from the galactic bulge due the decay N 1 → e ± ν is given by [79]: Here, Γ N 1 →e ± ν is the decay width of N 1 → e ± ν and ρ DM (s, ∆Ω) is the dark matter density profile in the galaxy.During our analysis, we have taken Einasto profile [82] with α einasto = 0.17 for the computation of gamma-ray flux.The angular integration over the solid angle ∆Ω is performed within the 2 • angular resolution of the spectrometer while the spacial integration is over the line of sight (l.o.s) distance of galactic bulge from the position of solar system.The extra factor of 2 appearing in Eq. ( 32) is due to the production of two photons per decay of N 1 .Using Eq. (32) we have computed the photon flux for different values of M N 1 and α.In Fig. (9) the red band shows the correct combination of M N 1 and α which is needed to explain the INTEGRAL observed flux.The dark cyan region is for those values of M N 1 and α which satisfies only the relic density constraint of N 1 .From Fig. 9, we can see that in the chosen range of M N 1 (∼ 1 -10 MeV) the active-sterile mixing angle α required to explain both the relic density as well as the INTEGRAL anomaly is ∼ 10 −12 − 10 −14 rad.

Conclusion
In this work we have shown that non-thermal sterile neutrino in U(1) B−L model can be a viable dark matter candidate.But the formalism developed here is in general applicable to any U(1) X extension of the Standard Model.Any such model trying to describe a non-thermal dark matter scenario (through IR Freeze-in) will have in general a very weakly coupled Z ′ as well as a feebly coupled dark matter candidate.Under such circumstances (i.e. when the mother particle responsible for most of the production of the dark matter has itself gone out of thermal equilibrium), we have shown how to solve a set of coupled Boltzmann equations to calculate the final relic density.We have seen that the sterile neutrinos are mostly produced from the decay of Z BL and W ± .We have also shown that though the contribution from W ± was neglected in the previous works (under the assumption of smallness of the active-sterile mixing angle), it can actually be sizeable (even dominating over the production from Z BL decay for some values of α and g BL ) depending on the parameter space we are focussing on.Note that for generic values of α we use in this work, one of the active neutrinos has to be very light.However this is allowed by the data available from the various present day neutrino experiments.Finally for completeness, we have also checked that such an O(MeV) mass non-thermal sterile neutrino can explain the 511 keV line observed by INTEGRAL/SPI.The α required to explain this signal falls in the region where the dark matter production is mostly dominated by Z BL decay.The decay width required has a corresponding life time ∼ 10 25 s, which much larger than the present age of the Universe (∼ 10 17 s).

Acknowledgement
Authors would like to thank Raj Gandhi, Alexander Merle, Sourav Mitra, Tirthankar Roy Choudhury, Osamu Seto and Bibhushan Shakya for many useful suggestions and discussions.Authors would also like to thank Mehedi Masud, Avirup Shaw and Taushif Ahmed for helping out with various aspects of numerical calculations.Authors also acknowledge Department of Atomic Energy (DAE), Govt. of India for their financial assistance and the cluster computing facility at HRI (http://cluster.hri.res.in).

A Analytical expressions for cross sections and decay widths
A.1 Production processes of N 1 from the decays of SM and BSM particles In this section, we give the expressions of all the relevant decay widths which are needed to solve the coupled Boltzmann equation for N 1 (Eqs.( 28)).
where M x denotes the mass of particle x, α is the active-sterile mixing angle of first generation (i.e.mixing of ν 1 with N 1 ) while g HN 1 N 1 and g hN 1 N 1 are vertex factors corresponding to the vertices HN 1 N 1 and hN 1 N 1 respectively.With respect to our chosen set of independent parameters, these vertex factors are given as: A.2 Production processes of N 1 from annihilation In this section, we present the expressions of all the relevant annihilation cross sections i.e. the production processes of N 1 through the annihilations of SM as well as BSM particles.In all the expressions given below, M X and Γ X denote the mass and total decay width of the particle X while g ijk denotes the coupling of the vertex involving fields i, j, k.Further, √ s is the center of mass energy of a particular annihilation process.All the annihilation cross sections given below are written in terms of our chosen set of independent parameters.
In this annihilation process, three s-channel diagrams mediated by h, H and Z and one electron mediated t-channel diagram are possible.However, Z boson and electron mediated diagrams are suppressed respectively by the fourth and second power of the active-sterile mixing angle α.Therefore we have considered other two s-channel diagrams only.
There are two s-channel diagrams and two t-channel diagrams for Z Z → N 1 N1 annihilation process.The t-channel diagrams mediated by active and sterile neutrinos are suppressed by fourth and eighth power of α respectively.Hence we have considered only two s-channel diagrams mediated by h and H.
f f → N 1 N1 (where f denotes any SM quarks or leptons) In this annihilation process four s-channel diagrams, mediated by Z, Z BL , h and H, are possible.However, the Z boson mediated diagram is suppressed by α 4 and consequently we have neglected it.
where n c is the colour charge of the corresponding fermion (f ).
Although, the annihilation processes HH → N 1 N1 , hh → N 1 N1 consist of two t-channel and two s-channel diagrams, we consider only the two dominant s-channel diagrams mediated by H and h.
This annihilation process is also mediated by two s-channel diagrams and two t-channel diagrams.However the t-channel diagram mediated by active neutrino is suppressed by α 4 .Therefore we have considered two s-channel diagrams and one t-channel diagram mediated by H, h and N 1 respectively.For simplicity, due to smallness of α, we have considered cos α ≃ 1 and consequently sin α ≃ 0 in the following expressions (Eq.( 44)-( 51)).
N x Nx → N 1 N1 (where N x s are the other two sterile neutrinos with x = 2, 3) In this annihilation process we have considered the only Z BL mediated s-channel diagram as it contributes dominantly over the other possible diagrams.
where, δV = 2 for W boson and 1 for Z boson and g HNxNx = 2 cos α x sin θ sin α x M Nx m νx v − cos θ cos α x g BL M Nx M Z BL , g Hνxνx = 2 sin α x sin θ cos α x M Nx m νx v + cos θ sin α x g BL M Nx M Z BL .
Total Decay width of h A.4 Decay width of N 1 → e ± ν i (i = 1 to 3) In this section, we have calculated the three body decay width of sterile neutrino dark matter N 1 into e ± and ν i .In this calculation we have considered only W ± and Z bosons mediated diagrams as these three diagrams contribute dominantly to this decay process of N 1 .Also in our calculation, for simplicity, we have neglected terms involving neutrino masses as these are extremely tiny compared to the masses of other particles.Further, we have neglected the intergenerational mixing between active and sterile neutrinos.We have define two quantities X and Y in terms of four momentums of ν i , e + and N 1 as where P is four momentum of N 1 while that of ν i and e − are p 1 and p 2 respectively.Now, and with Γ W and Γ Z are the total decay widths of W ± and Z bosons respectively.Therefore, where U e α s are the elements of PMNS matrix of neutrino mixing and in terms of neutrino mixing angles these are defined as (assuming Dirac CP phase δ = 0) U e 1 = cos θ 12 cos θ 13 , U e 2 = sin θ 12 cos θ 13 , U e 3 = sin θ 13 .
Finally, the expression of Matrix amplitude square for the process N 1 → e ± ν i is given by, The corresponding decay width in terms of |M 2 | is given by, The upper and lower limits of the quantities X and Y are given below with and (74) Y and U(1) B−L are denoted by g, g ′ and g BL .The corresponding gauge fields are W a µ (a = 1, 2, 3), B µ and Z BLµ .After spontaneous breaking of SU(2) L × U(1) Y × U(1) B−L symmetry by the VEVs of Φ and Ψ we get two physical neutral scalar fields h and H which can be expressed as a linear combinations of φ and ψ in the following way h H = cos θ − sin θ sin θ cos θ φ ψ ,

Figure 1 :
Figure 1: Evolution of comoving number density of Z BL with respect to z.

Figure 2 :
Figure 2: Evolution of comoving number densities of Z BL and N 1 .

Figure 3 :
Figure 3: Dependence of Y Z BL and Y N 1 on different sets of initial temperatures.

Figure 4 :
Figure 4: Comparison of comoving number densities of Z BL and N 1 with respect to different sets of chosen model parameters.

Fig. 5 (
(a)-(d)) we show the variation Y N 1 and Y Z BL with respect to different sets of independent parameters as before.Interesting feature of Fig. 5((a)-(d)) when contrasted with Fig. 4((a)-(d))

Figure 5 :
Figure 5: Comparison of Z BL and N 1 comoving number density with respect to different sets of chosen parameters.

Figure 6 :
Figure 6: Relic abundance of N 1 as function of z along with the relative contributions of W ± and Z BL decay channels.All the plots are drawn for M N 1 = 1 MeV.

Figure 7 :
Figure 7: Allowed region in g BL Vs α plane satisfying the relic density criteria.

Figure 8 :
Figure 8: Allowed region in M Z BL Vs g BL plane satisfying the relic density criteria.

Figure 9 :
Figure 9: Values of M N 1 and active-sterile mixing angle α allowed by the relic density of N 1 (0.1172 ≤ Ω N 1 h 2 ≤ 0.1226), are shown by the dark cyan points.The points lying within the red coloured band reproduced the flux observed by INTEGRAL/SPI.