The BACCO Simulation Project: Exploiting the full power of large-scale structure for cosmology

The identification of bound structures is a key aspect of $N$-body simulations, thus L-Gadget3 features a version of the SUBFIND algorithm with several improvements.

The first improvement is the ability to track subhalos on-the-fly across snapshots – defining progenitor and descendants –, computing various additional quantities such as peak halo mass, peak maximum circular velocity, infall subhalo mass, and mass accretion rate, among others. These properties become useful when modelling galaxy formation within gravity-only simulations (e.g. Chaves-Montero et al., 2016; Moster et al., 2018).

The second improvement of our updated version of SUBFIND is the use of the subhalo catalogue in the previous snapshot to better identify structures. In the original algorithm, particles are first sorted according to the local density, then when a saddle point is detected, the most massive group at that point is considered as the primary structure. This, however, can cause inconsistencies across time, as small changes can lead to fluctuations in the structure considered as primary. In our version of SUBFIND instead, when a saddle point is detected, we consider as primary the substructure whose main progenitor is the most massive. This has proven to return more stable merger trees and quantities that are not local in time (e.g. peak maximum circular velocity).

Finally, during their evolution, substructures can disappear from a simulation due to artifacts in the structure finder, finite numerical resolution, or because its mass falls below the resolution limit owing to tidal stripping. The last of our improvements to SUBFIND is the ability to track all subhalos with no recognizable descendant, keeping the position and velocity of their most bound particle. This is a indispensable feature to correctly model the small-scale galaxy clustering in dark matter simulations (Guo & White, 2014).

At every output time, we store FoF groups and SUBFIND subhalos with more than 10 particles. In total, there are approximately $129$ billion groups and $214$ billion subhalos in our outputs. This means our simulations are able to resolve halos with mass $5\times 10^{10}h^{-1}{\rm M_{\odot}}$, and subhalos with a number density of $0.1\,h^{3}{\rm Mpc}^{-3}$ at $z=0$.

In order to quantify the accuracy of our results, we have carried out a suite of small, $L=64\,h^{-1}{\rm Mpc}$ box, simulations where we systematically vary the numerical parameters around those adopted in the BACCO simulations. Specifically, we consider the time-integration accuracy, force calculation accuracy, softening length, and mass resolution. We adopt the parameters of the nenya simulation as the reference point for these tests.

In Fig. 2 we compare the nonlinear power spectrum between these test runs and one that adopts the same numerical parameters and mean interparticle separation as our main BACCO simulations. Solid and dashed lines denote the results at $z=0$ and $z=1$, respectively. In all panels the grey shaded region indicates $\pm 1\%$ agreement.

In the top panel we display simulations with different mass resolutions, $N^{1/3}/L=[2,3,4]$. The main effect of mass resolution is how well small non-linear structures are resolved. We see no significant effect up to $k\sim 5h\,{\rm Mpc}^{-1}$ and a mild increase, of about 1% in the power at $k\sim 10h\,{\rm Mpc}^{-1}$ when improving the mass resolution.

The main source of inaccuracies in the force calculation are the terms neglected in the oct-tree multipole expansion. Specifically, Gadget considers only the monopole contribution down to tree-nodes of mass $M$ and size $\ell$ that fulfill $\frac{GM^{2}}{r^{2}}\left(\frac{\ell}{r}\right)>\alpha|\vec{a}|$ for a particle at a distance $r$ and acceleration $|\vec{a}|$. Thus, the accuracy in the force calculation is controlled by the parameter $\alpha$. In the second panel we vary this parameter and see that the power spectrum is converged at a sub-percent level.

As in previous versions of Gadget, time-steps in L-Gadget3 are computed individually for each particle as $\sqrt{2\eta\epsilon/|\vec{a}|}$, where $\epsilon$ is the softening and $|\vec{a}|$ is the magnitude of the acceleration in the previous timestep. The parameter $\eta$ therefore controls how accurately orbits are integrated. In the third panel we vary this parameter increasing/decreasing it by a factor of 4/2 with respect to our fiducial value, $\eta=0.05$. We see that the power spectrum varies almost negligibly with less than a $1\%$ impact up to $k\sim 10h\,{\rm Mpc}^{-1}$.

Perhaps the most important degree of freedom in a numerical simulation is the softening length, $\epsilon$, a parameter that smooths two-body gravitational interactions (note that in Gadget, forces become Newtonian at a distance $2.7\epsilon$). In the fourth panel we compare three simulations with 50% higher and lower values of $\epsilon$ to values equal to $1/140$ and $1/35-th$ of the mean interparticle separation. On small scales we see that the amount of power increases systematically the lower the value of the softening length. In particular, our fiducial configuration underestimates the power by $1(2.5)$% at $k\sim 5(10)\,h\,{\rm Mpc}^{-1}$.

In summary, our results appear converged to better than 1% up to $k\sim 5h\,{\rm Mpc}^{-1}$, and to $\sim 3\%$ up to $k\sim 10\,h\,{\rm Mpc}^{-1}$. The main numerical parameter preventing better convergence is the softening length. To mitigate its impact, we have adopted the following empirical correction to our power spectrum results

To further validate our $N$-body code and quantify the accuracy of our numerical simulations, we have carried out the main simulation of the “Euclid code comparison project”, presented in Schneider et al. (2016). This simulation consists of $2048^{3}$ particles of mass $1.2\times 10^{9}\,h^{-1}{\rm M_{\odot}}$ in a $500\,h^{-1}{\rm Mpc}$ box, and has been carried out with several $N$-body codes: RAMSES (Teyssier, 2002), PkdGrav3 (Potter et al., 2017), Gagdet3 , and recently with ABACUS by Garrison et al. (2019).

Our realisation of this simulation adopts the same force and time-integration accuracy parameters as of main BACCO simulations, and the same softening length as nenya, $\epsilon=5\,h^{-1}{\rm kpc}$. The full calculation required $6230$ timesteps and took $3\times 10^{5}$ CPU hours employing $1024$ MPI Tasks.

In Fig. 3 we compare the resulting power spectra at $z=0$. We display the ratio with respect to our L-Gadget3 results including the correction provided in Eq. 1. For comparison, we display the uncorrected measurement as a blue dashed line. Note the spikes on large scales are caused by noise due to a slightly different $k$-binning in the spectra.

Our results, PkdGrav, RAMSES, and that of ABACUS agree to a remarkable level – they differ by less than 1% up to $k\sim 5h\,{\rm Mpc}^{-1}$. This is an important verification of the absolute accuracy of our results, but it is also an important cross-validation of all these 4 $N$-body codes. Interestingly, the Gadget3 result presented in Schneider et al. (2016) is clearly in tension with the other $4$ codes. Since the underlying core algorithms in our code and in Gadget3 are the same, the difference is likely a consequence of numerical parameters adopted by Schneider et al. (2016) not being as accurate as those for the other runs. In the future, it will be important to conduct code comparison projects were numerical parameters are chosen so that each code provides results converged to a given degree.

Having validated our numerical setup, we now present an overview of the results of our BACCO simulations at $z=0$ in Fig. 4.

In the left panel we show the nonlinear power spectrum in real space. Firstly, we see the low level of random noise in our predictions owing to the “Paired-&-Fixed”  method. On large scales, we have checked that our results agree at the $0.5\%$ level with respect to the linear theory solution, which we compute for each cosmology using the Boltzmann code CLASS (Lesgourgues, 2011). Also on large scales, we see that the three simulations display significantly different power spectra, despite the three of them having identical values for $\sigma_{8}$. This is mostly a consequence of their different primordial spectral index, $n_{s}$. In contrast, on scales smaller than $k\sim 0.2h\,{\rm Mpc}^{-1}$, the spectra become much more similar as their linear spectra also do (shown by dotted lines).

In the middle panel we show the monopole of the redshift-space correlation functions for subhalos with a spatial number density equal to $10^{-3}\,h^{3}{\rm Mpc}^{-3}$ selected according to their peak maximum circular velocity. This is roughly analogous to a stellar mass selection above $5\times 10^{10}h^{-1}{\rm M_{\odot}}$. We see a significant difference in the correlation amplitude among simulations. This hints at the potential constraining power of LSS if a predictive model for the galaxy bias is available. We display each of the two “Paired-&-Fixed” simulations in each cosmology as solid and dashed lines. Unlike in the power spectrum plot, the effect of pairing the initial phase fields is visible.

Finally, in the right panel of Fig. 4 we display the mass function of halos identified by the FoF algorithm. As in previous cases, there are clear differences among cosmologies. The minimum halo mass in our simulations is $\sim 3-4\times 10^{10}\,h^{-1}{\rm M_{\odot}}$, which should suffice to model galaxies with star formation rates above $10\,h^{-1}{\rm M_{\odot}}$/year at $z\sim 1$ (Orsi & Angulo, 2018), as expected to be observed by surveys like EUCLID, DESI, or J-PAS (Favole et al., 2017).

In the next section we will employ our BACCO simulations to predict the nonlinear mass power spectrum as a function of cosmology.

In this paper we aim to cover 8 parameters of the $\Lambda$CDM model extended with massive neutrinos and dynamical dark energy. Specifically, we consider the parameter range:

We assume a flat geometry, i.e. $\Omega_{k}=0$ and $\Omega_{m}+\Omega_{w}+\Omega_{\nu}=1$. We keep fixed the effective number of relativistic species to $N_{\rm eff}=3.046$, and the temperature of the CMB $T_{\rm CMB}=2.7255$ K and neglect the impact of radiation (i.e. $\Omega_{r}=0$). Note, however, that it is relatively straightforward to relax these assumptions and include in our framework varying curvature, relativistic degrees of freedom, or other ingredients.

For comparison, we notice that the range of parameters covered by our emulator is approximately twice as large as that of the Euclid emulator (with the exception of $w_{0}$, which is similar), which implies a parameter space volume $\sim 200$ times larger, and that it covers cosmological parameters beyond the minimal $\Lambda$CDM: dynamical dark energy and massive neutrinos. In contrast, our parameter space is similar to that considered by the Mira-Titan Universe project (Heitmann et al., 2016; Lawrence et al., 2017), however, as we will discuss next, we cover the space with approximately 20 times the number of sampling points.

We now define the cosmologies with which we will train our emulator. This is usually done by sampling the desired space with a Latin-Hypercube (e.g. Heitmann et al., 2006). We adopt a slightly different strategy based on the idea of iterative emulation of Pellejero-Ibañez et al. (2020) (see also Rogers et al., 2019), where we preferentially select training points in regions of high emulator uncertainty.

Let us first define the quantity we will emulate:

Emulating the quantity $Q(k)$, rather than the full nonlinear power spectra, reduces significantly the dynamical range of the emulation, simplifying the problem and thus delivering more accurate results (see e.g. Euclid Collaboration et al., 2019; Giblin et al., 2019). We emphasise, however, that these transformations are simply designed to improve the numerical stability of our predictions – our results regarding large scales and BAO are still uniquely provided by our simulation results.

To build our training set, we first construct a Latin-Hypercube with 400 points. We then build a Gaussian Process (GP) emulator for $Q$, which provides the expectation value and variance for the emulated quantity (c.f. §3.1.7 for more details). We then evaluate the GP emulator over a Latin-Hypercube built with 2,000 points, and select the 100 points expected to have the largest uncertainty in their predictions. We add these points to our training set and re-build the GP emulator. We repeat this procedure 4 times.

We display the final set of training cosmologies, comprised of 800 points, in Fig. 5. Blue symbols indicate the initial training cosmologies, whereas orange and green symbols do so for the 2nd and 4th iteration, respectively. We can appreciate that most points in the iterations are located near the boundaries of the cosmological space, which minimise extrapolation.

It is worth mentioning that we sample our space significantly better than typical emulators. For instance, the Mira-Titan and Euclid Emulator projects employ $36$ and $100$ simulations, respectively. This implies that we can keep the emulator errors under any desired level, which could be problematic otherwise, as parameter-dependent uncertainties in models, could in principle bias parameter estimates. On the other hand, sampling points can be designed optimally so that emulation errors are much smaller than a naive Latin Hypercube (Heitmann et al., 2014; Rogers et al., 2019), however, this can be only be done optimally for a single summary statistic: what could be optimal for the matter power spectrum, is not necessarily the best for, e.g, the quadrupole of the galaxy power spectrum.

Our next step is to predict nonlinear structure for each of our training cosmologies. For this, we employ the latest incarnation of the cosmology-rescaling algorithm, originally introduced by Angulo & White (2010).

In short, the cosmology-rescaling algorithm seeks a coordinate and time transformation such that the linear rms variance coincides in the original and target cosmologies. This transformation is motivated by extended Press-Schechter arguments, and returns highly accurate predictions for the mass function and properties of collapsed objects. On large scales, the algorithm uses 2nd-order Lagrangian Perturbation Theory to modify the amplitude of Fourier modes as consistent with the change of cosmology. On small scales, the internal structure of halos is modified using physically-motivated models for the concentration-mass-redshift relation (e.g. Ludlow et al., 2016).

The accuracy of the cosmology algorithm has been tested by multiple authors (Ruiz et al., 2011; Angulo & Hilbert, 2015; Mead & Peacock, 2014; Mead et al., 2015; Renneby et al., 2018; Contreras et al., 2020), and it has been recently extended to the case of massive neutrinos by Zennaro et al. (2019). Specifically, by comparing against a suite of $N$-body simulations, Contreras et al. (2020) explicitly showed that the cosmology rescaling achieves an accuracy of $\lesssim 3\%$ up to $k=5h\,{\rm Mpc}^{-1}$ over the same range of cosmological parameter values we consider here (c.f. Eq. 3.1.1). The largest errors appear on small scales and for dynamical dark energy parameters but, when restricted to the 6 parameters of the minimal $\Lambda$CDM model, the rescaling returns $1\%$ accurate predictions. Note that the level of accuracy is set by the performance of current models for the concentration-mass-redshift relation (which usually are not calibrated for beyond $\Lambda$CDM parameters), and future progress along those lines should feedback into higher accuracy for our emulator.

For our task at hand, we first split the parameter space into three disjoint regions where nenya, narya, or vilya will be employed (Contreras et al., 2020). Then, we load in memory a given snapshot of a given simulation and then rescale it to the corresponding subset of the 800 cosmologies. We employ 10 snapshots per simulation.

The full rescaling algorithm takes approximately 2 minutes (on 12 threads using OpenMP parallelization) per cosmology and redshift. Thus, all the required computations for building the emulator required approximately $10,000$ CPU hours – a negligible amount compared to that employed to run a single state-of-the-art $N$-body simulation.

We estimate the power spectrum in our rescaled simulation outputs using Fast Fourier Transforms. We employ a cloud-in-cells assignment scheme over two interlaced grids (Sefusatti et al., 2016) of $N=1024^{3}$ points. To achieve efficient measurements at higher wavenumbers, we repeat this procedure after folding the simulation box $4$ times in each dimension. We measure the power spectrum in $160$ evenly-spaced bins in $\log{k}$, from $k=0.01$ to $5h\,{\rm Mpc}^{-1}$, and transition from the full to the folded measurement at half the Nyquist frequency. We have explicitly checked that this procedure returns sub-percent accurate power spectrum measurements over the full range of scales considered.

Although our BACCO simulations have high mass resolution, for computational efficiency, hereafter we will consider a subset of $1080^{3}$ particles (uniformly selected in Lagrangian space) as our dark matter catalogue. This catalogue, however, is affected by discreteness noise. For a Poisson sampling of $N$ points over a box of side-length $L$, the power spectrum receives a contribution equal to $(L/N)^{3}$. However, this might not be an accurate estimate for our discreteness noise since our sampling is homogeneous in Lagrangian coordinates. We have estimated its actual contribution as a third-order polynomial by comparing the spectra of the full and diluted samples at different redshifts. We found that the amplitude is proportional to $\sim 0.6$ times the Poisson noise at $z=0$, and that it progressively decreases in amplitude at higher redshifts, to reach $\sim 0.2$ times the Poisson noise at $z=2$.

All our power spectrum measurements will be corrected by discreteness noise by subtracting the term described above, and further corrected for finite numerical accuracy following section §2.3.1. However, this procedure is not perfect, and we still detect $\sim 2$% residuals at $z\sim 1$ at $k\sim 5h\,{\rm Mpc}^{-1}$. This will contribute to uncertainties in our emulator, which, however, are smaller than systematic uncertainties induced by the cosmology scaling.

In total, we employ more than 16,000 power spectrum measurements; $400+4\times 100$ training cosmologies at approximately 10 different cosmic times for two paired simulations. Grey shaded regions in Fig. 6 show the range covered by this data at $z=0$ and $z=1$, with the inset focusing on the range of wavemodes where baryonic acoustic oscillations are found. Blue and red lines show two randomly-chosen cosmologies.

On the largest scales, we can see that our rescaled simulations agree almost perfectly with linear theory. Although expected, it provides further validation of the dataset. On intermediate scales, we can see a lack of any oscillatory behaviour – better appreciated in the figure insets. This is a consequence of our $P_{\rm linear}^{\rm smeared-BAO}$ correctly capturing the nonlinear smearing of baryonic acoustic oscillations caused mostly by large-scale flows. This is true at both $z=0$ and $z=1$, which confirms the correct redshift dependence of $P_{\rm linear}^{\rm smeared-BAO}$. On small and nonlinear scales, we see an increase of more than one order of magnitude with different cosmologies differing by even factors of 3-4.

To reduce the dimensionality of our power spectrum measurements, we have performed a principal components (PC) analysis over our whole dataset, after subtracting the mean. We have kept the 6 $k$-vectors with the highest eigenvalues, which together can explain all but $10^{-3}$ of the data variance.

In the top panel of Fig. 7 we show the amplitude of these PCs, as indicated by the legend. The most important vector is a smooth function of wavenumber and roughly captures the overall impact of nonlinear evolution. Subsequent vectors describe further modifications, but none of them shows significant oscillations, which indicates that $P_{\rm linear}^{\rm smeared-BAO}$ accurately models the BAO as a function of cosmology and redshift. It is worth noting that only the 6th vector displays any noticeable noise, owing to the highly precise input dataset.

In the bottom panel of the Fig. 7 we display the ratio of the full $Q\equiv\log(P/P_{\rm linear})$ over that reconstructed using the first PCs aforementioned. We show the results for a random 2% of the power spectrum in our training data. We can see that almost all of them are recovered to better than $0.5\%$. It is interesting to note that the residuals, although increase on intermediate scales, are mostly devoid of structure, which suggests that including additional PCs will simply recover more accurately the intrinsic noise in our dataset rather than systematic dependencies on cosmology.

To confirm this idea, we have repeated our emulation but keeping twice as many PCs. Although the description of the values of Q in our dataset gained accuracy (decreasing the differences down to $\sim 0.25\%$), the performance of the emulator at predicting other cosmologies did not increase. This supports the idea that PCs beyond the 6-th are simply capturing particular statistical fluctuations in the training set, rather than cosmology-induced features. It is interesting to note that this also means that all nonlinear information of the power spectrum beyond the BAO can be described by only $6$ numbers.

To interpolate between our training data, and thus predict $Q(k,z)$ for any cosmology, we will construct two emulators: one built by employing Gaussian Processes Regression, and another one built by training a Neural Network.

In short, Gaussian processes assume that every subset of points in a given space is jointly Gaussian distributed with zero mean and covariance $K$. The covariance is a priori unknown but it can be estimated from a set of observations (e.g. our training set). Once the covariance is specified, the Gaussian process can predict the full probability distribution function anywhere in the parameter space.

In our case, we measure the amplitude associated to each PC in each training set. We then build a separate Gaussian Process for each of our PCs, using the package GPy (GPy, 2012). We assume the covariance kernel to be described by a squared exponential, with correlation length and variance set to the values that best describe the correlation among our data, found by maximising the marginal likelihood.

For regression tasks, neural networks perform remarkably well, serving as a highly-precise approximation to almost any continuous function. Thus, they provide an alternative for emulating $Q(k,z)$ (see e.g. Agarwal et al., 2014; Kobayashi et al., 2020).

In our case, we found that a feed-forward neural network with a simple architecture allowed us to obtain an accurate emulation. Specifically, we employed a fully-connected network with 2 hidden layers, each with 2000 neurons. For the activation function – which takes our input vectors and transforms them non-linearly – we use of Rectified Linear Units (ReLUs) which are commonly used in machine learning. The training of the network is done by using the adaptive stochastic optimization algorithm Adam (Kingma & Ba, 2014), with a default learning rate of $10^{-3}$. For the implementation of the neural network we made use of the open-source neural-network python library, Keras (Chollet et al., 2015).

We employ the PCA-reconstructed $Q(k,z)$ as our training data for the neural network (NN), so we can compare its results to those of the Gaussian process in equal terms. We randomly selected 5% of our power spectra as our validation set, and employ the remaining 95% as our training data. We define the loss function as the mean squared error function (MSE), and employed 2000 epochs for the training, after which our results appeared converged. At this point, the value of the loss function evaluated in the training and validation sets differed by less than $\sim 50\%$, which suggest no overfitting. We note that we also tried the so-called batch-normalization and data-dropouts, without finding significant improvements in our predictions.

For a given target cosmology in our parameter space, we can predict the amplitude associated to each PC and then reconstruct the full $Q(k,z)$ vector. We found both emulators displaying similar performances: we compared their predictions at 1000 random locations within our parameter space, and differences were typically smaller than $1-2$ percent. Both emulators display similar computational performances, about $200$ milliseconds for predicting $Q$ on any cosmology.

We present a first look at our emulator results in Fig. 8, where we show the nonlinear power spectrum at $z=0$ for a cosmology consistent with a recent analysis of Planck data. Specifically, we consider: $\sigma_{8}=0.8102$, $\Omega_{m}=0.30964$, $\Omega_{b}=0.04897$, $n_{s}=0.9665$, $h=0.6766$, $M_{\nu}=0.0$, $w_{0}=-1.0$, $w_{a}=0.0$. The left panel shows the full power spectrum, the middle panel do so relative to linear theory, and the rightmost panel relative to HaloFit calibrated by Takahashi et al. (2012).

For comparison, we include as orange circles, the results of a full $N$-body simulation with the same mass resolution and numerical parameters as those of our main BACCO suite, but in a smaller box, $L=512\,h^{-1}{\rm Mpc}$. Also as in our BACCO suite, the initial conditions have been “Paired-&-Fixed” for this simulation.

Overall, we can see that the BACCO emulator and the $N$-body simulation agree to a remarkable level, being indistinguishable by eye in the leftmost and middle panels. In particular, on large scales, both agree with linear theory (which only can be appreciated thanks to the “Paired-&-Fixed” initial conditions); on intermediate scales, both also predict a BAO featured smeared out compared to linear theory.

In the rightmost panel we can see these aspects in more detail. Firstly, we note that the simulation and emulator results agree to about 1% on all the scales considered. This is consistent with the expected accuracy of the cosmology rescaling method, but also note that due to its somewhat small volume, the cosmic variance in the $N$-body simulation results are not negligible. It is also interesting to note the systematic disagreement with the predictions from halofit. In following subsections we will explore these differences further.

To start testing the accuracy of our emulator, we have defined 10 cosmologies distributed over the target parameter space (c.f.§3.1.1) using a Latin hypercube. We then rescale our BACCO outputs to those parameters and compare the results against our emulation predictions. This essentially tests how accurate the PCA decomposition and emulation via Gaussian Process Regression are.

In Fig. 9 we show the ratio of the emulated to the rescaled nonlinear power spectrum at $z=0$ and $z=1$ for the 10 cosmologies mentioned before (we recall that neither $z=0$ nor $z=1$ were explicitly included in the emulator). We can see that on all scales our results are better than $\pm 2\%$ percent, and that most of the cosmologies are predicted to better that 1% at $z=0$. Although this is already a high accuracy, we note that more and more rescaled results can be added over time to progressively improve the quality of the emulation. We recall that this level of uncertainty is comparable to the accuracy expected by the cosmology rescaling algorithm: $1\%$ for parameters of the minimal $\Lambda$CDM and $\sim 3\%$ when considering massive neutrinos and dynamical dark energy.

We now compare our emulation results against four widely-used methods to predict the nonlinear evolution: HaloFit (Takahashi et al., 2012), the Euclid Emulator (Euclid Collaboration et al., 2019), NGenHalofit (Smith & Angulo, 2019), and Mira-Titan (Lawrence et al., 2017). Since not all of them have been calibrated over the whole parameter space covered by our emulator, we have restricted the comparison to the volume covered by the Euclid Emulator. We note that of our 800 training cosmologies, only 2 of them fall within this parameter volume.

In Fig. 10 we display the ratio of those predictions to that of our BACCO emulator. Coloured lines show 10 cosmologies set by a latin hypercube inside the Euclid Emulator parameter space. In addition, we show as a heavy line the Euclid reference cosmology: $\Omega_{\rm cdm}=0.26067$, $\sigma_{8}=0.8102$, $\Omega_{\rm b}=0.04897$, $n_{s}=0.9665$, $h=0.6766$, $M_{\nu}=0$, $w_{0}=-1$, $w_{a}=0$, employed by (Euclid Collaboration et al., 2019).

On large scales, $k<0.08h\,{\rm Mpc}^{-1}$, our emulator agrees almost perfectly with NGenHalofit and the Euclid Emulator. Halofit, on the other hand shows a small constant power deficit, whereas the Mira-Titan displays a weakly scale-dependent offset. Over this range of scales, Mira-Titan is given by TimeRG perturbation theory (Pietroni, 2008; Upadhye et al., 2014), which might indicate inaccuracies in that approach.

On intermediate scales, $0.08<k/[h\,{\rm Mpc}^{-1}]<0.5$, our results agree very well with those of NGenHalofit; for almost all cosmologies the differences are within the expected accuracy of our emulator. In contrast, we see that Halofit overestimates the amount of power by about 2% at $z=0$ and it does not correctly captures the BAO nonlinear smearing (as was already noted by Euclid Collaboration et al., 2019). On the other hand, MiraTitan underestimates the power by up to 2.5% at $z=0$, and displays a strong feature at $k\sim 0.2h\,{\rm Mpc}^{-1}$ at $z=1$. Over this range of scales MiraTitan employs a suite of low-resolution, large volume simulations, thus the disagreement could be due to an insufficient numerical accuracy in their simulations.

On small scales, $k>0.5h\,{\rm Mpc}^{-1}$, we continue to see differences among the three methods. For halofit, NGenHaloFit and Mira-Titan, they are systematic, roughly independent of cosmology, and decrease at higher redshift. A possibility is that this could be originated by the different numerical accuracy of the underlying simulations. Specifically, the simulations used by Smith & Angulo (2019) employ a softening length $\epsilon=8h^{-1}{\rm kpc}$, which according to Eq. 1 is expected to produce an underestimation of 1.5% at $k=5h\,{\rm Mpc}^{-1}$. In addition, the transition between their high and low-resolution runs occurs at $k\sim 0.6h\,{\rm Mpc}^{-1}$ which might be related to the deficit we observe at $k\sim 0.7-0.7h\,{\rm Mpc}^{-1}$ at $z=0$.

In contrast, the differences with respect to Euclid Emulator vary significantly for different cosmologies. Specifically, for the Euclid Reference Cosmology, the agreement is sub-percent on all scales. However, for our other test cosmologies, differences have a spread of $\sim 5$%, even at $k<1h\,{\rm Mpc}^{-1}$. We do not see this behavior with other methods, which might suggest that there are significant uncertainties in the emulated power spectra of the Euclid Emulator beyond their quoted precision.

Finally, we note that there is a 1% “bump” at $k\sim 3h\,{\rm Mpc}^{-1}$ in all $z=1$ panels, which is originated by our imperfect shot-noise correction.

The previous comparison was done in a rather restricted cosmological parameter volume, which served as a strong test of our accuracy. However, one of the biggest advantages of our method is the ability to predict much more extreme cosmologies even with non-standard ingredients. We provide an example of this next.

In Fig. 11 we show the predictions for the effects of massive neutrinos on the $z=0$ cold matter (baryons plus dark matter) power spectrum. We display the ratio between cases with various neutrino masses, $P_{c}^{\nu}$, relative to that without massive neutrinos $P_{c}^{\nu}$. In all cases we use $A_{s}=2.1\times 10^{-9}$ and fix all the other cosmological parameters to those of the “Euclid reference cosmology”. We only display the predictions of HaloFit and Mira-Titan, since NGenHaloFit nor the Euclid Emulator have been calibrated in the case of massive neutrinos.

We can see that on large scales, our predictions agree with both linear theory and HaloFit. For very massive neutrino cases, there is also a good agreement with MiraTitan, but there is a significant disagreement for $M_{\nu}=0.1$, this might be caused by the scale-dependent features on large scales seen in the previous figure. On intermediate scales, all predictions also agree on the broadband shape of the neutrino-induced suppression, but our emulator is also able to capture the slightly different BAO suppression expected when massive neutrinos are present. On small scales, all three methods describe the well-known neutrino-induced spoon-like suppression, disagreeing slightly on the magnitude of the maximum suppression. Note however, we expect our emulator to predict more precisely the full shape of the nonlinear power spectrum, owing to the systematic uncertainties in HaloFit and Mira-Titan discussed before.