Constraints on sterile neutrino models from strong gravitational lensing, Milky Way satellites, and Lyman-$\alpha$ forest

The nature of dark matter is one of the most important unsolved questions in science. Some dark matter candidates do not have sufficient nongravitational interactions to be probed in laboratory or accelerator experiments. It is thus important to develop astrophysical probes which can constrain or lead to a discovery of such candidates. We illustrate this using state-of-the-art measurements of strong gravitationally-lensed quasars to constrain four of the most popular sterile neutrino models, and also report the constraints for other independent methods that are comparable in procedure. First, we derive effective relations to describe the correspondence between the mass of a thermal relic warm dark matter particle and the mass of sterile neutrinos produced via Higgs decay and GUT-scale scenarios, in terms of large-scale structure and galaxy formation astrophysical effects. Second, we show that sterile neutrinos produced through the Higgs decay mechanism are allowed only for mass $>26$ keV, and GUT-scale scenario $>5.3$ keV. Third, we show that the single sterile neutrino model produced through active neutrino oscillations is allowed for mass $>92$ keV, and the 3 sterile neutrino minimal standard model ($\nu$MSM) for mass $>16$ keV. These are the most stringent experimental limits on these models.


INTRODUCTION
The nature of dark matter (DM) is one of the most important questions in modern physics, with implications spanning from particle physics to astrophysics and cosmology. This unknown particle contributes 25% of the total energy of the universe [1], but is not made of ordinary matter and has no electromagnetic interaction.
Many DM models have been proposed. A number of candidates fall into the class of cold dark matter (CDM) [2], made of collisionless particles considered "cold" due to their small velocity dispersion relative to the speed of light. This model is extremely successful on supergalactic scales but there are open challenges at subgalactic scales [3]. CDM predicts more satellites than are observed around galaxies of Milky Way (MW) mass, "cuspy" dark matter density profiles in contrast to the flatter cores observed in dwarf galaxies and clusters, and predicts that subhalos hosting the largest MW satellites are either underdense or too small. It is still unclear whether these challenges can be solved by a better understanding of baryonic processes, or whether alternative dark matter models are needed [e.g. 4, and references therein].
Plenty of DM models have been proposed to eliminate these small-scale tensions between observations and CDM [5]. DM particles which are generated with higher velocity dispersions erase fluctuations in the matter power spectrum at scales smaller than a characteristic 'free-streaming length', suppressing structures below this scale. So-called 'hot' dark matter candidates such as standard neutrinos are ruled out by observations [6][7][8][9], as the main DM component. However, a broad range of "warm" DM (WDM) with smaller but non-negligible free streaming lengths are viable. One popular class of WDM models are sterile neutrinos (SNs).
Sterile neutrinos [10][11][12][13][14] are particles with righthanded chirality, no charge, and no color charge, and therefore do not interact with standard model particles except via mixing with neutrinos, or via some nonstandard-model interactions. They were first introduced for the purpose of explaining the masses of active (lefthanded) neutrinos, and can have masses in the range from eV to the Planck scale. In the early universe they can decouple from the plasma before electron-positron annihilation, when they are still relativistic [12,[15][16][17][18].
The exact production mechanism for a given SN model determines the clustering of dark matter [19][20][21][22]. Dodelson and Widrow [23] (hereafter the DW) model adds a single SN with coupling to active neutrinos. DW SNs are produced in neutrino oscillations at temperatures below a few GeV [24].
Alternatively, there may be multiple SN species produced in neutrino oscillations. In the seesaw theory of neutrino masses, one needs at least two right-handed states [25], but more right-handed states are allowed. In the presence of a sizeable lepton asymmetry, the active to SNs conversions select a lower-momentum part of the thermal distribution, leading to somewhat colder dark matter [Shi and Fuller (SF),26]. Thus, the popular neutrino minimal standard model [26, 27, hereafter νMSM or SF] postulates three right-handed neutrinos with masses below the electroweak scale.
SNs could also be produced through mechanisms other than active neutrino oscillations, including "freezein" production from decays of the inflaton [28] or an SU(2)×U(1) singlet Higgs boson [PK, [19][20][21][22]. Most of the SNs production from oscillations (DW or SF) takes place at temperature ∼ 0.1 GeV. In contrast, Higgs boson decays can produce a population of SNs at a temperature ∼100 GeV. Subsequent cooling and entropy production dilutes and redshifts this population, making the resulting DM colder than the WDM produced by DW or SF. This model, which produces particles in the 1-10 keV range, has been shown to produce the correct dark matter abundance, resulting in the so-called "keV Miracle Model" (hereafter PK).
Another possible production mechanism for SNs is the split seesaw mechanism [29, hereafter KTY]. The model predicts two large and one small Majorana masses due to a natural separation of scales. The large Majorana masses allow for thermal leptogenesis, while the keV mass produces a DM candidate. The model can be embedded into an SO(10) Grand Unified Theory (GUT), or some other theory containing a gauge U(1) B−L symmetry. The resulting DM is colder than the models described above, due to dilution and redshifting of SNs as the plasma cools from the GUT-scale production temperature. The 4 SN models described above could explain the unidentified 3.5 keV X-ray line found in observations of galaxies and clusters [30,31], even if they do not account for 100% of DM. Furthermore, keV SNs can have a dramatic effect on supernovae, e.g., explaining the pulsar kick velocities in excess of 1000 km/s which so far has evaded other explanations [32,33].
We use state-of-the-art measurements of gravitationally lensed quasars [34,35], MW satellites [36][37][38][39], and Lyα [40,41] , to constrain the 4 popular SN models described above. We take their limits in terms of the thermal relic WDM, and compute the equivalent limits for the four SN models described above, from the point of view of cosmological structure formation and the halo mass function. As we will show, our analysis provides the most stringent limits to date on these four SN models. the GUT-scale scenario (KTY), the 'keV miracle model' Higgs production mechanism (PK), and single particle neutrino oscillation production mechanism (DW). The posteriors do not go to 0 on the right limit, so we cannot impose upper constraints on the particle masses; however, since they do go to 0 on the small limit, we can derive a lower limit. The vertical dashed line marks the 95% lower boundary interval, corresponding to 4.6/2.1/11/34 keV for thWDM/KTY/PK/DW. The limits for the νMSM model depend on lepton asymmetry, and are discussed in the text. The case where the assumption for the average background density of the universe includes (Case I), or does not include (Case II), baryonic matter in addition to dark matter (see Appendix) is shown.

THERMAL RELIC WARM DARK MATTER CONSTRAINTS FROM STRONG GRAVITATIONAL LENSING
Strong gravitational lensing depends only upon gravity and is thus sensitive to the abundance of halos irrespective of their ability to emit or absorb light.It can thus determine the halo mass function directly, avoiding uncertainties related to the physics of star formation in low mass galaxies that affect traditional methods.
[35] used eight quadruply-imaged quasars systems to constrain the amplitude of the halo mass function and the free-streaming length of dark matter. For each system, many realizations of dark matter structure are drawn from analytic dark matter halo mass functions flexible enough to describe mass functions produced by a broad range of thermal relic dark matter masses.
The predicted flux ratios for each realization are compared with the observed flux ratios to estimate the likelihood using Approximate Bayesian Computing. The likelihoods from each system are multiplied together to infer the parameters common to all systems. A key parameter is the 'half-model mass' m hm which is the mass at which there are half as many halos as there would have been in the case of CDM.
We marginalize over all the other parameters to obtain  the posterior distribution for m hm . As always in Bayesian statistics, the posterior depends on the choice of priors.
Since we do not know the order of magnitude of the SN mass we adopt a uniform prior in log 10(m hm ), within the range 10 4.8 − 10 10 M (reported in units of Sun mass following existing literature on the halo mass function). The constraints on m hm can be tied into constraints on the mass of the DM particle, given a model. A traditional reference model is the thermal relic WDM. It does not refer to a physical particle in particular, but serves as a standard tie-in model for the properties of WDM models with thermal relics. [42], [43] derived a one-to-one mapping between the half-mode mass and the mass of the thermal relic WDM. The general form of the conversion is where A has two possible values (Case I and II), depending on assumptions about the background density of the universe, as detailed in the Appendix. Using Eq. 1, we obtain the posterior shown in Fig. 1.

RELATION BETWEEN THERMAL RELIC WDM AND STERILE NEUTRINO TRANSFER FUNCTIONS
Our goal is to use the constraints discussed as an illustration in §, and those obtained by [39] and [40,41] to derive constraints on the four SN dark matter candidates discussed in §: the GUT-scale scenario (KTY), the 'keV miracle model' Higgs production mechanism (PK),  TABLE I. Coefficients for the fits for the relation between the msn and m thWDM , for the cases of the Higgs production mechanism (PK) and GUT scale (KTY). The power law fit, as well as the first 2 degree polynomials are shown.
the single particle neutrino oscillation production mechanism (DW), and the Shi-Fuller mechanism within the neutrino minimal standard model (νMSM).
The key quantity is the transfer function, T , which describes the effect of free-streaming on matter distribution. Given the power spectrum of initial density fluctuations P i , T describes its evolution as a function of scale k and cosmic time, with respect to a standard CDM model: Transfer functions for SNs (Fig. 2) for the Higgs production mechanism and the GUT-scale scenarios have been previously obtained by [22]. We calculate the transfer functions of the several models by using the momentum-space distribution functions as tables that are provided to CLASS. Including these with the proper effective temperature of the dark matter models allows for accurate calculation of the transfer functions relative to CDM. The ones given here supersede the published ones by adopting more up to date cosmological parameters. We also corrected a mismatch between expected and provided dilution factors, using the non-CDM features of CLASS [44]. For the initial conditions, we fix the cosmology to the mean results from [1] and the CMB temperature from [45].
The thermal relic WDM transfer functions for a given m thWDM can be obtained from the analytical form presented in equation A9 of [42] (with more recent numbers from [46]), and it is given by with α, µ given in the Appendix. We tested two methods to find the relation between SN and thermal relic WDM transfer functions: first, we fit the SN transfer functions using the mass of the thermal relic WDM as a free parameter (Eq. 2); second, we match the 'half-mode' wavenumber where the transfer functions decrease to 0.5. The half-mode k hm and the mass of the particle m thWDM can be related analytically from Eq. 2, to obtain: The masses of the thermal relic WDM particles and SN particles corresponding to the PK transfer functions shown in Fig. 2, as well as those for the KTY model, are shown as scattered points. Polynomial and power-law fits are shown as solid lines. These relations allow us to map constraints on the mass of thermal relic WDM to the corresponding mass of a SN, for the Higgs production mechanism (PK), and the GUT-scale scenario (KTY).
In the first method, the fits might depend too much on how the numerical T (k) for the SNs was sampled, and the errors that arise at high k. The 1/2 mode matching gets rid of those issues, so we use this one for the rest of the paper. The difference in the end results of the two methods is of order 1%. We obtain the fits shown in Fig.  2. For root-mean-square goodness of fit estimator, the PK model KTY model were comparable.
The relation between the mass of the SN which generates an equivalent transfer function to a thermal relic particle with a given mass can be approximated by a polynomial of degree deg, m sn = f (m thWDM ) = deg 0 a i · m i thWDM , or by a power law m sn = a · m b thWDM . Results of the polynomials fits of orders 1, 2 (higher orders resulted in over fitting), as well as the power law, can be seen in Table I. Using these coefficients, a relation can also be derived between m sn and m hm , seen below in the case of the power law fit: where A takes the values described in the Appendix. These relations are calibrated on the m thWDM mass interval of (0. 75,22) keV, and then extended to 60 keV when applied to data. Fig. 3 shows the results of these fits. We use the results of the second order polynomial but we note that using a linear approximation would change the inferred bounds on SN mass by only 11%.
This mapping allows one to convert any results obtained for thermal relic WDM to SN in post-processing, without redoing the experiment or the analysis.
For SNs produced through the oscillation mechanism for the Dodelson-Widrow model, the relation between m sn and m thWDM is taken from [46]. For the νMSM model, [47] derives the model connections to the half-mode mass m hm , and we use the m hm posterior to constrain them.

MAPPING THERMAL RELIC WARM DARK MATTER CONSTRAINTS ONTO STERILE NEUTRINOS
We obtain the most stringent experimental limits on four SN models: PK, KTY, νMSM, and DW.
To recap, our starting point are the following 95% confidence limits on thermal relic WDM. The lensingonly analysis described in § gives m thWDM > 4.6keV [35]. Combination with satellite counts extends it to m thWDM > 9.7 keV [39]. Independent work on the Lyman-α forest yields 3.3 keV [40] and 5.3 keV [41], the latter using additional assumptions for the relevant thermodynamics.
For the PK and KTY models, the relations derived in § can now be used together with the ones for νMSM and DW to translate limits from thermal relic WDM into limits for the SNs masses (Fig. 1). The 95% limits for the four models are given in Tab. II.
We note that the posterior shown in Fig. 1 vanishes at the lower bound but not on the upper bound, as expected because the warmest models are ruled out by a number of observations [8,9,39]. This puts a lot of weight on the choice of priors, which can influence the limits reported. We thus convert the posteriors to lower limits, although of course the full posterior is more informative. Likelihood ratios can be obtained from Fig. 1.
For the νMSM model, we use the posterior on m hm to eliminate the model space as shown in Fig. 2 of [47] , which presents the expected half-mode mass as a function of lepton asymmetry (L6) for different neutrino masses. Our upper limit on log 10 (m hm [M ]) from strong lensing alone is 8.1. This rules out masses under 7.0 keV for all lepton asymmetries. For higher masses, only limited ranges of lepton asymmetries are allowed: 7keV: L6 ∈ (6.8, 7.6), 9keV: L6 ∈ (5.2,7.8), 11keV:L6 ∈ (4.3, 7.6), 14keV: L6 ∈ (1.7, 7.9), 16 keV: L6 ∈ (1.6, 11.5). After incorporating the MW satellite counts constraints, log 10 (m hm [M ]) > 7.0, which corresponds to m νM SM >16keV. These limits are improved compared to existing work ( [48]). These results are contingent on the assumption that DM is made from a single component, the SN model of choice. However, DM could be a mixture of different components, such as the "mixed cold+warm DM model" [49,50]. Adding a parameter to TABLE II. 95% lower limits for four sterile neutrino models: Higgs production mechanism (PK), GUT scale scenario (KTY), 3 sterile neutrino minimal standard model (νMSM), single model sterile neutrino produced through active neutrino oscillations (DW). The limits are derived from 4 datasets: gravitational strong lensing [35], strong lensing combined with Milky Way galaxy counts [39], Lyman-α forest [40], and Lyman-α forest combined with thermodynamic assumptions [41]. The case I and case II labels correspond to different assumption cases about the average background density of the universe, as described in the Appendix. In the last two rows of the table, we also presents the limits in terms of the half-mode mass, and the thermal relic WDM mass.
control the abundance ratio could make the constraints weaker. Future work would aim to combine the limits from strong lensing, galaxy counts, and Lyman-α forest in a joint analysis. Work has already been done in this direction ( [51]), however their dataset obtained less stringent limits than the strong lensing combined with galaxy counts obtained by [39]. Future analysis combining the data sets used in Table II will be useful. In addition, work exploring analytical connections in the non-linear regime (as was done at the level of the transfer function by [52]) at the sub-halo level would be useful. CONCLUSION We used flux ratios of strong gravitationally-lensed quasars, MW satellites, and Lyman-α forest to constrain four of the most popular SN models.
First, we derive effective relations to describe the correspondence between the mass of a thermal relic WDM particle and the mass of SNs produced via Higgs decay and GUT-scale scenarios, in terms of astrophysical effects. We take advantage of the similarity between the transfer functions of the SNs mechanism presented by [22], to that of thermal relic WDM.
We note that our derived equivalence relations are of general importance, and can be used to put limits on SN models for any thermal relic WDM measurement, not just the ones we present here.
The limits on the PK, KTY, νMSM and DW models summarized in Table II are the most stringent experimental limits on these four models. We note that the limits from lensing and MW satellites are independent of and agree with those from the Lyα forest. We have effectively ruled out part of the parameter space for SN generated through these 4 models.
IZ and TT acknowledge support by the National Science Foundation grant NSF-1836016, by the Gordon and Betty Moore Foundation Grant 8548 and by NASA grant HST-GO-15177. Part of the data used in this paper were obtained as part of HST- GO-15177 [58], NumPy [59], Python [60,61], scikit-learn [62] Appendix The thermal relic WDM transfer functions can be approximated by an analytical function with a dependence on m thWDM , as show in equation A8 of [42], and Eq. 2. We use the more recent fit to the characteristic length scale factor α(m thWDM ) obtained by [46,63], with µ = 1.12:  (5) in units of Mpc h −1 . By convention, the number of degrees of freedom is taken to be 1.5 for the warm particle, based on the equivalent contribution of a light neutrino species [42]. Following [43], we assume that the characteristic length scale α can be related to an effective free-streaming length scale λ eff fs . The 'half-mode' length scale λ hm corresponds to the scale at which the thermal relic WDM transfer function (2) decreases to 1/2: λ hm = 2πλ eff fs (2 µ/5 − 1) − 1 2µ .
We define the corresponding 'half-mode' mass scale and obtain a relation between m hm and m thWDM : whereρ is the background density of the universe. The A constant thus depends on the assumedρ. We consider two cases, where the background density of the universe can be taken as: where Ω M includes the contributions of both baryonic and dark matter (as done in [39], or II.ρ = Ω DM × ρ critic , where only the contributions of dark matter are included (as done in [35]).
We thus obtain two equations between the 'half-mode' mass and the thermal relic mass: I. m hm = 3.81 · 10 8 m thWDM 3.3keV