Lyman-$\alpha$ forest constraints on Primordial Black Holes as Dark Matter

The renewed interest in the possibility that primordial black holes (PBHs) may constitute a significant part of the dark matter has motivated revisiting old observational constraints, as well as developing new ones. We present new limits on the PBH abundance, from a comprehensive analysis of high resolution, high redshift Lyman-$\alpha$ forest data. Poisson fluctuations in the PBH number density induce a small-scale power enhancement which departs from the standard cold dark matter prediction. Using a grid of hydrodynamic simulations exploring different values of astrophysical parameters, we obtain a marginalized upper limit on the PBH mass of $f_{\rm PBH}M_{\rm PBH} \sim 60~M_{\odot}\,(170~M_{\odot})$ at $2\sigma$ (depending on priors), which significantly improves previous constraints from the same physical observable. We also extend our predictions to non-monochromatic PBH mass distributions, ruling out large regions of the parameter space for some of the most viable PBH extended mass functions.

Another suitable and mostly unexplored method for constraining the mass and abundance of PBHs is offered by the Lyman-α forest, which is the main manifestation of the InterGalactic Medium (IGM).Such observable represents a very powerful tool for tracing the DM distribution at (sub-) galactic scales (see, e.g., [49][50][51]).
Lyman-α forest data was used almost two decades ago to set an upper limit of few 10 4 M on PBH masses, in the simple case in which all DM is made by PBHs with the same mass [52].In this work, we update and improve such limit, by using the highest resolution Lyman-α forest data set up-to-date [51], and a new set of high resolution hydrodynamic simulations that allow a more precise modeling of the full 1D flux power.Furthermore, we generalize our results to different PBH fractions with respect to the total DM amount, and to non-monochromatic mass distributions.
Poisson noise and impact on the matter power spectrum.Stellar-mass PBHs would cause observable effects on the matter power spectrum; due to discreetness, a small-scale plateau in the linear matter power spectrum is induced by a Poisson noise contribution [47,[52][53][54].
Let us firstly consider the simple case where PBHs are characterized by a Monochromatic Mass Distribution (MMD), i.e., they all have the same mass M PBH .In this scenario, PBHs are parameterized only by their mass and abundance, so that the fraction parameter f PBH ≡ Ω PBH /Ω DM = 1 in the case where the whole of the DM is made of PBHs.
Assuming that PBHs are randomly distributed, their number follows a Poisson distribution, and each wavenumber k can be associated to an overdensity δ PBH (k), due to Poisson noise.The PBH contribution to the power spectrum is thus defined as where n PBH is the comoving PBH number density, i.e. with ρ cr being the critical density of the universe.Since n PBH is a k-independent quantity, P PBH is scale-invariant.
In other words, as anticipated above, the existence of PBHs induces a small-scale plateau departing from the standard ΛCDM prediction.
A crucial concept pointed out in [52,54] is that the PBH overdensity δ PBH can be interpreted as an isocurvature perturbation (see also [55] for the case of a Lognormal spatial distribution of PBHs).Hence, the CDM power spectrum can be written as: where D(z) is the growth factor, P iso is the isocurvature power spectrum, and P ad ∝ A s k n s is the primordial adiabatic power spectrum.T ad and T iso are the adiabatic and isocurvature transfer functions, respectively (see, e.g., [56] for their analytic expressions).The PBH linear power spectrum is thus defined by: where we set the pivot scale k * = 0.05/Mpc, and the primordial isocurvature tilt n iso = 4 in order to ensure the scale-invariance of the power spectrum.Given that the adiabatic power spectrum evolves as k −3 at large k's, the isocurvature contribution is expected to become important only at the scales probed by the Lyman-α forest.
The parameter A iso sets the amplitude of the isocurvature modes, depending on the PBH mass of the considered model.Indeed, by combining the equations reported in this Section, we can eventually express the isocurvatureto-adiabatic amplitude ratio as: where the last equality holds only for MMDs.In this case it is immediate to see a degeneracy between f PBH and M PBH : different combinations of PBH mass and abundance correspond to the same isocurvature-to-adiabatic amplitude ratio if the quantity f PBH M PBH is the same.In our framework, the effect on the linear matter power spectrum due to the presence of isocurvature modes consists of a power enhancement with respect to the standard ΛCDM spectrum, in the form of a small-scale plateau.As a straightforward consequence of Equation (1), the larger is the PBH mass, the stronger is the small-scale power enhancement.
In Figure 1 we provide the relative differences with respect to a pure ΛCDM scenario for the 3D linear and non-linear matter power spectra, at redshift z = 5, for ΛPBH models with M PBH = {10 2 , 10 3 }M , assuming f PBH = 1.We also show the 1D flux power spectra, which are the Lyman-α forest observables, associated to the same ΛPBH models.The gray shaded area refers to the scales covered by our Lyman-α data set, obtained from MIKE/HIRES spectrographs.The nonlinear power spectra have been extracted from the snapshots of cosmological simulations, thus they include both the linear contribution (encoded in the initial conditions) and, on top of that, the effects of the non-linear evolution computed by the numerical simulation itself.The PBH contribution is thereby included in the initial conditions, whereas during the non-linear evolution both the isocurvature and adiabatic DM modes are treated as cold and collisionless (see Data set and methods Section for further details).It can be easily seen how non-linearities in the 3D matter power spectrum wash out the differences induced by the presence of PBHs.On the other hand, the 1D flux spectra are a much more effective observable to probe the small-scale power.
Extended Mass Distributions.The PBH formation is, in the most standard case, a consequence of large perturbations in the primordial power spectrum; while the exact details of the peak required to form PBH and how this is linked to the real-space overdensities are still unclear [57], it seems not unlikely that the most realistic scenario involves the presence of a PBH population with an extended mass function.Moreover, a non-monochromatic mass distribution would be created by different merger and accretion history of each PBH.Focusing on EMDs is also intriguing due to the possibility that PBHs in the high mass tail of the distribution may account for the seeds of supermassive BHs [47,58].A general method to convert MMD constraints to limits on Extended Mass Distributions (EMDs) has been developed in [44].
The extension to EMDs of the observable considered in this work arises naturally from the second equality in Equation ( 5), by directly considering the PBH number density corresponding to a given EMD.Let us then consider EMDs in the form: where the function dΦ PBH /dM PBH describes the shape of the EMD, and ρ DM = Ω DM ρ cr .Given a certain EMD, one can define the so-called Equivalent Mass M eq , which is the mass of a MMD providing the same observational effect (i.e., in our case, the same A iso ∝ f 2 PBH /n PBH ) [44].Therefore, the conversion from MMD(M eq ) to EMD(M PBH ) is given by (7) where we assume that the PBH abundances are the same for both the MMD and EMD cases.We finally have: Let us now estimate M eq for two popular EMDs: Lognormal and Powerlaw.
The Lognormal EMD is defined by where σ and µ are the standard deviation and the mean of the PBH mass, respectively.Such function is representative of a wide family of EMDs, since it describes well the scenario of PBHs forming from a smooth symmetric peak in the inflationary power spectrum [59,60].
The Powerlaw EMD is given by characterized by an exponent γ ∈ (−1, +1), a mass interval (M min , M max ), and a normalization factor N PL ; Θ is the Heaviside step function.The case with γ = 0 corresponds to the scenario of PBH formation from collapsing cosmic strings or scale-invariant density fluctuations [61].
Data set and methods.In order to extract limits on the properties of PBHs as DM candidates from the Lyman-α forest, we adapted the method proposed in [62] for testing suppressed linear matter power spectra with structure formation data.To do so, we have built a new grid of hydrodynamic simulations in terms of the properties of PBHs, corresponding to initial linear matter power spectra featuring a small-scale plateau.Beside that, our analyses rely on a pre-computed multidimensional grid of hydrodynamic simulations, associated to several values of the astrophysical and cosmological parameters affecting the Lyman-α flux power spectrum, sampling all the viable volume of the corresponding parameter space.The simulations have been performed with GADGET-III, a modified version of the publicly available numerical code GADGET-II [63,64].The initial conditions have been produced with the numerical code 2LPTic [65], at redshift z = 199, with input linear matter power spectra for the ΛPBH models obtained by turning on the isocurvature mode in the numerical Boltzmann solver CLASS [66].
For the cosmological parameters to be varied, we sample different values of σ 8 , i.e., the normalization of the linear matter power spectrum, and n eff = d ln P m (k)/d ln k | k α , namely the slope of the matter power spectrum evaluated at the scale probed by the Lyman-α forest (k α = 0.009 s/km) [69][70][71].Hence, we have included five different simulations for both σ 8 (in the range [0.754, 0.904]) and n eff (in the interval [−2.3474, −2.2674]).Additionally, we included simulations corresponding to different values for the instantaneous reionization redshift, i.e., z reio = {7, 9, 15}.
Regarding the astrophysical parameters, we modeled the thermal history of the IGM with amplitude T 0 and slope γ of its temperature-density relation, parameterized as T = T 0 (1 + δ IGM ) γ−1 , with δ IGM being the IGM overdensity [72].We ran simulations with temperatures at mean density T 0 (z = 4.2) = {6000, 9200, 12600} K, evolving with redshift, as well as a set of three values for the slope of the temperature-density relation, γ(z = 4.2) = {0.88,1.24, 1.47}.The redshift evolution of both T 0 and γ are parameterized as power laws, such that T 0 (z) = , where the pivot redshift z p is the redshift at which most of the Lyman-α forest pixels are coming from (z p = 4.5 for MIKE/HIRES).The reference thermal history is defined by T 0 (z = 4.2) = 9200 and γ(z = 4.2) = 1.47, since such values provide a good fit to observations [73].Furthermore, we have considered the effect of ultraviolet (UV) fluctuations of the ionizing background, the impact of which is controlled by the parameter f UV .Its template is built from a set of three simulations with f UV = {0, 0.5, 1}, where f UV = 0 corresponds to a spatially uniform UV background [67].We have also included 9 grid points obtained by rescaling the mean Lyman-α forest flux F(z), namely {0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4} × FREF , with the reference values being the ones of the SDSS-III/BOSS measurements [74].With the goal to have a more refined grid in terms of mean fluxes, we have also included 8 additional values, obtained by rescaling the optical depth τ = − ln F, i.e. {0.6, 0.7, 0.8, 0.9, 1.1, 1.2, 1.3, 1.4} × τ REF .
Concerning the PBH properties, we have extracted the flux power spectra from 12 hydrodynamic simulations (512 3 particles; 20 comoving Mpc/h box length) corresponding to the following PBH mass and fraction products: log(M PBH f PBH ) = {1.0,1.5, 2.0, 2.2, 2.3, 2.5, 2.6, 2.7, 3.0, 3.5, 4.0}1, where M PBH is expressed in units of M .For this set of simulations, all the astrophysical and cosmological parameters have been fixed to their reference values, and the equivalent ΛCDM flux power spectrum has also been determined.
Our analysis is based on an advanced interpolation method, i.e., the so-called Ordinary Kriging method [75], particularly suitable to deal with the sparse, non-regular grid defined by our simulations.The interpolation is done in terms of ratios between the flux power spectra of the ΛPBH models and the reference ΛCDM one.Along the lines of [62], we first interpolate in the astrophysical and cosmological parameter space for the standard ΛCDM case; we then correct all the (M PBH f PBH )-grid points accordingly, and we finally interpolate in the (M PBH f PBH )space.This procedure relies on the assumption that the corrections due to non-reference astrophysical or cosmological parameters are universal, so that we can apply the same corrections computed for the ΛCDM case to the ΛPBH models as well (see [62] for a detailed discussion on the accuracy of this approach).
Our data set is constituted by the MIKE/HIRES samples of quasar spectra, obtained with the MIKE/Magellan spectrographs and the HIRES/KECK, at redshift z = {4.2,4.6, 5.0, 5.4}, in 10 k-bins in the range [0.001 − 0.08] s/km, with spectral resolution of 6.7 and 13.6 s/km, respectively [51].As in the analyses of [62,67], we applied a cut on the flux power spectra, by using only the measurements at k > 0.005 s/km, to avoid large-scale systematic uncertainties due to continuum fitting.Moreover, we did not include in our analyses the highest redshift bin for MIKE data, for which the errors on the flux power spectra are very large [51].We have thus used a total of 49 (k, z) data points.
Results and Discussion.We have determined the constraints on both astrophysical, cosmological, and PBH properties, by maximising a Gaussian likelihood with a Monte Carlo Markov Chain (MCMC) approach, using the publicly available affine-invariant MCMC sampler emcee [76].We adopted Gaussian priors on the mean fluxes F(z), centered on their reference values, with standard deviation σ = 0.04 [67], and on σ 8 and n eff , centered on their Planck values [68], with standard deviation σ = 0.05, given that the latter two parameters, whereas very well constrained by CMB data, are poorly constrained by Lyman-α forest data alone [62].We adopt a logarithmic prior on f PBH M PBH .However, let us stress that our results are not affected by this choice.By rerunning our analyses assuming a flat prior on it, the constraints are indeed unmodified, provided that a conservative Gaussian prior on the IGM temperature is also imposed, namely T A 0 = 7500 ± 1500 (1σ) (motivated by the most up-to-date IGM studies by [77]).All the other  Flat prior on z reio Gaussian prior on z reio FIG.2: 1-and 2-σ contour plots for some of the parameters of our analyses, for the two different prior choices on z reio that we adopted.The values for M PBH are expressed in units of M .
priors are flat within the boundaries defined by our grid of simulations.
Let us firstly focus on the simple case of PBHs featuring a MMD.In Table I we report the marginalized 2σ constraints and the best fit values for all the parameters considered in our analyses.The first two columns refer to the case in which a flat prior is applied to the reion- ization redshift.The limit on the PBH abundance under this assumption corresponds to However, both Planck CMB results and measurements by [77] clearly favour z reio values around 8.5.For this reason, we have repeated our analysis by imposing on the reionization redshift a Gaussian prior centered around z reio = 8.5, with standard deviation σ = 1.0.The results obtained under such assumption are shown in the last two columns of Table I, and in this case we have In the simplest case where all DM is made by PBHs ( f PBH = 1), these constraints can be interpreted as absolute limits on the PBH mass.On the other hand, such absolute upper bounds on M PBH weaken linearly when a smaller PBH abundance is assumed (0.05 The degeneracy between z reio and the PBH mass can be understood as follows: a higher reionization redshift, having fixed the thermal history at late times, corresponds to a more effective (i.e., larger) filtering scale, and thus FIG.4: Overview on the present experimental constraints on the abundance of PBH for a monocromatic mass distribution (from [31,40,78]), in shades of blue, and in red are the limits from this work.
Patterned areas show limits that are the most dependent on astrophysical assumptions [79,80].The bottom panel is zoomed in the LIGO mass range.
to a suppression of power that is compensated by a larger value of the PBH mass.The fact that the degeneracies are much more prominent for this parameter, unlike the warm dark matter case in which the cutoff is more degenerate with the thermal cutoff, is telling us that the increase of power at small scales is a distinctive feature whose effect is more likely to be degenerate with a different gas filtering scale.
In Figure 2 we show the 1 and 2σ contour plots for some of the parameters considered in our analyses, for both the aforementioned prior choices on z reio .The mild degeneracy between the amplitude of the IGM temperature T A 0 (z = z p ) and the PBH mass derives from the opposite effects on the flux power spectra due to the increase of the two parameters.A hotter IGM implies a small-scale power suppression which can be balanced by increasing M PBH f PBH .Slightly larger values for the mean fluxes F(z) are also required for accommodating the power enhancement induced by relatively large values of the PBH mass.
In order to test the stability of our results, we also performed an analysis with flat both on σ 8 and n eff .Under these assumptions, the constraint (using a Gaussian prior for z reio ) on f PBH M PBH is mildly weakened, up to 100 M .However, the largest values for the PBH mass are allowed only in combination with extremely low values for n eff , allowed in turn by our data set due to its poor constraining power on such parameter.As we have already stated, this is the main reason to impose a (still conservative) Gaussian prior motivated by CMB measurements on n eff .
In Figure 3 we show the adimensional 1D flux power spectra for our best fit model (last column in Table I), together with the spectra corresponding to different values of M PBH .The latter ones are obtained by keeping fixed to their best fit values all other cosmological and astrophysical parameters.Symbols refer to MIKE/HIRES data.
In Figure 4 we report the updated plot with the constraints on the fraction of DM in PBHs, as functions of their masses, in the monochromatic case.The "LIGO window" of PBHs between ∼ 20 − 80M initially suggested in [15] has been probed and tentatively closed by constraints from Ultra-Faint Dwarf galaxies [35] and lensing of Type Ia Supernovae [31]; these constraints have been questioned because of astrophysics uncertainties (see, e.g., [79][80][81]3), and so we show those constraints in a patterned area.In this work we robustly close the higher mass part of that remaining window.There remains however, an interesting possibility in the very low mass range, 10 −10 M , see e.g., [78,82].
By defining an equivalent mass M eq one can easily convert the aforementioned marginalized 2σ limits on f PBH M PBH for the MMD case, to bounds on the parameters of a given EMD.In Figure 5 we provide such bounds, by plotting the parameters of the two EMDs discussed in this paper, as functions of M eq .In other words, each of the panels maps the limits on MMDs peaked at M eq to constraints on EMDs.The left panel shows the Powerlaw EMD, with γ = 0, focusing on the following mass range: M min , M max ∈ [10 −2 , 10 7 ] M .In the right panel we focus on the Lognormal EMD, scanning the parameter space defined by µ ∈ [10 −2 , 10 7 ] M , and σ ∈ [0, 5].The two black lines correspond to the marginalized 2σ contraints quoted above, i.e.M eq = 60 M (solid), and M eq = 170 M (dashed).The blue regions are therefore admitted by our analyses, while the red areas are ruled out.
Conclusions.In this work we have presented new bounds on the fraction of DM in PBHs, using an extensive analysis of high redshift Lyman-α forest data, improving over previous similar analyses in three different ways: 1) we used the high resolution MIKE/HIRES data, exploring better the high redshift range where primordial differences in the linear matter power are more prominent; 2) we relied on a very accurate set of high resolution hydrodynamic simulations which expands over a thermal history suggested by data; 3) we used the full shape of the 1D flux power rather than a single amplitude parameter.Our results improve previous constraints by nearly 3 orders of magnitude; furthermore, we have generalized our results to non-monochromatic PBH mass distributions, and ruled out a large part of the parameter space for two of the most popular EMDs, i.e., the Powerlaw and Lognormal EMDs.
In the near future, it is expected that a larger number of high redshift, high resolution and signal-to-noise quasar spectra collected with the ESPRESSO spectrograph [83] or at the E-ELT could further allow to use the small scale flux power spectrum to achieve tighter constraints.Another relevant aspect would be an accurate modeling of the heating and ionization due to accretion effects around the PBHs, to quantify how and if they could impact on the (much larger) scales of the Lyman-α forest.
The accumulation of limits on the PBHs as DM model in the mass range probed by LIGO seems to suggest that the hypothesis of 30M PBHs being the DM is less and less likely to be true.It has however become clear that these studies brought a plethora of astrophysical information, and even the exclusion of certain PBH mass ranges will bring information on some of the processes happening in the very early Universe.

TABLE I :
2σ limits and best fit values for all the parameters of our analyses, for the two different prior choices on z reio that we adopted.The values for M PBH are expressed in units of M .
FIG. 3: 1D flux spectra for the ΛCDM and ΛPBH cases, for different values of the PBH masses.Symbols are data from HIRES and MIKE, lines are obtained by interpolating in the (M PBH f PBH )-space defined by our grid of simulations; while the best fit is technically for M PBH 0, it is indistinguishable from the ΛCDM case.Red, blue, black and green indicate z = 4.2, 4.6, 5.0, 5.4, respectively.