keV Sterile Neutrino Dark Matter from Singlet Scalar Decays: The Most General Case

We investigate the early Universe production of sterile neutrino Dark Matter by the decays of singlet scalars. All previous studies applied simplifying assumptions and/or studied the process only on the level of number densities, which makes it impossible to give statements about cosmic structure formation. We overcome these issues by dropping all simplifying assumptions (except for one we showed earlier to work perfectly) and by computing the full course of Dark Matter production on the level of non-thermal momentum distribution functions. We are thus in the position to study all aspects of the resulting settings and apply all relevant bounds in a reliable manner. We have a particular focus on how to incorporate bounds from structure formation on the level of the linear power spectrum, since the simplistic estimate using the free-streaming horizon clearly fails for highly non-thermal distributions. Our work comprises the most detailed and comprehensive study of sterile neutrino Dark Matter production by scalar decays presented so far.


Introduction
Invoking something invisible that even our most sensitive detectors cannot detect does not sound like science. Yet, in contemporary cosmology, we have so much indirect evidence for Dark Matter (DM) that hardly any scientist doubts its existence. DM is a non-luminous form of matter, i.e., it does not interact with light -unlike everyday objects around us. Nevertheless, having entered an era of precision cosmology, we have been able to determine that DM outweighs ordinary matter by a factor of about five in the energy balance of the Universe [1]. Further observational evidence such as galaxy [2] or cluster dynamics [3] and the Bullet Cluster [4] allow us to constrain the properties of DM: we are sure that it is a form of matter (and not, e.g., modified gravity [5]), and our best guess for its identity is a new elementary particle that is electrically neutral and massive [6]. Importantly, it must have been produced in the early Universe in the right amounts and with a momentum spectrum suitable not to spoil the formation of cosmic structures.
It is generally accepted that DM was responsible for the emergence of structures in the Universe [7], as seen in elaborate N -body simulations on high performance computers [8]. Historically, the most generic DM candidates were WIMPs (Weakly Interacting Massive Particles), which appear in popular scenarios like supersymmetry. Such particles typically form cold DM (CDM), i.e., particles produced by thermal freeze-out [9,10] which have non-relativistic velocities. The velocity of the particles strongly impacts structure formation. For example, hot DM (HDM), which is highly relativistic, is excluded as it would have wiped out all small structures in the Universe [11,12]. But CDM may have problems, too: it possibly forms "too many" small halos, which might even have the "wrong" structure (known small scale issues include the missing satellite problem [13,14], the abundance of isolated small halos [15], the too-big-to-fail problem [16,17], and the cusp-core problem [18,19]). While these may be cured once baryons are correctly included in the simulations (see, e.g., [20][21][22]), an alternative attempt is to consider non-cold DM settings. In the literature, these ideas have triggered the slightly unfortunate terminology of warm DM (WDM), i.e., DM particles with a thermal spectrum (i.e., Bose-Einstein or Fermi-Dirac) of a temperature close to their mass. This is also assumed in many astrophysical studies, from Lyman-α bounds [23][24][25][26] over galaxy formation [27][28][29] to N -body simulations [30][31][32]. However, looking at realistic settings, many DM-models feature a non-thermal spectrum, which one cannot associate any temperature with.
Among the most popular non-cold DM candidates is a sterile neutrino with a mass of a few keV, see Ref. [33] for a recent collection of information on this topic. Sterile neutrinos are motivated from a particle theory point of view, as they are related to light (active) neutrinos and possibly even involved in their mass generation, see e.g. Refs. [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51] for concrete models. Furthermore, the topic has gained some attention in recent years, because a sterile neutrino N would -with a very small rate -decay like N → νγ, thereby producing a nearly monoenergetic X-ray photon. A detection of the corresponding line signal has been claimed by two groups in 2014 [52,53], but it was very actively disputed -see Ref. [33] for a detailed discussion and Ref. [54] for the (still ambiguous) data that the Hitomi satellite could take before enduring its unfortunate fate. It may also be worth mentioning that there could be a "dip" in the cluster data [55], tending to make a nonobservation of the line in stacked data more consistent.
Sterile neutrinos with a sufficiently large lifetime are excellent DM candidates provided that 1.) an efficient production mechanism exists in the early Universe which 2.) produces a spectrum in accordance with cosmic structure formation. The most generic idea is to produce sterile neutrinos by their small admixtures to active neutrinos, by non-resonant active-sterile transitions, a mechanism first proposed by Langacker [56] and related to DM by Dodelson and Widrow (DW) [57]. In modern terms, one would refer to these sterile neutrinos as FIMPs (Feebly Interacting Massive Particles [58]), which do not thermalise but are instead produced gradually in the early Universe via freeze-in by their feeble couplings to Standard Model (SM) particles. While this mechanism is rather simple, it is excluded by data, due to the particles being too close to the HDM limit [26,59]; still, it is an unavoidable addendum to the spectrum produced by any other mechanism, which can modify the DM distribution function at least for sterile neutrino masses below 3 keV [60].
A popular alternative is based on a resonant enhancement of the active-sterile neutrino transitions by the presence of a primordial lepton number asymmetry. First proposed by Enqvist and collaborators [61] and related to DM by Shi and Fuller (SF) [62], this mechanism has been studied actively [63][64][65][66][67][68][69], and it does indeed yield a spectrum colder than that produced by DW. However, while there is no perfect agreement between the results of different groups, at least for the results that are publicly available [63,67,69] it seems unclear whether they are in full agreement with cosmic structure formation [70][71][72].
On the other hand, one could produce sterile neutrinos thermally via freeze-out, if they had non-trivial charges beyond the SM, as long as the resulting overabundance is diluted by a sufficiently efficient production of additional entropy [73][74][75]; this, however, is constrained by big bang nucleosynthesis [76].
A completely different direction is to produce sterile neutrinos by the decays of other particles. A generic possibility is a decaying singlet scalar S, which can easily couple to sterile neutrinos after having been produced itself in the first place. Apart from this scalar being the inflaton [77][78][79], in which case it had been present all along in the early Universe, one can tune the Higgs portal coupling of S in such a way that it is either produced like a WIMP via freeze-out [80][81][82] or like a FIMP via freeze-in [83,84]. Variants have been presented in [85][86][87][88][89][90][91][92][93][94], some also featuring other parent particles such as vectors [95][96][97], Dirac fermions [98], or pions [99,100]; related aspects such as influences from inflation [101] or thermal corrections [102] have been discussed, too.
Our main goal is to close a gap in the treatment of sterile neutrino production from scalar decays. In earlier works [80,81,83,84], simplifying assumptions have been applied to at all arrive at a result, such as neglecting active-sterile mixing, taking the number of relativistic degrees of freedom to be constant, and assuming a heavy scalar. This is also true for a previous paper by two of us (AM & MT) [84], which was the first to show how to numerically compute momentum distribution functions of sterile neutrinos from scalar decay. While it has been proven that neglecting active-sterile mixing is in fact a very good assumption [60], going to the low-mass region of the singlet scalar -where the number of degrees of freedom is not constant -poses new technical challenges. The only treatment available for this regime was put forward in [88]. However, the authors only used rate equations, such that an actual computation of the momentum distribution function and the confrontation with structure formation data cannot be reliably performed. We will in this work present the full numerical computation on the level of momentum distribution functions instead, including a detailed treatment of all subtleties related to the many new technical aspects that arise in this regime. We will give an a-posteriori justification of some of the assumptions made in [88] (such as the Higgs staying equilibrated during S-production), while we will also reveal that others (like the Higgs degrees of freedom in the unbroken phase) are maybe less justified. We will furthermore show in great detail how to constrain the resulting spectra by cosmic structure formation.
Our study presents the first complete treatment of scalar decay production in the whole parameter space. It involves a fully general solution of the production equations on the level of distribution functions. Apart from our techniques being transferable to virtually any type of decay production, our study will guide the particle physics community towards using limits from cosmic structure formation without having to leave their comfort zone completely. We thus contribute to closing the gap between particle physics models, early DM-production, and their phenomenological consequences.
This text is structured as follows. We first present a qualitative discussion of decay production of sterile neutrinos in Sec. 2, before explaining how to solve the evolution equations and to apply the relevant bounds in Sec. 3. Our main results are presented in Sec. 4, which features a thorough discussion of all relevant aspects from the production process over the distribution functions to structure formation. We conclude in Sec. 5. Technical aspects on the Boltzmann equation, on the evolution of the Higgs distribution, on the failure of the free-streaming horizon as a reliable estimator, and on the robustness of our half-mode analysis are discussed in Appendices A, B, C, and D, respectively.

The basic idea: Qualitative discussion
This section is intended to introduce the particle physics model used in our analysis and to give a qualitative understanding of its different regimes. Furthermore, we will justify the assumptions that go into the numerical computations.
Our setting introduces one real scalar singlet S (with mass m S ) and one right-handed neutrino N (with mass m N ) beyond the particle content of the SM. 1 The right-handed neutrino is coupled to the singlet scalar via a Yukawa-type interaction with coupling y, while the new scalar singlet is coupled to the Higgs doublet Φ via the most generic potential symmetric under a global Z 4 -symmetry: 2 Accordingly, the complete Lagrangian of the model reads where L ν is the part of the Lagrangian that can give mass to the active neutrinos. We do not assume any vacuum expectation value (VEV) for S in our analysis, however, we will include the constraints which would arise of the VEV of S gave the sterile neutrinos their mass, so that the case of S = 0 can be qualitatively recovered from our analysis. Note that we assume active-sterile mixing to be so small that the contribution from the DW mechanism to the production of sterile neutrinos is negligible compared to the part produced by scalar decay, based on taking into account X-ray limits on the decay of sterile neutrino DM. In Ref. [60] it has been shown that the DW modification to any previously produced population of sterile neutrinos (from whichever main production mechanism) 1 In general, any number of right-handed neutrinos can be assumed in this model. In the case of more than one generation, the scalar will then decay into all kinematically accessible right-handed states N i with branching fractions given by y 2 i / k y 2 k . If the mixing between the different right-handed states is large enough, all right-handed neutrinos will decay into the lightest state (N 1 by convention) quickly and the results of this paper stay unaltered under the substitution y 2 → k y 2 k . If, however, the mixing inside the sterile sector is small, there might be additional complications due to late injection of highly energetic DM particles by the decay of the heavier states into the lighter ones.
2 Suitable charge assignments are given by S → −S and N → ±iN . For further comments on this assumption, see [83,84] and references therein. The (rather mild) consequences of giving up this simplification are discussed in [81]. regime production channels is of a few percent at most, and completely negligible for masses m N larger than about 4 keV. Our setup mainly allows us to produce scalars from their coupling to the SM with interactions of the type X SM X SM ↔ SS. Subsequently, the scalars decay into sterile neutrinos via S → N N . We have to distinguish three different regimes: I Production before the electroweak phase transition (EWPT): all four degrees of freedom of the SU (2) L -doublet Higgs Φ contribute equally to the production/depletion of scalars from/into the thermal bath.
II Production after EWPT with m S > m h /2, where m h is the mass of the Higgs after electroweak symmetry breaking. Now the Higgs and the massive gauge bosons interact with the scalar different ways.
III Production after EWPT with m S < m h /2. This is similar to case II, the difference being that the Higgses present in the thermal plasma are now kinematically allowed to decay into pairs of scalars.
The channels corresponding to regimes I-III are listed in Tab. 1. Note that we have omitted the diagrams with a scalar S in the tor the u-channels, since they will contribute to the production only at order O (λ 3 ), i.e., suppressed by another small factor of λ compared to the leading order. In regimes II and III, the fermionic initial state f only receives a contribution from the top quark in practice, since all other channels are suppressed by their small Yukawa couplings. Couplings to lighter fermions may appear relevant at first sight, since in that case the Higgs could potentially be on-shell, which may enhance the cross section by orders of magnitude. However, this case of an on-shell Higgs needs to be subtracted adequately to avoid double-counting of decay events of thermal Higgses [103], and they ultimately do not contribute much. Note that, at precision level, even thermal corrections to the parameters can play a role [102].
In general, production will always start in regime I and then, unless finished in that regime, proceed either in regime II or in regime III after EWPT, depending on the value of m S . In the case of small Higgs portal couplings λ, the scalar will freeze in. This production mechanism is most efficient when T ∼ m S . In the case of larger λ, the scalar first equilibrates and then freezes out. In this scenario, the freeze-out time will also be linked to m S . Accordingly, there is a finite span in cosmic time (or in the temperature T , equivalently) in which the production is effectively taking place.
This time span is illustrated in Fig. 1 for four different masses m S . Red arrows correspond to small Higgs portal couplings λ, implying that the scalar never equilibrates, but freezes in instead. In this case, we have defined the time span of production starting when 10% of the final yield Y of a would-be stable scalar is produced and ending when 90% are produced. Blue arrows indicate the time spans for scalars freezing out before decaying into sterile neutrinos. In this case, there is no physically preferred initial time. Note that, even when assuming a vanishing initial abundance, thermalisation is fast compared to the time scale of the decay. The late-end boundary (i.e., at low temperatures) of the time span is, however, defined by Y /Y eq = 10, where Y eq denotes the equilibrium yield.
One can clearly see that, in cases where m S < m h /2, the freeze-in of scalars occurs only after EWPT. This can be explained by the fact that the decay of thermal Higgses after EWPT completely dominates the production as soon as this channel is kinematically accessible [88]. For m S = 65 GeV, production sets in before EWPT and stops at temperatures of some tens of GeV. In the case of a much heavier scalar (say, m S = 500 GeV), production even ends before EWPT, since the kinetic energies of the Higgses after EWPT are insufficient to produce the heavy scalar.
Note that we need the particle distribution function of all sourcing species in the plasma, like gauge bosons or Higgses, to account for the production of scalars induced by the channels shown in Tab. 1. In general, the dynamics of each of these species is determined by their own Boltzmann equation. Fortunately, it is not necessary to incorporate more than two equations into the system, since we can safely assume these particles to be in thermal equilibrium, an assumption that has been readily made in [88] even for small T , but not proven to work. We have explicitly checked that the Higgs Figure 1: Temperature ranges of the plasma temperature T ≡ T γ at which the scalar singlet S is efficiently produced in the freeze-in scenario (red) and in the freeze-out scenario (light blue). We present four different scalar masses (30,60,65, and 500 GeV) that span the variety of possible production times. T EWPT indicates the temperature when electroweak phase transition happens, while T PT is supposed to give a rough idea of the temperature when the Higgs and gauge bosons become strongly Boltzmann-suppressed, even if still equilibrated with the remaining SM degrees of freedom.
indeed stays in equilibrium at least until T ≈ 1 GeV. Note that the usual estimate for WIMP-like particles of mass m is that freeze-out occurs at T freeze-out ∼ m/20. The reason that the Higgs stays in equilibrium longer are inverse decay processes, and similar arguments apply to the gauge bosons and to the top-quark, cf. App. B. Still, the number densities of all SM particles become strongly Boltzmann-suppressed for temperatures of, say, 10-20 GeV. Hence, production always ceases at these temperatures, simply because there are hardly any particle left in the bath, even if they are still in equilibrium.

The Boltzmann equation and its solution
In order to compute the distribution functions of the scalar and the sterile neutrino, we employ the Boltzmann equation,L[f ] = C[f ], where f is the distribution function we want to determine,L = ∂ ∂t − Hp ∂ ∂p (with the Hubble function H) is the Liouville operator in a homogeneous and isotropic Friedman-Robertson-Walker universe, and C contains the collision terms describing interactions of the particles. 3 In our case, we actually need to solve two coupled Boltzmann equations, for the two distribution functions f N and f S of the sterile neutrino N and the scalar S, respectively: The collision terms are themselves sums of several individual terms which encode the details of the respective processes as well as information on whether they contribute to production or depletion of the species under consideration. Depending on the regime, various diagrams listed in Tab. 1 can contribute to the production of scalars and hence, ultimately sterile neutrinos. Explicitly, the scalar collision term C S consists of the following: • Regime I (T > T EWPT ): • Regime II (T < T EWPT and m S > m h /2): • Regime III (T < T EWPT and m S < m h /2): The labels should be quite intuitive. For example, in regime I, the term C S φφ↔SS encodes both the annihilation of two Higgs doublets Φ into two scalars S and vice versa, with opposite signs for the two cases. On the other hand, C S S→N N will always come with a negative sign, since it only describes the decay of the scalar into sterile neutrinos. Similarly, the only difference between regimes II and III is the appearance of the Higgs decays into two scalars, C S h↔SS , for sufficiently small scalar masses.
The sterile neutrino collision term is comparatively simple, and it is given by: but in this case with a positive sign as it creates sterile neutrinos. Here we neglected the inverse decay, i.e. the process N N → S, because the sterile neutrino is a FIMP, such that its abundance is always far below its would-be equilibrium abundance for large enough temperatures T . Hence the process N N → S is suppressed by the small number of sterile neutrinos present in the Universe. Note that the collision terms will in general depend on the distribution functions f i of all particle species i. However, as argued, we only need to compute the distributions f S and f N , since all SM species can be assumed to be in thermal equilibrium such that their distributions are well approximated by simple Boltzmann factors. 4 To anticipate one particular example collision term, we quote from Eq. (A-23) the one describing the production of sterile neutrinos of momentum p at temperature T by the decaying scalars, where the lower boundary of the integral over scalar momenta p is given by . This term is indeed positive, since it produces sterile neutrinos. As expected, this equation contains the distribution function f S of the scalar at temperature T , making it explicit that Eqs. (4) and (5) are coupled and have to be solved as a system of equations. For a detailed discussion of all collision terms we refer the reader to App. A.1.

How to crack a coupled system of non-linear partial integro-differential equations in two variables
Glancing at the explicit forms of the collision terms in App. A.1, we can see that the task of computing sterile neutrino production from scalar decay is a highly non-trivial one: it necessitates the solution of a coupled system of partial integro-differential equations in two variables. In an abstract, yet intuitive manner, the system of equations we have to solve turns out to have the following form: 5 where K S/N and D are known functions directly related to the collision terms C, and f eq i denotes the (possibly hypothetical, in case the interactions are too feeble) equilibrium distribution of particle i. The question is now how to tackle such a problem. Before starting the computation we note that, while the sterile neutrino distribution function f N does depend on the distribution function f S of the scalar, the converse is not true. The reason is that, as we had argued, we can neglect any inverse processes involving sterile neutrinos. This already yields to a considerable simplification, allowing us to solve the first Eq. (11) independently and to then insert the result obtained into the second equation. However, even then, we are still left with a non-linear partial integro-differential equation in f S which is highly non-trivial to solve numerically (and impossible to solve analytically). We will thus need to play some tricks, in order to tame this equation.
The first step is to perform a transformation of variables, (t, p) → (r, ξ), such that the left-hand side of the first Eq. (11) only involves a derivative with respect to a single variable. As shown in App. A.2, this is possible for r = f (t) being an arbitrary function f independent of the scalar momentum p, and ξ(p, t) = g a(t) a(t 0 ) p simultaneously being an arbitrary function g of a specific combination of time t and momentum p, which involves the scale factor a(t) and an arbitrary reference time t 0 . Thus, in fact, there exists a whole family of functions f and g which would simplify our complicated equation, however, choosing them in a smart manner may even lead to further simplifications. In our view, the following choice of variables is particularly convenient: 6 where g s (T ) is the number of effective entropy degrees of freedom (for which we have used the numerical fit developed in Ref. [105]). The arbitrary reference mass m 0 and temperature T 0 have been chosen by us to both equal the Higgs mass, m 0 = T 0 = m h , which will later prevent us from working with extremely small or large numbers in our numerical computation. This choice of variables indeed implies certain simplifications, e.g., the Liouville operator now reads: where denotes a derivative with respect to the temperature T . A valid interpretation of the new variables is to view r as "time" variable and ξ as "momentum" variable. Indeed, r increases as the temperature decreases, just as the cosmic time t, so it can really be thought of as a rescaled or stretched time. In turn, ξ can be thought of as something like a comoving rescaled momentum. Alternatively, one could think of it as a rescaled momentum that is red-or blue-shifted with respect to the reference temperature T 0 . Indeed, the big advantage of the variable ξ is that the distributions remain constant in time r as soon as the production of particles is finished. For example, if the scalar S was stable, its distribution for late time would be given by f S (ξ, r) ≡ f S (ξ, r prod ), ∀r ≥ r prod , where r prod marks the time when no S-particles are produced anymore (e.g. the freeze-out time in case S was a stable WIMP). This effect translates into the sterile neutrino distribution, i.e., f N will also be constant in r once the scalar decays cease to be efficient. All redshifts of the momenta are automatically included in the definition of ξ.
The resulting form of the Boltzmann equations used in our implementation is where i = S, N . Indeed, we have by our substitution transformed the partial differential equation in two variables into an equation containing a differential with respect to one variable only. Of course, the collision terms in their general form of Eq. (11) have to be transformed according to this change of variables. Note that this general form can be written down for both the scalar and the sterile neutrino, but it is only necessary for the former, as the sterile neutrino evolution equation is comparatively simple once f S is known. We will therefore concentrate on the case i = S in the following. Still, the main issue remains: the collision term is nonlinear in f S , and it furthermore includes an integral over ξ, which renders this equation to be of a partial integrodifferential type, which are very difficult to solve. While several strategies for certain types of integro-differential equations (such as the Volterra-type) are available in the literature, the equation under consideration turns out not to be of any such type, so that we have to find a dedicated workaround. We apply a trick to get at least an (arbitrarily good) approximation by discretising the equations and substituting the integral over ξ by a finite sum over M discrete momentum values ξ i , i ∈ {1, ..., M }. This transforms the original non-linear integro-differential equation (14) into a system of coupled non-linear ordinary differential equations for the different modes f i S (r): Note that the integral results in a "non-local" coupling, since it couples any mode f i S not only to some "neighbouring" modes, but to all modes instead.
Finally, we have to numerically solve the resulting system of differential equations. This step involves one more difficulty, since the most simple strategies (such as the Runge-Kutta method) fail for large parts of the parameter space due to the stiffness of the system. Another more technical issue is that packages such as Mathematica natively employ list-based algorithms, while matrix-based ones are much more appropriate for the equations under consideration. We have thus found it more convenient to use the built-in Matlab solver ode15s [106,107], which is particularly efficient when dealing with stiff problems. Nonetheless, one should try to tweak any solver algorithm to take advantage of our a-priori knowledge that a distribution function must never attain negative values. This way, we have been able to solve the equations for the scalar distribution functions f S efficiently. The resulting numerical functions can then be inserted into Eq. (14) for i = N or, alternatively, one can use the "master formula" as given in Eq. (16) of Ref. [84], once the variables are transformed according to Eqs. (12).
Obtaining the functions f N (r, ξ) was the goal we wanted to reach. These distribution functions form the basis to compute all properties of the DM species for a given point in the parameter space, see Ref. [84] for details. Not only can we compute the DM abundance Ω N h 2 , we can also use f N (r, ξ) to derive predictions for cosmic structure formation, which can then be matched to observations. We will explain what to do with f N in the next subsection. For instance, the particle number density of sterile neutrinos is computed as expected, where g N = 2 counts the spin degrees of freedom of the sterile neutrino N . From the particle number density, the Dark Matter abundance follows as where n (r prod ) and s (r prod ) are the number and entropy desnsities at r = r prod , respectively, s 0 = 2891.2 cm −3 [108] is today's entropy density, and ρ crit /h 2 = 1.054 · 10 −2 MeV cm −3 [108] is the critical density in units of the squared reduced Hubble constant h.

Relevant bounds
By solving the set of Boltzmann equations for the scalar S and for the sterile neutrino N , we obtain the momentum distribution function of the sterile neutrino, which contains all relevant information and must therefore be in accordance with a number of (partly model-dependent) observational or experimental constraints. In the following, we will present how this is tested and which assumptions go into these tests.
Mass of the sterile neutrino Since we neglected active-sterile mixing (as argued in Sec. 2 and in [60]), the momentum distribution function f N does not yet contain any information on m N , but only on the number density of steriles present in the Universe. We can, however, fix the mass of the sterile neutrinos by demanding that they make up all (or a certain part) of the cosmic DM energy density, cf. Eq. (17). This will later on result into the allowed regions marked in our plots.
Tremaine-Gunn bound The mass of the sterile neutrinos is also restricted by an absolute lower limit, the so-called Tremaine-Gunn (TG) bound [109], which arises from fermions having a maximal phase space density. When applied to astrophysical ob-jects such as the central regions of galaxies, the resulting mass turns out to be m N 0.5 keV [25]. This bound will result into a relatively large excluded region in our plots, and it will in practice separate the regions of the scalars freezing out and freezing in.
Overclosure In principle, even for the lowest possible mass of m N = 0.5 keV, a too large sterile neutrino number density will overclose the Universe (i.e., lead to Ω tot > 1). Since this is not in agreement with observations, a certain patch of the parameter space is excluded by this bound. However, the overclosure bound turns out to be irrelevant in practice, since it is inferior compared to the TG bound.
Structure formation A simple estimator for the predictions of structure formation is the so-called free-streaming (FS) horizon λ fs [25], which yields an estimate of the average length scale (usually given at redshift z = 0) that a DM particle would have travelled from the time of production t prod until today (t 0 ) -had it not been gravitationally trapped. It can be computed via A free-streaming horizon below something like 0.01 Mpc indicates a scenario that -in terms of structure formation -cannot be distinguished from cold Dark Matter (CDM), although the border is arbitrary to some extend (but reasonably well motivated to be about one order of magnitude below the border to HDM). A value of 0.1 Mpc, the typical scale of dwarf galaxies, conventionally indicates the cross-over from models that suppress power on small scales (compared to CDM) in a way that is still in accordance with observational data from HDM models that suppress too much power. As already mentioned, in the literature this regime is mostly referred to by the very unfortunate term warm DM, which seems to indicate a thermal spectrum. However, as we will see in the following, realistic distribution functions are in many cases (highly) non-thermal, which may be very badly reflected by such an over-simplistic quantity like the FS horizon. It is essential to bear in mind that λ fs is -at best -to be understood as an order-ofmagnitude estimator. In fact, depending on the case under consideration, it might have to be modified by an O (1) factor to give more accurate classifications of models. Moreover, by its definition in Eq. (18), it does not take into account the full spectral form of the DM distribution and will therefore never accurately treat non-thermal distributions. We have nevertheless exemplified this analysis (and its failure) for one selected mass of the scalar for the sake of a comparison, cf. App. C. A more robust analysis uses the linear matter power spectrum P (k), which reflects the full spectral information. For every point in the parameter space, we therefore compute the linear power spectrum using the CLASS code [110,111] for non-cold DM species, and then normalise it to the linear power spectrum of a perfectly cold distribution. This ratio, P (k) /P CDM (k) is also known as squared transfer function T 2 (k), which gives information on how strongly the power (i.e., the number of halos formed) is suppressed compared to a perfectly cold distribution at the scale k. Another advantage is that this quantity can then be compared to limits on the squared transfer function obtained from Lyman-α data, displayed by a limiting squared transfer function T 2 lim (k). Limiting functions T 2 lim (k) are usually obtained assuming a thermal DM distribution (in that case truly warm DM). Accordingly, it is for principle reasons difficult to compare a given model to the observational limits. Ideally, one would re-analyse the Lyman-α data using the exact spectral shape of the DM distribution and the respective mass of the sterile neutrino. However, this is impossible given the virtually infinite number of possible distribution functions. Nonetheless, if T 2 (k) ≥ T 2 lim (k) for all k, we can safely categorise the model as in agreement with the Lyman-α bound applied. If, on the other hand, T 2 (k) < T 2 lim (k) for all k, the model clearly has to be discarded. There is a problem with this procedure, though. In some cases, T 2 (k) ≥ T 2 lim (k) might only hold true for a certain interval of wave numbers k, but not for the entire range. For example, T 2 (k) ≥ T 2 lim (k) could be true for all k smaller than a particular valuek, but not anymore for k >k. This complication originates from the very fact that the limiting transfer functions are derived from thermal spectra, while realistic spectra can be highly non-thermal, which results in a difference in the exact evolution of the squared transfer function around its cutoff (e.g., the slope may be different).
In order to still use the limiting squared transfer functions, which are very good estimators accounting for observational relevance, we have developed the following procedure: 1. First, we compute the half-mode k 1/2 , i.e., the wave number at which the squared transfer function has dropped to 1/2: 2. Subsequently, we check whether If the condition is met within this range of small wave numbers (i.e., larger length scales), we consider the model as being consistent with a given bound. This is clearly an approximate method, since we intrinsically disregard some relevant information (namely the power spectrum below the half-mode). Furthermore, the value of 1/2 is, ultimately, pure choice and one could equally well justify other values. However, it is clear that the border has to be somewhere between the two extreme cases of a certain squared transfer functions being completely untouched by the Lyman-α bounds (restrictive view) or by at least in small parts being allowed by the bound (conservation view), which is why a value of 1/2 appears to be a fair compromise.
We have in addition shown that our results are, in fact, quite robust with respect to this arbitrary choice: due to the power spectra considered still having a slope somewhat similar to that of a thermal relic, even a choice of 0.05 instead of 0.5 would not alter our results very significantly, cf. App. D. We are currently working on more advanced analyses that will result into even more reliable procedures.
This procedure is in fact rather simple, as can be seen from the cartoon-like illustration in Fig. 2. For the Lyman-α bounds, we use the squared transfer functions derived from thermal spectra with thermal masses of m lim = 2.0 keV (m lim = 3.3 keV) in a conservative (more restrictive) scenario. Note that only the restrictive bound is displayed in Fig. 2, as the purpose of this figure is only to illustrate the principles behind the procedure, while later on we will display both bounds for completeness. The values quoted are motivated in [24], where the authors also provide an analytical fit formula for the transfer function of a given thermally distributed species of mass m lim , which is assumed to make up all the DM.
Bounds from the observed amount of dark radiation The combination of the momentum distribution function of the sterile neutrino and its mass is also constrained by the observation of the cosmic microwave background (CMB) and by the light element abundances related to big bang nucleosynthesis (BBN). The corresponding observables allow to constrain the effective number N eff of neutrinos and its deviation from the SM value of 3.046 [112]. The effective number of neutrinos is a measure of the radiation present in anything else than photons at a given epoch. It is therefore a measure of the Universe's expansion rate, which in turn leaves its imprint on both the CMB and the abundances of the elements produced during BBN. A similar analysis was shown in [84, fig. 9], however, with a small error in the numerical computation of the contribution to ∆N eff from the non-cold sterile neutrinos. While this error led to an overestimation of their impact, these bounds would in any case only be relevant in the region that is already excluded by the requirement of DM not being classified as "hot". Of course the precise numbers have to be computed numerically, but it can easily be understood that the bounds from structure formation and from N eff must be closely related, since they both critically depend on how long DM stays ultra-relativistic.
Note that, apart from possibly impacting the region in the parameter space where the DM is produced extremely late, the dark radiation bound could in fact also be relevant for very small DM masses. However, as we have seen, any value of m N below 0.5 keV is already excluded by the TG bound.
Model-dependent bounds on the most minimal particle physics setting So far, we have fixed the mass of the sterile neutrino by matching its number density to the observed DM energy density. In the most minimal setting, we could simultaneously demand that the mass of the sterile neutrino is generated by a VEV of the singlet scalar S, i.e., m N = y S . This would allow us to use bounds from perturbative unitarity and bounds derived on the scalar mixing angle from its contribution to the W -mass. In principle, there are also direct collider bounds [113], but they are off our plots and do in practice not constrain the relevant region of the parameter space.
Following [113], we can bound the VEV of the scalar by: Using S = m N y enables us to conclude an upper bound on the Yukawa coupling: A singlet scalar mixing with the Higgs via its VEV can also yield a radiative correction to the W -boson mass. Combining Eqs. (8), (9), and (11) from Ref. [113], we can deduce an upper limit on the Higgs portal coupling: where v EW is the VEV of the SM-Higgs and α is related to the mixing of the scalars.
This completes our list of bounds which can possibly be relevant. As we will see, apart from the Lyman-α bound, the most relevant constraint arises from the TG bound. In the most minimal setting, collider-related constraints can become relevant, too, however these are strongly model-dependent and in fact very easy to circumvent. Finally, the overclosure and dark radiation bounds play no role in practice.

Results
In this section, we present our main results for sterile neutrino Dark Matter. However, before discussing the sterile neutrinos themselves, it is very useful to understand the evolution of the singlet scalar in the early Universe, which will make many of our results much easier to grasp.
So let us, just for a moment, assume that the scalar is stable and understand what happens to it. As we had already mentioned, the scalar can either freeze in or freeze out in the early Universe, depending on its interaction strength λ. We have illustrated all the different cases in Fig. 3, where we display the ratio between interaction rates Γ int of the singlet scalar and the Universe's expansion rate H, as functions of the time parameter r = m h /T . The interaction rates Γ int are computed using all diagrams available in the different regimes, cf. Sec. 2, which depending on the processes may receive contributions from different scattering or decay diagrams. As can be concluded from Tab. 1, all diagrams are proportional to λ, which is why the evolution of the ratio Γ int /H with r looks identical for every plot in Fig. 3, up to a rescaling. However, it is an important information whether Γ int /H < 1, in which case the scalar freezes in, or whether Γ int /H > 1 and falls off only later on, in which case the scalar first thermalises and then freezes out.
Freeze-in Let us start with the freeze-in case. This regime does depend somewhat sensitively on the scalar mass, which is why we have in Fig. 3 depicted the two cases of λ = 10 −9 (left panel : for m S = 30, 60 GeV) and λ = 10 −8 (central panel : for m S = 65, 100, 500, 1000 GeV). In all cases, we can clearly see that the scalars are always far Figure 3: Interaction rates Γ int compared to the Universe's expansion rate H, as functions of the time parameter r = m h /T . We illustrate two freeze-in cases -depending on the scalar mass we show either λ = 10 −9 (left) or λ = 10 −8 (center ) -as well as one freeze-out case, with λ = 10 −5 (right). Notably, in the latter, all scalars undergo a cold freeze-out a temperatures within [m S /20, m S /4], which explicitly proves that it is a good approximation to consider the scalar to be at rest at the point of freeze-out [70], contrary to the claim made in [114]. from equilibrium. 7 We can furthermore see the change in the interaction rates due to the new diagrams available after the EWPT. In particular for very small scalar masses, m S = 30 or 60 GeV, we can see that the interaction rate jumps up by some two orders of magnitude. This is because for such small masses the Higgs decay h → SS is available, which is the dominant contribution as noted in Ref. [88]. For larger scalar masses, in turn, there is hardly any difference visible (if at all, the interaction rates seem to decrease a little); however, the main reason is that larger scalar masses are thermally suppressed at this late stage, as most clearly visible for the largest scalar masses of m S = 500 or 1000 GeV. But even for medium-scale masses of m S = 65 or 100 GeV, a small drop in the rate is visible. This is in parts due to kinematics, since the W -, Z-, and Higgs bosons as well as the top quark are already thermally suppressed at this stage, cf. diagrams in Tab. 1, but it happens also because of the different structures of the diagrams involved. In any case, the freeze-in of the scalar ceases at a temperature that is roughly equal to its mass, and indeed all interaction rates start to drop shortly after (or even before) T = m S is reached. A larger coupling λ translates into a larger number density, as expected.
Freeze-out For the larger coupling, λ = 10 −5 , scalars of all masses equilibrate. As visible in the right panel of Fig. 3, all interaction rates are larger than the expansion rate already for very early times. As long as the scalars are equilibrated, the actual interaction rate does not matter much. However, a boost in the rate -as due to the EWPT for m S = 30 and 60 GeV -can "delay" the freeze-out. This is particularly visible for a m S = 60 GeV, which freezes out at r ≈ 41 (or at T ≈ m S /20), and thus even later than the scalar with m S = 30 GeV, for which freeze-out happens at r ≈ 30 (or T ≈ m S /7). Thus, as we will later see to be correct, we can expect a much larger number density of scalars (and hence of sterile neutrinos) for m S = 30 GeV than for 60 GeV. Similarly, the case of m S = 65 GeV should yield a larger abundance than the one for m S = 100 GeV, which also turns out to be correct.
Sterile neutrinos Let us now come to sterile neutrinos as DM. In order to include as much information as possible in the figures, without however making them too busy, we illustrate the allowed ranges and all relevant constraints as "spaghetti plots" displayed in Figs. 4, 7, 10. In these plots, we depict the regions in the λ-y parameter space where the correct DM abundance is met (where only part of the abundance is produced) by the dark coloured solid lines (by the lightly coloured bands), for different values of the sterile neutrino mass: m N = 2, 7.1, 20, 50, 100 keV. To obtain these regions, we have for each combination of the parameters (λ, y, m S ) computed the resulting distribution function along the lines described in Sec. 3, and then integrated the result to obtain the DM abundance according to Eq. (17).
In addition, we display the TG phase space bound, as explained in Sec. 3.2. The overclosure bound (gray area) marks the part of the parameter space where we would obtain overclosure of the Universe even for the minimum sterile neutrino mass of 0.5 keV.
We also display the model-dependent bounds, which only hold in the most minimal setting, cf. Sec. 3.2. While they may not even exist in more involved settings, they do serve the purpose of representing the maximal effect collider-related bounds could possibly have. To indicate some numbers, the effect of the unitarity bound, Eq. (21), for a scalar mass 100 GeV is that the Yukawa coupling y has to be smaller than about 8.  22), is somewhat more subtle, since the scalar mixing angle α depends on the Higgs portal λ which in turn influences the DM abundance and is thus related non-trivially to the Yukawa coupling y and m N . As discussed, these bounds do only exist in the most minimal model, which is why we have marked them by thin dashed lines (see areas left of the TG bound and down in the lower right corners). For some cases, they are even off the plot.
The main bound however arises from structure formation, where we take into account the two Lyman-α bounds derived in Ref. [115] and perform the half-mode analysis described in Sec. 3.2 to determine whether a given distribution function is consistent with the data, or not. We have for each scalar mass m S numerically computed the linear power spectrum for each pair (λ, y). Reproducing the correct abundance results in a condition on the sterile neutrino mass m N , so that every point with the correct abundance can be characterised by a point (m S , λ, y, m N ) in the parameter space. Depending on which Lyman-α bound we used, we have in our plots marked the following regions: • forbidden: If the upper half of the squared transfer function is forbidden by both the conservative and restrictive Lyman-α bounds, the respective point is red.
• constrained: If the upper half of the squared transfer function is only forbidden by the restrictive Lyman-α bound but allowed by the conservative one, the respective point is purple.
• allowed: If the upper half of the squared transfer function is allowed by both the conservative and restrictive Lyman-α bounds, the respective point is blue.
This colour code is slightly reminiscent of the historically grown terms "hot", "warm", and "cold" DM, however, we would like to stress once more that the distributions obtained are non-thermal, and one thus cannot associate any temperatures with them; see App. C for an explicit counterexample showing that the classification drawn from thermal spectra is bound to fail for non-thermal distributions. Nevertheless it is correct to say that, roughly, the red regime is associated with rather high momenta compared to the DM mass, while the blue one tends to feature smaller momenta. By the light colours, we have indicated regions with a sizable but insufficient abundance for a given mass. Technically, these scenarios correspond to a slightly larger mass which yields the same effect as replacing a fraction of, say 10% of the Dark Matter by a perfectly cold Dark Matter component. We have checked that this simplified procedure yields correct results up to the level of accuracy limiting analyses using Boltzmann equations anyway. A final comment on the unavoidable but tiny contribution by Dodelson-Widrow production at temperatures of a few 100 MeV: as had been shown in Ref. [60], this modification does not only hardly contribute to the abundance for sterile neutrino masses of 4 keV and higher, also the modification of the actual spectrum is completely negligible. Thus, in our plots, there would be no effect visible even if we did take into account the largest modification allowed by data. This may only be different for a sterile neutrino mass of 2 keV, but even then the error by not considering the modification is of about 10%, which is in practice still negligible. For even smaller masses, the correction could be sizable, but in such cases a mixed scenario would be excluded by structure formation anyway.

Very light scalars: m S < m h /2
Let us begin with the case of very light scalar masses, i.e., m S < m h /2, displayed in Fig. 4. As discussed in Sec. 2, this leaves us with an interval of roughly 5 GeV m S 62 GeV. Here, the lower limit arises from the temperature where the relevant interactions are frozen out and/or not efficient anymore, while the upper limit is obtained from m h 125 GeV. Here, depending on the coupling λ, we may enter different regimes. First of all, for very small couplings λ, we will be in the FIMP region, i.e., in the left part of the plots. Freezein generically is most efficient at temperatures T around the mass m S of the FIMP, i.e., in regime III (see Fig. 1 and Tab. 1). Of course the production of scalars has started in regime I, at high T , but in any case -even with all diagrams of regime III at work -the interaction remains feeble. Once the scalar is produced, it decays with a certain lifetime proportional to y −2 into sterile neutrinos. Thus, for y large, the decays happen very close to the freeze-in of S and thus the resulting DM particles will be more strongly redshifted and "cooler". For too small y, in turn, the scalars remain present in the Universe and decay only very late, thus injecting too much energy into the sterile neutrinos and rendering them HDM -and thus excluded. Notably, as can be seen in the plots, the actual mass of the scalar within the allowed interval does not make much of a difference. The decay rate is proportional to y 2 m S , such that a larger scalar mass m S should translate into more red-shift, both because the correct abundance is obtained earlier and and on top of that the decay happens faster. However, for a fixed y, the unboosted decay rate also becomes smaller, such that the lifetime effectively increases, leaving less time for the steriles to redshift. To a first approximation, these two effects compensate each other, but of course there will be corrections due to scalars of different masses being boosted more or less (and thus having bigger or smaller lifetmes in the cosmic frame). Apart from these minor effects, the FIMP regions for m S = 30 GeV and m S = 60 GeV appear rather similar, up to a shift in the Higgs portal λ.
On the other hand, for larger values of λ, the scalar can enter thermal equilibrium and then freeze out, see the lower right regions of the plots. Depending on the exact combination of (λ, y), the decay into sterile neutrinos happens while the scalar is in equilibrium (horizontal part of the isoabundance-lines -which is independent of λ from a certain value on), after freeze-out (vertical part -which is similar to freeze-in in the sense that the scalars cool down before their decay and thus "forget" about their cosmological history for a long enough lifetime or small enough coupling y), or in between (kink-region -where a strong dependence on both λ and y is present). For this part of the plot, the question of whether there is at all any allowed region depends strongly on m S . This observation can be understood by realising that the scalar freezes out earlier and thus with a much larger abundance for m S = 30 GeV, which is due to the dependence of the interaction rates on the scalar masses. Hence, given many more scalars than for m S = 60 GeV, a too large abundance can only be avoided for sufficiently small sterile neutrino masses -which is why only one strip is visible in the lower right corner of the left plot, and this one is forbidden by structure formation due to the small DM mass favouring a "hotter" spectrum, i.e., with a stronger tendency for large momenta.
Two representative example evolutions of the yield can be seen in the top row of Fig. 5, both for the scalar freezing in (left) or out (right), which correspond to the crosses marked in the right Fig. 4. If the scalar was stable (gray dashed line), it would compare trivially to the hypothetical equilibrium curve (gray dotted line): as a FIMP (left), it would gradually be produced but never reach a thermal abundance before freezing in, while as WIMP (right) it would quickly thermalise before freezing out. Once the decay into sterile neutrinos is switched on, the true abundance of the scalar (black solid line) again increases initially but ultimately decreases to zero. Note that, in both plots, one can see the increase in the interaction rate due to the EWPT, cf. Fig. 3. The sterile neutrino   Let us make a quick comparison to the results obtained in Ref. [88]. Looking at the top two plots in their Fig. 2, which feature the same scalar masses as our top two plots in Fig. 5, our yields look qualitatively very similar except for a missing "bump" in the yield of the scalar for temperatures around 10 GeV, which appears in Ref. [88] but not in our result. While we cannot explain this feature, we do however have good agreement with the overall abundance obtained in [88], 8 up to some 20% difference, which can be attributed to applying different numerical procedures and also due to Ref. [88] using some approximations such as rate equations assuming a (suppressed) thermal shape.
The evolution of the resulting sterile neutrino distributions for the same two points is displayed in the bottom row of Fig. 5, with all spectra being functions of ξ = gs(T 0 ) 1/3 p T , cf. the second Eq. (12). We can see a certain deformation in both curves, although the one on the left (for the scalar being a FIMP), it is probably a small bump rather than an actual double peak. In both cases, these shapes come from two different production phases. In the FIMP case (left), the two phases correspond to the production before and after the EWPT, whereas the visible bump for larger ξ for the WIMP case (right) rather comes from production while the scalar is still in equilibrium and after its freeze-out, respectively. 9 Finally, in Fig. 6, we display the squared transfer functions for the two example points. Glancing at the right Fig. 4, we can see that both points lie inside "purple" regions, i.e., they are constrained but not excluded by the Lyman-α bounds. According to the procedure detailed in Sec. 3.2, this would imply that the parts of the squared transfer functions with k ≤ k 1/2 should both be consistent with the conservative Lyman-α bound, but inconsistent with the restrictive one. This is just what we can see in Fig. 6. It is also visible that, in particular for the FIMP case depicted on the left, the slope of the functions is different from that of the bounds, for which the distributions were thermal by construction. Furthermore, note that the sterile neutrino mass on the left plot, 7.1 keV, is significantly different from the mass of 2 keV corresponding to the conservative bound, which originates from the non-thermal shape of the scalar-decay produced DM spectrum. Obviously, according to the left Fig. 6, one would classify the case displayed as "warmer" than the restrictive bound corresponding to a mass of 3.3 keV, even though the true DM mass is even larger than that. This is one more reflection of conclusions which do hold for thermal spectra being invalidated once the shape of the momentum distribution function is more complicated. Thus, any "translation" of the structure formation properties to a hypothetical spectrum of thermal shape, as attempted for DW-produced sterile neutrinos on several occasions in the literature (e.g. Refs. [26,98,[115][116][117][118][119]), must be treated with extreme care and should not be regarded as a precision statement. We can also see that, although both of the points that we had explicitly illustrated above fall into the same category in what concerns structure formation, the FIMP case tends to be slightly "warmer", i.e., closer to being excluded than the WIMP example. This may come as a small surprise, given that the average momentum of the WIMP-produced sterile neutrino is nearly two-times larger than the FIMP-produced one. 10 However, when dividing the average momenta by the sterile neutrino masses, hence looking at the average velocities, the factor of nearly two remains present -just that in this case the FIMPproduced particle yields the larger value.
Looking back at Fig. 4, to get a more global picture, it is visible in the freeze-in (left) regions that for larger Yukawa couplings y the DM particles get "colder" (i.e., shifted to smaller momenta). This is due to the parent particles decaying earlier, thus leaving more time for the DM particles to redshift. Of course, this effect is more pronounced when the DM particles are heavier, see upper left corners of the plots. Interestingly, the unitarity bound would cut off part of that otherwise allowed parameter space (but only if m N = y S ). For the freeze-out (right) regions, instead, the case of small Yukawa coupling corresponds to a very late decay, always far after the freeze-out. This leaves the DM particles no time to cool down, so that they will generically be threatened by structure formation. The second limit of large λ keeps the scalars in equilibrium for a very long time, so that practically all sterile neutrinos are produced while the scalar is tiny that it does not have a real influence on the spectrum. still equilibrated. That leads to a "colder" spectrum, because the DM particles have more time to redshift. The turnover, corresponding to the double peak structure observed in Ref. [84] for the first time, can be somewhere in between. Depending on the (non-trivial) interplay between the different parameters and in particular on the relative strength of the two peaks, it can be ruled in or out by structure formation.

Light scalars: m h /2 < m S < m h
The next case to be discussed, cf. Fig. 7, is that of "light" scalars, in the sense that their mass is still less than the mass m h of the Higgs boson, but larger than half the Higgs mass. This corresponds to regime II in Tab. 1 for low temperatures T , while the production nevertheless starts at high T and thus in regime I, cf. Fig. 1. Here, the main characteristic is that, at the EWPT, all of a sudden lots of new channels open up. For the freeze-in case, i.e. small λ, this means that it is easier to produce scalars at a temperature just above their mass, which will result in an increase in their abundance and thus also in that of sterile neutrinos -however this increase is much less dramatic than for the very light scalars, and hence hardly visible in the plots, due to the Higgs decay into two scalars not anymore being accessible in this mass range. For the freeze-out case, i.e. large λ, this means that the scalars will be kept in equilibrium for a longer time in regime II, resulting into a large overall number of sterile neutrinos produced while the scalar is still equilibrated.
The spaghetti plots for two different scalar masses, m S = 65 and 100 GeV, are depicted in Fig. 7. At first sight, these plots appear qualitatively similar to those depicted in Fig. 4, except for an overall shift towards larger values of λ. However, an interesting observation is that, for the FIMP case, the spectrum seems to become "warmer" when going from m S = 65 to 100 GeV, while the trend was opposite when going from m S = 30 to 60 GeV. The reason is that, for the case at hand, Higgs decay does not play a role in the production of the scalar. Instead, many scalars are already produced rather early (see upper left panel of Fig. 8), such that they can indeed decay earlier (and are less boosted), so that in this case a smaller scalar mass results into smaller initial sterile neutrino velocities and this effect dominates the production.
As already anticipated, the lighter scalar (m S = 65 GeV) freezes out earlier, namely at T ∼ 17 GeV, than the heavier one (m S = 100 GeV), which does so at T ∼ 15 GeV, cf. Fig. 3. This may look surprising at first, since usually heavier particles tend to freeze out earlier due to the thermal suppression for non-relativistic particles. But the thermal suppression is similar for both cases, e −65/100 ∼ 0.5, so that the actual interactions are decisive. Already at T ∼ 50 GeV, all heavy bosons and the top are thermally suppressed, so that they cannot easily be produced anymore. However, a scalar with a mass of 100 GeV can even at rest annihilate into both W + W − and Z 0 Z 0 , while a lighter scalar cannot. Thus, the interaction rate of the heavier scalar is considerably larger, thereby keeping it in equilibrium significantly longer. It follows that the number density of the lighter scalar, and thus of the sterile neutrinos originating from its decay, is higher and hence must be compensated by a smaller sterile neutrino mass -which is the origin of only small masses being allowed in the freeze-out region of the left Fig. 7.
Again, two example points are depicted in Fig. 8, with the order of the plots as in Fig. 5. It is clearly visible that, in the yields of the scalars, the effect of the EWPT is much less dramatic than for the very light scalars. In the resulting sterile neutrino spectra, the effect of the EWPT is thus much less pronounced, although little bumps can be spotted on the left sides of both spectra in the bottom row of the figure. Qualitatively, the course of the DM production is similar to case of very light scalars: again, if the scalar freezes in, it gradually decays into sterile neutrinos; if it freezes out, then it may decay either in or out of equilibrium, which a gradual transition between the two regimes. Note that, for the masses displayed, the comparison to Ref. [88] is in fact much more difficult. While the same scalar masses are used in that reference, we strongly suspect that paper to have a small error in the degrees of freedom of the Higgs field above the EWPT. While the Higgs field does possess four physical components in the unbroken phase, in Ref. [88], only one seems to have been used. In other words, at high T , while the gauge boson 10 -5 10 -5  contributions are switched off, the contributions of the would-be Goldstone bosons are not switched on, as one would need to do. We will confirm this claim later in Sec. 4.3, where it results in a factor of four between our results and those obtained in Ref. [88]. However, this factor is only close to four for very heavy scalar masses (even when taking into account the difference in the normalisation of λ, which also happens to be a factor of four), where the production almost entirely happens in regime I. But for the comparatively light scalars treated in this subsection, the production is spread out over regimes I and II, cf. Fig. 1, such that a discrepancy by less than a factor of four but no perfect agreement can be expected. This is confirmed by our numerical results, making us confident that our treatment is the correct one for high temperatures. The squared transfer functions (i.e., the ratios of the resulting linear power spectra to the CDM case) are depicted in Fig. 9. As to be anticipated from the left Fig. 7, the FIMP point is located in a blue region, and it should thus be perfectly allowed by all bounds. The WIMP point, in contrast, is in a region that is entirely red -and it should thus correspond to DM that is by far too hot. Indeed, the squared transfer functions displayed in Fig. 9 confirm this conclusion. The whole curve of the FIMP case is in the allowed region, whereas the whole curve of the WIMP case is forbidden by far, even by the conservative bound. For such clear cases, an elaborate analysis would not even be necessary. Again we observe that, although both points feature a sterile neutrino mass of 7.1 keV, the predictions for cosmic structure formation are different by far. This clearly illustrates that the DM mass itself is not at all a good way to discriminate between two different DM spectra, even if viewed at one and the same temperature.

Heavy scalars: m h < m S
Finally, the scalars may also be heavier than the Higgs, see Fig. 10, which is the case that had already been discussed in our earlier reference [84]. In this case, given that both freeze-in and freeze-out are correlated to the mass of the scalar, most of the production of scalars happens above the EWPT -we basically remain in regime I in Fig. 1 during the whole production. Only for somewhat light scalar masses, such that m S /20 T EWPT , a small amount may be produced after the EWPT. However, for most of the production, only the simple 4-scalar interaction in the top row of Tab. 1 plays a role. In particular, one can greatly simplify the equations involved for (m h /m S ) 2 1, see Ref. [84] for details. As we can see from the central and right panels of Fig. 3, it is in fact a good approximation to assume that the production of sterile neutrinos from very heavy scalars happens entirely before the EWPT -thereby providing a clear justification of the assumptions made previously in [84].
Understanding the results for this case is somewhat less subtle than for the previous two. Comparing the two plots depicted in Fig. 10, one can see that there are no overly drastic changes. The reason is that, as argued in Ref. [84], there is a trade-off between a heavier scalar and an earlier decay: a larger scalar mass means that the production of scalars ceases at higher temperatures. In the freeze-in case, this indicates a lower overall production, which is what the vertical lines in that regime slightly shift to the right when going from m S = 500 GeV to 1000 GeV. However, given that no new production channels open up while staying in regime I, the change is rather mild. In the other limiting case, where the scalar freezes out, a larger scalar mass means a larger abundance for the case where the scalar decays mainly after freeze-out (i.e., the vertical lines on the bottom right of the plots). Thus, these vertical lines should slightly shift to the left to compensate for that, which they indeed do, as visible in the plots. For the horizontal parts of the lines, where the scalar decays while in equilibrium, however, freezing out earlier means less sterile neutrinos unless this is compensated by a larger decay rate: thus, these lines should shift to slightly larger values of the Yukawa coupling y when going from m S = 500 GeV to 1000 GeV, which is just what is visible in the plot.
Two representative example evolutions of the yield can be seen in the upper row of Fig. 11, both for the scalar freezing in (left) or out (right). Just as before, in the FIMP case, the scalar is gradually produced but never reaches a thermal abundance before freezing in (while constantly decaying), whereas for the WIMP-case it quickly thermalises before freezing out. Also here, the decay can happen earlier or later, depending on the value of y. As already hinted in Sec. 4.2, we expect a factor four difference compared to the results obtained in Ref. [88], as DM production happens nearly completely in regime I. This approximation is best for very large scalar masses, however, the most extreme case discussed in [88] involves a scalar with m S = 500 GeV. Taking the same couplings as in the reference and taking into account the different normalisation of λ, the discrepancy obtained numerically is already larger than a factor of three, and should converge to four for even larger scalar masses (which we however cannot explicitly check, given the absence of this case in [88]). Yet, given that no distribution functions had been computed in Ref. [88], their results can be rectified by simply including the missing factor.
The corresponding evolutions of the distribution functions are depicted on the bottom row of Fig. 11, with similar characteristics as for smaller masses. One big difference, though, is that the DM production is finished much earlier, which comes from heavier scalars decaying faster. This allows more time to redshift, but it is compensated to some extend by a larger initial velocity.
In what concerns structure formation, however, these shifts do not change very much. In the freeze-out case, while the sterile neutrinos are produced by a heavier scalar and are thus faster, they are also produced earlier (since the decay rate of the scalar is proportional to its mass), so that hardly any change in the colouring can be seen for the freeze-in case when comparing the two plots in Fig. 10. For the case of the scalar freezing out, only the in-equilibrium decays are phenomenologically relevant, as the late decays in any case produce too "hot" Dark Matter. However, for scalars decaying while in equilibrium, a larger mass translates into larger Yukawa couplings and thus into earlier decays so that, again, the resulting sterile neutrinos have a larger momentum at production but they have also more time to redshift, such that the prediction for structure formation remains basically unaffected -hence the similarity between both plots. Note that, for large scalar masses, none of the "kink" regions on the bottom right of the plot is allowed, which would 10 -5 which analysed just that case, and it is also backed up by the example squared transfer functions, cf. Fig. 12.

Conclusions & Outlook
In this work, we have completed the investigation of sterile neutrino Dark Matter production from scalar decays. We have shown how to perform the full computation using Boltzmann equations on the level of momentum distribution functions. While this would in principle amount to numerically solving a coupled system of partial integro-differential equations of a rather uncommon (and possibly unknown) type, we have argued why it is allowed to decouple the equations, how to transform the partial into an ordinary integrodifferential equation, and which steps to take to finally solve it using an algorithm based on a discretisation in one variable. All that comprised pioneering work, which had previously not been available without drastic assumptions or simplifications. While our findings are of course somewhat specific to the examples studied, our methods carry over to more general settings such as decays of non-scalar particles or multiple generations of sterile neutrinos: any reader inclined to investigate such scenarios can adopt our methods and will be able to apply the same strategies to crack the equations involved. The importance of investigating sterile neutrino Dark Matter on the level of momentum distribution functions lies in these quantities carrying all the information about the Dark Matter model at hand. Not only can they be used to extract binary information such as the Dark Matter number or energy densities, but having the distribution functions is vital to derive predictions for and apply constraints from cosmic structure formation. Many previous works tried to go around this by solving oversimplified rate equations and using a simplistic estimate of the free-streaming horizon to conclude about whether or not a setting is consistent. However, we have shown that this procedure can fail once nonthermal distribution functions are involved. In these cases, it is just not a very reliable strategy. We have instead proposed to compute the linear power spectrum (from which we derive the squared transfer function), for which task dedicated tools are available, and to use the half-mode analysis developed in this work to obtain a fairly good estimate of the compatibility of a certain distribution function with data [120].
Having investigated the production of sterile neutrino Dark Matter from scalar decay in full generality, the final question is where to go from here. While we have shown how to obtain a reliable constraint from Lyman-α data using the half-mode analysis presented in Sec. 3.2, in principle, one could refine the analysis presented, as we intend to do, and/or study further halo properties of settings like the one presented, as done e.g. for mixed Dark Matter settings [121,122]. The distribution functions obtained from scalar decay are conceptually not very different, however, as we have seen, they are highly non-thermal in both shape and number of momentum scales, so that even a description by a superposition of two thermally-shaped distributions may in fact not be sufficient. Various constraints can be applied, with maybe the most generic apart from the Lyman-α forest being dwarf satellite galaxy counts [123]. 11 The main limitation of our method is that, currently, we are using only the linear power spectrum [124], while more considerable differences of decay-produced Dark Matter to thermal spectra may be visible in the non-linear regime. Ideally, one would perform and N -body simulation for each distribution presented here, however, this is clearly not doable given that these computations are numerically very expensive. Thus, a better strategy is to first obtain a pre-selection of cases which could be potentially interesting for a full N -body simulation, while discarding those which are clearly ruled out, or in practice no different from cold Dark Matter. One way to do such an analysis is to estimate the effect of the non-linearities, by the extended Press-Schechter approach [125,126], modified by incorporating a collapse model that describes how the free-streaming Dark Matter particles get bound to form structures [123,127]. This is the next step to be done, which we will leave for future work [120]. Given that present-day data is already sufficient to yield vastly different constraints on different sterile neutrino Dark Matter production mechanisms [70] it can be expected that -with more detailed computations -we may be able to clearly distinguish various production mechanisms by observational data.
general collision term. The most general form possible for a collision term describing the reaction ψ + a + b + ... ↔ α + β + ... is [129]: withp being the 4-momentum and E p the energy of the particle ψ under consideration. The symbol |M| 2 denotes the (initial and final) spin-averaged matrix element, which also contains symmetry factors for identical final and initial states, 12 as well as multiplicity factors for reactions involving multiple particles per process. Otherwise, the standard definitions apply: the phase space element is dP X = g X d 3 p X 2E X (2π) 3 , for a particle X with g X internal degrees of freedom, momentum p X , and energy E X = m 2 X + p 2 X with mass m X . To give one concrete example, we will now explicitly derive the collision term describing the decay of a SM-like Higgs into two singlet scalars S, C S h↔SS [f S ](p, T ). The following quantities label the 4-and 3-momenta, as well as the absolute value of the latter: •p, p, p (≡ |p|): for either one of the scalars, •p , p , p (≡ |p |): for the remaining scalar, •q, q, q (≡ |q|): for the Higgs boson.
We furthermore abbreviate E i p j ≡ m 2 i + p 2 j . Starting from Eq. (A-1), we obtain: where |M| 2 = 16λ 2 v 2 , cf. Eq. (A-21). The 3-dimensional δ-function eliminates q: Due to the remaining δ-function, the term in parentheses can be cast into: As can readily be seen, these amount not to four but to only two distinct values: (A-10) (A-11) Employing arguments on continuity and limits of cos α 0 for p → ∞ and p → 0, these two have to be the boundaries of the p -integral, so the final form of this collision term is: (A-12) Having seen one concrete example, we will now list the full set of collision terms needed to reproduce our results: • 2 → 2-scattering processes: All scattering processes are of the form where ii = φφ, 13 hh, tt, W + W − , ZZ and f eq S (p) = exp − p 2 + m 2 S /T is the (would-be) equilibrium distribution of the scalar S, which is of Boltzmann-shape by virtue of the principle of detailed balance. Furthermore,ŝ is the square of the centre-of-mass energy, explicitly given by: s(p, p , cos α, m S ) = 2(m 2 S + (m 2 S + p 2 )(m 2 S + p 2 ) − pp cos α).
(A- 14) In addition, p is the momentum of the second scalar participating in the process and α is the angle between p and p , whose maximum value is given by cos α max = min{max{cos α im , −1}, 1}, defined by 4m 2 i =ŝ(p, p , cos α im , m i , m S ). This upper integration boundary excludes all values of cos α for which the integral becomes imaginary due to the square root contained in Eq. (A-13).
The spin-averaged matrix elements including a factor of 2 in all of the following because two scalars are annihilated or produced, respectively, in each process, are given by (assuming CP -invariance): 14 • 1 → 2-decay processes: Several different decays have to be taken into account. Starting with the decay of a SM-Higgs into two singlet scalars, the corresponding collision term is given by , as well as the matrix element including a factor 2 because of annihilation/production of two scalars: Furthermore, there are two collision terms related to the decay of a singlet scalar S into two sterile neutrinos N . The first is the one used in the Boltzmann equation for S, with the decay widths Γ S→N N = y 2 m S /(16π). The second version is the one used in the Boltzmann equation for N , 4p . Note that the two collision terms C S S→N N and C N S→N N appear to be somewhat different, although one may very naively expect one to be just the negative of each other. However, we should keep in mind that we are working on the level of momentum distribution functions, which implies that the collision terms look rather different depending on whether or not the desired distribution function is integrated over.
In these equations it is understood that p, p , and T are substituted in favour of ξ, ξ , and r, as specified in Eq. (12).

A.2 Transformation of Variables
As stated in the main text, we perform a transformation of variables in order to bring the Liouville operatorL = ∂ ∂t − Hp ∂ ∂p into a more convenient form. To see how this results into Eq. (13), consider a general transformation into new variables r and ξ: These new variables can be inserted into the Liouville operatorL: To simplifyL, we need to get rid of one of the two differential operators. In other words, the new variables would be most useful if we could choose them in such a way that they transform the partial differential equation into an effective ordinary differential equation. The first step is to demand that r does not depend on p, which eliminates one term: Next, we demand This is a rather simple partial differential equation. Fixing the initial condition where ξ 0 is some arbitrary C 1 -function, it has the simple solution This implies that, if we fulfill the requirements that r only depends on t and the dependence of ξ on p and t is given by Eq. (A-29), the Liouville operator in terms of the new coordinates will have the simple formL = ∂r ∂t ∂ ∂r . (A-30) We can make our life even easier by a smart choice for the functions r(t) and ξ 0 . Exploiting the one-to-one correspondence between temperature T and time t, one possible choice is: for some reference mass m 0 and some reference temperature T 0 , both of which we choose to equal the Higgs mass: For the last equality in Eq. (A-31), we have used the fact that the comoving entropy density s is constant,

Appendix B The freeze-out of the Higgs boson
We argued that we do not have to solve a system of Boltzmann equations for all species in the early Universe, but only for the scalar singlet S and the sterile neutrino N , since we rely on the SM particles sourcing the production of S being in thermal equilibrium. This assumption is for sure good for m S m H , but we should assess its quality if we want to proceed to smaller m S . The smaller m S the later (in cosmic time) the scalar will be produced in general, and hence the less reliable the assumption of particles like the Higgs or gauge bosons being in equilibrium could be. In fact, all species relevant for sourcing S in the broken phase (W ± , Z, h, t) have similar masses and are therefore expected to decouple from the thermal plasma at a similar time.
Still, if we restrict ourselves to m S > 30 GeV, it is enough to know whether the source particles are still in equilibrium at the corresponding temperature. Even if they are in equilibrium much longer, all the particles mentioned will by then have disappeared due to their masses resulting in strong Boltzmann-suppressions. Since we have not found any detailed source to cross-check this assumption, we have assessed it ourselves with the following analysis.
Let us assume that the Higgs boson h and the top quark t decouple from the plasma before the gauge-bosons, since their mass is larger by a factor of O (1). We therefore looked at the following system of coupled Boltzmann equations: where SM denotes all SM degrees of freedom except for the t and the h. The matrix elements going into all the collision terms were computed at tree-level only. The Higgs decay width was taken from Ref. [108]. Solving them on the level of rate equations, we find that in this system, both the Higgs and the top closely track their equilibrium abundance all the way down to temperatures of about 1 GeV, where their abundances are exponentially suppressed already. When switching off inverse decays and taking into account only two-to-two-processes, the top would freeze out around 5 GeV while the Higgs would freeze out around 3 GeV. Note that the massive gauge bosons can also be produced via inverse decays, which is another argument to assume that they will decouple only after the Higgs and the top.
Appendix C Failure of the free-streaming horizon as measure for non-thermal Dark Matter Here, we would like to discuss one important point which we had stressed on several occasions throughout the manuscript. Given that we are dealing with highly non-thermal DM spectra, we cannot expect paradigms developed for thermal relics to carry over to this much more general case. In particular, a non-thermal spectrum cannot be described very well by a single number, such as a temperature, an average momentum, or a velocity. Given that, we have in fact no reason to believe that a quantity based on an average velocity, such as the free-streaming (FS) horizon λ fs as defined in Eq. (18), should give us any countable information on the DM spectrum. Yet, it is incorrectly used in many occasions in the literature.
In order to clearly illustrate how the FS horizon fails, we compare two versions of an example spaghetti plot, which are depicted in Fig. 13. Here, on the left panel, we can see the plot for m S = 100 GeV how we obtained it in Sec. 4.2 (this figure is basically Figure 13: Comparison of half-mode analysis (left; figure identical to the right Fig. 7, except for the model-dependent bounds not being displayed here for simplicity) with the computation of the free-streaming horizon (right).
identical to the right Fig. 7, apart from the missing collider-related bounds which we skipped here in order not to distract the reader). The identical patch of the parameter space is displayed on the right panel of Fig. 13, however, this time with an analysis based on the FS horizon, classifying the different points as "cold", "warm", or "hot"even though, as already explained, these categories do not really suit any non-thermal spectrum. Comparing both plots, one can immedately see that the analysis based on the FS horizon 15 is much more pessimistic than the result based on the half-mode analysis described in Sec. 3.2. While there is no perfect correspondence of the colours red, purple, and blue between the two plots, it is nevertheless evident that some points which are not even constrained by current data in the left plot, seem to be excluded completely when looking at the right plot. This is particularly true for the two points marked in the plots: Obviously, both these points are unconstrained (blue) in the half-mode analysis, while point A would be classified as "hot" (red) in the analyses based on the FS horizon, and even point B would still be labeled as "warm" (purple). This second classification is based on the results obtained for the FS horizon, which turn out to be: A : λ fs = 0.117 MPc ⇒ "hot", B : λ fs = 0.089 MPc ⇒ "warm", (C-2) which perfectly matches the classification from the right panel of Fig. 13. But how can we be sure that this classification is insufficient? The simplest way is to confront the DM distribution functions with the actual data, which we can do by explicitly displaying the corresponding squared transfer functions in comparison to the Lyman-α data, as done in Fig. 14. Looking at the curves for both points A & B, we can immediately see that both of them are not at all constrained by the data. Thus, both these points are, in fact, even indistinguishable from cold DM, from a structure formation point of view.
Given that this is truly obvious for the two example cases, it serves as a clear example for the FS horizon analysis leading to a conclusion that would be completely incorrect (namely discarding point A all along, while viewing point B as borderline case). Thus, as should now be obvious, the free-streaming horizon is not at all a suitable measure when applied to non-thermal DM distributions. Figure 15: Compatability of the scalar decay model with m S = 60 GeV for different threshold wave numbers. The left panel contains the analysis with k x = k 1/2 (as used throughout Sec. 4; the plot displayed here is basically identical to the right Fig. 4, except for the model-dependent bounds not being displayed to enable a better comparison), while the right panel displays k x = k 0.05 . The comparison shows that the changes are minor, but they will be finally specified in an upcoming work of us aiming to compared several advances methods to derive bounds from structure formation.
At first sight, the definition of the halfmode in Eq. (19) might seem just as arbitrary as the free-streaming boundary between "cold" and "warm" ("warm" and "hot") at 0.01 Mpc (0.1 Mpc). Even though the transfer function falls off steeply around k 1/2 , it is a priori unclear at which value of the squared transfer function the discrimination becomes indeed negligible. While a more fundamental analysis of structure formation (e.g. by rederiving Lyman-α bounds for non-thermal spectra) is beyond the scope of this paper and is projected for future work, we can nonetheless test the robustness of our analysis against changes in the threshold given in Eq. (19). More specifically, we can make our analysis more restrictive by comparing not only wavenumbers smaller than k 1/2 but wavenumbers smaller than some k x with x < 0.5, which translates to k 1/2 < k x . In order to demsonstrate that the analysis is rather robust even against large changes in this threshold, we have reanalysed the case of a scalar with m S = 60 GeV with a k x = k 0.05 . Fig. 15 shows that even this rather drastic change in the threshold power (by one whole order of magnitude) inflicts rather mild changes on the results. This would be dramatically different when changing the values put in as boundaries in a free-streaming analysis (cf. Fig. 13).
Even though being approximate itself, our method of incorporating bounds from structure formation, using Lyman-α data derived for thermal spectra, clearly proves more reliable than the free-streaming approach.