Dodelson-Widrow Production of Sterile Neutrino Dark Matter with Non-Trivial Initial Abundance

The simplest way to create sterile neutrinos in the early Universe is by their admixture to active neutrinos. However, this mechanism, connected to the Dark Matter (DM) problem by Dodelson and Widrow (DW), cannot simultaneously meet the relic abundance constraint as well as bounds from structure formation and X-rays. Nonetheless, unless a symmetry forces active-sterile mixing to vanish exactly, the DW mechanism will unavoidably affect the sterile neutrino DM population created by any other production mechanism. We present a semi-analytic approach to the DW mechanism acting on an arbitrary initial abundance of sterile neutrinos, allowing to combine DW with any other preceding production mechanism in a physical and precise way. While previous analyses usually assumed that the spectra produced by DW and another mechanism can simply be added, we use our semi-analytic results to discuss the validity of this assumption and to quantify its accurateness, thereby also scrutinising the DW spectrum and the derived mass bounds. We then map our results to the case of sterile neutrino DM from the decay of a real SM singlet coupled to the Higgs. Finally, we will investigate aspects of structure formation beyond the usual simple free-streaming estimates in order to judge on the effects of the DW modification on the sterile neutrino DM spectra generated by scalar decay.


Introduction
The biggest mystery of modern cosmology is the identity of Dark Matter (DM). This unknown substance amounts to ∼ 80% of the matter content in the Universe [1], and it is clear that it was important at the time when galaxies and other large structures have formed. Our best guess for the identity of DM is a new and electrically neutral particle, which is present in sufficient amounts in the Universe. Historically, the most plausible particle physics candidate was a Weakly Interacting Massive Particle (WIMP) [2,3]. However, up to now we have no unambiguous detection, neither in direct [4][5][6] or indirect [7][8][9] detection attempts nor by trying to directly produce them at colliders [10][11][12][13]. Instead, the limits are pushed towards smaller and smaller couplings.
On the DM side, N -body simulations using ordinary cold (i.e., non-relativistic) DM do not perfectly reproduce the Universe at small scales. While there is a rather obvious excess of dwarf galaxy numbers (historically referred to as the missing satellite problem [14,15], but also present in the field outside of galaxies [16]), more subtle discrepancies related to the internal kinematics of dwarfs, the cusp-core [17,18] and too-big-too-fail [19,20] problems, have been reported as well. Whether these problems can altogether be solved with a proper modelling of gas and stars in cosmological simulations or whether they hint towards a shift of the DM paradigm is still an open question (see, e.g., [21] for a review).
A generic non-cold DM candidate is a sterile neutrino with a mass of a few keVa heavier and much more feebly interacting version of the active neutrino. It can act as -1 -

JCAP04(2016)003
DM if produced in the early Universe in the right amount and with a suitable spectrum. The natural production for this DM is by freeze-in [22,23] (see also ref. [24] for a recent review) driven by the small active components of the predominantly sterile mass eigenstates. It was proposed by Dodelson and Widrow [25] in the context of DM ("DW mechanism"), based on refs. [26,27]. However, given the observational bounds on active-sterile mixing from not observing DM decay 1 as well as the relatively hot spectrum resulting from the DW mechanism, this simple way of producing sterile neutrino DM is excluded [52,53].
In this paper, we will re-investigate the DW mechanism. There are several reasons to do so. First, one often encounters incorrect statements about DW, e.g., that it would produce a spectrum of thermal shape [25,85]. This assumption is even used when deriving consequences for structure formation [86][87][88][89], and it also enters the exclusion itself [53,90]. Second, while not all DM can be produced by DW, the mechanism cannot be switched off as long as there is mixing between active and sterile neutrinos (which is generic, see e.g. refs. [91][92][93][94][95][96][97][98], or ref. [99] for a review). Thus, virtually any production mechanism will experience a modulation by the DW mechanism. This is the case for many settings investigated recently [69,71,75,79], yet the effect has not been taken into account. Third, we present the explicit example of scalar decay production being affected by the DW contribution. Decay production has been studied earlier by two of us (AM & MT) [70], so that the consequences of the DW modification do not only serve as an illustration for the current paper, but also as a-posteriori justification of our previous results. This paper is structured as follows. In section 2 we present a formal but general solution for the DM-distribution produced from or modulated by the DW mechanism, valid for any initial spectrum. We then use this formal solution in section 3 to derive an absolute upper bound on the strength of the DW-modification, which is strongly constrained by not observing X-rays from sterile neutrino decay. In section 4 we discuss quite generally the implications of 1 Summary of the status of the 3.5 keV line: evidence for this X-ray signal has initially been reported by [28,29] using XMM-Newton and Chandra data for several clusters as well as for Andromeda. Three main criticisms have been raised: 1.) The signal does not show up in all data sets [30][31][32][33][34][35] -although some groups found new versions of the line [33,36]. 2.) The emission lines of chemical elements are not treated correctly [37][38][39]. 3.) The signal does not follow the generic DM profiles at the centres of galaxies [33,40]. All these arguments have also been criticised [31,[41][42][43], and at the moment it is not clear what the final outcome will be. On top, also the original authors scrutinise their observations [44,45], which will hopefully clarify some aspects. In general, however, we should keep in mind that the "signal" is suffering from low statistics, which implies that both its "detection" and "exclusion" depend on the data sets used and analysis applied. Finally, future observations from Astro-H [46,47], LOFT [48,49], eROSITA [50], or searches using X-ray microcalorimeter sounding rockets [51] could soon clear up the situation [45]. 2 For parent particles with gauge charges, thermal corrections may be non-negligible [78].
the number of satellite galaxies on the maximal allowed fraction that DW produced sterile neutrinos can contribute to DM, irrespective of what makes up the rest. As a concrete example for a modulation of the DM distribution function by DW, we use section 5 to compute the effect on singlet scalar decay production. We will also exemplify some structure formation aspects. 3 We conclude in section 6. Technical details, such as a formal proof of the DW-consistency relation can be found in appendix A.
2 DW mechanism with non-zero initial abundance In this section, we present the most general solution to the Boltzmann equation describing the DW mechanism. Although the solution is formal, as it contains integrals which are not yet evaluated as they depend on the exact parameters and functions inserted, our solution allows to get an intuitive picture of the workings behind the mechanism.
In the generic setup used for DW production, the SM is extended by n R right-handed (RH) fermion singlets N i , where i = 1, . . . , n R . In fact, given that at least two light neutrinos are massive [100,101], we require at least two such RH neutrinos, n R ≥ 2 [102,103], or possibly even more such as the three fields considered in the neutrino-minimal SM (νMSM) [104]. However, when aiming at DW production only, we can work with n R = 1 for simplicity. This will not affect the results but allow us to concentrate on the actual mechanism. The inclined reader can find a more general discussion in appendix A.1.
With n R = 1, the following new terms can appear in the Lagrangian: We have introduced four new parameters: the Majorana mass M 1 and three Yukawa couplings y α with α ∈ {e, µ, τ }. Instead of the Yukawas, however, we will use the active-sterile mixing angles (θ e , θ µ , θ τ ): where v EW is the vacuum expectation value of the SM Higgs field. Note that eq. (2.2) is in principle only valid for small mixing, i.e. y α v EW /M 1 1. However, in practice this condition is always fulfilled due to the X-ray bound [105]. The three mixing angles (θ e , θ µ , θ τ ), although small, are the driving forces behind N 1 -production: their values decide about how often a process involving SM particles produces a sterile neutrino.

The full DW Boltzmann equation and its formal solution
The fundamental quantity describing the properties of a particle species in the Universe is its momentum distribution function f (p, t), abbreviated as MDF. It is governed by a set of Boltzmann equations. In the epoch we are interested in, well before the decoupling of the cosmic microwave background, the Universe is homogeneous and isotropic. Therefore, the MDF will only depend on the modulus p of the momentum and on time t or, equivalently, on the temperature T : f = f (p, T ). Note that, however, the temperature T denotes the temperature of the photons that stay in equilibrium until CMB decoupling.

JCAP04(2016)003
In a Friedman-Robertson-Walker Universe, the Boltzmann equation describing the dynamics of the sterile neutrino -and thus yielding the distribution f N (p, T ) -reads: 3) by the temperature-time derivative and inserting the DW collision term in an abstract form, we arrive at: where have defined the redshift integrand κ (T ) to be: In eq. (2.4), f th denotes a thermal Fermi-Dirac distribution and h, explicitly displayed in eq. (A.3), encodes all details of the conversion of active to sterile neutrinos. This latter quantity is where all the physics enters, and it will be discussed in detail in appendix A. It is also discussed at length in the literature, see e.g. [56,61]. Note that h = h[T, p, M 1 , (θ e , θ µ , θ τ )] also depends on the mass M 1 of the sterile neutrino and on all three mixing angles (θ e , θ µ , θ τ ). For the sake of clarity, we will however suppress these arguments whenever there is no risk of confusion. The compact notation allows to grasp the essential parts of the Boltzmann equation for the sterile neutrino. Since the active-to-sterile conversion is a 1 ↔ 1 process, it is clear that the Boltzmann equation must depend linearly on both the MDFs of the sterile and active neutrinos (the latter being given by f th as long as active neutrinos are in thermal equilibrium). The right-hand side of eq. (2.4) can be interpreted as the sum of a gain term: h (T, p) f th (T, p) , and a loss term: − h (T, p) f N (T, p) . (2.6) In the standard DW scenario without any initial abundance, the loss term is usually neglected, 5 which greatly simplifies the computation. We will keep the term, though, thus 4 In principle, also scattering channels have to be accounted for in the collision terms. While not changing the number density of a species, scatterings may change the distribution of the momentum modulus p. In cases where eq. (2.3) will finally be integrated to obtain a Boltzmann equation on the level of particle number densities, scatterings are usually neglected, though they can have an effect in theory. 5 It can, however, not be neglected in the case of resonant active-sterile conversion, even for vanishing initial abundance, see ref. [56]. allowing for an arbitrary distribution function of sterile neutrinos produced from another production mechanism before the onset of DW, which becomes efficient only at temperatures of O (100 MeV) [56].
Before advancing to the full solution of eq. (2.4), let us investigate the properties of κ (T ). Conservation of the comoving entropy density can be cast as g S (T ) T 3 a 3 (T ) = const., where g S (T ) is the number of effective entropy degrees of freedom (d.o.f.). Differentiating this equation with respect to time and changing variables from t to T yields: where g S denotes the derivative of g S with respect to T . By virtue of eq. (2.7), we can immediately derive: which will be used to compactify our final solution. The above result justifies the name "redshift integrand" for κ: in a collisionless Boltzmann equation, the solution of eq. (2.3) has to account for redshift only, and for an initial distribution f ini (p) it reads 6 indeed fulfilling the boundary condition f collisionless (p, T ini ) = f ini (p). The exponentiated integral of κ turns out to be precisely the factor that accounts for the redshift of a collisionless species. If the number of d.o.f. stays constant, the expected approximate redshift proportionality to T −1 is recovered, too.
Let us now proceed to the solution of eq. (2.4). In its most condensed form, the solution at some final temperature T f reads where we have made use of the following abbreviations: In eq. (2.12), we have introduced the notation (hf th ) (T, p) ≡ h (T, p) f th (T, p). As before, we have suppressed the arguments M 1 and (θ e , θ µ , θ τ ). The prefactor S can be interepreted as damping factor converting part of the initial distribution f ini into active neutrinos by activesterile conversion in addition to redshifting the distribution. The product of Sf DW can be interpreted as the pure DW contribution, resulting from a vanishing initial abundance.

JCAP04(2016)003
Consistency check. The form of eq. (2.10) is most suitable to discuss an important consistency check of the solution: the MDF at some temperature T 1 has no memory of the dynamics that shaped it at temperatures T > T 1 . Accordingly, the temperature T ini can be chosen arbitrarily if the DW contribution produced at T > T ini is included into f ini . To be more concrete, we expect the following relation to hold: . This equation states that sterile neutrinos produced via DW until some temperature T 3 can just be re-interpreted as initial abundance f ini present at T 3 ; in fact, they could have been produced by any mechanism. But, "stopping the clock" at T 3 and re-starting it, they can be used as a starting condition for DW production from temperature T 3 onwards. While the physics of this relation should be intuitively clear, we will present a formal-analytic proof in appendix A.3.

A note on the DW case without initial abundance
Though being ruled out as sole production mechanism [53], we have argued before that DW can contribute subdominantly to the relic abundance, e.g. in a mixed DM setting. Some features attributed to pure DW production in the literature are in fact not correct. In particular the close-to-thermal shape noted in [25,85] has been used to exclude the DW mechanism in the first place. But, while the DW mechanism is in any case excluded as only source of sterile neutrino DM, the current bounds can actually be enhanced if one takes into account the correct spectral shape. It is worthwhile to use our semi-analytical results to discuss the case of DW production without initial abundance in a precise and physical way, in order to probe the quality of the approximations used elsewhere. We will in particular discuss why a suppressed thermal shape of the DW spectrum, as often adopted in the literature, is in fact not a very accurate estimate, especially if the high momentum part of the distribution is important -like in analyses concerning cosmological structure formation where precisely that part puts the DW mechanism into trouble.
To do so, let us solve eq. (2.4) neglecting the term −h (T, p) f N (T, p) on the right-hand side. Then, the solution at temperature T f , as derived from eq. (2.10), reads: Since a thermal distribution of a nearly massless species just depends on the ratio p/T of momentum and temperature, f th in eq. (2.14) depends on T 2 only via the term g S (T 2 ). Thus, only if g S varied sufficiently slowly with T 2 , one could replace g S (T 2 ) by some average value g S and pull the thermal part f th in front of the integral, hence resembling the shape of a thermal distribution (only multiplied by a suppression factor). If one can do that, the solution to eq. (2.4) is indeed given by: Even if that was the case, which is however certainly not true at least during the QCD transition as it generically happens around the peak of DW production [25,26], for the resulting distribution to be of thermal shape the function h would in addition need to vary only very slowly with the momentum p. But, plotting h as a function of p for different values of T 2 , cf. left panel of figure 1, we see that this is not at all the case. Accordingly, the statement about the distribution being of thermal shape (only redshifted and with a suppression factor) is just not correct.
To quantify the discrepancy between the approximation and the numerical result, we have chosen three benchmark cases: M 1 = 3 keV, M 1 = 7 keV, and M 1 = 25 keV. We have then computed the ratios between the approximate result in eq. (2.15) and the exact one in eq. (2.14) for pure e-mixing (the results for pure µ-or τ -mixing are very similar). The freedom of choosing g S was eliminated by fixing this parameter such that the particle number density is equal to the numerical result, which -of course -is unknown in case one uses the approximation only! In this sense, our ex-post choice of g S yields the best approximation possible, and even this is not very accurate for high momenta. The right panel of figure 1 illustrates the deviation for the three benchmark cases as a function of rescaled momentum x = p/T . Note that these benchmark cases are valid for any choice of the mixing angle, since by virtue of eqs. (2.14) and (2.15) both the approximation and the numerical result with vanishing initial abundance are exactly proportional to sin 2 (2θ).
The right panel of figure 1 reveals that the approximation is perfect for x → 0, which is due to the fact that the Fermi-Dirac distribution tends to 1/2, irrespective of the choice DW parameter space (before Ly-α) X-ray (hyp) Tremaine-Gunn X-ray (Suzaku) Tremaine-Gunn X-ray (Suzaku) Parameter space for pure DW production before applying limits from structure formation with isoabundance lines for numerical limits and approximative results.
of g S . However, the thermal shape approximation systematically underestimates the high momentum modes, which are in fact the most decisive ones when excluding the DW mechanism by cosmic structure formation. We want to emphasise once more that the best choice of g S is a priori unclear. To give an impression of the significance of different meaningful choices, we show in figure 2 numerical and estimated isoabundance lines in the plane spanned by M 1 and sin 2 (2θ) for the cases of e-, µ-, and τ -mixing. The blue curves represent the contour where a pure DW production yields the correct abundance if calculated numerically, while the magenta lines use two different a priori choices of g S , namely as arithmetic mean and (2.16) (2.17) The figure contains limits from X-ray observations (for a detailed explanation of the most conservative bound dubbed hyp, see section 3) as well as the Tremaine-Gunn bound. In all three cases, using a meaningful average g S can lead to an overestimate of the square of the mixing angle by about half an order of magnitude when fixing the abundance to the current best-fit value from Planck [1]. To complete this discussion, we show the numerical distribution function as compared to the estimated ones in figure 3. We have chosen a mass of M 1 = 2 keV and pure e-mixing since, according to figure 2a, this is about the maximum mass that can reproduce the observed relic abundance without violating the most conservative hyp X-ray limit. Of course we anticipate that all these spectra will not be in agreement with bounds from Ly-α observations. Nonetheless, we will show the effect of only estimating the distribution for the sake of completeness. Before advancing to more sophisticated analysis methods in section 4, let us present a simple way to estimate the compatability of DW production and structure formation. We will derive a translation between the mass of a non-thermally produced DM particle with mass M 1 and a thermal relic with mass m w , such that the implications for structure formation are roughly the same. Such a translation can already be found in [106, eq. (12)]. While we could  not reproduce the exact numerical coefficients of [106, eq. (12)], because the assumptions made to obtain it are not given explicitly in that paper, we present the most relevant steps of the derivation in a fully parametric way, i.e. keeping all model dependent parameters in the final result.
The rationale behind our translation is the idea that the average velocity of the DM particles at a certain time after production is a good indicator for structure formation. Let us thus first analyse a thermally produced species of mass m w and g w internal d.o.f. that has a temperature T w (which may deviate from the plasma temperature T = T γ ). Let us assume that this species decouples at a temperature where the background plasma counts g dec entropy d.o.f.
Calculating the relic abundance for this species, we find the relation: for a fermionic species, 1 for a bosonic species, (2.18) where h = H 0 / 100 kms −1 Mpc −1 , ρ crit is the critical density of the Universe, and s 0 is today's entropy density. For a thermal species in the mass range of some keV to MeV, we can safely assume that it has become non-relativistic at matter-radiation equaility, such that the average velocity at this epoch is given by:

JCAP04(2016)003
In contrast to the species that we have discussed now, we do not assume that the sterile neutrino has any thermal history, i.e., we assume no spectral form for the momentum distribution function. We only assume that the production ceases at some temperature T prod when there are g prod entropy d.o.f. in the plasma. We parametrise the average momentum p at production temperature by p = αT prod , where α is usually of O (1). Then, at equality, the average velocity of the sterile neutrino with mass M 1 -again assuming it is non-relativistic -is given by: Equating (2.20) and (2.21), we can solve for M 1 : for a fermionic species. Choosing suitable values for the d.o.f. and setting α to a value in the region 2 − 3, as appropriate for DW production, we find a result that is in the same ballpark as [106, eq. (12)], but tends to yield a slightly smaller sterile neutrino mass for a fixed thermal mass m w . Also the recent conversion formula between thermal masses and sterile neutrino masses in [107] yields slightly smaller masses than [106, eq. (12)]. This tendency shows that the conversion formula in [106, eq. (12)] is slightly too agressive, resulting in somewhat too tight limits as we will show in section 4, together with the analysis of the approximations for a pure DW spectrum presented in eqs. (2.15)-(2.17). On the other hand, we will also show in this section that using the exact spectrum could even improve the limits.

Quantitative analysis of the mechanism
Having seen the formal-analytical solution of eq. (2.4), let us have a look at the numerical properties of this solution, using numerical fits for g S presented in [108]. We start by discussing whether or not it is a good approximation to just add the DW-component by hand to any (correctly redshifted) initial distribution of sterile neutrinos, i.e., to neglect potential damping or distortion effects of the active-sterile conversion on the spectrum. This has been done in several references [69,71,75,79], however, in none of those works the reliability of this assumption has actually been checked. Given our general solution, we will now show that the approximation applied in the earlier refereces does indeed work. We thus provide an a-posteriori justification of the treatments used in previous works.
Looking again at the formal solution reported in eq. (2.10), we realise that this question about the reliability of the approximation described above can easily be answered by analysing how much the factor S (T f , T ini , T f , p) can deviate from unity. To this end, we plot the contours of S in the plane spanned by x f = p/T f and M 1 , both for a (quite large) example value of 5 · 10 −5 for the different mixing angles θ e,µ,τ (figure 4) and for the maximal mixing    angle allowed by X-ray observations dubbed θ = θ max (M 1 ), see figure 5. We chose T ini = 10 GeV and T f = 10 MeV, thereby spanning the range relevant for DW production [56]. Extending this temperature range does not alter the reults. The numerical integration has been performed using the CUBA library [109].
In the case of θ = θ max , we show the results of the analysis for two different versions of the X-ray bounds. The top panel uses an analysis based on data from Suzaku [110] (dubbed Suzaku), which update the combined limits obtained in [111], while the bottom panel relaxes the limits on θ 2 max by a factor of 5 (dubbed hyp). The latter option illustrates what would change if our current face value bounds were in fact too strong, which is not completely unrealistic given the intrinsic uncertainties of X-ray observations. 7 We further indicate by a dashed line the mass below which the maximally allowed mixing angles overproduces sterile neutrinos from DW alone. This means that all masses below this line are excluded, irrespective of the initial distribution before onset of DW production.
As evident from figure 5, the maximal suppression of any momentum mode in the initial spectrum can be a few percent at most. No matter which scenario is taken, it is either the X-ray bound itself which forces the mixing angles to be small and thus implies a small effect (for sterile neutrino masses of roughly 4 keV and larger), or it is the constraint of not producing too much DM which does not allow for large active-sterile mixing (below a mass of, say, 2 keV). Only around sterile neutrino masses of M 1 ∼ 3 keV the DW modification could possibly be significant. Even in that region, however, the effect turns out to be very marginal, and it never amounts to more than 5%.
This shows that simply adding a DW component to any previously produced sterile neutrino population is indeed a pretty good approximation. No big modifications of the spectrum are expected. While the DW modification turns out to be negligible in practice, this result was not a priori clear beforehand. However, one can see it from our equations, as the small values of h force the factor S to be close to unity.

Constraints from structure formation
It is well known that structure formation imposes stringent constraints on the minimal sterile neutrino mass produced by the DW mechanism. In combination with the bounds from X-ray data, they rule out DW sterile neutrinos as the dominant DM component [52,112]. Different sterile neutrino production mechanisms may still have subdominant DW contributions, though, as explained in the sections above. It is therefore useful to identify the maximum allowed contribution of DW sterile neutrinos with respect to the total DM budget.
To constrain the DW sterile neutrino abundance, we assume the limiting case of a mixed DM model where the remaining DM component is considered to be perfectly cold. The structure formation of such a scenario has been investigated in various other studies (see e.g. [113,114]) and there are constraints based on both Lyman-α forest [115] and on dwarf galaxy counts [116]. In practice, the remaining DM component does not have to be perfectly cold (as for example in the case of resonantly produced sterile neutrino DM where the colder component still consists of luke-warm DM [115]), which would lead to even stronger constraints on the abundance of the DW component.

JCAP04(2016)003
In this paper we apply the method presented in ref. [116], which consists of comparing the expected number of sub-haloes to the observed abundance of satellite galaxies in the Milky-Way. A good estimate of the sub-halo abundance is given by the relation: where variances S i and masses M i are defined as with sh and hh standing for subhalo and host-halo, respectively. Eq. (4.1) is based on a modification of the extended Press-Schechter recipe (using a sharp-k window function) described in [117] (see also [118]) and a normalisation to N -body simulations of [89] to account for tidal stripping. The only input it requires is the linear power spectrum P (k) for a given cosmology and DM model. For the estimate of Milky-Way satellite numbers, we follow the method of [119] (see also [120]) which consists of adding up 11 classical and 52 ultra-faint dwarf galaxies (the latter is an estimate based on the 15 ultra-faint dwarfs observed by SDSS multiplied by a factor of 3.5 to account for the limited SDSS sky coverage). We refer to the refs. [116,119,120] for details about this method including discussions about statistical and systematical uncertainties.
Based on stellar kinematics it is possible to estimate the total mass of the observed Milky-Way satellites, albeit with rather large systematic uncertainties. We therefore only assume that all known satellites have total halo masses above 10 8 h −1 M , which is a conservative estimate based on results from ref. [121]. For the mass of the Milky-Way, we consider the range 5.5 × 10 11 < M hh < 3.2 × 10 12 h −1 M given by ref. [122].
The bounds on the maximally allowed abundance of DW sterile neutrinos are shown in figure 6. Assuming an average Milky-Way mass of 1.2 × 10 12 M /h (pink exclusion area) the allowed fraction of DW sterile neutrinos (i.e. f = Ω DW /Ω DM ) never exceeds f ∼ 0.3. It is strongly suppressed beyond M 1 = 3 keV due to the X-ray limit from Suzaku. The allowed fraction increases to a maximum of f = 0.5 if a very large Milky-Way mass of 3.2×10 12 M /h is assumed instead (red exclusion area). Note that our analysis rules out a sizable fraction of parameter space that was found to be valid in a similar analysis done in refs. [75,123,124]. Especially in the phenomenologically interesting region around 7 keV, the maximal fraction of sterile neutrinos produced from the DW mechanism in ref. [75] turns out to be far too optimistic and is not in agreement with our results using the new Suzaku data, as to be expected from the tight constraints on the mixing angle for such masses.
As announced in section 2.2, we complete this session by showing the relative power spectra computed from the different versions of the MDF in figure 7. The blue solid line (to the very left) describes the full numerical solution, while the red and the green lines correspond to the different choices of g S in the approximation of eq. (2.15). The dotted black line is the approximation used in [86]. As already discussed, in the case of the two approximations, the relative power spectra fall off only at higher wave-numbers as compared to the numerical case and are therefore too conservative. On the other hand, the approximation used in ref. [86] is slightly too restrictive, showing that the correct treatment of a DW component in a mixed model for sterile neutrino DM might reallow former borderline cases.

Example: initial abundance from scalar decay
Having laid the technical foundations for the computation of DW production combined with an initial abundance, let us now proceed to a well-motivated mechanism that could produce such an initial spectrum: one real scalar singlet S, coupled to the Higgs via some portal coupling λ, and one sterile neutrino N 1 , coupled to the scalar via a Yukawa-type interaction with coupling strength y. This setting was extensively discussed in the literature (see e.g. [67,68,125]). A detailed analysis on the level of distribution functions was presented in an earlier work by two of us [70], however, with the assumption of zero active-sterile mixing. We now use the opportunity to kill two birds with one stone: the shortcoming of the previous paper will be rectified simultaneously to delivering an explicit example for the current paper.
In order not to unnecessarily prolong this text, we will only sketch the most important aspects of scalar decay production. We refer to ref. [70] for details.

Singlet scalar decay production of sterile neutrinos
Our setup extends the SM by two gauge singlets, one real scalar S and one sterile neutrino N 1 . The Lagrangian is given by:  Figure 7. Relative power spectra for sterile neutrinos of a mass of M 1 = 2 keV being produced by the DW mechanism exclusively. Using the exact numerical version of the MDF will yield the tightest constraints on a model with a non-nigligible fraction of sterile neutrinos produced via the DW mechanism. The light-and dark-grey shaded areas illustrate the 2-σ exclusion limits from Lyman-α observations [53] and dwarf galaxy counts [119].
where L SM is the SM Lagrangian, L ν is the part of the Lagrangian that gives mass to the active neutrinos [and contains the mixing angles (θ e , θ µ , θ τ )], and V scalar is the scalar potential containing a Higgs portal coupling ∝ λ H † H S 2 between the SM-like Higgs H and the new scalar. With this Lagrangian, two basic processes are relevant for the production of sterile neutrinos: As in [70], we will apply the assumption that the scalar decay occurs well before the QCD transition, such that g S ∼ O(100) can be taken to be roughly constant during DM production. Denoting the scalar mass as m S , we can express our results as functions of two effective parameters [70], the effective decay width:  • The effective Higgs portal C HP determines whether or not the scalar equilibrates in the early Universe. If C HP is small enough, the scalar undergoes a freeze-in process [68] before decaying into sterile neutrinos, while it instead freezes out if C HP is large enough [67]. This translates directly into different sterile neutrino MDFs.  Table 1. Overview of the benchmark cases. All parameter sets are chosen such that the relic abundance of sterile neutrinos is in accordance with the best-fit value [1]. Note that we do not assume maximal mixing in agreement with X-ray constraints in case 3β since this would violate the overabundance bound. Still, the parameter C Γ in case 3γ is chosen to reproduce the observed relic abundance with the same mass M 1 as in case 3β.
• The effective decay width C Γ determines how fast the scalar decays. It thus dictates when the energy stored in the mass of the scalar is injected into the Universe in the form of highly relativistic sterile neutrinos (recall that M 1 m S ). Therefore it ultimately decides about the spectrum being warmer or cooler, and it is very relevant to determine the consequences for structure formation. Also, if the scalar equilibrates, there is an upper limit for C Γ beyond which the particle number of sterile neutrinos becomes so large that any mass not violating the Tremaine-Gunn bound [126] will always overclose the Universe, see ref. [70] for details.

Benchmark cases with and without a DW component
With the preceeding discussion we can motivate our choice of three benchmark cases that will be taken as initial abundance for the DW production: figure 5 suggests that the sterile neutrino mass M 1 should not exceed a few keV in order not to violate the overabundance bound, while simultaneously not having very small mixing angles due to X-ray constraints. First estimates on structure formation simultaneously demand a very early production for sterile neutrinos in this mass range (see [70, figure 8]). Recall that a very early production is equivalent to a large C Γ , which restricts us to the freeze-in region of the scalar decay parameter space. The exact parameters defining these cases will be summarised in table 1. Each case will in turn be subdivided into three subcases. The first subcase -dubbed α -corresponds to a fixed choice of (C HP , C Γ , (θ e , θ µ , θ τ ) = (0, 0, 0)), i.e., DW production switched off completely.
Since we assumed M 1 m S , the mass of the sterile neutrino is completely irrelevant when calculating the distribution function f N (p, T ) from scalar decay. Combining this fact with the parametrisation from eq. (5.2), it is clear that the distribution from scalar decay will solely depend on our effective coupling parameters C HP and C Γ . The mass scales M 1 and m S only enter when calculating dimensionful quantities like the free-streaming horizon or the relic density. Furthermore, the discussion of section 3 shows that the suppression of the initial abundance, which does depend on M 1 , cf. eq. (2.10), is very weak -such that the -16 -

JCAP04(2016)003
Subcase Value of C HP Value of C Γ θ Mass α C HP C Γ θ = 0 matched to Planck data β C HP C Γ θ = 0 matched to Planck data γ C HP C Γ θ = 0 same value as in β relic abundance of the combined production mechanism can be cast into the following form, to a very good approximation: Here, dΩ SD dM 1 is the relic abundance per sterile neutrino mass for the case of pure scalar decay (cf. [70, figure 3]).
If we choose, for instance, (θ e , θ µ , θ τ ) = (θ max (M 1 ) , 0, 0) in eq. (5.3), we can numerically solve for M 1 for a fixed set of (C HP , C Γ ), in such a way that Ω SD+DW (M 1 ) matches the current best-fit values from Planck [1]. These are the subcases β, corresponding to scenarios where the correct abundance is achieved for the same couplings (C HP , C Γ ) as in subcase α, even when including the DW contribution, however paying the price of having to rescale the sterile neutrino mass. 8 This case is useful whenever the parameter space is scanned in terms of the effective couplings, e.g., when aiming at producing new versions of figures 3, 8(a), or 8(b) in ref. [70].
The value for the mass obtained via this procedure can then be mapped back to a case with (C HP , C Γ ), such that SD without DW yields the same mass when applying the relic abundance constraint (subcases γ). Keeping C Γ constant guarantees that the production time of both scenarios is comparable, which is the most meaningful comparison when it comes to structure formation. The defining characteristics of subcases α, β, and γ are briefly summarised in table 2.
In figure 8, we present the distribution functions for all (sub-)cases. In each panel, the solid blue curve corresponds to subcase α, while the dashed blue curve corresponds to subcase β and the green one to subcase γ. We also quote the average momenta x of the distribution. Note that the areas under the curves of subcase β and under the ones for subcase γ are identical within one and the same case by definition.
In order to assess the compatability of the distribution functions shown in figure 8 beyond the simple approach using the average quantity of the free-streaming horizon, we have computed the linear power spectra P for all 9 (sub-) cases and normalised them to the linear power spectra of a CDM analogue (P CDM ). The results are shown in figure 9, where we have adopted the colour coding of figure 8. The figure clearly shows that, in all three benchmark cases, the spectrum from pure SD is always colder than the one with a DW component added. For masses around 7 keV (left panel) the difference is, however, very small, and all subcases are in agreement with the bounds from structure formation. For intermediate masses (central panel), the combined spectrum is only borderline compatible with structure formation, while the pure SD cases lie just a bit beyond the exclusion region. For smaller masses, around 1 keV, the difference between the different cases is largest, as expected from the fact that JCAP04(2016)003  . Relative power spectra of the benchmark-cases from table 1. The light-and dark-grey shaded areas illustrate the 2σ exclusion limits from Lyman-α observations [53] and dwarf galaxy counts [119]. The black dotted curves illustrate the pure DW case with the lower of the two particle masses given in table 1.
X-ray bounds are very weak in this regime. However, the corresponding subcases of case 3 (right panel) are clearly not compatible with Lyman-α observations. This benchmark analysis shows that the additional DW component for a SD spectrum can usually be neglected, either because the effect is tiny to start with or because the model will be excluded with or without the additional DW component. Only in a very narrow range of parameters the additional DW component can possibly affect the validity of formerly borderline cases. 9 In the phenomenologically interesting region around 7 keV, it is completely irrelevant whether the DW component is taken into account or not.

JCAP04(2016)003 6 Conclusions
In this work, we have presented for the first time a fully comprehensive semi-analytical method for calculating the effect of the DW mechanism on any initial distribution of sterile neutrinos. Treating the combination of DW and any other production mechanism in the most abstract way has allowed us to find a formalism where the interference effects between the initial spectrum and the DW effect can be put into an intuitive form. The formalism yields an analytical expression of the combined momentum distribution function for every point in time (or temperature), the evaluation of which is however only feasable numerically.
Taking into account constraints from X-ray observations, we have used our method to show that the spectra of the preceeding production mechanism and DW production can simply be added (after being correctly redshifted) to a very good approximation: interference effects are of the order of a few percent at most. We also used our approach to assess the quality of the common approximation of a suppressed thermal shape of the DW spectrum. We have seen that the estimates with suppressed thermal shape are quite sensitive to the definition of the average number of degrees of freedom during production. Two meaningful choices of this average clearly underestimate the abundance for a given mass and mixing angle. They also underestimate the weight of high momenta in the sterile neutrino spectrum, which implies that the limits previously obtained from structure formation considerations are in fact too weak.
Furthermore, we have set limits on the maximal fraction of sterile neutrino DM produced from the DW mechanism using our numerical results and astrophysical constraints. In a conservative analysis, DW can contribute about 50% of the total DM abundance, however, only for a mass of about 2 keV. This results sets tight constraints on earlier and possibly too optimistic interpretations of the data.
We finally applied our method to the example case of an initial sterile neutrino abundance produced from the decay of a scalar singlet particle. Since a pure scalar decay spectrum of sterile neutrinos is not perfectly cold, the additional DW contribution must in any case be smaller than the 50% quoted above. In this scenario, we find that the maximal contributions of DW to the initial spectrum lie in a negligible parameter range, where the resulting combined spectrum would be excluded by structure formation anyway. In regions where the remaining combined spectrum does not violate structure formation bounds, the effect of DW is negligibly small, which justifies the assumption of neglecting the DW contribution for scalar decay production as made in [70]. In the phenomenologically interesting region of a sterile neutrino of M 1 = 7.1 keV we find that it is indeed irrelevant whether DW is accounted for or not if the scalar decays fast enough not to violate bounds from structure formation. Comparing the linear transfer function of this case to current bounds from Lyman-α observations, we showed that scalar decay is indeed a viable production mechanism for keV sterile neutrinos.

JCAP04(2016)003
A Technical details on the formal solution of eq. (4) A. 1 On the number of sterile neutrinos n R In the general case with arbitrary n R , eq. (2.1) takes the form where the Majorana mass matrix M ij must be symmetric to fulfill the Majorana condition and can therefore always be diagonalised unitarily. There is no theoretical upper limit on n R . If active-neutrino masses are to be explained by a seesaw mechanism, n R ≥ 2 is needed to accommodate the two distinct mass scales ∆m 2 atm and ∆m 2 sol observed by neutrino osciallation experiments [102,103].
Before discussing the physical aspects of DW, it is worthwhile to work out the number of parameters in a general 3 + n R framework: choosing a basis where M ij is diagonal, we are left with the n R real eigenvalues (if CP is conserved) and 6n R real parameters in the complex Yukawa matrix y αi , which can be subdivided into 3n R mixing angles and equally many phases. Three of these phases can be absorbed into the charged lepton fields and hence we are left with n R Majorana masses, 3n R mixing angles and (3n R − 3) phases, i.e. (7n R − 3) free parameters in total. In this case, the mixing angles are given by again assuming them to be sufficiently small.

A.2 The DW mechanism
To recall the Dodelson-Widrow (DW) mechanism (see [25] or also [26,27] for earlier works) in a nutshell, it is basically some freeze-in [22,23] type of DM-production. This means that the DM particles are interacting so feebly that they never enter thermal equilibrium. Instead, they start of with a certain zero or non-zero abundance and are gradually produced in the early Universe as long as they are kinematically accessible. In the case of sterile neutrinos, their small admixtures θ α with the active-neutrino sector cause them to be produced in the small fraction |θ α | 2 of reactions where a vertex of flavour α happens to produce the keV-scale mass eigenstate N 1 instead of one of the three light mass eigenstates ν 1,2,3 . The Boltzmann equation describing this type of production has been given in eq. (2.4), and the part of the right-hand side which we abbreviated as h(p, T ) is explicitly given by 10 h(p, T ) = 1 8 Γ α (p, T )∆ 2 (p) sin 2 (2θ α ) ∆ 2 (p) sin 2 (2θ α ) + [Γ α (p, T )/2] 2 (p) + [∆(p) cos(2θ α ) − V α (p, T )] 2 , where Γ α (p, T ) are the interaction rates of active (anti-) neutrinos of flavour α, V α (p, T ) is the background potential for active (anti-) neutrinos of flavour α, ∆(p) = (M 2 1 − m 2 ν )/(2p) M 2 1 /(2p) (where m ν denotes any light neutrino mass which is always negligible compared to M 1 ), and θ α is the active-sterile mixing angle for neutrinos of flavour α, i.e., the fraction of the sterile neutrino mass eigenstate contained in the active flavour α. The basic difficulty is to compute the interaction rates Γ α (p, T ) accurately, and to have a reliable expression for the potential V α (p, T ). 10 In the case of a primoridal lepton number asymmetry there would be further potentials [56]. A.2.1 The interaction rates Γ α (p, T ) The basic form of the interaction rate Γ α (p, T ) is given by [56,127,128]: where G F = 1.166 · 10 −5 GeV −2 is the Fermi constant, p is the sterile neutrino momentum, and C α (T ) are temperature-dependent functions. These functions depend on the details of the dynamics in the plasma in the early Universe, and while first computations have been presented quite a while ago [129], a detailed 2-loop calculation was performed only later on [57]. Most importantly, the results obtained in this reference have been recently made available in machine-readable numerical data files. 11 The values of C α (T ) as functions of the temperature are depicted in figure 10. Note that the functions as displayed here include the QCD contributions which are neglected in some references, like in ref. [128]. In particular this information complements the treatment presented in ref. [56], where the interaction rates were only presented for a relatively narrow temperature range.
A.2.2 The potentials V α (p, T ) Also the potentials V α (p, T ) are not often discussed in the generality required and, partially, the literature is plagued by unfortunate typos. We thus display the potentials, although known in principle, in full generality (see the discussions in [56,131]): where the upper (lower) sign holds for neutrinos (anti-neutrinos), ζ(x) is the Riemann ζfunction, and η B = 6.05 · 10 −10 is the baryon asymmetry. 12 Here, the number densities and In the case of f DW , we can only find an approximate relation: f DW T a , T b , p, ∆m 2 , (θ e , θ µ , θ τ ) ≈ i=e,µ,τ This approximation is valid as long as the term converting sterile neutrinos into active ones is small, which is certainly the case for pure (i.e. non-resonant) DW production without initial abundance. In this case, we get the simple scaling behaviour expected for the solution, being directly proportional to the sine-square of the active-sterile mixing angle.