Constraints from Ly-$\alpha$ forests on non-thermal dark matter including resonantly-produced sterile neutrinos

TL;DR

This paper uses Ly-$\\alpha$ forests from BOSS data to constrain non-thermal dark matter models, including sterile neutrinos, by analyzing flux power spectra and hydrodynamical simulations, tightening bounds on particle masses and abundances.

Contribution

It provides new constraints on sterile neutrino and warm dark matter models using Ly-$\alpha$ forest data and hydrodynamical simulations, including effects of lepton asymmetry and higher-resolution datasets.

Findings

Sterile neutrinos as light as 6.5 keV are consistent with BOSS data at 2 sigma.

Constraints on warm dark matter relic mass follow F ~ 0.35 (keV/m_x)^-1.37.

A 7 keV sterile neutrino with specific lepton asymmetry is compatible at 1.9 sigma with BOSS and at 3.1 sigma with extended data.

Abstract

We use BOSS DR9 quasars to constrain 2 cases of dark matter models: cold-plus-warm (C+WDM) where the warm component is a thermal relic, and sterile neutrinos resonantly produced in the presence of a lepton asymmetry (RPSN). We establish constraints on the relic mass m_x and its relative abundance $F=\Omega_{wdm}/\Omega_{dm}$ using a suite of hydrodynamical simulations in 28 C+WDM configurations. We find that the 3 sigma bounds approximately follow F ~ $0.35 (keV/m_x)^{-1.37}$ from BOSS alone. We also establish constraints on sterile neutrino mass and mixing angle by producing the non-linear flux power spectrum of 8 RPSN models, where the input linear power spectrum is computed directly from the particles distribution functions. We find values of lepton asymmetries for which sterile neutrinos as light as 6.5 keV (resp. 3.5 keV) are consistent with BOSS at the 2 sigma (resp. 3sigma) level. These limits tighten by close to a factor of 2 for lepton asymmetries departing from those yielding the coolest distribution functions. Our Ly-a forest bounds can be strengthened if we include higher-resolution data from XQ-100, HIRES and MIKE. At these smaller scales, the flux power spectrum exhibits a suppression that can be due to Doppler broadening, IGM pressure smoothing or free-streaming of WDM particles. In the current work, we show that if one extrapolates temperatures from lower redshifts via broken power laws in T_0 and gamma, then our 3 sigma C+WDM bounds strengthen to F~ $0.20 (keV/m_x)^{-1.37}$, and the lightest RPSN consistent with our extended data set have masses of 7.0 keV at the 3 sigma level. Using dedicated hydrodynamical simulations, we show that a 7 keV sterile neutrino produced in a lepton asymmetry $L = 8 \times 10^{-6}$ is consistent at 1.9 sigma (resp. 3.1 sigma) with BOSS (resp. BOSS + higher-resolution), for the thermal history models tested in this work.

Tables

Table 1. Table 1 : Hydrodynamical simulation parameter grid. Our best-guess is run on the central values. First and second order derivatives are computed using the stepped values in the right column.

parameter	central value & step
$H_0$	$\left(0.675\pm 0.05\right)\times {10}^2\mathrm{km}{s}^{-1}{\mathrm{Mpc}}^{-1}$
${\Omega }_m$	$0.31\pm 0.05$
$n_s$	$0.96\pm 0.05$
${\sigma }_8$	$0.83\pm 0.05$
$T_0^{z=3}$	$\left(1.4\pm 0.7\right)\times {10}^{4}K$
${\gamma }^{z=3}$	$1.3\pm 0.3$
$A^{\tau }$	$\left(2.5\pm 2.0\right)\times {10}^{-3}$
${\eta }^{\tau }$	$3.7\pm 0.4$

Table 2. Table 2 : Summary of Ly- α 𝛼 \alpha constraints on NRP neutrino mass according to the data set used. Quoted lower bounds are the 95% confidence level. Tension with the upper bound from X-rays is expressed in standard deviations in the right-most column. Simulation resolution refers to the quantity N L 𝑁 𝐿 \frac{N}{L} in Mpc − 1 superscript Mpc 1 \rm{Mpc^{-1}} .

Reference	QSO spectra	Data	Simultion	Bounds on $m_{{\nu }_s}^{\mathrm{nrp}}$	Tension with
study	resolution	set	resolution	from Ly-$\alpha$ forests	X-ray bounds
[19]	low only	SDSS-I	12.8 (hydro)	$\ge 12$ keV	$6\sigma$
[22]	low + high	SDSS-I + UVES	6.7 (N-body)	$\ge 10$ keV	$5\sigma$
[21]	high only	HIRES + MIKE	25.6 (hydro)	$\ge 19$ keV	$9.5\sigma$
[25]	low only	SDSS-III	30.72 (hydro)	$\ge 24$ keV	$12\sigma$
[44]	low + medium	SDSS-III + XQ100	30.72 (hydro)	$\ge 25$ keV	$12.5\sigma$
[75]	medium + high	XQ100 + MIKE	38.4 (hydro)	$\ge 34$ keV	$17\sigma$

Table 3. Table 3 : Best-fit values and 68% C.L. of the most relevant parameters. We take 1 ​ keV / m x = 0.05 1 keV subscript 𝑚 𝑥 0.05 \mathrm{1\,keV}/m_{x}=0.05 , F wdm = 0.30 subscript 𝐹 wdm 0.30 F_{\rm{wdm}}=0.30 as the best-fit C+WDM model, although values for the parameters of the table change by < 1 ​ σ absent 1 𝜎 <1\;\sigma along the degeneracy F wdm = α ​ ( 1 ​ keV / m x ) β subscript 𝐹 wdm 𝛼 superscript 1 keV subscript 𝑚 𝑥 𝛽 F_{\rm{wdm}}=\alpha(1\,{\rm keV}/m_{x})^{\beta} curve as well as from a pure WDM fit.

Parameter	SDSS	SDSS+XQ+HR
${\sigma }_8$	$0.855\pm 0.020$	$0.815\pm 0.020$
$n_s$	$0.935\pm 0.010$	$0.950\pm 0.010$
$T_0\left(z=3\right)$ (K)	$9600\pm 4000$	$14500\pm 3000$
$\gamma$	$0.9\pm 0.2$	$0.9\pm 0.2$
	$-2.7\pm 0.7$	$-1.9\pm 0.4$
${\eta }^{T_0}\left(z>3\right)$	$-4.1\pm 1.4$	$-2.1\pm 0.6$
${\eta }^{\gamma }$	$0.7\pm 0.5$	$-0.3\pm 0.4$

Table 4. Table 4 : Tension with respect to measured P φ ​ ( k ) subscript 𝑃 𝜑 𝑘 P_{\varphi}(k) for the hydrodynamical simulations that fall in the range between 1 and 4  σ 𝜎 \sigma . RPSN configuration is given in column 1, comparison to ‘SDSS only’ in column 2 and to ‘SDSS+XQ+HR’ in column 3.

RPSN model	SDSS	SDSS+XQ+HR
M4L12⋆	$3.06\sigma$	$>4\sigma$
M6L6	$3.13\sigma$	$>4\sigma$
M6L9⋆	$2.0\sigma$	$3.3\sigma$
M7L8⋆	$1.9\sigma$	$3.1\sigma$
M8L4	$2.7\sigma$	$>4\sigma$
M8L8⋆	$1.5\sigma$	$2.5\sigma$

Equations

$$n_{\text{textsc}\mathrm{Hi}}=\left(\frac{{\sigma }_{\mathrm{comb}}\sqrt{2k_b/m_e}}{c{\sigma }_{\mathrm{ion}}}\right)\frac{n_b^2\sqrt{T}}{n_{{\gamma }_{\mathrm{UV}}}}$$

$${\phi }_{\mathrm{obs}}\left(\lambda \right)=e^{-\tau (\lambda )}{\phi }_{\mathrm{qso}}\left(\lambda \right)$$

$$\tau \left(v\right)={\int }_0^v\mathrm{d}{v}_{\parallel }^{\prime }\frac{n_{{\mathrm{H}}_{\mathrm{I}}}{\sigma }_{\mathrm{Ly}\alpha }\left(z\right)}{\nabla v_{\parallel }^{\prime }}$$

$${\delta }_{\phi }(\lambda )=\frac{\phi (\lambda )-⟨\phi ⟩}{⟨\phi ⟩}=\frac{e^{-\tau }}{e^{-{\tau }_{\mathrm{eff}}}}-1$$

$$P_{\mathrm{Ly}\alpha }(k)={⟨\frac{|\widetilde{{\delta }_{\phi }}(k){|}^2-P_{\mathrm{noise}}(k)}{W^2(k)}⟩}_{QSOs}$$

$$T(\delta ,z)=T_0(z)\times (1+\delta {)}^{\gamma (z)-1}$$

$$-\ln ⟨\phi ⟩={\tau }_{\mathrm{eff}}=A^{\tau }\times (1+z{)}^{{\eta }^{\tau }}$$

$$P_{\mathrm{cdm}}(k,z)={\mathcal{T}}^2(k,z)\times \mathcal{P}(k)\times {\mathcal{D}}^2(z)$$

$$\mathcal{P}(k)=\frac{2{\pi }^2}{k^3}{\mathcal{A}}_s{\left(\frac{k}{k_{\star }}\right)}^{n_s-1}$$

$$2\pi P^{1d}(k)={\int }_k^{\infty }d{k}^{\prime }{k}^{\prime }{P}^{3d}(k^{\prime })$$

$${\lambda }_{\mathrm{FSH}}(t)=a(t){\int }_0^{a(t)}da\frac{⟨v⟩}{a^{2}H}$$

$$\frac{{\lambda }_{\mathrm{FSH}}(t)}{a(t)}=\frac{2\pi }{k_{\mathrm{FSH}}(t)}$$

$$\frac{{\lambda }_{\mathrm{FSH}}^0}{a_0}=\frac{2\pi }{k_{\mathrm{FSH}}^0}\simeq {\int }_0^{a_{\mathrm{nr}}}\frac{da}{a^{2}H}+{\int }_{a_{\mathrm{nr}}}^{a_0}\frac{a_{\mathrm{nr}}da}{a^{3}H},$$

$${\Omega }_{\mathrm{wdm}}=F_{\mathrm{wdm}}\times {\Omega }_{\mathrm{dm}}=F_{\mathrm{wdm}}\times \left({\Omega }_{\mathrm{wdm}}+{\Omega }_{\mathrm{cdm}}\right)$$

$$m_{{\nu }_s}=\kappa m_x^{\mu }/{\omega }_{\mathrm{wdm}}^{1/3}$$

$$\mathrm{CL}({\theta }_i,{\theta }_j,...,{\theta }_n)=1-{\int }_{\Delta {\chi }^2({\theta }_i,{\theta }_j,...,{\theta }_n)}^{\infty }\mathrm{d}x{f}_{N_{\mathrm{dof}}}(x)$$

$$f_{N_{\mathrm{dof}}}(x)=\frac{e^{-x/2}{x}^{\frac{N_{\mathrm{dof}}}{2}-1}}{\sqrt{2^{N_{\mathrm{dof}}}}\Gamma (N_{\mathrm{dof}}/2)}$$

Full text

Our Lyman-$\alpha$ forest bounds can be additionally strengthened if we include higher-resolution data from XQ-100, HIRES and MIKE that allow us to probe smaller scales. At these scales, the measured flux power spectrum exhibits a suppression that can be due to Doppler broadening, IGM pressure smoothing or free-streaming of WDM particles. In order to distinguish between these mechanisms, thermal history at redshifts $z\geq 5$ should be determined. In the current work, we show that if one extrapolates temperatures from lower redshifts via broken power laws in $T_{0}$ and $\gamma$, then our $3\sigma$ C+WDM bounds strengthen to $F_{\rm{wdm}}\sim 0.20(\mathrm{keV}/m_{x})^{-1.37}$, and the lightest resonantly-produced sterile neutrinos consistent with our extended data set have masses of $\sim 7.0~{}\mathrm{keV}$ at the $3\sigma$ level. In particular, using dedicated hydrodynamical simulations, we show that a hypothetical 7 keV sterile neutrino produced in a lepton asymmetry of $\mathcal{L}=|n_{\nu_{e}}-n_{\bar{\nu}_{e}}|/s~{}=~{}8\times 10^{-6}$ is consistent at $1.9~{}\sigma$ (resp. $3.1~{}\sigma$) with BOSS (resp. BOSS + higher-resolution) data, for the thermal history models tested in this work. More information about the state of the IGM at redshifts 5–6 will allow one to conclude whether the small-scale suppression of the flux power spectrum is due to such sterile neutrino or to thermal effects.

1 Introduction

It appears with increasingly robust evidence from gravitational lensing, non-Newtonian galaxy rotation curves, baryon acoustic oscillations (BAO) and the cosmic microwave background (CMB) that 80% of all gravitating matter in the Universe is dark. Even under the assumption that dark matter (DM) is an additional component to baryonic matter, the standard model of particle physics currently has no viable DM candidate, i.e., a massive, neutral, long-lived particle insensitive to the electromagnetic interaction. The discovery of neutrino flavor oscillations (see, e.g., [1]) proved our current picture of the standard model to be incomplete by providing direct evidence for non-zero neutrino masses. Neutrinos are, however, too light to account for all of dark matter. Therefore, a particle origin of dark matter necessarily implies extensions of the standard model of particle physics, many of which provide dark matter candidates (for a review see, e.g., [2, 3]). One such candidate, similar to neutrino dark matter while being free of its drawbacks, is the right-handed (a.k.a. sterile) neutrino [4]. Right-handed neutrinos carry no lepton charge and only directly couple (mix at quadratic level) to left-handed (active) neutrinos. As such, sterile neutrinos never thermalize in the early Universe. Therefore their number density is much lower than that of the cosmic neutrino background. Nevertheless, enough may be produced to account for the observed total matter density, and are thus an inherent DM candidate [4, 5, 6, 7, 8, 9]. Although produced while relativistic, these particles become non-relativistic deeply in the radiation dominated epoch (RDE) — as any dark matter should — and therefore evade the cosmic microwave background constraint on the effective number of stable, relativistic species in thermal equilibrium in the early Universe, $N_{\rm{eff}}=2.99\pm 0.20$ (from the combined TT+TE+EE+lowP Planck analysis [10]).

Such particles, becoming non-relativistic deep in the RDE, are generically known as warm dark matter (WDM). WDM particles free-stream out of gravitational potential wells, delaying structure growth on scales below a characteristic distance $\lambda_{\rm{FSH}}$ that scales with the expansion rate $H=\dot{a}/a$. WDM is expected to produce fewer dark matter halo satellites than CDM [11]. However, there is still no direct way to probe the dark matter distribution in the Local Group other than by gravitational lensing. The fewer-than-expected galactic satellite count around the Milky Way does not constitute conclusive evidence in favor of a warm dark matter since there are several well understood mechanisms which lead to poorly-populated dark matter haloes such as re-ionization and feedbacks from active galactic nuclei (AGNs) and supernovæ [12]. Incorporating these and other baryonic effects on sub-Mpc scales in N-body simulations can also account for the observed less cuspy inner density profiles of galaxy halos. Distinguishing cold from non-cold dark matter (NCDM) requires either directly probing the dark matter distribution in galaxy clusters such as with strong gravitational lensing [13, 14, 15], or probing the scales at which WDM is expected to produce no structures beyond a certain mass. Lyman-$\alpha$ (Ly-$\alpha$) forests in the spectra of distant $z\in\left[2.1,4.5\right]$ quasars (QSO) probe the intergalactic $\sim 10^{0-2}~{}h^{-1}~{}\rm{Mpc}$ scales encompassed by the integrated free-streaming of keV particles while they remained relativistic. They are thus a formidable probe for establishing lower bounds on NCDM particle mass [16, 17, 18, 19, 20, 21, 22, 23, 24, 25].

Since they mix with left-handed neutrinos, sterile neutrinos can be copiously produced in the early Universe  [4, 6, 7, 8]. A net lepton asymmetry (an excess of leptons over anti-leptons) present in the primordial plasma during the production epoch boosts active-sterile neutrino mixing and thus can account for the observed abundance of dark matter [5, 7, 9]. This requires lepton asymmetries of several orders of magnitude larger than the baryon asymmetry of the Universe [9] estimated from light element abundances from big bang nucleosynthesis. These resonantly-produced sterile neutrinos (RPSN) are relevant in that they have much cooler distribution functions than is assumed in the so-called non-resonant production (NRP) mechanism, which occurs in the absence of lepton asymmetry. Through their cooler distribution functions on the one hand and their smaller mixing angles on the other hand, a RPSN of mass $m_{\nu_{s}}^{\rm{rp}}$ is more compatible with clustering, Ly-$\alpha$ and X-ray data than a NRP neutrino of the same mass $m_{\nu_{s}}^{\rm{nrp}}$ [23]. Additional production mechanisms exist (see [26]) including production from the decay of a scalar field [27, 28, 29] and diluted thermal overproduction [30, 31] to name just a few. They all have in common a cooler momentum distribution than NRP neutrinos which translate into shorter free-streaming horizon scales.

In this work we investigate several NCDM dark matter model, different from thermal relic WDM. Only a handful of similar studies have been conducted thus far. Ref. [22] issued bounds on the cold+warm dark matter mixtures, using SDSS-II data. Ref. [23] have issued bounds on the mass of RPSN by rescaling the Ly-$\alpha$ bounds on the mixture of cold+warm dark matter [22]. Several other works [24, 32, 33] used the bounds on thermal relic WDM and extrapolated the expected flux power spectrum in the Ly-$\alpha$ forest using a correction from the linear matter power spectrum (in [33], the matching was also performed at the level of the non-linear power spectrum obtained from DM-only simulations).

The current work issues bounds on the mass of RPSN using dedicated hydrodynamical simulations, a approach which had not been implemented directly before for Ly-$\alpha$ forest constraints. The paper is outlined as follows. Sec. 2 recaps how Ly-$\alpha$ forests are used to probe the $10^{0-2}~{}h^{-1}~{}\rm{Mpc}$ scales through the flux power spectrum with the help of state-of-the-art hydrodynamical simulations. We mention how NCDM cosmologies are incorporated in those simulations in Sec. 2.2. We present the Ly-$\alpha$ forest power spectrum constructed from 3 samples of QSOs, from the Sloan Digital Sky Survey (SDSS), the Very Large Telescope (VLT) and W. M. Keck Observatory, that enable us to constrain NCDM cosmologies in Sec. 2.1. Sec. 3 discusses different neutrino production mechanisms and their impact on the initial matter power spectrum, used in our simulations pipeline. We also investigate mixed cold plus warm dark matter (C+WDM) models in Sec.3.2 and use them to put constraints on the mass and fraction of the warm component. We also use these C+WDM models as a means to cross-check our bounds on RPSN mass by issuing a mapping between the parameters involved in both cases. We recap our results in Sec. 4 and conclude in Sec. 5.

2 Probing the power spectrum with Ly-$\alpha$ Forests of QSOs

Neutral Hydrogen in the intergalactic medium (IGM) constitutes a biased tracer for the (total) matter density fluctuations at $1-100~{}h^{-1}~{}\rm{Mpc}$ scales. It is present wherever photo-ionization by the ultraviolet (UV) background is balanced by electron recombination $\gamma+\rm{H}\rightleftharpoons e^{-}+\rm{H}^{+}$. The density of neutral Hydrogen (Hi) can be obtained by equating the interaction rates of both reactions and by setting global neutrality ($n_{e}=n_{p}$):

[TABLE]

where $\sigma_{\rm{comb,ion}}$ are respectively the recombination and photo-ionization cross-sections111It should be noted that the interaction rates involved are implicit functions of temperature $T$. The dependence on $T$ is rather $n_{\rm{\textsc{Hi}}}\propto T^{-0.7}$ than $n_{\rm{\textsc{Hi}}}\propto T^{0.5}$ as is implied in Eq. 2.1. and $n_{b,\gamma_{\rm{UV}}}$ the baryon and UV photon densities.

Photons of wavelength $\lambda=121.6~{}\rm{nm}$ interacting with a neutral Hydrogen atom have a probability $\propto e^{-\tau}$ of exciting the bound electron from the fundamental to the first excited state, known as the Ly-$\alpha$ transition. The Ly-$\alpha$ optical depth in the IGM $\tau$ thus relates the observed and intrinsic flux from a background quasar with a Hydrogen-rich IGM in the foreground via

[TABLE]

in the (rest-frame) wavelengths that span between the Ly-$\alpha$ $\lambda$1216 and Ly-$\beta$ $\lambda$1026 emission lines. The optical depth is given by the Hydrogen density times the Ly-$\alpha$ absorption cross section $\sigma_{\rm{Ly}\alpha}$ integrated over distance. Since we probe density fluctuations in the Hubble flow, we measure quantities in velocity-space and thus the Ly-$\alpha$ optical depth at velocity $v$ (with respect to Earth at $v=0$) is the aforementioned quantity integrated along the line-of-sight

[TABLE]

where $\nabla v_{\parallel}$ is the velocity gradient parallel to the line-of-sight.

Because of cosmological redshift, the rest-frame wavelength of the Ly-$\alpha$ absorption gets redshifted from its observed wavelength as the quasar signal travels to Earth in an expanding Universe. Hence a series of absorption features between the Ly-$\alpha$ and Ly-$\beta$ emission lines in the spectra of high-redshift QSOs that entails the distribution of neutral Hydrogen along the line-of-sight. This Ly-$\alpha$ forest is a widely-used tool to probe density fluctuations at intergalactic scales.

2.1 From Ly-$\alpha$ forest to the flux power spectrum

Perhaps the most convenient and widely-used statistical tool to compare observations with theoretical predictions is the one-dimensional flux power spectrum [34]. It is obtained by the Fourier transform of the QSO’s transmitted flux fraction

[TABLE]

normalized by a mean transmitted flux which defines an effective optical depth $\langle\varphi\rangle=e^{-\tau_{\rm{eff}}}$.

The flux power spectrum in the Ly-$\alpha$ forest is obtained by deconvolving the power spectrum of the transmitted flux fraction defined in Eq. 2.4 from the instrumental window function $W$ after substracting the (white) noise power spectrum, and averaging over the entire set of quasars:

[TABLE]

2.1.1 SDSS / BOSS

Our sample consists of $13,821$ out of a parent sample of $\sim 60,000$ quasar spectra from the DR9 of SDSS-III/BOSS [35, 36, 37, 38, 39, 40]. They are selected for the following criteria on their Ly-$\alpha$ forest: a signal-to-noise ratio per $\Delta\lambda/\lambda=10^{-4}$ pixel greater than 2, absence of broad absorption line features, absence of damped or detectable Lyman-limit systems, and an average resolution in the Ly-$\alpha$ forest of at most $85~{}\rm{km}~{}s^{-1}$.

The spectra in this sample are used to measure the transmitted flux power spectrum in 12 redshift bins from $\langle z\rangle=4.4$ to $2.2$ (each bin spanning $\Delta z=0.2$) and in 35 equally-spaced spatial modes ranging from $k=10^{-3}$ to $2.10^{-2}~{}s~{}\rm{km}^{-1}$ (cf. left panel of Fig. 1). To reduce correlations between neighboring $z$-bins, we split the Ly-$\alpha$ forest of each quasar spectrum into up to three distinct redshift sectors. Each sector has a maximum extent of $\Delta z<0.2$. The transmitted flux power spectrum $|\tilde{\delta_{\varphi}}(k)|^{2}$ is computed separately in each $z$-sector. We checked that the resulting power spectrum agreed with that derived from a likelihood approach.

Complete details on our selection procedure as well as calibrations, computation of the flux power spectrum and determination of both statistical and systematic uncertainties are extensively described in [41].

2.1.2 VLT / X Shooter

In addition to the SDSS/BOSS data described above, we make use of the Ly-$\alpha$ forest power spectrum from the XQ-100 Legacy Survey [42], which consists of a sample of 100 medium-resolution QSO spectra observed with the VLT/XShooter instrument [43]. We measure the Ly-$\alpha$ power spectra in 70 $k$ bins for 3 redshifts centered on $\langle z\rangle=3.200,3.555,3.925$ following the methodology described in [44] and shown in the right panel of Fig. 1. We overlay the power spectra constructed with our hydrodynamical simulations in their central configuration (see Sec. 2.2) extracted from the 3 nearest redshift bins ($z=3.2,3.6,4.0$). The raw power spectrum measured in these bins asymptotically approaches the noise power spectrum computed as a white noise at small scales. Because of an uncertainty on the correction of the spectrograph resolution, we limit our study to 50, 60 and 70 $k$ bins respectively (corresponding to $k\leq 0.05,0.06,0.07~{}s~{}\mathrm{km}^{-1}$) for the 3 aforementionned redshifts. This screening ensures the raw power spectrum is always dominant over the noise power spectrum at these small scales.

2.1.3 HIRES and MIKE

Finally, we add the Ly-$\alpha$ power spectrum measured by [45] in high-resolution QSO spectra taken with the Magellan Inamori Kyocera Echelle (MIKE [46]) instrument at the Las Campanas Observatory and the High Resolution Echelle Spectrometer (HIRES [47]) at the W.M. Keck Observatory. This high-resolution power spectrum is measured in 4 redshift bins ($\langle z\rangle=4.2,4.6,5.0,5.4$) and 9 $k$ bins spanning down to $k\leq 0.08~{}s~{}\mathrm{km}^{-1}$. Since we extract our simulated power spectra from $z=4.6$ downwards, we only make use of the two lower $z$ bins for this data set, which we plot on the right-hand panel of Fig. 1 along with the simulated power spectrum extracted at the corresponding redshifts.

Our combined Ly-$\alpha$ power spectrum data thus consist of $35\times 12$ high-statistic low-resolution measurements with BOSS, $50+60+70$ low-statistic medium-resolution measurements with XQ-100, and $9\times 2$ low-statistic high-resolution measurements with HIRES/MIKE for a total of 618 measurements that span various redshifts and scales.

2.2 Probing the non-linear power spectrum at galactic scales

We compare the Ly-$\alpha$ flux power spectrum ($P_{\varphi}(k)$ herein) described in [41] and in Sec. 2.1 above with a prediction model obtained by simulating the distribution of matter in the photo-ionized IGM using the hydrodynamics code Gadget-3222http://www.mpa-garching.mpg.de/gadget/ [48, 49]. We extract the position, density and velocity fields of a baryon and a DM population in 13 equidistant $\Delta z=0.2$ redshift snapshots in $z\in\left[2.2,4.6\right]$. In each of these redshift bins, we compute the Ly-$\alpha$ optical depth from the underlying Hydrogen density distribution and IGM temperature (see Eq. 2.1) from a sample of $10^{5}$ randomly traced lines of sight. The IGM temperature is extracted from a sample of $10^{6}$ particles from which we derive the temperature-density power law slope and intercept:

[TABLE]

The power spectrum of neutral Hydrogen density fluctuations in this velocity field is derived from these computed quantities. The photo-ionization rate in the IGM is fixed in all 13 redshift bins to be in agreement with typical measurements (see for instance [50]), which is equivalent to re-normalizing the $P_{\varphi}(k)$ amplitude through the effective optical depth, which we model by a simple redshift power law:

[TABLE]

We derive the $P_{\varphi}(k)$ of $N^{3}=3072^{3}$ particles (for each species) in a $L^{3}=(100~{}\rm{Mpc})^{3}$ co-volume with periodic boundary conditions from a subset of 3 smaller simulations. We correct for their lack of resolution or size using a splicing technique described in [51]. Residuals between the exact and the spliced $P_{\varphi}(k)$ have been tested and modelled using $N=$ 1024, 1600 and 2048 simulations [52]. The simulations were run at the French TGCC333Très Grand Centre de Calcul supercluster Curie under three PRACE444Partnership for Advanced Computing in Europe and a GENCI555Grand Équipement National de Calcul Intensif allocations totalling $19\times 10^{6}$ CPU hours (see Acknowledgements).

To derive constraints on cosmological parameters, we compute the $P_{\varphi}(k)$ from a central benchmark (“best-guess” herein) model which sets our chosen free parameters at their central values in Tab. 1. Our cosmological parameters ($h$, $\Omega_{m}$, $\sigma_{8}$, $n_{s}$) are centered on the Planck 2013 [53] best fitted values. The baryon density parameter $\omega_{b}=\Omega_{b}h^{2}$ has a negligible impact on the flux power spectrum, and is therefore kept fixed to the Planck 2013 best-fit value $\omega_{b}=0.04858$. Astrophysical parameters ($T_{0}^{z=3}$, $\gamma^{z=3}$, $A^{\tau}$, $\eta^{\tau}$) describing IGM thermodynamics666the $T_{0}$ and $\gamma$ mentionned here are the logarithmic $y$-intercepts of a redshift-dependant power law for $T(z),\gamma(z)~{}\propto\left[(1+z)/4\right]^{\eta}$, taken at a pivot $z=3$. The logarithmic slopes $\eta_{\gamma}$, $\eta_{T}(z<3)$ and $\eta_{T}(z>3)$ are let free and fitted in addition to $T_{0}(z=3)$ and $\gamma(z=3)$. (Eqs. 2.6,2.7) are consistent with measurements in [54, 55, 56, 57, 50]. We then compute the first and second order derivatives of $P_{\varphi}(k)$ and use them to derive a second-order Taylor expansion around our best-guess model, with our derivatives taken at the steps in Tab. 1.

The Gadget code solves for the hydrodynamics at the onset of the non-linear regime of structure formation, which we arbitrarily set at $z=30$. We solve the Boltzmann equations in the expanding Universe numerically with the CAMB777http://camb.info [58] and CLASS888http://class-code.net/ [59] softwares and compute the linear power spectrum of matter density fluctuations at $z=30$. When implementing different NCDM cosmologies, it is useful to define a linear transfer function $T(k)$ defined such that $P_{\rm{ncdm}}(k)=T^{2}(k)P_{\rm{cdm}}(k)$, where

[TABLE]

is the linear (three-dimensional) power spectrum of total (baryon+DM) matter density fluctuations in the benchmark $\Lambda$CDM model and

[TABLE]

is the primordial scalar power spectrum. Because the Ly-$\alpha$ flux power spectrum is a unidimensional probe, we make the distinction between 1D and 3D in the above definition of the transfer function, where

[TABLE]

A complete description of our simulations and its technical aspects are extensively detailed in [60]. The implementation of our different NCDM cosmologies is detailed in Sec. 3.

2.3 Systematics and caveats

Although the Ly-$\alpha$ forest is a formidable albeit biased tracer for the matter distribution at intergalactic scales, several factors limit the accuracy of measurement and modeling of the $P_{\varphi}(k)$. These include feedback processes, such as active galactic nuclei (AGN) and supernovæ, which can potentially correlate very disparate scales. Other baryonic effects can intervene at small scales [61, 62, 63, 24, 64]. Moreover, the scarce existing measurements of the IGM thermal state and UV photo-ionization background at the relevant $z\geq 3$ redshifts compel us to have a very conservative apprehension of the astrophysical parameters. We therefore model the redshift dependence of the IGM density-temperature normalization and power index as simple power-laws and we use large step sizes in the computation of the first and second order derivatives for the Taylor expansion, in order to encompass the large observational uncertainties. We also introduce in the fitting procedure several nuisance parameters which account for uncertainties related to our estimate of the instrumental noise in the data or of the spectrograph resolution, to our modeling of the IGM, to residual biases in the splicing technique used to achieve high-resolution simulations over 100 $h^{-1}\,\rm Mpc$ scales, to supernova and AGN feedbacks, and to the redshift of reionization. Details on the scale and redshift dependence of the nuisance parameters are extensively described in [65, 52].

It should be noted that at the velocity scales probed by the medium-resolution BOSS DR9 data ($k\leq 2\times 10^{-2}~{}s~{}\rm{km}^{-1}$), these small-scale baryonic effects and the warmth of the IGM are not a limiting factor. For this reason, throughout this paper we quote bounds on NCDM masses using BOSS data only (’SDSS’), in addition to our full set of BOSS + XShooter + MIKE + HIRES data (’SDSS+XQ+HR’).

3 NCDM signatures on the matter power spectrum

When particles have relativistic enough velocities, they free-stream to a horizon scale $\lambda_{\rm{FSH}}$ effectively unaffected by gravitational potentials. Hence the matter power spectrum is suppressed below the free-streaming horizon, which is given by

[TABLE]

where the velocity dispersion $\langle v\rangle$ is given by the speed of light during the relativistic regime, and by ${\langle p\rangle}/{m}$ afterwards. For warm and cool dark matter, this transition takes place during the radiation dominated era. The associated comoving scale

[TABLE]

grows like $t^{1/2}$ during the relativistic regime, like $\ln(t)$ during the non-relativistic regime as long as radiation dominates, and remains asymptotically constant during matter domination. Finally the comoving free-streaming horizon today can be estimated from

[TABLE]

where $a_{\rm{nr}}\simeq{\langle p\rangle_{0}}/{m}$ is the scale factor at the time of the non-relativistic transition, and $\langle p\rangle_{0}$ is the momentum dispersion today. Distant QSOs probe the power spectrum at scales of several Mpc (see Sec 2.1). They enable putting upper bounds on $k_{\rm{FSH}}^{0}$ of keV NCDM particles, which translate into lower bounds on their mass. The velocity and momentum dispersion requires knowledge of the explicit distribution function of the particle, which differs from a Fermi-Dirac (thermal) distribution function depending on the production mechanism. Thus the mass bounds are different for each production mechanism.

Below we discuss several non-CDM (NCDM) models, inspired by the sterile neutrino dark matter, that are the subject of the current investigation.

3.1 Non-resonantly produced sterile neutrinos as WDM

Sterile neutrino were originally proposed as dark matter candidates by Dodelson and Widrow [4] (DW herein). In the framework of DW, sterile neutrinos are predominantly produced at $T\sim 150~{}\mathrm{MeV}(m_{\nu_{s}}/\mathrm{keV})^{1/3}$ when the oscillation production rate is most efficient while not reaching thermal equilibrium [4, 6, 66, 8]. The resulting distribution function can be roughly approximated by a rescaled Fermi-Dirac distribution [6], in which case the average momentum $\langle p\rangle$ would be identical to that of active neutrinos. The proper treatment however, based on the quantum Liouville equation [67] shows that sterile neutrinos produced in (non-resonant) oscillations do not feature a re-scaled thermal distribution and their average momentum is about 10–40% colder depending on sterile neutrino mass (see Fig. 8 in [8] or Fig. 6 in [68]). To distinguish them from the idealized DW case, we refer to sterile neutrinos produced via this mechanism as non-resonant (see Sec. 3.3 for distinction), or NRP sterile neutrinos.

The requirement $\Omega_{\nu_{s}}=\Omega_{\rm{dm}}$ fixes the $\theta$ – $m_{\nu_{s}}^{\rm{nrp}}$ relationship, represented as the upper black solid line in Fig. 8. The flux of photons from the radiative decay channel $\nu_{s}\rightarrow\gamma\nu_{\alpha}$ is a function of $\theta$ and $m_{\nu_{s}}$ [69]. Decay lines in astrophysical spectra, or the lack thereof, thus establishes constraints on these parameters [6, 70, 71, 72]. Comparing the upper bounds on the putative dark matter decay flux with currently measured DM abundance, [73] and [74] yield an upper limit of $m_{\nu_{s}}^{\rm{nrp}}\leq 4~{}\rm{keV}$ for the NRP mechanism. The non-detection of small-scale damping in the flux power spectrum of the Ly-$\alpha$ forest due to $\nu_{s}$ free-streaming has yielded lower bounds consistently above the $4~{}\rm{keV}$ limit with $\geq 5\sigma$ (see Tab. 2 below, with first constraints dating back to 2006 [19]). If right-handed neutrinos constitute all of dark matter, a growing consensus suggests they cannot be produced in this oscillation mechanism in absence of a net lepton asymmetry.

3.2 Mixed C+WDM models

A wide class of NCDM dark matter models can be approximated by adding a CDM component in addition to, e.g., thermal WDM. Such models start to deviate from CDM at scales determined by the mass of WDM component, but the overall amount of suppression is controlled by the warm DM fraction, $F_{\rm{wdm}}$, of the total dark matter. The warm-to-total DM fraction can thus be defined such that

[TABLE]

The solid and dashed transfer functions ($T(k)=\sqrt{P_{\rm{ncdm}}/P_{\rm{cdm}}}(k)$) at $z=0$ in Fig. 2 feature the free-streaming cutoff scale for different NCDM masses. As discussed above, heavier DM particles damp power on smaller scales, which makes them more consistent with the benchmark $\Lambda$CDM model, at least in the linear regime. In a cold plus warm dark matter model (C+WDM), the smaller is the fraction of the warm component, the colder is the overall transfer function as is shown by the lighter colored lines in Fig. 2 reaching an asymptotical plateau when $k\rightarrow\infty$ whose height is a function of the warm-to-cold fraction. For low $F_{\mathrm{wdm}}\lesssim 5\%$ ratios, the plateau is well approximated by $1-T(k\rightarrow\infty)\sim(1-F_{\mathrm{wdm}})$ (see [22]). For masses of a few keV, the thermal velocities of WDM particles can be neglected and thus the total dark matter distribution can be treated as a mono-species collisionless fluid in our hydrodynamical simulations whose linear transfer function is obtained by setting the DM abundance as $F_{\rm{wdm}}$ warm and $1-F_{\rm{wdm}}$ cold. The warm-to-total DM fraction $0\leq F_{\rm{wdm}}\leq 1$ is therefore an additional free parameter that interpolates between the pure WDM limit ($F_{\rm{wdm}}=1$) described above and the benchmark CDM ($F_{\rm{wdm}}=0$) limit.

We run 28 C+WDM models on a non-uniform ($\mathrm{keV}/m_{x}$, $F_{\rm{wdm}}$) grid referenced as the black dots in Fig. 6, where $m_{x}$ is the mass of the WDM component assuming it is an early-decoupled thermal relic. NRP sterile neutrinos and thermal relics feature some nearly identical transfer functions, such that one can establish a correspondance between $m_{x}$ and $m_{\nu_{s}}$ from matching their relic density and average velocity. This gives the following mass mapping:

[TABLE]

where $\kappa=4.43~{}\rm{keV}$, $\mu=4/3$, and $\omega_{\rm{wdm}}=F_{\rm{wdm}}\times\Omega_{\rm{dm}}h^{2}$ is expressed in units of $0.25\times 0.7^{2}$ [17]. In reality, this mapping is only approximate, since the actual phase-space distribution of NRP sterile neutrinos departs slightly from a rescaled Fermi-Dirac distribution [67, 8], and the authors of [76] suggest a mapping with $\kappa=3.90~{}\rm{keV}$ and $\mu=1.294$. Since our set of hydrodynamical simulations used to establish constraints on NRP neutrinos as pure WDM were run with $\mathrm{keV}/m_{x}$ as the variable, we chose to do the same for our set of C+WDM models for the sake of consistency. Sterile neutrino masses can be established simply by using the mapping in Eq. (3.5) without running additional simulations.

3.3 Resonantly-produced Sterile Neutrinos (RPSN) as Cool DM

The presence of a net lepton asymmetry at temperatures $T\sim 0.1~{}\rm{GeV}$ can significantly enhance the production of sterile neutrinos from active neutrinos through forward scattering in dense media [5]. In a mechanism similar to the MSW999Mikheyev-Smirnov-Wolfenstein effect [77, 78] accounting for the solar neutrino deficit, the excessive abundance of leptons and their conjugate neutrinos with respect to anti-leptons can yield the correct DM density $\Omega_{\rm{dm}}$ with weaker mixing angles $\theta$. The authors of [5] showed that this resonant production (RP) yields sterile neutrinos with significantly cooler momenta than the NRP ones. The resonant momenta depend on the sterile neutrino mass $m_{\nu_{s}}^{\rm{rp}}$ and net leptonic (assumed electronic) asymmetry $\mathcal{L}=(n_{\nu_{e}}-n_{\bar{\nu}_{e}})/s$ in units of entropy density where $s\propto g_{\star}T^{3}$. If the resonance occurs before the QCD phase transition, only the low momenta states are populated from the quasi thermally-distributed active neutrinos ($\langle q=p/T_{\nu}\rangle\simeq 3.15$), resulting in cooler neutrino and anti-neutrino distribution functions with $\langle q\rangle\simeq 1.6$.

The right panel of Fig. 3 displays the value of $\langle q\rangle/m_{\nu_{s}}$ for RPSN distribution functions computed in [9, 79]. The coolest distribution functions occur for given values of $\mathcal{L}$ and $m_{\nu_{s}}^{\rm{rp}}$, which we denote $\mathcal{L}^{\star}(m_{\nu_{s}})$. For larger asymmetries than $\mathcal{L}^{\star}$ for a given mass, the resonantly boosted forward scattering occurs later than the QCD phase transition, which yield quasi-Fermi populated momenta states (with weaker mixing angles). The left panel of Fig. 3 features the distribution functions of several leptonic asymmetries for a $m_{\nu_{s}}^{\rm{rp}}=8~{}\rm{keV}$ RPSN.

We run 8 hydrodynamics simulations for the following RPSN models: M3L16⋆, M4L12⋆, M6L6, M6L9⋆, M7L8⋆, M8L4, M8L8⋆ and M13L6⋆, where M is the RPSN mass in keV and L the asymmetry parameter in units of $10^{-6}$. Their linear matter power spectra at $z=30$, input to the non-linear code Gadget, are computed from the distribution functions using the CLASS software.

4 Constraints on NCDM mass from Ly-$\alpha$ forests

4.1 Mapping between C+WDM and RPSN as cool DM

As we discussed in the previous section, there are two ways one can make the dark matter distribution cooler than the Dodelson-Widrow pure WDM case. When one considers a sterile neutrino of mass $m_{\nu_{s}}$, one can thus explore the following possibilities:

• $\bullet$

the sterile neutrino constitutes the entirety of dark matter ($F_{\rm{wdm}}=1$) and is produced in absence of a net leptonic asymmetry ($\mathcal{L}=0$);

• $\bullet$

the sterile neutrino constitutes $0\leq F_{\rm{wdm}}<1$ of the total dark matter and is produced in absence of a net leptonic asymmetry ($\mathcal{L}=0$);

• $\bullet$

the sterile neutrino constitutes the entirety of dark matter ($F_{\rm{wdm}}=1$) and is produced in presence of a net leptonic asymmetry in the early Universe $\mathcal{L}>0$;

• $\bullet$

the sterile neutrino constitutes $0\leq F_{\rm{wdm}}<1$ of the total dark matter and is produced in presence of a net leptonic asymmetry in the early Universe $\mathcal{L}>0$.

In [25], we explored the first of these 4 listed scenarios and have concluded that to be consistent with Ly-$\alpha$ forest data from SDSS-III, the pure WDM sterile neutrino has to be more massive than $24.4~{}\rm{keV}$ with 95% likelihood101010the lower bound is relaxed to $16.0~{}\rm{keV}$ when adding CMB data. This lower bound is at 12$\sigma$ tension with the upper bound issued by X-ray data, set at $4~{}\rm{keV}$. Scenario 1 has thus been disfavored. Search for keV sterile neutrino DM has shifted to scenarios 2 or 3, which we investigate in this paper. We have run 36 hydrodynamical simulations with a resolution of $2\times 3072^{3}$ particles in a $(100~{}h^{-1}\rm{Mpc})^{3}$ co-volume, 28 of which explored C+WDM models ( Sec. 3.2) while the remaining 8 were dedicated to RPSN as pure cool dark matter (Sec. 3.3). To extend the sterile neutrino parameter space we sample, we use the 28 C+WDM configurations in addition to the 8 hydrodynamical simulations modelling RPSN. We take advantage of the fact that both yield similar transfer functions up to some $k$ scale to map out a correspondance between ($\mathcal{L}>0$, $F_{\rm{wdm}}=1$) and ($\mathcal{L}=0$, $0\leq F_{\rm{wdm}}<1$) assuming particles of the same mass $m_{\nu_{s}}$.

Fig. 4 illustrates the correspondance between the coolest RPSN model of 4 keV ($\mathcal{L}^{\star}=1.2\times 10^{-5}$) and the 4 keV neutrino produced in absence of a lepton asymmetry that constitutes $\sim 35\%$ of the total dark matter. The correspondance is obtained with a least-square method out to $k_{\rm max}=1.35\,h^{-1}{\rm Mpc}$ on the linear transfer function. Because Ly-$\alpha$ forests are a unidimensional probe for the matter distribution, we perform our mapping on the 1D transfer function, obtained using Eq. 2.10. This $\mathcal{L}-F_{\mathrm{wdm}}$ mapping enables us to convert our bounds on ($m_{\nu_{s}}$, $F_{\rm{wdm}}$) obtained with our 28 C+WDM simulations to bounds on ($m_{\nu_{s}}$, $\mathcal{L}$).

For illustration purposes, the top left panel in Fig. 5 shows the linear power spectra of the M8L8⋆ (in thick black) along with the closest matching linear matter $T_{\mathrm{1d}}(k)$ of a C+WDM model assuming $m_{\nu_{s}}=8~{}\rm{keV}$, which occurs for a warm-to-total DM fraction of $F_{\rm{wdm}}=25\%$ (in thin black). As illustrated in Fig. 4, this linear transfer function correspondance is adequate up to some $k$ scale, beyond which the corresponding C+WDM model $T_{\mathrm{1d}}(k)$ breaks away from its comparative RPSN model to an asymptotical plateau ($T_{\mathrm{1d}}(k\rightarrow\infty)\propto(1-F_{\mathrm{wdm}})\geq 0$). For most values of $m_{\nu_{s}}$ explored in this work, this breakaway $k$ is beyond the scales probed by our Ly-$\alpha$ forest data set, which we’ve materialized on Fig. 5. We illustrate the negligible impact of differences in the linear 1D transfer function beyond the breakaway scale by overplotting three C+WDM models that exhibit similar $T_{\mathrm{1d}}$ on large scales. The differences in the non-linear regime measured by the flux power spectra are within the statistical uncertainties of the simulations, almost an order of magnitude smaller than data uncertainties on similar scales.

While we derive constraints in the sterile neutrino parameter space through the mapping from our C+WDM grid described above, we also use the eight exact RPSN hydrodynamical simulations to determine the accuracy of our mapping procedure. To this end, we compare the $\chi^{2}$ obtained using the exact (non-linear) RPSN simulation with that of its corresponding C+WDM model obtained with our mapping procedure in the linear regime (see Fig. 5 for illustration with M8L8*⋆*). This cross-check is done for both sets of data, SDSS/BOSS alone (’SDSS’) and combined with VLT/XShooter, Keck/HIRES and LCO/MIKE (’SDSS+XQ+HR’). We observe a systematic shift at the level of $0.2\,\sigma$ on average in the first case, and of $0.5\,\sigma$ in the second, with the RPSN simulation showing a smaller $\chi^{2}$ (better agreement with the data) than its C+WDM match. The better agreement for SDSS only is consistent with a better match of the transfer function on large scales. The results we present hereafter for the RPSN models are corrected for this systematic shift.

4.2 Constraints on C+WDM mass and fraction

The probability distribution in the $(F_{\rm{wdm}},m_{x})$ plane being strongly non-Gaussian, a Taylor expansion in either of these two parameters would not provide accurate results. We therefore extend the method described in [60, 65, 52, 25] in the following way. We use the likelihood of previous work, based on a second-order Taylor expansion, to capture the dependence of the Ly-$\alpha$ flux power spectrum with the cosmological and astrophysical variables of Table 1, and to model identified nuisance parameters accounting for IGM thermal state modeling, re-ionization redshift, spectrometer noise, and simulation uncertainties. To capture the dependence with $F_{\rm{wdm}}$ and $m_{x}$, we produced a grid of 28 C+WDM simulations with non-zero $(F_{\rm{wdm}},m_{x})$ while all other parameters are set to their best-guess value in Table 1. For each of the C+WDM simulations, we compute a $\chi^{2}$ with respect to our $35\times 12$ $P_{\varphi}(k,z)$ data points from the BOSS DR9 described in Sec. 2.1, assuming a Gaussian distribution for the Hubble parameter, $h=0.673\pm 0.010$ (Planck 2015), and minimizing over all the parameters of the likelihood described above. We interpolate within this grid to predict the probability distribution at any point in the $(F_{\rm{wdm}},m_{x})$ plane.

We identify the minimal value of $\chi^{2}$ using the MINUIT package [80], letting all parameters free. We set a confidence level (CL) on any parameter $\theta_{i}$ (out of $n$) by minimizing the $\chi^{2}$ function on all remaining $n-1$ parameters for each scanned value of $\theta_{i}$. To set confidence levels on a hypersurface of 2 parameters ($\theta_{i}$, $\theta_{j}$), the $\chi^{2}$ minimization is performed on the $n-2$ remaining parameters. Assuming all experimental errors are normally distributed,

[TABLE]

[TABLE]

where $\Gamma(z)=\displaystyle\int_{0}^{\infty}\mathrm{d}x~{}x^{z-1}e^{-x}$ is the gamma function. $1\sigma$, $2\sigma$ and $3\sigma$ confidence levels we quote on parameter $\theta_{i}$ or ($\theta_{i}$, $\theta_{j}$) correspond to a $\chi^{2}$ value with respect to the minimum value of $\Delta\chi^{2}(\theta_{i})=1,4$ and 9 and $\Delta\chi^{2}(\theta_{i},\theta_{j})=2.30,6.18$ and 11.83 respectively.

Fig. 6 displays a cubic-interpolated map of the $\Delta\chi^{2}$ probability in the ($\mathrm{1\,keV}/m_{x}$, $F_{\rm{wdm}}$) plane using the 28 C+WDM flux power spectra computed using our hydrodynamical simulations. The $\chi^{2}$ values we obtain along the $F=1$ axis are in excellent agreement with the ones we reported in [25] and [44] for 95% confidence level limits in the pure WDM case, which led to $m_{x}\geq 4.09~{}\rm{keV}$ and $m_{\nu_{s}}^{\rm{nrp}}\geq 24.4~{}\rm{keV}$ using our SDSS Ly-$\alpha$ sample only and to $m_{x}\geq 4.65~{}\rm{keV}$ and $m_{\nu_{s}}^{\rm{nrp}}\geq 28.8~{}\rm{keV}$ using in addition the XQ-100, HIRES and MIKE data. The bounds shown in Fig. 6 at $2\sigma$ are of course weaker than the ones reported previously since they are here computed for a full two-dimensional analysis in the ($\mathrm{1\,keV}/m_{x}$, $F_{\rm{wdm}}$) plane, thus for $\Delta\chi^{2}=6.18$ instead of 4.

Independently of the redshift range considered for SDSS and of the use or not of the higher-resolution data, the contours in the ($\mathrm{1\,keV}/m_{x}$, $F_{\rm{wdm}}$) plane can be approximated by $F_{\rm{wdm}}=\alpha(1\,{\rm keV}/m_{x})^{\beta}$, with $\beta=-1.37$ in all cases. For the 95% C.L. limits, the normalization is $\alpha=0.24$ for the full (12-bin) SDSS data, and $\alpha=0.14$ for SDSS+XQ+HR. For the $3~{}\sigma$ limits, the normalization is $\alpha=0.35$ for SDSS and $\alpha=0.20$ for SDSS+XQ+HR. Therefore, thermally decoupled relics as light as $m_{x}\geq 0.7~{}\rm{keV}$ are consistent with Ly-$\alpha$+$H_{0}$ data (95% CL) if they constitute 15% of the total dark matter or less. The contribution of warm to total dark matter is reduced to $\sim 10\%$ when including higher-resolution data. The bounds are derived under the assumption that the IGM thermal history can be modeled with 5 parameters: a broken power law for $T_{0}$ and a simple power law for $\gamma$. Non-monotonic thermal histories as discussed in [24], for instance, are not included in the marginalization. Although the impact of such a hypothesis is small for SDSS-only bounds, which can therefore be considered as conservative lower bounds, it can lead to looser limits in the case of SDSS+XQ+HR.

In the pure WDM case, we showed that the highest two redshift bins of the SDSS data significantly tightened the 95% CL limit despite their low statistical significance [25]. The same is true here for the study in the full ($\mathrm{1\,keV}/m_{x}$, $F_{\rm{wdm}}$) plane. Considering the lowest ten redshift bins only (i.e., redshifts in $2.1<z<4.1$) loosens the bound by about 25% on $\mathrm{1\,keV}/m_{x}$ for a given $F_{\rm{wdm}}$. The $\chi^{2}$ of the best fit increases by 71.7 when including the 70 SDSS data points at $z>4.1$, indicating that the highest two redshift bins are consistent with the bins at lower redshift.

The best-fit parameters for either a pure WDM or a C+WDM scenario, shown in Table 3, are all in excellent agreement with our previous WDM analyses [25, 44]. The IGM temperature history is compatible with other recent estimates [75] and the cosmological parameters are compatible with the latest Planck results, except for a $\sim 2\sigma$ tension on $n_{s}$ when fitting BOSS data alone. As was shown in [25], the lower preferred value of $n_{s}$ in BOSS Ly$\alpha$ data compared to CMB has an impact on the constraint one can set on the mass of a pure WDM particles. The use of BOSS Ly$\alpha$ data alone, or, equivalently, allowing for a running of $n_{s}$ that accommodates for the different values of $n_{s}$ on large (CMB regime) and small (Ly$\alpha$ regime) scales, leads to tighter constrains than when fitting BOSS and Planck data together in the absence of running. A similar effect is true here. Approximating the 95% CL contour by $F=\alpha(1\,{\rm keV}/m_{x})^{\beta}$ as we did above, we obtain a constraint on $m_{x}$ that is looser by about 35% for fixed $F_{\rm{wdm}}$ when adding Planck to BOSS data. The situation is different with the extended SDSS+XQ+HR data, for which it was shown in [44] that the tension on $n_{s}$ was mostly resolved. We thus expect similar constraints on $m_{x}$ whether or not we include CMB data in addition to this extended set. For fixed $F_{\rm{wdm}}$, the 95% CL constraint on $m_{x}$ indeed shifts by less than 12% between the two configurations (extended Ly$\alpha$ alone or with the addition of Planck).

4.3 Constraints on RPSN mass

We convert the bounds on ($m_{x}$, $F_{\mathrm{wdm}}$) established in the previous subsection into bounds on ($m_{\nu_{s}}$, $\mathcal{L}$) using the mapping procedure described in Sec. 4.1, corrected for systematic shifts using the $P_{\varphi}$ from our 8 RPSN hydrodynamical simulations. Fig. 7 displays the 1$\sigma$, 2$\sigma$ and 3$\sigma$ confidence level contours in the ($m_{\nu_{s}}$, $\mathcal{L}$) plane using our SDSS/BOSS Ly-$\alpha$ data (left panel) or the full SDSS+XQ+HR data (right panel). The tension in standard deviations with these two data sets is reported in Table 4 for the relevant RPSN models for which we ran hydrodynamical simulations. A lepton asymmetry during the era of RPSN production in the early Universe boosts the oscillation frequency from active to sterile neutrinos, thus enabling ample production of dark matter sterile neutrinos with weaker mixing angles $\theta$. We run a Boltzmann code that computes the mixing angle for a given dark matter density, mass and lepton asymmetry parameter. Fig. 8 displays the quantity $\sin^{2}2\theta$ as a function of $m_{\nu_{s}}$ assuming $\Omega_{\mathrm{dm}}h^{2}=0.26142\times 0.675^{2}\simeq 0.119$ and values of the lepton asymmetry parameter in the range $0\leq\mathcal{L}\leq 7\times 10^{-4}$. The black lines in Fig. 8 display the relation between mass and mixing angle for eight values of the primordial lepton asymmetry shown along each black curve (in units of $10^{-6}$). The $\mathcal{L}_{6}=0$ thick line corresponds to non-resonant production. Values above $\mathcal{L}\gtrsim 10^{-3}$ are inconsistent with Big Bang nucleosynthesis (BBN). For $\sin^{2}2\theta\gtrsim 10^{-7}$, dark matter is overproduced. The region excluded at the $3\sigma$ CL by SDSS/BOSS (resp. by SDSS + XQ + HR) data is shaded in blue (resp. in red). The dashed blue line indicates the limit using the 10 lowest redshift bins only, i.e., excluding $z=4.2$ and 4.4.

As expected from the white stripe visible on the right panel of Fig. 3, the coolest RPSN models, which occur for $\mathcal{L}=\mathcal{L}^{\star}(m_{\nu_{s}})$, feature the longest free-streaming length and are more consistent with Ly-$\alpha$ forest data than other values of $\mathcal{L}$. This is visible on Fig. 7 as a horn-like valley in the $\chi^{2}$ maps, which extends to sterile neutrino masses around $\sim 7~{}\mathrm{keV}$ in the right-hand panel. This area is of particular interest since it matches the range of masses and mixing angles for which the $3.55~{}\mathrm{keV}$ X-ray signal reported in [82, 81, 83] can be interpreted as photons emitted by the decay of a $7.1~{}\mathrm{keV}$ right-handed neutrino. Although this region exhibits a $\sim 3\sigma$ tension with the SDSS+XQ+HR Ly-$\alpha$ data, two caveats should be considered. First, IGM thermal histories impact the small scales ($0.02~{}\leq~{}k/s~{}\mathrm{km}^{-1}~{}\leq~{}0.07$) probed by these high-resolution data. Although we marginalize over 5 parameters to describe the thermal history (as explained above), more general models (non-monotonic temperature evolution for instance) could loosen our constraint. Second, the flux power spectrum exhibits large gradients with respect to the RPSN parameters around the “horn” region, where the interpolation procedure is thus more delicate. Therefore, because of the interest of this region for RPSN constraints, we located our eight RPSN simulations in that area: six correspond to the coolest models for their mass (M3L16⋆, M4L12⋆, M6L9⋆, M7L8⋆, M8L8⋆ and M13L6⋆), the remaining two (M6L6 and M8L4) being slightly warmer than their corresponding counterparts at the same mass (M6L9⋆ and M8L8⋆). The results presented in Table 4 show that the horn is a real feature, although the exact location of its boundaries might require additional hydrodynamical simulations to assess. Hence the shape of the blue and red contours on Fig. 8 may be less accurate in the regions around the 6 bottom-most gold squares that correspond to our coolest RPSN models. The neutrino decay origin of the $3.55~{}\mathrm{keV}$ X-ray line, shown as the black dot with error bars, is located in this region.

For the reasons just stated, we suggest scanning the area around $m_{\nu_{s}}=7.1~{}\mathrm{keV}$ and $\sin^{2}2\theta=4.9\times 10^{-11}$ with a set of dedicated hydrodynamical simulations in order to properly account for the strong dependence of the power spectrum on model parameters in that region. These simulations should also implement the different IGM thermal histories prognosticated in [24]. We leave this for future work.

5 Conclusion

We used the SDSS-III/BOSS DR9 Ly-$\alpha$ forest data to constrain warm dark matter models. In the previous study [25], these data were used to put bounds on the mass of warm dark matter in the form of thermal relics. In this paper, we extend the previous results to two important classes of non-thermal WDM models: cold-plus-warm dark matter and sterile neutrinos produced via mixing with active neutrinos in the presence of a net lepton asymmetry (known as resonantly produced sterile neutrinos). While many works have used Ly-$\alpha$ forest to constrain thermal relic WDM (and non-resonantly produced sterile neutrinos) [16, 17, 19, 18, 21, 44, 75, 25], in the current study, for the first time, RPSN are constrained by running a set of hydrodynamical simulations with realistic initial power spectra predicted by this scenario.

The size of the SDSS/BOSS DR9 quasar dataset significantly reduces statistical uncertainties compared to previous SDSS-II data. The systematic uncertainties associated with feedback and IGM modeling, instrument noise etc, now become comparable to the statistical ones. Uncertainties in the IGM thermal history, in particular, mostly affect the highest redshift bins. To provide both conservative (although weaker) and more ambitious (although more prone to systematics) bounds, we computed all limits with and without the two highest z = 4.2, 4.4 redshift bins. These bounds can be additionally strengthened by including the higher-resolution XQ-100, MIKE and HIRES (XQ + HR) Lyman-$\alpha$ datasets that probe the power spectrum at smaller scales and higher redshifts (see Fig. 1). In this regime, however, the flux power spectrum exhibits a suppression on the smallest scales that a number of astrophysical effects can lead to: Doppler broadening, pressure smoothing, expansion of the filaments along the line of sight etc., [21, 24]. Dark matter free-streaming could also contribute to this suppression.111111A possible way to identify the mechanism responsible for the small-scale suppression of the flux power spectrum is to measure the thermal history of the IGM independently of Ly-$\alpha$ forest data. To this end, a novel method — Gaussian optical depth decomposition — was proposed in [84]. It is based on the idea that for high-resolution spectra it is possible to identify individual absorption lines and to infer gas temperature directly via measurement of their broadening. This project is currently underway. The detailed analysis of the influence of the thermal history on the flux power spectrum at small scales is beyond the scope of this study. In this work, we focus on thermal histories that can be modeled with 5 parameters: a broken power law for $T_{0}$ and a simple power law for $\gamma$.

RPSNs lighter than $m_{\nu_{s}}\sim 3.5$ keV produce a Ly-$\alpha$ flux power spectrum inconsistent with the one measured from SDSS/BOSS at more than $3\sigma$. This bound tightens to $m_{\nu_{s}}\sim 7$ keV when including the higher resolution data. Heavier sterile neutrinos can be consistent with data if one assumes lepton asymmetries of about $\mathcal{L}\sim 10^{-6}-10^{-5}$ at the time of RPSN production. Outside this range of lepton asymmetries, both limits tighten by about a factor 2. Models that would potentially explain the observed $3.5~{}\mathrm{keV}$ line are compatible with SDSS/BOSS 10-bin and 12-bin data at a level between $2\sigma$ and $3\sigma$, and with the extended SDSS+XQ+HR data at slightly above $\sim 3\sigma$.121212Curiously, the 7 keV RPSN with the mixing angle within the range reported in [81, 82] produces the suppression of the flux power spectrum at small scales that would fit the XQ + HR data if the intergalactic medium is cold at redshifts $z\geq 5$ [24]. We notice that the distribution function of RPSN for the coolest models (those in the region of the white stripe of Fig. 3) varies drastically with $m_{\nu_{s}}$ and $\mathcal{L}$. It is therefore challenging to obtain precise constraints in this region of the RPSN parameter space. We therefore also run 8 hydrodynamical simulations with exact input for RPSN models, among which one corresponding to a 7 keV neutrino produced in a lepton asymmetry $\mathcal{L}=8\times 10^{-6}$. This corresponds to the coolest, and thus most conservative, model for $m_{\nu_{s}}=7~{}\mathrm{keV}$. Its power spectrum is consistent at $1.9~{}\sigma$ with BOSS data, and at $3.1~{}\sigma$ with BOSS + XQ + HR. We leave further dedicated investigation of the region of interest for future work.

The combined BOSS and XQ + HR Ly-$\alpha$ data sets are only compatible with a very restricted region in parameter space of RPSNs. This result is consistent with the conclusion of [85] based on the clustering of SDSS galaxies.

Finally, in the case of C+WDM, we find that thermally decoupled relics as light as $m_{x}\geq 0.7$ keV are consistent at 95% C.L. with SDSS Ly-$\alpha$ data if they constitute less than 15% of the total dark matter. The bounds on $F_{\rm{wdm}}$ tighten to 10% when derived from the SDSS + XQ + HR data set. More generally, we find that the limits can be well approximated by $F_{\rm{wdm}}=\alpha(1\,{\rm keV}/m_{x})^{\beta}$, with $\beta\sim-1.37$. At 95% C.L., $\alpha=0.24$ when using SDSS only and $\alpha=0.14$ with SDSS+XQ+HR. At $3~{}\sigma$, $\alpha=0.35$ for SDSS only and $0.20$ for SDSS+XQ+HR.

RPSNs require significant lepton asymmetry to be present in the primordial plasma during the epoch of dark matter production [5, 9]. Such a lepton asymmetry can be generated e.g. in the $\nu$MSM [86], minimal extension of the Standard Model with 3 right-handed neutrinos that would explain both the dark matter but also the baryon asymmetry of the Universe as well as neutrino masses and oscillations [87, 88, 89].

Sterile neutrino dark matter remains to this day an active area of research. We focused on sterile neutrinos produced by the mixing with active neutrinos, yet other production mechanisms can yield cool transfer functions as well (see [26] for a review). Our bounds do not apply to these alternative production scenarios. Our work shows that Ly-$\alpha$ forest constraints are consistent with those of other studies. Considering the possibility that sterile neutrinos are a fraction of dark matter (just like left-handed neutrinos are), any sterile neutrino mass is still allowed, given that it constitutes at most $\sim 10\%$ of the total dark matter density and is consistent with upper bounds from X-ray data.

Acknowledgments

We thank Matteo Viel for providing us high resolution quasar spectra from the MIKE and HIRES, as well as Volker Springel for making GADGET-3 available to our team. JB thanks Alexander Merle and Kevork Abazajian for their helpful contributions and discussions on the resonant and non-resonant production mechanisms. This project received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (GA 694896).

We acknowledge PRACE (Partnership for Advanced Computing in Europe) for access to thin and xlarge nodes on the Curie cluster based in France at the TGCC (Très Grand Centre de Calcul) under allocation numbers 2010PA2777, 2014102371 and 2012071264. We also acknowledge the French national access to high-performance computing GENCI (Grand Équipement National de Calcul Intensif) for access to the Curie cluster under allocation t2016047706.