Production of keV Sterile Neutrinos in Supernovae: New Constraints and Gamma Ray Observables

We study the production of sterile neutrinos in supernovae, focusing in particular on the keV--MeV mass range, which is the most interesting range if sterile neutrinos are to account for the dark matter in the Universe. Focusing on the simplest scenario in which sterile neutrinos mixes only with muon or tau neutrino, we argue that the production of keV--MeV sterile neutrinos can be strongly enhanced by a Mikheyev--Smirnov--Wolfenstein (MSW) resonance, so that a substantial flux is expected to emerge from a supernova, even if vacuum mixing angles between active and sterile neutrinos are tiny. Using energetics arguments, this yields limits on the sterile neutrino parameter space that reach down to mixing angles of the order of $\sin^2 2\theta \lesssim 10^{-14}$ and are up to an order of magnitude stronger than those from X-ray observations. While supernova limits suffer from larger systematic uncertainties than X-ray limits they apply also to scenarios in which sterile neutrinos are not abundantly produced in the early Universe. We also compute the flux of $\mathcal{O}(\text{MeV})$ photons expected from the decay of sterile neutrinos produced in supernovae, but find that it is beyond current observational reach even for a nearby supernova.

One of the most auspicious candidate particles for the dark matter in the Universe is the sterile neutrino-an electrically neutral fermion with a mass on the order of keV-MeV that couples to ordinary matter only through a tiny mass mixing with Standard Model (SM) neutrinos [1,2]. In the simplest sterile neutrino scenarios, it is assumed that the abundance of sterile neutrinos ν s is zero at the end of inflation, and they are later produced through their mixing with SM (active) neutrinos ν a [3,4] (see also [5]). Experimental constraints on the mass of keV sterile neutrinos and their mixing with SM neutrinos arise from the measured DM relic density [6], from Pauli blocking (the Tremaine-Gunn bound) [7,8], from Lyman-α forests [9], and from X-ray searches for radiative decays of sterile neutrinos ν s → ν a +γ [10][11][12][13]. From the combination of these constraints, one concludes that the ν s mass should be m s 4 keV, and its mixing angle with the SM neutrinos should be sin 2 2θ 10 −6 in the simple two-flavor approximation.
In this letter, we add a new limit to this inventory of constraints by considering sterile neutrino production in core-collapse supernovae (SN) [14][15][16][17][18][19][20][21][22][23]. A supernova develops when a 9M star runs out of nuclear fuel. The thermal pressure that normally counteracts gravity disappears, and the core of the star collapses into a neutron star. The temperature in the nascent neutron star is MeV, so that a thermal population of (active) neutrinos is produced. These ν a can oscillate into ν s , which escape the exploding star unhindered and may carry away significant amounts of energy [16,20,24]. Constraints on anomalous energy loss from SN 1987A will thus allow us to constrain the sterile neutrino parameters. Since ν a → ν s conversion can be resonant thanks to the ultrahigh matter density ∼ 10 14 g/cm 3 in the SN core, these limits will be very strong. The flux of sterile neutrinos with O(MeV) energies escaping from a supernova leads to a flux of secondary gamma rays when they decay, and we study this flux as well.
Our main results are summarized in fig. 1. The solid orange exclusion region shows that, in the mass range m s ∼ 2-80 keV, limits from energy loss in supernovae surpass previous limits by up to two orders of magnitude in sin 2 2θ. Note that, unlike the other limits shown in fig. 1, our bounds would still hold if sterile neutrinos are not part of the DM in the Universe. In the following, we discuss in detail how we have obtained our new limits and sensitivity estimates.
Sterile neutrino production in supernovae. We consider a simplified two-flavor oscillation picture with mixing between a sterile neutrino ν s and one species of active neutrinos ν x , x = µ or τ . We do not consider mixing between ν s and electron neutrinos to avoid complications arising from charged current interactions between ν e and electrons/positrons. The flavor basis Hamiltonian describing neutrino propagation in matter includes the vacuum oscillation term and a Mikheyev-Smirnov-Wolfenstein (MSW) potential V x that describes coherent forward scattering of ν x on the background matter via Z exchange [29][30][31]: Here, θ 1 is the ν s -ν x mixing angle in vacuum and ∆m 2 m 2 s is the mass squared difference between the two mass eigenstates [32]. The MSW potential is V where G F is the Fermi constant and N n , N νe , and Nν e are the neutron, electron neutrino, and electron antineutrino number densities. The + (−) sign corresponds to the potential experienced by neutrinos (antineutrinos). Since |N νe − Nν e | N n /2, the potential is positive for antineutrinos and negative for neutrinos. Note that we neglect terms that would arise from differences between the ν x andν x number densities because such differences are expected to be small in the  Supernova bounds on the sterile neutrino mass ms and mixing sin 2 2θ (orange, this work) compared to previous constraints [12,13,25] from the Tremaine-Gunn bound [7,8] (yellow), from X-ray searches in the Andromeda Galaxy M31 [10] (green), from the Fermi GBM all-sky analysis [26] (light blue), and from the galactic [25] (dark blue) and extragalactic [11][12][13] (gray) diffuse X-ray background. For the latter constraint, we show also how it is modified if sterile neutrinos account for only a fraction of the DM in the Universe. It is understood that the other X-ray limits would shift in a similar way if Ωh 2 < 0.11. We also show in pink the region preferred by the properties of the core of the Fornax dwarf galaxy [27], and as a yellow star the parameter point that could explain the 3.5 keV X-ray line hinted at by refs. [13,28]. Black curves illustrate the parameter regions in which the Dodelson-Widrow mechanism [3] (L = 0) or the Shi-Fuller mechanism [4] with lepton assymetry L > 0 would yield the correct DM relic density Ωh 2 = 0.11. parameter region where our limits will lie [24]. We have checked that even adding an MSW potential corresponding to a maximal asymmetry does not alter our results significantly. We also neglect the momentum and angular dependence of neutrino self-interactions [33,34] and instead restrict ourselves to the simplified formalism used in [16,17,20]. Since N n , N νe , and Nν e are extremely large in the supernova core and gradually decrease with radius, most antineutrinos will encounter an MSW resonance on their way out of the exploding star. At the resonance, and the effective mixing angle θ m in matter becomes maximal [35]. We consider two physically different mechanisms for sterile neutrino production in supernovae: (i) Adiabatic flavor conversion at an MSW resonance.
A ν x of energy E streaming away from the supernova core can convert to ν s when the matter density, and thus the MSW potential V x , has reached the value satisfying the resonance condition, eq. (2). Each point (t, R) in time and space corresponds to a specific value of V x ; therefore, neutrinos of a specific energy E res (t, R) are resonantly converted at this point. We take the local neutrino luminosities and spectra and the local matter densities from the simulation of an 8.8M supernova by the Garching group [36] (see also [37][38][39][40]).
Hard scattering processes must be rare to give neutrinos sufficient time to convert adiabatically. We therefore require the spatial width of the MSW resonance regions, to be smaller than the mean free path λ mfp . The number dN MSW s /dE of ν s in an energy interval [E, E + dE] produced by adiabatic flavor conversion at the MSW resonance is given by [16] dN MSW Here, the time integral runs from the time of core bounce (t = 0) until ∼ 9 sec later, and R res (t , E) is the radius at which the resonance energy is E at time t . The quantity n ν (t , R res ) is the active neutrino number density at time t and radius R res , f x (E) is the energy distribution of active neutrinos,Ē is the average neutrino energy, P res (E) is the flavor conversion probability at the resonance, and the Heaviside Θ function implements the condition that neutrinos must have enough time between collisions to convert adiabatically. It is crucial that neutrinos do not encounter more than one MSW resonance on their way out of the supernova. (This would be different if we considered mixing between ν s and ν e instead of ν µ,τ [18,19].) Note that we can make the strongly simplifying assumption of radial symmetry, and we also neglect the depletion of active neutrinos by conversion into ν s . Moreover, we do not need to consider ν x streaming inwards. They would convert to ν s at the resonance, then travel through the core, and convert back to ν x on its far side. We parameterize f x (E) as [41] f with normalization (This relation definesĒ.) The "pinching parameter" α desribes the degree to which f x (E) differs from a Maxwell-Boltzmann distribution.  [36]. Contributions from adiabatic flavor conversion and collisional production are shown separately. The νs spectrum extends to much higher energies than the νx spectrum because flavor conversion occurs in a region where temperatures are much higher than at the last scattering surface of νx.
The ν x → ν s conversion probability at the resonance is given by the Landau-Zener formula [16,17,35] with the oscillation length at the resonance, L res osc 2π/(V x sin 2θ). We find that adiabatic flavor conversion occurs mostly at radii ∼ 10-15 km, still inside the neutrino sphere at ∼ 20-30 km. In fig. 2, we compare dN MSW s /dE to the spectrum of active neutrinos. We see that the ν s spectrum extends to higher energies because most of the flavor conversion happens in a high-temperature region from which ν s can stream out freely, while ν x are still trapped.
(ii) Collisional production. Sterile neutrinos can also be produced from a mixed ν x -ν s state in hard scatterings on nucleons, electrons, and positrons. We take the interaction rate Γ x to be approximately equal to the dominant scattering rate on neutrons, and the cross section for this process is computed following [42,43]. The physical picture for collisional ν s production is as follows: starting with an ensemble of only ν x at t = 0, each of them soon acquires a small ν s admixture ∝ sin 2 2θ m by oscillation. A collision causes the collapse of the resulting mixed state into either a pure ν x or a pure ν s . Afterwards, oscillations start anew, and ν x are quickly replenished. After many collisions, the ν s abundance is proportional to sin 2 2θ m and Γ x .
Up to a factor of 1/2, this intuitive picture leads to the correct quantum mechanical Boltzmann equation [1,3,15,44,45] ∂ ∂t Here, dn s /dE and dn x /dE are the energy spectra of sterile and active neutrinos, respectively. At any given spacetime point, dn x (t, R, E)/dE is related to the distribution function in eq. ( The averaged oscillation probability P νx→νs is [1] P νx→νs = 1 2 The extra term (Γ x E/m 2 s ) 2 in the denominator compared to the usual expression for the mixing angle in matter [35] accounts for the suppression of ν s production when λ mfp is much smaller than the oscillation length, so that oscillations do not have time to develop between collisions (quantum Zeno effect).
Integrating eq. (7) over time and radius leads to the energy spectrum of sterile neutrinos produced collisionally, We again evaluate dN s /dE numerically using the data from [36]. The resulting ν s spectra, shown in fig. 2, can be harder than the ones from adiabatic production because the collisional production rate depends on Γ x , which grows proportional to E 2 . Constraints from supernova luminosity. We can constrain the energy output in sterile neutrinos from SN 1987A by comparing the observed energy output in active neutrinos of E a = few × 10 53 ergs [46,47] to the gravitational energy released in the collapse of a stellar core at the Chandrasekhar mass, which is also on the order of E tot = few × 10 53 ergs [36,48]. If a substantial fraction of E tot was carried away by sterile neutrinos, the observed E a could not be explained [49]. We therefore consider the ratio R(sin 2 2θ, m s ) ≡ E s (sin 2 2θ, m s )/E tot . We assume that R depends only weakly on the mass and type of the progenitor star, so that the values obtained for the supernova simulated in [36] are a good proxy for SN 1987A. Our computation of ν s production is only selfconsistent for R 1 because we neglect depletion of active neutrinos. We nevertheless extrapolate it to larger values and set a limit by requiring R < 1. The justification for this is that the associated uncertainty in our results is on the same order as the uncertainty of the predicted and measured energy outputs from SN 1987A. We have verified this assumption by rescaling in each time step of our calculation the active neutrino number densities to account for the energy carried away by sterile . 3. Constraints on the sterile neutrino parameter space from energy loss in supernovae, considering (a) adiabatic νx → νs conversion only or (b) collisional νs production only. The solid and dotted orange curves correspond to different assumptions on the maximum allowed energy loss, expressed here in terms of the ratio R of the energy output in sterile neutrinos and the total energy output. For comparison, we also show in green the parameters for which the Dodelson-Widrow mechanism [3] predicts the correct relic abundance of νs.
neutrinos up to this time. Doing so modifies our limit by 30%.
Our limits on the ν s parameter space are shown in fig. 3 for adiabatically and collisionally produced ν s separately, and in fig. 1 for the combination of both production mechanisms (requiring R < 1). Thanks to MSW enhancement, our constraints reach down to sin 2 2θ ∼ 10 −14 at m s ∼ 10-100 keV, surpassing all other limits in this mass range.
The shape of the exclusion regions in fig. 3 can be understood as follows. For adiabatic conversion at small m s , the oscillation length at the MSW resonance L res osc is large, making flavor conversion non-adiabatic according to eq. (6). At large m s 100 keV, the resonance condition of eq. (2) cannot be satistifed. At m s ∼ few × 10 keV, adiabatic flavor conversion is most effective at low sin 2 2θ ∼ 10 −14 . For larger sin 2 2θ, the radial width R width of the resonance region becomes too large, so that neutrinos scatter before having a chance to convert. Also collisional production is most effective when the mixing angle in matter is MSW-enhanced. At m s 10 keV, the resonance condition of eq. (2) is fulfilled only in the outer regions of the supernova, where the scattering rate is too low for effective collisional production. At m s 100 keV, the resonance condition is never satisfied.
Thus, if ν s are abundantly produced in a supernova, we expect the explosion to be accompanied by a flux of energetic (O(1-100 MeV)) secondary gamma rays. Photons in this energy range are not normally expected from a supernova or supernova remnant, both of which emit Xrays only at energies 10 keV [51,52]. The arrival times of the gamma rays from ν s decay are spread out over a time interval ∆t 3.6 hrs × d 1 kpc m s 1 keV where d is the distance to the supernova and E γ is the gamma ray energy (which is just half the ν s energy, E γ = E ν /2). Sterile neutrinos decaying immediately after their production lead to gamma rays that reach the Earth at the same time as the active neutrino burst. Gamma rays from ν s decaying only after travelling a distance d are delayed by ∆t.
The photon flux from ν s decay in units of cm −2 sec −1 MeV −1 is given by where γ = 2E γ /m s is the Lorentz boost factor. The factor in parenthesis accounts for the ν s decay probability. Note that the sterile neutrino spectrum dN νs (E ν )/dE ν is evaluated at E ν = 2E γ . We plot dφ γ (E γ )/dE γ for two sets of benchmark parameters in fig. 4, and the total gamma ray flux as a function of sin 2 2θ and m s is shown in fig. 5. In both figures, we have chosen d = 1 kpc. As expected from  fig. 1 (orange region) and to constraints from X-ray searches [10-13, 25, 26] (green, dark blue, gray, and light blue regions).
fig. 2, and taking into account E γ = E ν /2, the photon spectrum peaks at E ∼ 10-100 MeV. Comparing the expected photon spectrum to the diffuse astrophysical gamma ray background measured by COMPTEL [53][54][55], EGRET [56][57][58] and Fermi-LAT [59], we find that the signal is still several orders of magnitude below the uncertainty on the background. Most likely, a direct observation of the ν s -induced photon signal is therefore beyond the scope of even the next generation of Compton telescopes [61][62][63] and would require a factor ∼ 10 5 improvement compaired to the projected sensitivity of ComPair [62]. Summary. In conclusion, we have computed the flux of hypothetical keV-MeV sterile neutrinos ν s from a supernova. We have constrained the ν s parameter space using energy loss arguments ( fig. 1), and we have estimated the gamma ray flux from ν s produced in a nearby supernova. Directions for future work include a more detailed calculation of the ν s flux from a supernova, going beyond the two flavor approximation and with a more detailed treatment of collective effects. We moreover plan to investigate how our results depend on the mass of the progenitor star [64].