A model for mixed warm and hot right-handed neutrino dark matter

We discuss a model where a mixed warm and hot keV neutrino dark matter rises naturally. We arrange active and sterile neutrinos in the same $SU(3)_L$ multiplet, with the lightest sterile neutrino being dark matter. The other two heavy sterile neutrinos, through their out-of-equilibrium decay, contribute both to the dilution of dark matter density and its population, after freeze-out. We show that this model features all ingredients to overcome the overproduction of keV neutrino dark matter, and explore the phenomenological implications for Big Bang Nucleosynthesis and the number of relativistic degrees of freedom.


Introduction
Although there is no doubt about the existence of dark matter (DM) [1,2], we have no idea about its nature. There are compelling pieces of evidence that dark matter may be composed of elementary particles, all based on its gravitational interaction with ordinary matter. Such particles must be electrically neutral (at least effectively) and cosmologically stable. Dark matter is also crucial for evolution of structure formation as we observe today. In general, dark matter candidates are classified as hot dark matter (HDM), warm dark matter (WDM), or cold dark matter (CDM) depending on their free-streaming around the period of structure formation. Structure formation requirements do not allow for dark matter to be comprised of mostly HDM [3].
Weakly interacting massive particles (WIMPs) are by far the most extensively studied class of CDM as the correct dark matter abundance is easily reproduced with cross sections around the weak scale [4]. They have been extensively searched for by many experiments (direct and indirect detection, and colliders) with no success [5]. The null results reported thus far motivate us to explore alternative candidates. As CDM faces problems at small-scale astrophysical scales, mixed populations of dark matter are well motivated [6]. Right-handed neutrinos (from now on sterile neutrinos) with mass around keV scale are suitable candidates. Sterile neutrinos arise in many popular extensions of the standard model (SM) such as leftright model [7][8][9], B-L model [10][11][12][13][14] and in a particular version of gauge models based in the SU (3) C × SU (3) L × U (1) N symmetry [15][16][17].
Although interesting alternatives to CDM candidates, keV sterile neutrinos are usually overproduced in simplified models [18]. In this case, the most plausible way to dilute this dark matter population and obtain the correct abundance is through entropy injection [19]. In general, entropy can be injected into the early universe when a long-lived particle, a diluton, that decouples while relativistic, dominates the energy density of the universe and decays once non-relativistic [19]. Thus, a successful keV dark matter model should feature a longlived particle that plays the role of such diluton. The natural diluton candidates are the right-handed neutrino themselves. The lightest right-handed neutrino is stable, while the other two act as the diluton [20][21][22]. This apparent easy solution faces solid constraints from Big Bang Nucleosynthesis (BBN) and the Cosmic Microwave Background (CMB).
The study of keV neutrino dark matter has been discussed elsewhere in simplified models, our goal here is to embed this mechanism in a UV complete model, which is well motivated for other theoretical reasons as family replication [23] and electric charge quantization [24,25]. We will discuss this keV neutrino dark matter in a model based on the SU (3) C × SU (3) L × U (1) N gauge symmetry, 331 for short. There are several ways to arrange the fermion generations in this gauge symmetry. The different ways give rise to different models. Here we will focus on the 331ν R [15][16][17], which features right-handed neutrinos in the same SU (3) L multiplet where resides the active neutrinos [26].
In this work, we check under which conditions the 331ν R accommodates a successful keV dark matter candidate. To do so, we invoke a discrete symmetry to guarantee the stability of the lightest sterile neutrino. We calculate the relativistic freeze-out of the lightest sterile neutrino and consider the heavier ones as dilutons. As we are dealing with an extended gauge sector there are new interactions that impact the sterile neutrino abundance, differing from past keV neutrino dark matter studies. We highlight that the two heavy sterile neutrinos dilute and produce dark matter. To test our model, we estimate the free-streaming length of our dark matter candidate and find that we have a mixed WDM+HDM scenario, which successfully obeys the constraints stemming from CMB and BBN.
This work is organized as follows. In Section 2 we present the key aspects of the model; in Section 3 and Section 4 we address the production mechanisms; in Section 5 we outline the viable parameter space; lastly in we draw our conclusions in Section 6.
2 The essence of the 331ν R model

Particle content
In the 331ν R model the lepton generations are arranged as, where l = e, µ, τ refers to the three generations.
In the quark sector, anomaly cancellation requires the first two generations of quarks to come in the anti-triplet while the third generation in triplet representation as, The new quarks (u , d ) have the usual electric charges. The gauge sector of the model is formed by the standard gauge bosons, A, W ± and Z, and five others called W ± , U 0 , U 0 † and Z . The interactions of these gauge bosons with matter can be found in Ref. [27].
The scalar sector of the original version involves three scalar triplets, namely With such scalar content, when the symmetry SU (3) C ×SU (3) L ×U (1) N is spontaneously broken to SU (3) C × U (1) em , all particles acquire masses at tree level, except the neutrinos.
We adopt the following vacuum structure, To further simplify the model, we assume the following discrete symmetry transformation over the full Lagrangian, ( η, ρ, e lR , u aR , d aR ) → − ( η, ρ, e lR , u aR , d aR ) , where a = 1, 2, 3. This discrete symmetry (Z 2 ) will play an important role in the dark sector of the model, as we shall see. With such matter and scalar content we build the following Yukawa interactions invariant under the gauge symmetry, (2.5) which generate masses for all fermions, with the exception of neutrinos. We remark that while the SM does not contain any dark matter candidate, the 331ν R model poses three candidates, namely, η 0 , U 0 or ν R , which are mutually exclusive. U 0 is underabundant, while η 0 and ν R are viable multi-TeV dark matter candidates [28][29][30][31][32][33][34][35][36][37]. Concerning U 0 , despite of being an interesting candidate, unfortunately it does not provide the correct abundance (it is under-abundant), while ν R did not receive any attention until now. In other words, this is the first time that right-handed neutrino is being treated as dark matter in the 331ν R .
Before considering ν R as dark matter candidate we discuss in the next section how to generate masses for the neutrinos in the model.

Neutrino Masses
Right-handed neutrinos are hypothetical particles and their masses are free parameters that may take a wide range of values varying from eV up to GUT scale. In the original version of the 331ν R model neither ν L nor ν R gain masses. The most immediate way of providing masses for them is through effective dimension-5 operators [26]. In the case of left-handed neutrinos, this operator is constructed with the scalar triplet η and the lepton triplet L, According to this operator, when η 0 develops a VEV, v η , the left-handed neutrinos develop Majorana mass terms, Regarding right-handed neutrinos, the dimension-5 operator that give them mass is constructed with the scalar triplet χ and the lepton triplet L, When χ 0 develops a VEV, v χ , this effective operator provides Majorana masses for the right-handed neutrinos, Once v χ > v η , then m ν R > m ν L . Thus, light right-handed neutrinos is a natural result of the model. Observe that the discrete symmetry discussed above avoids the operator ∼ 1 Λ (Lη)(χL) which would generate mixing among active and sterile neutrinos. This is a particularly interesting result, as it renders the lightest right-handed neutrino automatically stable and suitable to be our dark matter candidate. The ways these operators can be realized are not relevant for the issue we are going to address here. However, for the sake of completeness, by invoking the existence of a sextet of scalar with mass belonging to the GUT scale, we can realize such operators according to type II seesaw mechanism [38], as done in Refs. [39][40][41].

Main interactions
To work with the physical neutrinos we have to diagonalize m ν L and m ν R . From now on we refer to the active physical neutrinos as ν L = (ν 1 , ν 2 , ν 3 ) T L , and to the physical sterile ones as To simplify even more we assume ν R in a diagonal basis. We now present the interactions involving sterile neutrinos that matter for us here: where C W = cos(θ W ) and S W = sin(θ W ) with θ W being the Weinberg angle and l = (e , µ , τ ) T . The dominant interactions involving quarks that matter for the calculation of the dark matter abundance are

Relic abundance of a light sterile neutrino
We assume that N 1 is the lightest of the sterile neutrinos, being in principle a dark matter candidate. We therefore check under which conditions it thermalizes with species of the standard model bath in the early universe and calculate its relic abundance accordingly.
In the context of the 331ν R model, the sterile neutrinos are able to thermalize with the standard fermions through exchanges of W , Z and U 0 , as depicted in Fig. 1. This happens whenever their interaction rates Γ N i (T ) are faster than the Hubble rate H(T ) at a given temperature T : The rate at which N i self-annihilate into species 3 and 4, with masses m N 3 and m N 4 , is given by Γ N i (T ) = n eq N i (T ) σv , with n eq N i (T ) their equilibrium number density and a thermally averaged annihilation cross-section given by 2) where S is the symmetrization factor, s is the Mandelstam variable, λ (x, y, z) is the Källen function, K i is the modified Bessel function of the second kind of order i, Ω is the solid angle between initial and final states in the center of mass frame, and |M | 2 the (not averaged) squared amplitude of the process.
As usual, we compute the freeze-out temperatures T f at which the sterile neutrinos decouple from the thermal bath by equaling n eq N i (T f ) σv (T f ) = H(T f ) for the main processes. For simplicity but without loss of generality for our purposes, we will assume the hierarchy m N 3 m N 2 m N 1 . As a consequence, we can neglect co-annihilation processes [42]. In Fig. 2 we show the ratio between the annihilation rates of N 1 (in yellow) and of N 2 (in red) and the Hubble rate in the standard case of a radiation-dominated era as a function of the inverse of temperature 1 for m W = 1.0 × 10 6 GeV (left panel) and m W = 1.0 × 10 10 GeV (right panel). To obtain the curves of Fig. 2 we consider all annihilation processes shown in Fig. 1. This analysis is also applicable to the case of N 3 , but for the purpose of this paper we do not need to take it into account, as we will see later. As shown in Fig. 2, the freeze-out temperature T f increases with the mass of the mediator. In our model, the gauge bosons have similar masses, and from now on we assume m W = m U 0 0.72 × m Z [27]. We can see in Fig. 2 that the maximum of the ratio happens when the gauge bosons are produced on-shell, making the decoupling to take place close to that scale.
We have chosen in Fig. 2 masses of N 1 , N 2 , and W close to the values which will provide the right amount of relic abundance for N 1 , as we will show later. In such a case, we can see that N 1 and N 2 freeze-out almost at the same time (when n N i σv ∼ H, as indicated by the dashed horizontal line), with T f 1.0 × 10 3 GeV for m W = 1.0 × 10 6 GeV (left panel) and T f 2.2 × 10 8 GeV for m W = 1.0 × 10 10 GeV (right panel). Thus, both N 1 and N 2 decouple while still relativistic. As shown in Fig. 2 the freeze-out temperature T f strongly depends of mediator mass, for that benchmark value of mass. When the masses of N 1 and N 2 increase (indicated by the yellow and red dashed curves in Fig. 2, respectively) the Boltzmann suppression (which occur when m N i ∼ T ) occurs earlier, as shown in Fig. 2. Given that in our model the couplings between sterile neutrinos, gauge bosons, and standard fields are all of order O(10 −1 ), the only way of avoiding thermalization of sterile neutrinos is by invoking gauge bosons heavier than about O(10 16 ) GeV.
We remark that among all processes contributing to the freeze-out of N 1 and N 2 , the most relevant are annihilations into SM charged leptons and active neutrinos. In the relevant limit of m N i T m W , they happen at the following rate: We have found the typical freeze-out temperature to be given by where g e (T f ) is the number of energetic degrees of freedom at freeze-out, with g e (T f ) ∼ 100 for T f > 100 GeV. We have checked that our numerical solution shown in Fig. 2 is in a good agreement with the estimation above.
We can therefore conclude that N 1 would be a cosmic relic which was once thermalized, and proceed with the computation of their final abundance in order to determine whether it can constitute the cosmological dark matter.

Relativistic freeze-out
The final relic abundance of N 1 is defined to be where the label " 0 " indicates quantities as measured today, with Ω 0 DM h 2 0.12 being inferred by the Planck satellite [2], and Y i ≡ n i /s is the yield of a species i, with s the entropy density in a comoving volume.
In this work, we are interested in the case of a light sterile neutrino dark matter. For N 1 thermally produced the lower limit on its mass is m N 1 2.0 keV [44]. As we have just seen, for the mass hierarchy m N 1 , m N 2 m W , both N 1 and N 2 decouple from the thermal bath almost simultaneously and at high temperatures (T ∼ m W /100). The yield of a Majorana neutrino that decouples while ultra-relativistic is given by , the agreement with the relic abundance constraint, would require an unreasonable amount of relativistic degrees of freedom by the time of freezeout in order for a light sterile neutrino to not overclose the universe.
To the best of our knowledge, the only way of depleting the yield of a decoupled species is by considering entropy production after freeze-out. As it is well known [45], a long-lived particle that decoupled while ultra-relativistic can dominate the cosmic expansion before decaying, thus injecting a sizable amount of entropy into the thermal bath. It is therefore interesting to notice that the heavier sterile neutrinos are natural candidates to deplete the yield of N 1 . For simplicity, we will investigate the out-of-equilibrium decay of N 2 as the source of entropy production, while assuming that N 3 is heavy enough as to not affect our analysis.
It is straightforward to see that an increase of total entropy S in a comoving volume after the freeze-out of any relic, by a factor of ∆ ∼ S(T 0 )/S(T f ), will dilute its yield by the same factor: For freeze-out happening at TeV scale, our model provides g s (T f ) ∼ 100. Thus, N 1 with mass in the keV-MeV range requires an entropy injection of ∆ = 10 − 100 in order to be a viable dark matter candidate. As we show in what follows, though, the out-of-equilibrium decay of N 2 contributes in a non-trivial way to the final abundance of N 1 in the context of the 331ν R model.

Non-thermal production
As we have just discussed, the out-of-equilibrium decay of N 2 into species of the thermal bath dilutes the abundance of our dark matter candidate, N 1 . However, in the 331ν R model, treebody decays of N 2 into N 1 can be sizable, which could potentially repopulate (and overclose) the universe with dark matter. Therefore, the final relic abundance of N 1 will have a (thermal) contribution from the relativistic freeze-out and also a (non-thermal) contribution from the tree-body decays of N 2 into N 1 : Let us parametrize the total decay width of N 2 in terms of the partial width into N 1 , Γ (N 1 ) The dimensionless parameter N 2 contains all other channels which do not involve N 1 as final product 2 . Of course, N 2 must decay into species which thermalize with the SM bath in order to dilute N 1 . Decay channels into 331ν R states do not necessarily thermalize, but here we will treat α as a free parameter encoding only decay channels which instantaneously thermalize with the SM bath.
We have found that the leading contributions to Γ (N 1 ) N 2 are the three-body decays into charged leptons and neutrinos. In the limit m N 1 , m l m N 2 m W , we have (3.10) In the next section, we develop the tools needed to properly dealing with the competing effects of the N 2 out-of-equilibrium decays.

Coupled evolution of sterile neutrinos
Let us now discuss how to properly find the final relic abundance of our dark matter candidate, the keV scale sterile neutrino N 1 . The relativistic freeze-out of both N 1 and N 2 , as well as the non-thermal production and dilution of N 1 due to the out-of-equilibrium decay of N 2 , can be taken into account by solving the following coupled Boltzmann fluid equations for the yields of N 1 and N 2 3 : where the scale factor a is used as a time parameter.
The relativistic freeze-out and the non-thermal production of N 1 are accounted for by the first term in the right hand side of the equation above, whereas the second term accounts for the dilution of Y N 1 after the entropy production.
The reaction rate densities R N 1 ,N 2 contain all processes that can change the number of N 1 and N 2 in a comoving volume. We have found the following leading contributions: The terms proportional to σv N i N i represent the annihilations into SM leptons (see Eq. (3.3)) and their backreactions. We can represent σv N i N i as the sum of the channels that contribute for it: The other terms represent the contribution of the N 2 decays and inverse decays. We recall that α encodes the channels without N 1 as a decay product, assuming that they all thermalize instantaneously. We therefore see that the decay of N 2 into N 1 couples their evolution in the early universe, even if the entropy production were negligible.
In order to inject a significant amount of entropy into the thermal bath after decaying, N 2 must be significantly long-lived and have dominated the total energy density of the universe. This is indeed a natural consequence of our framework, since N 2 decouples while ultra-relativistic -which means that ρ N 2 /ρ R ∝ m N 2 a once it becomes non-relativistic, with ρ R the energy density of radiation. The Hubble rate in Eq. (4.1) will be therefore given by where M P l 2.4 × 10 18 GeV is the reduced Planck mass. From the temperature at which N 2 starts dominating the energy density, T i , until its complete decay, at the so-defined reheat temperature T RH , the universe would have therefore undergone an early matter-dominated era. The duration of such an era is determined by the amount of entropy produced, T i /T RH ∝ ∆, and the reheat temperature is found to be given by [46] T RH = 5π 2 72 g e (T RH ) In order to not jeopardize the BBN predictions [47], we must ensure T RH 4 MeV. This guarantees that N 2 decays before the weak decoupling of active neutrinos and all the standard leptons thermalize.
The rate of injection of entropy due to the decay of N 2 is given by [45], where f T represents the fraction of the decay products of N 2 that thermalize in the plasma [48,49]. The fraction f N T of decay products that do not thermalize will populate the sea of decoupled relativistic species (contributing to ∆N ef f = 0, see Section 5.1). It is therefore convenient to rewrite the energy density of N 2 in terms of f T and f N T : where Br(N 2 → others) = 1/(1 + 1/α) (see Eq. (3.9)) is the branching ratio into all decay channels without N 1 in final states. Since T RH < T f , the N 1 produced via decay will not be able to thermalize anymore. This is why such dark matter population is said to be non-thermal. On the other hand, above 4 MeV, all SM leptons are able to thermalize. It is then easy to see that f l T = f ν T = 2/3, whereas f l N T = f ν N T = 1/3, so that (4.8) Finally, since ρ N 2 and ρ R evolve non-trivially during the evolution of N 1 and N 2 , the set of Eq. (4.1) must be solved together with the following Boltzmann fluid equations: (4.9)

Numerical results
We numerically solve the set of equations (4.1), (4.4), (4.6), and (4.9). For numerical convenience, we re-scale the scale factor by A ≡ aT RH , the energy density of N 2 by Φ N 2 ≡ ρ N 2 a 4 , and the energy density of radiation by Φ R ≡ ρ R a 4 . Regarding the initial conditions, at A = A I 1 (the actual value does not change results), the yields of N 1 and N 2 follow their equilibrium values. The freeze-out of the sterile neutrinos take place while they are still relativistic, during the radiation era. From A I until the moment when N 2 becomes non-relativistic, at A = A N R , there is no entropy production and we just need to solve the coupled set of Eq. (4.1) without the last terms. In this case the Hubble rate is the usual one, H ∝ T 2 /M P l .
For A A N R , we consider the full set of equations. We follow Ref. [46] and assume an inflationary model that yield an inflationary reheat temperature of T Inf RH 7 × 10 15 GeV, which implies S I = S N R 897 and Φ I R = Φ N R R 1790. At A N R , both Y N 1 and Y N 2 are given by Eq. (3.6). Since at this point N 2 is non-relativistic, we have (4.10) In Fig. 3 we present the evolution of the set of equations (4.1), (4.6), and (4.9) for α = 0, 10 and 100 (continuous, dashed and dotted curves, respectively). On the left panel of Fig. 3 we show the full solutions for Y N 1 (blue curve) and Y N 2 (red curve), as well as the solution for Y N 1 in the absence of non-thermal contribution, Y T hermal N 1 (yellow curve). On the right panel we show the solutions for entropy S (pink curve) and for the quantities ρ N 2 a 4 (red curve) and ρ R a 4 (green curve) and we can observe that when N 2 becomes non-relativistic (at A N R , as shown) the energy density of radiation is still greater than the energy density of N 2 . However the ratio between the energy density of non-relativistic N 2 and the energy density of radiation evolves like ρ N 2 /ρ R ∝ am N 2 , then ρ N 2 can dominate the energy density of the universe if N 2 is sufficiently long lived, as shown in Fig. 3. The complete decay of N 2 correspond in Fig. 3 to the abrupt decrease of Y N 2 on the left panel (or ρ N 2 a 4 on the right panel). (yellow) and Y N 2 (red). Right panel: Curves of entropy S (pink), ρ N 2 × a 4 (red) and ρ R × a 4 (green), also shown are how the value of ∆ changes as we increase α (∆ = 100, 41 and 14, for respective α).
We can observe in the right panel of Fig. 3 that when N 2 decays completely, the entropy is increased by a factor ∆ and, as expected, it keeps constant before and after the decaying of N 2 . The increase of entropy dilutes the abundance of N 1 , as we have discussed, and we can observe this behavior from the evolution of the Y N 1 (blue curves) on the left panel in the Fig. 3. When N 2 decays completely, the injection of entropy ceases and Y N 1 , S and ρ R a 4 levels off.
As shown in Fig. 3 the free parameter α plays an important role in the dilution of N 1 . Naively, one would expect that the dilution would increase with α, since more thermalized decay channels are allowed. However, since Γ N 2 increases with α, increasing α makes N 2 to decay earlier, such that it does not dominate the evolution of the universe long enough for a significant entropy injection to take place. We can observe this behavior on the right panel of Fig. 3: when α = 0 (solid curves), ρ N 2 > ρ R for a much longer period and inject more entropy (∆ = 100) than when we take α = 10 2 (dotted curves) and obtain ∆ = 14. As the injection of entropy (parameterized by ∆) is responsible for diluting N 1 , this explains why the final value of Y N 1 is higher as we increase α.
The free parameter α has also an important role on the non-thermal contribution to the abundance of N 1 . As shown in the left panel of Fig. 3 there is a gap between Y N 1 (thermal and non-thermal contribution) and Y T hermal N 1 (only thermal contribution) which means that the N 1 produced via N 2 decays is responsible for increasing the abundance of N 1 from Y T hermal N 1 to Y N 1 . We can observe that this gap decreases as we increase α. This behavior is explained by the fact that when α increases, it allows N 2 to decay in another particles beyond N 1 .
According to our computation, for α = 10 2 we have that the non-thermal contribution is approximately 1%, while for α = 0 the non-thermal contribution is approximately 50%.

Viable parameter space
In this section we obtain the constraints on the parameter space (m W , m N 2 ), that makes N 1 produced by means of relativistic freeze-out and diluted from N 2 decay a realistic DM candidate. In Fig. 4 we present these constraints for α = 0 (solid curves) and α = 100 (dashed curves).
Our first constraint is on the reheating temperature, given in Eq. (4.5), which is the temperature soon after the decay of N 2 . As we have pointed out, N 2 must decay prior to the active neutrinos decoupling as to not disturb the BBN predictions, such that we need to ensure T RH 4 MeV. In Fig. 4 the region in pink is excluded because it gives T RH < 4 MeV. As T RH ∝ Γ N 2 and Γ N 2 increases with α, this bound, translated to an upper limit on M W , is weakened as we increase α.
As we discussed above, N 2 must decouple ultra-relativistic in order to produce enough entropy and sufficiently dilute N 1 . This criterion rule out the green-shaded region in Fig. 4, in which m N 1 T f with T f given by Eq. (3.4). This bound is independent of the parameter α. Finally, the LHC constraint over the mediator mass, m W 3 T eV [50][51][52][53], is indicated by the blue-shaded region.
The correct relic abundance of N 1 today, Ω N 1 h 2 0.12, as constrained by the Planck satellite [2], is ensured by demanding In order to take into account the thermal and non-thermal contributions, we find Y 0 N 1 by solving Eq. (4.1) numerically. For this we developed a Python algorithm which searches the values of m N 2 and m W obeying Eq. (5.1). The contours of correct relic density are shown for m N 1 = 2keV (light-grey curve) and m N 1 = 1MeV (grey curve), for α = 0 (continuous) and α = 100 (dashed). The region below the light-grey (grey) contours are excluded by Planck as they overclose the universe. As the purpose of this paper is to study N 1 as a viable DM candidate, we are interested in the light-grey (grey) curves themselves. As we can see, the heavier the dark matter, the heavier the mediators need to be for the agreement with Planck.

Contribution to ∆N ef f
The amount of relativistic particles at matter-radiation equality epoch contributes to the number of relativistic degrees of freedom and directly affects the CMB power spectrum. An important parameter in this regard is the effective number of neutrino species defined as N ef f = 8 7 ( 11 4 ) 4/3 ( ρ rad −ργ ργ ), where ρ rad and ρ γ are respectively the total radiation and photon energy density. The current value for N ef f from Planck 2018+BAO [2] is N ef f = 2.99 ± 0.17  which is in agreement with the standard model prediction N SM ef f = 3.0440 [54]. The change in the amount of radiation is quantified by means of the parameter ∆N ef f ≡ N ef f − N SM ef f . Since the fraction of N 1 produced non-thermally becomes non-relativistic at temperatures of the order O (eV ) (see Appendix A), it increases the value of N ef f . We remark that the N 1 produced thermally and then diluted does not contribute to ∆N ef f because it becomes non-relativistic before equality (see Appendix A).
The contribution to ∆N ef f in our model comes therefore from the non-thermal population of N 1 and is mainly controlled by the parameter α: . According to this, α = 0 yields ∆N ef f 2.05 which is excluded by the bounds discussed above. In other words, we need α > 0. Future experiments as CMB-S4 will have a sensitivity to constraint ∆N ef f = 0.060 (at 95% C.L.) [55]. In the absence of evidence for ∆N ef f = 0, we will have the limit ∆N ef f < 0.06 which requires α 33.

Structure formation (free-streaming)
As dark matter can be classified as cold (CDM), warm (WDM) or hot (HDM) according to its free-streaming λ f s [56], in this section we will roughly estimate the free-streaming of N 1 .
The free-streaming at the epoch of matter-radiation equality is an important parameter that allows us to understand how the first structures were formed. It is given by [57] where t prod is the production time, t e the equality time, and v(t) the mean velocity of N 1 .
The N 1 population thermally produced and then diluted due to an injection of entropy is "colder" than what it would be if only thermally produced. The relation between the temperature of the decoupled N 1 and that of the plasma is affected by the injection of entropy, parametrized by ∆, and is given by where the subscript "b" refers to any epoch prior to the dilution. The thermally produced fraction of a keV N 1 , which is relativistic at production, has its free-streaming suppressed to the scale of O (0.1Mpc) due the dilution ∆, being therefore classified as warm or even cold dark matter (see for instance Refs. [48,56]).
On the other hand, the free-streaming of N 1 produced non-thermally is given by [58,59] where v 0 is the initial velocity of N 1 and z e is the red-shift at equality. Since the non-thermal population of N 1 is produced at T RH through 3-body decay of N 2 , and in the limit of massless final states, we have p N 1 prod m N 2 /3. For temperatures below T RH , the universe is radiation-dominated and the temperature red-shifts as a = a 0 × T 0 /T . We have therefore For m N 1 = 2.0 keV, m N 2 = 60.0 GeV and m W = 4.7 × 10 6 GeV we obtain λ f s 2.4M pc. The non-thermal fraction of N 1 is therefore classified as HDM.
Our non-thermal component can compose up to 50% of the total abundance (for α = 0) and 1% (for α = 10 2 ). As the CMB constraint requires α > 0, we have a scenario of mixed warm-hot DM. Such a scenario has important implications for structure formation and can address some small-scale problems such as the core-cusp [60]. Warm dark matter candidates can delay the formation of structures, with lighter dark matter increasing the delay. The recent measurement of the 21-cm absortion signal due to the light of the first stars, reported by the EDGES Collaboration [61], can provide lower bounds on the mass of non-cold dark matter candidates. If N 1 is mainly warm, this observation implies roughly m N 1 3 keV [62]. A rigorous approach to the implications of mixed DM in the structure formation may also put bounds on the parameter α [63,64]. However this requires the computation of the fluctuation of the power spectrum for our scenario and is beyond the scope of this paper.

Conclusions
We investigated the possibility of having a viable mixed warm and hot keV neutrino dark matter in a model based on the SU (3) L gauge group. Active and sterile neutrinos are arranged in the same SU (3) L multiplet, with the lightest sterile neutrino being dark matter. Its abundance is set by interactions with Standard Model particles controlled by the new gauge bosons rising from the extended gauge sector. Its stability is warranted by a discrete symmetry that prevents mixing between active and sterile neutrinos, making the model safe from otherwise stringent bounds, such as X-rays.
We have shown that the sterile neutrinos easily thermalize with the Standard Model bath via exchanges of these heavy new gauge bosons, unless their masses are at the GUT scale. When the sterile neutrinos are much lighter than the gauge bosons, they decouple nearly at temperatures much larger than their masses, rendering them ultra-relativistic, and therefore overproduced, at the thermal freeze-out.
In this work, we have shown that this common issue with keV right-handed neutrino is solved within this SU (3) L gauge group. Since the heavier sterile neutrinos decouple while ultra-relativistic, and are long-lived due to the heavy gauge bosons mediating their decays, they eventually dominate the cosmic expansion after freeze-out. Their out-of-equilibrium decay into SM bath species and also into keV right handed neutrino dark matter. Hence they contribute both to the dilution of warm dark matter population and as well as to its non-thermal HDM population.
Our scenario of mixed warm-hot dark matter is amenable cosmological constraints. As this entropy injection episode should take place before BBN, constraints are derived on the masses of the heavy sterile neutrino states, and gauge bosons. We found that they must be larger than the TeV scale. Moreover, the hot dark matter can also increase N ef f in a detectable way, and for this reason we imposed α 33 to make sure that the heavy sterile neutrino decay mostly into particles of the SU (3) L gauge group other than dark matter.
In summary, we conclude that our model can successfully host a mixed warm plus hot dark matter setup in agreement with existing bounds. A N 1 temperature at the matter-radiation equality Here we will check when N 1 produced thermally and non-thermally becomes non-relativistic and if they contribute to ∆N ef f at CMB.
Following [58,65], the momentum of N 1 produced non-thermally p N T N 1 by N 2 decay gets red-shifted where a D is the scale factor at the moment of decay and p N T N 1 (T D ) is the momentum of non-thermal N 1 when it is produced. Assuming that m N 2 m N 1 , m e,µ,νe,νµ , which imply in We assume that N 1 produced via decay becomes non-relativistic at some temperature T nr (with the respective scale factor a nr ). This happen when p N T N 1 (T nr ) = m N 1 [65]. Then, from Eq. (A.3) we have We can represent the red-shift when non-thermal N 1 becomes non-relativistic as For a universe dominated by radiation the temperature red-shifting as T = T 0 × a 0 /a and using Eq. (A.4) we obtain where, As a result of Eq. (A.6) we obtain that N 1 produced non-thermally becomes non-relativistic at temperatures O (eV ) for a long range of α and then contribute to N ef f . On the other hand the red-shift when thermal N 1 becomes non-relativistic is given by: where a f is the scale factor when N 1 freeze-out at the temperature T f . The thermal N 1 momentum gets red-shifted However we should note that as the thermal decoupling of N 1 happens before the N 2 decay, the thermal N 1 will be cooler than the bath particles [48], and a f a where a and T represent the scale factor and temperature for some epoch after N 2 decays. As non-thermal N 1 is produced during the injection of entropy due decay of N 2 we did not consider this contribution for non-thermal N 1 treatment.
As the thermal N 1 becomes non-relativistic when p T N 1 ∼ m N 1 , therefore from Eq. (A.9) and Eq. (A.10), we can rewrite Eq. (A.8) as The parameter ∆ plays an important role in Z T nr . For m N 1 = 2.0 keV and for some α we only have a unique ∆ that can reproduce N 1 as DM, such as for α = 0 we obtain numerically that we need of ∆ 100 to obtain Ω N 1 h 2 = 0.12. As we discussed in Section 4.1 if we increase α we decrease the abundance of non-thermal N 1 which implies that we need to decrease ∆ in order to obtain N 1 as DM. Then we can obtain the lower value to Z T nr for m N 1 = 2.0 keV if we take a scenario of the minimum allowed ∆ value. This scenario can be reproduced if we assume that N 1 is only produced via freeze-out or if we assume that α is so large that nonthermal contribution becomes irrelevant. In these cases we can assume that the lower value to ∆ to provide the right amount of relic abundance for N 1 is ∆ 19.0, given by Eq. where we take g s (T f ) = 100. For ∆ = 19.0 we obtain Z T N T = 6.9×10 7 Z e , where Z e 3365 is the red-shift at matter-radiation equality. This is an important result because it shows us that the N 1 produced via freeze-out and consequently diluted by the decay of N 2 becomes non-relativistic before matter-radiation equality epoch.

B Evaluation of ∆N ef f
Here we derive Eq. (5.2). At temperature T 0.5 MeV only the photon (γ), SM neutrinos (ν) and non-thermal N 1 contribute to the radiation energy density of the universe. Then we can represent the energy density of radiation in that epoch as ρ R = ρ γ + N SM ν ρ ν + ρ N 1 The energy density for a ultra-relativistic particle (radiation) is given by: where g accounts for its spin degeneracy. As we discussed in Section 4, we can write the energy density of non-thermal N 1 as function of the energy density of N 2 at the time that N 1 was produced (at T RH ), such that ρ N 1 = f N T · ρ N 2 (T RH ), with f N T given by Eq. (4.8). According to [45,57], between the epoch that ρ N 2 starts to dominate the energy density of Universe until it decays (at T RH ) the radiation produced from decaying of N 2 is the dominant radiation component. Then we can assume that ρ N 2 (T RH ) ρ R (T RH ). As non-thermal N 1 does not thermalize anymore, its energy density only gets red-shifted. Then we can write the energy density of radiation as where T ν represents the temperature of SM neutrinos and T RH N 1 is the temperature of N 1 at the time of production (where T RH N 1 = T RH , after that T RH N 1 gets red-shifted). As g ν = 2, T /T ν = (11/4) 1/3 [57] and using the value of f N T given by Eq. However we would like to know the ∆N ef f at matter-radiation equality epoch T e . As the temperature of N 1 gets red-shifted after its production, we represent the temperature of N 1 at equality by T e N 1 = T RH N 1 × a RH /a e . From entropy conservation we have g s (T RH )T 3 RH a RH = g s (T e )T 3 e a 3 e , then