Accelerating Large-Scale-Structure data analyses by emulating Boltzmann solvers and Lagrangian Perturbation Theory

The linear matter power spectrum is an essential ingredient in all theoretical models for interpreting large-scale-structure observables. Although Boltzmann codes such as CLASS or CAMB are very efficient at computing the linear spectrum, the analysis of data usually requires 10 4-10 6 evaluations, which means this task can be the most computationally expensive aspect of data analysis. Here, we address this problem by building a neural network emulator that provides the linear theory (total and cold) matter power spectrum in about one millisecond with ≈0.2%(0.5%) accuracy over redshifts z ≤ 3 (z ≤ 9), and scales10 -4 ≤ k [ h Mpc -1] < 50. We train this emulator with more than 200,000 measurements, spanning a broad cosmological parameter space that includes massive neutrinos and dynamical dark energy. We show that the parameter range and accuracy of our emulator is enough to get unbiased cosmological constraints in the analysis of a Euclid-like weak lensing survey. Complementing this emulator, we train 15 other emulators for the cross-spectra of various linear fields in Eulerian space, as predicted by 2nd-order Lagrangian Perturbation theory, which can be used to accelerate perturbative bias descriptions of galaxy clustering. Our emulators are specially designed to be used in combination with emulators for the nonlinear matter power spectrum and for baryonic effects, all of which are publicly available at http://www.dipc.org/bacco.


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 10 4 -10 6 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.
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.

Amendments from Version 1
In this version, we have included in the text the discussions raised by the comments of the referees.We have added a Figure (Figure 5) that shows the accuracy of our emulation in the galaxy and galaxy-matter power spectra with realistic priors on the galaxy bias.
Any further responses from the reviewers can be found at the end of the article

REVISED
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 available 1 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.

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).

Parameter space
We consider 8 cosmological parameters: the primordial spectral amplitude A s and index n s ; 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 ≡ H 0 /(100 kms −1 Mpc −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) = w 0 + (1a) w a ; 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 A s and n s is straightforward at linear level, we can avoid to train the neural network for these two parameters.Thus, we train the network with fixed s A′ and s n′ , and then rescale the power spectrum as follow:

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., w 0 = -1, and w a = 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.

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 P NN is the power spectrum predicted by the neural network, P CLASS 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 10 5 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 0.2 epoch N − .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 N epochs , 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.

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.
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.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.
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, P NN , to that of the Boltzmann solver CLASS, P CLASS , 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.
In the left panel of Figure 3 we show the linear growth factor, D ( ) ( , ) / ( , 0) D k P k z P k z ≡ = 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/P no-wiggle , which highlights the contribution of the BAO to the power spectrum.We compute P no-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 A s 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.

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.

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, P gg , and the matter-galaxy cross-power spectrum, P gm , at second order are given as: where, b δ , b δ 2, b s 2, and b ∇ 2 δ are free "bias" parameters; δ is the linear density field, s 2 is the shear field defined as In other words, the clustering of galaxies is given as a weighted sum of 15 cross-spectra of five Lagrangian fields, δ L ∈ {1, δ, δ 2 , s 2 , ∇ 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 P 11 , P 1δ , 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 k s = 0.75hMpc −1 .

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 , 1s 2 , δδ, δs 2 , δ 2 ∇ 2 δ, s 2 ∇ 2 δ, and ∇ 2 δ∇ 2 δ) more accurate than 0.5%.
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 .
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.

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.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; s 2 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.
Our mock data corresponds to the auto and cross power spectra of weak lensing shear measurements for 10 equipopulated redshift bins z i ∈ [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 g i (χ) 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 with z i 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)Δℓf sky is the number of independent multipoles falling in a bin centered in ℓ with width Δℓ, f sky is the fraction of the sky covered by the survey, σ e and eff i n 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 f sky = 0.36, σ e = 0.37, eff i n = 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, w 0 , and w a .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.

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 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.
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 project 3 , 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.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.
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.

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 Y He = 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 N r = 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 & 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 w 0 w a dynamical dark energy, and finally in our full extended cosmological space.In Figure 7 we show the results of this comparison.
In the minimal ΛCDM scenario and in the w 0 w a 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 −1 hMpc −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 −3 h 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.Although BAO peaks are well recovered by emulators, there are still residual errors for example in the suppression due to massive neutrinos.Can this be improved further by sampling k-bins more densely around BAO peaks?
For the LPT emulator, it is good to see the accuracy of the final combined power spectrum with a set of reasonable bias parameters.Not all components will contribute equally, and it is good to know the accuracy of the final galaxy (and cross) power spectrum.This can be done by randomly sampling bias parameters with reasonable priors for example.
For the mock cosmic shear analysis, it is useful to know a typical improvement on the MCMC run times using emulators as opposed to using CLASS.The difference between CLASS and CAMB for massive neutrinos is larger than emulation errors.It will be good to add a short discussion on why using CLASS is preferred to construct emulators.

Is the work clearly and accurately presented and does it cite the current literature? Yes
Is the study design appropriate and does the work have academic merit?Yes

Are sufficient details of methods and analysis provided to allow replication by others? Yes
If applicable, is the statistical analysis and its interpretation appropriate?Yes Are all the source data underlying the results available to ensure full reproducibility?Yes

Are the conclusions drawn adequately supported by the results? Yes
Competing Interests: No competing interests were disclosed.
I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
first PCs describe the broad shape of the power spectrum, higher PCs are connected to the BAO oscillations.In this work, we have retained the first 64 PCs, which are sufficient to describe the broad shape of the power spectrum at 0.03%, but leave residuals at the BAO scales at 0.1% level.We have indeed chosen the number of PCs used in this work as a compromise between the accuracy of the decomposition and the complexity of the network emulation.Moreover, by adding more PCs we also include small-scale noise that the network will likely learn.In order to get accuracy below 0.1%, more PCs have to be included, at the price of a more complex training of the network (and noisier predictions).Nonetheless, if such precision will be required in the future, a better emulation can be achieved in this way.We have added this discussion to the paper.
We thank the referee for raising this point, it is indeed important to show a final accuracy for the combined cross spectra.We have added a Figure to our paper (Fig. 5) that shows the accuracy of the galaxy auto power spectrum and galaxy-matter cross power spectrum when drawing the bias parameters from realistic Lagrangian galaxy priors presented in Zennaro et al. 2021.We find a percent accuracy in both the power spectra at most of the scales considered. 2.
In our case, when considering massless neutrinos (where the NN gain is effectively lower because Boltzmann solvers are faster), by using the NN emulator we reduced the analysis times by a factor between 10 and 100 with respect to CLASS.We stress that these factors are highly dependent on the setup (the algorithm used for sampling, machine employed, number of cores).

3.
We showed that emulation techniques can be more accurate than the typical difference between CLASS and CAMB, arguably the two most used Boltzmann solvers.We highlight that even when dropping the fluid approximation of CLASS, we find a difference between the two codes of about 1% at the smallest scales considered (k = 50 h/Mpc ).Although such small scales are beyond the reach of the next generation LSS surveys, and therefore using either of the two Boltzmann solvers will likely produce the same outcome, we think that a more in-depth investigation of the differences between CLASS and CAMB within such a large cosmological space and wide range of scales is necessary.We leave such investigation for future works.We have added this discussion to the paper.

4.
Competing Interests: No competing interests were disclosed.therefore all the differences between CAMB and CLASS in this case should come from their different dynamical dark energy implementations.We plan to wait for the second referee report before uploading the new version of the text, to avoid confusion among different versions of the paper.
this procedure is practically equivalent to the emulation of the power spectrum transfer function.

Figure 1 .
Figure 1.The mean absolute fractional error of our neural network,

Figure 3 .
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.
and ∇ 2 δ the Laplacian of the linear density; and et al., 2021, for details).

Figure 4 .
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; s 2 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.

Figure 6 .
Figure6.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.

Figure 7 .
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.

Cosmology A s n s Ω c Ω b h M ν [eV] w 0 w a z
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 A s and n s 2.

of galaxy clustering and lensing from a hybrid N-body- perturbation theory model
. arXiv e-prints.art.arXiv: 2101.11014.2021.
This is an open access peer review report distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.