Constraining the mass of light bosonic dark matter using SDSS Lyman-$\alpha$ forest

TL;DR

This study constrains the mass of ultra-light fuzzy dark matter using Lyman-alpha forest data from SDSS, employing cosmological simulations to compare theoretical predictions with observations, and excludes certain mass ranges at 95% confidence.

Contribution

It provides the first robust constraints on fuzzy dark matter mass using Lyman-alpha forest data, incorporating hydrodynamical simulations and a novel scaling law with warm dark matter models.

Findings

Excluded fuzzy dark matter masses between 10^{-22} and 2.3×10^{-21} eV at 95% CL.

Extended exclusion up to 2.9×10^{-21} eV with higher-resolution data.

Quantum effects are negligible for masses above 10^{-22} eV in the relevant scales.

Abstract

If a significant fraction of the dark matter in the Universe is made of an ultra-light scalar field, named fuzzy dark matter (FDM) with a mass $m_a$ of the order of $10^{-22}-10^{-21}$ eV, then its de Broglie wavelength is large enough to impact the physics of large scale structure formation. In particular, the associated cut-off in the linear matter power spectrum modifies the structure of the intergalactic medium (IGM) at the scales probed by the Lyman-$\alpha$ forest of distant quasars. We study this effect by making use of dedicated cosmological simulations which take into account the hydrodynamics of the IGM. We explore heuristically the amplitude of quantum pressure for the FDM masses considered here and conclude that quantum effects should not modify significantly the non-linear evolution of matter density at the scales relevant to the measured Lyman-$\alpha$ flux power, and for $m_a \geq 10^{-22}$ eV. We derive a scaling law between $m_a$ and the mass of the well-studied thermal warm dark matter (WDM) model that is best adapted to the Lyman-$\alpha$ forest data, and differs significantly from the one infered by a simple linear extrapolation. By comparing FDM simulations with the Lyman-$\alpha$ flux power spectra determined from the BOSS survey, and marginalizing over relevant nuisance parameters, we exclude FDM masses in the range $10^{-22} \leq m_a < 2.3\times 10^{-21}$ eV at 95 % CL. Adding higher-resolution Lyman-$\alpha$ spectra extends the exclusion range up to $2.9\times 10^{-21}$ eV. This provides a significant constraint on FDM models tailored to solve the "small-scale problems" of $\Lambda$CDM.

Figures

Equations

$$k_c=6.46{\left(\frac{m_X}{\mathrm{keV}}\right)}^{1.11}{\mathrm{Mpc}}^{-1}\mathrm{where}\sqrt{\frac{P_{\mathrm{WFM}}(k_c)}{P_{\mathrm{CDM}}(k_c)}}=0.5$$

$$k_J(\rho )$$

$$=$$

$$\frac{m_X}{\mathrm{keV}}\sim \frac{\sqrt{m_a{M}_{\mathrm{Pl}}}}{\mathrm{keV}}\sim 0.5{\left(\frac{m_a}{10^{-22}\mathrm{eV}}\right)}^{0.5}$$

$${\ddot{\Psi }}_0+3H(z){\dot{\Psi }}_0+m_a^2{\Psi }_0=0$$

$$c_{\mathrm{s},\mathrm{eff}}^2=\frac{k^2/4m_a^2}{1+k^2/4m_a^2}$$

$$k_c=4.97{\left(\frac{m_a}{10^{-22}\mathrm{eV}}\right)}^{0.46}{\mathrm{Mpc}}^{-1}$$

$$m_X=0.79{\left(\frac{m_a}{10^{-22}\mathrm{eV}}\right)}^{0.42}\mathrm{keV},$$

$$P_{1D}(k)=\frac{1}{2\pi }{\int }_k^{\infty }d\tilde{k}\tilde{k}{P}_{3D}(\tilde{k}),$$

$${\partial }_t\vec{v}+H\vec{v}+\frac{1}{a}(\vec{v}\cdot \nabla )\vec{v}=-\frac{1}{a}\nabla \left[\varphi -\frac{{\hslash }^2}{2m_a^2{a}^2}\left(\frac{{\nabla }^2\sqrt{\rho }}{\sqrt{\rho }}\right)\right]$$

$$T_{\mathrm{IGM}}=T_0(z)(1+\delta \rho /\rho {)}^{\gamma (z)-1}$$

$${\tau }_{\mathrm{eff}}=\alpha \times (1+z{)}^{\beta }\mathrm{with}\alpha =0.0025\mathrm{and}\beta =3.7$$

$$r_Q\equiv \frac{|\nabla Q/m_a|}{|\nabla \varphi |}$$

$$m_X=0.715\times {\left(\frac{m_a}{10^{-22}\mathrm{eV}}\right)}^{0.558}\mathrm{keV}$$

Full text

††pubyear: 2017††pagerange: Constraining the mass of light bosonic dark matter using SDSS Lyman-$\alpha$ forest.–Constraining the mass of light bosonic dark matter using SDSS Lyman-$\alpha$ forest.

1 Introduction

While non-baryonic dark matter is a cornerstone in our current understanding of the observable Universe, its nature remains completely unknown after decades of theoretical and experimental investigations. Cold dark matter (CDM) is necessary to understand at the same time the global dynamics of the Universe and the gravitational formation of large scale structures from an initial, almost scale-invariant, power spectrum. The current-time dark matter density derived from CMB observations, $\Omega_{c}h^{2}=0.118\pm 0.001$ (Ade et al., 2016) is six times larger than the baryon density $\Omega_{b}h^{2}=0.0222\pm 0.0002$, which can also be independently derived from the measurements of light element abundances. Dark matter also seems to be necessary to understand the dynamics of gravitationally bound systems such as galaxies (Persic et al., 1996) and clusters of galaxies (Clowe et al., 2006).

To understand the nature of dark matter, two complementary strategies are possible. A first option is to test specific models inspired by particle physics, using more or less direct detection strategies. Two outstanding candidates predicted in extensions of the Standard Model of particle physics are the supersymmetric neutralino, a thermal relic whose mass is in the GeV$-$TeV range, and the QCD axion, a relic pseudo-scalar field with a mass around $1-100\,\mu\rm{eV}$. They are tested with dedicated experiments or observations, without a positive outcome as of now (Olive et al., 2016) although future experiments will explore much more of their parameter space, in particular in the case of the currently poorly tested QCD axion. A second option consists in measuring or constraining generic properties of dark matter from astrophysical or cosmological observations. For example, assuming that the dark matter is made of fermions, then estimations of their phase-space density in dwarf galaxies imply the Tremaine-Gunn bound (Tremaine & Gunn, 1979) on their mass, $m\gtrsim 400$ eV. Such a bound does not apply for bosonic particles.

The small-scale clustering features of visible matter may also be used in different ways to constrain dark matter properties. In fact, at the galactic or sub-galactic scale it has been argued that CDM predictions do not match observations, concerning the number density of low-mass and dwarf galaxies (Klypin et al., 1999; Moore et al., 1999; Boylan-Kolchin et al., 2012), and the absence of an observed central dark matter cusp in these objects (Oh et al., 2011). While models which take into account the full richness of baryon physics at these scales may surpass those issues, it is still possible that these features result from specific DM properties, not included in the most naive CDM model.

A popular scenario to explain the apparent lack of structures at small scales is that the DM velocity dispersion could be large enough at the time of equality, when structures start to form, so that the related free-streaming of these warm DM (WDM) particles would suppress small-scale structures. For thermal WDM of mass $m_{X}$, using cosmological parameters from Ade et al. (2016), the associated cutoff in the linear density power spectrum $P(k)$ is at the scale (Viel et al., 2005) :

[TABLE]

In this equation, $k_{c}$ does not depend on $h$. Recent observations of a 3.5 keV X-ray line in several DM-rich objects (Bulbul et al., 2014; Cappelluti et al., 2017) were also interpreted as a hint for a 7 keV sterile neutrino which would constitue WDM. An issue with WDM is that it seems difficult to solve the halo abundance and the cusp-core problems at the same time, if both are to be taken seriously. This is the "Catch 22" problem (Macciò et al., 2012) : for WDM with a mass of a few keV, the halo core size, which originates from the above-mentioned Tremaine-Gunn phase-space density argument, is too small to solve the apparent cusp-core issue.

Other DM properties were suggested and studied to alleviate the small-scale CDM issues. One possibility is that DM is self-interacting with large cross-sections of the order of $0.1-1$ cm2/g, thus producing constant-density cores inside halos (Spergel & Steinhardt, 2000; Rocha et al., 2013). The free-streaming effect characteristic of WDM could also be produced by decaying DM with a small mass splitting (Wang et al., 2013).

In this article, we now focus on the alternative scenario according to which DM consists in extremely light bosons, with a mass scale $m_{a}\sim 10^{-22}$ eV, so that the quantum wave properties of these particles are large enough to smooth the associated density fluctuations on the relevant small scales. The basic physics involved at large scales in this Fuzzy Dark Matter (FDM) scenario is presented in eg. Hu et al. (2000); Marsh (2016b); Hui et al. (2017). The effective "quantum pressure" associated with FDM generates an effective Jeans wavenumber, which is the following for a fluid of density $\rho$ :

[TABLE]

Structure formation is suppressed for comoving scales smaller than the corresponding Jeans scale at equality time, which is roughly $\sim 10$ Mpc*-1* for $m_{a}=10^{-22}$ eV. This results in a linear power spectrum with a small-scale cutoff produced at high redshift. The related phenomenology is therefore close to the WDM one. On the other hand, in dense environments such as the center of DM halos, where the effective Jeans scale is small due to the $\rho^{1/4}$ scaling, simulations which take explicitly into account these wave properties (Woo & Chiueh, 2009; Schive et al., 2014; Veltmaat & Niemeyer, 2016; Zhang et al., 2016) demonstrate a rich phenomenology including the development of solitonic cores inside Navarro-Frenk-White halos. In particular, these solitonic cores can fit the observed kinematic properties of some dwarf galaxies, with FDM masses ranging from $\sim 0.4\times 10^{-22}$ eV up to $\sim 4-6\times 10^{-22}$ eV depending on the considered objects and analysis procedure (Schive et al., 2014; Marsh & Pop, 2015; Calabrese & Spergel, 2016; Chen et al., 2017; Gonzáles-Morales et al., 2016).

There is to our knowledge no strong argument to explain why the DM mass would be such that their de Broglie wavelength is of the order of galactic scales. On the other hand, light bosonic pseudo-scalars are a generic prediction in string models, for example as the zero-mode of the Kaluza-Klein tower resulting from the compactification of antisymmetric tensor fields, with the mass scale arising due to the modulus fixing the GUT scale coupling (Svrcek & Witten, 2006). It is plausible that DM is made of several such species, one of which would have the right mass to produce the FDM phenomenology. However for simplicity, we consider in this article only the case of a single (fuzzy) DM species.

Existing bounds on the mass of FDM particles arise from both linear and non-linear cosmological probes. The most robust contraints are based on linear probes of structure formation, mostly the CMB, which show that FDM is not a dominant component of dark matter for $m_{a}\lesssim 10^{-25}$ eV (Hlozek et al., 2015). Models with higher values for $m_{a}$ cannot be distinguished from CDM given that the scales for which structure formation is damped with respect to CDM are too small to be observed with the primary CMB anisotropies. Pushing into the mildly non-linear regime, CMB lensing at high $\ell$ could explore FDM masses up to $\sim 10^{-23}$ eV (Hložek et al., 2017). For even higher-mass FDM, the halo mass function (HMF) was semi-analytically computed from the linear power spectrum, using formalisms similar to the case of pure CDM (Marsh & Silk, 2014; Marsh, 2016a; Du et al., 2017) : for the benchmark mass $m_{a}=10^{-22}$ eV, the HMF starts to be suppressed for $M\la 10^{10}\,M_{\sun}$, while it is fully cut off at $M\sim 10^{7}-10^{8}\,M_{\sun}$ depending on the calculation. The galaxy UV-luminosity function at high redshift $z\sim 4-10$ is therefore expected to be strongly attenuated for $M_{\rm UV}\ga-16$, and accordingly reionization is also expected to be delayed (Bozek et al., 2015; Sarkar et al., 2016). Recent estimations of the high-redshift UV-luminosity function (Bouwens et al., 2015; Livermore et al., 2017), as derived from deep observations such as the Hubble Frontier Fields, are still subject to large uncertainties (Bouwens et al., 2017) but they do not clearly point towards a deviation from the CDM Schechter law. FDM bounds could therefore be derived, ranging from $m_{a}>1.2\times 10^{-22}$ eV to $5\times 10^{-22}$ eV or even $1.2\times 10^{-21}$ eV, depending on the data used and analysis strategy (Schive et al., 2016; Corasaniti et al., 2017; Menci et al., 2017).

The absorption of the Lyman-$\alpha$ emission from distant quasars by neutral H i in the intergalactic medium (IGM) constitutes a powerful, high-redshift probe for spatial fluctuations of the matter density at comoving scales going down to $\sim 1-0.1$ Mpc depending on the instrumental wavelength resolutions. As in the case of WDM, it is therefore expected that Lyman-$\alpha$ forest observations can constrain FDM models, with very different sources of uncertainties with respect to the aforementioned measurements. In fact, Amendola & Barbieri (2006) disfavors pure FDM models with $m_{a}<10^{-22}$ eV, by making use essentially of high-resolution Lyman-$\alpha$ forest spectra. Given the similar phenomenologies of FDM and WDM, it is also possible to apply existing WDM bounds on $m_{X}$ (eg. Viel et al. (2013a); Baur et al. (2016)) to the FDM scenario, by using the correspondance between $m_{a}$ and $m_{X}$ which matches the cutoff position $k_{c}$ in the linear matter power spectrum (Marsh & Silk, 2014). This correspondance can roughly be estimated as highlighted in Marsh (2016b) : the power spectrum suppression for WDM takes place for modes inside the horizon when $T\sim m_{X}$, while for FDM the suppressed modes are inside the horizon when the scalar field starts oscillating, i.e. for $H(T)\sim m_{a}$. Since $T\sim\sqrt{HM_{\rm Pl}}$ (radiation era), the order-of-magnitude mass matching between both models is thus

[TABLE]

This relation is confirmed by more sophisticated calculations. The bound $m_{X}>4.09$ keV (Baur et al., 2016) would therefore naively translate into $m_{a}\ga 7\times 10^{-21}$ eV. However the significant differences of power spectrum shapes between FDM and WDM at the cutoff scale, as computed in Marsh & Silk (2014); Schive et al. (2016), make this interpretation uncertain and, as we will see, too optimistic.

Given the potential of the Lyman-$\alpha$ forest to constrain FDM, it is desirable to dedicate more in-depth studies on this scenario. In this article we examine the impact of the exact FDM linear power spectrum on the Lyman-$\alpha$ forest. We present the results of hydrodynamical simulations which make use explicitly of FDM transfer functions. We highlight the impact of the difference between WDM and FDM linear transfer functions on the predicted Lyman-$\alpha$ flux spectra. While we do not take into account the full wave behavior of the FDM fluid in the numerical simulations, we present quantitative arguments which suggest that this approximation should be valid at the scales probed by our simulations for $m_{a}\ga 10^{-22}$ eV. We then combine our models with existing high-statistics Lyman-$\alpha$ power spectra obtained from the BOSS survey, as published in Palanque-Delabrouille et al. (2013), to derive constraints on $m_{a}$, by applying a procedure similar to that used in e.g. Baur et al. (2016) to take into account various sources of uncertainties. Finally, we also add higher-resolution Lyman-$\alpha$ spectra from López et al. (2016) and Viel et al. (2013a) to the analysis in order to improve the sensitivity to higher values of $m_{a}$.

2 Calculation of the Lyman-$\alpha$ power spectrum

2.1 Linear matter spectrum

Throughout this paper we consider for simplicity a pure-FDM model, in which all of the dark matter is made of a (pseudo) scalar field $\Psi$, of mass $m_{a}$. We also do not consider possible self-interactions. In the scenario when such interactions are given by the expansion of a periodic axion-like potential $V\sim f_{a}^{2}\,m_{a}^{2}(1-\cos\Psi/f_{a})$, their effects are negligible at the cosmological scales of interest for values of $f_{a}$ that do not require fine-tuning to give the relic abundance. The homogeneous component $\Psi_{0}$ of the field obeys then the Klein-Gordon equation in the Friedman-Robertson-Walker metric:

[TABLE]

For the masses of interest here, $m_{a}\gg H(z)$ during all the epoch of structure formation so that the field oscillates and behaves on large scales as CDM. On small scales, the effect of the pressure associated with $\Psi$ gives to the FDM fluid an effective sound velocity (eg. Hwang & Noh (2009))

[TABLE]

The related Jeans scale damps structure formation in the linear regime at small scales. Therefore the transfer function defined from the linear matter power spectrum $P_{\rm FDM}(k)=T^{2}(k)\times P_{\rm CDM}(k)$ presents a strong attenuation at the comoving Jeans scale during equality. We calculate $P_{\rm FDM}(k)$ using the AxionCAMB (Hložek et al., 2017) modification of the CAMB Boltzmann solver (Lewis et al., 2000), with all cosmological parameters from Ade et al. (2016).

Fig. 1 illustrates the corresponding cutoff calculated for three axion masses in the range of interest for this study, namely $3\times 10^{-22}\la m_{a}\la 5\times 10^{-21}$ eV. We also show the FDM transfer function as parametrized by Hu et al. (2000). For identical values of $m_{a}$, we observe a slight, but significant shape difference between this parametrization and AxionCAMB : the power spectrum cutoff from the full AxionCAMB calculation is smoother than the one from Hu et al. (2000). This trend is confirmed by Ureña-López & Gonzalez-Morales (2016), a CLASS-based (Lesgourgues, 2011) calculation of the FDM transfer function. Based on AxionCAMB calculations for several values of $m_{a}$ in the range $10^{-23}-10^{-20}$ eV, we find the following scaling relation for the cutoff mode $k_{c}$ defined by $T(k_{c})=0.5$ :

[TABLE]

This relation is in good agreement with those found in Hu et al. (2000); Marsh & Silk (2014). From Eqn. 1, the mass-matching relation between WDM and FDM becomes

[TABLE]

again in reasonable agreement with previous estimates, including the order-of-magnitude value given in Eqn. 4.

The WDM transfer function, also represented on Fig. 1 for $m_{X}=2.5$ keV, is much smoother than the FDM one. Therefore it is expected that SDSS-based constraints on FDM will differ significantly from those obtained by the mapping of Eqn. 8, which justifies the need for dedicated simulations. Given that the Lyman-$\alpha$ forest essentially probes the one-dimensional matter power spectrum,

[TABLE]

the cutoff in the 3-D spectrum modeled by a given $T(k)$ will generically produce a cutoff in the 1-D spectrum with a different, smoother shape, and at much smaller values of $k$ (which explains why SDSS spectra, with $k\la 2\,h\,{\rm Mpc}^{-1}$, are sensitive to models with $k_{c}$ of several tens $h\,{\rm Mpc}^{-1}$). Non-linear and hydrodynamical effects make the picture more complicated, so that it is difficult to predict without full calculations how Lyman-$\alpha$ flux power spectra in the SDSS range are modified by changes in the shape of $T(k)$. It is nevertheless expected that SDSS Lyman-$\alpha$ data probe essentially the low-$k$ part of $T(k)$, and therefore constraints on $m_{a}$ should be less severe than those anticipated by making use of WDM bounds with relations such as Eqn. 8. Adding higher-resolution Lyman-$\alpha$ data will naturally extend the range of explored $k$, therefore improving the sensitivity to higher values of $m_{a}$.

2.2 Non linear matter spectrum

In the non-linear, newtonian and non-relativistic regime i.e. at low redshifts and reasonably small scale, the evolution of the FDM fluid can be rigorously followed using two equivalent methods : either by explicitly solving the Poisson-Schrödinger equations as in Schive et al. (2014), or by solving the corresponding Madelung fluid equations. Given the orders of magnitude for the DM density, velocity and de Broglie wavelength, the particle occupancy numbers are very high so that these equations can be interpreted in a classical way. For a FDM fluid whose wave function writes $\Psi=\sqrt{\rho}\,e^{i\theta}$, with $m\vec{v}=\hbar\nabla\theta$, the dynamical equation in a background geometry characterized by the scale factor $a$ and expansion rate $H$, with newtonian potential $\phi$, is (see e.g. Chavanis (2011)):

[TABLE]

This is equivalent to the Euler equation as implemented in standard N-body simulations, with an additional "quantum pressure" term $Q\equiv-(\hbar^{2}/2m_{a})\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}$. Note that, as Veltmaat & Niemeyer (2016) show, solving this fluid equation gives identical results to the Poisson-Schrödinger system away from interference fringes. Implementing this term into N-body simulations is difficult (Marsh, 2015; Mocz & Succi, 2015) mostly because it is a second order derivative of the density field, which means that most of its power, in terms of Fourier transform, is on very small scales. In Veltmaat & Niemeyer (2016), the effect of quantum pressure in the non-linear regime was recently studied at cosmological scales using a particle-mesh scheme. For $m_{a}=2.5\times 10^{-22}$ eV, significant modifications of the matter power spectrum are visible for $k\ga 100$ Mpc*-1*. This is consistent with arguments from Schive et al. (2016) : simulations with $m_{a}\sim 10^{-22}$ eV have shown that the $Q$ term has the most severe impact on the internal structures of DM halos, and much less on larger scales. However, the SDSS Lyman-$\alpha$ spectra do not probe comoving scales smaller than a few hundreds of kpc. As another argument, we computed with AxionCAMB the linear growth rate ratio $\xi$ defined in Schive et al. (2016) : for $m_{a}=3.4\times 10^{-22}$ eV, we find $\xi=1$ within 1 % for $k\la 20\,h\,{\rm Mpc}^{-1}$. We therefore expect the $Q$ term to have a negligible influence on the DM dynamics with respect to gravity, at the scale relevant for the observable Lyman-$\alpha$ forest.

As a consequence, for this study it is most probably appropriate to use standard N-body simulations. In Section 3.1, we will show a posteriori distributions for the $Q$ term which support this choice.

2.3 Hydrodynamical simulations

The methodology and numerical tools used here are very similar to and were extensively described in Borde et al. (2014); Rossi et al. (2014); Palanque-Delabrouille et al. (2015b); Baur et al. (2016). We use the Gadget-3 N-body simulation, an update of the Gadget-2 program (Springel et al., 2001; Springel, 2005). In this exploratory work, which aims to assess specifically the impact of the shape of the linear power spectrum on the resulting Lyman-$\alpha$ flux, we use a single set of $\Lambda$CDM cosmological parameters from Ade et al. (2014). They are in accordance to the central values used in Baur et al. (2016) : $h=0.675$, $\Omega_{M}=0.31$, $n_{s}=0.96$ and $\sigma_{8}=0.83$. Note that we checked that, as for WDM, the cutoff in the linear matter power spectrum induced by the FDM models considered here does not change significantly the value of $\sigma_{8}$, for a given primordial scalar perturbation amplitude.

The initial conditions are set at $z=30$ with the 2LPTIC software, starting from the linear matter power spectrum as computed by AxionCAMB for this redshift. The (fuzzy) dark matter fluid is then treated as a collection of fixed-mass point particles. The baryon fluid is evolved using the Smoothed Particle Hydrodynamics technique, with stars created in cold and dense baryon environments. Of importance for the Lyman-$\alpha$ forest, we model the IGM heating by the UV background light, using internal Gadget heating rate parameters which result in the redshift-dependent IGM temperature-density relation

[TABLE]

As for the case of cosmological parameters, we adopt here fixed benchmark parameters such that $T_{0}(z=3)=14000$ K and $\gamma(z=3)=1.3$, which are in agreement with the measurements from Becker et al. (2011).

From Gadget snapshot files in the redshift range $z=4.6-2.2$ adapted to SDSS Lyman-$\alpha$ data, we infer the line-of-sight-averaged one-dimensional Lyman-$\alpha$ flux power spectrum. This observable is defined from the fluctuations of quasar’s transmitted flux fraction, $\delta_{\varphi}(\lambda)=\varphi(\lambda)/\bar{\varphi}-1$, where $\bar{\varphi}$ is the mean transmitted flux fraction at the Hi absorber redshift, computed over the entire sample. We use here again a single H i optical depth model, which is known to roughly match existing data:

[TABLE]

In practice, the Lyman-$\alpha$ flux power spectrum is inferred in the range adapted either to published SDSS spectra, $k=10^{-3}-2\times 10^{-2}$ s km*-1*, or to higher-resolution spectra, $k=10^{-3}-0.1$ s km*-1*. We exploit the splicing technique as described in Borde et al. (2014). It consists in combining the results of two N-body simulation outputs, one of high resolution with $768^{3}$ DM particles in a 25 $h^{-1}\,$Mpc box, and one on larger scales, with $768^{3}$ particles in a $100\,h^{-1}\,$Mpc box, making use of a third low-resolution simulation with $128^{3}$ particles in a 25 $h^{-1}\,$Mpc box.

Simulations were computed for four different values of $m_{a}$ between $3.4\times 10^{-22}$ and $4.1\times 10^{-21}$ eV using AxionCAMB. To assess FDM-related systematic effects, an additional simulation was run with $T(k)$ given by the formula from Hu et al. (2000).

3 Results

Before comparing the predictions for the one-dimensional Lyman-$\alpha$ flux power spectrum with measured spectra, we first provide a discussion of the quantum pressure term, which was ignored in the N-body simulations. All the calculations presented here are therefore a posteriori and only hold if the dynamical impact of quantum pressure in the non-linear regime is negligible.

3.1 Quantum pressure

Fig. 2 illustrates the properties of the DM fluid, derived from the Gadget snapshots for $z=2.6$ at scales of a few Mpc, which correspond to the median redshift and smallest comoving scales of relevance for the SDSS Lyman-$\alpha$ flux. The top panel provides a by-eye comparison of the DM mass density for CDM with respect to the lowest-mass FDM used in the N-body simulations, $m_{a}=3.4\times 10^{-22}$ eV. The severe attenuation of small-scale structures due to the linear cutoff in the FDM scenario is evident.

In order to assess the relative importance of quantum pressure, the bottom panel compares the gravitational potential $\phi$ (left), as calculated explicitly with Gadget, with the quantum pressure $Q/m_{a}$ (right). Both are expressed in $(m/s)^{2}$. The $Q$ term is estimated from the numerical laplacian of the density field $\rho$, which is itself obtained by smoothing the DM point particle distribution with a kernel adapted to the local density of DM particles in the simulation, so that higher resolutions are obtained in higher density regions. We checked that the resulting $Q$ distributions are stable with respect to the kernel size parameter, which means that our estimation for $Q$ is not severely biased by shot-noise related fluctuations. On the other hand, we find that the gradient estimator is limited by the simulation resolution in regions where the DM density fluctuates rapidly, such as inside halos of course, but also along the filament edges : this means that in these regions higher-resolution simulations would necessarily predict larger values for $Q$. With this caveat in mind, we can observe that $Q$ is significant only in small regions of space, again along filaments and in the vicinity of DM halos. The fact that gradients of $Q$ look concentrated along filaments is consistent with the results of Veltmaat & Niemeyer (2016), who find that the main differences in DM density maps when fully taking into account the $Q$ term are in filaments (appart from halo cores which they don’t resolve). Note that, although it could be expected that $Q$ would blow up in deep voids due to the $1/\sqrt{\rho}$ factor (Marsh, 2015), we do not find any hint of this behavior from our calculations.

To quantify the validity of the N-body approximation used, we compute from Gadget snapshots the ratio between the two "force" terms which appear in the right-hand-side of Eqn. 10, based on our available estimate for $Q$:

[TABLE]

The mass-weighted distribution of this dimensionless ratio is represented on Fig. 3 for several values of $m_{a}$, using high-resolution snapshots at $z=2.6$.

For $m_{a}=1.8\times 10^{-21}$ eV (which is close to the limits that we will set in Section 3.3), and even for the lowest FDM mass implemented in our simulations, $m_{a}=3.4\times 10^{-22}$ eV, the quantum force at the scales resolved by our simulations is always much smaller than the gravitational force. For $m_{a}=3.4\times 10^{-22}$ eV, we find that $r_{Q}>0.01$ for much less than 1 % of the DM particles. The impact of $Q$ on the DM dynamics is therefore probably negligible. We observe a strong increase of $r_{Q}$ for higher redshifts due to the variation of gradients with the scale factor $a$, but that does not change the conclusion. However we remind again the caveat that our estimator for $Q$ is such that the quantum force is significantly underestimated in unresolved, high-gradient regions - and we know indeed from eg. Schive et al. (2014), that the quantum force does play a major role especially inside DM halos for such values of $m_{a}$.

For values of $m_{a}$ significantly lower than $\sim 10^{-22}$ eV, such as $m_{a}=5.0\times 10^{-23}$ eV, Fig. 3 shows the distribution of $r_{Q}$ obtained by a simple $1/m_{a}^{2}$ rescaling with respect to the case $m_{a}=3.4\times 10^{-22}$ eV. We see that a significant population of FDM particles is expected to undergo a "quantum force" more than 10 % times the gravitational force, mostly along filaments. This implies that the use of standard N-body simulations may be inappropriate in this mass regime, if one wants to predict density power spectra with high precision.

This discussion suggests that it is relevant to use classical N-body simulations in the context of this work, for scales resolved by our simulations (ie. larger than $\sim 100$ kpc), for FDM masses large enough, roughly $m_{a}\ga 10^{-22}$ eV. This is consistent with the results of Veltmaat & Niemeyer (2016), who find more explicitly that the effect of quantum pressure in the non-linear regime does not affect the power spectrum for modes $k\lesssim 100$ Mpc*-1* in the case $m_{a}=2.5\times 10^{-22}$ eV. Importantly, it is clear that the validity of this approximation is directly linked to the explored scales : the use of the relatively low-resolution SDSS Lyman-$\alpha$ flux power spectrum is therefore more robust against wave effects not accounted for, with respect to higher-resolution probes such as the XQ-100, HIRES and MIKE spectra.

3.2 Lyman-$\alpha$ flux power spectra

The Lyman-$\alpha$ flux power spectrum as calculated according to Section 2.3 at scales measured by SDSS is illustrated in Fig. 4, for several redshifts of interest, both in the case of CDM and FDM with $m_{a}=3.4\times 10^{-22}$ eV. The difference between the predicted fluxes is rather subtle, of the order of 5 %. This is however sufficiently large so that SDSS data can discriminate between models, using both the $k$ and $z$ signal dependences. On the contrary, for the highest FDM mass probed by our simulations, $m_{a}=4.1\times 10^{-21}$ eV, the typical differences in flux power with respect to CDM are below 1 % and therefore hard to test with only SDSS spectra.

Fig. 5 illustrates the impact of shape modifications in the linear transfer function $T(k)$ on the final Lyman-$\alpha$ flux prediction, for transfer functions with the same cutoff value $k_{c}$. While the differences between the linear power spectra predicted by AxionCAMB and the Hu formula are not negligible at all, as discussed in Section 2.1, the top panel of Fig. 5 shows that the related difference for the SDSS Lyman-$\alpha$ spectra is well below 1 %. Therefore we conclude that the uncertainties related to the modeling of the scalar field linear evolution are negligible for Lyman-$\alpha$ forest predictions.

On the contrary, the shape difference between FDM and WDM transfer functions do propagate significantly to the Lyman-$\alpha$ forest, as Fig. 5 (bottom) demonstrates. This is due to the fact that the WDM transfer function is smaller than the FDM one for values of $k$ much smaller than the cutoff $k_{c}$. Finally, let us observe from Fig. 5 that the flux power spectrum differences are more pronounced at high redshift, while they are reduced at low $z$ where non-linear effects tend to erase the original shape differences.

The Lyman-$\alpha$ flux power spectrum is also attenuated by the thermal pressure which smooths the spatial distribution of IGM over cosmological times. The relevant Jeans smoothing scale is given by (Schaye, 2001) $L_{J}\simeq 0.52\,{\rm kpc}\times(T/10^{4}{\rm K})^{1/2}\,(n_{H}/{\rm cm}^{-3})^{-1/2}$. Using the values for $T$ and $n_{H}$ calculated in our CDM simulation, the spatial distribution of $L_{J}$ is peaked around $L_{J}\sim 500$ kpc. The associated mode $2\pi/L_{J}$ corresponds to the value $m_{a}\simeq 7\times 10^{-22}$ eV, using Eqn. 7. The effect of thermal Jeans smoothing happens therefore at a scale comparable to the one associated to the FDM effect we’re looking for. There is therefore a degeneracy between IGM temperature parameters $(T_{0},\gamma)$ and $m_{a}$.

3.3 Constraining $m_{a}$ from the SDSS/BOSS Lyman-$\alpha$ forest

The SDSS Lyman-$\alpha$ forest data used here is the same as in Palanque-Delabrouille et al. (2013), in which all the computation details for the flux power spectrum and its related statistical and systematic uncertainties are extensively described. From a parent sample of $\sim 60,000$ SDSS-III/BOSS DR9 quasars (Ahn et al., 2012; Eisenstein et al., 2011; Ross et al., 2012), the data consist in a selection of $13,821$ spectra that have signal-to-noise ratios per pixel greater than 2, no broad absorption line features, no damped or detectable Lyman-limit systems, and an average resolution in the Lyman-$\alpha$ forest of at most $85\>\rm km~{}s^{-1}$. The Lyman-$\alpha$ forest is defined as the region spanning $1050<\lambda_{RF}/\r{A}<1180$. 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\times 10^{-2}~{}\rm s~{}km^{-1}$. The flux power spectrum is obtained from the Fourier Transform of the fractional flux transmission $\delta_{\varphi}$ defined in Section 2.3, computed separately for each $z$-sector.

It has been shown (eg. in Palanque-Delabrouille et al. (2013) that this SDSS Lyman-$\alpha$ forest matches well predictions based on the CDM model, by appropriately choosing relevant astrophysical and instrumental parameters such as the IGM thermal state modeling or the spectrometer resolution. We can therefore constrain the $m_{a}$ parameter of FDM in a way similar to constraints previously set on the sum of active neutrino masses, or the mass $m_{X}$ of WDM, using the same data. We use the same likelihood as used in Baur et al. (2016); Palanque-Delabrouille et al. (2015b), in order to handle the influence of relevant cosmological and astrophysical parameters : $H_{0}$, $\Omega_{M}$, $n_{s}$ and $\sigma_{8}$ on the one hand, and on the other hand the previously mentioned parameters $(T_{0},\gamma)$ (Eqn. 11) and $(\alpha,\beta)$ (Eqn. 12), which provide a simple description of the IGM temperature and mean Lyman-$\alpha$ optical depth. No CMB-based prior is set for the cosmological parameters, except for $H_{0}$, to which we apply the Gaussian constraint $H_{0}=67.3\pm 1$. No (ie. flat) priors are used for $(\alpha,\beta)$, while we constrain $\gamma=1\pm 0.3$ and $T_{0}\in[0;25000]$ K.

Other physical processes may impact significantly our FDM constraints, but are poorly known and not explicitly modeled in our simulations. To capture their effects, we introduce additional nuisance parameters, with simple analytical modifications of the predicted flux power spectra, based on other published simulations or models. The description of these corrections is detailed in Baur et al. (2016), and they include feedback processes from galactic outflows (Viel et al., 2013b), fluctuations of the UV background from discrete sources, and homogeneous reionization history (McDonald et al., 2005). Other identified nuisance parameters are taken into account, to model the Si iii-H i cross-correlation, the spectrometer noise, the presence of residual damped Lyman-$\alpha$ systems in the data, and simulation uncertainties related to the splicing algorithm. The signal dependence with $m_{a}$ is captured using the grid of simulations with $m_{a}=3.4\times 10^{-22},7.9\times 10^{-22},1.8\times 10^{-21}$ and $4.1\times 10^{-21}$ eV, and for which all other parameters are fixed to the benchmark values given in Section 2. For each FDM simulation a $\chi^{2}$ is computed with respect to the $35\times 12$ SDSS data points in $(k,z)$ space, assuming $h=0.673\pm 0.010$ (Ade et al., 2016) while floating all other likelihood parameters.

We then explored two methods to derive a bound on $m_{a}$ from these $\chi^{2}$s. In a first approach, we use already existing CDM simulations which vary all the aforementioned astrophysical and cosmological parameters while keeping $1/m_{a}=0$ in order to extrapolate the grid at any point of the whole parameter space. Since no FDM simulation including cross terms between $m_{a}$ and other parameters was carried out, this is done by assuming that the variation of the flux power spectrum with all parameters at any given value of $m_{a}$ is identical to the one calculated for $1/m_{a}=0$. We then identify the minimal $\chi^{2}$ value, letting all parameters free. It occurs for $\sigma_{8}=0.85$, $n_{s}=0.94$, $\Omega_{M}=0.29$, $(T_{0},\gamma)=(8723,0.93)$, $(\alpha,\beta)=(0.0025,3.7)$, and $1/m_{a}=0$ - meaning that the best fit is the pure CDM scenario. The minimal $\chi^{2}$ is 408.1, comparable to the one for the highest tabulated FDM mass, $\chi^{2}(m_{a}=4.1\times 10^{21}~{}{\rm eV})=409.0$. For the tabulated FDM model with lowest mass we find $\chi^{2}(m_{a}=3.4\times 10^{-22}~{}{\rm eV})=431.6$, which shows that such a value for $m_{a}$ is clearly excluded. A frequentist 95 % lower limit on $m_{a}$ is set, following eg. Palanque-Delabrouille et al. (2015a) : $m_{a}>2.0\times 10^{-21}$ eV. As already mentioned in Baur et al. (2016), all adjusted best-fit IGM nuisance parameters are compatible with no correction within $1\sigma$.

In order to by-pass the lack of FDM cross term calculations, and also to highlight the correspondance with WDM models, we also used a second approach which consists in using the WDM bound published in Baur et al. (2016) in the following way. We determine again a relationship between $m_{X}$ and $m_{a}$, this time so that the same $\chi^{2}$ is obtained in the SDSS adjustment, $\chi^{2}(m_{X})=\chi^{2}(m_{a})$. This requirement is satisfied with the following scaling law:

[TABLE]

This scaling is significantly different from Eqn. 8 which was derived from the cutoff position in linear matter spectra. However, it is the best adapted for our goal in the sense that the WDM and FDM models associated by this relation do predict Lyman-$\alpha$ flux power spectra that are indiscernible using current SDSS data. The 95 % CL limit $m_{X}>4.09$ keV derived by Baur et al. (2016) then converts to $m_{a}>2.3\times 10^{-21}$ eV. Given that this second approach takes into account cross terms, calculated in Baur et al. (2016), we consider it to be more accurate than the first one, and therefore conclude that the SDSS Lyman-$\alpha$ forest exclude FDM models with $m_{a}<2.3\times 10^{-21}$ eV at 95 % CL.

The above-mentioned matter spectral index adjusted with this SDSS data only, $n_{s}=0.94$, is in slight tension with the CMB-derived value $n_{s}=0.97$ (Ade et al., 2016). As was noticed in Baur et al. (2016), using SDSS only without Planck-based priors (except for $H_{0}$) is similar to allowing for a running of $n_{s}$ between the large scales of CMB and the smaller Lyman-$\alpha$ regime. Applying Planck-based priors on all our cosmological parameters therefore leads to a less tight constraint on the FDM or WDM mass. The 28 % change in $m_{X}$ derived in Baur et al. (2016) translates into the looser bound $m_{a}>1.3\times 10^{-21}$ eV if one includes Ade et al. (2016) into the adjustment.

Uncertainties linked to the modeling of FDM itself were already extensively discussed. We noticed that uncertainties in the linear matter spectrum do not modify the result. The impact of quantum pressure in the N-body simulations is more an open question, but we argued that at least for $m_{a}\ga 10^{-22}$ eV it should be negligible at the scale and precision probed by the SDSS Lyman-$\alpha$ forest.

Finally, errors related to the simplified treatment of IGM physics in our model were taken into account by allowing to float "nuisance" parameters, from $(T_{0},\gamma)$ to the analytical corrections for feedback effects or UV fluctuations. However, these corrections are derived from simulations or models with specific underlying assumptions, so that it cannot be certain that they fully take into account all processes relevant for our study. In addition, other phenomena that may impact the H i gas properties were not taken into account, such as residual fluctuations due to the inhomogeneity of the reionization process (Cen et al., 2009). This specific process should however affect more the high-redshift IGM properties.

3.4 Adding higher-resolution Lyman-$\alpha$ forest data

The analysis presented in Section 3.3 can easily be extended by including additional Lyman-$\alpha$ forest data with higher resolution. For that purpose, in this Section we add to the SDSS flux power spectrum the data from both the XQ-100 survey (López et al., 2016) and HIRES/MIKE data as presented in Viel et al. (2013a).

The XQ-100 survey, obtained with the XSHOOTER spectrograph at VLT, consists of a homogeneous and high-quality sample of 100 quasar spectra at $z=3.5-4.5$. The Lyman-$\alpha$ flux power spectrum was first derived from this survey in Iršič et al. (2017b), and the XQ-100 publicly available raw data was recently reanalyzed by Yèche et al. (2017), taking into account in particular an improved determination of the spectrograph resolution. Here we use the Lyman-$\alpha$ flux power spectrum as tabulated in Yèche et al. (2017), for three values of average redshifts in the range $z=3.2-3.93$, and spatial modes reaching values as high as $k=0.05-0.07$ s km*-1*, depending on the redshift.

In addition to XQ-100, we also include the Lyman-$\alpha$ flux power spectrum as presented in Viel et al. (2013a) (Fig. 7), which is based on a collection of quasar spectra from Keck/HIRES and Magellan/MIKE. These high-resolution spectrometers permit to reach spatial modes up to $k=0.09$ s km*-1*. Here we use only the flux power spectra determined at $z\leq 4.6$ due to the lack of simulation snapshots available at higher redshifts.

With the same methodology as in Section 3.3, we combine the SDSS Lyman-$\alpha$ spectra with those higher-resolution data, as was done in Yèche et al. (2017). We make use of the same simulated FDM flux power spectra, described in Section 2. With the procedure based on the $m_{X}-m_{a}$ mass matching, the limit $m_{X}>4.65$ keV from Yèche et al. (2017) translates into the following bound $m_{a}>2.9\times 10^{-21}$ eV at 95 % CL. This corresponds to a 30 % improvement with respect to the limit based on SDSS only : the higher resolutions of XSHOOTER, HIRES and MIKE permit to explore fluctuations at smaller comoving lengths.

If we include Planck priors on the cosmological parameters as was done with SDSS alone, the 95 % CL bound becomes $m_{a}>2.4\times 10^{-21}$ eV. The relative change in the constraint is less severe than with SDSS data alone. This is because the inclusion of higher-resolution spectra to the original SDSS data also reduces the above-mentioned tension on the $n_{s}$ parameter, as can be see in Yèche et al. (2017).

The uncertainties which may impact this bound are similar to those discussed in the end of Section 2.3. The impact of potential small-scale fluctuations in the IGM properties related to the distribution of ionizing sources and patchy reionization is increased, although our measurements do not rely on observations at $z\geq 5$. Also, for a given value of $m_{a}$, the effect of quantum pressure is stronger at the scales probed by high-resolution spectrographs : the "quantum force", being a third-order derivative, scales as $\sim k^{3}$. However, while this effect could be problematic for $m_{a}\sim 10^{-22}$ eV, it should not be a concern for the mass range $m_{a}=2.3-2.9\times 10^{-21}$ eV which is not covered by SDSS alone.

4 Conclusions

In this article, we studied the Lyman-$\alpha$ forest phenomenology within the Fuzzy Dark Matter scenario. The non-linear physics and hydrodynamics were modeled using standard N-body simulations with initial conditions which reproduce the FDM linear physics. This approach is an approximation which neglects the "quantum force" present in the Madelung equation. We presented semi-quantitative arguments, based on the simulation outputs, which strongly suggest that at least for $m_{a}\ga 10^{-22}$ eV the quantum force at the scales resolved by our simulations is always much smaller than the gravitational force. In agreement with the conclusions of Veltmaat & Niemeyer (2016), we find indications that for $m_{a}\la 10^{-22}$ eV the quantum term may impact the non-linear predictions for $P(k)$, and should therefore be included in the simulations in order to provide more reliable predictions.

Our main result is obtained by comparing the SDSS Lyman-$\alpha$ flux spectra and simulations, which permits to derived the limit $m_{a}>2.3\times 10^{-21}$ eV (95 % CL). Including higher-resolution spectra from the XSHOOTER, HIRES and MIKE spectrographs extends the exclusion range up to $m_{a}=2.9\times 10^{-21}$ eV. This analysis was done with a similar method, and similar assumptions, as for the case of WDM models (Baur et al., 2016; Yèche et al., 2017). For this study we also derived the most appropriate $m_{a}-m_{X}$ scaling law adapted to the SDSS Lyman-$\alpha$ data, which is different from the one which was previously found using the linear transfer function cutoffs.

These constraints confirm that the Lyman-$\alpha$ forest provides sensitivity to FDM masses as large as or larger than other observables such as high-redshift galaxy luminosity functions (Menci et al., 2017). More importantly, the Lyman-$\alpha$ forest is an independent probe with different sources of observational and astrophysical uncertainties. Together with others, our result now provides severe constraints on models tailored to solve small-scale problems of the CDM paradigm, for which $m_{a}\la 10^{-21}$ eV seems to be necessary.

While in the process of final redaction of this article, we were informed of the work presented in Iršič et al. (2017a). In their study, the authors of this work focus on the high-resolution Lyman-$\alpha$ forest to derive a bound on $m_{a}$. Their conclusions are very similar to ours in terms of constraints for $m_{a}$. In particular, our $m_{X}-m_{a}$ scaling law given by Eqn. 14 lies in between the curves quoted "$k_{1/2}$" and "$k_{0.75}$" in their Fig. 3, and interpolates reasonably their result points given by an MCMC analysis. We point out the complementarity in terms of systematics between both works. Whereas the work of Iršič et al. (2017a) uses exclusively high-resolution and high signal-to-noise Lyman-$\alpha$ spectra, including at redshifts above $z=5$, their statistical uncertainties are larger, and we note that some astrophysical uncertainties associated with the physics of the IGM and the effect of reionization are more important at high redshift. We also stress that the results of Iršič et al. (2017a), like ours, could be impacted by wave effects in the non-linear regime for $m_{a}\la 10^{-22}$ eV.

In addition to progress in the modeling of IGM physics, which is needed to improve the robustness of all Lyman-$\alpha$-based studies, future work should also be devoted to study the impact of non-linear wave effects for $m_{a}\la 10^{-22}$ eV. Baring this issue, our results provide among the best absolute lower bounds for the mass of bosonic dark matter. Even given this constraint, a number of interesting phenomena are predicted in the FDM framework that are distinctive from CDM and deserve to be explored in the future. More elaborate, but plausible scenarios such as mixed CDM-FDM models or models with mass-varying FDM are also worth investigating.

Acknowledgements

We acknowledge fruitful discussions with Jim Rich on the physics of FDM, and with Christophe Magneville on the N-body simulation pipeline. We thank Daniel Grin for exchanges on the AxionCAMB program. The authors are also grateful to Vid Irsic and Matteo Viel for providing the HIRES/MIKE-derived flux power spectra.

We thank Volker Springel for providing us with the Gadget-3 program. Gadget snapshots were exploited using the pynbody software (Pontzen et al., 2013). 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 Equipement National de Calcul Intensif) for access to the Curie cluster under allocation t2016047706. DJEM acknowledges support of a Royal Astronomical Society Postdoctoral Fellowship, hosted at King’s College London.