General framework for cosmological dark matter bounds using $N$-body simulations

We present a general framework for obtaining robust bounds on the nature of dark matter using cosmological $N$-body simulations and Lyman-alpha forest data. We construct an emulator of hydrodynamical simulations, which is a flexible, accurate and computationally efficient model for predicting the response of the Lyman-alpha forest flux power spectrum to different dark matter models, the state of the intergalactic medium (IGM) and the primordial power spectrum. The emulator combines a flexible parametrization for the small-scale suppression in the matter power spectrum arising in “noncold” dark matter models, with an improved IGM model. We then demonstrate how to optimize the emulator for the case of ultralight axion dark matter, presenting tests of convergence. We also carry out cross-validation tests of the accuracy of flux power spectrum prediction. This framework can be optimized for the analysis of many other dark matter candidates, e.g., warm or interacting dark matter. Our work demonstrates that a combination of an optimized emulator and cosmological “effective theories,” where many models are described by a single set of equations, is a powerful approach for robust and computationally efficient inference from the cosmic large-scale structure.

Figure 1

Transfer functions $T(k)$ [as defined in Eq. (1) at $z=99$] for various configurations of the noncold dark matter (nCDM) model as a function of comoving distance wave number $k$. The colored lines show some extremal values of our prior range, while the black lines show nCDM fits for cold dark matter ($\alpha =0$), 2 keV warm dark matter ($[\alpha ,\beta ,\gamma ]=[0.02,2.2,-4.5]$) and ${10}^{-22}\mathrm{eV}$ ultralight axion dark matter $[\alpha ,\beta ,\gamma ]=[0.1,5.2,-3.4]$.

Figure 2

The transfer function $T(k)$ [as defined in Eq. (1) at $z=99$] for the ultralight axion dark matter model as a function of comoving distance wave number $k$ and axion mass $m_a$. The dashed lines show $T(k)$ as calculated by the Boltzmann code axionCAMB [89]; the solid lines show our nCDM model fit using the $[\alpha ,\beta ,\gamma ]$ parameters. The nCDM model cannot characterize the small-scale oscillations in the axionCAMB transfer functions but in the power spectrum, these are strongly suppressed relative to the initial cutoff scale.

Figure 3

The nCDM model fit for ultralight axion dark matter transfer functions as a function of axion mass $m_a$. From top to bottom, the three nCDM parameters $[\alpha ,\beta ,\gamma ]$, showing a polynomial model fit (see text for more details).

Figure 4

The positions of the Bayesian-optimized ULA emulator training simulations projected onto each parameter plane. The shade of the crosses indicates the iteration at which the training point is added, such that the initial set is shaded lightest and the final optimization point is shaded darkest. We show only the planes where parameters are emulated together (i.e., not the planes with parameters from different redshift bins; or the ${\tau }_0(z=z_i)$ planes which are densely sampled by simulation postprocessing, see Sec. 2b). It follows that on the $x$ axis, $z_i$ corresponds to the redshift indicated on the $y$ axis. We note that for $\log (m_a[\mathrm{eV}])$, the initial set fully spans the $[\alpha ,\beta ,\gamma ]$ volume and does not project onto that axis; however, the initial set contributes to the final emulator model as we always emulate in $[\alpha ,\beta ,\gamma ]$. It can be seen how the darker-shaded optimization points concentrate in a smaller subvolume, which corresponds to the peak of the posterior distribution (see Fig. 6).

Figure 5

The convergence of the exploration term of the acquisition function (see Sec. 3b) as a function of optimization simulation number for the Bayesian-optimized ULA emulator.

Figure 6

The convergence of the posterior distribution (set of marginalized 1D and 2D distributions) using the Bayesian-optimized ULA emulator. Each set of colored contours shows the estimate of the posterior at various iterations of the emulator. The darker and lighter shaded regions, respectively, indicate the 68% and 95% credible regions. In each plane [except the ${\tau }_0(z=4.2)$ planes which are densely sampled by simulation postprocessing, see Sec. 2b], crosses indicate the projected positions of emulator training simulations. The crosses are color-coded by the iteration of the emulator by which the training points are added. We note that for $\log (m_a[\mathrm{eV}])$, the initial set fully spans the $[\alpha ,\beta ,\gamma ]$ volume and does not project onto that axis; however, the initial set contributes to the final emulator model as we always emulate in $[\alpha ,\beta ,\gamma ]$. $T_0(z=4.2)$ is in units of ${10}^4\mathrm{K}$; $u_0(z=4.2)$ is in units of $\mathrm{eV}{m}_p^{-1}$, where $m_p$ is the proton mass; and $m_a$ is in units of eV. For clarity, we show only the projections onto the IGM parameters at $z=4.2$; we observe similar convergence at the other redshifts we consider.

Figure 7

The convergence of summary statistics of the posterior distribution using the Bayesian-optimized ULA emulator. From top to bottom, the number of sigma shift (defined by the marginalized posteriors at a given optimization epoch) between emulator iterations for the 1D marginalized posterior means, $1\sigma$ and $2\sigma$ constraints. Each colored line shows the convergence for each model parameter.

Figure 8

A violin plot (set of kernel density estimation histograms) showing the ratio of empirical (emulator mean—truth) to predicted emulator error for leave-one-out cross-validation simulations, (from left to right) as a function of binned wave number $k$ and (from top to bottom) redshift $z$. The left set of PDFs considers only the simulations added by Bayesian optimization; the right set considers all simulations including the initial set. For comparison, we indicate $1\sigma$ and $3\sigma$ limits, respectively, by the darker and lighter shaded bands. If the emulator model was fitting perfectly, each distribution would be a unit Gaussian; the distribution being fatter indicates underfitting (underestimating emulator error) and the distribution being narrower indicates overfitting (overestimating emulator error). The wave number bins span the following ranges (where $k$ is in $h{\mathrm{Mpc}}^{-1}$): 1 [$0.6<k<10$]; 2 [$10<k<19$]; 3 [$19<k<28.3$].

Figure 9

A violin plot (set of kernel density estimation histograms) showing (the logarithm of) the ratio of theory to data error in the flux power spectrum, (from left to right) as a function of binned wave number $k$ and (from top to bottom) redshift $z$. For each element, the left PDF shows the ratio of emulator error to data error; for the majority of training data in the optimization set (which constitutes the lower-valued mode in each case), this is $<1$. Since Fig. 8 shows through cross-validation a tendency to overestimate the prediction uncertainty, the right PDFs show the ratio of the empirical error (for validation simulations only in the optimization set for clarity) to data error. This ratio is often $\ll 1$ (and less than the theoretical emulator-to-data error ratio, see also Fig. 8), indicating no statistically significant bias in emulation of the flux power spectrum. See more details in the text. The wave number bins span the following ranges (where $k$ is in $h{\mathrm{Mpc}}^{-1}$): 1 [$0.6<k<10$]; 2 [$10<k<19$]; 3 [$19<k<28.3$].

Figure 10

Cross-validation test on the inferred posterior distribution using the Bayesian-optimized ULA emulator: set of marginalized 1D and 2D distributions given mock simulated data removed from the training set. The darker and lighter shaded regions, respectively, indicate the 68% and 95% credible regions. The blue set of contours (obscured by the red set) uses a Gaussian prior on $T_0(z=z_i)$ with means matched to the true values of the mock data (and standard deviations of 3000 K): respectively at $z=[5.0,4.6,4.2]$, [10208, 9471, 10262] K. As a test of sensitivity to the prior distribution, the red set of contours instead uses $T_0(z=z_i)$ prior means set to the fiducial simulation values from our data analysis in Ref. [84]: respectively, at $z=[5.0,4.6,4.2]$, [8022, 7651, 8673] K. Dashed lines indicate the true parameter values of the mock data. $T_0(z=4.2)$ is in units of ${10}^4\mathrm{K}$; $u_0(z=4.2)$ is in units of $\mathrm{eV}{m}_p^{-1}$, where $m_p$ is the proton mass; and $m_a$ is in units of eV. For clarity, we show only the projections onto the IGM parameters at $z=4.2$; we observe similarly unbiased inference at the other redshifts we consider.

Figure 11

From top to bottom: a test of convergence in the simulated flux power spectrum at $z=5$ with respect to the number of simulation particles (for fixed comoving box volume); a test of convergence with respect to the comoving volume of the simulation box [for fixed particle mass resolution; the emulator is built with simulations with a comoving volume of $(10h^{-1}\mathrm{Mpc}{)}^3$ and $2\times 512^3$ particles]; a test of the systematic error arising from misestimating the spectral resolution of the data (up to $\pm 20\%$; the uncertainty on the spectral resolution is generally estimated to be $\sim 10\%$, e.g., [65]); and a test of the systematic error arising from ignoring the redshift evolution of the mean flux in sections of Lyman-alpha forest (the residual error in our analysis will be much less since data are corrected for this effect by the “rolling mean” estimator; e.g., [65]). For comparison, uncertainty in flux power spectrum data are $\sim 10\%$ to 25% depending on wave number (e.g., [65]).