A field-level emulator for modeling baryonic effects across hydrodynamic simulations

Weak gravitational lensing is a powerful tool for measuring the clustering of matter in our universe, and thus obtaining information about matter content and initial conditions of our universe (Kilbinger, 2015) through various summary statistics. However, deriving precise cosmological constraints from weak lensing observations necessitates highly accurate theoretical models that account for baryonic physics, which redistributes matter on small scales via processes such as Active Galactic Nuclei (AGN) feedback. These processes remain poorly understood and inadequately constrained by current observations, leading to challenges in formulating a predictive theory.

While causality ensures that baryonic effects are negligible on large scales (Lewis & Challinor, 2011), baryonic effects gain significance on smaller scales, affecting structure formation through hydrodynamic processes. These processes, including AGN and stellar feedback, can heat gas and inject large amounts of energy into galaxies and the surrounding halo. Specifically, AGN feedback can eject gas to very large distances, which can further modify the dark matter distribution through gravitational interactions.

Studies have shown that probing the small scales contains a wealth of information, leading to stronger parameter constraints (Lu & Haiman, 2021), while ignoring small scale information leads to significant deterioration of these constraints (Köhlinger et al., 2016, 2017; Hikage et al., 2019). This small-scale information is heavily influenced by baryonic effects, and the large uncertainty associated with these effects makes them one of the primary sources of systematic error in weak lensing analyses. Hence, accurate modeling of baryonic effects is crucial when probing the information-rich small scales for an unbiased cosmological analysis.

At present, numerical simulations remain the sole comprehensive method for precise simulation of baryonic effects and the deeply non-linear evolution of cosmic structures. Hydrodynamic simulations provide a detailed and accurate representation of the behavior of baryonic matter by modeling complex physical processes such as gas dynamics, star formation, and feedback mechanisms (Teyssier, 2002; Di Matteo et al., 2005; Jenkins et al., 1998). However, these simulations are not ab initio parameter free, but instead must parametrize the lack of physics understanding via free parameters that can be varied. Furthermore, hydrodynamic simulations require substantial computational resources as compared to dark matter only N-body simulations.

The development of emulators has emerged as a powerful technique to overcome this computational challenge, enabling rapid and accurate predictions of physical properties without the need for running costly simulations. These emulators interpolate simulation results and have been shown to be remarkably accurate. Various emulators have been developed for cosmology, catering to various observables, encompassing the matter power spectrum (Heitmann et al., 2014; Euclid Collaboration et al., 2019; Winther et al., 2019; Angulo et al., 2021; Knabenhans et al., 2021), mass function (McClintock et al., 2019; Bocquet et al., 2020), and the galaxy correlation function and Lyman-$\alpha$ Forest (Zhai et al., 2019; Bird et al., 2019).

While most previous work has focused on modeling the baryonic effects on the matter power spectrum (Huang et al., 2019; Mead et al., 2021a; Aricò et al., 2021a; Schneider et al., 2019, 2020; Giri & Schneider, 2023), there is an increasing need for developing fast baryon models at the field level for analysis beyond two-point statistics. For example, simulation-based inference methods (Cranmer et al., 2020) show great promise in extracting rich non-Gaussian information either through high-order statistics (e.g., Hahn et al., 2023), or directly from the fields (e.g., Dai & Seljak, 2021, 2022; Villaescusa-Navarro et al., 2021a, b). These approaches rely on fast and accurate cosmological predictions from numerical simulations. Previous field-level baryon models, such as Baryon Correction Model (Schneider & Teyssier, 2015) and Enthalpy Gradient Descent (Dai et al., 2018), move the dark matter particles from N-body simulations to mimic the baryonic effects. While they have been shown to accurately predict the power spectrum from hydrodynamical simulations (Schneider et al., 2019), they can be computationally expensive when the particle resolution is high.

By analyzing a diverse range of baryonic feedback hydrodynamics simulations across multiple redshifts, we will show that adding baryons to N-body simulations can be achieved using a field-level transfer function to augment N-body fields with a Fourier mode amplitude, k, dependent transfer function correction. We develop a transfer function emulator using Gaussian process for modeling the baryonic effects in terms of $P_{\rm hydro}(k)/P_{\rm nbody}(k)$, where $P_{\rm hydro}(k)$ is the total matter power spectrum, and $P_{\rm nbody}(k)$ is the dark matter power spectrum.

In this paper, we develop the emulator and show that it is accurate at a percent level over our whole parameter space, which covers scales $0.01\leq k\leq 10$ h/Mpc and redshifts $0\leq z\leq 1.5$. We validate the performance of our emulator against thousands of hydrodynamical simulations and their respective gravity-only counterparts. In particular, we make use of CAMELS-Astrid (Ni et al., 2023), CAMELS-IllustrisTNG, CAMELS-SIMBA(Villaescusa-Navarro et al., 2021c, 2023), BAHAMAS (McCarthy et al., 2017, 2018), Horizon AGN (Dubois et al., 2014), Owls (Schaye et al., 2010; van Daalen et al., 2011), and Eagle (Schaye et al., 2015; Crain et al., 2015; McAlpine et al., 2016; Hellwing et al., 2016) simulations. We also compare against commonly utilized emulators like BACCO (Aricò et al., 2021a), HMcode (Mead et al., 2021a, b), and BCemu (Giri & Schneider, 2023) on all simulations. Finally, we show the improvement of the field-level baryon model against hydrodynamic fields at varying redshifts. This emulator is fast as it only requires a single FFT and its inverse, which enables large-volume N-body simulations for generating realistic weak lensing mock data for cosmological analysis at the field level.

This paper is organized as follows: in section 2, we describe the suite of simulations that are used for training and testing our emulator. In section 3 we explain how our emulator can be used to emulate baryonic effects at the field level. In section 4 we describe the methods and construction of the Gaussian process emulator. In section 5 we test the emulator’s robustness on multiple hydrodynamic test simulations, compare it with currently available emulators, and show field-level improvements using our emulator. We summarize and conclude in section 6.

In this section, we describe the simulations that we employ throughout this paper. Our main suites of simulations are part of the Cosmology and Astrophysics with MachinE Learning Simulations (CAMELS) (Villaescusa-Navarro et al., 2021c, 2023; Ni et al., 2023). CAMELS is a suite of 10,421 cosmological simulations each with a comoving volume of $(25~{}h^{-1}\text{Mpc})^{3}$ evolved from $z=127$ to $z=0$ with $256^{3}$ dark matter particles and $256^{3}$ gas particles in the initial conditions. These contain 5,097 N-body simulations and 5,324 hydrodynamic simulations. Notably, each hydrodynamic simulation in CAMELS pairs with an N-body counterpart, sharing identical cosmological parameters and initial random seeds.

Simulations in CAMELS are categorized into various suites (Astrid, IllustrisTNG, and SIMBA) and sets based on the employed code for running the simulations and the arrangement of cosmological and astrophysical parameters ($\Omega_{m},\sigma_{8},A_{\rm SN1},A_{\rm AGN},A_{\rm SN2},A_{\rm AGN2}$), as well as the initial random seeds. The Astrid suite comprises 1,092 hydrodynamic simulations executed using the MP-Gadget simulation code (Feng et al., 2018), employing analogous subgrid physics as the original Astrid simulations (Ni et al., 2022; Bird et al., 2022). Additionally, the IllustrisTNG suite (based on Vogelsberger et al., 2013; Torrey et al., 2014) and the SIMBA suite (based on Davé et al., 2016), with 1,092 hydrodynamic simulations each, are executed using the AREPO code (Springel, 2010; Weinberger et al., 2020) and the GIZMO code (Hopkins, 2015), respectively.

Each simulation is characterized by its cosmology (given by $\Omega_{m}$ and $\sigma_{8}$) and its astrophysical feedback (given by $A_{\rm SN1},A_{\rm AGN},A_{\rm SN2},A_{\rm AGN2}$). In particular, throughout CAMELS’ suites, the astrophysical parameters represent the value of subgrid physics parameters that influence stellar and Active Galactic Nuclei (AGN) feedback mechanisms.

We made use of the Latin Hypercube, LH, set within each suite of CAMELS, which contains 1,000 simulations whose cosmological and astrophysical parameters are arranged in a latin-hypercube within a very broad range¹¹1We note that in the case of Astrid, the parameter $A_{\rm AGN2}$ varies between 0.25 and 4.

and every simulation has a different value of the initial random seed. Within each LH set, CAMELS provides 1,000 simulations for each redshift that span a wide range of cosmologies and baryonic feedbacks, perfect for our purposes of capturing the underlying physics.

Importantly, each suite has been run with a different code and therefore the subgrid physics model is completely different. So, notably, while the range of variation of the above parameters remains consistent across all CAMELS suites, the precise definitions and overall impact of these astrophysical parameters vary significantly across suites. These simulations are designed to train and set machine learning algorithms given the way their cosmological, astrophysical, and initial random seed parameters are set. We utilize this set for each suite throughout this paper to train and test our methodology on simulations with significantly different cosmologies and astrophysics.

Figure 1 shows the baryonic effect in the matter power spectrum, $P_{\rm hydro}(k)/P_{\rm nbody}(k)$, across different redshift values in all CAMELS suites utilized in this study. From the figure, it is clear that baryonic feedback can have very diverse and strong effects on the matter power spectrum, especially on small scales. SIMBA, with its aggressive AGN feedback, produces the most prominent suppression of the matter power on large scales. On the other hand, IllustrisTNG exhibits a more moderate impact on the matter power spectrum as a consequence of having milder AGN feedback, and Astrid spans the widest range of baryonic feedback, encompassing effects seen in both SIMBA and IllustrisTNG.

To cover the broadest range of baryonic feedback, based on Figure 1, we used simulations from the Astrid suite at $z=0.0$ to train our emulator. For selecting the training simulations, we performed a random sampling of 800 simulations (out of the 1000 available) from the Astrid suite within CAMELS at $z=0.0$.

Post-training, we test the emulator using the remaining 200 simulations from Astrid at $z=0.0$ alongside all other available hydrodynamic simulation suites in CAMELS. Figure 2 illustrates the matter power spectrum ratio in the Astrid $z=0.0$ test data. Additionally, the IllustrisTNG suite and the SIMBA suite are part of the test dataset.

Previous studies (Heitmann et al., 2014; Smith et al., 2014; Rasera et al., 2014; Heitmann et al., 2013) have shown that both high physical resolution and large box sizes are required to guarantee convergence of the power spectrum. Schneider & Teyssier (2015) showed that deviations of the power spectrum ratio using small boxes, like in CAMELS, are at the $5\%$ level. However, these could also be due to cosmic variance which affects the small CAMELS-like box volumes. To study the effects of cosmic variance on the matter power spectrum ratio, we show the ratios for simulations with the same cosmology and astrophysics from CAMELS-Astrid’s CV set in figure 3. From this, we can see that the matter power spectrum ratio is affected by cosmic variance up to $\sim 10\%$ on small scales, suggesting that the deviations for small boxes are due to cosmic variance.

In addition to the diverse array of simulations in CAMELS, we extended the validation of our emulator by testing it against simulations outside the CAMELS database. These external simulations, including BAHAMAS (McCarthy et al., 2017, 2018), Horizon AGN (Dubois et al., 2014), Owls (Schaye et al., 2010; van Daalen et al., 2011), and Eagle (Schaye et al., 2015; Crain et al., 2015; McAlpine et al., 2016; Hellwing et al., 2016), serve as crucial benchmarks to assess the robustness and generalizability of the emulator’s predictions beyond CAMELS. These simulations encompass various diverse physical processes including AGN feedback, supernovae feedback, mass loss from Asymptotic Giant Branch stars, radiative cooling, stellar winds, and stellar initial mass function, among others. Additionally, these simulations differ from those in CAMELS in volume, resolution, and subgrid physics code (van Daalen et al., 2020), enabling us to test the robustness of our emulator to these effects. The solid lines in the left panel of Figure 6 illustrate the matter power spectrum ratio derived from these external simulations, allowing us to assess how the emulator performs across varied simulations beyond the scope of the CAMELS.

In this section, we present our methodology to model baryonic effects at the field level for the total matter density field. We start by computing the cross-correlation coefficient between the matter field in full hydrodynamic simulations and their N-body counterparts. The cross-correlation coefficient, $r(k)$, is defined as:

Here, $P_{\rm cross}(k)$, $P_{\rm nbody}(k)$, and $P_{\rm hydro}(k)$ represent the cross-power spectrum, the N-body power spectrum, and the power spectrum of the hydrodynamic fields respectively. This coefficient’s range spans from $-1$ to $1$, where values closer to $1$ signify a strong positive linear relationship, $-1$ indicates a strong negative linear relationship and $0$ implies no linear relationship between the datasets.

Figure 4 shows the cross-correlation coefficients derived from all CAMELS suites at different redshift values pertinent to this study. These coefficients serve as indicators of the correlation strength between N-body and hydrodynamic fields within the simulations.

We can see that the calculated cross-correlation coefficients are very close to 1, down to $k\sim 10$ h/Mpc, for all the simulations, with most deviations being within $5-10\%$. On the other hand, we see in Figure 1 that baryonic effects can cause deviations of up to $\sim 50\%$ on the matter power spectrum, with the effects getting more dominant at smaller scales. This suggests that the baryonic effects predominantly impact the amplitude of the Fourier modes (given by the power spectrum) rather than their phases (given by cross-correlation coefficients).

Since the amplitude and phase of the Fourier modes completely describe the fields, with baryonic effects mainly changing the amplitudes, aligning the power spectra of N-body fields with their hydrodynamic counterparts would also effectively align them at the field level, facilitating a cost-effective field-level analysis using just N-body simulations.

We can achieve this power spectra alignment by applying a transfer function to the N-body fields (Bond & Szalay, 1983; Seljak & Zaldarriaga, 1996). Transfer functions operate by performing specific modifications to the field data. In Fourier space, each mode of a field is represented by its amplitude, typically denoted by $|\mathbf{k}|$, and its phase. Transfer functions act on $\mathbf{k}$ to alter their amplitudes according to certain criteria (Peacock, 1999; Dodelson, 2003).

The field-level transformation using a transfer function, $T(k)$, is mathematically defined as:

Here, $F(\mathbf{k})$ symbolizes the original field in Fourier space, and $F^{\prime}(\mathbf{k})$ represents the transformed field in Fourier space after the element-wise application of the transfer function $T(k)$ to the original field $F(\mathbf{k})$. This transformation enables the adjustment of simulated fields to match desired characteristics or observational data, enhancing the accuracy or realism of the simulation results.

In our context of incorporating baryonic effects in N-body simulations, we apply a transfer function to the simulated field to adjust its power spectrum. Defining baryonic suppression as:

Our field transformation on N-body fields is then:

This could increase or suppress the power of certain scales and correct discrepancies arising from missing baryonic physics effects in N-body simulations, aligning the power spectra of N-body fields with the full hydrodynamic fields.

To fulfill the promise of field-level modeling of baryonic effects, we need to characterize $P_{\rm hydro}(k)/P_{\rm nbody}(k)$ for our transfer function. In this section, we describe the main numerical methods we use to create our emulator of $P_{\rm hydro}(k)/P_{\rm nbody}(k)$ using Gaussian Processes.

Gaussian Processes (GPs) (Rasmussen et al., 2004) are a versatile tool within machine learning and statistics, renowned for their efficacy in regression, interpolation, and uncertainty quantification. They provide a flexible framework for modeling functions along with their associated uncertainties (e.g. Bird et al., 2019; Rogers et al., 2019; Rogers & Peiris, 2021; Pedersen et al., 2021). Moreover, a Gaussian process emulator is computationally efficient, enabling its use within standard inference methodologies such as Markov Chain Monte Carlo (MCMC) for evaluations.

For the length scales we want to model, the CAMELS simulations have 39 linearly-spaced $k$ bins spanning the range $0.36<k<9.93$ h/Mpc. Hence, to model baryonic effects at small scales down to $k\sim 10$, our emulator delineates these 39 linearly-spaced $k$ bins. We treat the baryonic effects in each $k$ bin as an individual Gaussian process, enabling independent training for each bin, with every simulation serving as a training point for these $k$ bins.

At a specific point $X=[\Omega_{m},\sigma_{8},A_{\rm SN1},A_{\rm AGN},A_{\rm SN2},A_{\rm AGN2}]$ in the parameter space, a Gaussian process models the target function — $S(\mathbf{k}|X)\coloneqq P_{\rm hydro}(\mathbf{k}|X)/P_{\rm nbody}(\mathbf{k}|X)$ in our case — as an assembly of random variables that form a joint Gaussian distribution. This model is defined via $S(\mathbf{k}|X)\sim\mathcal{N}(0,K(X,X_{i}))$, where $X_{i}$ signifies the parameter values at the training simulations, and $K(X,X^{\prime})$ represents a covariance kernel.

The choice of the kernel function $K(X,X_{i})$ plays a pivotal role in characterizing the correlation or similarity between data points $X$ and $X_{i}$. This function serves as a prior, encapsulating the expected behavior of the underlying function — baryonic effects in our case — to be modeled. For our emulator, we adopt a Matérn 5/2 kernel, a generalized form of the radial basis function (RBF), defined by:

Here, $r=||X-X^{\prime}||_{2}$ represents the L2 distance between two data points, $\sigma_{0}^{2}$ denotes the variance parameter, and $\ell_{i}$ signifies the length scale of each input dimension, influencing the smoothness and range of correlations between data points. Our choice of this covariance kernel is motivated by the need for flexibility, achieved through the squared exponential kernel, the efficacy of linear interpolation, and allowing for noise in the training data.

At a test point $X_{*}$, the joint distribution of the test data $S(\mathbf{k}|X_{*})$ and training data $S(\mathbf{k}|X_{i})$ can be expressed as:

Here, $\sigma_{n}^{2}$ serves as a hyperparameter signifying Gaussian noise within the training data.

Consequently, our Gaussian process involves a total of 8 hyperparameters: 6 correlation lengths, $\sigma_{0}^{2}$, and $\sigma_{n}^{2}$. The hyperparameters are optimized by maximizing the marginal log-likelihood of the training data (Rasmussen et al., 2004).

The posterior predictive distribution over test data is obtained through:

The mean and variance derived from the posterior predictive distribution, using the training information at $X_{i}$, serve as estimators for the value and interpolation uncertainty associated with $S(X_{*})$.

We implement our emulator using tinygp (Foreman-Mackey et al., 2022), a Python library for GP Regression (GPR) built on top of the JAX library for numerical computing (Bradbury et al., 2018).

The GP model offers a broad prior across function space, enabling the modeling of the diverse baryonic effects we see in figure 1 without imposing strong prior constraints on its parameter dependencies. Since it is stochastic, this model provides predictions for the baryonic suppression beyond the training points, accompanied by associated uncertainties that can be integrated into our statistical model. This emulation methodology provides a robust approach for modeling and predicting the baryonic effects in simulations, enabling efficient and accurate interpolation, and quantification of uncertainties.

In this section, we use the trained Gaussian Process emulator to generate predictions for $P_{\rm hydro}(k)/P_{\rm nbody}(k)$ as a function of four astrophysical parameters - $A_{\rm SN1},A_{\rm AGN},A_{\rm SN2},A_{\rm AGN2}$ within their respective ranges. We emphasize that our emulator was trained on Astrid simulations at $z=0$, and therefore, the meaning of these astrophysical parameters is, in principle, associated with the Astrid subgrid physics model. However, to make our emulator generic and robust, from now on, we will consider these four astrophysical parameters as nuisance parameters that one needs to tune to reproduce the result of one particular hydrodynamic simulation.

We employ the differential evolution global optimizer from the SciPy library (Virtanen et al., 2020) to obtain the best-fit value of these nuisance parameters. This optimization technique is adept at exploring the parameter space to seek optimal solutions, especially in scenarios with complex, multi-dimensional parameter spaces. The differential evolution (Storn, 1996; Storn & Price, 1997; Price et al., 2005) method operates stochastically, offering a non-gradient approach to locating the minimum and can search through large volumes in parameter space.

We now show the accuracy of our emulator for simulations within and outside CAMELS, showing its precision to changes in simulation cosmology, feedback, subgrid physics, resolution, volume, and redshift. On top of that, we compare the accuracy of our emulator against other emulators in the literature. Finally, we demonstrate the emulator’s efficacy in creating field-level improvements when applied to the N-body fields of simulations within CAMELS, validating its potential for advancing field-level weak-lensing analyses using cosmological simulations.

Field-level approaches have the potential to extract all the available information from cosmological surveys. Modeling and marginalizing over baryonic effects at the field-level becomes a key ingredient in these efforts. In this work we have developed a new method to model baryonic effects for the total matter density field, the relevant quantity for weak lensing analyses.

The key finding in this work is that by computing the cross-correlation between the total matter density field in hydrodynamic and N-body simulations from thousands of simulations of the CAMELS project (see Figure 4) we conclude that baryonic effects weakly affect the phases of Fourier modes of the total matter density field down to scales as small as $k\sim 10~{}h{\rm Mpc}^{-1}$. This finding implies that baryonic effects will predominantly modify Fourier mode amplitudes. Thus, we can baryonify the total matter field of an N-body simulation by rescaling its Fourier mode amplitudes.

In this work we have built an emulator using Gaussian processes for the total to dark matter power spectrum ratio $S(k)$ that takes as input 2 cosmological parameters ($\Omega_{\rm m}$ and $\sigma_{8}$) and 4 astrophysical parameters ($A_{\rm SN1}$, $A_{\rm SN2}$, $A_{\rm AGN1}$, $A_{\rm AGN2}$). We have trained our emulator using 800 state-of-the-art hydrodynamic simulations from the Astrid suite of CAMELS. We then show that our emulator is able to reproduce the baryonic effects of thousands of hydrodynamic simulations that have different cosmologies, astrophysics, subgrid physics, resolutions, volumes, and redshifts within a few percent precision.

We have compared our emulator against others in the literature, such as BACCO, HMCode, and BCemu. We find that our emulator shares a similar level of accuracy with those, but it has a wider range of validity given that it has been trained on CAMELS, where variations in cosmology and astrophysics are very large. We also showed explicitly how using our method reduces the residuals when working at the field level by comparing the results of hydrodynamic simulations against baryonified N-body simulations. A limitation of using CAMELS is that the box size is very small, and baryonic effects may not be fully captured due to the absence of larger halos in these simulation boxes. This will need to be investigated in more detail using a suite of simulations varying box size.

Our emulator enables robust, cost-effective field-level weak lensing modeling and facilitates precise power spectra analyses at the two-point level. The versatility and accuracy of our GP baryonification emulator underscore its potential as a powerful tool in cosmological simulations, offering opportunities for enhanced analyses and deeper insights into baryonic effects in large-scale structures. However, whether this emulator suffices at the field level depends on the specifics of the observational program. For example, for weak lensing, this will require making weak lensing maps using ray-tracing techniques. The overall detectability of the effects that go beyond our field level emulator in the the weak lensing depends on the density of background galaxies and the observed area of the sky. This analysis goes beyond the purpose of this paper, and will be presented elsewhere.

We thank Raul Angulo and Aurel Schneider for their comments on the usage of the BACCO and BCemu emulators. DS thanks James Sullivan for helpful discussions on Gaussian processes. The work of FVN is supported by the Simons Foundation. The CAMELS project is supported by the Simons Foundation and the NSF grant AST 2108078.