Sterile neutrino dark matter: Impact of active-neutrino opacities

The resonant production of keV sterile-neutrino dark matter mainly takes place during the QCD epoch of the early universe. It has been argued that it could be strongly affected by the opacities (or damping rates) of active neutrinos, which receive non-perturbative QCD-contributions. We find that for lepton asymmetries $ n _ { L _ \alpha } /s $ below $ 10 ^ { -6 } $ the opacities significantly affect the sterile-neutrino yield, but that for larger asymmetries, which are necessary for producing a significant fraction of the dark matter, the yield is insensitive to changes of the opacities. Thus non-perturbative QCD contributions to the opacities at temperatures around 160 MeV will not affect this dark matter scenario. We obtain larger sterile-neutrino yields than previous studies, and thus weaker lower limits on the active-sterile mixing angle from Big Bang Nucleosynthesis.


Introduction
Sterile neutrinos with mass in the keV range have long been discussed a dark matter candidate [1]. Through their mixing with active neutrinos they can be produced from the Standard Model (SM) plasma in the early universe. They are warm dark matter, giving rise to less structure on small scales than cold dark matter, which may solve some problems of the standard ΛCDM model. The detection of small scale structure in the Lyman-α forest results in lower bounds on the mass of the sterile neutrinos [2]. The sterile neutrinos can decay into an active neutrino and a photon, giving rise to monochromatic X-rays. Non-observations of X-ray lines imposes upper limits to the active-sterile mixing angle which depend on the sterile-neutrino mass. The original production scenario proposed by Dodelson and Widrow [1], has already been ruled out as the sole source of dark matter production by combining Lyman-α and X-ray constraints [3][4][5]. One way to circumvent these constraints was suggested by Shi and Fuller [6]: They assume lepton asymmetries much larger than the observed baryon asymmetry. They lead to resonant production, resulting generally in a non-thermal spectrum, which can be colder than in the Dodelson-Widrow scenario (see, however, [7]), and so evade Lyman-α constraints. Furthermore, resonant enhancement makes the production much more efficient, requiring smaller mixing angles and thus escaping the upper limits from X-ray constraints [8][9][10].
In order to calculate the final abundances, one has to track a set of coupled evolution equations for the sterile-neutrino phase space densities and the lepton asymmetries, the latter of which get depleted during the resonant conversion process. One difficulty in solving these equations comes from high lepton asymmetries resulting in sharp resonances, which necessitate high numerical precision. Furthermore, the production process typically starts at temperatures of a few GeV and ends at a few MeV prior to the onset of Big Bang Nucleosynthesis (BBN), introducing uncertainties from the QCD epoch T ∼ 160 MeV, when QCD interactions are strong and thus non-perturbative.
The cosmological expansion depends on the pressure-energy relation (equation of state) of the matter in the universe. Lattice QCD results have been incorporated into an equation of state [11] which is used to compute sterile neutrino production [12]. The sterileneutrino reaction rate depends both on the real and the imaginary part of the self-energy of the active neutrinos. The latter receives contributions from interactions with leptons and with quarks. The quark contribution is determined by susceptibilities of the QCD plasma's conserved charges. Lattice determinations of susceptibilities [13,14] have been included into sterile-neutrino evolution equations in [15]. The most difficult to compute are hadronic contributions to the imaginary part of the active-neutrino self-energy, or opacity, which can be written as a momentum integral over mesonic spectral functions [16]. Determining these spectral functions on the lattice is challenging. Previous works on resonant sterile-neutrino production [15,17] have treated the opacity in different approximations, but it still remains an open question how important the non-perturbative contributions are. In this work we address it by computing the sterile-neutrino dark matter production for various lepton asymmetries and by carefully treating the resonances which appear during the time evolution.
After introducing our setup in sec. 2.1, we state in sec. 2.2 evolution equations for the sterile-neutrino phase space densities and lepton asymmetries. We will include various leptonic and hadronic active-neutrino self-energy contributions in sec. 2.3. We solve the coupled system of equations for various values of initial lepton asymmetries in sec. 3. Motivated by our findings we update the lower limit on the active-sterile mixing angle in the two-flavor (one active and one sterile flavor) scenario by using maximal values of lepton asymmetries allowed by BBN. Conclusions are in sec. 4. Appendix A contains additional figures showing how the resonances emerge when the lepton asymmetry is increased.
2 Non-equilibrium evolution equations

Setup
We consider the Standard Model augmented by one family of sterile Majorana neutrinos N with Majorana mass M and non-zero Yukawa couplings h α to all active neutrino flavors, where ϕ = iσ 2 ϕ * is the conjugate Higgs doublet and α = (ν Lα , e Lα ) the left-handed lepton doublet. The sterile-neutrino field in the interaction picture reads with the energy k 0 = (k 2 + M 2 ) 1/2 . The spinors u and v satisfy the Majorana condition u = v c , where c denotes charge conjugation. Furthermore, the creation-/annihilation operators fulfill {a kλ , a † qλ } = δ k,q δ λλ with helicities λ = ±1/2. V denotes the volume of our system. With these operators we define the sterile-neutrino phase space density operators as 3) The number operator for left-handed leptons of flavor α reads We will start the time evolution at temperatures of a few GeV, where sphaleron processes have long terminated. Then baryon number B is conserved, and its tiny value can be well approximated by zero for our purposes. Furthermore electric charge Q is conserved and exactly zero.

Equations of motion
The expectation values of the operators introduced above are conserved by Standard Model interactions and thus evolve much more slowly than the other degrees of freedom. Their deviations from equilibrium characterize the non-equilibrium state, and are assumed to be small. We introduce chemical potentials µ Lα , µ B and µ Q and denote them collectively by µ in the following. In the infinite volume limit the phase space densities and the lepton number densities n Lα ≡ L α /V satisfy the evolution equations [18,19] with the Fermi-Dirac distribution f F (x) = 1/(e x/T + 1) and the spinors u and v which appear in (2.2). Furthermore, is the spectral function which is determined by the retarded and advanced 2-point function of the operator J α ≡ ϕ † h α α , which couples toN in (2.1).
Oftentimes the kinetic equations are expanded in µ. However, for a proper treatment of resonances, which show up in the spectral functions ρ α , one has to include all orders in µ, which will become more apparent below. The relation of the chemical potentials to the charge densities, to be discussed in sec. 2.3, can still be assumed to be linear.
In the broken phase, where ϕ = (v/ √ 2, 0) with v = 246 GeV, the 2-point function (2.8) is proportional to the active-neutrino propagator, Here we have factored the chiral projectors P R,L = 1 2 (1 ± γ 5 ) out of the propagator and furthermore introduced the active-sterile mixing angle The active-neutrino self-energy in the plasma rest frame can be approximated as [17,20] Σ ret with real b α , c α . Γ α is the imaginary part of the refractive index for the active neutrinos [21], and is also referred to as neutrino opacity. Γ α /2 is the neutrino damping rate (see e.g. [22]). The function c α is odd in µ. For the chemical potentials under consideration we only need to keep the linear order in µ for c α , and we can neglect the µ-dependence of b α and Γ α . The advanced self-energy is obtained from (2.11) by replacing iΓ α → −iΓ α . Then one obtains the spectral function A quick calculation yields (2.14) Resonances occur when a square bracket in the denominator of (2.13) or (2.14) vanishes. This can happen when c α is large enough and has the appropriate sign. For a given sign of an initial lepton asymmetry, only one of the expressions (2.13) and (2.14) can lead to resonances. Moreover, at leading order in M/|k| only the terms containingū k− ρ α u k− orv k+ ρ α v k+ will contribute in the evolution equations. The subleading terms will be dropped in the following. Due to the isotropy of the universe the phase space density only depends on |k|. Then the Hubble expansion is taken into account by replacinġ The Hubble parameter is given by where ρ is the energy density and M Pl 1.22 · 10 19 GeV is the Planck mass. By using K ≡ |k|a(t)/a(t end ) as an independent variable, where a is the scale factor, (2.15) turns into ∂ t f Kλ and the equation for f becomes an ordinary differential equation. a(t end ) is the scale factor at the time corresponding to the temperature T end = 10 MeV at which we compute the final abundances. For lepton number densities the Hubble expansion is taken into account through the replacemenṫ (2.17) The term proportional to H is eliminated by considering the differential equation for n Lα /s where s is the entropy density. Finally, the time derivatives are replaced by temperature derivatives via dT /dt = −T H(T )3c 2 s (T ), with the speed of sound c s .

Active-neutrino self-energy
For vanishing chemical potentials the real part of the active-neutrino self-energy arises at O G F /m 2 W , where G F is the Fermi constant and m W is the W -boson mass [12,21], where θ W is the weak mixing angle, E α ≡ (p 2 + m 2 α ) 1/2 , and m α is the mass of the charged lepton of flavor α. We have neglected the masses of active neutrinos. b α is positive, which corresponds to an index of refraction greater than 1, or a negative thermal mass squared.
The leading contribution due to non-zero chemical potentials can have either sign. It arises at O G F [15,17,21], without the 1/m 2 W suppression of (2.18). Therefore it can be of similar size as (2.18) when the chemical potentials are small. n να , n eα are particle minus anti-particle number densities of neutrinos and charged leptons. They can be written in terms of the particle chemical potentials µ i , where the lepton susceptibilities χ i can be evaluated in the ideal gas limit, with g να = 1, g eα = 2. The hadronic contribution to the electric charge density n had Q can be written as where χ had QQ is the hadronic contribution to the electric-charge susceptibility. The particle chemical potentials in (2.20) can be written in terms of the chemical potentials of the slowly varying and of the conserved charges, The latter can be expressed through the lepton number densities n Lα by inverting The susceptibilities χ had QQ , χ BQ , and χ BB have been determined on the lattice for temperatures near the QCD crossover [13,14]. Reference [15] has used a hadron resonance gas model below and perturbation theory above and connected all three regions via spline interpolations, which we are going to use. 3 The dominant contribution to the active-neutrino opacity Γ α appears at O G 2 F since the O G F contributions are suppressed by exp(−m W /T ). It can be split into a leptonic and a hadronic piece, Over a large part of the temperature range which is relevant for sterile-neutrino production, Γ had α is non-perturbative. Two different approaches have been taken to calculate this function. In [12] the free-quark approximation is used for the whole temperature range, but, in order to account for the strong interaction, the number of colors N c is replaced by a temperature dependent N c,eff (T ) which vanishes at low temperatures, and equals 3 at the highest temperature. In [15], on the other hand, the free-quark approximation at high temperatures is connected to chiral perturbation theory at low temperatures via spline interpolations. We refer to these two approximations for (2.32) as Γ N c,eff and Γ spline . We will also consider the approximation Γ had α = 0 for which we write Γ lep .

Numerical results
We integrate the coupled eqs. (2.5), (2.6) from T = 4 GeV down to 10 MeV, using the parameterizations of the energy ρ(T ) and entropy density s(T ) as well as the speed of sound c s (T ) from [12], 4 based on calculations in [11].
Depending on the value of the lepton asymmetry, high numerical precision is needed to handle sharp resonances. In practice, the numerical solution requires an implicit multistep method to deal with possible equation stiffness. We assume that at the starting temperature there are no sterile neutrinos, i.e., f kλ = 0. We recast the mixing angles into the total active-sterile mixing angle through which is the quantity that one can constrain experimentally from X-ray observations.

Influence of the opacity on sterile-neutrino production
In this section we make the simplifying assumption that only one Yukawa coupling is nonzero, namely h µ . This allows us to compare the effect of the opacities Γ N c,eff α , which are available for 3 lepton flavors, and Γ spline α , which is currently only available for α = µ. Then only the term with α = µ contributes in (2.5), and only the muon-flavor asymmetry will be dynamical. In principle, non-zero electron and tau-flavor asymmetries can influence the evolution equations as they appear in the functions c α , but we assume these to be zero. We choose M = 7.1 keV, motivated by the tentative signal reported in [23,24], and θ 2 µ = 2.5 · 10 −13 as a representative point in the available parameter space. We compare the resulting sterile-neutrino phase space densities at T = 10 MeV obtained with the different opacities Γ spline is available for 0.03 ≤ |k|/T ≤ 12.5. We use the latter range when solving the evolution equations, which is sufficient for our purposes. For positive lepton asymmetries, resonances mainly contribute to the production of sterile neutrinos with negative helicity, while for positive helicity the resonant contribution is suppressed with M/|k|. If there are resonances, then there are usually two resonance frequencies for each |k| [8,15,17]. For most of the relevant temperatures, the two resonances lie in the momentum range we consider. In practice, the smaller resonance frequency dominates the sterile-neutrino dark matter production, and the larger one plays a negligible role [15]. We show results for three different initial values of n Lµ /s in figs. 1, 2 and 3, additional ones can be found in appendix A. Generally we observe that the higher the initial lepton asymmetry, the larger the phase space densities become. In fig. 1 the lepton asymmetry is so low that resonances are outside the displayed momentum range (the dominant one leads to the slight increase at small momenta) and only give a small contribution to the production. In contrast, the initial lepton asymmetry in fig. 2 is high enough so that each momentum mode in the shown range passes through a resonance, giving much larger phase space densities. The same is true for fig. 3, where we chose the initial asymmetry such that we obtain the complete dark matter abundance. In fig. 1 we see how the different approximations for the opacity influence the sterile-neutrino production and in parallel the depletion of n Lµ /s. The purely leptonic contribution is the smallest one, resulting in the least efficient production. In figs. 2 and 3 one can see that for a larger initial lepton asymmetry, there is only a sub-percent difference in the final abundance of sterile neutrinos between using the full opacity and using only the leptonic contribution. The resulting phase space densities have become indistinguishable.
In the limit Γ α → 0 (2.13) turns into a delta function [17]. This indicates that the dominant resonance in figs. 2 and 3 is so sharply peaked, that the differences in the active-neutrino opacities become irrelevant. The same is true for the lepton asymmetry evolution. No matter what opacity is used, the depletion is almost identical. We find that this behavior occurs in all of the allowed (white) parameter space shown in fig. 4, if we  fig. 1 but with an initial lepton asymmetry tuned such that the sterile neutrino energy density gives the complete relic dark matter abundance, Ω s = Ω DM .
tune the lepton asymmetry such that the resulting sterile neutrino energy density gives the correct dark matter abundance. The differences in energy densities obtained with the different opacities are typically below the 2%-level, for very low masses and high mixing angles at most 5%. The transition from quite different to basically equivalent solutions by increasing the lepton asymmetry can be followed in smaller steps in appendix A. We have used the publicly available code of [17] to check our calculation, and what we find is mostly in agreement with our results described above. For very high asymmetries we find that the resulting phase space densities suffer from sporadic kinks, hinting at numerical instabilities which we could not get rid of by naively increasing the desired precision. Nevertheless the resulting figures resemble ours quite well.
Our findings partly disagree with the ones in [15], which were calculated using steriledm, a publicly available code created by the authors of [15]. It uses 1,000 momentum bins as a default, which apparently misses parts of the resonances in the sterile-neutrino production. While this problem is absent for non-resonant production, it becomes more and more severe for increasing asymmetry. We have explicitly checked that increasing the number of momentum bins to 30,000 gives results which mainly agree with ours. This problem could be the cause of the rather large differences in the phase space densities using either Γ N c,eff α or Γ spline α which was observed in [15]. The lepton asymmetries needed to produce the complete dark matter abundance (Ω s /Ω DM = 1) are quite large compared to the baryon asymmetry. Most baryogenesis mechanisms produce comparable amounts of lepton and baryon asymmetries before electroweak sphaleron freeze-out. In the νMSM [25], which contains two additional heavier sterile neutrinos, a larger lepton asymmetry can be produced thereafter [9,26]. However, it turns out this can boost the lepton asymmetry by at most a factor 1,000 [27], and that one can reach at most Ω s /Ω DM = 1/10 [18]. In this scenario, improving on calculations for Table 1: Mixing angle that leads to the complete relic DM abundance for various masses, with total initial asymmetry n L /s = 2.5 · 10 −3 . Left: asymmetry only in the muon flavor, right: all three asymmetries initially equal, n Lα /s = n L /3s. M/keV sin 2 (2θ) · 10 13  hadronic contributions to active-neutrino opacities can be important [18], as the lepton asymmetries are even smaller than in fig. 1.

Limits from Big Bang Nucleosynthesis
Measuring the helium abundances as a result of BBN one obtains the upper bound [4] n L s ≤ 2.5 · 10 −3 (3.2) on the total lepton asymmetry. We now compute the resulting lower limit on the mixing angle, for which one can obtain Ω s /Ω DM = 1. For simplicity, we take the maximal value in (3.2) as an initial condition at T = 4 GeV. The lepton asymmetry at times prior to BBN could be higher, but the depletion turns out to be only on the level of a few percent for the low mixing angles considered here. We calculate, for various masses, the mixing angle that leads to the complete relic DM abundance. The results are given in table 1, displaying the scenario where all asymmetry is in the muon flavor (left) and the scenario where the asymmetry is split equally onto all three flavors (right). The only non-zero neutrino Yukawa coupling is h µ .
Combining our lower limit on the mixing angle with current X-ray constraints closes the parameter space for masses M 40 keV. A lower mass limit of M 1 keV can be deduced from DM phase space density restrictions within dense galaxies [28,29]. In the low mass range, very high mixing angles in the Dodelson-Widrow scenario are excluded because they result in dark matter overproduction [12]. Therefore the parameter space is bounded from all sides. We show the combined constraints in fig. 4. The applicability of lower mass bounds from Milky Way satellite counts or Lyman-α methods in the case of resonant production, which typically leads to colder than thermal spectra, is not clear at this stage. Hence we did not include them in our parameter plot. Previous calculations of X-Ray  Figure 4: Combined constraints for keV sterile-neutrino dark matter. The X-ray constraints are taken from [30], see also [10,[31][32][33]. Phase space density constraints are from [29]. The BBN limit given by the solid black line holds if all of the input lepton asymmetry is only in the muon flavor. The dashed line corresponds to the BBN limit if the input lepton asymmetry is split equally onto all three flavors.
such mass bounds have been performed e.g. in [7,34] with the code sterile-dm. Generally we find colder spectra than [7], with a mean momentum that is typically between 25% and 50% lower, depending on the region in the parameter space.
As we have seen in sec. 3.1, our calculations generally give larger phase space densities than sterile-dm, if it is used for resonant production "as is" with 1,000 default momentum bins. This code was used to calculate the BBN limit in [7,[33][34][35], giving much stronger limits than the ones we find, especially for the lower end of the mass range shown in fig. 4. Again, increasing the number of momentum bins, sterile-dm gives better agreement with our results. On the other hand we note that our BBN limits are in closer agreement with the only slightly lower ones in [4], which are also displayed in [32,36], and also with the ones shown in [37], which are based on [38].
One has to keep in mind that the BBN bound only applies to the total lepton asymmetry. In fact, n Lµ could be larger than (3.2) if it is partly compensated by the other lepton flavor asymmetries. But the same compensation would not take place in (2.19) where the different flavors enter with different coefficients. Therefore c µ would increase leading to a larger production rate and to a weaker bound on sin 2 (2θ).

Summary and conclusions
In this work we have traced the evolution of keV sterile-neutrino phase space densities and lepton number densities in the early universe. Lepton asymmetries much higher than the baryon asymmetry significantly influence the active-neutrino spectral function and resonantly boost sterile-neutrino production.
Standard Model input enters our calculation in several places: through susceptibilities, which relate charges to chemical potentials, and through spectral functions of various currents, which determine the opacities of active neutrinos. The sterile-neutrino dark matter production mainly happens during the cosmic QCD epoch when quarks and gluons are strongly coupled. Then both the hadronic susceptibilities and hadronic spectral functions are non-perturbative. For hadronic susceptibilities, lattice QCD results are available and used in our calculation. Lattice calculations of hadronic spectral functions are notoriously difficult, which in principle could lead to large theoretical uncertainties.
We have studied how different approximations for the hadronic contributions to the active-neutrino opacity affect sterile-neutrino production. For initial conditions with vanishing sterile neutrino density, and lepton asymmetries large enough so that the produced sterile neutrinos make up all of today's dark matter, we find that the production is dominated by very sharp resonances for which the effect of the hadronic opacities are negligible. This implies that more precise, non-perturbative determinations will not be needed for this scenario.
For such large lepton asymmetries we found much larger dark matter yields than previous studies. Therefore we could weaken lower bounds for the mixing angle derived from the upper bound on the total lepton asymmetry from BBN. This opens up the available parameter space of this sterile-neutrino dark matter scenario.
it becomes blind to the choice of the active-neutrino opacity.