Improved determination of sterile neutrino dark matter spectrum

The putative recent indication of an unidentified 3.55 keV X-ray line in certain astrophysical sources is taken as a motivation for an improved theoretical computation of the cosmological abundance of 7.1 keV sterile neutrinos. If the line is interpreted as resulting from the decay of Warm Dark Matter, the mass and mixing angle of the sterile neutrino are known. Our computation then permits for a determination of the lepton asymmetry that is needed for producing the correct abundance via the Shi-Fuller mechanism, as well as for an estimate of the non-equilibrium spectrum of the sterile neutrinos. The latter plays a role in structure formation simulations. Results are presented for different flavour structures of the neutrino Yukawa couplings and for different types of pre-existing lepton asymmetries, accounting properly for the charge neutrality of the plasma and incorporating approximately hadronic contributions.


Introduction
Despite considerable phenomenological success, the Standard Model of particle physics suffers from a number of shortcomings: it can explain neither the presence of neutrino oscillations, nor of dark matter, nor of a cosmological baryon asymmetry.It is remarkable that all of these problems can in principle be cured by a simple enlargement of the theory through three generations of right-handed neutrinos [1]- [4], without changing the underlying theoretical principles of gauge invariance and renormalizability.Despite its simplicity it remains unclear, of course, whether nature makes use of this possibility.
In the present paper we are concerned with the dark matter aspect (for reviews see, e.g., refs.[5,6]).For the parameter values that are phenomenologically viable, the dark matter sterile neutrinos do not contribute to the two observed active neutrino mass differences, so that the dark matter aspect is partly decoupled from the neutrino oscillation and baryon asymmetry aspects.The decoupling is not complete, however: it turns out that the sterile neutrino dark matter scenario is tightly constrained [7], and only works if the dynamics responsible for generating a baryon asymmetry also generates a lepton asymmetry much larger than the baryon asymmetry [8], which subsequently boosts sterile neutrino dark matter production through an efficient resonant mechanism first proposed by Shi and Fuller [9].Despite the tight constraints, indications of a possible observation [10,11] demand us to take this scenario seriously.For our purposes, the complicated dynamics of ref. [8] simply amounts to the fact that the initial state at a temperature of a few GeV possesses certain lepton asymmetries.
The goal of the present paper is to refine the analysis of ref. [7] in a number of ways, both theoretically and as far as the numerical solution of the rate equations is concerned.The resulting spectra could be used as starting points in structure formation simulations (cf.e.g.refs.[12]- [14] and references therein).
Our plan is the following.In sec. 2 we review and refine the theoretical description of sterile neutrino dark matter production from a thermalized Standard Model plasma with preexisting lepton asymmetries.The practical implementation of the corresponding equations and a strategy for their solution are discussed in sec.3. Numerical solutions are presented in sec.4, and we conclude with a discussion in sec. 5.

General derivation of the rate equations
Our first goal is to derive a closed set of rate equations for the sterile neutrino distribution function and for lepton number densities, 1 valid both near and far from equilibrium.In ref. [7] such equations were postulated but the argument involved phenomenological input, in order to account for two types of "back reaction".Our derivation yields an outcome that differs slightly from ref. [7].Moreover, the influence of electric charge neutrality of the plasma on the relation between lepton number densities and lepton chemical potentials was not properly accounted for in ref. [7].Finally, in order to obtain a maximal effect, the lepton asymmetries of the different flavours were treated as equilibrated in ref. [7], even though flavour equilibrium (through active neutrino oscillations) is expected to be reached only at temperatures below 10 MeV or so [15,16].We eliminate all of these simplifications in the present paper.
We start by defining, in sec.2.1, the quantities appearing in the equations.In sec.2.2 an evolution equation is obtained for the sterile neutrino distribution function.In sec.2.3 the same is achieved for lepton number densities, with a right-hand side expressed in terms of lepton chemical potentials and sterile neutrino distribution functions.The system is completed in sec.2.4, where we relate lepton chemical potentials to lepton densities, taking into account electric charge neutrality of the Standard Model plasma.

Definitions
We consider a system consisting of right-handed (sterile) neutrinos and Standard Model (SM) particles.The SM particles are assumed to be in thermal equilibrium at a temperature T .The initial state is characterized by non-zero lepton chemical potentials, µ a , associated with different lepton flavours.(At low temperatures T < ∼ 10 MeV the lepton flavours equilibrate through active neutrino oscillations, so that there is only a single lepton chemical potential, denoted by µ L .)The initial state may also contain an ensemble of sterile neutrinos.The two sectors communicate through Yukawa interactions, parametrized by neutrino Yukawa couplings.As a result, the distribution function of the sterile neutrinos evolves towards its equilibrium form, and the lepton densities decrease, because lepton number is violated in the presence of neutrino Yukawa interactions and Majorana masses.(However, in practice neither process gets completed within the lifetime of the Universe.) Let us denote by h the matrix of neutrino Yukawa couplings.We work out a set of rate equations to O(h 2 ).This means that, as soon as we have factorized a coefficient of O(h 2 ), the sterile neutrinos can be treated as mass eigenstates and free particles, with masses given by Majorana masses.Sterile neutrinos can then be represented by on-shell field operators, where The index I enumerates the sterile neutrino "flavours", and M I is their mass, which we assume to be real and positive.The creation and annihilation operators satisfy the anticommutation relations {â Ikσ , â † The on-shell spinors u, v are normalized in a usual way, for instance The Majorana nature of the spinors requires that u = C vT , v = C ūT , where C is the charge conjugation matrix.Now, let us define the ensemble occupied by the sterile neutrinos.Their distribution is described by a density matrix, denoted by ρN .We take a rather general ansatz for ρN , assumed however to be diagonal in flavour, momentum, and spin indices: where the function µ Ik is left unspecified except for being spin-independent, Z N takes care of overall normalization, and k ≡ d 3 k/(2π) 3 .The density matrix has also been assumed to be x-independent.Note that even though the function µ Ik bears some resemblance to a chemical potential, it is not identical to one (Majorana fermions are their own antiparticles and no chemical potential can be assigned to them).The property originating from ρN that we need in practice is a phase space distribution function, denoted by f Ik , which can be defined as Setting the indices equal, this amounts to where V denotes the volume.We remark that the normalization in eq.(2.5) differs from that in ref. [7] by a factor (2π) 3 , and is such that the total number density of sterile neutrinos, summed over the two spin states, reads Let us then turn to the Standard Model (SM) part.For its density matrix we adopt the ansatz ρSM where B is the baryon number operator (at T ≪ 160 GeV baryon and lepton numbers are separately conserved within the SM).The lepton number operator associated with flavour a reads La ≡ Here e a denotes a charged lepton of flavour a, and a L ≡ (1 − γ 5 )/2, a R ≡ (1 + γ 5 )/2 are chiral projectors.The left-handed doublet is denoted by a L ℓ a = a L (ν a , e a ) T .The interaction between SM degrees of freedom and sterile neutrinos is described by the interaction Hamiltonian Ĥint ≡ where h is a 3 × 3 Yukawa matrix with complex elements h Ia .The gauge-invariant operators to which the sterile neutrinos couple are where φ ≡ iσ 2 φ * is the conjugate Higgs doublet.For the moment the only property that is needed from SM dynamics is the spectral function corresponding to these operators, In practice we assume this function to be flavour-diagonal, i.e. ∝ δ ab , however the weight is flavour-dependent because charged lepton masses play an important role. 2or describing the dynamics of the coupled system, we start from an "instantaneous" initial state, ρ(0) = ρSM ⊗ ρN , (2.12) where the SM part is from eq. (2.7) and the sterile neutrino part is from eq. (2.3).In an interaction picture, the density matrix evolves as where ĤI is the Hamiltonian corresponding to eq. (2.9) in the interaction picture.This evolution causes both the SM and sterile neutrino density matrices to evolve.For the Standard Model this change is almost negligible whereas for sterile neutrinos it is of O(1).
In order to obtain differential equations for the set of observables to be defined, the time t is considered large compared with time scales of the SM plasma, t ≫ 1/(α 2 T ), where α is a generic fine-structure constant.This guarantees that decoherence takes place in the SM part of the Hilbert space.At the same time t should be small enough to avoid secular terms.In practice, through an appropriate choice of ρSM and ρN , we can arrange things in a way that the limit t → ∞ can be taken and equilibration is correctly accounted for without secular terms (see below).An implicit assumption we make is that since the sterile neutrinos are being produced from a statistical plasma, they are decoherent, and their density matrix is assumed to retain the diagonal form in eq.(2.3); this assumption would surely be violated by terms of higher order in h.However, given that we do find a system evolving towards equilibrium, the ansatz of eq.(2.3) can be considered self-consistent.

Sterile neutrino production rate
As a first ingredient, we determine the production rate of sterile neutrinos from an asymmetric plasma.The computation can be carried out with the formalism of ref. [17]. 3From eq. (2.13) the density matrix of the system evolves as where terms containing odd powers of h have been suppressed because they are projected out later on.The density matrix of the initial state is given by eq.(2.12).The observable we are interested in corresponds to the time derivative of the expectation value of eq.(2.5), normalized like in eq.(2.4) (and summed over spin states): Inserting eq.(2.14), we are faced with 4-point functions of the creation and annihilation operators, evaluated in an ensemble defined by ρN .Given that the density matrix is assumed diagonal in flavour, momentum and spin indices and has an effectively Gaussian appearance, the 4-point functions can be reduced to 2-point functions:4 ) where the indices incorporate all dependences.Subsequently a somewhat tedious analysis, tracing the steps outlined in ref. [17], shows that all terms quadratic in the distribution functions drop out.The final result reads where n F is the Fermi distribution.On the right-hand side we have replaced f Ik (0) with f Ik (t), i.e. ρI (0) with ρI (t) in eq.(2.14).This is consistent with our goal of an O(h 2 ) determination of f Ik , since the difference between ρI (0) and ρI (t) is of higher order in h.
As was anticipated at the end of the previous section, this prescription removes any secular terms from the evolution.Finally, we remark that, in the limit f Ik → 0, the result agrees with ref. [7], whereas in equilibrium, corresponding to µ a = 0, f Ik = n F (E I ), the production rate vanishes as must be the case.

Lepton number washout rate
The lepton number defined by eq.(2.8) is not conserved because of the interactions in eq.(2.9).The operator equation of motion reads [19] La = where (T a ) ij ≡ δ ai δ aj .We propose to evaluate the expectation value of this operator (taken in the interaction picture) in the time-dependent ensemble described by eq.(2.13).This time the leading contribution comes from the term linear in ĤI (n a ≡ L a /V ): where X ≡ (t, x), Y ≡ (t ′ , y), and ... ≡ Tr [(...)ρ(0)].The SM part of this expectation value can be expressed in terms of Wightman correlators, which in turn can be expressed in terms of the spectral function in eq.(2.11).For the sterile neutrino part, we can insert the field operators from eqs. (2.1), (2.2), and eliminate the annihilation and creation operators through eq.(2.4).A tedious but straightforward analysis yields where we have employed the prescription described below eq.(2.21).In the limit f Ik → 0, the result corresponds to eq. (2.32) of ref. [7].On the other hand in equilibrium, i.e. µ a = 0, f Ik = n F (E I ), the washout rate vanishes as must be the case.An interesting crosscheck of eq. ( 2.24) can be obtained by putting sterile neutrinos in equilibrium (f Ik → n F (E I )) and by expanding lepton number densities to first order around equilibrium.In the limit of small µ a , we can re-express the µ a through the lepton densities through a susceptibility matrix describing how lepton asymmetries disappear (or are "washed out") near equilibrium.The coefficients read which after an appropriate adjustment of conventions agrees with eq.(4.7) / (29) of ref. [19].

Relation of lepton densities and chemical potentials
In order to close the set of equations, we need to determine the relations between lepton densities and chemical potentials.Even though the densities n a are separately conserved within the Standard Model with h = 0, their fluctuations are correlated.The reason is that, because of charge neutrality, an excess of electrons over positrons is compensated for by an excess of antimuons over muons.This implies that the relation of µ a and n a is nondiagonal.In the present section we work it out to leading order in Standard Model couplings, at T < ∼ 1 GeV.(Recently, similar computations have been extended up to higher orders in Standard Model couplings at T > ∼ 160 GeV [19,20].) The desired relations can most conveniently be obtained by first computing the pressure (i.e.minus the grand canonical free energy density, p = −Ω/V ) as a function of chemical potentials associated with all conserved charges.Apart from baryon and lepton numbers, chemical potentials need to be assigned to gauge charges, in our case the electromagnetic U(1) em charge (≡ µ Q ) and the weak SU(2) L charge (µ Z ) [21].The chemical potentials assigned to gauge charges correspond to the fact that the zero components of the associated gauge fields can develop an expectation value.In our case, the weak gauge bosons can be omitted, because their effective potential has a large tree-level term ∼ µ 2 Z v 2 , which implies that µ Z is suppressed by ∼ T 2 /v 2 with respect to µ Q .Therefore, only µ a , µ B , and µ Q need to be included.Given that the chemical potentials are very small compared with the temperature (µ a < ∼ 0.02T ), it is sufficient to determine p up to O(µ 2 ) in the chemical potentials.
Omitting exponentially suppressed contributions from the W ± gauge bosons and from the top quark, the pressure can be expressed as where N c = 3 is a book-keeping variable for hadronic effects, and χ is a (diagonal) "susceptibility".In the free limit the susceptibility reads For vanishing mass this evaluates to χ(0) = T 2 /6, whereas for a non-vanishing mass it can be expressed as where K 2 is a modified Bessel function.(Radiative corrections to diagonal quark susceptibilities have been computed up to a high order in the massless limit [22], and the same quantities can also be measured on the lattice, cf.e.g.refs.[23,24] and references therein.)Given eq.(2.28), µ B is eliminated in favour of the baryon density through n B = ∂p/∂µ B , yielding where we have denoted In the following we neglect n B in comparison with lepton densities, n B ≈ 0, which fixes µ B (if n B were kept non-zero, a Legendre transform should be carried out to an ensemble with fixed n B ). Charge neutrality is imposed by requiring ∂p/∂µ Q = 0. From the resulting expression, we can obtain lepton densities as functions of the chemical potentials as n a = ∂p/∂µ a .
The solution to the charge neutrality condition ∂p/∂µ Q = 0 reads The total lepton density of flavour a is The first term accounts for the neutrino density whereas the latter term represents the density n ea of charged leptons of flavour a.
It is important to note that because the three lepton asymmetries are independent of each other, charge conservation can be balanced by the other lepton flavours.For instance, even if n ντ = 0, so that µ τ = 0, it is still possible to have n τ = 0 because tau-leptons couple to µ Q and can thereby neutralize some of the charges carried by electrons and muons.
If, however, the different lepton densities are assumed to equilibrate (which is physically unlikely at T > 10 MeV [15,16]), the situation changes.We consider this case as well, given that it can serve as a useful test case leading to the largest possible resonance effect [7].For equilibrated lepton densities we set µ a = µ L , a ∈ {1, 2, 3}, so that eq.(2.33) simplifies to and the lepton density is The neutrino and charged lepton asymmetries read n ν a = χ(0) µ L and n ea = 2χ(m a )(µ L − µ Q ), respectively.Note that if all quark masses are made large compared with temperature, or we set N c → 0, then eq. ( 2.36) implies that µ Q = µ L , and lepton asymmetry is only carried by neutrinos.This is because in the absence of hadronic plasma constituents, charge neutrality would be violated if charged leptons were chemically equilibrated and had non-vanishing asymmetries.

Summary of the basic setup
For simplicity, we assume that during the period under consideration only one sterile neutrino is active, i.e. has an interaction rate comparable with the Hubble rate. 5Then only one sterile neutrino flavour contributes in eq.(2.24).We denote this flavour by I = 1.With K being the associated on-shell four-momentum, two rates are defined by where ρ aa is the Standard Model spectral function from eq. (2.11).Then the basic equations to be solved, (2.21) and (2.24), take the forms ( with µ a 's and n a 's related through eq.(2.34).The Fermi distributions in eq.(3.2), which are always below unity, take care of Pauli blocking.The two terms on the right-hand sides of these equations can be interpreted as originating from reactions involving SM "particles" and "antiparticles". 6In the case of equilibrated active flavours, eq.(3.3) is replaced with Our goal is to solve these equations (generalized to an expanding background) in the regime where v denotes the Higgs vacuum expectation value.

Active neutrino properties
Given that we are deep in the Higgs phase, the Higgs doublet φ can to a good approximation be replaced by its (gauge-fixed) expectation value, cf.eq.(3.5).Then the neutrino Yukawa couplings only appear through neutrino Dirac masses, defined as The spectral function ρ aa is proportional to the spectral function of an active neutrino of flavour a.This assignment refers to the weak interaction eigenbasis.The (retarded) propagator of an active neutrino is of the form [25] and the spectral function is given by ρ(K) = Im S(K + iu 0 + ), where u ≡ (1, 0) is the fourvelocity of the heat bath.The function a can to a good approximation be replaced by its tree-level value a = 1 [26].The functions b and c are often called (thermal and asymmetry induced) matter potentials, and Γ can be called (up to trivial factors) a thermal width, a damping rate, an interaction rate, or an opacity.The function b is real, even in µ a , and odd in K.The function c is real, odd in µ a , and even in K.The imaginary part Γ can to a good approximation be evaluated at µ a = 0 and is then even in K. Inserting the form of eq.(3.7) into eq.(3.1), accounting properly for various sign conventions, and dropping terms which are subleading in the regime of eq.(3.5), we obtain 6 More precisely, the rate R − a describes a transition lepton ↔ N 1 and the rate R + a a transition antilepton ↔ N 1 .Time can run in either direction.Eq. (3.2) states that sterile neutrinos can be produced from leptons and antileptons.Eq. (3.3) states that asymmetries between lepton and antilepton densities decrease by the difference of the rates felt by leptons and antileptons.
Let us now discuss the explicit forms of the functions appearing in eq.(3.7).In the regime of eq.(3.5) the imaginary part Γ originates at 2-loop level, and has been computed with account of all SM reactions in ref. [27] (there it was represented by the combination where W ) is the Fermi constant, and the dimensionless function ÎQ ∼ 1 has been tabulated for various momenta, temperatures, and lepton flavours on the web page related to ref. [27].
The function b originates at 1-loop level and was determined in ref. [26] for where α w = g 2 w /(4π) and θ w is the weak mixing angle.In analogy with eq.(3.9) we can write Because of the 1/α w factor, the dimensionless function b is much larger than ÎQ , b ∼ 80.It is plotted in fig. 1 of ref. [27] for different flavours and temperatures.The last ingredient is the function c, which incorporates effects from charge asymmetries.This function was also determined in ref. [26].Assuming chiral equilibrium and including all light SM particles, we obtain where the second line encodes the contribution of quarks, coming from tadpole diagrams mediated by the Z boson.Because the latter couples differently to up-and down-type quarks, the hadronic contribution contains a part that is not proportional to n B = 1 3 i=u,d,s,c,b n i and thus survives even if, as we assume, the baryon density vanishes.The quark densities read and correspondingly for other up-and down-type quarks, respectively.Eq. (2.31) can be used for expressing µ B in terms of µ Q , yielding (for n B → 0) We again consider two cases, corresponding to those in sec.2.4.For unequilibrated lepton asymmetries, the neutrino asymmetries are different, as given by eq.(2.35).Then µ Q can be read off from eq. (2.33).The charged lepton densities originate from the second term in eq.(2.34).In contrast, for equilibrated lepton asymmetries, all neutrino densities are equal, n ν a = χ(0) µ L , whereas µ Q can be read off from eq. (2.36).

Expanding background
In an expanding background, the left-hand sides of eqs.
where H is the Hubble rate, H = 8πe/(3m 2 Pl ), and e denotes the energy density.The inhomogeneous term can be eliminated from the equation of motion for f k by integrating along a trajectory of redshifting momentum, where s is the total entropy density, and from that for n a by normalizing by s, It is also convenient to integrate in terms of the temperature T rather than the time t.
Denoting the final moment of integration by T * ≡ 1 MeV, we get where c 2 s is the speed of sound squared.Numerical values for the thermodynamic functions appearing in these equations (e, s, c 2 s ) have been tabulated in ref. [28].Note that the righthand side of eq.(3.19) is even in charge conjugation, whereas that of eq.(3.20) is odd.For the case of equilibrated active flavours, eq.(3.20) . (3.21)

Resonance contributions
For small Γ the rates in eq.(3.8) resemble Dirac delta-functions, if the first factor in the denominator has a zero.For positive c, this can be the case with the term R + a .Denoting we can then approximate and zeros exist if The zeros are located at and the "Jacobian" reads In practice, of course, I Q is not infinitesimally small, and the resonance is not arbitrarily narrow.In this situation resonance effects interfere with non-resonant contributions.One way to account for this in a practical numerical solution is sketched in appendix A.

Relic density
Once the final spectrum f k * has been obtained through the integration of eqs.(3.19)-(3.21),we need to relate it to the present-day dark matter energy density.In practice we choose the lowest temperature of the integration to be T * ≡ 1 MeV, by which time all the source terms have switched off, whereas T 0 denotes the present-day temperature of the cosmic microwave background.Today, the sterile neutrinos are non-relativistic, so that the energy density carried by them reads The dark matter energy density can be written as (Ω dm ≡ ρ dm /ρ cr ) where ρ cr is the critical energy density and s(T 0 ) = 2 891/cm 3 is the current entropy density.Making use of the known value of ρ cr yields ρ cr /[h 2 s(T 0 )] = 3.65 eV.Recalling that Ω dm h 2 = 0.12 according to Planck data [29], and dividing eq.(3.26) by eq.(3.27), we get where we made use of the facts that d 3 k T f k T /s(T ) is temperature-independent at T ≤ T * and that s(T * ) ≈ 4.67T 3 * .So, given the known f k * , eq. (3.28) allows us to determine Ω 1 /Ω dm .

Numerical results
We have integrated eqs.In case (a) Y e , Y µ grow initially, even though the source terms R ± a are not active yet, because charged τ -leptons cannot carry their share of the asymmetry when T ≪ m τ (Y L ≡ a Y a is constant).In case (e) such a re-distribution is not possible and Y µ and Y τ are exactly conserved.The values of the initial neutrino asymmetries n νa are given in table 1; the values of the corresponding lepton asymmetries n a = n νa + n ea follow from eqs. (2.33)-(2.35).Lepton asymmetries would be expected to equilibrate below T = 10 MeV [15,16], in the region shown by a grey band, however the rates R ± a have switched off by then so this has no effect on sterile neutrino distributions.
Let us reiterate that in the case of equilibrated active flavours, one would have to assume active neutrino oscillations to proceed much faster than the processes considered in the present paper, which is unlikely to happen at T > 10 MeV [15,16].Nevertheless we display the results in order to allow for a comparison with ref. [7], to be performed in sec. 5.The initial state is parametrized by the neutrino asymmetry normalized to the entropy density, n νa /s.The mixing angles are parametrized through which is the combination that appears in the (inclusive) decay rate of sterile neutrinos to an active neutrino and a photon.We consider the value sin 2 (2θ) ≈ 7 × 10 −11 mentioned in ref. [10] and the limits of sin 2 (2θ) ∼ (2 − 20) × 10 −11 from ref. [11].Confining effects are modelled through the phenomenological replacement N c → N c,eff as suggested in ref. [27].(In ref. [27] it was checked that this recipe is consistent with Chiral Perturbation Theory at low T ; unfortunately Chiral Perturbation Theory is not applicable at T > ∼ 100 MeV.)In fig. 1, the two resonance locations (in each channel) are shown for the cases (a) and (e).In fig.2, the evolution of the densities Y a is shown, and in fig. 3 the same is done for the distribution function f k T .The ratio Ω 1 /Ω DM from eq. (3.28) is illustrated in fig.4, whereas the differential shape of f k T at T = T * = 1 MeV can be inferred from figs. 5

Conclusions
In view of the exciting (if unconsolidated) prospect of accounting for dark matter through 7.1 keV sterile neutrinos [10,11], the purpose of this paper has been to promote a previously proposed quantum field theoretic framework [7] from a qualitative towards a more quantitative level.In order to reach this goal, two types of "back reactions" (i.e.non-linearities) entering the basic equations have been derived from stated assumptions (cf.secs.2.2, 2.3).The relation of the lepton densities and lepton chemical potentials entering these equations has been systematically worked out to leading order in small chemical potentials (cf.sec.2.4).
The equations have been written in a form which separately tracks three different flavours of non-equilibrated lepton densities (cf.eqs.(3.19), (3.20)).Finally, the equations have been numerically solved "as is", without imposing further model assumptions at this stage.
In a previous study [7], which relied otherwise on similar approximations as the present one, it was assumed that all active flavours are in chemical equilibrium, and that both the charged and the neutral leptons carry the same asymmetry, so that the total initial lepton asymmetry is effectively n L = 9n νe .This leads to a large coefficient c and correspondingly to a maximally efficient resonant contribution.In reality, as we have discussed, charged leptons cannot carry such large asymmetries because of electric charge neutrality.Concretely, this implies that we need a larger active-sterile mixing angle or initial asymmetry for a comparable effect.For instance, for case (a), where we find the initial asymmetry n νe /s = 12.25 × 10 −6 for sin 2 (2θ) = 7 × 10 −11 , the analysis of ref. [7] would have produced the correct dark matter abundance already with n νe /s = 9.05×10 −6 for sin 2 (2θ) = 7×10 −11 , or already for sin 2 (2θ) = 1.5 × 10 −11 with n νe /s = 12.25 × 10 −6 .In other words, the difference between ref. [7] and the present work is of order unity.On the logarithmic scale of figs.5, 6 the distributions of ref. [7] do however bear some similarity with ours, if considered at the same value of sin 2 (2θ).
Our numerical results have been presented in sec.4. The final spectra for all the cases considered can also be downloaded from http://www.laine.itp.unibe.ch/dmpheno/.It remains a theoretical challenge to confirm whether some of the pre-existing neutrino asymmetries in table 1 can indeed be produced by mechanisms such as the one described in ref. [8].
Despite the improvements of the present paper, it should be acknowledged that the solution still contains theoretical uncertainties.The reason is that most of the sterile neutrino production takes place at temperatures of a few hundred MeV (cf.fig.3), where hadronic effects play a significant role.In our work, hadronic effects have been handled through a phenomenological recipe introduced in ref. [27], which does correctly incorporate the fact that QCD displays a rapid but smooth crossover rather than an actual phase transition.Then hadronic uncertainties remain on a level of some tens of percent as discussed previously [7].Initial neutrino number density n νx /s in units of 10 Eventually, if the sterile neutrino dark matter scenario establishes itself, many of these uncertainties can be reduced through lattice Monte Carlo measurements.As has been outlined in ref. [17] and in the present paper, lattice input is needed for the equation of state, for quark number susceptibilities, and for mesonic correlation functions in various quantum number channels.The first two ingredients would already be available but it is not clear whether including lattice input in some places and not in others would consistently improve on our results. 7Nevertheless their gradual inclusion seems to present an interesting challenge for future work.
In general, the resonances cannot be considered as infinitely narrow, and their treatment cannot be separated from that of non-resonant contributions.Consider the vicinity of a structure like in eq.(3.23).Numerical integrations take place on a discrete grid.Suppose that between two grid points, x i−1 and x i , where x can be chosen as ln(T * /T ) in eq.(3.19) and as E 1 on the right-hand side of eq.(3.20) ), we find that F has changed its sign; let x 0 denote the location of the zero.In the case of the integral on the right-hand side of eq.(3.20), the correct contribution would be ∆ ≈ whereas naively we could have estimated this through In order to correct for the error, we may subtract δ and add ∆ to the result.If x 0 is close to x i or x i−1 , neighbouring cells need to be corrected as well.
In the case of eq.where a ′ enumerates non-resonant terms; r is the resonant contribution; φ ≡ φ I Q /(F 2 + I 2 Q ); and χ ≡ χ I Q /(F 2 + I 2 Q ).Assuming the non-resonant contribution to be subleading, which can be arranged by choosing x i − x i−1 sufficiently small, the unknown value f (x 0 ) can be estimated from Inserting this into eq.(A.3) we can solve for f (x i ) in terms of the known f (x i−1 ), including now the non-resonant contributions as well.

. 23 )
If c < 0, a similar term exists in R − a ; it can only exist in one of the terms at a time.In general, there are two resonances to be considered.The simplest way to see this is to fix T and consider F as a function of E 1 .Indeed the energy dependence of the variables appearing in eq.(3.22) is simple: b is linear in E 1 whereas c is constant (cf.eqs.(3.10), (3.15)).If we write b = b E 1 , where b ≪ 1, then

Figure 1 :
Figure 1: Resonance locations in the different flavour channels for case (a) (left) and case (e) (right)[only a = e has a physical effect in these cases because only h 1e = 0].In each case we have set sin 2 (2θ) = 7 × 10 −11 and tuned the initial asymmetry to the value producing the correct dark matter abundance, cf.table 1.We have considered comoving momenta k T * ≤ 12.5T * at T * = 1 MeV; the continuous line indicates the upper edge of this range.

( 3 . 11 , 1 Figure 2 :
Figure 2: The evolution of the lepton asymmetries Y a for case (a) (left) and case (e) (right).The parameters are like in fig. 1.In case (a) Y e , Y µ grow initially, even though the source terms R ±a are not active yet, because charged τ -leptons cannot carry their share of the asymmetry when T ≪ m τ (Y L ≡ a Y a is constant).In case (e) such a re-distribution is not possible and Y µ and Y τ are exactly conserved.The values of the initial neutrino asymmetries n νa are given in table1; the values of the corresponding lepton asymmetries n a = n νa + n ea follow from eqs. (2.33)-(2.35).Lepton asymmetries would be expected to equilibrate below T = 10 MeV[15,16], in the region shown by a grey band, however the rates R ± a have switched off by then so this has no effect on sterile neutrino distributions.

11 , 1 Figure 3 :
Figure 3: Examples for the evolution of the right-handed neutrino distribution f kT for case (a) (left) and case (e) (right), assuming that f kT (T = 4 GeV) = 0.The final temperature is T * = 1 MeV, and k * ≡ k T * denotes momenta at this temperature.The parameters are like in figs. 1, 2. For smallish k * /T * most of the production takes place at the lower resonance temperature (cf.fig.1).

Figure 4 :
Figure 4: The total energy density carried by sterile neutrinos today, normalized to the dark matter density, as a function of the initial lepton asymmetry, for case (a) (left) and case (e) (right).The couplings span the range indicated by refs.[10, 11].

Ω 1 / Ω DM = 1 Figure 5 :
Figure 5: The sterile neutrino distribution function at T ≤ T * = 1 MeV, normalized to the Fermidistribution n F (k) = 1/[exp(k/T ) + 1],for the initial lepton asymmetry producing precisely the observed dark matter energy density, for case (a) (left) and case (e) (right).Given that f ≪ n F , sterile neutrinos are far below equilibrium despite their efficient resonant production.

Figure 6 :
Figure 6: Like fig. 5 but for the other cases [for ease of comparison, case (e) is reproduced here].
and 6.The initial neutrino densities yielding the correct dark matter abundances in all cases (a)-(j) are summarized in table 1.It is remarkable that despite quite different asymmetries (cf.table 1),