High-accuracy emulators for observables in $\Lambda$CDM, $N_\mathrm{eff}$, $\Sigma m_\nu$, and $w$ cosmologies

We use the emulation framework CosmoPower to construct and publicly release neural network emulators of cosmological observables, including the Cosmic Microwave Background (CMB) temperature and polarization power spectra, matter power spectrum, distance-redshift relation, baryon acoustic oscillation (BAO) and redshift-space distortion (RSD) observables, and derived parameters. We train our emulators on Einstein-Boltzmann calculations obtained with high-precision numerical convergence settings, for a wide range of cosmological models including $\Lambda$CDM, $w$CDM, $\Lambda$CDM+$N_\mathrm{eff}$, and $\Lambda$CDM+$\Sigma m_\nu$. Our CMB emulators are accurate to better than 0.5% out to $\ell=10^4$ which is sufficient for Stage-IV data analysis, and our $P(k)$ emulators reach the same accuracy level out to $k=50 \,\, \mathrm{Mpc}^{-1}$, which is sufficient for Stage-III data analysis. We release the emulators via an online repository CosmoPower Organisation, which will be continually updated with additional extended cosmological models. Our emulators accelerate cosmological data analysis by orders of magnitude, enabling cosmological parameter extraction analyses, using current survey data, to be performed on a laptop. We validate our emulators by comparing them to CLASS and CAMB and by reproducing cosmological parameter constraints derived from Planck TT, TE, EE, and CMB lensing data, as well as from the Atacama Cosmology Telescope Data Release 4 CMB data, Dark Energy Survey Year-1 galaxy lensing and clustering data, and Baryon Oscillation Spectroscopic Survey Data Release 12 BAO and RSD data.


INTRODUCTION
Next-generation Cosmic Microwave Background (CMB) surveys, such as the Simons Observatory (SO; Simons Observatory 2019), CMB-S4 (Abazajian et al. 2019;CMBS4 2022), and CMB-HD (Sehgal et al. 2019), will provide high signal-to-noise measurements of the CMB temperature and polarization anisotropy down to small scales, characterized by their power spectra. The final seasons of the Atacama Cosmology Telescope (ACT; Henderson et al. 2016) and of the South Pole Telescope (SPT; Benson et al. 2014) will also make progress in that direction in the forthcoming years. In addition, current and next-generation galaxy surveys including the Dark Energy Survey (DES; Dark Energy Survey Collaboration 2016), the Kilo-Degree Survey (KiDS; Kuĳken et al. 2015), the Hyper Suprime-Cam Survey (HSC; Aihara et al. 2018), the Vera C. Rubin Observatory (Ivezić et al. 2019), Euclid (Laureĳs et al. 2011), Roman (Spergel et al. 2015;Akeson et al. 2019) or SPHEREx (Doré et al. 2014) probe the matter power spectrum via weak gravitational lensing and galaxy clustering, and spectroscopic galaxy surveys such as the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013)  Einstein-Boltzmann codes such as (Lewis et al. 2000) and (Lesgourgues 2011a;Blas et al. 2011) are routinely used to accurately compute linear-theory cosmological power spectra, as well as the background cosmic evolution (e.g., the distance-redshift relation). These codes also implement various prescriptions for the calculation of non-linear corrections to the matter power spectrum and quantities derived therefrom, e.g., Halofit (Smith et al. 2003) and HMCode (Mead et al. 2015). The computation of these cosmological quantities through a Boltzmann code represents a significant computational requirement in Markov Chain Monte Carlo (MCMC) methods, traditionally used to extract constraints on cosmological parameters. MCMC methods require ∼ 10 4 − 10 6 calls to a Boltz-mann code as they achieve convergence through a similar number of likelihood evaluations. The problem is exacerbated if one requires high numerical accuracy in the computation of the cosmological quantities obtained from Boltzmann codes, as needed for upcoming surveys (e.g., McCarthy et al. 2022).
Several groups have developed emulators of cosmological quantities in order to accelerate the MCMC calculations by bypassing the call to a Boltzmann code with a faster algorithm (e.g. Fendt & Wandelt 2007;Auld et al. 2007;Albers et al. 2019;Mootoovaloo et al. 2022;Aricò et al. 2021;Günther et al. 2022). Recently, Spurio Mancini et al. (2022) developed C P , an emulation framework for cosmological quantities based on TensorFlow neural networks (Abadi et al. 2015). The accuracy of the C P emulators presented in Spurio Mancini et al. (2022) was tested on CMB and large-scale structure analysis of current and future data. More recently, Balkenhol et al. (2022) constructed C P emulators for the analysis of the SPT data, covering angular scales up to ℓ = 3000.
With this suite of emulators, we are able to accurately reproduce marginalized posterior distributions from current CMB and largescale structure (LSS) likelihood codes (based on CMB, matter power spectra, and BAO distances) within minutes on a laptop, whereas the original MCMC analyses took hours to days on large computing clusters. To accomplish this, we use the MCMC sampler and its implemented likelihoods (Torrado & Lewis 2021). Our C -P emulators are wrapped into via a simple wrapper that will be available online 1 . We reproduce parameter constraints from the ACT Data Release 4 (Aiola et al. 2020;Choi et al. 2020), Planck 2018 temperature, polarization, and lensing potential data (Planck Collaboration 2020b,a), DES-Y1 (Krause et al. 2017;Troxel et al. 2018;DES Collaboration 2018), and BAO distances from BOSS and other surveys (Beutler et al. 2011;Ross et al. 2015;Alam et al. 2017).
Our emulators are explicitly constructed with sufficient accuracy to remain valid for use in upcoming CMB analyses for the foreseeable future, including ACT, SPT, SO, CMB-S4, and CMB-HD, as well as ongoing galaxy surveys.
In addition to the emulators, we release the full pipeline to generate them, which can be used by the community to build emulators for other extensions to the standard cosmological model. The cosmological models covered here include Λ Cold Dark Matter (CDM), CDM, where is the dark energy equation of state; ΛCDM with varying eff , where eff is the effective number of relativistic species at recombination; and ΛCDM with varying neutrino mass Σ . In ΛCDM, ΛCDM+ eff , and CDM, we assume one massive and two massless neutrino states, as in the baseline Planck analyses (Planck Collaboration 2020a). In ΛCDM+Σ , we generate emulators for both the case of three species of neutrinos with degenerate mass, and the case of one massive and two massless neutrino states. This paper is organized as follows. In Section 2 we outline our 1 The wrapper will be shared in a forthcoming paper which will cover in detail how to embed these emulators in future experiments' likelihood pipelines. methodology, while in Section 3 we study and calibrate accuracy settings of Boltzmann solvers by comparing and . This sets the accuracy target for our emulators, which we validate in Section 4. In Section 5 we illustrate the validity and applicability of our emulators on a variety of likelihoods mentioned above. We conclude in Section 6.

METHOD
The general methodology used to construct training and testing sets for our neural network emulators follows that of Spurio Mancini et al. (2022). Here we provide a brief summary and refer the reader to that paper for further details. To compute all of the cosmological quantities on which our emulators are trained and tested we use v2.9.4. 2 We proceed in six steps, as follows.
(i) We specify the cosmological parameters upon which the emulators are built, i.e., the inputs for our neural networks. We choose the following six ΛCDM parameters, as defined in : • , the physical baryon density; • cdm , the physical dark matter density; • 0 , the Hubble parameter today; • , the reionization optical depth; • s , the scalar spectral index; • ln 10 10 s ; where s is the amplitude of primordial scalar perturbations.
In extensions, we also add , eff , or Σ , defined in the previous section.
Note that we opt for 0 rather than for the angular size of the sound horizon at recombination, 100 (which we obtain as a derived parameter), due to the slightly different definitions of the angular acoustic scale implemented in and 3 . For MCMC analyses, if one needs to sample over 100 rather than 0 , it is straightforward to do so, for example using a shooting method with the derived parameter emulators 4 .
(ii) We specify the range of the parameter values over which we train our emulators, reported in Table 1. These ranges should be conservative enough for all physically relevant analyses, assuming the use of data with similar or better constraining power than Planck for the CMB, and DES or BOSS for galaxies. We emphasize that the emulators should not be used outside the ranges of Table 1 or for different cosmological models, as there is no guarantee on their accuracy in such cases. For the ΛCDM+ eff model, we let eff vary between 1.5 and 5.5. In practice, we vary the parameter ur in between 0.49 and 4.49 and obtain eff as a derived parameter 5 . In order to generate the data necessary for the matter power spectrum 2 We initially attempted to use the latest version, i.e., v3; however, it appears that the numerical accuracy of the code has changed compared to v2.9.4 and it was not clear how to recover sufficiently high accuracy for our purposes at high ℓ (see https://github.com/lesgourg/class_public/issues/494). Thus, here we opted for v2.9.4. We note that the updates made in v3 compared to v2.9.4 are detailed on the repository webpage. 3 In , is defined at the maximum of the visibility function, while in it is defined when the optical depth equals unity, and these times are slightly different. The of our emulator is the same as that in . 4 This adds an extra 20 ms to each MCMC step and would therefore increase slightly the convergence time of the chains. Note that the shooting method is also the standard way to sample over 100 in . 5 ur = eff − ( ncdm /(4/11) 1/3 ) 4 , where ncdm is in units of cmb today. This is the convention used by .
at different redshifts, we add pk , the redshift at which we output the matter power spectra, as an extra varying input parameter to our grid. (iii) We generate a Latin hypercube (LHC) of the parameter space. We set S = 128, 000 as the number of samples in the LHC 6 , which is the same as the number of spectra we aim at computing. To generate the LHC, we use the python library pyDOE 7 .
(iv) We use a computing cluster to calculate the S training samples in parallel. Note that despite the total number of training samples being similar to that needed for an MCMC, the task of generating training samples is embarrassingly parallel, since the computation for a given set of input parameters does not depend on any otherunlike an MCMC, which is a path-dependent calculation in parameter space. Given this, we find that computing resources are more optimally used by running each computation on a single thread (i.e., setting OMP_NUM_THREADS = 1 before running ). The total amount of disk space taken by the training and testing samples is ≈ 150GB. Note that for each sample, in addition to the cosmological observables that we want to emulate, we also save 14 derived parameters, such as 8 and 100 s (see the end of this section for details). For each class of emulators (ΛCDM + extensions) the computation of the data requires ≈ 48 hours on a modern high-performance computing cluster. (The accuracy settings for the Boltzmann codes are given in Section 3.) (v) We process the data so that the training samples can be used in the C P training pipeline. For all quantities (TT, TE, EE, PP, PKL, PKNL, DA, H, S8, DER) we select 80% of the samples for training and leave the remaining ones for testing. The TT, EE, PP, PKL, PKNL, H, and DER data consists of positive numbers with large dynamic ranges; to ease training, we take the logarithm of all these quantities before passing them to C P . The TE power spectra are oscillatory, zero-crossing functions, and thus we cannot directly take the logarithm. In this case, we use Principal Component Analysis (PCA), which reduces the dimensionality of the dataset and its dynamic range. We retain only 64 principal components for each spectrum. We verify that adding more does not lead to any significant improvement (Spurio Mancini et al. 2022).
(vi) We generate the emulators using the C P functions and verify their accuracy on the test set. Generating each emulator involves training a neural network model. This operation takes O (1 hr) for each emulator and is best performed on a GPU, to fully take advantage of the acceleration provided by the T F library, which C P uses for its neural network implementations. We also stress that for a given cosmological model once the training is done and the emulator generated, this step does not need to be repeated. Note that we emulate TT, TE, EE, PP, PKL, PKNL, DA, H, S8, and DER separately. For all emulators we use dense neural networks with 4 hidden layers, each with 512 nodes. When building the emulators for TT, TE, EE, PP, DA, H, S8, and DER we do not include pk in the mapping from parameters to data, since pk has no effect on these quantities. Similarly we also remove from the parameter sets for the PP, DA, H, S8, PKL, and PKNL emulators.
The trained TT, TE, EE, PP, PKL, PKNL, DA, H, S8, and DER emulators are stored as pickle files. The size of each CMB power spectra emulator file is 25.9MB, except the TE emulator which is lighter because of the PCA and is 6.3MB. The PKL and PKNL emulators are 4.2MB; the H, DA, and S8 emulator files are 13.5MB; and the DER emulator is 3.2MB. (The size of each emulator file 6 This is roughly the same amount of samples as in a standard cosmological MCMC analysis. 7 https://pythonhosted.org/pyDOE/ is proportional to the number of data points that are saved: 11000 multipoles for the CMB TT and EE spectra, the number of PCA weights for the TE spectra, 500 wavenumbers for PKL and PKNL, and 5000 redshifts for DA, H, and S8.) For illustration, examples of CMB power spectra are shown in Figure 1; examples of matter power spectra calculations are shown in Figure 2; and examples of calculations of ( ), ( ), and 8 ( ), which are used in the BAO and RSD calculations, are shown in Figure 3.

Matter power spectrum
To construct the emulators for the linear power spectrum and its non-linear corrections, we add an extra input to the neural networks, namely, the redshift at which the power spectra are computed, pk . This represents an additional parameter, sampled between 0 and 5. For instance, the LHC used for our ΛCDM emulator has seven dimensions: the six ΛCDM parameters (cf. Table 1) augmented by pk . For each LHC sample we save the matter power spectra at pk on a logarithmically spaced -grid between min = 10 −4 Mpc −1 and max = 50 Mpc −1 with 500 points. (See Section 3 for our perturbation and non-linear settings.)

BAO and redshift-space distortions
Recent galaxy survey data allow us to constrain models based on various BAO distance measures, as well as using redshift-space distortions to constrain the parameter combination 8 , where is the growth rate of cosmic structures. To compute BAO distances and 8 , we need to save the angular diameter distance , the Hubble parameter , and 8 as a function of , as well as the comoving sound horizon at baryon drag , which is last in the list of derived parameters discussed above. We save ( ), ( ), and 8 ( ) on a linearly spaced -grid between min = 0 and max = 20 with 5000 points. The upper bound max = 20 is much higher than what is relevant to galaxy surveys; however we record the high-distances as they may be useful to other applications, such as studies of reionization or cosmic-dawn 21cm measurements. With these quantities, we compute BAO distances and 8 straightforwardly (see, e.g., Alam et al. 2017). Note that 8 = −(1 + ) 8 ( )/ where the derivative can be evaluated numerically. This is how we compute 8 , as is done in .

Derived parameters
We save 14 commonly used derived parameters, namely: • the angular size of the sound horizon at decoupling 100 s , • the amplitude of matter clustering 8 , • the primordial Helium fraction P , • the reionization redshift reio , • the number of effective relativistic degrees of freedom in the early Universe eff , • the conformal time at which the visibility function reaches its maximum (i.e., the recombination time) rec , • its associated redshift rec , • the comoving sound horizon at recombination s,rec , • the conformal angular diameter distance to recombination a,rec , • the conformal time at which the photon optical depth crosses unity ★ , • the redshift ★ at which the photon optical depth crosses unity, Our emulators should never be used outside these ranges. For the matter power spectrum emulators the LHC is supplemented by pk , the redshift at which the power spectrum is evaluated, which is varied between 0 and 5. For , the values in parentheses refer to the prior bounds that are used in the emulators for extensions, which are slightly broader than for the ΛCDM emulators (outside the parentheses).
• its associated comoving sound horizon s,★ , • and conformal angular diameter distance a,★ • and the comoving sound horizon at baryon drag .
These parameters are saved into a list for each sample in the LHC, with the ordering given above.

BOLTZMANN CODE ACCURACY AND SETTINGS
To ensure high precision in our calculations, we carefully study and calibrate the numerical precision parameters in and . 8 For the maximum multipole, we choose ℓ max = 11000 in both codes. We chose ℓ max = 11000 as some forecasts and calculations, e.g., in the context of CMB-HD (Sehgal et al. 2019), require CMB power spectra at such high multipoles. For non-linear corrections to the matter power spectrum, we use HM (Mead et al. 2021) with fiducial values min = 3.13 and 0 = 0.603 in both codes. These values correspond to a dark matter-only scenario and are the same as in the baseline Planck 2018 analyses (Planck Collaboration 2020a), justifying our choice. For reionization modelling we kept the default setting (i.e., the same as in ). We remark here that in future work it will be useful to train emulators in which these HM parameters are also varied, or for other non-linear regime prescriptions. Indeed, as shown by, e.g., McCarthy et al. (2022), future CMB power spectrum data will be sensitive to the impact of baryonic physics on the matter power spectrum (via CMB lensing); in addition, LSS statistics are sensitive to the nonlinear and baryonic prescriptions at high values of . With emulators that include these parameters, we will be able to marginalize over this uncertainty.
The TT, TE, EE, and PP spectra computed with 9 are fully converged when using the following parameters (see, e.g., McCarthy et al. 2022;Hill et al. 2022): set_matter_power: We use the spectra computed with these settings as a highaccuracy reference. This high-accuracy calculation takes O (1 min) per sample on 16 threads.
It is possible to match TT, TE, EE, TE, PKL, and PKNL spectra to to better than 0.1% precision (for ℓ < 3000 and < 10 Mpc −1 ), by using the precision parameters in the cl_ref.pre file in the GitHub repository (this file ensures convergence at the 0.01% level internally to for TT and EE) (Lesgourgues 2011b) and setting k_max_tau0_over_l_max: 15 to ensure high accuracy at ℓ > 4000 as discussed below). This is illustrated in Figures 4 and 6, where this calculation is shown as the dotted-dashed yellow lines, labeled "class prec". We refer to it as the ultra-high-precision prediction. Nonetheless, this computation takes O (1 h) which is prohibitively long to generate O (10 5 ) spectra for training the emulators.
By optimizing the tradeoff between precision and computation time, we find that the minimal settings required in order to match the high-precision calculation with v2.9.4 are obtained by setting the following three parameters: 10 • accurate_lensing: 1 • k_max_tau0_over_l_max: 15.
• perturb_sampling_stepsize: 0.05 These are the precision parameters that we use in order to generate our emulator training data. The first parameter ensures that the lensed TT, TE, and EE spectra are converged for ℓ > 3000. Without it, they have non-physical oscillatory features. The second parameter ensures convergence at high-ℓ for all the spectra, including the unlensed ones. The third parameter is critical to get high-accuracy spectra over the whole multipole range. It is particularly important to get a converged PP spectra for ℓ < l_switch_limber (default value in v2.9.4: l_switch_limber=10). Note that out of these three parameters, only the second one has an impact on the linear and non-linear matter power spectra. This high-precision calculation takes roughly the same amount of time as the one, i.e., O (1 min). 9 We use v1.3.6. 10 Note that these parameters are optimized here, whereas McCarthy et al.
(2022) used extremely high-precision parameters without optimizing the tradeoff between precision and computation time.

EMULATORS
Our trained TT, TE, EE, PP, PKL, PKNL, H, DA, S8, and DER emulators in ΛCDM, CDM, ΛCDM+ eff , ΛCDM+Σ are made publicly available on the C P GitHub repository. 11 To predict the power spectra, they need to be imported in Python, via the C P package. They allow for rapid computation of power spectra and derived parameters, and hence of parameter posterior probability distributions in MCMC analyses (see Section 5). Loading the emulators takes ≈ 0.1 sec. Computing one set of TT, TE, EE, PKL, PKNL power spectra and derived parameters takes ≈ 60 ms, compared to O (1 min) with class or camb with high-precision settings, i.e., we achieve a factor of 1000 speed-up. The GitHub repository also contains a set of notebooks showing how to use the emulators.
In Figure 1-6 we show the emulator's predictions for one set of cosmological parameters, namely, the central values of the right column of Table I in Planck Collaboration (2020a), with one massive neutrino with = 0.06 eV and in ΛCDM. In Figure 1, we compare the power spectra from , , and C P explicitly, The (i.e., the high-precision calculation discussed in Sec. 3), "class prec" (i.e., the ultra-high-precision calculation discussed in Sec. 3), high-precision , and C -P power spectra are indistinguishable for ℓ < 10000. For 10000 < ℓ < 11000, some differences can be seen in the TT and EE power spectra: the and C P predictions tend to fall off compared to . Since this is less significant for the "class prec" prediction, we attribute this to numerical precision settings, and do not investigate further as this is probably irrelevant for any CMB experiment in the foreseeable future. The shaded areas indicate the sensitivity of Advanced ACT (Henderson et al. 2016), SO (Simons Observatory 2019), and CMB-S4 (CMBS4 2022) forecast sensitivity. We refer to Section 4 of Bolliet et al. (2022) for details on the TT and PP sensitivity curves and the links where they can be downloaded from. The TE sensitivity can be computed from TT and EE 12 (see, e.g., Spurio Mancini et al. 2022, for the formula we use). Note that we use sky = 0.3 for Advanced ACT and sky = 0.4 for SO and CMB-S4. Also, note that all noise curves also assume the use of Planck data in combination with that from each ground-based experiment (see Sec. 2 of Simons Observatory (2019) for details); the use of Planck is important in extending the forecast measurements to low multipoles, where otherwise atmospheric noise in the ground-based data would be very large. Figure 2 shows the same as Figure 1 but for the linear and non-linear matter power spectra, computed at three redshifts ( = 0, 2.5, 5). In these cases, the , , and C P predictions are indistinguishable across the full range. The same is true for the H, DA, and S8 predictions, as can be seen in Figure 3.
In Figure 4, we show the relative difference between CMB angular power spectra (TT,TE,EE,PP) and , in percent. The spikes in TE are where the spectra cross zero and are not problematic. As discussed before, the agreement between the ultra-high precision prediction and is at the 0.1% level until ℓ ∼ 3000, and degrades at higher ℓ. Nonetheless, the agreement between C P and both predictions remains at the 0.1% level across the whole range of multipoles, including ℓ > 3000, except in two cases. First, 11 https://github.com/cosmopower-organization 12 The CMB-S4 EE noise is from S4_190604d_2LAT _pol_default_noisecurves_deproj0_SENS0 _mask_16000_ell_EE_BB.txt and the SO EE noise is from v3.1.0/SO_LAT_Nell_P_baseline_fsky0p4_ILC_CMB_E.txt in the low-ℓ part (ℓ < 10) of the lensing convergence power spectrum we see very small differences ( 0.5%) between the ultra-high precision prediction on one hand and the high-precision and C -P prediction on the other hand (the agreement between the C P and high-precision calculation remains better than 0.005%). This difference can be reduced further by decreasing the parameter perturb_sampling_stepsize (see Sec 3), at the cost of a longer runtime. Second, in the low-ℓ part (ℓ < 50) of the EE power spectrum we see differences between the C P and high-precision prediction at the 0.6% level. If needed, this could be solved by generating more training data.
To assess whether this level of difference is acceptable for CMB Stage-IV analyses, we compare the relative difference between the CMB spectra in terms of CMB-S4 sensitivity in Figure 5. We see that the power spectra agree with to better than 0.03 across all multipoles and with the high-precision prediction to better than 0.03 . Therefore, such level of agreement is sufficient.
In Figure 6 we show the relative difference between the linear and non-linear matter power spectra with respect to in percent, at three redshifts, = 0, 2.5, and 5. In the left panels, we see that the C P linear matter power spectra agree with the highprecision predictions at the 0.05% level across the whole range. The difference between these and the ultra-high precision prediction is at the 0.3% level. Finally, the difference between the ultra-high precision prediction and is at the 0.3% level up to ≈ 20 Mpc −1 and then becomes ≈ 0.8% at higher . In the right panels, showing the non-linear matter power spectra, we see that differences between on one hand and and on the other hand are more important for the non-linear matter power spectra, becoming larger in the non-linear regime, but remain at the 1-1.5% level, including at = 5.
Since there is no direct measurement of the matter power spectrum we cannot compare it with an instrumental sensitivity level. However, we can check the accuracy on observables based on the matter power spectrum. The lensing convergence power spectrum (PP, bottom right panel of Figure 1) is one such example. Other examples include galaxy clustering and galaxy weak lensing observables and to a lesser extent, the lensing of CMB angular power spectra. We perform Stage III posterior inference analyses based on these observables and recover nearly identical constraints to those obtained running the full Boltzmann code calculations (see Section 5), giving us confidence that our emulators are sufficiently precise. We defer to future work a more detailed quantitative analysis, including the release of high-accuracy emulators for the matter power spectrum tested against (forecast) Stage-IV configurations.
The results shown in Figures 1-6, are for one set of cosmological parameters. However, it is important to quantify the precision of our emulators over the whole prior range of parameters (see Table  1), in all the cosmological models. We do so in Figures 7-13 13 . These figures show the relative difference between the emulators, i.e., the C P prediction, and the exact high-precision prediction on the testing data sets 14 , containing ≈ 25000 predictions. The results for TT are in Figure 7, for EE in Figure 8, for TE in Figure 9, for PP in Figure 10, for PKL and PKNL in Figure 11, for H, DA, and S8 in Figure 12, and for the derived parameters in Figure 13. For TT, TE, EE, and PP we show the difference in 13 The "mnu" results on this figure are for one massive and two massless neutrino states. 14 The testing data is made of the exact computations which are not used for training the emulator neural networks. We use 20% of our dataset for testing. (for default precision as the solid red line and ultra-high precision settings as the dotted yellow line labeled "class-prec"), C P , and . We show the dimensionless ℓ (ℓ + 1) ℓ /2 for TT, TE and EE, and [ℓ (ℓ + 1) ] 2 ℓ /2 for PP (lensing convergence power spectrum) in a ΛCDM model with one set of cosmological parameters from the Planck 2018 results (see text for details). The shaded areas indicate the forecast 1 uncertainty for Advanced ACT, SO, and CMB-S4 (see Section 4 for details). This figure illustrates the emulated quantities, the relative difference between the spectra are shown in Figure 4 (in percent) and 5 (in units of CMB-S4 statistical error bars). On the EE plot, the thin black solid lines indicate the cosmic variance.  Figure 1 but for matter power spectra. We show the linear matter power spectra on the left and the non-linear spectra on the right (computed according to HMcode with min = 3.13 and 0 = 0.603, within ). In both panels, we show the power at three redshifts, = 0, 2.5, 5 from top to bottom, spanning the range of the emulators. This figure illustrates the emulated quantities, the relative difference between these spectra are shown in percentages in Figure 6. units of CMB-S4 sensitivity and in percentages, and for the other emulators we show the difference in percentages (since there is no simple way to compare with instrumental sensitivity; see comments in the previous paragraph). In all the figures, the three shades of red correspond to 68%, 95%, and 99% of the testing data set, from dark to light, respectively. For instance, for TT ( Figure 7) we find that our emulators predict the power spectra to better than 0.02 at CMB-S4 precision level, 68% of the time in all the four cosmological models considered here. The agreement for TE and EE is similar. For PP the agreement is at the 0.02 CMB-S4 precision level, 99% of the time in all the four cosmological models considered here.
The quantity that is the least accurately emulated is the non-linear matter power spectrum, which is reproduced at the 1 − 1.5% level 95% of the time in all models (see Figure 11). This level of agreement is satisfactory for application to current survey configurations, especially since non-linear modeling and baryonic uncertainty is at the ≈ 10% level (see, e.g., Mead et al. 2021) and even exceeds 30% at > 1 Mpc −1 (see, e.g., Amon & Efstathiou 2022).
The emulated redshift-dependent ( ), ( ), and 8 ( ) agree to better than 0.05% between = 0 and 20, for 99% of the testing set in all models (see Figure 12). We find a similar agreement between emulated derived parameters and the testing dataset (see Figure 13). , angular diameter distance (middle), and 8 (right) as computed with and C P , between = 0 and 20 (i.e., spanning the redshift range of our emulators). We use the same settings as Figure 1 and note that in this parameters configuration the relative difference between and C P (not shown here) is less than 0.2%. This figure illustrates the emulated quantities, see Figure 4 for the relative difference over the full testing set.  settings of (high precision as the solid red line and ultra-high precision settings as the dotted yellow line labeled "class-prec"), and C P with respect to in percent. We use the same parameters settings as in Figure 1. Note that the non-negligible fractional deviations seen in a few places for TE occur due to the zero-crossing of the spectrum at those multipoles, and do not reflect inaccuracy in the emulator.

ACCELERATED LIKELIHOOD ANALYSES
In this section we run posterior inference analysis and use the emulators presented above in MCMC extractions of cosmological parameters. A full Stage-IV analysis will require optmization of likelihood codes alongside emulators which is beyond the scope of this paper and will be covered in future work. Therefore, here we demonstrate accelerated likelihood analyses for current (Stage II-III) data and reproduce marginalized constraints for ACT DR4 (Subsection 5.1), Planck lensing + DES + BAO (Subsection 5.2) and Planck TT,TE,EE in CDM, ΛCDM+ eff and ΛCDM+Σ (Subsection 5.4). We use (Torrado & Lewis 2021) for MCMC sampling. To analyse the resulting chains we use GetDist (Lewis 2019).
We define a C P theory object within such that we are able to implement the C P call completely outside of the likelihoods. Our implementation of the C P theory object is available upon request. In a forthcoming paper, we plan to release an official C P wrapper for both and (Zuntz et al. 2015).

Accelerated ACT DR4 power spectrum analysis
We use the ACT Data Release 4 (DR4) "actpol_lite" likelihood (Choi et al. 2020;Aiola et al. 2020) with C P to reproduce the original DR4 results. This is the foreground-marginalized likelihood publicly available in the pyactlike repository 15 . For the C P runs, this likelihood involves the TT, TE, and EE emulators, relying on multipoles up to ℓ ≈ 4500. In Figure 14 we show the resulting 2D marginalized posterior probability distribution obtained in a few minutes on a laptop with C -P as red empty contours, and compare it with the publicly available chains 16 shown as the filled green contours. The reference chains were obtained in Hill et al. (2022) using high-accuracy settings and took several days to converge. The overlap of the C -P and reference contours is nearly perfect. We note that for the C P run we sample over 0 and obtain as a derived parameter using our DER emulator in post-processing of the chains, whereas the reference chains were computed with as an input parameter and 0 was saved as a derived parameter.
As a second test, we perform a maximum-likelihood analysis using the full ACT DR4 likelihood, including foregrounds. To do so we use a public python implementation bpllike 17 . The C P results (blue contours) are shown in Figure 15 and compared with the reference high-precision chains that we ran using the original full Fortran ACT DR4 likelihood 18 (red contours). Again, the overlap between the contours is nearly perfect. The foreground parameter contours (for the thermal SZ and kinematic SZ amplitudes and CIB SED power-law index and amplitude) match exactly. Cosmological parameter contours agree to better than 0.1 , except for . This is simply because of different definitions in and . Note that our chains reach − 1 = 0.05 in ≈20 minutes on a laptop with ten threads.

Accelerated Planck lensing + DES + BAO analysis
To further demonstrate the efficiency and wide range of our emulators, we reproduce cosmological parameters extracted from the 16 CLASS2p8_ACTPol_lite_DR4_leakfix_yp2_baseLCDM_taup_hip_R0p01 downloaded from LAMBDA. 17 https://github.com/ACTCollaboration/bplike 18 ACT DR4 Fortran likelihood, named actpolfull there.
Planck CMB lensing power spectrum + DES + BAO analysis. We compare with the publicly available Planck chains 19 . The chains of interest are labeled DES_lenspriors_lensing_BAO.
We use the Python-native re-implementation of the Planck reconstructed lensing power spectrum likelihood, assuming their fiducial (conservative) multipole range. This likelihood is available in as planck_lensing_2018 (Planck Collaboration 2020b). We apply the priors used in the Planck lensing analysis: Gaussian priors on s = 0.96 ± 0.02 and Ω b ℎ 2 = 0.0222 ± 0.0005, a flat prior 40 < 0 /(km/s/Mpc) < 100, and a fixed = 0.055 (see Planck Collaboration 2020b, for details).
Finally, we use the BAO likelihoods corresponding to BOSS DR12 (Alam et al. 2017), the SDSS Main Galaxy Sample (Ross et al. 2015), and the 6dF survey (Beutler et al. 2011). These are implemented in as sdss_dr12_consensus_bao, sdss_dr7_mgs, and sixdf_2011_bao.
For the C P runs, this likelihood combination involves the PKNL, H, DA, and PP emulators. Because this calculation requires the computation of the non-linear ( ) at many redshifts, it is slightly more time-consuming. Indeed, we need to evaluate the matter power spectrum on a 2D grid in ( , ), hence requiring us to call the emulator as many times as there are points in the dimension. We set the ( , ) grid to cover redshifts between 0 and 4 with 15 linearly spaced points. With this choice, performs ≈16 evaluations per second on an ARM64 MacBook Pro. We reach − 1 0.15 after ≈30 minutes on a laptop, with four chains taking a total of 56,000 accepted steps with acceptance rate of 0.2.
The results are shown in Figure 16. The overlap between the C - . Relative difference of the linear (left) and non-linear (right) matter power spectra between two different settings of (high precision as the solid red line and ultra-high precision settings as the dotted-dashed yellow line labeled "class-prec") and C P with respect to in percent. P and reference contours is perfect. As derived parameters, we show 8 ≡ 8 (Ω m /0.3) 0.5 and 0.25 8 ≡ 8 (Ω m /0.3) 0.25 , which are most constrained by galaxy weak lensing and CMB lensing data, respectively. Our C P runs recover the reference constraints on 8 and 0.25 8 exactly. This validates the matter power spectrum emulators for current galaxy WL surveys.

RSD Full shape analysis
To demonstrate that we can also compute 8 and use it in a likelihood analysis, using our S8 emulator (see Section 2), we perform a likelihood analysis and use sdss_dr12_consensus_fullshape with Planck lensing. We compare and C P chains in Figure 17.

Extensions
To validate our emulators in extended cosmology we run maximumlikelihood analyses in CDM, ΛCDM+ eff , and ΛCDM+Σ . For the data and likelihood we choose the official clik Planck 2018 plikHMTTTEEE+lowl+lowE (Planck Collaboration 2020a) implementation, which we call through . We show the CDM, ΛCDM+ eff , and ΛCDM+Σ constraints in Figure 18, 19, and 20, respectively. In each extension, we overplot the reference contours from chains (empty red contours) that we obtain online (see footnote 19). Our C P contours are in blue. In all cases the agreement with the reference Planck chains is excellent. To reach − 1 ≈ 0.1, C P needs approximately 25 minutes, while or chains typically take O (1 day) to converge.

CONCLUSIONS
We have extended the work of Spurio Mancini et al. (2022) by creating emulators for CMB temperature and polarization angular anisotropy power spectra, CMB lensing convergence power spectra, linear and non-linear matter power spectra, and the redshift evolution of ( ), ( ) and 8 ( ), in ΛCDM and extensions, namely, CDM, ΛCDM+ eff and ΛCDM+Σ with one or three massive states (see Section 2). All of these quantities are computed using high-precision settings (see Section 4). They are sufficiently precise for Stage-IV CMB analyses (see Section 3) and current galaxy weak lensing and clustering surveys. We have tested all of our emulators in maximum-likelihood analyses, involving high-ℓ CMB, galaxy weak lensing and clustering, BAO and RSD data (see Section 5).
The main outcome of this work is to open the door to fast param-    Figure 8. Relative difference between C P and the "true" high-precision prediction for CMB EE power spectra in percent (top) and in units of CMB-S4 sensitivity (bottom). See end of Section 4 for details. eter inference, via the widely used MCMC method, with accelerated Boltzmann computations thanks to our C P emulators. For instance, as we demonstrate here, most of Stage-III parameter inference is now feasible on a laptop.
In a forthcoming paper we plan release wrappers of our emulators so that they can be readily used within (Torrado & Lewis 2021) and (Zuntz et al. 2015). We created an online repository to store our emulators and will continue updating it. For CMB polarization, we have made emulators for E modes, which could be used to accelerate LiteBIRD (LiteBIRD Collaboration 2022) analyses dedicated on reionization history constraints (Zaldarriaga et al. 2008;Allys et al. 2022). We defer B-modes power spectra emulators to future work.
In addition to the applications considered here, these emulators can also be used to save time at the initial step of libraries for computations of cosmological LSS observables like _ (Bolliet et al.   Figure 9. Relative difference between C P and the "true" high-precision prediction for CMB TE power spectra in percent (top) and in units of CMB-S4 sensitivity (bottom). Note that the spikes are simply due to zero crossings of the TE power spectrum. See end of Section 4 for details.  Figure 10. Relative difference between C P and the "true" high-precision prediction for CMB lensing convergence power spectra in percent (top) and in units of CMB-S4 sensitivity (bottom). See end of Section 4 for details.
Recent works Cabass et al. 2022) have used the effective field theory of large-scale structure (Baumann et al. 2012;Carrasco et al. 2012) in order to derive constraints from spectroscopic surveys, based on the one-loop galaxy power spectrum and bispectrum. These calculations are time-consuming (MCMC convergence takes O(1day) using _ (Chudaykin et al. 2020). The bottleneck is the calculation of higher-order correlators. The work presented here does not directly help to accelerate such analyses, but nonetheless the method can directly be applied in this context as well.

AUTHOR CONTRIBUTIONS
BB led the work. ASM contributed to the development of the project and software pipeline, helped define the library of models and parameters to be studied, and contributed to the writing and editing of the manuscript. JCH contributed to studies of numerical precision settings, comparisons of and , contributed to the writing and editing of the manuscript, and coordinated the initiation and overall direction of this work. MM aided in investigations into precision settings and prior ranges required for various parameters, helped define the library of models to be produced, developed the bplike software for the ACT DR4 likelihood analysis and contributed to the writing and editing of the manuscript. HTJ is developing the wrappers that will import these emulators into future ACT and SO analyses, and reviewed the paper at the final stages and provided feedback. EC developed the ACT DR4 likelihood, contributed to discussions and to the writing and editing of the manuscript. JD coordinated the initiation of this work, and contributed to the editing of the manuscript.

DATA AVAILABILITY
The emulators are available online at CosmoPower Organisation, including tutorial notebooks. If you use the emulators, please cite our work as well as Spurio Mancini et al. (2022).   , et al., 2008Zuntz J., et al., 2015 This paper has been typeset from a T E X/L A T E X file prepared by the author. . The yellow contours are the constraints using the ΛCDM+Σ emulator with one massive and two massless neutrinos (unlike the Planck chains), rather than three massive neutrinos (blue and red contours). The model with three degenerate states is a better approximation to the mass splitting results. (We thank Antony Lewis for pointing this out to us.) See Subsection 5.4 for further details.