Minimal sterile neutrino dark matter

,

We propose a novel mechanism to generate sterile neutrinos νs in the early Universe, by converting ordinary neutrinos να in scattering processes νsνα → νsνs. After initial production by oscillations, this leads to an exponential growth in the νs abundance. We show that such a production regime naturally occurs for self-interacting νs, and that this opens up significant new parameter space where νs make up all of the observed dark matter. Our results provide strong motivation to further push the sensitivity of X-ray line searches, and to improve on constraints from structure formation.
Introduction.-The existence of sterile neutrinos, putative particles that are uncharged under the standard model (SM) gauge interactions, is extremely well motivated. For example, such sterile states provide natural candidates [1][2][3][4][5] to explain the observed tiny nonzero neutrino masses [6]. If one of the sterile neutrinos has a mass in the keV range and is stable on cosmological time scales, furthermore, it is an excellent candidate for the dark matter (DM) in our Universe [7]. A smoking-gun signal for this scenario would be an astrophysical X-ray line, resulting from DM decaying into an active neutrino and a photon [8]. Such X-ray signatures are very actively being searched for, leading to ever more stringent limits on how much sterile and active neutrinos can mix [9] (while Refs. [10,11] report a potential detection).
Sterile neutrinos can be produced by neutrino oscillations in the early Universe, which is known as the Dodelson-Widrow (DW) mechanism [12]. However, the region of parameter space where this mechanism produces an abundance of sterile neutrinos consistent with the entirety of the observed DM is excluded [13]. Alternative scenarios that remain viable include resonant production in the presence of a large lepton asymmetry [14], production by the decay of a scalar [15][16][17][18][19][20], via an extended gauge sector [21][22][23], and production by oscillations modified by new self-interactions of the SM neutrinos [24][25][26][27] or by interactions between the sterile neutrinos and a significantly heavier scalar in combination with a large lepton asymmetry [28].
All known viable mechanisms thus require new particles in addition to the sterile neutrino, either explicitly or implicitly (e.g., to ensure gauge invariance or to create a lepton asymmetry much larger than the observed baryon asymmetry). Our objective here is to propose a novel scenario that is minimal in the sense that it re-quires only a single real degree of freedom on top of the DM candidate.
Recently, some of us introduced a new generic DM production mechanism [29] that is characterized by DM particles transforming heat bath particles into more DM, thereby triggering an era of exponential growth of the DM abundance (see also Ref. [30]). As discussed there, one possibility to provide the necessary DM seed population for such a DM density evolution is through an initial freeze-in [31] period. In the scenario proposed in this article, instead, sterile neutrinos ν s are initially generated through oscillations like in the DW mechanism and subsequently transform active neutrinos ν α through the process ν s ν α → ν s ν s . We demonstrate that such a scenario generically emerges when sterile neutrinos feel the presence of a dark force, and that this opens up significant portions of parameter space for sterile neutrino DM that may be detectable with upcoming experiments.
Model setup.-A necessary requirement to realize DM production via exponential growth is σv tr σv fi [29], where σv tr is the thermally averaged interaction rate for the transmission process, i.e. the conversion of a heat bath particle to a DM particle, and σv fi is the corresponding quantity for the more traditional freeze-in process [31], where a pair of DM particles is produced from the collision of heat bath particles. Here we point out that a simple and generic way to realize this condition is a secluded dark sector [32][33][34][35] where DM particles interact among each other via some mediator φ, while interacting with the visible sector only through mixing by an angle θ. In that case, both processes dominantly proceed via the s-channel exchange of φ, resulting in σv tr ∝ sin 2 θ and σv fi ∝ sin 4 θ. We note that such 'secret interactions' of sterile neutrinos have been studied in different cosmological contexts before . Motivated by these general considerations, we concentrate in the following on a single sterile neutrino ν s interacting with a light scalar φ, both singlets under the SM gauge group. Assuming Majorana masses for ν s and the active neutrinos, ν α , the relevant Lagrangian terms are given by where repeated indices α are summed over and ν α (ν s ) are left-(right-) handed spinors. We will concentrate on the case of heavy mediators, m φ > 2 m s , for most of this article, but later also briefly discuss phenomenological consequences of lighter mediators. We checked that even for large mediator self-couplings, number-changing interactions like 2φ ↔ 4φ or 3φ ↔ 4φ do not qualitatively change our results and therefore neglect the scalar potential in all practical calculations. We further assume, for simplicity and concreteness, that ν s dominantly mixes only with the active neutrino species ν e , and that m s m α . Expressed in terms of mass eigenstates, which for ease of notation we denote by the same symbols as flavor eigenstates, the interactions of the mediator then take the form (1) with θ m αs /m s 1. The unsuppressed couplings among φ and ν s turn out to be sufficiently strong to equilibrate the dark sector during the new exponential production phase that we consider below. On the other hand, mass-mixing-induced interactions between ν s and electroweak gauge bosons are suppressed by the Fermi constant, G F , and will only be relevant in setting the initial sterile neutrino abundance.
Evolution of ν s number density.-For an initially vanishing abundance, in particular, the interactions in Eq. (1) only allow freeze-in production of ν s . While the corresponding rate scales as sin 4 θ, active-sterile neutrino oscillations at temperatures above and around Λ QCD ∼ 150 MeV, in combination with neutral and charged current interactions with the SM plasma, lead to a production rate scaling as sin 2 θ [12]. Adopting results from ms m φ sin 2 (2θ) y  Ref. [62], we use the ν s number density, n s , and energy density, ρ s ∼ p n s , that result from this DW production.
Once it is completed, and in the absence of dark sector interactions, the expansion of the Universe will decrease these quantities as n s ∝ a −3 and ρ s ∝ a −4 , respectively, where a is the scale factor. Some time later, various decay and scattering processes (cf. Fig. 1) become relevant due to the new interactions appearing in Eq. (1) and, for the parameter space we are interested in here, eventually thermalize the dark sector via the (inverse) decays ν s ν s ↔ φ. From that point on, the phase-space densities of ν s and φ follow Fermi-Dirac and Bose-Einstein distributions, respectively, that are described by a common dark sector temperature T d as well as chemical potentials µ s and µ φ . Similar to the situation of freeze-out in a dark sector [63], the evolution of these quantities is determined by a set of Boltzmann equations for the number densities n s,φ and total dark sector energy density ρ = ρ φ + ρ s : where H ≡ȧ/a is the Hubble rate, P = P s + P φ is the total dark sector pressure, and C i are the various collision operators (see Appendix for details). With φ ↔ ν s ν s in equilibrium, the chemical potentials are related by 2µ s = µ φ , allowing us to replace the first two of the above equations with a single differential equation for n ≡ n s + 2n φ . Noting that ρ ∝ a −4 andñ ∝ a −3 , both right before and after φ ↔ ν s ν s starts to dominate over the Hubble rate, the initial conditions to the coupled differential equations forñ and ρ can then be determined at the end of DW production. In order to illustrate the subsequent evolution of the system, let us consider two concrete benchmark points, cf. Tab. I, for which the sterile neutrinos obtain a relic density that matches the observed DM abundance of Ω DM h 2 0.12 [64], with a mixing angle too small to achieve this with standard DW production. As demonstrated in Fig. 2, with solid (dashed) lines for BP1 (BP2 ), this leads to qualitatively different behaviors: BP1 Here the only additional process (beyond φ ↔ ν s ν s ) where the rate becomes comparable to H, at m s /T ν ∼ 0.2 with T ν the active neutrino temperature, is ν s ν α → φ (left panel, light blue). This triggers exponential growth in the abundance for both

Rate [keV]
Dark Thermalization ν s and φ (right panel, green and orange) through ν s ν α → φ * → ν s ν s , with φ being (almost) on shell, cf. Fig. 1 (c). Once T ν m φ the transmission process becomes inefficient and the final ν s abundance is obtained. Afterwards, since both φ and ν s are non-relativistic, the dark sector temperature decreases with T d ∝ a −2 both before and after kinetic decoupling (right panel, dark gray).
BP2 In this case the larger coupling y (needed to compensate for the smaller θ) leads to another process impacting the evolution of the system: at m s /T ν ∼ 0.01, the rate for ν s ν s → φφ (left panel, cyan, and Fig. 1 (b) starts to be comparable to H. As φ predominantly decays into ν s ν s (left panel, orchid), this effectively transforms kinetic energy to rest mass by turning 2ν s to 4ν s -very similar to the reproductive freeze-in mechanism described by Refs. [28,65,66]. As expected, this leads to a significant drop in the temperature T d (right panel, black). This process becomes inefficient for T d m φ , due to the Boltzmann suppression of φ. Subsequently, the rate for ν s ν α → φ (left panel, light blue) becomes comparable to H, leading to a phase of exponential growth in the same way as for BP1.
Observational constraints.-Due to the mixing with active neutrinos, ν s is not completely stable and subject to the same decays as in the standard scenario for keVmass sterile neutrino DM. The strongest constraints on these decays come from a variety of X-ray line searches. We take the compilation of limits from Ref. [9], but only consider the overall envelope of constraints from Refs. [67][68][69][70][71][72]. Furthermore, we consider projections for eROSITA [73], Athena [74], and eXTP [75].
Observations of the Lyman-α forest using light from distant quasars place stringent limits on a potential cutoff in the matter-power spectrum at small scales, where the scale of this cutoff is related to the time of kinetic decoupling, t kd . In our scenario, this is determined by DM self-interactions and we estimate t kd from Hn s = C νsνs→νsνs [76], where the collision term C νsνs→νsνs is stated in the Appendix. A full evaluation of Lyman-α limits would require evolving cosmological perturbations into the non-linear regime, which is beyond the scope of this work. Instead, we recast existing limits on the two main mechanisms that generate such a cutoff. At times t < t kd , DM self-scatterings prevent overdensities from growing on scales below the sound horizon r s = t kd 0 dt c s /a, where c s = dP/dρ is the speed of sound in the dark sector [77]. We use the results from Ref. [78] for cold DM in kinetic equilibrium with dark radiation (with c s = 1/ √ 3) to recast the current Lyman-α constraint on the mass of a warm DM (WDM) thermal relic m WDM > 1.9 keV [79] to the bound r s < 0.34 Mpc.
Overdensities are also washed out by the free streaming of DM after decoupling. We evaluate the free-streaming length as λ fs = redshifts of roughly z ∼ 50. We translate the WDM constraint of m WDM > 1.9 keV to λ fs < 0.24 Mpc, which we will apply in the following to our scenario. We note that the WDM bounds from Ref. [79] are based on marginalizing over different reionization histories. Fixed reionization models tend to produce less conservative constraints, which however illustrate the future potential of Lyman-α probes once systematic errors are further reduced (m WDM > 5.3 keV [80], e.g., corresponds to r s < 0.09 Mpc and λ fs < 0.07 Mpc, respectively). DM self-interactions are also constrained by a variety of astrophysical observations at late times [81]. We adopt σ T /m s < 1 cm 2 /g as a rather conservative limit, where σ T is the momentum transfer cross-section as defined in Ref. [82]. Far away from the s-channel resonance, we find , largely independent of the DM velocity v. For such cross sections cluster observations [83,84], or the combination of halo surface densities over a large mass range [85], can be (at least) one order of magnitude more competitive than our reference limit of σ T /m s < 1 cm 2 /g.
Viable parameter space for sterile neutrino DM.-In Fig. 3 we show a slice of the overall available parameter space for our setup in the sin 2 (2θ) − m s plane, for a fixed mediator to DM mass ratio of m φ /m s = 3. For every point in parameter space, the dark sector Yukawa coupling y is chosen such that the sterile neutrinos ν s make up all of DM after the era of exponential growth. In the yellow band, DW production can give the correct relic abundance, including QCD and lepton flavor uncertainties [62]; the dashed brown line corresponds to the central prediction, which is the basis for our choice of initial conditions for number and energy densities of ν s . In the blue region DM will be overproduced, Ω s h 2 > 0.12, while the other filled regions are excluded by bounds from X-ray searches (gray), Lyman-α observations (orange), and DM self-interactions (violet). The white region corresponds to the presently allowed parameter space.
It is worth noting that, unlike in standard freeze-out scenarios [86], later kinetic decoupling implies a shorter free-streaming length in our case because the dark sector temperature scales as T d ∝ T 2 ν already before that point. At the same time, the sound horizon increases for later kinetic decoupling. The shape of the Lyman-α exclusion lines reflects this, as kinetic decoupling occurs later for larger values of y.
In Fig. 3 we also show the projected sensitivities of the future X-ray experiments eROSITA [73], Athena [74], and eXTP [75], which will probe smaller values of sin 2 (2θ). Similarly, observables related to structure formation will likely result in improved future bounds, or in fact reveal anomalies that are not easily reconcilable with a standard non-interacting cold DM scenario. While the precise reach is less clear here, we indicate with dashed orange and violet lines, respectively, the impact of choosing λ fs < 0.12 Mpc, r s < 0.15 Mpc, and σ T /m s < 0.1 cm 2 /g rather than the corresponding limits described above. Overall, prospects to probe a sizable region of the presently allowed parameter space appear very promising.
Discussion.-While an X-ray line would be the cleanest signature to claim DM discovery of the scenario suggested here, let us briefly mention other possible directions. For example, the power-spectrum of DM density perturbations at small, but only mildly non-linear, scales may be affected in a way that could be discriminated from alternative DM production scenarios by 21 cm and high-z Lyman-α observations [87][88][89][90]. Another possibility would be to search for a suppression of intense astrophysical neutrino fluxes due to φ production on ν s DM at rest. We leave an investigation of these interesting avenues for future work.
We stress that the parameter space is larger than the m φ /m s = 3 slice shown in Fig. 3. Larger mass ratios, in particular, have the effect of tightening (weakening) bounds on λ fs (r s ), because kinetic decoupling happens earlier, and weakening self-interaction constraints; this extends the viable parameter space shown in Fig. 3 to smaller mixing angles and allows for a larger range of m s (cf. Fig. 2 in the Appendix). Changing the interaction structure in the dark sector, e.g. by charging the sterile neutrinos under a gauge symmetry, is a further route for model building that will not qualitatively change the new production scenario suggested here.
For completeness, we finally mention that smaller mediator masses are yet another, though qualitatively different, route worthwhile to explore. For m s < m φ < 2m s the mediator is no longer dominantly produced on-shell in transmission processes, so the cross section for transmission, ν s ν α → ν s ν s , scales as y 4 rather than y 2 and larger Yukawa couplings are needed in order to obtain the correct relic density. This, in turn, implies that it may only be possible to satisfy the correspondingly tighter selfinteraction and Lyman-α constraints by adding a scalar potential for φ (because additional number-changing interactions would potentially allow an increase in the ν s abundance, similar to what happens for the dashed green curve in Fig. 2, right panel, at m s /T ν ∼ 0.02). For even lighter mediators, m φ < m s , extremely small Yukawa couplings or mixing angles would be required to prevent DM from decaying too early through ν s → ν α φ (while the decay ν s → 3ν α , present also in the scenario we focus on here, is automatically strongly suppressed as Γ ∝ y 4 sin 6 θ).
Conclusions.-Sterile neutrinos constitute an excellent DM candidate. However, X-ray observations rule out the possibility that these particles, in their simplest realization, could make up all of the observed DM. On the other hand, there has been a recent shift in focus in general DM theory, towards the possibility that DM may not just be a single, (almost) non-interacting particle. Indeed, it is perfectly conceivable that DM could belong to a more complex, secluded dark sector with its own interactions and, possibly, further particles.
By combining these ideas in the most economic way, a sterile neutrino coupled to a single additional dark sector degree of freedom allows for a qualitatively new DM production mechanism and thereby opens up ample parameter space where sterile neutrinos could still explain the entirety of DM. Excitingly, much of this parameter space is testable in the foreseeable future. In particular, our results provide a strong motivation for further pushing the sensitivity of X-ray line searches, beyond what would be expected from standard DW production.
Here we complement the discussion in the main text with additional, more technical information about the scenario that we study.
Chemical potentials.-In the main text, in Fig. 2, we showed and discussed the evolution of the dark sector temperature and number densities for the two benchmark points specified in Tab. I. Here we complement in particular the discussion of the number densities by directly showing, in Fig. 4, the evolution of the chemical potentials µ s and µ φ for these benchmark points. In both cases, DW production leads to an abundance of sterile neutrinos with average momentum of the same order of magnitude as the SM neutrino temperature at the time of production, but a highly suppressed number density compared to the SM neutrino density. This leads to a large negative chemical potential µ s = µ φ /2 −T d after thermalization in the dark sector. For BP1 (solid lines) this only changes when the exponential growth of the dark sector abundance starts.
For BP2 (dashed lines) the process ν s ν s → φφ becomes important around m s /T ν ∼ 0.01, which increases the abundance by effectively converting 2ν s → 4ν s (see main text). This transformation of kinetic energy to rest mass decreases the temperature and very efficiently increases the chemical potential. If equilibrium with the inverse reaction φφ → ν s ν s was to be established (enforcing µ φ = µ s ) this would eventually result in a vanishing chemical potential (since µ φ = 2µ s is still enforced by ν s ν s ↔ φφ). This point is never quite reached for BP2, however, because φ becomes Boltzmann suppressed due to T d m φ before that could happen.
Other mass ratios.-In the main text, we already discussed the general impact of changing the mass ratio of m φ /m s = 3 that we chose for defining our benchmark points, as well as for illustrating the newly opened parameter space for sterile neutrino DM in Fig. 3 (of the main text). In Fig. 5 we further complement this discussion by explicitly showing the situation for m φ /m s = 5. We note that, apart from the general strengthening or weakening of constraints discussed in the main text, the Lyman-α forest constraint on the sound horizon first becomes stronger when increasing from very small mixing angles, against what would be expected from an earlier kinetic decoupling due to a smaller Yukawa y. This feature can be explained by the fact that at larger y also the process ν s ν s → φφ is more important, which can decrease T d (after the decays of φ) and thereby lower the speed of sound.
Collision operators.-Here, we provide explicit expressions for the types of collision operators that we used in our analysis. As discussed in the main text, we can restrict ourselves to 1 → 2, 2 → 1, and 2 → 2 interactions. We always state the collision operators C n for the number density of some particle χ ∈ {ν s , φ} from a given reaction, assuming that χ appears only once in the initial state (for an appearance in the final state the expressions need to be multiplied by −1). Multiple occurrences can be treated by summation over the corresponding expressions, and the total collision operator is given by the sum over all possible processes. We do not state explicitly the collision operators C ρ for the energy density; these can be obtained in the same way as the C n , after multiplying the integrand with the energy E χ of χ.
We assume CP -conservation for all reactions and denote particles by numbers i ∈ {1, 2, 3, 4} (χ being one of the particles i), four-momenta by p i , three-momenta by p i , absolute values of three-momenta by p i , and energies by E i = m 2 i + p 2 i , where m i is the particle's mass. The matrix elements in our expressions are summed over the initial-and final-state spin degrees of freedom. As the dark sector is in general not non-relativistic (or satisfies µ −T d ), we fully keep spin-statistical factors and the complete Fermi-Dirac or Bose-Einstein distributions everywhere. This is particularly important for scenarios where the reaction ν s ν s → φφ is relevant (roughly y 4 × 10 −4 in Fig. 3 in the main text, cf. dashed lines in main text Fig. 2) and can lead to a change in the final sterile neutrino abundance by more than two orders of magnitude compared to a simplified treatment assuming Boltzmann distributions and no Pauli blocking or Bose enhancement. 1 Starting with inverse decays, 12 → 3, the collision operator for the number density takes the form ± p 1 m 4 3 + (m 2 1 −m 2 2 ) 2 − 2m 2 3 (m 2 1 +m 2 2 ) , where the symmetry factor κ 12→3 = 1/2 if particles 1 and 2 are identical and κ 12→3 = 1 if not, and the phasespace distribution functions are denoted by the f i (with ± denoting Bose enhancement/Pauli suppression if 3 is a boson/fermion, analogously for other particles henceforth). Similarly, for the reverse reaction 3 → 12, we find Turning to 2-to-2 reactions, with particle identifications 12 → 34, the collision operator for the number den-1 Due to the exponential dependence on y, this however only changes the required value of y for the correct DM abundance by a few percent. sity is given by with the symmetry factor κ 12→34 being 1 if neither the initial-state nor the final-state particles are identical, 1/2 if either the initial-state or the final-state particles are identical, and 1/4 if the initial-state and the final-state particles are identical. We follow a similar approach as discussed in Ref. [91] for the Boltzmann equation at the phase-space level. After integrating over p 4 using the spatial part of the δ-distribution, we can write the 3momenta in spherical coordinates as p 2 = p 2 (sin β, 0, cos β) T (11) where −1 ≤ cos β ≤ 1, 0 ≤ φ < 2π, and −1 ≤ cos θ ≤ 1, and hence p 2 4 = (p 1 + p 2 − p 3 ) 2 = p 2 1 + p 2 2 + p 2 3 + 2p 1 p 2 cos β − 2p 1 p 3 cos θ − 2p 2 p 3 (cos β cos θ + sin β sin θ cos φ) .
Matrix elements.-For completeness, we finally provide a full list of all matrix elements that are relevant for our scenario (summed over the spin degrees of freedom of all initial-and final-state particles). Note that we can safely assume CP -conservation for these. We start with the decay width Γ φ = Γ φ→νανα + Γ φ→νανs + Γ φ→νsνs (31) of the mediator φ, implicitly defining the matrix elements via the partial widths We note that, typically, Γ φ Γ φ→νsνs as sin θ 1.
In the main text, we have stressed that dark sector thermalization and the phase of exponential growth are mostly due to the 3-body interactions ν s ν s ↔ φ and ν s ν α → φ. These scale as y 2 and are equal to the onshell contribution to the s-channel part of ν 1 ν 2 → ν 3 ν 4 , cf. the first term in Eq. (36). From the full expression of the matrix element, see also Fig. 1 in the main text, it is however evident that there are also contributions from t/u-channel diagrams, as well as s-channel contributions with an off-shell mediator. The processes ν s ν s → ν s ν s and ν s ν α → ν s ν s thus have interaction rates scaling as y 4 far away from the s-channel resonance. Since ν s ν α → φ is at most barely larger than the Hubble rate (cf. left panel of Fig. 2 in the main text), the additional suppression by y 2 causes the off-shell contribution to ν s ν α → ν s ν s to be negligible.
As also discussed in the main text, it is sufficient to solve the Boltzmann equations forñ = n s + 2n φ and ρ = ρ s + ρ φ since ν s ν s ↔ φ thermalizes the dark sector. Self-scatterings ν s ν s ↔ ν s ν s do not change number or energy density, i.e. their contribution to the corresponding collision operators is zero. Similarly, ν s ν s ↔ φ does not changeñ or ρ and does not need to be calculated for the evolution. Note however that the full process ν s ν s → ν s ν s is relevant for kinetic decoupling, which can occur after φ becomes highly Boltzmann suppressed at T d m φ . * torsten.bringmann@fys.uio.no † frederik.depta@mpi-hd.mpg.de ‡ marco.hufnagel@ulb.be § joern.kersten@uib.no ¶ ruderman@nyu.edu * * kai.schmidt-hoberg@desy.de