New Production Mechanism for keV Sterile Neutrino Dark Matter by Decays of Frozen-In Scalars

We propose a new production mechanism for keV sterile neutrino Dark Matter. In our setting, we assume the existence of a scalar singlet particle which never entered thermal equilibrium in the early Universe, since it only couples to the Standard Model fields by a really small Higgs portal interaction. For suitable values of this coupling, the scalar can undergo the so-called freeze-in process, and in this way be efficiently produced in the early Universe. These scalars can then decay into keV sterile neutrinos and produce the correct Dark Matter abundance. While similar settings in which the scalar does enter thermal equilibrium and then freezes out have been studied previously, the mechanism proposed here is new and represents a versatile extension of the known case. We perform a detailed numerical calculation of the DM production using a set of coupled Boltzmann equations, and we illustrate the successful regions in the parameter space. Our production mechanism notably can even work in models where active-sterile mixing is completely absent.


Introduction
Our picture of the whole Universe has been strengthened first by the analysis of the WMAP 9-year data set [1] in combination with the data from the ground based telescopes SPT [2] and ACT [3], which had been supplemented in March 2013 by the release of the long awaited data obtained by the Planck satellite [4]. Still we are puzzled by the ingredients of our Universe, one of the biggest mysteries being the identity of the socalled Dark Matter (DM). Even if the ΛCDM model, involving a cosmological constant Λ and cold, i.e. non-relativistic, DM (CDM), provides a very good fit to the data [4], the intermediate case of warm Dark Matter (WDM) is still a valid possibility [5,6,7,8,9,10,10,11,12,13]. However, hot (i.e., highly relativistic) DM is clearly excluded by structure formation arguments [14,15], One particularly interesting candidate particle which in most settings turns out to be WDM would be a sterile [i.e., mainly a Standard Model (SM) singlet] neutrino with a mass of a few keV. If such a particle exists, in addition to two heavier (i.e., GeV) neutrinos which are nearly degenerate in mass, the resulting setting, called the νMSM [16], can indeed simultaneously accommodate for neutrino masses, for DM, and for the baryon asymmetry of the Universe [17,18]. However, while the νMSM can successfully accommodate for such a peculiar set of sterile neutrinos, it does unfortunately not yield an explanation for the required mass pattern. This fact has triggered the construction of a variety of models in the recent years, which try to give such an explanation. The ideas used to obtain light sterile neutrinos thereby range from the application of the Froggatt-Nielsen mechanism [19,20,21], over flavour symmetries [22,23,24], extra dimensions [25,26,27], extensions of the seesaw mechanism [21,28,29,30,31,32,33,34], composite neutrinos [35,36], global symmetries [37,38], loop suppressions [39], to gravitational effects [40] -see Ref. [41] for a recent review. Even lighter sterile neutrinos have also attracted considerable interest (see, e.g., Refs. [42,43] for two up-to-date reviews), but in general right-handed neutrinos have applications at various scales [44].
A frequent "problem" with non-standard DM candidates such as keV sterile neutrinos is that they cannot be produced easily via the generic process of thermal freeze-out. This is simple to understand, since this mechanism requires particles to be in thermal equilibrium with the plasma in the early Universe, which does not work for sterile neutrinos as their interactions are too weak. Nevertheless, sterile neutrinos will in general have slight admixtures to active neutrinos. Thus, they can be produced from time to time from the thermal plasma even though they never entered thermal equilibrium. For the case of keV sterile neutrinos, this simple scenario is called the Dodelson-Widrow (DW) mechanism [45] and it is nowadays known to be excluded by observations, in case that singlet scalar. For example, if λ 10 −6 , the scalar will enter the thermal equilibrium [60]. If it is unstable but its lifetime is large enough, it could then freeze-out as thermal relic and afterwards decay to produce keV sterile neutrino DM. On the other hand, for smaller values, λ ∼ 10 −10 , the scalar could freeze-in instead and, provided that it is heavy enough and stable (or at least very long lived), itself be the DM in the Universe. This case has been studied for the scalar either being a generic heavy singlet [66] or a light pseudo Nambu-Goldstone boson [67]. However, what has not been studied up to now in the literature is the combination of the two settings, namely the freeze-in of an unstable scalar σ which subsequently produces keV sterile neutrinos via its decay. The current study will close this gap and show that such a scalar FIMP (Feebly Interacting Massive Particle [65]) can also lead to an interesting and valid possibility to produce keV sterile neutrino DM.
The paper is structured as follows: We first give an illustrative description of the idea behind the mechanism in Sec. 2. The more technical details, such as a description of the model setting and of the Boltzmann equations to solve, as well as a discussion of the relevant bounds are provided in Sec. 3. Our actual results are presented in Sec. 4, before concluding in Sec. 5. The appendices provide further technical details, such as definitions of the effective degrees of freedom (Appendix A), remarks on the use of modified Bessel functions (Appendix B), as well as detailed analytical derivations of the Boltzmann equations (Appendix C) and of the free-streaming horizon (Appendix D).

2
The basic idea: keV neutrino production by the decays of scalar FIMPs Before entering the technical details, we would like to give an illustration of how the proposed mechanism works. The decisive point of the production of any particles through "late" decays of a metastable species is that the parent particle has to be produced in the first place. While this might seem as a disadvantage, since two production stages are needed, it can in many circumstances actually be advantageous, because different constraints may hold for the two particles involved. For example, the keV sterile neutrino cannot be produced from the thermal plasma only: in case it enters thermal equilibrium, it is typically overproduced since it is relativistic at freeze-out [48,49]. If it does not enter thermal equilibrium and is only produced non-resonantly by small admixtures, its spectrum is too warm if the correct abundance is produced. The constraints from structure formation [7] can then only be realised for relative large keV sterile neutrino masses which are in conflict with the constraints from the non-observation of the decay of the keV neutrino into a light neutrino and a photon [17,18]. The singlet scalar, instead, can be produced via thermal freeze-out, and then by its decay lead to a suitable abundance of keV sterile neutrino DM, while at the same time being compatible with all bounds [59,60]. In this paper we pursue a different path to produce the singlet scalar σ: if the Higgs portal coupling is small enough, λ ≪ 10 −6 [59,60], the scalar particle never enters thermal equilibrium (due to its feeble interactions), but it can still be produced by the plasma. This opens up a new region in the parameter space where a non-negligible abundance can be produced, which actually increases for increasing λ, contrary to what would happen in thermal freeze-out. This idea is not new, but it was recently summarised and systematised in Ref. [65], where also the term FIMP (feebly interacting massive particle) was introduced. Furthermore, a scalar that has been produced in this way had not been studied before for the case of keV sterile neutrinos production.
In our setup, the physical singlet scalar σ is produced via freeze-in from the thermal plasma. Approximately, it will have a spectrum with thermal shape, but with an overall suppression factor. This scalar σ will then fully decay into keV neutrinos N 1 via the reaction σ → N 1 N 1 . Note that, in principle, it also couples to the heavier sterile neutrinos N 2,3 . We assume this decay to be kinematically forbidden, since we consider M 2,3 ≫ m σ . Thus, the decisive decay mode is σ → N 1 N 1 . On the other hand, it might be possible to construct interesting scenarios with m σ /2 < M 2,3 < m σ , which could open up channels like σ → N 1 N 2,3 . For simplicity, we will discard these possibilities and always assume M 1 ≪ m σ ≪ M 2,3 , keeping in mind that there are several models and mechanisms which can indeed generate such a mass pattern for Majorana sterile neutrinos [19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,37,38,39,40,41].
An example evolution of the yields Y of σ and N 1 with decreasing temperature T is displayed in Fig. 1. As can be seen, we start essentially with a zero abundance of both particles (the precise value of the initial abundance plays no role as long as it is negligibly small), but with decreasing temperature the abundance of σ increases before reaching a plateau at the freeze-in temperature T ∼ m σ . However, this abundance decreases again later, due to the decays σ → N 1 N 1 . Since every scalar σ decays into exactly two N 1 's, this implies Y N 1 (late times) = 2Y σ (early times) as long as no N 1 's are produced from other sources. If the N 1 's are fully non-relativistic at late times, this also implies the relation mσ Ω σ h 2 between the final abundances, which makes it evident that this mechanism is useful to correct an overabundance of σ by a suitable mass ratio M 1 /m σ . However, since the N 1 's can also be semi-relativistic ("warm"), at least for times close to their production, the above relation could receive a correction factor in case the yield Y N 1 is not evaluated at a late enough time. Nevertheless, as an estimate, the above formula can be applied.  Figure 1: Example variation of yields Y N 1 and Y σ as a function of the temperature T , cf. Sec. 3.1 for details. As can be seen from the figure, a significant abundance of σ gradually builds up due to freeze-in, before the decays σ → N 1 N 1 set in and, at the same time, a significant amount of keV sterile neutrinos N 1 is produced.
Finally, note that the assumption that the keV neutrinos N 1 are produced exclusively by the described scalar decays does not always need to be true. In particular, in a setting where there is a non-negligible active-sterile mixing between N 1 and the light neutrinos ν i , a certain contribution to the abundance of N 1 's produced by the DW mechanism is unavoidable. We will take this contribution into account by estimating the maximal amount of keV neutrinos which can be produced by DW, without violating the X-ray bound or overproducing the DM. However, we would like to stress that our production mechanism does not need active-sterile mixing. While such a mixing may or may not be desirable from a phenomenological perspective, there are settings known in which it is exactly zero [37,62,63,64]. In such a scenario the production of keV neutrinos by a combination of the standard DW and SF mechanisms would fail, while our mechanism (as well as the version where σ does enter thermal equilibrium) could still be valid.
After having discussed the general idea behind our proposal, we will now present the more technical aspects of our work.

3
Details of our analysis 3.1 The Model The particle content of the SM is extended by three right-handed sterile neutrinos N a (a = 1, 2, 3) and one real scalar singlet S [60]. The Lagrangian is which consists of the SM, kinetic terms of the sterile neutrinos N a , Yukawa interactions f a of the singlet S with N a , and a scalar potential V scalar . Finally, L ν is the part of the Lagrangian giving mass to the light neutrinos. In the simplest setting, we would Then, a type I seesaw mechanism [68,69,70,71,72] could be at work using the right-handed Majorana masses for N a arising from a VEV f = S , at least if the Yukawa couplings respect the observational X-ray bound [73]. Alternatively, there could exist, e.g., more complicated seesaw-type mechanisms or radiative light neutrino mass generation [74,75,76,77,78]. Since we do not rely on a specific mechanism, we will leave the mass generation of light neutrinos unspecified. Any realistic setting must provide a way to generate a viable light neutrino mass and mixing pattern, but the details do not play a decisive role in our production mechanism.
We restrict our considerations to a potential V scalar which only depends on the absolute value of the SM Higgs field H and on even powers of the real scalar singlet S. Such a potential results from a global symmetry, e.g., lepton number, and does not impose a severe restriction. Assuming a global Z 4 = {±1, ±i} symmetry, 3 such that S → −S and N k → iN k (while all other fields transform trivially), the most general potential is: The SU(2) Higgs doublet H ∼ (2, +1) and the scalar singlet S (1, 0) are parametrized as and S = f +σ.
Note that the Goldstone bosons h ± are eaten by W ± to make them massive, similar to neutral boson ρ being eaten by the Z 0 . All other components are physical:h is the SMlike Higgs andσ is a physical singlet scalar. The VEVs are given by where v = 246 GeV (in our convention), and S = f . Note that f could potentially be large.
Inserting the VEVs, H † H → v 2 /2 and S 2 → f 2 , and differentiating the potential with respect to v 2 and f 2 , respectively, gives the minimum conditions The Higgs portal coupling λ results in mixing of the physical scalar fields. Concentrating on the potential terms which are proportional toσ 2 ,h 2 , andσh, and inserting the minimum conditions, Eq. (4), the mass matrix in the interaction basis (h,σ) T reads: In the basis (h, σ) T of mass eigenstates we have, in the limit of small λ, Interpreting the transition from the interaction to the mass basis as an abstract rotation, such that Eqs. (5) and (6) yield: The independent parameters are the singlet mass m σ , the Higgs portal λ, and the VEV f of the singlet. Since the Higgs portal λ is small, we can practically identify h withh and σ withσ. We will use the notation h and σ in the following. In our numerics, we have fixed the SM Higgs mass to 125 GeV in accordance with the experimental results by the ATLAS [85,86,87,88,89,90,91,92,93,94] and the CMS [95,96,97,98,99,100,101,102,103] collaborations. In addition we assume m σ > m h for definiteness and m σ < f in order to avoid being in danger to enter a non-perturbative regime, i.e., we vary m σ between the upper 1σ limit of m h < 126.4 GeV [104] and f .

Relic density
The relic density of our DM candidate particle N 1 is produced by the decays of a frozenin real scalar singlet particle σ. The Boltzmann equations for the annihilation and the decay processes are given in Eqs. (C-2) and (C-12), respectively. To calculate the relic density of N 1 , we have to solve a system of coupled equations describing simultaneously the annihilation and decay processes as it is done in, e.g, Ref. [105].
We have numerically solved the following two coupled Boltzmann equations: see Appendix C, in particular Eqs. (C-6) and (C-13), for detailed information. The equilibrium yield is given by with g σ = 1 being the spin degrees of freedom for the particle σ, x ≡ mσ T , and For the definitions of h eff , g eff and of the Bessel functions, see Appendices A and B. As already explained, the DM particle is the lightest sterile neutrino N 1 which is produced by the frozen-in real scalar singlet σ due to out-of-equilibrium decays, σ → N 1 N 1 . The thermally averaged cross section times relative velocity σ ann v for the real scalar singlet σ is calculated numerically using the micrOMEGAs package [106]. Γ(σ → N 1 N 1 ) is the thermally averaged decay rate for the decay σ → N 1 N 1 and the analytically determined decay width in the rest frame of the decaying particle σ is See Eqs. (C-9) and (C-11) for the definition of Γ(σ → N 1 N 1 ) . Finally, the relic density can be obtained using the following formula: with Y 0 = Y N 1 (T 0 ) being the yield of the DM particle at late times.

Existing constraints on the free streaming horizon
In order to determine whether the neutrinos generated act as CDM or WDM, one would actually need to determine the velocity profile and do a full simulation of the resulting structures in the Universe, see e.g. Ref. [107]. However, this is a big effort and far beyond the scope of this paper. Alternatively, one gets at least an indication by computing the so-called (co-moving) free-streaming horizon r FS , which can be interpreted as the mean distance which the DM particles would travel if they were not bound by gravitation at some point. The free-streaming horizon is defined as [7] r FS = where t in is the initial time at which the integration starts, t 0 is the current time, v(t) is the mean velocity of the DM particles, and a(t) is the scale factor. The free-streaming horizon is a co-moving quantity and, as we will see, one can define a free-streaming horizon of 0.1 Mpc [108], which is about the size of a dwarf galaxy, as the separation between HDM (λ FS > 0.1 Mpc) and WDM (λ FS < 0.1 Mpc). In turn, free-streaming horizons which are considerably smaller typically correspond to CDM. Note that this is in some sense an artificial definition, as we will explain in detail in Sec. 4, but it nevertheless gives a good orientation in practice. As we will see, the condition r FS < 0.1 Mpc will lead to a lower bound on the mass M 1 of the keV sterile neutrino. We will compute this bound by an approximate solution of the integral in Eq. (16). 4 Let us start with the integral boundaries. The production time of the DM particles can be approximated by t in ≡ t prod + τ , where t prod is the time of freeze-in (i.e., the time when the temperature equals the FIMP mass m σ [65]) and τ = 1/Γ is the lifetime of the scalar particle σ, cf. Eq. (14). The scale factor a(t) can be approximated as a(t) ∝ t 1/2 [a(t) ∝ t 2/3 ] for radiation [matter] dominance. Note that it is perfectly fine to neglect the vacuum-dominated part of the integral in Eq. (16) and to assume matter-dominance until t 0 , since very late times practically do not have any effect on the result [7]. This treatment is perfectly motivated and serves as an easy approximation. However, one still has to take into account the entropy dilution from the time of production, which happens at a very high temperature, to the current time. This amounts to a further factor of ξ −1/3 [60], with an entropy dilution factor given by where we have taken both the real scalar σ and the keV Majorana neutrino N 1 to contribute to radiation at high temperatures. 5 Since the scalar σ has been produced in a significant amount at the time of its decay, this should not be a bad approximation, and for the same reason also a significant amount of N 1 's should be around. The crucial question is how to determine the average velocity v(t) . For simplicity, we assume an instantaneous transition between the highly relativistic and the fully nonrelativistic regimes, where t nr is the time when the particle becomes non-relativistic, defined by the equality between its average momentum and its mass, p(t) = M 1 . This average momentum can be extracted from the distribution function of the DM particles. For non-relativistic parent particles σ (which is a good approximation [65]), this distribution function is given by [52,110,111,112] f (p, where β is a normalization factor that will turn out to be irrelevant for our purposes, p is the co-moving momentum, and the DM temperature is defined as ≃ mσ 2 is the DM momentum in the center-of-mass frame and the decay time t d is defined as H(t = t d ) = 1 2t in [112]. Since the particle production happens during radiation dominance, we know that H(t d ) = 1/(2t d ) and can thus identify t d ≡ t in .
Defining "early" and "late" production of the DM particles as t in < t eq and t in > t eq , respectively, where t eq is the time of matter-radiation equality, one can easily compute the free-streaming horizon by splitting the integral from Eq. (16) into different pieces for radiation/matter dominance and highly relativistic/non-relativistic DM particles. The result is see Appendix D for details and in particular Eqs. (D-4) and (D-5). Note that the two parts of Eq. (20) coincide for t nr → t eq . Furthermore, as to be expected, r FS always increases with increasing t in . We have used Eq. (20) to mark the excluded region of HDM (r FS > 0.1 Mpc) later on in Figs. 3 and 4. We will furthermore indicate the CDM regions (r FS < 0.01 Mpc), which are not excluded and instead reveal that also a keV-mass particle can act as CDM, depending on the details of its production.

Collider bounds on the production of Dark Matter
In colliders, a DM signal can be detected through monojet or cascade events. If the DM particle is stable, it does not decay inside the detector volume and thus leaves its track as missing energy, which can be reconstructed. Comparing simulated DM events with data analysis allows to constrain the DM interaction and its mass. Bounds exist on DM masses around 1 GeV. Specific collider constraints on Majorana fermion DM can be found in [113]; for Dirac fermion, complex scalar, and real scalar DM, see [114]. Note that these constraints are not relevant in our case, since the DM is a sterile neutrino N 1 with a mass in the keV range, produced by the decay of a frozen-in scalar singlet FIMP. The scalar singlet FIMP itself is produced via the Higgs portal. In Ref. [115], the allowed region for the Higgs portal coupling λ and the mass of the scalar singlet is presented. From that reference it follows that the strongest upper bound on the Higgs portal is λ < 0.01, but for scalar singlet masses m σ 60 GeV there is essentially no constraint. In the mechanism we propose, the Higgs portal coupling is of order λ ∼ O(10 −8 ), i.e., given the mass range of our scalar singlet and its feeble interactions the constraints in Ref. [115] for the allowed λ − m S region are not relevant for us. For completeness, see also [116] for LHC sensitivities on the Higgs portal.
In addition, Ref. [115] presents constraints on the invisible decay width of the SM Higgs doublet H; for m h = 125 GeV, the invisible decay width has to be smaller than approximately 0.0025 GeV. In our model, the SM Higgs doublet H decays invisibly into the DM particle N 1 and an active neutrino at tree-level, if light neutrino masses are generated by type I seesaw or by any other mechanism which allows for a neutrino Yukawa coupling given by L ν = −y α1 D L αH N 1 . Since H is of the order O(100) GeV, M 1 of order O(1 − 100) keV, and the light neutrino masses in the sub-eV range, the Yukawa couplings y α1 D must be tiny such that the decay width of H → N 1 ν α is much smaller than 0.0025 GeV.
To conclude, all existing collider bounds on production of DM are not relevant for our mechanism and do not constrain the parameters we are considering in our numerical analysis.

Results
We have numerically solved Eqs. (9) and (10) in order to determine the final abundance of keV sterile neutrinos N 1 . First of all we scanned over a range of values for the Higgs portal coupling λ in order to identify the successful region to obtain the correct relic abundance. The only requirement we impose on λ is that λ 10 −6 in order not to enter thermal equilibrium [60]. The result of this scan can be found in Fig. 2, where we plot the abundance regions for different values of the coupling, λ = 10 −7,8,9 , as a function of the keV neutrino mass M 1 . The broadening of the corresponding bands originates from the variation over the scalar mass m σ . For definiteness, we assume that the singlet scalar mass is always larger than the SM-like Higgs mass m h ≈ 125 GeV (corresponding to the upper end of the bands in the plot). Furthermore, in order to avoid entering a potentially non-perturbative regime, we also assume that m σ < f (corresponding to the lower ends of the bands in the plot), where f is the VEV of the singlet field S. In Fig. 2, we present the plots for the two example values f = 500 GeV and f = 1 TeV, which are perfectly compatible with all bounds. As can be seen from Fig. 2, the successful value of the Higgs portal coupling λ should be around 10 −8 , more or less independently of the value of the VEV f . Accordingly, we will focus on the region λ ≈ 10 −8 in what follows and investigate this region in greater detail in what concerns the relic abundance and in particular the experimental and observational bounds.  A more detailed investigation of the successful regions in parameter space can be found in Figs. 3 and 4, where we have indicated the region of the correct abundance (i.e., within the 3σ ranges of the Planck data [4]), as generated by scalar FIMP production only, by the orange band in the plot. The parameter values for the plots are chosen as f ∈ {500 GeV, 1 TeV}, with λ ∈ {1.0·10 −8 , 1.2·10 −8 } for Fig. 3 and {1.5·10 −8 , 2.0·10 −8 } for Fig. 4. As can be seen from the plots, the iso-abundance lines reveal a more or less linear dependence of the keV sterile neutrinos mass M 1 on the scalar singlet mass m σ . This feature can be understood easily by observing that the final DM energy density must be equal to the initial energy density in scalar σ particles, which can at most be redshifted. Since this initial energy density is non-relativistic, it can be written as ρ σ = m σ n σ , where n σ is the number density of σ-particles. Similarly, the energy density in N 1 can be computed by the non-relativistic expression for late times, cf. discussion in Sec. 2, since in the successful regions in the parameter space the DM particles become non-relativistic within the age of the Universe.
In addition, we have indicated some important bounds. As explained, we have assumed that m σ > m h , and by gray rectangles we indicate the corresponding regions left of the upper 1σ bound on m h of 126.4 GeV [104]. Furthermore, we know that HDM is excluded or, rather, bound to make up at most 1% of the DM in the Universe [14,15] by considerations of cosmological structure formation. A rough way to quantify when DM particles are HDM, WDM, or CDM is the co-moving free-streaming horizon r FS , cf. Sec. 3.3. Since it is a bit crude to attribute the property of a whole velocity spectrum of DM particles to one single number, it is to some extent a question of definition where to draw the lines between the three DM categories. A relatively common choice, which somewhat representatively reflects the use of the three terms in the literature [117,118] is to take the border between HDM and WDM at a free-streaming horizon of roughly r FS = 0.1 Mpc, where larger values signal HDM which is forbidden. Note that this value is physically motivated due to the size of dwarf satellite galaxies being in that range. However, between CDM and WDM there is not a very well-defined boundary, since it is not easy to unambiguously define at which value of r FS the structure formation on small scales starts to depart from the pure CDM case [7,119,120]. However, it is clear that the free-streaming horizon for CDM should be "significantly smaller" than the one for WDM. For definiteness, we have therefore decided to simply take a value that is by one order of magnitude smaller than the one for the HDM-WDM boundary. Keeping in mind that this distinction between WDM and CDM is also a matter of definition, the values of r FS which we used are: Thus, in the plots displayed in Figs. 3 and 4, the thick red line at the bottom of the plots marks the transition between WDM and HDM, and the light red region below this line is excluded by structure formation. The light blue region in the upper part of the plots, bounded by the thick blue line, corresponds to CDM and the white region marks the WDM sector. Here, it is worth to point out that in a considerable region of the parameter space, keV sterile neutrinos (with large enough masses) can be cold DM (or, more precisely, indistinguishable from CDM according to our definition), in contrast to most of the scenarios for keV sterile neutrino DM. To some extent, this is a simple reflection of the fact that our DM production happens in the early Universe, but the more crucial point is that the mass ratio m σ /M 1 happens to be in the correct region to allow for a sufficient cooling time. Figs. 3 and 4 reveal that the region of the correct DM relic abundance lies in the cold or warm DM parameter space, depending on the specific value of λ, respectively.
We also have the possibility to produce part of the DM in keV sterile neutrinos by the ordinary DW-mechanism [45], in addition to the production by the mechanism proposed here. This contribution depends on the active-sterile mixing angle θ 1 of the keV sterile neutrino N 1 , and it can be estimated by the approximate formula [121]: Note that, if the keV sterile neutrino makes up all the DM in the Universe and if it is unstable under N 1 → νγ, then there is a strong bound from the non-observation of the corresponding X-ray line (see Refs. [17,18,61] for recent collections of bounds). In the plots, we have represented the corresponding maximal (i.e., for the largest allowed value of sin 2 θ 1 ) addition of particle production due to the DW mechanism by the purple bands. As can be seen from the plots, this would shift the allowed regions (i.e., the regions where the total abundance of keV neutrinos, as produced by both mechanisms together, is within the 3σ regions of Planck) towards slightly larger values of m σ . For very low M 1 , there is a considerable DW-production resulting from the comparatively weak X-ray bound in this mass region. In this region, nearly all the DM can be produced by the DW-mechanism, which for these masses completely dominates the production by frozenin scalars, if the maximally possible value is taken for the active-sterile mixing. However, from studies of the Lyman-α forest, the corresponding lower bound on the keV sterile neutrino mass, when the DW-mechanism is at work, is between 8 and 10 keV [7] (note that this accidentally coincides with the light red HDM region in our plots). Thus, this region of the parameter space is excluded. On the other hand, depending on the exact value of the active-sterile mixing, the combined abundance of keV sterile neutrinos produced by both mechanisms could also lay in between the orange and purple bands, that we indicate in the plots. We want to stress that the orange bands correspond to production by scalar FIMPs only, i.e., this is the region of correct abundance for a vanishing active-sterile mixing, θ 1 ≡ 0. While this may not be desirable from a phenomenological point of view (e.g. for a possible detection of the X-ray line [122,123,124] or for a potential detection of modifications of neutrino-less double beta decay [61,125]), vanishing active-sterile mixing may be very natural in certain settings [37,62,63,64]. In such frameworks it would be impossible to produce keV sterile neutrinos via the DW and/or SF mechanisms, but our mechanism (as well as the version in which the scalar freezes out) could be easily used as an alternative. Finally, we have also indicated possible mass limits arising from the X-ray bound. For example, if the bound on the active-sterile mixing is taken to be sin 2 (2θ 1 ) < 10 −13 , this excludes keV sterile neutrino masses M 1 above 64.5 keV, while the values below are consistent with the X-ray bound (but not necessarily with the HDM bound). If activesterile mixing is not present, then there is no fixed upper bound on M 1 , and alternative scenarios with stable, e.g., MeV or GeV sterile neutrino DM could be found, too.
The general message of our plots is that there is considerable room for keV sterile neutrinos to be produced by scalar FIMPs and to be compatible with all bounds. Accordingly, if in a certain setting the Higgs portal coupling of a singlet scalar is bound to be very small, it could still be used to produce sterile neutrino DM.

Conclusions
In this paper, we have presented a new and successful mechanism for the production of keV sterile neutrino DM. The mechanism is based on the so-called freeze-in of a scalar particle σ, which has too feeble interactions with SM particles (and hence with the primordial thermal plasma) to ever enter thermal equilibrium, but whose mixing with the SM-like Higgs boson is nevertheless large enough to gradually produce a significant abundance of σ's in the early Universe. If these σ's are unstable, they can decay efficiently into pairs of keV sterile neutrinos N 1 , thereby generating the required DM abundance. We have numerically solved the corresponding system of coupled Boltzmann equations and we have also presented a discussion of all potentially relevant bounds.
Depending on the exact value of the first generation active-sterile mixing angle, the DM abundance generated by the mechanism proposed here must be corrected by a contribution from the (generic) Dodelson-Widrow mechanism. We have estimated the maximal effect of this additional amount of keV neutrinos, which alters the successful regions in the parameter space, without however spoiling the proposed production mechanism. On the other hand, it is worth to note that our mechanism does not at all rely on active-sterile mixing, and it could very well live even with a vanishing active-sterile mixing angle. This point could be particularly interesting for models which avoid the strong X-ray observational bound on the active-sterile mixing, by stabilising the keV neutrinos and at the same time forbidding any mixture of active and sterile states.
While similar mechanisms had been proposed previously for early and late freeze-out (and subsequent decay) of the scalar, our proposal opens up a new window in a region of the parameter space where freeze-out is not at all possible. This is particularly interesting for models which predict a very small Higgs portal coupling between the singlet scalar field σ and the SM-like Higgs. Apart from being applicable to many settings where a suitable scalar is available, the main advantage of our production mechanism is that it happens at relatively early times, thereby causing the DM particles to become non-relativistic already at high temperatures (which they do not feel due to their feeble interactions). Hence, depending on the exact values of the parameters, the keV neutrinos can be cold DM in a significant fraction of the parameter space. This is particularly interesting in case the X-ray bound can be circumvented in a concrete model, in which case a significant region of the M 1 -m σ parameter plane can lead to the correct relic abundance of DM and be consistent with all bounds.
A Appendix: Effective degrees of freedom We followed the notation of Ref. [126], where the energy density ρ i and the entropy density s i for a particle species i are defined as: with the temperature T i of the particle species i. The effective degrees of freedom g i eff and h i eff for energy and entropy density, respectively, are defined as with x i ≡ m i /T i , η i = 1 for Fermi-Dirac, η i = −1 for Bose-Einstein and η i = 0 for Maxwell-Boltzmann. The number of internal degrees of freedom is denoted by g i . In our numerics, we accounted for the contribution of the real scalar singlet σ and the sterile neutrino N 1 to the total energy and entropy effective degrees of freedom given as

B Appendix: Modified Bessel functions
The modified Bessel functions K n (x) of the second kind obey the identity For a Maxwell Boltzmann distribution with zero chemical potential, the equilibrium number density n eq of a particle with g internal degrees of freedom and mass m is where x = m/T . For the yield Y = n s with entropy density s = 2π 2 45 h eff T 3 follows:

C Appendix: Annihilation and decay reactions
In the usual Friedman-Robertson-Walker metric, the Boltzmann equation for the number density n of a particle species can be written as d dt n + 3Hn = C[n] , (C-1) where C is the collision operator expressing the number of particles per phase space volume that are lost or gained per unit time due to interactions with other particles. For the standard annihilaton process σ σ ⇋ SM SM of a real scalar singlets σ into Standard Model particles, the Boltzmann equation for the number density n σ reads d dt n σ + 3Hn σ = − σ ann v (n 2 σ − n 2 σ,eq ) ≃ σ ann v n 2 σ,eq , (C-2) where the latter approximation is valid for the freeze-in case, for which the initial number density and thus the initial abundance can be neglected [65]. Furthermore, σ ann v is the relativistic thermally averaged annihilation cross section and v is the Møller velocity. Following the discussion of Ref. [126], it is possible to write We have generated the correct Feynman rules using LanHEP [127] and we have used micrOMEGAs [106] for the calculation of Eq. (C-3).
In the radiation dominated era, the Hubble expansion rate can be expressed as where G N is Newton's gravitational constant. Furthermore, in the radiation dominated era, the expansion age t of the Universe with Ω tot = 1 equals: In terms of the abundance Y = n s with the entropy density s = 2π 2 45 h eff T 3 , it follows: The superscript A serves as indication of the annihilation process. The decay processes σ → N 1 N 1 of a real scalar singlet σ into two sterile neutrinos N 1 is described by the following phase space integration: Neglecting, the Pauli blocking and enhancing factors, we can define with Γ * (σ → N 1 N 1 ) the decay width for the particle at energy E σ . The above phase space integration yields: with Γ(σ → N 1 N 1 ) the decay width in the rest frame of the decaying particle σ, i.e., Γ(σ → N 1 N 1 ) = Eσ mσ Γ * (σ → N 1 N 1 ). Thus, for the decay process σ → N 1 N 1 of a real scalar singlet σ into two sterile neutrinos N 1 , the Boltzmann equation for the number density n N 1 reads as (again we define x = m σ /T ) d dt n N 1 + 3Hn N 1 = 2 K 1 (x) K 2 (x) Γ(σ → N 1 N 1 ) n σ . (C-12) The factor 2 accounts for the fact that two sterile neutrinos N 1 are produced per decay. In terms of the abundance Y = n s with the entropy density s = 2π 2 45 h eff T 3 , it follows: where the superscript D serves as indication of the decay process.

D Appendix: Free-streaming horizon
With the distribution function f (p, t) of the DM particle, given as (D-1)