1 Introduction
The optimal exploitation of large scale structure (LSS) data is one of the most important challenges in modern cosmology. For this, multiple theoretical models have been developed based on perturbation theory, the halo model, excursion set theory, or N-body simulations. Regardless of the nature of the modelling, essentially all approaches rely on predictions of the linear matter power spectrum as a function of cosmological parameters.
The linear matter power spectrum can be accurately computed by taking moments of the Bolztmann equation describing the co-evolution of all types of energy (neutrinos, cold dark matter, radiation, etc.) in the universe. This truncated Bolztmann hierarchy can be solved very efficiently by publicly available codes such as CMBFast (Seljak & Zaldarriaga, 1996), CAMB (Lewis et al., 2000), and CLASS (Blas et al., 2011; Lesgourgues, 2011a), which can provide the matter power spectrum in 1–10 seconds of computing, depending on the desired accuracy.
Estimating cosmological parameters from LSS observations typically requires 104–106 evaluations of the relevant theoretical model. Thus, the computational cost of estimating the linear power spectrum can reach thousands of CPU hours. Traditionally, this has not been an issue because other aspects of the modelling (e.g. the calculation of loops in perturbation theory) were significantly more expensive. This situation, however, in changing with the use of emulators. Emulators are a mathematical tool that allow for an efficient multidimensional interpolation of a given function. Popular choices in the field of LSS are Gaussian processes (e.g. Bocquet et al., 2020; Heitmann et al., 2014; McClintock & Rozo, 2019), polynomial chaos expansions (e.g. Euclid Collaboration et al., 2020; Knabenhans et al., 2019) and neural networks (e.g. Aricò et al., 2021; Kobayashi et al., 2020; Zennaro et al., 2021). In all cases, a number of (computationally expensive) predictions at different locations of a parameter space are used to build the emulator which can then be evaluated in any other point of the parameter space.
Multiple emulators have been built for LSS summary statistics (usually based on the results of N-body simulations), which include the nonlinear matter power spectrum, halo mass function, halo bias, and the effect of baryons (e.g. Angulo et al., 2021; Aricò et al., 2021; Bird et al., 2019; Bocquet et al., 2020; Euclid Collaboration et al., 2020; Heitmann et al., 2014; Knabenhans et al., 2019; McClintock et al., 2019; Schneider et al., 2020; Winther et al., 2019; Zennaro et al., 2021; Zhai et al., 2019).
Regardless the technique, the emulated quantity is usually the ratio with respect to the linear theory expectation (e.g. the nonlinear over the linear power spectrum). This reduces the dynamical range of the function and removes some of the cosmology-dependence, allowing a more robust and accurate emulation. Since the evaluation of an emulator takes a negligible amount of CPU time, the calculation of the linear power spectrum becomes the bottleneck of the analysis. Similarly to the emulators, also the widely used halofit algorithm (Smith et al., 2003; Takahashi et al., 2012) applies its interpolation functions to the linear power spectrum. In this case too, the evaluation of the linear power spectrum is by far the bottleneck of the calculation.
A possible way to tackle this is to reduce the number of evaluations needed in parameter constraints. Specifically, one technique proposed by several authors is to directly emulate the likelihood function (e.g. Leclercq, 2018; McClintock & Rozo, 2019; Pellejero-Ibañez et al., 2020). Another ideas, which we explore here, is to directly speed up the evaluation of the linear power spectrum.
Even for efficient parameter samplers, reducing the computational cost of Boltzmann solvers at fixed precision would be of great benefit. This would allow a more accurate determination of the respective likelihood functions and a faster exploration of different models and data. An option for speeding up these calculations was proposed by Albers et al. (2019) who trained a neural network to replace the most time-consuming parts of the code CLASS. Another alternative is to directly build an emulator for the linear power spectrum. This has been attempted by PICO (Fendt & Wandelt, 2007) using a fifth-order polynomial expansion to interpolate between pre-computed CMB temperature power spectra, and by CosmoNET (Auld et al., 2007; Auld et al., 2008) using neural networks to emulate CMB fluctuations and the matter power spectrum. Having been developed over a decade ago, these works were limited by the size of their training sets (which contained of the order of 1,000 measurements) and by the computational cost of the training itself.
Here we take advantage of all the recent developments in neural networks and computer architecture to develop a new emulator for the linear matter power spectrum with a focus on the analysis of forthcoming LSS experiments. Our emulator covers a eight-dimensional ΛCDM parameter space including massive neutrinos and dynamical dark energy, and reaches about 0.2% precision for 10−4 ≤ k[h Mpc−1] < 50 in a region around the parameter values preferred by current data analyses. The accuracy somewhat degrades to 0.5% when going to extreme cosmological models. Additionally, we build an emulator for the 15 different cross-spectra that enter a second-order Lagrangian bias expansion of galaxy clustering, which we compute using second-order Lagrangian perturbation theory.
Our neural networks provide predictions in about one millisecond on a single CPU core and have negligible memory requirements. Furthermore, they are designed to be employed together with our emulators for the nonlinear matter power spectrum (Angulo et al., 2021), baryonic effects (Aricò et al., 2021), and galaxy bias expansion (Zennaro et al., 2021). All these emulators are part of the baccoemu project, and they are publicly available1 and updated continuously to improve their accuracy.
The outline of this paper is the following: in section 2 we present our linear matter power spectrum and quantify its accuracy; in section 3 we introduce our Lagrangian perturbation theory spectra emulators; in section 4 we employ the linear matter power spectrum emulator to show that it provides unbiased results for the analysis of a mock stage-IV cosmic shear survey; we provide our conclusions in section 5.
2 Linear matter power spectrum emulator
In this section, we describe our emulator for the linear cold matter power spectrum. We start by defining our parameter space (subsection 2.1), then we describe our training and validation sets (subsection 2.2), as well as our neural network setup (subsection 2.3). We finish by validating our predictions for the growth function, baryonic acoustic oscillations, and the neutrino-induced suppression on the power spectrum (subsection 2.4).
2.1 Parameter space
We consider 8 cosmological parameters: the primordial spectral amplitude As and index ns; the density of cold matter and baryons in units of the critical density of the universe, Ωc and Ωb, respectively; the dimensionless Hubble parameter, h ≡ H0/(100 kms−1Mpc−1); the sum of neutrinos masses in units of eV, Mν; and a time-evolving dark energy equation of state defined through a CPL parameterisation (Chevallier & Polarski, 2001; Linder, 2003), w(z) = w0 + (1 – a) wa; the redshift z.
We define two separate cosmological hyperspaces. The first one, hereafter dubbed as standard, spans values roughly 10σ around Planck best-fitting parameters (Planck Collaboration et al., 2018). We note that this hyper-volume is very similar to that adopted in the nonlinear, baryonic, and bias emulators of Angulo et al. (2021); Aricò et al. (2021); Zennaro et al. (2021), respectively. The second cosmological parameter space, dubbed as extended, is defined by parameter ranges roughly twice as large, with which we aim at expanding the possible usage of the emulator to different kind of analyses and applications. The cosmological parameters and the respective ranges of both hyperspaces are provided in Table 1. We note that, since the dependence of the power spectrum on As and ns is straightforward at linear level, we can avoid to train the neural network for these two parameters. Thus, we train the network with fixed and , and then rescale the power spectrum as follow:
Table 1. The range of cosmological parameter values defining the two hyperspaces over which we train our neural network emulator.
As and ns are the primordial spectral amplitude and tilt, respectively; Ωc and Ωb are the density of cold matter and baryons in units of the critical density of the universe; h is the dimensionless Hubble parameter h = H0/(100 km s–1 Mpc–1);Mν is the sum of neutrinos masses in eV; and w0 and wa are parameters describing the time-evolving dark energy equation of state via w(z) = w0 + (1 – a)wa; z is the redshift.
| Cosmology | As | ns | Ωc | Ωb | h | Mν [eV] | w0 | wa | z |
|---|
| Standard | [-, -] | [-, -] | [0.23, 0.4] | [0.04, 0.06] | [0.6, 0.8] | [0.0, 0.4] | [-1.15, -0.85] | [-0.3, 0.3] | [0, 3] |
| Extended | [-, -] | [-, -] | [0.06, 0.7] | [0.03, 0.07] | [0.5, 0.9] | [0.0, 1.] | [-2., -0.5] | [-0.5, 0.5] | [0, 9] |
whit a pivot scale kp = 0.05Mpc. In this way, the neural network is facilitated by a lower dimensionality problem, we can in principle use less points in the training set, and additionally we are not limited by boundaries in As and ns2.
2.2 Training and validation sets
We define our training set by evenly sampling the parameter space with a Latin-hypercube (LH) algorithm, a statistical method which maximise the distance between the sampling points. Specifically, our neural networks are trained with the co-addition of various LH samples. The first one is a 50,000-point LH defined in the standard cosmological space. The second one is a LH of 100,000 points, defined in the extended cosmological space. We then added two LHs, of 10,000 and 20,000 points, defined respectively in the standard and extended space, but fixing the redshift at z = 0. Two extra LHs of 10,000 and 20,000 points are defined in the standard and extended space but within a ΛCDM cosmology, i.e. fixing Mν = 0., w0 = –1, and wa = 0. Finally, the last two LHs, also of 10,000 and 20,000 points, are built in the standard and extended space, within a ΛCDM cosmology and at z = 0. In this way, we improve our predictions for the ΛCDM case and at z = 0, both of which of particular cosmological interest but that might display degraded performances since they are located on the edge of the parameter hyperspace.
We note that the standard space is more densely sampled, and thus we expect our emulator to be more accurate in this region. We use 90% of the points in each LH in the training set, and add the remaining 10% to the validation set.
Overall, our training set is made of 216,000 models, whereas our validation set contains 24,000. We employ the Boltzmann solver CLASS to compute the cold-matter power spectrum in each of the 240,000 points. We make our predictions on a fixed grid in wavenumber, with 600 k–bins between 10−4 ≤ k[h Mpc−1] ≤ 50, and at the redshift indicated by the LH sampling. In Appendix A we provide details of our specific CLASS setup and a comparison with CAMB. We highlight that we emulate both the power spectrum of cold dark matter plus baryons (i.e. cold matter), and the total mass case (including massive neutrinos when present). We made this choice since the cold matter power spectrum has been shown to be a better prediction for the various LSS statistics (e.g. Castorina et al., 2015; Zennaro et al., 2019), but the total matter power spectrum can be useful, for example, for weak lensing analyses.
Before the training, we take the logarithm of the power spectra, normalising them by their mean in each k-bin. We find that, for the size of our training set, this is sufficient to reduce the variance of the data, and that using the ratio of the power spectra with some approximated method, e.g. Eisenstein & Hu (1999), does not improve sensibly the emulation. We have also tried to split the training set in two different components and emulating them separately, e.g. power spectrum at z = 0 and growth function, or smooth power spectrum and BAO oscillations, but we did not find any obvious advantage in any of these strategies.
Finally, we perform a principal component analysis (PCA) decomposition and retain the first 64 eigenvectors, which combined are enough to reproduce the broad band power spectrum below 0.03% and the BAO below 0.1%; this effectively filters out small scale noise, which aids the neural network training. We note that, in order to get the BAO accuracy below 0.1%, more PCs have to be included, at the price of a more complex training of the network.
2.3 Neural network setup
We build our emulator using a feed-forward neural network trained with 216,000 power spectrum measurements. Compared to Gaussian Processes, neural networks have the advantage of significantly better scaling of the computational and memory requirements with large training sets. We employ the publicly-available libraries Keras and TensorFlow (Abadi et al., 2016; Chollet et al., 2015) to build a relatively simple architecture with two hidden layers of 400 neurons each, and a rectified linear unit as an activation function. We use an Adam optimiser with an initial learning rate of 10−3, and define as loss function the mean absolute fractional error. This quantity, to be minimised during the training procedure, is defined as
where PNN is the power spectrum predicted by the neural network, PCLASS is the CLASS power spectrum, and the mean 〈…〉 runs over the training points.
To avoid overfitting, we monitor the loss function λ computed over both the training and the validation datasets. We stop the training when the loss function evaluated on the validation set does not decrease for more than 10,000 epochs (to avoid being stuck in local minima). We then reduce the learning rate by a factor of 10 and repeat the training until the loss function becomes flat again. This resulted into a training of approximately 105 epochs, being the final loss function of the order of 10−3. The whole training took approximately 24 hours. We apply this procedure twice, one for the total matter power spectrum and one for the cold matter one.
In Figure 1 we display the loss function λ computed in the training and validation sets, as a function of the number of epochs employed in the training of the total matter power spectrum. We see that the neural network return progressively more accurate results, with the loss function decreasing roughly as . We can also see that the accuracy is roughly identical in the training and validation sets, which suggests the network has learnt relevant features in the power spectrum data, rather than any specific source of noise. This trend is qualitatively similar for the cold matter power spectrum. We note that it appears that further training could yield to an even more accurate neural network without overfitting. However, as the improvement scales slowly with Nepochs, it becomes impractical to train for much longer. Furthermore, the accuracy of the network is comparable with the residuals of the PCA decomposition and with the level of agreement between CLASS and CAMB that we observe in Appendix A.

Figure 1. The mean absolute fractional error of our neural network, , as a function of the number of epochs employed for its training.
Black circles and red diamonds show the results when λ is evaluated on the training and the validation set, respectively. The vertical blue line marks the minimum of λ in the validation set and thus the training we will adopt thereafter in this paper.
2.4 Validation of the emulator
We test the accuracy of our emulator using the validation set described in subsection 2.2. We recall that this sample contains roughly 10% of the training data, distributed between the standard and extended cosmological spaces.
In Figure 2 we display the ratio of the total matter power spectrum computed with our emulator to that of CLASS. Shaded regions contain 68, 95, and 99.7% of the measurements in the standard and extended space (left and right panel, respectively). We display 100 randomly-selected ratios for comparison.

Figure 2. The accuracy of our neural network predictions for the linear matter power spectrum for multiple cosmologies and redshifts.
We display the ratio of the power spectra computed by our emulator, PNN, to that of the Boltzmann solver CLASS, PCLASS, in the standard (left panel) and in the extended (right panel) cosmological parameter spaces (see Table 1). In the extended space, all the cosmological parameters have values not included in the standard space. The shaded regions enclose 68%, 95%, and 99.7% of the cosmologies in our validation set, and the mean is shown as a thick black line. As an example, thin coloured lines show the results for 100 randomly selected cosmologies.
We can see that, in the standard cosmological space, our emulator is unbiased at the 0.01% level, and 68% of the validation set lies within 0.1%, while 95% within 0.2%. Outliers appears to be well within 0.4%. In the extended cosmological space, the accuracy is about ≈ 0.3% for most of the scales, although with more outliers, especially around the BAO scales. Nevetheless, we highlight the standard cosmological space is expected to contain the full range of currently allowed cosmologies. In any case, even in the extended space, our emulator accuracy is significantly higher than that of current predictions for the nonlinear matter power spectrum for which state-of-the-art N-body codes agree at ~ 2% for k ~ 10hMpc−1 (Angulo et al., 2021; Schneider et al., 2016; Springel et al., 2020). Moreover, in general we expect in the extended space an accuracy in-between the one found in the standard and extended space, being the latter a very particular case when all the parameters simultaneously have very extreme values. Therefore, any uncertainty of our emulator is arguably subdominant for LSS data analyses.
We further test our total matter power spectrum emulator in Figure 3 by examining its predictions for the redshift dependence of growth factor, the baryonic acoustic oscillations, and the neutrino-induced suppression of the matter power spectrum. In each case we illustrate the accuracy by displaying the emulator and CLASS predictions for a small number of selected cosmologies within our standard space.

Figure 3. Validation of the predictions of our neural-network emulator for the linear matter power spectrum.
Left panel: Growth factor at k = 0.2 h Mpc–1 as a function of expansion factor. Middle panel: The ratio of the linear power spectrum over its smooth (or de-wiggled) counterpart which isolates the contribution of baryonic acoustic oscillations to the power spectrum. In both of these panels we show 10 randomly-selected cosmological models within the standard space. Right panel: The ratio of the power spectrum computed including neutrinos of a various mass, as specified in the legend, over its respective neutrino massless case. In all three cases, we show the results obtained with CLASS as solid lines, and with our emulator as dashed lines. In the bottom panels, we display the compare these predictions indicating differences of 0.1% and 0.2% as shaded regions.
In the left panel of Figure 3 we show the linear growth factor, over the Einstein-de-Sitter solution, D ∝ a, computed using our emulator and CLASS at 10 different cosmological models. We see that the emulator predictions are very accurate, within 0.1%, as expected from the performance shown previously. Also, the accuracy does not depend on the expansion factor, which indicates that the performance of our emulator is the same at all redshifts.
Given their importance for LSS analyses, we explicitly check how well the baryonic acoustic oscillation (BAO) feature is recovered. In the middle panel of Figure 3 we display the ratio between linear power spectrum and its de-wiggled counterpart, P/Pno–wiggle, which highlights the contribution of the BAO to the power spectrum. We compute Pno–wiggle by performing a discrete sine transform of the linear theory power spectrum, smoothing the result, and returning to Fourier space with an inverse sine transform (e.g Baumann et al., 2018). We can see how each oscillation is remarkably well reproduced by our emulator, with small in-phase deviations of the order of 0.3%. This suggests the accuracy of our emulator is also high enough for BAO analyses.
To close this section, we examine the dependency of the power spectrum on massive neutrinos. The rightmost panel of Figure 3 shows the ratio of the power spectrum computed with increasingly massive neutrinos, Mν = [0.06, 0.1, 0.2, 0.3] eV, over the massless case, Mν = 0 eV, keeping fixed the primordial spectral amplitude As and the matter density Ωm. As with our previous tests, we find that the distortion caused by neutrinos, a suppression at small scales proportional to the neutrino mass fraction, is accurately recovered by our emulator at 0.2% level.
In this section, we have shown the tests carried out to validate the accuracy of our total matter power spectrum emulator. We have performed similar tests for the cold matter power spectrum emulator, which we do not include here for the sake of brevity, finding a similar level of accuracy.
3 Lagrangian perturbation theory emulator
Besides predicting the clustering of matter, another challenge for present-day and forthcoming cosmological surveys is to predict the clustering of galaxies. In this section we build emulators for several matter statistics as predicted by Lagrangian perturbation theory (LPT) that are typically employed in models for the power spectrum of galaxies.
3.1 Lagrangian bias expansion
Among the many possible models to describe the clustering of galaxies, a particularly promising one is a perturbative Lagrangian bias expansion (Matsubara, 2008). In this formalism, the power spectrum of galaxies, Pgg, and the matter-galaxy cross-power spectrum, Pgm, at second order are given as:
where, bδ, bδ2, bs2, and b∇2δ are free “bias” parameters; δ is the linear density field, s2 is the shear field defined as s2 ≡ sijsij where sij ≡ ∂i∂j∇–2δ–δD,ijδ(q), and ∇2δ the Laplacian of the linear density; and (see Zennaro et al., 2021, for details).
In other words, the clustering of galaxies is given as a weighted sum of 15 cross-spectra of five Lagrangian fields, δL ∈ {1, δ, δ2, s2, ∇2δ}, advected to Eulerian coordinates:
where δD is a Dirac’s delta. The displacement field Ψ(q) (and thus the cross-spectra) can be computed perturbatively at a given order, or, as proposed recently, measured directly from N-body simulations. The latter can accurately describe the power spectrum of galaxies down to much smaller scales than using perturbative solutions, which opens up the possibility of an accurate modelling of galaxy clustering (Modi et al., 2020).
In fact, Zennaro et al. (2021) and Kokron et al. (2021) have built emulators for these spectra as a function of cosmology from simulation suites. However, even in this case, large scales are described with perturbation theory owing to the noise (cosmic variance) in the N-body displacements. Additionally, the emulation is performed for the ratio of the N-body spectra over the perturbative solution, which, as for the matter power spectrum, improves the quality of the emulation.
The LPT terms can be computed with publicly available codes (such as velocileptors, Chen et al., 2020) or, as in Zennaro et al. (2021), by solving the three dimensional LPT integrals with a adaptive quadrature algorithm. Both implementations are very efficient and produce the LPT spectra in around 1 second of computing time, depending on accuracy parameters and machine architecture. However, this can become the most costly step in the context of obtaining the galaxy model from an emulator.
In this section, we build an emulator for each of the 15 LPT terms to speed up this process. Specifically, we emulate the terms P11, P1δ, and Pδδ where we expand δ at first order in Lagrangian perturbation theory; this is enough for the applications mentioned. All remaining terms are computed expanding densities at second order in Lagrangian perturbation theory, retaining only contributions of order (11) and (22). We refer the reader to Zennaro et al. (2021) for the explicit expression of each of these spectra in LPT. It is important to note that we assume the linear power spectrum entering our LPT calculations to be smoothed on a scale ks = 0.75hMpc−1.
3.2 LPT emulator
We build the LPT emulators in the standard hyper-parameter space (see Table 1) and over the expansion factor range: a ∈ [0.4, 1.] (note that this redshift range is slightly smaller than that employed in subsection 2.1). We note that our cosmological parameter space matches that used by Zennaro et al. (2021), where we presented non-linear boost-factors for these spectra. We then compute the LPT predictions in 46 logarithmic bins over the range 10−2 < k[hMpc−1] < 0.75. We sample the parameter space with a single LH with 20,000 points. In each of these points, we compute the 15 Lagrangian cross power spectra as described above. As noted by Zennaro et al. (2021), some of these spectra can be negative. In the case where a term has at least a negative value in the training set, we emulate the quantity:
so that the spectra are always positive. As in the case of the linear matter power spectrum, prior to the training we take the logarithm of the power spectra, and subtract the mean in each k-bin. We then perform a PCA decomposition prior to the training, retaining for each spectrum a number of PC sufficient to recover the power spectrum at 0.01%, which ranges from five to 20 PCA vectors.
We use 90% of the sample as our training dataset and the remaining 10% as our validation set. The architecture of the neural network used to emulate the LPT power spectra is the same as that described in subsection 2.3, i.e. 2 hidden layers with 400 neurons each.
We quantify the accuracy of our emulators in Figure 4, which shows the ratio of the emulation prediction in the validation set over the spectra computed with LPT. Overall, the accuracy of the 15 spectra is better than 1%, with some terms (1δ2, 1s2, δδ, δs2, δ2∇2δ, s2∇2δ, and ∇2δ∇2δ) more accurate than 0.5%.

Figure 4. The accuracy of our emulators for the cross-spectrum of linear fields in Eulerian coordinates predicted by Lagrangian Perturbation Theory.
The two fields defining the cross-spectra are indicated in the legend of each panel, where 1 is an homogeneous Lagrangian field; δ and δ2 are the linear density field and its square, respectively; s2 is the shear field; and ∇2δ is the Laplacian of the linear density field. In each panel we display the ratio of the emulator prediction over the same quantity computed by directly solving the relevant LPT expression. Shaded regions enclose 68% and 95% of the measurements in our validation set, and the mean is is marked by the thick black line. For comparison we show a randomly-selected set of cosmologies as coloured lines.
When computing the power spectrum of a biased tracer (e.g. galaxies), the 15 emulated cross-spectra are weighted by the respective bias factors, and thus have different contribution to the total spectrum. To quantify the accuracy of the emulation in the case of realistic galaxies, we randomly sample the priors on Lagrangian galaxy bias found in Zennaro et al. (2021). Figure 5 shows the accuracy we get in the galaxy auto power spectrum and in the galaxy - matter cross power spectrum, that is around 1% at most of the scales considered, and approaching 2% at scales k > 0.4 h Mpc−1.

Figure 5. The accuracy of our emulators for the galaxy auto power spectrum (left panel) and galaxy - matter cross-spectrum predicted by Lagrangian Perturbation Theory.
We used bias factors randomly drawn from the priors computed in Zennaro et al. (2021). In each panel we display the ratio of the emulator prediction over the same quantity computed by directly solving the relevant LPT expression. Shaded regions enclose 68% and 95% of the measurements in our validation set, and the mean is is marked by the thick black line. For comparison we show a randomly-selected set of cosmologies as coloured lines.
This accuracy is higher than that of the nonlinear emulators, thus the LPT emulation uncertainty should add a subdominant contribution to the global modelling error. Nevertheless, as with our linear emulator, it is straightforward to add further points to the training set, if higher accuracy is required.
4 Parameter constraints from mock cosmic shear power spectra
In this section, we illustrate the use of our linear matter power spectrum emulator with a simple application: we infer cosmological parameters from a mock lensing power spectrum using either CLASS or our emulator. In this way, we will confirm that the accuracy of our emulator is adequate in the context of LSS data analysis.
Our mock data corresponds to the auto and cross power spectra of weak lensing shear measurements for 10 equipopulated redshift bins zi ∈ [0.1, 2.5]. We adopt the specifics of a stage-IV survey, and in particular those of the Euclid mission (Amendola et al., 2018). The cross-spectrum of cosmic shear is given by:
where P(k,z) is the linear (total) matter power spectrum (given by CLASS for our mock data), and gi(χ) is the the lensing kernel of the i-th redshift bin:
where c is the light speed and χ(z) is the comoving distance to z. We assume a redshift distribution of galaxies in the i-th bin ni(z) ∝ z2 with zi roughly setting the redshift where the galaxy number density peaks. The galaxy distribution in each bin is normalised such that
For simplicity, we only consider the Gaussian contribution to the covariance:
where Nℓ = (2ℓ + 1)Δℓfsky is the number of independent multipoles falling in a bin centered in ℓ with width Δℓ, fsky is the fraction of the sky covered by the survey, σe and are the RMS ellipticity and the effective projected number density of the source galaxies in the i bin, respectively (Barreira et al., 2018). The presence of the Kronecker deltas δℓi,ℓj forces the Gaussian term to be diagonal, and δij to have shape noise terms only for matching redshift bins. We adopt a fiducial setup of fsky = 0.36, σe = 0.37, = 30 arcmin−2.
We fit the mock data over 20 logarithmically-spaced multipoles between ℓ ∈ [20, 5000], using the affine invariant MCMC sampler emcee (Foreman-Mackey et al., 2013) and employing 10 walkers of 10,000 steps each, considering a burn-in phase of 1,000 steps. Our data model is given by Equation 6 computed with a linear power spectrum provided by either CLASS or by our emulator, and letting free the cosmological parameters Ωc, σ8, h, w0, and wa. Using 1 CPU, a single call to CLASS takes ≈ 0.5 seconds (≈ 7 seconds when having massive neutrinos), whereas the emulator evaluation takes ≈ 1 milliseconds.
In Figure 6, we display the posterior distributions estimated with the MCMC analysis using CLASS and the emulator (black and red lines, respectively). We see that by using the linear emulator, we recover unbiased parameters at less than 0.05σ. Moreover, the 1D marginalised PDFs, the parameters degeneracy, and the contours are all almost indistinguishable from those using directly CLASS. Overall, the results show a remarkable agreement between both approaches, which suggests that the accuracy of our linear emulator is sufficient for the analyses of forthcoming LSS surveys.

Figure 6. The marginalised posterior distributions on cosmological parameters obtained from mock weak lensing data.
Our mock data is comprised by the shear power spectrum and cross-spectra using 10 tomographic bins over z ∈ [0., 2.5], and a redshift distribution of background galaxies, shape noise, and cosmic variance expected for a stage-IV Euclid-like survey. The contours in each panel show 1,2, and 3 σ levels obtained by performing an MCMC analysis using the linear matter power spectrum directly provided by the Boltzmann solver CLASS (black lines and contours) or by our emulator (red lines and contours). Blue lines indicate the cosmological parameters adopted in our mock data, whereas black and red lines show the best-fitting values obtained by our analysis.
5 Summary
In this paper, we have presented and validated a set of fast and accurate emulators aimed at speeding up the analysis of forthcoming LSS data. These emulators can provide their predictions in about one millisecond of computing time and cover a broad 8-dimensional cosmological parameter space that includes massive neutrinos and dynamical dark energy. First, we have built and validated an emulator for the linear cold matter power spectrum, which is faster than a typical Boltzmann solver by a factor of 1000. The accuracy of the emulator is subpercent over 10−4 < k < 50 (Figure 2) and it also accurately predicts the growth of fluctuations, the baryonic acoustic oscillations, and the suppression of clustering induced by massive neutrinos (Figure 3). Our second set of emulators predict multiple cross-spectra of Lagrangian fields relevant for a second-order perturbative model of galaxy bias. We compute these fields with second-order Lagrangian perturbation theory which can then be emulated with percent accuracy (Figure 4). Finally, we have shown that the accuracy of the linear emulator can be used to provide unbiased cosmological constraints for a tomographic analysis of a Euclid-like weak lensing survey (Figure 6).
The emulators presented in this work are part of the baccoemu project3, and have been specifically designed to be used with those we have previously built for the nonlinear matter power spectrum (Angulo et al., 2021), the modifications induced by baryonic physics (Aricò et al., 2021), and galaxy bias (Zennaro et al., 2021). All together, they can contribute towards a comprehensive, fast, and accurate cosmological exploitation of LSS data.
6 Data and software availability
6.1 Underlying data
Training and validation sets of the linear matter and Lagrangian bias power spectra presented.
This project contains the following underlying data:
linear_emulator_training_set.tar. Set of 240,000 linear matter power spectra obtained using the Boltzmann solver CLASS, and used as training and validation sets for the artificial neural network.
lpt_emulator_training_set.tar Set of 20,000 Lagrangian bias power spectra, used as training and validation sets for the artificial neural network.
Data are provided as pickle files readable in python, which contain dictionaries organised as described in the included README files. Data are available under the terms use of MIT license.
Acknowledgements
We warmly thank Jonás Chaves-Montero, Sergio Contreras, Marcos Pellejero-Ibañez, and Jens Stücker, for useful discussions.
A Appendix: Class setup
In this Appendix we provide some details on our CLASS setup and compare its predictions against those of CAMB.
In our CLASS calculations, we set the primordial helium fraction to YHe = 0.24, the optical depth at reionisation τreio = 0.0952, the number of (degenerate) massive neutrinos to Nν = 3, and the number of relativistic species Nr = 3.046. The neutrino to photon temperature is computed accounting for non-instantaneous neutrino decoupling and spectral distortions induced by the reheating, Tν = 0.71611. We assume the standard CLASS neutrino fluid approximation, and solve the neutrino integrals with an automatic choice of the best quadrature method, employing 150 momentum bins. To get an estimate of the absolute accuracy of the CLASS power spectra, we make a comparison against the predictions of another Boltzmann solver, CAMB (Lewis & Bridle, 2002). We use in CAMB, whenever possible, an analogous setup, specifying a degenerate neutrino hierarchy, and furthermore setting to True the accurate_massive_neutrino_transfer and Reionization options.
We compare both codes in 100 points distributed as a Latin hyper-cube, in a pure ΛCDM scenario, adding massive neutrinos, in a ΛCDM plus w0wa dynamical dark energy, and finally in our full extended cosmological space. In Figure 7 we show the results of this comparison.

Figure 7. Comparison of the linear power spectrum provided by CLASS and CAMB, two independent Bolztmann solvers.
Each panel displays the ratio for multiple cosmologies in the minimal ΛCDM model, ΛCDM plus neutrinos, ΛCDM plus dynamical dark energy, and in ΛCDM plus neutrinos and dynamical dark energy. Shaded regions denote 0.1 and 0.02% differences.
In the minimal ΛCDM scenario and in the w0wa cosmologies, the two codes agree at the 0.1% level, except for k > 10h Mpc−1 where CAMB underestimates the CLASS solution by 0.2%, likely due to CAMB neglecting the impact of reionisation on the baryon sound speed (Lesgourgues, 2011b). When considering the massive neutrinos, the agreement degrades to a 0.5%, a value larger than what found in Lesgourgues & Tram (2011) but in agreement with Zennaro et al. (2017), and approaches 1% at the smallest scaled employed. When dropping the CLASS fluid approximation, the agreement between the two codes improves to 0.2% at k < 10−1hMpc−1, but get increasingly worse at smaller scales, reaching 1% at scales k ≈ 40hMpc−1. Within the dynamical dark energy models, the two codes agree mostly at 0.1%, except for scales k < 10−3h Mpc−1 and k > 30h Mpc−1, where the difference can reach 1%. Although the scales probed by the next generation LSS surveys are well recovered by both the Boltzmann solvers considered here, a more in-depth investigation of the differences between CLASS and CAMB within a large cosmological space and wide range of scales is necessary, since we find percent differences that might potentially impact future analyses. We leave such investigation for future works.
Faculty Opinions recommendedReferences
- Abadi M, Barham P, Chen J, et al.:
Tensorflow: A system for large-scale machine learning. In 12th {USENIX} Symposium on Operating Systems Design and Implementation ({OSDI} 16). 2016; 265–283. Reference Source
- Albers J, Fidler C, Lesgourgues J, et al.:
CosmicNet. Part I. Physics-driven implementation of neural networks within Einstein-Boltzmann Solvers.
J Cosmol Astropart P.
2019; 2019(9): 028. Publisher Full Text
- Amendola L, Appleby S, Avgoustidis A, et al.:
Cosmology and fundamental physics with the Euclid satellite.
Living Rev Relativ.
2018; 21(1): 2. PubMed Abstract
| Publisher Full Text
| Free Full Text
- Angulo RE, Zennaro M, Contreras S, et al.:
The BACCO simulation project: exploiting the full power of large-scale structure for cosmology.
Mon Not R Astron Soc.
2021; 507(4): 5869–5881. Publisher Full Text
- Aricò G, Angulo RE, Contreras S, et al.:
The BACCO simulation project: a baryonification emulator with neural networks.
Mon Not R Astron Soc.
2021; 506(3): 4070–4082. Publisher Full Text
- Auld T, Bridges M, Hobson MP, et al.:
Fast cosmological parameter estimation using neural networks.
Mon Not R Astron Soc.
2007; 376(1): L11–L15. Publisher Full Text
- Auld T, Bridges M, Hobson MP:
COSMONET: fast cosmological parameter estimation in non-flat models using neural networks.
Mon Not R Astron Soc.
2008; 387(4): 1575–1582. Publisher Full Text
- Barreira A, Krause E, Schmidt F:
Accurate cosmic shear errors: do we need ensembles of simulations?
J Cosmol Astropart P.
2018; 2018(10): 053. Publisher Full Text
- Baumann D, Green D, Wallisch B:
Searching for light relics with large-scale structure.
J Cosmol Astropart P.
2018; 2018(8): 029. Publisher Full Text
- Bird S, Rogers KK, Peiris HV, et al.:
An emulator for the Lyman-α forest.
J Cosmol Astropart Phys.
2019; 2019(2): 050. Publisher Full Text
- Blas D, Lesgourgues J, Tram T:
The Cosmic Linear Anisotropy Solving System (CLASS). Part II: Approximation schemes.
J Cosmol Astropart Phys.
2011; 2011(7): 034. Publisher Full Text
- Bocquet S, Heitmann K, Habib S, et al.:
The Mira-Titan Universe. III. Emulation of the Halo Mass Function.
Astrophys J.
2020; 901(1): 5. Publisher Full Text
- Castorina E, Carbone C, Bel J, et al.:
DEMNUni: the clustering of large-scale structures in the presence of massive neutrinos.
J Cosmol Astropart P.
2015; 2015(7): 043. Publisher Full Text
- Chen SF, Vlah Z, White M:
Consistent modeling of velocity statistics and redshift-space distortions in one-loop perturbation theory.
J Cosmol Astropart P.
2020; 2020(7): 062. Publisher Full Text
- Chevallier M, Polarski D:
Accelerating Universes with Scaling Dark Matter.
Int J Mod Phys D.
2001; 10(2): 213–223. Publisher Full Text
- Chollet F, et al.:
Keras. 2015. Reference Source
- Eisenstein DJ, Hu W:
Power Spectra for Cold Dark Matter and Its Variants.
Astrophys J.
1999; 511(1): 5–15. Publisher Full Text
- Euclid Collaboration,
Knabenhans M, Stadel J, et al.:
Euclid preparation: IX. EuclidEmulator2 -- Power spectrum emulation with massive neutrinos and self-consistent dark energy perturbations.
arXiv e-prints.
art. arXiv: 2010.11288, 2020. Reference Source
- Fendt WA, Wandelt BD:
Pico: Parameters for the Impatient Cosmologist.
Astrophys J.
2007; 654(1): 2–11. Publisher Full Text
- Foreman-Mackey D, Hogg DW, Lang D, et al.:
emcee: The MCMC Hammer.
Publ Astron Soc Pac.
2013; 125(925): 306. Publisher Full Text
- Heitmann K, Lawrence E, Kwan J, et al.:
The Coyote Universe Extended: Precision Emulation of the Matter Power Spectrum.
Astrophys J.
2014; 780(1): 111. Publisher Full Text
- Knabenhans M, Stadel J, Marelli S, et al.:
Euclid preparation: II. The EUCLIDEMULATOR - a tool to compute the cosmology dependence of the nonlinear matter power spectrum.
Mon Not R Astron Soc.
2019; 484(4): 5509–5529. Publisher Full Text
- Kobayashi Y, Nishimichi T, Takada M, et al.:
Accurate emulator for the redshift-space power spectrum of dark matter halos and its application to galaxy power spectrum.
Phys Rev D.
2020; 102(6): 063504. Publisher Full Text
- Kokron N, DeRose J, Chen SF, et al.:
The cosmology dependence of galaxy clustering and lensing from a hybrid N-body-perturbation theory model. arXiv e-prints. art. arXiv: 2101.11014. 2021. Reference Source
- Leclercq F:
Bayesian optimization for likelihood-free cosmological inference.
Phys Rev D.
2018; 98(6): 063511. Publisher Full Text
- Lesgourgues J:
The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview.
arXiv e-prints.
art. arXiv: 1104.2932, 2011a. Reference Source
- Lesgourgues J:
The Cosmic Linear Anisotropy Solving System (CLASS) III: Comparision with CAMB for LambdaCDM. arXiv e-prints. art. arXiv: 1104.2934. 2011b. Reference Source
- Lesgourgues J, Tram T:
The Cosmic Linear Anisotropy Solving System (CLASS) IV: efficient implementation of non-cold relics.
J Cosmol Astropart Phys.
2011; 2011(9): 032. Publisher Full Text
- Lewis A, Bridle S:
Cosmological parameters from CMB and other data: A Monte Carlo approach.
Phys Rev D.
2002; 66: 103511. Publisher Full Text
- Lewis A, Challinor A, Lasenby A:
Efficient Computation of Cosmic Microwave Background Anisotropies in Closed Friedmann-Robertson-Walker Models.
Astrophys J.
2000; 538(2): 473–476. Publisher Full Text
- Linder EV:
Exploring the Expansion History of the Universe.
Phys Rev Lett.
2003; 90(9): 091301. PubMed Abstract
| Publisher Full Text
- Matsubara T:
Nonlinear perturbation theory with halo bias and redshift-space distortions via the Lagrangian picture.
Phys Rev D.
2008; 78(8): 083519. Publisher Full Text
- McClintock T, Rozo E:
Reconstructing probability distributions with Gaussian processes.
Mon Not R Astron Soc.
2019; 489(3): 4155–4160. Publisher Full Text
- McClintock T, Rozo E, Becker MR, et al.:
The Aemulus Project. II. Emulating the Halo Mass Function.
Astrophys J.
2019; 872(1): 53. Publisher Full Text
- Modi C, Chen SF, White M:
Simulations and symmetries.
Mon Not R Astron Soc.
2020; 492(4): 5754–5763. Publisher Full Text
- Pellejero-Ibañez M, Angulo RE, Aricó G, et al.:
Cosmological parameter estimation via iterative emulation of likelihoods.
Mon Not R Astron Soc.
2020; 499(4): 5257–5268. Publisher Full Text
- Planck Collaboration, Aghanim N, Akrami Y, et al.:
Planck 2018 results. VI. Cosmological parameters. arXiv e-prints. 2018. Publisher Full Text
- Schneider A, Stoira N, Refregier A, et al.:
Baryonic effects for weak lensing. Part I. Power spectrum and covariance matrix.
J Cosmol Astropart P.
2020; 2020(4): 019. Publisher Full Text
- Schneider A, Teyssier R, Potter D, et al.:
Matter power spectrum and the challenge of percent accuracy.
J Cosmol Astropart P.
2016; 4: 047. Publisher Full Text
- Seljak U, Zaldarriaga M:
A Line-of-Sight Integration Approach to Cosmic Microwave Background Anisotropies.
Astrophys J.
1996; 469: 437. Publisher Full Text
- Smith RE, Peacock JA, Jenkins A, et al.:
Stable clustering, the halo model and non-linear cosmological power spectra.
Mon Not R Astron Soc.
2003; 341(4): 1311–1332. Publisher Full Text
- Springel V, Pakmor R, Zier O, et al.:
Simulating cosmic structure formation with the GADGET-4 code. arXiv e-prints. art. arXiv: 2010.03567. 2020. Reference Source
- Takahashi R, Sato M, Nishimichi T, et al.:
Revising the Halofit Model for the Nonlinear Matter Power Spectrum.
Astrophys J.
2012; 761(2): 152. Publisher Full Text
- Winther HA, Casas S, Baldi M, et al.:
Emulators for the nonlinear matter power spectrum beyond Λ CDM.
Phys Rev D.
2019; 100(12): 123540. Publisher Full Text
- Zennaro M, Angulo RE, Aricò G, et al.:
How to add massive neutrinos to your ΛCDM simulation – extending cosmology rescaling algorithms.
Mon Not R Astron Soc.
2019; 489(4): 5938–5951. Publisher Full Text
- Zennaro M, Angulo RE, Pellejero-Ibáñez M, et al.:
The BACCO simulation project: biased tracers in real space.
arXiv e-prints.
art. arXiv: 2101.12187, 2021. Reference Source
- Zennaro M, Bel J, Villaescusa-Navarro F, et al.:
Initial conditions for accurate N-body simulations of massive neutrino cosmologies.
Mon Not R Astron Soc.
2017; 466(3): 3244–3258. Publisher Full Text
- Zhai Z, Tinker JL, Becker MR, et al.:
The Aemulus Project. III. Emulation of the Galaxy Correlation Function.
Astrophys J.
2019; 874(1): 95. Publisher Full Text