On the road to percent accuracy II: calibration of the non-linear matter power spectrum for arbitrary cosmologies

TL;DR

This paper presents an emulator that predicts the non-linear matter power spectrum for various cosmologies beyond $\\Lambda$CDM using only standard $\\Lambda$CDM simulations, achieving high accuracy and enabling efficient testing of extended models.

Contribution

It introduces a novel emulator approach that leverages standard simulations and analytical models to accurately predict non-linear spectra for extended cosmologies, reducing computational costs.

Findings

Emulator achieves errors less than 1% in non-linear regime for tested cosmologies.

Method effectively models $f(R)$ gravity and massive neutrinos.

Reduces need for expensive beyond-$\\Lambda$CDM$ simulations.

Abstract

We introduce an emulator approach to predict the non-linear matter power spectrum for broad classes of beyond-$\Lambda$CDM cosmologies, using only a suite of $\Lambda$CDM $N$-body simulations. By including a range of suitably modified initial conditions in the simulations, and rescaling the resulting emulator predictions with analytical `halo model reactions', accurate non-linear matter power spectra for general extensions to the standard $\Lambda$CDM model can be calculated. We optimise the emulator design by substituting the simulation suite with non-linear predictions from the standard {\sc halofit} tool. We review the performance of the emulator for artificially generated departures from the standard cosmology as well as for theoretically motivated models, such as $f (R)$ gravity and massive neutrinos. For the majority of cosmologies we have tested, the emulator can reproduce the matter power spectrum with errors $\lesssim 1\%$ deep into the highly non-linear regime. This work demonstrates that with a well-designed suite of $\Lambda$CDM simulations, extensions to the standard cosmological model can be tested in the non-linear regime without any reliance on expensive beyond-$\Lambda$CDM simulations.

Tables

Table 1. Table 1 : The massive neutrino, f ​ ( R ) 𝑓 𝑅 f(R) gravity and hybrid models used in this work to design and test our emulator scheme, each labelled according to the sum of neutrino masses ( ∑ m ν subscript 𝑚 𝜈 \sum m_{\nu} ) and/or the strength of the deviation from standard gravity ( | f R ​ 0 | subscript 𝑓 𝑅 0 |f_{R0}| ). We consider two base Λ Λ \Lambda CDM cosmologies: (A) Ω b ​ h 2 = 0.0225 , Ω m ​ h 2 = 0.1382 , h = 0.6898 , n s = 0.969 , ln ⁡ ( 10 10 ​ A s ) = 3.195 formulae-sequence subscript Ω b superscript ℎ 2 0.0225 formulae-sequence subscript Ω m superscript ℎ 2 0.1382 formulae-sequence ℎ 0.6898 formulae-sequence subscript 𝑛 s 0.969 superscript 10 10 subscript 𝐴 s 3.195 \Omega_{\rm{b}}h^{2}=0.0225,\Omega_{\rm{m}}h^{2}=0.1382,h=0.6898,n_{\rm s}=0.969,\ln(10^{10}A_{\rm s})=3.195 ; and (B) Ω b ​ h 2 = 0.0173 , Ω m ​ h 2 = 0.0864 , h = 0.6 , n s = 0.969 , ln ⁡ ( 10 10 ​ A s ) = 3.807 formulae-sequence subscript Ω b superscript ℎ 2 0.0173 formulae-sequence subscript Ω m superscript ℎ 2 0.0864 formulae-sequence ℎ 0.6 formulae-sequence subscript 𝑛 s 0.969 superscript 10 10 subscript 𝐴 s 3.807 \Omega_{\rm{b}}h^{2}=0.0173,\Omega_{\rm{m}}h^{2}=0.0864,h=0.6,n_{\rm s}=0.969,\ln(10^{10}A_{\rm s})=3.807 . Note that (B) sits several standard deviations away from the  Planck Collaboration et al. ( 2018 ) best fit cosmology, thus working as stress test for our method.

model	$\sum m_{\nu }$	$\left|f_{R0}\right|$	base parameters
MNU_0.05	$0.05\mathrm{eV}$	–	(A)
MNU_0.1	$0.1\mathrm{eV}$	–	(A)
MNU_0.2	$0.2\mathrm{eV}$	–	(A)
MNU_0.4	$0.4\mathrm{eV}$	–	(A)
xMNU_0.05	$0.05\mathrm{eV}$	–	(B)
xMNU_0.1	$0.1\mathrm{eV}$	–	(B)
xMNU_0.2	$0.2\mathrm{eV}$	–	(B)
xMNU_0.4	$0.4\mathrm{eV}$	–	(B)
F4	–	${10}^{-4}$	(A)
F5	–	${10}^{-5}$	(A)
F6	–	${10}^{-6}$	(A)
F4-MNU_0.4	$0.4\mathrm{eV}$	${10}^{-4}$	(A)
F5-MNU_0.2	$0.2\mathrm{eV}$	${10}^{-5}$	(A)
F6-MNU_0.1	$0.1\mathrm{eV}$	${10}^{-6}$	(A)

Equations

$${\lambda }_{\mathrm{C}0}\approx 30\sqrt{\frac{(n+1)}{4-3{\Omega }_{\mathrm{m}}}\frac{|f_{R0}|}{10^{-4}}}{h}^{-1}\mathrm{Mpc},$$

$${\lambda }_{\mathrm{FS}}^{\max }\approx 350\sqrt{\frac{1}{{\Omega }_{\mathrm{m}}}\frac{1\mathrm{eV}}{m_{\nu }}}{h}^{-1}\mathrm{Mpc}.$$

$${\Omega }_{\nu }{h}^2=\frac{\sum m_{\nu }}{93.14\mathrm{eV}}.$$

$$P_{\mathrm{NL}}^{\mathrm{real}}(k,z)=\mathcal{R}(k,z)\times P_{\mathrm{NL}}^{\mathrm{pseudo}}(k,z),$$

$$P_{\mathrm{L}}^{\mathrm{pseudo}}(k,z)=P_{\mathrm{L}}^{\mathrm{real}}(k,z).$$

$$\mathcal{R}(k,z)\equiv {\frac{P_{\mathrm{NL}}^{\mathrm{real}}(k,z)}{P_{\mathrm{NL}}^{\mathrm{pseudo}}(k,z)}|}_{\mathrm{HM}},$$

$$B(k,z)\equiv \frac{P_{\mathrm{NL}}^{\mathrm{pseudo}}(k,z)}{P_{\mathrm{L}}^{\mathrm{pseudo}}(k,z)}.$$

$$\ln B(k,z;\mathbf{\pi })=\mu (k,z)+\sum _{j=1}^{n_{\psi }}{\psi }_j(k,z)w_j(\mathbf{\pi })+ϵ,$$

$$\begin{array}{rl}S(k,z;{\mathbf{\pi }}_{\Lambda },\Delta \mathbf{\alpha })& \equiv \frac{P_{\mathrm{L}}^{\mathrm{pseudo}}(k,z;{\mathbf{\pi }}_{\Lambda },{\mathbf{\pi }}_{\mathrm{X}})}{P_{\mathrm{L}}^{\Lambda }(k,z;{\mathbf{\pi }}_{\Lambda })}\\ & \approx 1+\sum _{i=1}^{n_{\Phi }}{\Phi }_i(k,z)\Delta {\alpha }_i,\end{array}$$

$$S(k_m)\approx \nu (k_m)+\sum _{i=1}^{n_{\Phi }}{\Phi }_i(k_m){\alpha }_i,$$

$$1-\nu (k_m)\approx \sum _{i=1}^{n_{\Phi }}{\Phi }_i(k_m){\alpha }_i^{\Lambda },$$

$${\alpha }_i^{\Lambda }=\sum _{m=1}^{n_k}{\Phi }_i(k_m)[1-\nu (k_m)].$$

$$\underset{\mathbf{D}}{\min }E(\mathbf{D})={\left(\frac{1}{\left(\genfrac{}{}{0ex}{}{N}{2}\right)}\sum _{i=1}^{N-1}\sum _{j=i+1}^N\frac{1}{{\prod }_{l=1}^d(x_{il}-x_{jl}{)}^2}\right)}^{\frac{1}{d}},$$

$$\frac{\partial E^d(\mathbf{D})}{\partial x_{rs}}=\frac{2}{\left(\genfrac{}{}{0ex}{}{N}{2}\right)}\sum _{i\ne r}\frac{1}{(x_{is}-x_{rs})}\frac{1}{{\prod }_{l=1}^d(x_{il}-x_{rl}{)}^2}.$$

$$\begin{array}{c}0.0215<{\Omega }_{\mathrm{b}}{h}^2<0.0235,\\ 0.12<{\Omega }_{\mathrm{m}}{h}^2<0.155,\\ 0.6<h<0.8,\\ 0.9<n_s<1.05,\\ 2.92<\log (10^{10}{A}_s)<\mathrm{3.16.}\end{array}$$

$${ϵ}^{\mathrm{emu}}\equiv \underset{0.01<k<10}{\max }\left|\frac{P_{\mathrm{NL}}^{\mathrm{emu}}(k)-P_{\mathrm{NL}}^{\mathrm{True}}(k)}{P_{\mathrm{NL}}^{\mathrm{True}}(k)}\right|$$

$$y(\mathbf{\pi })=f(\mathbf{\pi })+{ϵ}_{\mathrm{n}}(\mathbf{\pi }),$$

$$K(f,f^{\ast };\mathbf{h})\equiv \mathrm{cov}\left(f(\mathbf{\pi }),f({\mathbf{\pi }}^{\ast });\mathbf{h}\right)=A\prod _{l=1}^d\exp \left[\frac{({\pi }_l-{\pi }_l^{\ast }{)}^2}{p_l^2}\right].$$

$$S^{\mathrm{dw}}(k,z)=\frac{P_{\mathrm{L}}^{\mathrm{pseudo},\mathrm{dw}}(k,z)}{P_{\mathrm{L}}^{\Lambda ,\mathrm{dw}}(k,z)},$$

$$\mathcal{W}(k,z)=S(k,z)-S^{\mathrm{dw}}(k,z).$$

$$\tilde{\mathcal{W}}(k;a,b)=[a+b{\log }_{10}(k)]{\mathcal{W}}_{\mathrm{fid}}(k),$$

$$S(k;{\mathbf{\pi }}_{\Lambda },\Delta \mathbf{\alpha },a,b)=S^{\mathrm{dw}}(k;{\mathbf{\pi }}_{\Lambda },\Delta \mathbf{\alpha })+\tilde{\mathcal{W}}(k;a,b),$$

$$N_{\mathrm{eff}}=3.046+\Delta N_{\mathrm{eff}}.$$

Full text

1 Introduction

Large cosmological data sets have ushered in the era of “precision cosmology”, in which parameters describing the $\Lambda$CDM model are known with uncertainties of only a few percent. This increased precision has recently brought to light some level of discordance between early- and late-time cosmological probes. In particular, all three major weak lensing surveys (Hildebrandt et al., 2017; Hikage et al., 2018; Troxel et al., 2018) infer lower values of $\sigma_{8}\sqrt{\Omega_{\rm{m}}}$ compared to those obtained from observations of the temperature and polarisation of the cosmic microwave background (CMB; Planck Collaboration et al., 2018). Significant tensions have also arisen between the value of the Hubble constant derived from the CMB and its direct local measurement via the distance-redshift relation calibrated with Cepheid variables and Type Ia supernovae (Riess et al., 2016, 2018). If not the product of an unlikely statistical fluctuation, then these discrepancies may be ascribed either to unaccounted systematic uncertainties in the analyses or to an incompleteness of the concordance cosmology (see, e.g., Verde et al., 2013; Joudaki et al., 2017; Mörtsell & Dhawan, 2018; DES Collaboration et al., 2018).

In the event that the data prefer a revision of the standard cosmological model, accurate predictions of large-scale structure statistics for alternative cosmologies to $\Lambda$CDM become very valuable. However, we currently lack a method to predict these observables for general extensions to $\Lambda$CDM and over the range of scales relevant for ongoing and future galaxy surveys (Laureijs et al., 2011; LSST Dark Energy Science Collaboration, 2012; Hildebrandt et al., 2017; Abbott et al., 2018; Taylor et al., 2018). This means we cannot take advantage of the wealth of information contained in the non-linear scales without resorting to computationally expensive numerical simulations.

It is however infeasible to simulate all possible combinations of cosmological parameters, even within the $\Lambda$CDM paradigm. For current applications, fitting functions calibrated on a small number of $N$-body simulations are good enough to obtain unbiased constraints (Smith et al., 2003; Takahashi et al., 2012; Mead et al., 2015). On the other hand, future all-sky surveys will place more stringent requirements on our theoretical predictions, which makes emulator techniques a preferable strategy as far as accuracy is concerned. The basic idea is to interpolate between the output of $\sim$10–100 $N$-body simulations of cosmologies associated with carefully-chosen input parameters (Habib et al., 2007; Schneider et al., 2008; Heitmann et al., 2009, 2016; Lawrence et al., 2017; Euclid Collaboration et al., 2018; Winther et al., 2019; Harnois-Deraps et al., 2019). Previous implementations of this methodology have been applied either to flat, massless neutrinos (i.e. vanilla) $\Lambda$CDM, or to a few specific extensions, that is models with massive neutrinos, evolving dark energy with a linearly parameterised equation of state or $f(R)$ gravity. The computational expense of simulating modified gravity theories, however, can be orders of magnitude larger than that of the standard model (Llinares, 2018). Consequently, there currently exists no emulator capable of making non-linear predictions accurate at the level of $\sim$1% across a broad range of beyond-$\Lambda$CDM cosmologies. This severely limits our capacity for constraining these alternative models with data probing the growth of structure on cosmic scales.

In this paper we complement the work of Cataneo et al. (2019, referred to as ‘C19’ henceforth), who formulated an efficient method to compute the non-linear matter power spectrum in a broad range of extensions to the standard cosmology with the desired percent accuracy. The C19 strategy consists of rescaling the non-linear power spectrum of a tailored $\Lambda$CDM-like model – the pseudo cosmology – by an analytical function encapsulating the non-linear physics beyond the vanilla cosmology – the halo model reaction. By construction the pseudo cosmology shares the linear power spectrum of the real beyond-$\Lambda$CDM cosmology.

For $\Lambda$CDM extensions exhibiting a scale-independent linear growth of structure, the corresponding pseudo non-linear power spectrum can be obtained from a $\Lambda$CDM emulator by matching the amplitude of mass fluctuations, $\sigma_{8}$, to that of the real cosmology at the redshift of interest. In general, however, structures in alternative cosmological scenarios grow at different rates on different scales, which makes traditional emulators unfit for this task. Here, we present a new Gaussian process emulator designed to compute the pseudo non-linear matter power spectrum of arbitrary cosmologies, including both scale-independent and scale-dependent linear growth. We show that a suite of $\sim$1000 $\Lambda$CDM-like simulations is enough to predict with percent accuracy the pseudo non-linear power spectrum of popular cosmologies which have scale-dependent linear growth, such as $f(R)$ gravity and massive neutrino models. We also obtain similar performance for synthetic cosmological models created by rescaling the linear $\Lambda$CDM power spectra by smooth arbitrary functions.

The combination of the C19 halo model reactions with the emulator developed in this work will hence provide accurate predictions of the non-linear matter power spectrum in a broad class of cosmological models, including modified gravity, dark energy and massive neutrinos, for use in cosmic shear analyses (see, e.g., Troxel et al., 2018; Hikage et al., 2018; Hildebrandt et al., 2018). The fact that galaxies are biased tracers of the underlying matter distribution means that additional emulator complexity, including parameters to describe the halo occupation distribution, would be required for this methodology to be used in future clustering analyses (see, e.g., Zhai et al., 2019). There is also the potential to use this technique to model the impact of baryonic feedback.

This paper is organised as follows. In Section 2 we examine the expected range of deviations in the linear matter power spectrum of viable extensions to the $\Lambda$CDM cosmology, and outline the C19 reaction framework for modelling the non-linear power spectrum. In Section 3 we detail the emulator methodology for arbitrary input linear power spectra, where we define ‘arbitrary’ to mean any power spectra which deviates smoothly from the standard model approximately within the range allowed by observational constraints. We also describe the construction of our suite of surrogate numerical simulations used as a training set for the emulator, which we design to encompass a reasonable range of $\Lambda$CDM extensions. In Section 4 we present the performance of our emulator and determine the minimum number of simulations we need to reach the accuracy required by the high-quality data from the next generation of all-sky imaging surveys. Finally, we conclude and discuss the impact of our findings in Section 5.

2 Beyond vanilla $\Lambda$CDM

The late-time phenomenology of the concordance model – a globally flat Universe with a matter-energy content dominated by cold dark matter, baryonic matter and the cosmological constant – is entirely described by five base parameters, $\boldsymbol{\pi}_{\Lambda}\equiv\{\Omega_{\rm{b}}h^{2},\Omega_{\rm{m}}h^{2},h,n_{\rm s},\ln(10^{10}A_{\rm s})\}$, where $\Omega_{\rm{b}}$ is the current fraction of energy-density in baryonic matter, $\Omega_{\rm{m}}$ is the present-day total matter energy-density fraction, $h$ is the dimensionless Hubble constant, and $n_{\rm s}$ and $A_{\rm s}$ are, respectively, the slope and amplitude of the primordial scalar power spectrum. Under the assumption of the standard cosmology, measurements of the CMB constrain the radiation content (i.e. photons and three massless neutrino species) to better than one part in a thousand (Fixsen, 2009).

Changes to the dark sector and to the law of gravity have profound implications for the formation of structures in the Universe over a wide range of scales. Departures from the standard cosmology, therefore, leave distinctive features on the statistical quantities measuring the clustering of matter, such as the two-point correlation function, or its Fourier transform, the power spectrum $P(k)$. The landscape of $\Lambda$CDM extensions is exceptionally vast, and in this work we only consider a family of $f(R)$ gravity theories (Hu & Sawicki, 2007) and massive neutrino cosmologies (Lesgourgues & Pastor, 2006), both with a $\Lambda$CDM background. However, with little or no modification our methodology (see Section 3) can be applied to a much broader class of models, including non-standard dark matter candidates  (Marsh & Silk, 2014; Schneider, 2015; Cyr-Racine et al., 2016; Marsh, 2016; Poulin et al., 2016; Hložek et al., 2017; Dakin et al., 2019; Thomas et al., 2019), extra relativistic degrees of freedom (Baumann et al., 2018), and Horndeski theories (Zumalacárregui et al., 2017).

What makes models like $f(R)$ gravity and massive neutrino cosmologies interesting is their scale-dependent linear growth of structure. More specifically, $f(R)$ theories enhance the growth on scales comparable to, or smaller than, the Compton wavelength, which at $z=0$ reads (Hu & Sawicki, 2007)

[TABLE]

where $f_{R0}$ defines the extent of the departure from General Relativity, recovered for $f_{R0}=0$. We set $n=1$ and use $|f_{R0}|=10^{-4},10^{-5},10^{-6}$ (referred to as F4, F5 and F6 in the following), which brackets the range of values giving rise to interesting cosmological behaviours, with F4 already strongly disfavoured by current data (Terukina et al., 2014; Lombriser, 2014; Cataneo et al., 2015; Liu et al., 2016; Alam et al., 2016).

On the other hand, the high thermal velocities of massive neutrinos prevent them from clustering on scales smaller than the free-streaming length, $\lambda_{\rm FS}$. This results in the suppression of the growth of structure on scales smaller than the maximum free-streaming length, which is defined at the time of the non-relativistic transition (after recombination for $m_{\nu}\lesssim 0.5\,{\rm eV}$), and is given by (Lesgourgues & Pastor, 2006)

[TABLE]

For simplicity, we will consider cosmologies including two massless neutrinos and one massive neutrino with $m_{\nu}=0.05,0.1,0.2,0.4\,{\rm eV}$, where the smallest value is close to the minimum sum of neutrino masses consistent with neutrino oscillation experiments (Forero et al., 2014), and the largest value is outside the 95% confidence region found by Planck Collaboration et al. (2018). In addition, for all these extensions we keep the base cosmological parameters fixed to their $\Lambda$CDM values, and compensate for the presence of massive neutrinos by reducing the cold dark matter density as $\Omega_{\rm c}\rightarrow\Omega_{\rm m}-\Omega_{\rm b}-\Omega_{\nu}$, where

[TABLE]

Table 1 lists all models, collectively referred to as physical models, used to test our emulation scheme described in Section 3, including hybrid cosmologies with massive neutrinos and modified gravity. For these extensions, the linear power spectrum deviations from the standard predictions are obtained with mgcamb111http://aliojjati.github.io/MGCAMB/mgcamb.html(Zhao et al., 2009; Hojjati et al., 2011). Figure 1 shows a $z=0$ example of beyond-$\Lambda$CDM to $\Lambda$CDM linear matter power spectrum ratios.

2.1 The halo model reaction framework

By unlocking the information stored in the statistical distribution of matter on small scales we can substantially increase our sensitivity to physics beyond the vanilla cosmology, be it new particles/fluids, modifications to the theory of gravity, or astrophysics (see, e.g. Huterer et al., 2015; Copeland et al., 2018; Nori & Baldi, 2018; Heymans & Zhao, 2018; Schneider et al., 2019). Not surprisingly, then, the non-linear matter power spectrum is found at the core of all cosmological analyses hinged on galaxy survey data  (Parkinson et al., 2012; Kitching et al., 2014; Gil-Marín et al., 2016; Hildebrandt et al., 2017; van Uitert et al., 2018; Joudaki et al., 2018; Abbott et al., 2018; Gil-Marín et al., 2018; Hikage et al., 2018). To take full advantage of the exquisite observations from forthcoming experiments (Laureijs et al., 2011; LSST Dark Energy Science Collaboration, 2012; Levi et al., 2013), theoretical predictions must reach percent level accuracy deep in the non-linear regime of structure formation (Taylor et al., 2018). At present, however, no method is at the same time computationally efficient, accurate and flexible enough to be employed in future analyses aimed at constraining physics beyond the standard paradigm  (Lesgourgues et al., 2009; Bird et al., 2012; Heitmann et al., 2014; Brax & Valageas, 2013; Zhao, 2014; Blas et al., 2014; Massara et al., 2014; Levi & Vlah, 2016; Lawrence et al., 2017; Aviles & Cervantes-Cota, 2017; Senatore & Zaldarriaga, 2017; Bose et al., 2018; Euclid Collaboration et al., 2018; Heisenberg & Bartelmann, 2019; Winther et al., 2019).

The halo model reaction framework first proposed by C19, and based on the work of Mead (2017), stands out as a promising, practical solution to the long-standing problem of predicting the non-linear matter power spectrum with the required accuracy and speed, while also being readily applicable to a variety of cosmological scenarios. Within this approach, the non-linear power spectrum of a given real cosmology takes the form

[TABLE]

where the pseudo cosmology is defined as a $\Lambda$CDM cosmology sharing the same linear matter power spectrum of the real cosmology, that is

[TABLE]

The reaction,

[TABLE]

captures the non-linearities sourced by late-time physics beyond $\Lambda$CDM, and is predicted with the halo model (HM) and perturbation theory (see Cataneo et al., 2019, for details). In taking the ratio of the two halo model predictions, Mead (2017) showed that the impact of the known inaccuracies in the individual power spectra is significantly reduced. On linear scales and for real $\Lambda$CDM cosmologies one has the trivial relation $\mathcal{R}=1$. C19 found that the modelling approach summarised by equation 4 is accurate at the percent level for a broad range of general extensions to the standard model.

Here, we formulate a method to provide predictions for $P_{\rm NL}^{\rm pseudo}$ that are more accurate than any semi-analytical prescription based on the halo model. It is this accurate estimate that is then used to calibrate the reaction, $\mathcal{R}$, in equation 4 in order to derive the non-linear matter power spectrum of our target beyond-$\Lambda$CDM cosmology, $P_{\rm NL}^{\rm real}$. To this end, we develop a novel emulation scheme where in addition to the base $\Lambda$CDM parameters we use the shape of the linear matter power spectrum as input. This work is intended as a first feasibility study, and as such we approximate the output of otherwise expensive pseudo $N$-body simulations with the halofit fitting functions (Takahashi et al., 2012) implemented in a modified version of the public Einstein-Boltzmann solver camb222https://github.com/cmbant/CAMB (Lewis et al., 2000). In practice, for a model $X$, camb reads in both $\boldsymbol{\pi}_{\Lambda}$ and the externally generated $P_{\rm L}^{X}/P_{\rm L}^{\Lambda}$ ratio, evaluates $P_{\rm L}^{\rm pseudo}$ and gives it to halofit, which finally provides the desired non-linear quantity. The grey curves in the upper panel of Figure 2 show the pseudo non-linear power spectra for the physical models in Table 1, where we used the set (A) of base $\Lambda$CDM parameters also for the xMNU cosmologies. Note that despite halofit predictions being approximate, this strategy is theoretically consistent, in that pseudo cosmologies are simply $\Lambda$CDM cosmologies with non-standard initial conditions. Our approach enables a rapid design, construction and performance assessment of the emulator (Heitmann et al., 2009; Euclid Collaboration et al., 2018).

3 Methodology

In this Section, we describe the basics of Gaussian process emulation of cosmological statistics. We also detail our pipeline, which starts from the linear matter power spectrum of an arbitrary cosmology and maps this onto a point in the multi-dimensional parameter space over which the emulator performs the interpolation. This gives the emulator high versatility, in that it will be able to make predictions for a broad range of non-standard cosmologies absent from the training set.

3.1 Emulation strategy

An in-depth discussion of our emulation strategy, employing Gaussian process (GP) regression, can be found in Appendix A. The accuracy of this method is sensitive to the particular observable being emulated. Our aim is to predict the pseudo non-linear matter power spectrum discussed in Section 2.1, a quantity which has a dynamic range of many orders of magnitude and features a change in gradient sign between small and large scales, as shown in the upper panel of Figure 2 for the extensions to $\Lambda$CDM listed in Table 1. GP emulators on the other hand perform better at predicting smooth monotonic functions with a narrow dynamic range. Hence, following Euclid Collaboration et al. (2018), we consider the natural logarithm of the boost factor, $B(k,z)$, as our observable of interest, where

[TABLE]

The lower panel of Figure 2 shows that the natural logarithm of the boost factor indeed matches the desired properties of a GP observable.

An essential step in the emulator design is data compression. For any particular cosmological model, the power spectrum is sampled both in spatial frequencies (wavenumbers) and time (redshifts), for a total of $\sim 10^{4}$ support points. Training the emulator on each of these individual measurements requires prohibitively large computational resources. To overcome this issue, we perform a principal component analysis (PCA) following previous work (see, e.g, Habib et al., 2007; Schneider et al., 2008; Heitmann et al., 2016). For a given parameter vector, $\boldsymbol{\pi}$, and redshift, $z$, we apply the following linear decomposition

[TABLE]

where $\mu$ and $\{\psi_{1},\cdots,\psi_{n_{\psi}}\}$ are, respectively, the mean and basis functions of the training set of boosts at $z$, and the weights $w_{j}$ are the projections of the mean-subtracted observable onto $\psi_{i}$. The error term $\epsilon$ has two contributions: one component, arising from the negligible noise on the halofit training predictions (given by $\epsilon_{\rm n}(\boldsymbol{\pi})$ in equation 17); and a second component coming from missing relevant basis functions. We find that $n_{\psi}=12$ is sufficient to satisfactorily reconstruct the training set, and including more basis functions has only a marginal effect on the emulator accuracy. Ultimately, for a given $\boldsymbol{\pi}^{\ast}$ our emulator infers the boost factor by predicting the weights, $w_{j}(\boldsymbol{\pi}^{\ast})$, conditioned on the training set $w_{j}(\boldsymbol{\pi})$.

3.2 Model-independent parameterisation of the pseudo cosmologies

Extensions to $\Lambda$CDM are described by an arbitrary number of additional physical parameters. Consider, for example, $w$CDM, $f(R)$ gravity and massive neutrino cosmologies: matter clustering in all these models depends on five standard cosmological parameters, e.g. $\boldsymbol{\pi}_{\Lambda}=\{\Omega_{\rm{b}}h^{2},\Omega_{\rm{m}}h^{2},h,n_{s},\ln(10^{10}A_{s})\}$, along with one or more parameters describing model-specific features, $\boldsymbol{\pi}_{\rm X}=\{f_{R0},\sum m_{\nu},\cdots\}$. To construct an emulator capable of predicting generic pseudo non-linear matter power spectra, we opt for a model-independent parameterisation mapping an arbitrary model to a certain coordinate, $\boldsymbol{\pi^{*}}=\{\boldsymbol{\pi}_{\Lambda},\Delta\boldsymbol{\alpha}\}$, in our emulation parameter space, where $\Delta\boldsymbol{\alpha}$ is a vector of effective parameters that quantifies smooth deviations from the linear matter power spectrum of the standard cosmology. Thus, for $\Delta\boldsymbol{\alpha}=0$ we recover the late-time physics of the vanilla $\Lambda$CDM model. For this, we model the ‘shape’, namely the ratio of the linear pseudo matter power spectrum to that of the $\Lambda$CDM cosmology sharing the same standard parameters, as

[TABLE]

where $\{\Phi_{1},\cdots,\Phi_{n_{\Phi}}\}$ are the principal components of a training set, $\mathbb{T}$, consisting of smooth curves capturing a range of ‘reasonable’ deviations from $\Lambda$CDM 333Eq. (9) can be understood in terms of a principal component reconstruction based on the training set $\mathbb{T}$, that is

$S(k_{m})\approx\nu(k_{m})+\sum_{i=1}^{n_{\Phi}}\Phi_{i}(k_{m})\alpha_{i}\,,$

(10)

where for simplicity we dropped the dependence on redshift and cosmology, and made explicit that we sample the shape $S$ at $k_{m}\in\{k_{1},\cdots,k_{n_{k}}\}$. The function $\nu$ is the mean across $\mathbb{T}$, and in general we expect $\nu\neq 1$. Then, for any $\Lambda$CDM cosmology we must have

$1-\nu(k_{m})\approx\sum_{i=1}^{n_{\Phi}}\Phi_{i}(k_{m})\alpha_{i}^{\Lambda}\,,$

(11)

with $\alpha_{i}^{\Lambda}$ being the weights that exactly compensate the departure of the mean from unity, and upon using the orthonormality of the basis functions they read

$\alpha_{i}^{\Lambda}=\sum_{m=1}^{n_{k}}\Phi_{i}(k_{m})[1-\nu(k_{m})]\,.$

(12)

Finally, the combination of Eqs. (10) and (11) gives Eq. (9), with $\Delta\alpha_{i}\equiv\alpha_{i}-\alpha_{i}^{\Lambda}$. Note that for a given training set, $\mathbb{T}$, the $\Lambda$CDM weights need to be computed only once, and can subsequently be hard-coded in the emulator.. Use of PCA to parameterise deviations from vanilla $\Lambda$CDM cosmology has been previously implemented in a number of studies (see, e.g., Zhao et al., 2009; Hojjati et al., 2012; Eifler et al., 2015).

3.2.1 Generation of the basis set

We assemble a set of orthogonal basis functions to reconstruct a broad class of scale-dependent linear departures from the $\Lambda$CDM power spectrum. For this purpose, we generate random curves in Fourier space over the range $10^{-4}<k\,{\rm Mpc}/h<10$, with the constraint that all curves must converge to unity on large enough scales. This last requirement is motivated by the physics of well-known $\Lambda$CDM extensions, such as Generalised Brans-Dicke theories (De Felice & Tsujikawa, 2010; Park et al., 2010; Hinterbichler et al., 2011; Pogosian & Silvestri, 2016; Quiros et al., 2016; Joyce et al., 2016) and massive neutrino cosmologies (Lesgourgues & Pastor, 2006). Both these models possess a characteristic scale, $\bar{k}$, such that on scales $k\ll\bar{k}$ the linear matter clustering matches that of the vanilla $\Lambda$CDM cosmology. For $k\gg\bar{k}$ the growth of structure is either suppressed ($S<1$) or enhanced ($S>1$). Note that models with $\bar{k}\rightarrow 0$ exhibit in practice a scale-independent linear growth (see, e.g., Pogosian & Silvestri, 2016), and a simple rescaling of the amplitude of the vanilla $\Lambda$CDM power spectrum is sufficient to produce the corresponding pseudo power spectrum. Viable models in these two classes of theories above have $\bar{k}\gtrsim 10^{-3}\,h/\rm{Mpc}$ (Brax et al., 2012; Wang et al., 2012; Joyce et al., 2015; Planck Collaboration et al., 2018). We therefore adopt the value $\tilde{k}=10^{-3}\,h/\rm{Mpc}$ for the scale at which the shapes converge to unity, and allow for generous positive smooth deviations in the range $k\in[\tilde{k},10]\,h/{\rm Mpc}$. Moreover, since the basis functions derived from our generated shapes are sufficiently general to describe the power spectrum ratios at various redshifts, we drop the $z$-dependence, $\Phi_{i}(k,z)\rightarrow\Phi_{i}(k)$, and use one basis set for all redshifts.

In order to reduce the sensitivity of the emulator to the particular method used to generate the set of random shapes, we employ two different and independent generating strategies. The first of these methods, which we refer to as GPCurves, draws shapes from a 1-dimensional Gaussian process with a squared exponential kernel of the form given in equation 18 (where $d=1$), conditioned on a ‘training set’ of densely sampled points at $k<10^{-3}\,h/\rm{Mpc}$, all with values equal to unity, helping to fix the random shapes to this value on large scales. Here, for the kernel amplitude, we use $A=5$, which corresponds to the maximum value that the random shapes can have on scales $k>10^{-3}\,h/\rm{Mpc}$. We tune the only length-scale parameter, $p$, so that the artificial shapes exhibit at most one stationary point for $k>10^{-3}\,h/\rm{Mpc}$, thus resembling features of the physical shapes in Figure 1. There is no intrinsic restriction in the GP preventing generated shapes from taking unphysical negative values. We therefore generate an initial set of curves from which we sample 1000 viable positive shapes.

As a second method we use Smurves (Moews et al., 2019), a novel random curve generator. The functionality of the generator is illustrated by considering a particle travelling in Fourier space from low-to-high wavenumbers. Before $k=10^{-3}\,h/\rm{Mpc}$, the particle experiences no forces and hence travels undisturbed at a vertical coordinate of unity. After this point, forces randomly generated from a uniform distribution perturb the particle in the vertical direction, altering its trajectory. Also randomly generated are the locations of gradient sign changes. In this manner, we produce an additional 1000 smooth curves with a maximum of one stationary point and vertical intervals in the range [0, 8]. In order to facilitate extra flexibility in the deviations from $\Lambda$CDM, the parameters specifying the curve generation in Smurves are set purposefully to create a sample of curves which are slightly broader than GPCurves.

A comparison of 10 curves randomly selected from either samples are shown in different colours in the upper panel of Figure 3. The upper-middle panel shows the shapes corresponding to the physical models in Table 1 (orange) relative to the collective sample of 2000 random curves from GPCurves and Smurves (grey). The lower-middle panel shows the orthogonal basis functions, $\Phi_{i}$, obtained via a PCA of the 2000 curves. Note that the resulting basis set allows for a broader range of departures from $\Lambda$CDM than our test physical models.

Each basis function, $\Phi_{i}$, appearing in the shape reconstruction (equation 9) contributes one extra dimension to the parameter space of the emulator, with the additional parameters corresponding to the PCA weights, $\Delta\alpha_{i}$. Therefore, in order for our emulator to be computationally efficient whilst still achieving accuracies of $\lesssim$1%, we keep the smallest number of basis functions that guarantees sub-per-cent reconstruction errors on all physical test shapes (see Figure 1). We find that $n_{\Phi}=8$ basis functions are sufficient to reconstruct $f(R)$ gravity shapes with negligible errors. For cosmologies with massive neutrinos, however, changes to the BAO dynamics at early times produce rapid oscillations in their late-time power spectrum shapes. These features require one additional step in the reconstruction, but no extra parameters need be included in the emulation. Further details on this can be found in Appendix B. The resultant shape reconstruction accuracy for all physical models is shown in the lower panel of Figure 3. Overall our emulation volume is therefore carved out of a 13-dimensional space: 5 $\Lambda$CDM parameters, $\boldsymbol{\pi}_{\Lambda}$, and 8 shape parameters, $\Delta\boldsymbol{\alpha}$.

3.3 Building the training and trial sets

Having established a model-independent description for arbitrary cosmological models, we now discuss how we generate training and trial sets for our pseudo matter power spectrum emulator. In Section 3.3.1 we give details on the Latin hypercube (LH) space-filling strategies we adopt to efficiently sample our 13-dimensional space. In Section 3.3.2, we then define suitable ranges for the cosmological and shape parameters, $\{\boldsymbol{\pi}_{\Lambda},\Delta\boldsymbol{\alpha}\}$, such that the emulated boost factors are representative of viable extensions to $\Lambda$CDM. The nodes of the LHs take values in the range $[0,1]$, which we then properly rescale to obtain the “physical” coordinates $\boldsymbol{\pi}$ used to generate the training and trial sets.

3.3.1 Latin hypercube sampling

The use of LH-based distributions of points to evenly sample a parameter space was first introduced by McKay et al. (1979). The LH is a generalisation of the concept of a Latin square in which each row and column of a two-dimensional grid features exactly one sample, taking the form of a chess board with a number of rooks that do not threaten each other. Owing to the desirable space-filling properties, LH sampling of the input parameters for simulations has become a standard practice (Sacks et al., 1989; Currin et al., 1991; Heitmann et al., 2006; Schneider et al., 2008; Agarwal et al., 2012; Liu et al., 2015; Garrison et al., 2018; Harnois-Deraps et al., 2019).

Assembling 13-dimensional LH designs therefore presents an ideal way to construct predictions to both train the emulator and test its accuracy for different regions of the emulation space. The spatial distribution of the trial cosmologies does not require optimisation. Therefore we sample 300 trial coordinates from a 13-dimensional LH generated with the pyDOE Python package444https://pythonhosted.org/pyDOE/. The training nodes, however, impact significantly on the emulator accuracy, and sampling methods able to cover the high-dimensional volume as uniformly as possible are thus preferable. For this purpose, space-filling criteria based on the distance between samples are commonly used. One of the most prominent approaches is the Maximin LH design, which maximises the minimum distance between samples (Morris & Mitchell, 1995). While this approach leads to the sought-after space-filling properties in the bulk of the parameter space, most points cluster in the corners or at the boundaries of the lower-dimensional projection spaces. Since the matter power spectrum is most sensitive to subsets of cosmological and shape parameters, space-filling strategies also optimising projection space sampling would give an obvious advantage.

To this end, Joseph et al. (2015) introduced the maximum projection design (MaxPro), in which a weighted Euclidean distance, $E$, is minimised, resulting in the new criterion

[TABLE]

where $\boldsymbol{D}=\{\boldsymbol{x}_{1},\dots\boldsymbol{x}_{N}\}$ is an $N\times d$ design matrix containing $N$ samples occupying the $d$-dimensional parameter space. $x_{il}$ denotes the value of the $i$-th sample’s coordinate in the $l$-th dimension. The MaxPro criterion maximises the bulk and projection space-filling properties, since for any dimension $l$, if $x_{il}=x_{jl}$ and $i\neq j$, $E(\boldsymbol{D})=\infty$. This ensures that in the design minimising $E(\boldsymbol{D})$ no two points can be close to each other in any of the projections.

The MaxPro design optimisation implemented in this work is based on the simulated annealing algorithm described in Morris & Mitchell (1995)555https://CRAN.R-project.org/package=MaxPro. This consists of selecting an initial configuration before searching for progressively better design choices by randomly perturbing the current design, keeping the new matrix should it reduce the cost function, $E(\boldsymbol{D})$, or else do so with a given probability. Finally this implementation arrives at a locally optimal MaxPro design by ascending the gradient given by

[TABLE]

To investigate the performance of our emulator we generate sets of training nodes with sizes $N=100,300,500$ and 600, with the additional property of having an optimal space-filling pattern in each projection. Furthermore, to investigate the sensitivity of our results to the specific configuration of nodes in each training set, we produce 10 optimised designs of each size and measure the range of achieved accuracies across them.

3.3.2 Setting parameter ranges

We shall now map our dimensionless training and trial nodes belonging to the unit LHs onto our parameters, $\boldsymbol{\pi}\in\{\boldsymbol{\pi}_{\Lambda},\Delta\boldsymbol{\alpha}\}$. To do this, we must identify appropriate ranges for these parameters. We begin with the $\Lambda$CDM parameters, and note that our emulator is intended for cosmological analyses employing high-quality data from the next-generation galaxy surveys. Hence, we set the boundaries on $\boldsymbol{\pi}_{\Lambda}$ to be roughly consistent with the 95% marginalised constraints from the combined TT, TE, EE $+$ lowE analysis in Planck Collaboration et al. (2018), namely,

[TABLE]

These ranges are illustrated in Figure 4, which shows the 2D projections of the $\Lambda$CDM sub-space, with the grey points corresponding to one of the $N=500$ MaxPro distributed training designs.

Turning to the shape parameters, $\Delta\boldsymbol{\alpha}$, we take inspiration from the physical models listed in Table 1, and simply set the range of values in our training set to be slightly broader than those defined by the $\Delta\boldsymbol{\alpha}$’s associated with $f(R)$ gravity, massive neutrino and combined cosmologies. The upper panel of Figure 5 shows the $\Delta\alpha_{i}$ for the physical models at $z=0$ in orange, and those related to one of the designs with 500 optimally spaced training nodes (see Section 3.3.1) in grey. The convergence of the $\Delta\alpha_{i}$ values towards zero as the PCA rank, $i$, increases, is an indication that we are using a sufficient number of basis functions to reconstruct our test models.

The lower panel of Figure 5 shows the shapes corresponding to the $\Delta\boldsymbol{\alpha}$ vectors shown in the upper panel, both for the 500-node training set (in grey) and for the physical models (in orange), reproduced from Figure 1. We see that the training shapes cover a complex range of behaviours, whilst still encompassing the shapes of the physical models.

For the trial set, we opt for narrower parameter ranges, such that in each dimension the trial coordinates have values in the central 75% of the training set intervals. In doing this we test the emulator accuracy only in the bulk of the parameter space, so that our conclusions are less influenced by edge effects, where emulator accuracy degrades due to a shortage of nodes near the boundaries. This is visualised in Figure 4 by the distribution of trial nodes (in magenta). We also test the performance of our emulator for strictly $\Lambda$CDM cosmologies, where all shape parameters are set to zero. This $\Lambda$CDM trial ensemble is used to test whether the emulator is able to make accurate predictions for the standard model in the non-linear regime, even when $\Lambda$CDM cosmologies are completely absent from the training set.

3.3.3 Generating boost factors

With the training and trial nodes now assembled, the final task is to produce the corresponding linear as well as non-linear matter power spectra, and hence the logarithm of the boost factor, $B(k,z)$ (see equation 7), which the emulator is trained on and predicts. Conventionally, the non-linear statistics are measured from cosmological simulations, but for the purposes of this proof-of-concept study halofit can serve as a proxy for a suite of $N$-body simulations. Note that this approach is self-consistent as long as the growth of structure and background evolution of the ‘simulated’ cosmologies matches those of the $\Lambda$CDM model, which is indeed the case for the pseudo cosmologies. To estimate the surrogate $P^{\rm{pseudo}}_{\rm{NL}}$ we repeat the process described in Section 2.1 for the physical models: for a given pseudo cosmology, $\boldsymbol{\pi}=\{\boldsymbol{\pi}_{\Lambda},\Delta\boldsymbol{\alpha}\}$, our modified version of camb accepts as input the random shape described by $\Delta\boldsymbol{\alpha}$, and internally multiplies this by the $\Lambda$CDM linear power spectrum associated with $\boldsymbol{\pi}_{\Lambda}$. This generates the pseudo linear power spectrum which together with $\boldsymbol{\pi}_{\Lambda}$ serves as input for halofit non-linear predictions. The boost factor is then produced by simply taking the ratio of the non-linear to the linear power spectrum. We generate boost factors for the training and trial sets of our emulator at redshifts 0 and 1. The $z=0$ pseudo linear and non-linear matter power spectra, along with the corresponding boost factors, are shown in Figure 6 for a 500-node training set (in grey), as well as for the physical models (in orange) also illustrated in Figure 2.

4 Results

The non-linear matter power spectrum is key to the derivation of a number of cosmological observables, therefore we use this rather than the emulated statistic, $\ln B(k,z)$, to evaluate the accuracy of our emulator. We consider wavenumbers in the range $k\in[0.01,10]\,h/\rm{Mpc}$, excluding larger scales where the boost factor goes to unity in all models, presenting no challenge to the emulator. For each trial cosmology and redshift separately, we quantify the emulator accuracy as

[TABLE]

where $P^{\rm emu}_{\rm NL}$ and $P^{\rm True}_{\rm NL}$ refer to the emulated and true pseudo non-linear power spectra respectively, and for the purposes of this study we take halofit non-linear predictions as the truth. Figure 7 presents the fraction of the trial cosmologies with emulation errors smaller than a threshold value, $\bar{\epsilon}$, for various sizes of the training sets at both $z=0$ (solid lines) and $z=1$ (dashed). The fraction shown is the average across the 10 realisations per training set size. The left panel shows the results for the pseudo cosmology trials, whilst the right panel corresponds to the pure $\Lambda$CDM trial ensemble.

At $z=0$, for the pseudo cosmology trials and training sets of size $N=500$ or 600, the emulation is accurate to better than 1% for about $85\%$ of the trials, and the full ensemble is recovered to better than 2% accuracy. As expected, the emulator performs noticeably worse when trained on a smaller number of nodes, with only 13% (60%) of the pseudo non-linear spectra emulated to better than 1% when $N=100$ (300). This is because the emulator has less information on the complex relationship between input parameters, $\boldsymbol{\pi}$, and the boost factors. However, we find that the largest inaccuracy across all trial models is still $<6\%$ ($<3\%$) for $N=100$ (300).

At $z=1$, the fraction of pseudo trials with $\epsilon^{\rm emu}<1\%$ are similar to those at $z=0$. However, the reduced steepness of the cumulative distributions at $z=1$, compared to those at $z=0$, suggests an increasing fraction of outliers with redshift. This is caused by the training and trial boost factors displaying a broader dynamic range at the higher redshift, which poses a greater challenge to emulation. Such broadening can be understood in terms of differences in the boost factor evolution induced by the different growth rates and non-linear mode couplings for the various base $\Lambda$CDM parameters and shapes.

The right panel of Figure 7 shows similar trends for the $\Lambda$CDM trial ensemble. For the two largest training sets we find that over 90% of the trial cosmologies are reproduced with better than percent accuracy at both redshifts analysed. The larger emulation space dimensionality compared to similar emulator schemes designed for $w$CDM cosmologies only (see, e.g. Heitmann et al., 2014) demands a three-fold increase in the number of nodes to achieve equivalent accuracies for $\Lambda$CDM cosmologies. In fact, for our smallest training set ($N=100$) the emulator predictions can have errors as large as 4% (7%) at $z=0$ ($z=1$), comparable to those obtained by Heitmann et al. (2014). Although our emulation strategy requires a few hundred training nodes to reach percent accuracy, it is important to realise that it provides non-linear power spectrum predictions for a much broader class of cosmologies, including the $\Lambda$CDM model. Moreover, with the aid of pairing and fixing techniques (Angulo & Pontzen, 2016; Euclid Collaboration et al., 2018) the total number of training simulations required is analogous to that of other GP emulators (Lawrence et al., 2010, 2017).

Finally we assess the accuracy of the emulator for the physical models in Table 1. Motivated by our results for the synthetic trial cosmologies, for this second test we use $N=500$ training nodes, since a larger training set yields no significant improvement in the emulator performance.

The full range of the fractional errors of the emulated $P^{\rm{pseudo}}_{\rm{NL}}$ relative to the halofit predictions, across the 10 training set realisations and all physical models, are shown in Figure 8 for $z=0$ (upper panel) and $z=1$ (lower panel). The outer (pale) and inner (dark) bands show the ranges including or excluding the xMNU_0.4 cosmology, respectively. Comparison of these reveals that much of the observed error is associated with this extreme model, which resides outside of the central 75%-per-side hypervolume used in evaluating the accuracies with the synthetic pseudo cosmologies. We find that the remaining models are recovered to better than the target percent accuracy for $k<10\,h/\rm{Mpc}$ at $z=0$, with the majority of the training set realisations. We note that the range plotted, showing the total variation obtained across the training set realisations, presents a conservative estimate of the emulator performance; the $2\sigma$ range of absolute emulation accuracies, for example, is considerably tighter. We observe a slightly higher outlier rate at $z=1$, consistent with our observations with the synthetic pseudo cosmologies. Of the models included in the dark band, we find, once again, that the largest inaccuracies occur for cosmologies near the boundary of the sub-volume used to test the synthetic shapes, MNU_0.4 and xMNU_0.2, both lying many standard deviations away from the Planck best fit cosmology (Planck Collaboration et al., 2018).

5 Discussion and Conclusions

In this paper, we have designed and tested a methodology essential for predicting accurate non-linear matter power spectra in a broad range of extensions to the $\Lambda$CDM cosmology in the context of the halo model reaction framework (Cataneo et al., 2019, ‘C19’). The effectiveness of the C19 strategy rests on the availability of accurate baseline $\Lambda$CDM non-linear power spectra evolved from non-standard initial conditions, the pseudo power spectra (see Section 2.1). We showed that such quantities can be readily predicted with a Gaussian process emulator whilst simultaneously satisfying the stringent accuracy requirements. We demonstrated the power of this technique for theoretically motivated models, such as $f(R)$ gravity and massive neutrino cosmologies, as well as for phenomenological models, where in both cases structures on different scales grow at different rates.

Based on the results of this work, we advocate the following implementation of our methodology in the future:

A 13-dimensional optimised Latin hypercube is defined in terms of: $\boldsymbol{\pi}_{\Lambda}\,(\rm{dim}=5)$, the cosmological parameters that define the linear $\Lambda$CDM power spectrum, and $\Delta\boldsymbol{\alpha}\,(\rm{dim}=8)$, the shape parameters which modify the linear $\Lambda$CDM power spectrum to create the linear power spectrum for a given beyond-$\Lambda$CDM model. 2. 2.

Two standard $\Lambda$CDM $N$-body simulations, with initial conditions paired and fixed following Angulo & Pontzen (2016), are run for each node in the LH, i.e. with the input cosmology $\boldsymbol{\pi}_{\Lambda}$, and rescaling of the primordial power spectrum with the shape defined by $\Delta\boldsymbol{\alpha}$. 3. 3.

From the output of each $N$-body simulation, the pseudo non-linear matter power spectrum, $P^{\rm pseudo}_{\rm NL}(k,z)$, and hence the natural logarithm of the boost factor, $\ln B(k,z)$ (see equation 7), are measured. An emulator is subsequently trained to predict this quantity for any $\{\boldsymbol{\pi}_{\Lambda},\Delta\boldsymbol{\alpha}\}$ combination. 4. 4.

A user will provide the analytically computed linear power spectrum in their chosen beyond-$\Lambda$CDM model, $P_{\rm L}^{\rm real}(k,z;\boldsymbol{\pi}_{\Lambda},\boldsymbol{\pi}_{X})$, where $\boldsymbol{\pi}_{X}$ is a vector of parameters describing the particular extension to the standard cosmology, and the corresponding vector of base cosmological parameters, $\boldsymbol{\pi}_{\Lambda}$. We then perform a PCA decomposition of the ratio, $P_{\rm L}^{\rm real}(k,z;\boldsymbol{\pi}_{\Lambda},\boldsymbol{\pi}_{X})/P_{\rm L}^{\Lambda}(k,z;\boldsymbol{\pi}_{\Lambda})$, where the denominator is the standard $\Lambda$CDM linear spectrum, to map $\boldsymbol{\pi}_{X}$ to the eight PCA weights (see equation 9). 5. 5.

Given the set of values of the base cosmological parameters and PCA weights for the queried model, the emulator predicts $\ln B(k,z)$. 6. 6.

The pseudo non-linear matter power spectrum is then simply obtained as $P_{\rm NL}^{\rm pseudo}(k,z)=B(k,z)\times P_{\rm L}^{\rm real}(k,z)$. 7. 7.

Finally, a rescaling of the pseudo non-linear power spectrum by the C19 halo model reaction, in equation 4, produces the real non-linear power spectrum including the effects of physics for the chosen beyond-$\Lambda$CDM model.

The advantage of our emulator scheme over existing approaches (Heitmann et al., 2014; Lawrence et al., 2017; Euclid Collaboration et al., 2018; Winther et al., 2019) is twofold: the model-independent parameterisation enables predictions for beyond-one-parameter extensions to $\Lambda$CDM, effectively expanding the domain of applicability to a much broader class of cosmologically interesting models; and the simulations required for the training phase are all based on the standard $\Lambda$CDM cosmology, which makes them easier and faster to run.

This work is a first step towards a proper and accurate non-linear matter power spectrum emulator, one with a training set built from the output of $N$-body simulations. Here, instead, our goal was to perform a feasibility study, in which the approximate semi-analytical predictions of halofit act as surrogate simulations (replacing step (ii) above). Taking inspiration from the linear power spectrum shape of well-studied extensions to the standard cosmology, we computed the training pseudo non-linear power spectra by feeding smooth, random modifications of the $\Lambda$CDM linear power spectrum to halofit. Despite its intrinsic 5-10% inaccuracies, previous studies showed that halofit provides a quick and robust way to design, construct and test the emulator (Heitmann et al., 2009; Euclid Collaboration et al., 2018).

With the aid of pairing and fixing techniques for the initial conditions (Angulo & Pontzen, 2016), we estimate that a total of 1000 $\Lambda$CDM-like simulations (two realisations for each of the 500 training nodes) is enough to build a per-cent-level accurate pseudo non-linear matter power spectrum emulator. This setup also allows for predictions of the $\Lambda$CDM non-linear matter power spectrum with $\lesssim 1\%$ accuracy, even with the standard cosmology being completely absent from the training set. Despite the remarkable amount of computational resources required, the number of simulations is not too dissimilar from that used for the training of previous, more limited emulators (Lawrence et al., 2010, 2017).

Although not shown here, our methodology can also potentially be extended to include the effects on baryonic physics on non-linear scales (see, e.g., Schneider et al., 2019; Debackere et al., 2019). This correction, which could be incorporated into the halo model reaction, is left for future investigation.

The results presented here should be regarded as conservative, as other emulation schemes, such as sparse polynomial chaos expansion (Blatman & Sudret, 2011; Euclid Collaboration et al., 2018), and further design optimisation (Rogers et al., 2019; Caron et al., 2019) could help overcome the shortcomings of Gaussian process regression. Furthermore, we have only considered departures from $\Lambda$CDM on scales of $k>10^{-3}\,h/\rm{Mpc}$. This requirement could be relaxed by incorporating only one additional emulation parameter, $\Delta\bar{k}$, describing translations of the scale at which departures occur about this point. Since the scope of this preliminary study was to show the feasibility of pseudo cosmology emulation, these topics are also left to future work.

Our emulation method, coupled to the halo model reactions, enables future cosmic shear analyses by providing fast, accurate, and flexible predictions of the matter power spectrum deep in the non-linear regime, where we are more likely to find imprints of new physics beyond the concordance cosmology, if any (see, e.g., Heymans & Zhao, 2018). Moreover, the approach developed in this work can be extended to predict other pseudo cosmology statistics, for instance galaxy clustering (Zhai et al., 2019), as well as mean halo properties, such as their abundance and concentration (Kwan et al., 2013; McClintock et al., 2019).

6 Acknowledgements

BG, MC and CH acknowledge support from the European Research Council under grant number 647112. BM acknowledges the support of a Principal’s Career Development Scholarship from the University of Edinburgh. MC thanks Joey Faulkner for stimulating discussions in the early stages of this work.

Appendix A Gaussian process regression method

The mathematics behind Gaussian process (GP) regression emulators have been covered extensively in previous work; we refer the interested reader to Rasmussen & Williams (2006) for a general discussion of GP, and to Habib et al. (2007) and Schneider et al. (2008) for its first applications to cosmology. Here we summarise its key points, which we implement using the publicly available code scikit-learn666https://scikit-learn.org (Pedregosa et al., 2011).

GP regression is a non-parametric Bayesian machine learning algorithm for constraining the distribution of functions consistent with observed data. Let us consider a data set, $\mathcal{D}$, consisting of $n$ realisations of an observable, $\boldsymbol{y}$, each corresponding to different values of the $d$-dimensional input parameters, $\boldsymbol{\pi}\in\mathbb{R}^{d}$, i.e. $\mathcal{D}=\{(\boldsymbol{\pi}_{i},y_{i})|i=1,...,n\}$. The ensemble of matter power spectra measured from a suite of $N$-body simulations is a typical example of such a data set, and constitutes a training set for the emulator. The task of the GP emulator is to learn the distribution of functions, $f(\boldsymbol{\pi})$, which are consistent with the mapping between the training set input parameters - the ‘nodes’ - and output, via

[TABLE]

where $\epsilon_{\rm n}(\boldsymbol{\pi})$ is a noise term sampled from a mean-zero Gaussian distribution with a standard deviation given by the error on the training set observable, $y(\boldsymbol{\pi})$. The prediction, $y^{*}$, corresponding to an arbitrary coordinate $\boldsymbol{\pi^{*}}$, is then sampled from a generalisation of a Gaussian posterior probability distribution over the range of consistent functions777In this work we consider halofit training predictions as a proxy for hypothetical numerical simulations with negligible error. Nevertheless, to prevent unwanted numerical instabilities we use the arbitrarily small constant default $\epsilon_{\rm n}$ in scikit-learn for all $y(\boldsymbol{\pi})$.. In other words, the GP emulator interpolates the observable from the input coordinates of the training set to the trial coordinates across a $d$-dimensional parameter space.

A key ingredient of our posterior is the Gaussian prior distribution of functions deemed to reasonably map between input and output. The prior is determined by a mean, conventionally taken to be zero, and a covariance function, known as the ‘kernel’. The kernel can take various functional forms, each described by a vector of hyperparameters, $\boldsymbol{h}=\{A,p_{1},\cdots,p_{d}\}$, governing the kernel’s behaviour. Following Heitmann et al. (2009), in this work we adopt the squared exponential form, which specifies the covariance between the functions $f(\boldsymbol{\pi})$ and $f(\boldsymbol{\pi^{*}})$ as

[TABLE]

This kernel has the following properties: (1) the covariance varies smoothly within the parameter space; (2) it depends only on the Euclidean distance between points, such that $K(f,f^{*};\boldsymbol{h})=K(f^{*},f;\boldsymbol{h})$; (3) predictions become maximally correlated when $\boldsymbol{\pi}=\boldsymbol{\pi^{*}}$; (4) the correlation is large for points in relative proximity and small for largely separated points; (5) each $p_{l}$ corresponds to the functions’ characteristic length-scale of variation in each of the $d$ dimensions, whilst $A$ is the kernel amplitude.

The emulator is trained by finding values for the hyperparameters such that a distribution of functions which are optimally consistent with all realisations in the training set can be defined. In this work, we fit for these using the method built-in to scikit-learn, which employs a gradient ascent optimisation of the marginal likelihood conditioned on the training set.

Appendix B Modelling BAO residuals

As shown in Figure 1, models with massive neutrinos have shapes characterised by a rapidly oscillating component superimposed over a smooth function. These wiggles are associated with changes in the baryon acoustic oscillations (BAO) produced by differences in the sound horizon at the end of the drag epoch and shallower gravitational potentials before recombination (Hu & Sugiyama, 1996; Hu & White, 1996; Eisenstein & Hu, 1998; Lattanzi & Gerbino, 2017) compared to a massless neutrino $\Lambda$CDM cosmology with the same standard parameters, i.e. $\boldsymbol{\pi}_{\Lambda}^{m_{\nu}=0}=\boldsymbol{\pi}_{\Lambda}^{m_{\nu}\neq 0}$, and $\Omega_{\rm c}^{m_{\nu}=0}=\Omega_{\rm c}^{m_{\nu}\neq 0}+\Omega_{\nu}^{m_{\nu}\neq 0}$. Depending on the value of the standard cosmological parameters and neutrino mass, the amplitude of the BAO residual in the shape can be up to $\sim$1% of the slowly varying component. The PCA reconstruction outlined in Section 3.1 requires a large number of components to simultaneously capture long- and short-scale variations. Hence we find errors $\gtrsim 1\%$ in reconstructing the $\sum m_{\nu}$-induced oscillations with only $n_{\Phi}=8$ basis functions, as shown by the grey curves in the upper panel of Figure 9.

In order to improve the reconstruction accuracy of the PCA, we isolate and remove the BAO residuals by applying the de-wiggling algorithm presented in Appendix C. In short, the BAO in Fourier space maps to a localised bump in the discrete sine transform of the power spectrum. By removing this bump and performing the inverse transform, we obtain a smooth, “de-wiggled” version of the power spectrum. This is performed on both the model featuring massive neutrinos and on the $\Lambda$CDM baseline, before calculating the de-wiggled shape given by

[TABLE]

where the superscript ‘$\rm{dw}$’ on the quantities on the right-hand side denote power spectra obtained from applying the de-wiggling algorithm. This smoothed component of the shape is captured very accurately by the basis functions. Reconstruction of the full shape then rests on modelling the remaining oscillatory component, given by

[TABLE]

Notice that all quantities in equations 19 and 20 depend implicitly on the cosmological parameters. We experimented with fitting these oscillations with the function

[TABLE]

where $a$ and $b$ are free parameters controlling the amplitude and modulation of a fiducial wiggle template, $\mathcal{W}_{\rm fid}$, which we identify with the wiggle contribution of the MNU_0.4 cosmology in Table 1, represented by the red curve in the lower panel of Figure 9. We checked that our fiducial template choice in equation 21 is robust against changes in redshift, standard cosmological parameters, sum of neutrino masses and neutrino mass eigenstates. This ensures that the parameters $\{a,b\}$ are sufficient to capture BAO variations in the family of pseudo cosmologies treated in this work. The blue curves in the lower panel of Figure 9 show the behaviour of our fitting formula, Eq. (21), for all of the test cosmologies in Table 1 containing massive neutrinos.

In summary, for each redshift and cosmology separately, we reconstructed the total shapes of the physical models as

[TABLE]

where only the smooth, de-wiggled term is modelled via the PCA reconstruction given in equation 9. The magenta curves in the upper panel of Figure 9 show the overall accuracy of the two-step reconstruction, originally presented in the lower panel of Figure 3. We find that the inaccuracies in the BAO regime are greatly suppressed compared to performing the PCA directly on the full shape (grey curves).

The fact that the cosmologies with massive neutrinos require two extra parameters to summarise their shapes to sub-per-cent accuracy, might imply that the emulator need be informed of $\{a,b\}$ to achieve this target accuracy in its predictions, thereby increasing the dimensionality from 13 to 15. To test this, we produced MaxPro training sets optimised in 15 dimensions, with a range of wiggle parameters, $a$ and $b$, slightly broader than those obtained for the physical models with massive neutrinos. We then compared the emulator’s performance with the physical models, when trained on these 15-dimensional sets, and when trained on the same sets but in 13-dimensions, omitting $a$ and $b$. We found, in fact, that the emulator is mostly insensitive to these two parameters, which therefore contribute only marginally to improving its accuracy.

Figure 10 holds the key to this somewhat surprising result, where we adopt the xMNU_0.4 cosmology as a worst case scenario among our physical models. Plotted in solid orange is the $S/S^{\rm dw}$ ratio, which is equivalent to the ratio of pseudo linear power spectra with or without changes to the BAO scale compared to $\Lambda$CDM (see equation 19). In dashed dark blue is the ratio of their non-linear counterparts obtained with halofit. Solid magenta shows the ratio between the boost factors derived from using the full or de-wiggled shape, which is also equal to the ratio of the dark blue and orange curves. We see that although the BAO residual amounts to $\sim$1% of the non-linear power spectrum, its contribution to the boost factor is much smaller. In turn, this means that the emulator can reconstruct the boost factor with the required accuracy irrespective of whether the information on the BAO residual is provided or not. In other words, the boost factor effectively damps down the correction associated with changes to the BAO caused by massive neutrinos, eliminating the need for the additional wiggle modelling (equation 21)888Although we verified this statement with halofit we expect it to remain valid in simulations as well, assuming that halofit provides accurate predictions for the ratio $P_{\rm NL}/P_{\rm NL}^{\rm dw}$..

The fact that early-time modifications to the BAO physics have such a weak effect on the boost factor is a strength of our emulation scheme; an emulator designed to directly predict the pseudo non-linear power spectra is likely to incur $\gtrsim 1\%$ errors unless detailed modelling of the BAO residual is employed. The unimportance of the wiggle parameters means that the 13-dimensional emulator is sufficient to achieve our target accuracy. In Appendix D, we investigate the performance of the emulator with cosmologies featuring extra relativistic degrees of freedom, for which changes in the acoustic oscillations prior to the epoch of last scattering leave especially prominent wiggles in the late-time shape ratios. We find the emulator achieves sub-per-cent accuracy at $z=0$, but less consistent results at higher redshift for this model.

Although modelling of the BAO residuals proves to be unnecessary, their extraction through the de-wiggling algorithm remains important. By acting as a low-pass filter, the PCA reconstruction can in principle remove rapid oscillations and capture the broadband shape. However, the PC weights thus derived are somewhat different from those obtained after de-wiggling, generating an error that ultimately propagates to the boost factor. We compare the emulator performance with and without de-wiggling in Appendix E.

Appendix C De-wiggling methodology

Here we outline in more detail the method, discussed briefly in Appendix B, for de-wiggling linear matter power spectra following Baumann et al. (2018). This procedure is performed on the pseudo and $\Lambda$CDM spectra separately, and therefore we denote both with $P_{\rm{L}}(k)$ in the following:

We firstly interpolate the linear matter power spectra featuring BAO wiggles, $P_{\rm{L}}(k)$, from the logarithmically-spaced $k$ values onto a linear-spaced array, $k_{n}$, sampled at $2^{n}$ points, where $n$ is an integer. This is so as to improve the computational efficiency of the Fast Discrete Sine Transform (DST), which is then performed on $\log_{10}[k_{n}P_{\rm{L}}(k_{n})]$ using the orthonormalised type-II sine transform. The resulting array, which has indices denoted by $i_{\rm dst}$ and length $2^{n}$, is split into two separate arrays, one with the even $i_{\rm dst}$ indices and the other with odd $i_{\rm dst}$ indices. 2. 2.

The DST of the even and odd arrays each feature a bump, which is a localised manifestation of the BAO in the $k$-space of the matter power spectrum. This is shown by the solid lines for one of the physical models in the upper panel of Figure 11. In order to de-wiggle the power spectrum, we must remove these bumps. To this end, we take the second derivative of the DST, converting the observed bump in each array into a prominent single-period oscillation, as shown by the lower panel of Figure 11. We then identify $i_{\rm dst}^{\rm{min}}$ and $i_{\rm dst}^{\rm{max}}$ values bracketing the oscillation, which are close to the local minimum on the left-hand side of the approximate centre of this feature and the local maximum on the right-hand side respectively. 3. 3.

We remove the corresponding $i_{\rm dst}$-range from the even and odd DST arrays, multiply each by $(i_{\rm dst}+1)^{2}$ and interpolate across the gap between $i_{\rm dst}^{\rm{min}}$ and $i_{\rm dst}^{\rm{max}}$ with cubic splines. The $(i_{\rm dst}+1)^{2}$ factors are then divided out, generating smooth “de-bumped” even and odd DST arrays, where the information corresponding to the BAO has been removed. This is shown by the dashed lines in the upper panel of Figure 11. 4. 4.

We recombine the even and odd DST arrays into a single array. Performing an inverse type-II DST produces the logarithmic de-wiggled linear matter power spectra weighted by the linear-spaced $k_{n}$ array, $\log_{10}[k_{n}P^{\rm{dw}}_{\rm{L}}(k_{n})]$. From this we obtain the de-wiggled linear spectra itself, $P^{\rm{dw}}_{\rm{L}}(k_{n})$. 5. 5.

After performing this on both the pseudo and $\Lambda$CDM power spectra, we take the ratio of the respective $P^{\rm{dw}}_{\rm{L}}(k_{n})$ to obtain $S^{\rm{dw}}(k_{n})$, the de-wiggled shape sampled at the linearly-spaced wavenumbers. Finally we interpolate the de-wiggled shape back to the logarithmically-spaced $k$ sampling, obtaining $S^{\rm{dw}}(k)$. We find that there is less numerical noise at the upper and lower wavenumber bounds of the de-wiggled shape when we take the ratio of the linearly-sampled de-wiggled power spectra before interpolating, rather than the other way around.

Whereas Baumann et al. (2018) chose fixed values of the $i_{\rm dst}^{\rm{min}}$ and $i_{\rm dst}^{\rm{max}}$ relative to the oscillation in the second derivative of the DST arrays, we optimise for these parameters in our analysis with the physical models. Specifically, we cycle through a grid of $(i_{\rm dst}^{\rm{min}},i_{\rm dst}^{\rm{max}})$ combinations which are centred on the local minimum and maximum on either side of the oscillation, and identify those which return the smoothest de-wiggled shape. The smoothness is measured by integrating the second derivative of $S^{\rm{dw}}(k)$, convolved with a Gaussian filter, in the range of the BAO wiggles. The optimal $(i_{\rm dst}^{\rm{min}},i_{\rm dst}^{\rm{max}})$ combination is that which minimises this integral. We find the de-wiggled shapes are insensitive to the width of the Gaussian filter, which is employed so as to ensure numerical noise does not bias identification of the optimal $(i_{\rm dst}^{\rm{min}},i_{\rm dst}^{\rm{max}})$ range.

Appendix D Extension to extra relativistic degrees of freedom

Appendix B demonstrated the unimportance of modelling the BAO residuals in the shapes of the physical cosmologies, justifying the 13-dimensional emulation scheme. In order to further verify the validity of this, we test our emulator scheme against a theoretically motivated cosmology with an especially prominent BAO residuals in the shape ratio.

The existence of a relic sea of neutrinos with effective relativistic degrees of freedom, $N_{\rm eff}=3.046$, is a general prediction of the standard model (Mangano et al., 2005). However, physics beyond this paradigm includes scenarios with additional relativistic particles (or ‘dark radiation’) at the epoch of decoupling, typically quantified by the variation $\Delta N_{\rm eff}$, such that

[TABLE]

This has the effect of reducing the expansion rate and increasing the acoustic scale (Archidiacono et al., 2013), thereby acting to ease the tension between the early- and late-time measurements of the Hubble constant (Riess et al., 2019). We consider a model with $\Delta N_{\rm eff}=1$ extra relativistic degrees of freedom, which lies several standard deviations away from the best-fit result of the Planck Collaboration et al. (2018) analysis. The BAO residuals in the shape of this cosmology at $z=0$ are shown in the upper panel of Figure 12.

Following the same procedure employed for the physical models in Table 1, we obtain a smoothed version of this shape using the de-wiggling algorithm (see Appendix C), perform a PCA decomposition to obtain the PCA weights, $\Delta\boldsymbol{\alpha}$, but refrain from modelling the BAO residual component as described in Appendix B. We then infer the $P^{\rm pseudo}_{\rm NL}$ for this model using the emulator in the fiducial 13-dimensional setup and compare to the prediction from halofit.

The results using the 500-node training sets are shown in the lower panel of Figure 12. We find that despite the BAO residual in the shape being as large as $2\%$ for this cosmology999For comparison, this feature never exceeds 1% in our other test cosmologies., the emulator reproduces the $P^{\rm pseudo}_{\rm NL}$ to better than 1% accuracy at $z=0$ (pink) in almost all of the training set realisations. The reason that these accuracies are achievable, as with the physical models, is due to the boost factor suppressing the early-time modifications of the BAO to sub-per-cent levels. We also note that changes in the linear matter power spectrum for models with $\Delta N_{\rm eff}\neq 0$ (see Figure 12) are caused by processes happening well before the matter dominated epoch. At redshifts $\gtrsim 100$ the background and growth evolution in these extensions are indistinguishable from those in $\Lambda$CDM. Therefore, the full non-linear matter power is entirely captured by the pseudo cosmology (i.e. the reaction is unity for all scales, see Section 2.1). At $z=1$, however, we find that this model cannot be reliably predicted to within the target accuracy across the training set realisations. This is perhaps an indication that models with such marked BAO residuals are beyond the capabilities of our emulator as specified in this work. This may be improved upon given an alternative emulation method, such as sparse polynomial chaos expansion (Blatman & Sudret, 2011; Euclid Collaboration et al., 2018), or with careful optimisation of the training set configuration (Rogers et al., 2019; Caron et al., 2019).

Appendix E Emulation without de-wiggling

Although modelling of the BAO residuals in the shapes to obtain the wiggle parameters, $a$ and $b$, has negligible effect on the results (see Appendix B), we find that shape de-wiggling remains important for improving emulation accuracy. This is because the PCA is a low-pass filter, effective at capturing only the low-frequency component of the shapes. This property is imparted by the basis set, $\Phi_{i}$, derived from smooth curves (see Section 3.2.1), such that the corresponding weights, $\Delta\alpha_{i}$, are designed to contain smooth component information. Consequently, performing the PCA directly on the full shape leads to misestimation of the $\Delta\alpha_{i}$ values, and bias in the emulation.

Figure 13 presents the effect of this for the xMNU_0.4 cosmology (see Table 1) at $z=0$ (upper panel) and $z=1$ (lower panel). The pink band shows the full range of emulation accuracies achieved with the 10 realisations of the 500-node training set for our fiducial emulation method incorporating the de-wiggling, reproduced from Figure 8. The grey band shows this quantity for the case with no de-wiggling, where the PCA decomposition is performed directly on the full shape. Indeed, excess inaccuracies of $\sim$0.5% are observed in the absence of de-wiggling, for the offset PCA weights map to a boost factor related to a slightly different shape compared to the halofit predictions. Although currently subdominant, these errors will become more relevant for future emulators designed to have target performance of $\sim$0.5%, a requirement imposed by irreducible inaccuracies in the halo model reactions.