# Emulation of baryonic effects on the matter power spectrum and constraints
from galaxy cluster data

Sambit K. Giri  and Aurel Schneider

###### Abstract

Baryonic feedback effects consist of a major systematic for upcoming weak-
lensing and galaxy-clustering surveys. In this paper we present an emulator
for the baryonic suppression of the matter power spectrum. The emulator is
based on the baryonification model, containing seven free parameters that are
connected to the gas profiles and stellar abundances in haloes. We show that
with the baryonic emulator, we can not only recover the power spectra of
hydro-dynamical simulations at sub-percent precision, but also establish a
connection between the baryonic suppression of the power spectrum and the gas
and stellar fractions in haloes. This connection allows us to predict the
expected deviation from a dark-matter-only power spectrum using measured X-ray
gas fractions of galaxy groups and clusters. With these measurements, we
constrain the suppression to exceed the percent-level at $k=0.1-0.4$ h/Mpc and
to reach a maximum of 20-28 percent at around $k\sim 7$ h/Mpc (68 percent
confidence level). As a further step, we also perform a detailed parameter
study and we present a minimum set of four baryonic parameters that are
required to recover the scale and redshift dependence observed in hydro-
dynamical simulations. The baryonic emulator can can be found at
https://github.com/sambit-giri/BCMemu.

##  1 Introduction

The next generation of cosmological surveys will not only lead to a massive
increase of data but also allow us to push towards smaller cosmological
scales. While these scales may carry new information regarding the underlying
cosmology, they are at the same time strongly nonlinear and affected by
astrophysical phenomena, making them very challenging to model. One of the
main concerns in that respect are active galactic nuclei (AGN) expelling large
amounts of gas into the intergalactic medium (IGM). Such AGN feedback
processes cannot be modelled from first principles, and they have been shown
to alter the clustering signal at scales relevant for cosmology [1, 2, 3, 4].

Based on various cosmological simulations, it has been shown that baryonic
(AGN) feedback leads to a suppression of the matter power spectrum of order 10
percent at non-linear scales, i.e. wave-modes beyond $k\sim 0.3$ h/Mpc [5].
While most simulations predict a similar general trend, they do not agree
between each other at the quantitative level. This should, however, not come
as a surprise, since AGN feedback processes are included via different sub-
grid models into the simulations. On top of that, the model parameters are
usually tuned to reproduce key galaxy properties without specific focus on
effects relevant for cosmology.

The specific connection between AGN feedback parameters and cosmology was
first investigated using the suites of hydro-dynamical simulations OWLS [6]
and cosmo-OWLS [7, 8]. These simulations revealed that the strength of the
feedback ejection has a direct influence on the matter power spectrum and on
the amount of gas within haloes. The latter can be directly measured using
X-ray observations of individual galaxy groups and clusters. As a consequence,
the AGN sub-grid model parameters in simulations can be directly tuned to
X-ray observations, a strategy pursued for the cosmological simulation suite
BAHAMAS [9, 10].

Recently it has been shown that, instead of running expensive hydro-dynamical
simulations, it is possible to include baryonic effects in gravity-only
$N$-body simulation outputs by perturbatively modifying the positions of
simulation particles around haloes [11]. This baryonification method is based
on parametrised gas and stellar profiles with shapes motivated by
observations. It comes with a handful of well-defined, physically meaningful
baryonic parameters, which can be systematically varied in a cosmological
context. Note furthermore, that the baryonification method has shown to
provide predictions that are accurate enough for the next-generation
cosmological surveys, at least at the level of the matter power spectrum [12,
13].

Although the baryonification method is many orders of magnitude faster than
running full hydro-dynamical simulations, it still requires the modification
of large $N$-body output files, making it unsuitable to be used directly in
Bayesian parameter inference pipelines. The main goal of the present paper is
therefore to construct an emulator for the baryonic suppression effect on the
matter power spectrum. This emulator is based on the same set of parameters
that the underlying baryonification model and it allows for a more efficient
exploration of the vast parameter space. Since there is no hope of determining
the effects of baryonic feedback from first principles, this is the only way
of constraining baryonic nuisance parameters in a cosmological context.

After building the emulator based on a training set from a slightly modified
version of the baryonification model of Ref. [12], we compare it to results
from a variety of hydro-dynamical simulations available in the literature. The
goal is to not only show that the simulated power spectra can be recovered at
percent accuracy, but also to highlight that the baryonification method is
able to predict the simulated power spectra based on information about the gas
and stellar abundance in haloes. If both these aspects are confirmed, we can
use the baryonification model for cross-correlation analyses [14].

Other questions we want to address in the present paper are: what is the
minimum amount of baryonic nuisance parameters that is able to recover hydro-
dynamical simulation results at acceptable accuracy? Can we reproduce the
redshift dependence visible in hydro-dynamical simulations? What can we learn
from X-ray observations regarding the expected baryonic suppression of the
power spectrum?

The paper is structured as follows: in Sec. 2 we provide a summary of the
baryonification model and describe the set-up of the emulator and quantify its
accuracy. Sec. 3 is dedicated to a comparison with hydro-dynamical
simulations. In Sec. 4 we introduce a case study where we use the measured gas
fractions of X-ray measurements of galaxy clusters to predict the power
suppression due to baryonic feedback. In Sec. 5 we present a way to reduce the
number of free model parameters, while maintaining an acceptable agreement
with simulations. We next study the redshift dependence of the model
parameters in Sec. 6. Finally we present our conclusions in Sec. 7.

##  2 Baryonification model and emulator

The method of baryonifying $N$-body simulations is based on the principle of
modifying halo profiles in simulation outputs by slightly displacing particles
around halo centres. Although this displacement is performed in a spherical
symmetric way, the resulting $N$-body outputs still contain the full
complexity of the nonlinear cosmological density field with non-spherical
haloes arranged along the cosmic web. The result of the baryonification
process is primarily visible in statistical measures such as the matter power
spectrum.

In this section we first summarise the important ingredients of the
baryonification model with specific emphasis on the free model parameters. We
then present the emulator of the baryonic power suppression. A specific
emphasis is placed on the production of the data set, the emulation method,
and the testing phase.

###  2.1 Model ingredients

The guiding principle of the baryonification method is the transformation of
gravity-only profiles into modified, dark-matter-baryon (dmb) profiles that
include the effects of dark matter, gas, and stars [11]. This process can be
described as

$\rho_{\rm nfw}(r)\,\,\,\rightarrow\,\,\,\rho_{\rm dmb}(r)=\rho_{\rm
clm}(r)+\rho_{\rm gas}(r)+\rho_{\rm cga}(r),$ (2.1)

where $\rho_{\rm nfw}$ denotes the truncated NFW profile [15, 16], while
$\rho_{\rm dmb}$ is composed of a collision-less matter ($\rho_{\rm clm}$),
gas ($\rho_{\rm gas}$), and central galaxy profile ($\rho_{\rm cga}$). The
collision-less profile mainly consists of dark matter but also includes all
stars in satellite galaxies and the stellar halo.

The stellar profile of the central galaxy is described by an exponentially
truncated power law

$\rho_{\mathrm{cga}}(r)=\frac{f_{\mathrm{cga}}(M)}{4\pi^{3/2}R_{\mathrm{h}}r^{2}}\exp\left[-\left(\frac{r}{2R_{\mathrm{h}}}\right)^{2}\right],$
(2.2)

where $f_{\rm cga}(M)$ describes the fraction of stars in the central galaxy
(defined in Eq. 2.7 below) and $R_{\mathrm{h}}$ is the stellar half-light
radius which is set to 0.015 times the virial radius $r_{\mathrm{vir}}$ [a
value in agreement with observations, see e.g. Ref. 17]. Note that for the
purpose of the matter power spectrum at cosmological scales, the details of
the stellar profiles are not important and the simplistic profile of Eq. (2.2)
is therefore sufficient.

The distribution of gas in the halo is parametrised by a cored double-power
law of the form

$\rho_{\mathrm{gas}}(r)\propto\frac{\left[\Omega_{b}/\Omega_{m}-f_{\mathrm{star}}(M)\right]}{\left[1+10\left(\frac{r}{r_{\mathrm{vir}}}\right)\right]^{\beta(M)}\left[1+\left(\frac{r}{\theta_{\mathrm{ej}}r_{\mathrm{vir}}}\right)^{\gamma}\right]^{\frac{\delta-\beta(M)}{\gamma}}},$
(2.3)

where $f_{\rm star}(M)$ is the total fraction of stars in haloes (see Eq. 2.7
below). The $\beta$-parameter has an additional halo-mass dependence which is
given by

$\beta(M_{c},\mu)=\frac{3(M/M_{c})^{\mu}}{1+(M/M_{c})^{\mu}}\ ,$ (2.4)

allowing the gas profile to become shallower than the NFW profile towards
smaller halo masses ($M\lesssim M_{c}$). Note, furthermore that the mass-
dependence of $\beta(M)$ differs slightly from the definition in Ref. [14,
12], in the sense that it is now strictly limited to positive values [see e.g.
Ref. 18].

The gas profile of Eqs. (2.3, 2.4) contains the five free model parameters
$M_{c}$, $\mu$, $\theta_{\rm ej}$, $\gamma$, and $\delta$. These parameters
are crucial for the modelling of baryonic feedback effects, as they determine
how the gas is distributed in and around haloes.

The collision-less matter profile is dominated by the dark matter component,
but also contains satellite galaxies and intra-cluster stars. At first, it is
given given by

$\rho_{\mathrm{clm}}(r)=\left[\frac{\Omega_{\mathrm{dm}}}{\Omega_{\mathrm{m}}}+f_{\mathrm{sga}}(M)\right]\rho_{\mathrm{nfw}}(r)$
(2.5)

where $\rho_{\rm nfw}$ is the truncated NFW profile and $f_{\rm sga}(M)$ is
the stellar fraction attributed to satellite galaxies and the intra-cluster
light. Note, that the collision-less profile may react to the other components
via the adiabatic relaxation prescription of Ref. [19]. As a consequence, it
does not retain its original NFW shape. A more detailed description of this
procedure can be found in Ref. [12].

After having discussed the three components making up the dark-matter-baryon
profile in Eq. (2.1), we will now introduce the stellar-to-halo fractions that
appear as part of these components. The stellar fraction of satellites (and
halo stars) is obtained by subtracting the central galactic fraction from the
total fraction of stars in haloes, i.e.,

$f_{\rm sga}(M)=f_{\rm star}(M)-f_{\rm cga}(M),$ (2.6)

where $M$ refers throughout the paper to the halo mass at 200 times the
critical density. The total stellar fraction and the fraction of stars in the
central galaxy are each given by a power law

$f_{i}(M)=0.055\left(\frac{M}{M_{s}}\right)^{-\eta_{i}}$ (2.7)

with $i=\\{{\rm star},{\rm cga}\\}$ and $M_{s}=2.5\times
10^{11}h^{-1}M_{\odot}$. Note that this simple description is influenced by
the well-known [20] relation, but only contains the large-scale decrease
(since haloes below $\sim 10^{12}$ M⊙ do not have to be modelled for our
purposes). We furthermore redefine $\eta_{\rm star}\equiv\eta$ and $\eta_{\rm
cga}=\eta+\eta_{\delta}$, ending up with 2 additional parameters
${\\{\eta,\eta_{\delta}\\}}$. All seven model parameters
${\\{M_{c},\mu,\theta_{\rm ej},\gamma,\delta,\eta,\eta_{\delta}\\}}$
describing the gas distribution and stellar abundances, are summarized in
Table 1.

###  2.2 Baryonic emulator

After summarising the baryonification model, we now introduce the baryonic
emulator for the matter power spectrum. We start by presenting the training
set, before discussing the emulation and testing method.

####  2.2.1 Dataset

To simulate the distribution of dark matter in our Universe, we run a $N$-body
simulation with $512^{3}$ particles in cubic volume of
$(256~{}h^{-1}\mathrm{Mpc})^{3}$ using the code Pkdgrav3 [21]. This setup has
been shown in Ref. [11] to be sufficient to study baryonic suppression effects
of the power spectrum at sub-percent accuracy. We assume a standard
5-parameter $\Lambda$CDM cosmology where $\Omega_{m}=0.315$,
$\Omega_{b}=0.049$, $h_{0}=0.674$, $n_{s}=0.96$ and $\sigma_{8}=0.811$. This
cosmology is within the constraints given by Planck [22]. Previous studies
have shown that the baryonic effects are insensitive to (reasonable) changes
of cosmological parameter values except for the mean baryon fraction
$f_{b}=\Omega_{b}/\Omega_{m}$ [14, 23]. Therefore, our baryonic emulator will
only depend on cosmology via $f_{b}$ which greatly simplifies the analysis.

Parameter name | Range
---|---
$\mathrm{log_{10}}M_{\mathrm{c}}$ | [11, 15]
$\mu$ | [0, 2]
$\theta_{\mathrm{ej}}$ | [2, 8]
$\gamma$ | [1, 4]
$\delta$ | [3, 11]
$\eta$ | [0.05, 0.4]
$\eta_{\delta}$ | [0.05, 0.4]
$f_{b}$ | [0.10, 0.25]
Table 1: List of parameter abbreviations and ranges used for the baryonic
emulator. The top five parameters ($\log_{\mathrm{10}}M_{c}$, $\mu$,
$\theta_{\rm ej}$, $\gamma$, $\delta$) are related to the gas profiles,
parameters 6 and 7 determine the stellar fractions ($\eta$, $\eta_{\delta}$),
while the bottom parameter ($f_{b}$) describes the relevant changes in
cosmology.

We create the training set by varying seven baryonification parameters along
with $f_{b}$. In Table 1, we list the ranges of these parameters
$\boldsymbol{\theta}_{\mathrm{bar}}$ on which our emulator will depend as well
as their respective ranges. This eight-parameter space is randomly sampled
with 2700 points. We then apply our baryonification algorithm at each point
before measuring the power spectra of the 2700 modified $N$-body outputs. The
power spectra are calculated between wave-modes $0.03~{}\mathrm{h/Mpc}\lesssim
k\lesssim 12.52~{}\mathrm{h/Mpc}$ in 232 bins. We finally divide each absolute
power spectrum by the dark-matter-only power spectrum of the original
simulation output in order to obtain the relative baryonic power suppression
($\boldsymbol{\mathcal{S}}$). This training set is used to create the emulator
(described in later sections).

Next to the training set, we independently create another data-set with 300
randomly selected points in our eight-dimensional parameter space. This
corresponds to the testing set that will be used to estimate the accuracy of
our emulator.

####  2.2.2 Method

We build our emulator by fitting a regression model to our data-set
$\\{\tilde{\boldsymbol{\theta}}_{\mathrm{bar}},\tilde{\boldsymbol{\mathcal{S}}}\\}$,
where $\tilde{\boldsymbol{\theta}}_{\mathrm{bar}}$ and
$\tilde{\boldsymbol{\mathcal{S}}}$ are baryonification parameter and power
spectrum suppression vector, respectively. We use kriging, which is an
interpolation based on Gaussian process (GP), to determine the regression
model. In the most generic form, the kriging model can be written as [e.g.
24],

$\displaystyle\boldsymbol{\mathcal{S}}=\sum_{i=1}^{p}f_{i}(\boldsymbol{\theta}_{\mathrm{bar}})+\mathcal{Z}(\boldsymbol{\theta}_{\mathrm{bar}})\
,$ (2.8)

where the first term is a linear combination of a known function $f_{i}$ and
the second term is a realisation of a stochastic process or GP,
$\mathcal{Z}(\boldsymbol{\theta}_{\mathrm{bar}})\sim\mathcal{N}(0,\mathcal{K})$.
We assume a linear polynomial for $f_{i}$. $\mathcal{K}$ is the kernel that
quantifies the similarity between the points in the dataset. In this work, we
use in the squared-exponential kernel defined as,

$\displaystyle\mathcal{K}(\boldsymbol{\theta}_{\mathrm{bar}},\boldsymbol{\theta}_{\mathrm{bar}}^{\prime})=\sigma^{2}\prod_{l=1}^{d}\exp\left(-A_{l}({\theta}_{\mathrm{bar},l}-{\theta}_{\mathrm{bar},l}^{\prime})^{2}\right)\
,$ (2.9)

where $d$ is the dimension of each data point in
$\boldsymbol{\theta}_{\mathrm{bar}}$. $A_{l}$’s and $\sigma$ are the hyper-
parameters that are determined during training. Once the hyper-parameters are
known, Eq. 2.9 can be solved at any new point of the parameter space. However,
this calculation can be computationally expensive especially for large data
sets. We therefore use a modified version of kriging, which is known as kernel
partial least squares (KPLS) [25]. Emulators built with KPLS are much faster
and require less memory without loss in accuracy compared to the standard
kriging method111This speed-up is obtained due to inclusion of partial least
square (PLS) method [26]. The PLS method is a tool for high dimensional
problems that searches the direction that maximises the variance between input
and output variables. The PLS information is integrated into the kriging
correlation matrix to scale down computational complexity.. This makes them
more suitable for Bayesian inference pipelines. We use the KPLS implementation
in the SMT222This package is available at https://smt.readthedocs.io. [27].

In this work, we build independent emulators at each redshift
($z=0,0.5,1,1.5,2$). Note that Ref. [23] took a different approach by adding
$z$ as an additional emulation parameter. We follow the strategy of multiple
emulators at different redshifts for two main reasons: First, this approach
allows us to implement different models for the redshift evolution of
individual parameters after the emulators are constructed. Indeed, any
redhsift evolution of individual parameters can be trivially retracted from
the emulators at different redshifts as long as the parameter values (that are
now evolved with redshift) remain in the range of the emulator. Second, not
emulating redshift allows us to reduce the dimension of the regression model.
Due to the curse of dimensionality, the required size of the training set
increases substantially due to increase in the dimension [28, 29].

In Sec. 5, we will specifically study the redshift evolution of emulation
parameters. More particularly, we investigate what parameters need to change
with redshift in order to reproduce the observed evolution of the baryonic
power suppression in hydro-dynamical simulation.

Our baryonic emulator is available online and can be used by installing
BCMemu333This package is available at https://github.com/sambit-giri/BCMemu..
Note that this emulator not only predicts the power spectrum suppression for a
given set of baryonic parameters but it also provides the uncertainty.

####  2.2.3 Testing emulation accuracy

Figure 1: The percentage error of predicting the $\mathcal{S}(k)$ in the
testing set with our emulators at various redshifts. The dark and light shaded
regions represent the 1-$\sigma$ and 2-$\sigma$ regions respectively.

In this section, we investigate the performance of the emulator using the
testing set. We quantify the accuracy of our emulator with the coefficient of
determination (also known as $R^{2}$ score) given as,

$\displaystyle
R^{2}=\sum_{j=1}^{n_{\mathrm{bins}}}\left[1-\frac{\sum_{i=1}^{n_{\mathrm{test}}}\left(y_{i}(k_{j})-\hat{y}_{i}(k_{j})\right)^{2}}{\sum_{i=1}^{n_{\mathrm{test}}}\left(y_{i}(k_{j})-\bar{y}(k_{j})\right)^{2}}\right]\
,$ (2.10)

where
$\bar{y}(k_{j})=\frac{1}{n_{\mathrm{test}}}\sum_{i=1}^{n_{\mathrm{test}}}y_{i}(k_{j})$.
$n_{\mathrm{test}}$ and $n_{\mathrm{bins}}$ are the size of the testing set
and number of $k$ bins, respectively. The $R^{2}$ score measures of how well
the model or emulator replicates the outcomes $y$. In order to determine the
accuracy of our emulator in predicting the testing set, $y_{i}$ is replaced
with $\mathcal{S}_{i}(k)$ in the above equation. We find that the $R^{2}$
score for the testing set is more than 0.99 at every redshift, which
corresponds to a very high accuracy level.

In Figure 1, we show the relative percentage error between the power spectra
suppression in the the testing set $\mathcal{S}_{\mathrm{test}}$ and the
corresponding predicted power spectra suppression
$\mathcal{S}_{\mathrm{pred}}$. The dark-black and light-black regions
correspond to the 1$\sigma$ and 2$\sigma$ regions, respectively. The 1$\sigma$
regions is within 1 percent and the 2$\sigma$ region within 3 percent compared
the true answer at all wave-modes.

##  3 Studying baryonic effects in hydro-dynamical simulations

In this section we test the validity of the baryonification method by
comparing it to hydro-dynamical simulations. This comparison is important not
only to assure that the matter power spectrum can be reproduced at percent
level but also to investigate the connection between halo properties and the
power spectrum. The latter is a prerequisite for combining cosmological data
with individual observations of gas around haloes, which consists of a
potentially powerful way to constrain baryonic feedback properties [30].

In this paper we consider the hydro-dynamical simulations OWLS [6], cosmo-OWLS
[7], BAHAMAS [9, 10], Illustris-TNG [31, 32, 33, 34], and Horizon-AGN [35].
The reason for selecting these simulations is that we could find published gas
and stellar fractions next to the matter power spectrum. This is a requirement
for the type of comparison analysis carried out below. A list of the
simulations used for the comparison is provided in Table 2.

Simulation name | cosmology | code | AGN feedback
---|---|---|---
BAHAMAS [9, 10] |  WMAP9 [36] |  Gadget3 [37] | thermal
cosmo-OWLS [7] |  WMAP7 [38] |  Gadget3 [37] | thermal
OWLS [6] |  WMAP3 [37] |  Gadget3 [37] | thermal
Illustris-TNG [31, 32, 33, 34] |  Planck [39] |  AREPO [40] | thermal and kinetic
Horizon-AGN [35] |  WMAP7 [38] |  RAMSES [41] | thermal and kinetic
Table 2: List of hydro-dynamical simulations considered in this work along
with the cosmology, the simulation code, and the implemented AGN feedback
mode.

###  3.1 Reproducing the power suppression signal from simulations

The first and most important requirement of our baryonification model is its
ability to reproduce the baryonic power suppression predicted by simulations.
We use the 7-parameter baryonification model to fit the emulated relative
power spectra $\mathcal{S}$ to the results from the hydro-dynamical
simulations. The best fits are obtained using the maximum likelihood
estimation (MLE) technique444We use the optimization methods implemented in
scipy [42] python package to maximize the likelihood. [43, 44]. For the
likelihood, we assume a Gaussian error of 1 percent at each wave-mode of the
$\mathcal{S}(k)$ we want to fit.

The best-fitting curves together with the results from the different hydro-
dynamical simulations are shown in Fig. 2. The agreement between hydro-
dynamical simulations and the baryonification model is better than one percent
over all wave-modes considered at low redshifts ($z<2$). Only for the $z=2$
results, the agreement becomes slightly worse but remains at the 1 percent
level until $k\sim 8$ h/Mpc.

The level of agreement presented in Fig. 2 is comparable to the results from
Ref. [13] that are based on an alternative implementation of the
baryonification method [45]. Note however, that both these results are based
on a 7-parameter model implementation. For practical purposes, it will be
important to further reduce the amount of free parameters. We will show
results for reduced sets of parameters in Sec. 5.

Figure 2: The relative matter power spectrum $\mathcal{S}=P/P_{\rm N-body}$
from all the hydro-dynamical simulations along with the best fitting curves
from the baryonification model (black-dashed lines) at three redshifts. The
agreement between model and simulations remains better than one percent over
all $k$-modes considered (see bottom panel) at low redshifts ($z<2$). At
$z=2$, the agreement degrades at small scales ($k\gtrsim 8$ h/Mpc).

###  3.2 Predicting the power spectrum from halo properties

(a) cosmo-OWLS 8.0

(b) cosmo-OWLS 8.5

(c) BAHAMAS

(d) OWLS

(e) Illustris-TNG

(f) Horizon-AGN

Figure 3: Comparison between the baryonification model (blue bands) and the
hydro-dynamical simulations cosmo-OWLS, BAHAMAS, OWLS, Illustris-TNG, and
Horizon-AGN (black lines). For each sub-figure, we show both the gas and
stellar fractions ($f_{\rm gas,500}$, $f_{\rm star,500}$) as a function of
$M_{500}$ (left-hand panel) as well as the baryonic power suppression
($\mathcal{S}$) as a function of $k$ (right-hand panel). Note that the
baryonification model is only fitted to the gas and stellar fractions from
simulations. The power suppression is then predicted by the baryonification
model (based on the fitted parameter values) and compared to the simulation
results. This procedure validates the ability of the baryonification model to
recover $\mathcal{S}(k)$ based solely on information of the fractions $f_{\rm
gas,500}$ and $f_{\rm star,500}$.

Another important feature of the baryonification model is its ability to
provide a direct link between individual halo profiles and the total matter
power spectrum. In this section we will investigate this connection by
comparing the baryonifiaction model to the simulations.

The question we want to examine is to what extent the baryonification model is
able to reproduce the baryonic power suppression from simulation based solely
on the simulated mean gas fractions ($f_{\mathrm{gas,500}}$) and stellar
fractions ($f_{\mathrm{star,500}}$) within haloes (more precisely within the
radius $r_{500}$). We calculate the stellar and gas fractions using the
baryonification profiles (see Eq, 2.3) and fit them to the simulated fractions
from OWLS, cosmo-OWLS, BAHAMAS, IllustrisTNG, and HorizonAGN. For the fitting
procedure, we consider 20 bins between $10^{13}h^{-1}M_{\odot}$ and $2\times
10^{14}h^{-1}M_{\odot}$ and assume a standard deviation of 1 percent on both
$f_{\mathrm{gas,500}}$ and $f_{\mathrm{star,500}}$. We ignore the bins where
the data from simulations are not available. We then run an MCMC to fit all
seven parameters simultaneously. From the resulting chain, we take parameter
values within the 68 percent credible regions of the posterior distributions
to determine the corresponding $\mathcal{S}(k)$ using our emulator.

The gas and stellar fractions fitted to each simulation (left-hand panels) as
well as the corresponding prediction for the power spectrum suppression
(right-hand panel) are shown in Fig 3. While the black lines represent the
results from simulations, the coloured bands indicate the 68-percent contours
of the baryonification model. Fig. 3 shows that the baryonification model is
able to reproduce the shapes and amplitudes of the simulated baryonic power
suppression to good accuracy when fitted to the gas and stellar fractions from
simulations. This is an important result validating the methodology of the
baryonification method. Note that the results also agree with the analysis of
Ref. [12], where a similar exercise has been carried out, albeit using a much
less general and robust fitting procedure.

A closer look at Fig. 3 reveals that the agreement is best for the OWLS and
cosmo-OWLS simulations. The BAHAMAS simulations can also be recovered well,
albeit only when correcting for resolution effects at small scales. Indeed, it
is shown in Ref. [46] that the BAHAMAS power spectrum is biased due to
resolution effects above $k\sim 5$ h/Mpc. We correct for this bias in Fig.
4(c) by multiplying the original curve with the resolution-corrected
difference provided in appendix A of Ref. [46].

For the Illustris-TNG simulation shown in Fig. 3(e), the constraints on the
power suppression remain rather poor. This is a direct result of the fact that
we only have information about gas and stellar fractions from relatively small
halo masses (of $M_{500}\lesssim 2\times 10^{13}~{}h^{-1}M_{\odot}$). Note,
furthermore, that while the power spectrum and the gas fraction are obtained
from the Illustris-TNG simulation with box-length of $L=100~{}h^{-1}$ Mpc, the
stellar fractions are form the lower-resolution simulation with
$L=300~{}h^{-1}$ Mpc. Although there remains some uncertainty regarding the
stellar fraction, we still find the $\mathrm{S}(k)$ of Illustris-TNG to be
within the 68 percent credible region predicted by the baryonification model.

For the Horizon-AGN simulation illustrated in Fig. 3(f), we use the
$f_{\mathrm{gas,500}}$ and $f_{\mathrm{star,500}}$ values of individual haloes
instead of binned quantities (see coloured points in the left-hand panel). The
resulting constraints on the power suppression are rather poor. We speculate
that this could be a consequence of the large scatter in the individual gas
and stellar fractions from the simulation.

###  3.3 Predicting halo properties from the power spectrum

In the previous section we have shown that the baryonification model is able
to reproduce the baryonic power suppression based on the stellar and gas
properties of haloes. We now want to investigate the inverse procedure, i.e.
to test if the model can reproduce the gas and stellar fractions from
simulations when fitted to the baryonic power suppression. For the fitting
procedure, we use a setup similar to the one above, allowing for a 1 percent
error on the wave-modes between 0.1 and 12.5 h/Mpc.

(a) cosmo-OWLS 8.0

(b) cosmo-OWLS 8.5

(c) BAHAMAS

(d) OWLS

(e) Illustris-TNG

(f) Horizon-AGN

Figure 4: Comparison between the baryonification model (green bands) and the
hydro-dynamical simulations cosmo-OWLS, BAHAMAS, OWLS, Illustris-TNG, and
Horizon-AGN (black lines). In contrast to Fig. 3, we now fit the
baryonification model to the power suppression $\mathcal{S}$ of each
simulation (left-hand panel of each sub-figure) with the goal to recover the
gas ($f_{\rm gas,500}$) and stellar fractions ($f_{\rm star,500}$) as a
function of $M_{500}$ (right-hand panel of each sub-figure). The figure
confirms that the baryonification method is able to recover simulated halo
properties based on information from the power spectrum.

The fits of the power suppression ($\mathcal{S}$) and the resulting
predictions for the gas ($f_{\mathrm{gas,500}}$) and stellar fractions
($f_{\mathrm{star,500}}$) are shown in Fig. 4 as green bands. The width of the
bands are obtained by drawing baryonification parameters from 68 percent
credible regions. The corresponding measures from hydro-dynamical simulations
are provided as black lines. As before, each simulation is accounted for in a
separate sub-figure (a-e).

The best results are again obtained for cosmo-OWLS (Fig. 4(a), 4(b)), BAHAMAS
(Fig. 4(c)), and OWLS (Fig. 4(d)), where the baryonification model is able to
recover the simulated gas and stellar fractions over the full mass range. Note
that for the BAHAMAS simulation, we fit to the resolution-corrected
suppression signal of the power spectrum (see Sec. 3.2 and Appendix A of Ref.
[46] for more details). Ignoring this correction would lead to a visible bias
between the gas and stellar fractions from the simulation and the
baryonification model.

For the Illustris-TNG and the Horizon-AGN simulations shown in Fig. 4(e) and
4(f), we find good agreement for the gas fraction predicted by the
baryonificatin model. The stellar fraction, on the other hand, is predicted to
be somewhat smaller than the simulation results. At least for Illustris-TNG,
this could be due to the fact that the stellar fraction come from a lower-
resolution simulation (see discussion in Sec. 3.2).

###  3.4 Predicting the cluster baryon fraction from the power spectrum

The comparison tests carried out in the two previous sections are limited to
the six simulations with published data for $f_{\rm gas,500}$ and $f_{\rm
star,500}$ as function of $M_{500}$ (presented in Table 2). However, if we
restrict ourselves to the mean baryon fraction of clusters, then there are
more simulations to incestigate. Indeed, Ref. [46] published the baryonic
fractions at $M_{500}=10^{14}$ M⊙, $f_{b}(\sim 10^{14}M_{\odot})$, for the six
simulations introduced above plus the additional simulations BAHAMAS low-AGN,
BAHAMAS high-AGN, and Illustris555The BAHAMAS low-AGN and BAHAMAS high-AGN
simulations assume 0.2 dex lower and higher AGN heating temperatures compared
to the original BAHAMAS. For details about the Illustris simulations, see
Refs. [47, 48, 49]). We will now use this extended data-set to check how well
the baryonification model is able to predict the value of $f_{b}(\sim
10^{14}M_{\odot})$ when fitted to the baryonic power suppression
($\mathcal{S}$).

Figure 5: Mean cluster baryon fractions for haloes with $M_{500}\sim
10^{14}M_{\odot}$ in various hydrodynamical simulations (coloured symbols).
While the stars represent the values published in Ref. [46], the pink pentagon
is obtained from Ref. [50]. The black data points with error bars show
predictions from our baryonification model obtained by fitting the model to
the baryonic power suppression from each simulation. The size of the error-
bars corresponds to the uncertainty-range of the fitting procedure (see text
for more details). Note that all simulation results are well recovered by the
baryonification model, with only the Illustris results lying outside the
1$\sigma$ uncertainty range. The empty squares give the values estimated in
Ref. [13] based on a different implementation of the baryonification method.

In Fig. 5, we show the mean baryonic fractions $f_{b}(\sim 10^{14}M_{\odot})$
provided by Ref. [46] (coloured stars) and we compare them to the results from
the baryonification model (black data points). The latter are obtained by
fitting the baryonification model to the power suppression ($\mathcal{S}$) of
the corresponding simulations, before predicting $f_{b}(\sim
10^{14}M_{\odot})$ based on the obtained best-fitting parameters. The shown
error-bars correspond to the 68 percent confidence level based on the fits to
$\mathcal{S}$ assuming 1 percent errors on each $k$-mode in the range $k\sim
0.1-12.5$ h/Mpc (same as in Sec. 3.3).

Fig. 5 highlights the very good agreement between the cluster baryon fractions
predicted by the simulations and obtained with the baryonification model.
Nearly all of the simulation results comfortably lie within the 68 percent
uncertainty regime. Only Illustris has a cluster baryon fraction that is about
$2\sigma$ higher than the prediction by the baryonification model. Note,
however, that the difference is reduced to a little more than 1$\sigma$ when
comparing to the result from Ref. [50] (pink pentagon) instead of Ref. [46]
(pink star). In Ref. [50] it is furthermore shown that, due to the small
simulation box ($75~{}h^{-1}$ Mpc) of the Illustris simulation, the recovered
function $f_{b}(M_{500})$ becomes very noisy around the relevant scale of
$10^{14}$ $M_{\odot}$. This points towards the underlying reason for the
difference between [46] and [50], letting us conclude that the cluster baryon
fraction from Illustris remains somewhat uncertain.

Note that for the Illustris and the three BAHAMAS simulations, we have added
the published results from Ref. [13] (empty squares) that are based on a
different implementation of the baryonification method. While we obtain very
similar predictions for BAHAMAS and BAHAMAS low-AGN, there are notable
differences for the BAHAMAS high-AGN and Illustris simulations, with the
results from Ref. [13] significantly deviating from the simulations and from
our results.

##  4 Constraints on baryonic effects with X-ray observations

In the previous section we have validated the baryonifcation method as a tool
to recover the baryonic power suppression ($\mathcal{S}$) based on information
about the properties of haloes around galaxy groups and clusters. We will now
use this connection to predict $\mathcal{S}$ based on observed gas and stellar
fractions from X-ray and optical observations. In Sec. 4.1 we first quantify
the constraints on the baryonic parameters, before predicting the shape of the
power spectrum in Sec. 4.2. Throughout this analysis, we will fix all
cosmological parameters to the values from Planck [39].

Figure 6: Gas and stellar fractions determined from X-ray and optical
observations of galaxy groups and clusters [51, 52, 53, 54, 17]. We have
binned the gas and stellar fractions into four and two bins, respectively. The
error bars on the binned data includes both the observational errors on the
points as well as their scatter.

###  4.1 Constraints on baryonic parameters

The parameters of the baryonification model will be constrained using gas and
stellar fractions of individual galaxy groups and clusters. We thereby follow
a similar strategy to that in Ref. [12], but we will perform a more systematic
multi-parameter analysis using Bayesian inference.

Fig. 6 shows all the observational data in terms of gas (blue points) and
stellar fractions (orange points) that is used in this work. The data is taken
from Refs. [51, 52, 53, 54, 17] and contains estimates for the observational
errors. We furthermore expect the scatter of the data points to depend on the
stochastic nature of structure formation, where haloes of the same mass may
contain different amounts of gas and stars.

For our analysis it is convenient to bin the gas and stellar fractions of
individual clusters into four and two halo mass bins, respectively. These bins
are shown as black data points in Fig. 6. The error on the binned data
accounts for both the observational error and scatter of the points inside
each bin. In our analysis, we assume the errors on the binned data to be
uncorrelated.

In order to obtain gas and stellar fraction with X-ray observations, the total
halo mass has to be estimated assuming the gas to be in hydrostatic
equilibrium. This assumption is, however, known to be sub-optimal because of
non-thermal pressure contributions from in-falling gas. Various studies have
shown that mass estimates based on hydrostatic equilibrium are likely to
underestimate the true cluster mass by 5-30 percent [e.g. 55, 56, 57].
Following the usual practice, we quantify this uncertainty using a hydrostatic
mass bias ($b_{\mathrm{hse}}$) that connects the estimated to the true mass
via the relation

$\displaystyle M_{\mathrm{500,hse}}=(1-b_{\mathrm{hse}})M_{\mathrm{500}}\ .$
(4.1)

Here $M_{\mathrm{500,hse}}$ is the total mass estimated assuming hydrostatic
equilibrium and $M_{\mathrm{500}}$ is the true total mass inside of
$r_{\mathrm{500}}$.

Figure 7: Posterior contours of the seven parameters from the baryonification
model fitted to X-ray and optical data of gas and stellar fraction within the
halo radius $r_{500}$. Here we have fixed all cosmological parameters to the
best-fitting values from Planck and we have marginalised over the hydrostatic
mass bias parameter ($b_{\mathrm{hse}}$).

We infer constraints in the baryonic feedback processes by running MCMC using
our seven-parameter baryonification model. For all baryonic parameters we
consider flat priors within the range given in Table 1. For the hydrostatic
mass bias ($b_{\mathrm{hse}}$) we use a Gaussian prior with mean and standard
deviation of 0.26 and 0.07, respectively. These values are taken from the
analysis of Ref. [58] that is based on CMB lensing and Sunyaev-Zeldovich data.
Note that these values are in good agreement with other direct estimates based
on X-ray and lensing data [e.g. 59].

In Fig. 7, we plot the resulting posterior distribution from our MCMC run
(marginalising over $b_{\mathrm{hse}}$). The parameter
log${}_{\mathrm{10}}M_{c}$ is constrained to a value around 12.6, which is
similar to the standard value used in previous studies [12, 14]. For the
second parameter $\mu$, the data prefers low values, which is a consequence of
the fact that the X-ray gas fraction is only weakly dependent on halo mass.
The parameters $\theta_{\mathrm{ej}}$ and $\gamma$ are not well constrained by
the data. The parameter $\delta$, on the other hand, is pushed to rather low
values, corresponding to a slowly falling gas profile beyond the virial radius
(see Eq. 2.3). This means that large amounts of gas are pushed outside of the
halo, providing evidence for strong feedback processes.

Regarding the stellar parameters, Fig. 7 shows a clear trend between the
parameters $\eta$ and $\eta_{\delta}$. For large values of $\eta$,
$\eta_{\delta}$ is pushed to small values and vice versa. This means that
small total stellar fractions are only allowed if the satellite fraction is
suppressed. Large total stellar fractions, on the other hand, require a large
fraction of satellite galaxies.

Figure 8: Constraints on $f_{\mathrm{gas}}$, $f_{\mathrm{star}}$ and
$\mathcal{S}(k)$ based on the binned X-ray observations. Left panel gives the
binned X-ray observations of gas (red points) and stellar (green points)
fractions for different cluster masses. The light and dark shaded regions show
the 95 percent and 68 percent credible regions respectively determined from
our MCMC run. Right panel shows the 95 percent and 68 percent credible region
of our constraint on $\mathcal{S}(k)$ with light and dark shaded regions.

###  4.2 Constraints on the power spectrum

Using the constraints on the baryonification parameters (shown in Fig. 7), we
can now investigate the expected shape of the baryonic power suppression
($\mathcal{S}$). The result can then be compared to the predictions from
hydro-dynamical simulations, providing an indirect test for their validity
regarding cosmology.

In the left-hand panel of Figure 8, we show the $f_{\mathrm{gas,500}}$ and
$f_{\mathrm{star,500}}$ as blue and orange bands. They include the 68 percent
(dark band) and 95 percent (light band) credible regions obtained via the
posterior distribution of the baryonification parameters (see Fig. 7). The
plot confirms that the binned X-ray and optical data used for the fitting
procedure (shown as red and green symbols) is within the credible region
obtain by the MCMC run.

In the right-hand panel of Figure 8, we show the predicted baryonic power
suppression ($\mathcal{S}$) obtained by applying the baryonification model to
the X-ray and optical data from individual galaxy-groups and clusters. The
dark and light bands again correspond to the 68 and 95 percent credible
regions. At the 68 percent confidence level, the baryonic power suppression
$\mathcal{S}$ is predicted to exceeds the percent level between $k\sim 0.1$
and $0.4$ h/Mpc growing to a maximum 21-28 percent at around $k\sim 7$ h/Mpc.
At the 95 percent level, suppression exceeds one percent somewhere between
$k\sim 0.1$ h/Mpc and $1.3$ h/Mpc reaching a maximum of 14-38 percent beyond
$k\sim 7$ h/Mpc.

The right-hand panel of Fig. 8 also shows the predicted baryonic power
suppression signal from various hydro-dynamical simulations. We find that
cosmo-OWLS 8.0, cosmo-OWLS 8.5, BAHAMAS, and OWLS are consistent with the
predictions from X-ray data at the 95 percent confidence level666Note that for
cosmo-OWLS 8.0, cosmo-OWLS 8.5, and BAHAMAS this should not come as a
surprise, since these simulations have been calibrated to a similar collection
of X-ray observations.. The Illustris-TNG and Horizon-AGN simulations, on the
other hand, do not agree with the X-ray informed predictions from the
baryonification model. They show a suppression signal that is shifted to
higher $k$-values and is in apparent disagreement with the signal suggested by
X-ray observations.

##  5 Reducing the number of baryonic parameters

(a) Five parameters

(b) Three parameters

(c) Single parameter

Figure 9: Fits to the simulated baryonic power suppression signal
($\mathcal{S}$) for a reduced number of five (a), three (b), and one (c)
parameter. From left to right we show results from redshift 0, 1 ,and 2, which
are fitted independently. coloured lines correspond to simulation results,
dashed black lines to the emulator. For the one parameter model, we only fit
up to $k=3$ h/Mpc, beyond this range the emulator results are shown as dotted
lines. For the same plot with the original 7 free parameters, see Fig. 2, for
details about the parameter reduction procedure, see text and Appendix A.

Until now we allowed all the seven parameters in our baryonification model to
vary. An efficient analysis of cosmological survey data requires, however, a
minimum amount of baryonic nuisance parameters. The objective of this section
is therefore to find a model with the smallest number of free parameters that
still provides a decent fit to the results from hydro-dynamical simulations.

In a first step, we reduce the number of parameters by two, ending up with a
five-parameter model. We determine which parameters to fix, by applying a
model selection procedure based on the Akaike information criterion (AIC) [60,
61] and the Bayesian information criterion (BIC) [62, 61]. A detailed
description of the procedure is provided in Appendix A.

Our analysis guides us to fix $\delta$ and $\eta$ to 7 and 0.2, respectively.
With this five-parameter model, we fit the baryonic power suppression
($\mathcal{S}$) to the predicted signal from the hydro-dynamical simulations
introduced in in Sec. 3. The resulting fits (black dashed lines) and
simulation results (coloured lines) are shown in Figure 9(a). We find that the
five-parameter model provides fits that are nearly as accurate as the original
seven-parameter model (see Fig 2). The agreement is about one percent or
better over all scales, except for the $z=2$ case, where we get larger than
percent deviations at $k\gtrsim 3$ h/Mpc.

As a further step, we reduce the number of parameters from five to three. We
do this by fixing $\gamma$ and $\mu$ to 2.5 and 1.0, respectively. The
resulting fits for the three-parameter model are shown in Figure 9(b). For
redshift 0 and 1, the fits agree to better than 1.5 percent up to $k\sim 10$
h/Mpc. For the $z=2$ case, the agreement is better than 3 percent for
$k\lesssim 7$ h/Mpc.

Finally, we further reduce the number of free parameters to 1 by fixing
$\theta_{\mathrm{ej}}$ and $\eta_{\delta}$ to 3.5 and 0.2, respectively. The
resulting fits from this 1-parameter model are shown in Figure 9(c). They are
significantly less accurate than the ones from the 3-, 5-, and 7-parameter
models, showing discrepancies of order 5 percent below $k\sim 8$ h/Mpc
(degrading further at higher $k$-modes). We therefore conclude that a model
with a single free parameter is insufficient to accurately describe the
simulation results for scales up to $k\sim 10$ h/Mpc. Note, however, that a
single parameter model remains a valid option for large scales below $k=2$
h/Mpc. In Appendix B we show that for a limited fitting range, 3-percent
agreement up to $k\sim 2$ h/Mpc can be obtained with a one-parameter model.

##  6 Redshift dependence of baryonic effects

So far we have ignored the potential redshift dependence of our model
parameters. Instead, all models have been fitted to the simulation results at
each redshift independently. However, a cohesive framework, where observations
at multiple redshifts can be studied simultaneously, may require the
baryonification parameters to become functions of redshift [see e.g. Ref. 18].

In order to study potential redshift dependencies, we will consider the three-
parameter model introduced above, as it yields almost same level of accuracy
compared to the original seven-parameter model in reproducing the baryonic
power suppression signal. In terms of simulations, we focus on BAMAHAS, cosmo-
OWLS, and OWLS, which all provide published power spectra at redshifts 0.0,
0.5, 1.0, 1.5, and 2.0.

Figure 10: Redshift dependence of the free model parameters (assuming the
3-parameter model), when fitted to each redshift bin individually. Different
colours correspond to different fits to the simulations cosmo-OWLS 8.0 (blue
upward-triangles), cosmo-OWLS 8.5 (orange downward-triangles), BAMAHAS (green
stars), and OWLS (brown diamonds).

In Fig. 10, we show the parameter values (for the three model parameters
$\log_{10}M_{c}$, $\theta_{\rm ej}$, $\eta_{\delta}$) producing the best
independent fits at different redshifts. For the parameter $\log_{10}M_{c}$
(left-hand panel) a clear trend is visible, the values decreasing with
increasing redshifts. The decrease is most pronounced for cosmo-OWLS 8.0,
cosmo-OWLS 8.5 and BAHAMAS, while it is only week for OWLS. For the parameter
$\theta_{\rm ej}$ (middle panel), no clear and coherent trend is discernible.
The values seem to oscillate and then increase, except for OWLS where they
slightly decrease towards higher redshifts. Regarding the third parameter
$\eta_{\delta}$ (right-hand panel), a decreasing trend can be observed again.
The slope of this decrease is similar between the simulations except for the
cases where the parameter is pushed towards the edge of the allowed prior
range (which is at $\eta_{\delta}=0.05$).

Motivated by the trends illustrated in Fig. 10, we allow the parameters to
vary with redshift, following the simple relation

$\displaystyle X(z)=X_{0}(1+z)^{-\nu_{X}}.$ (6.1)

Here $X$ stands for either $\log_{10}M_{c}$, $\theta_{\rm ej}$, or
$\eta_{\delta}$. With this extended parametrisation, we have formally doubled
the number of parameters, each former parameter now being described by its
value at $z=0$ (given by $X_{0}$) and its power-law index (given by $\nu_{X}$)
describing the redshift evolution. However, based on the trends visible in
Fig. 10 and in order to avoid an excessive growth of the parameter space, we
fix the parameters $\nu_{\theta_{\rm ej}}$ and $\nu_{\eta_{\delta}}$ to 0 and
0.6, respectively. This means that we are left with a four-parameter model
with $\log_{10}M_{c,0}$, $\theta_{\rm ej,0}$, $\eta_{\delta,0}$, and
$\nu_{\log_{10}M_{c}}$ as freely variable parameters.

In Fig. 11, we compare the new 4-parameter model (black dashed lines) to the
results from simulations (coloured lines). This time, we do not fit the signal
at different redshifts independently, but we perform a fit to all three
redshift bins, allowing all four parameters ($\log_{10}M_{c,0}$, $\theta_{\rm
ej,0}$, $\eta_{\delta,0}$, $\nu_{\log_{10}M_{c}}$) to vary simultaneously. We
obtain a 3-percent agreement for all redshifts up to $k\sim 5$ h/Mpc. For
lower redshifts ($z\leq 1$) the agreement is at the level of 3 percent up to
$k\sim 10$ h/Mpc.

Next to the 4-parameter model with redshift dependence, we also show the
results for the redshift independent 3-parameter model as dotted black lines.
This time, however (and unlike Fig. 9(b)), the model parameters are only
fitted to the simulation results at $z=0$. This means that the black dotted
lines indicate the expected redshift dependence of $\mathcal{S}$ provided none
of the model parameters are allowed to vary with redshift. A comparison
between the dashed and dotted lines also illustrates the obtained gain when
adding one additional parameter describing the explicit redshift dependence of
$\log_{10}M_{c}$.

Figure 11: Fits to the simulated baryonic power suppression ($\mathcal{S}$)
assuming an explicit redshift dependence of the model. The coloured lines show
the simulation results, while the dashed black lines correspond to the
4-parameter model which includes a power-law redshift dependence of the
$\log_{10}M_{c}$ parameter (see text). Note that the model is fitted to the
simulation results from all redshifts simultaneously. For comparison, we also
show the results from the redshift-independent 3-parameters model fitted
solely to the $z=0$ results.

At the end of this section we want to emphasize that, while at least one
redshift-dependent model parameter is clearly required to describe the
simulation, there is currently no direct evidence from observations that model
parameters have to vary with redshift. Indeed, Refs. [63, 64] find no
significant redshift dependence in the X-ray cluster gas fraction between
$z=0$ and 1. Note, however, that the uncertainties on these measurements
remain rather large.

##  7 Conclusions

In this work we built an emulator for the power suppression signal
$\mathcal{S}=P/P_{N-\rm body}$, which is primarily caused by baryonic feedback
processes pushing gas out of galaxies into the intergalactic space. The
baryonic effects are modeled using a slightly modified version of the
baryonification method described in Ref. [12]. The baryonification method
relies on parametrised halo profiles and is implemented as a post-processing
step into gravity-only $N$-body simulations. While being many orders of
magnitude faster than running full hydro-dynamical simulations,
baryonification is not directly suitable for Bayesian inference using MCMC, as
it involves reading and writing large $N$-body output files. Constructing an
emulator based on a large training set of baryonified simulations is therefore
a natural and efficient approach.

In a first step, we build five independent emulators at different redshifts
($z=0$, 0.5, 1.0, 1.5, 2.0) assuming eight free parameters. Five parameters
are related to the gas profile ($\log M_{c}$, $\mu$, $\theta_{\rm ej}$,
$\gamma$, $\delta$), two to the stellar fractions ($\eta$, $\eta_{\delta}$),
and one parameter ($f_{b}=\Omega_{b}/\Omega_{m}$) captures effects from
cosmology. The emulators are constructed using a modified kriging method,
known as _kernel partial least squares_ (KPLS), which is a fast and memory-
efficient emulation technique well suited for MCMC analyses. We show that the
emulators have an accuracy of better than 1 percent at $z=0$, better than 2
percent at $z=0.5$, 1.0, 1.5, and better than 2.5 percent at $z=2$ for 95
percent of the cases and over all $k$-modes considered ($k\leq 12.5$ h/Mpc).

The emulators (and the underlying baryionfication method) are validated using
suites of hydro-dynamical simulations, such as cosmo-OWLS 8.0, cosmo-OWLS 8.5,
BAHAMAS, OWLS, Illustris-TNG, and Horizon-AGN. We find that the emulators can
fit all simulation results to better than one percent for $z=0$, 0.5, 1, and
1.5 over all scales ($k\leq 12.5$ h/Mpc). Only for $z=2$, the agreement
degrades somewhat but stays at the 3 percent level until $k\lesssim 5$ h/Mpc.

One of the important features of the baryonification model is that its
parameters are connected to halo profiles and are therefore physically
meaningful. This makes it possible to constrain baryonic parameters using
observational information from the gas and stellar distribution around
clusters. We investigate this connection by testing the ability of the
baryonification model to recover the power suppression signal from simulations
when provided with information on the gas and stellar fraction of haloes
alone. We show that the model is indeed capable to recover the power spectrum
based solely on halo gas and stellar fractions of galaxy groups and clusters.
Inversely, it can also predict the gas and stellar fractions based on the
matter power spectrum alone. This important cross-check opens the door towards
cross-correlation studies between e.g. gas and weak-lensing observations.

As a first step in that direction, we perform a case study where we use X-ray
and optical cluster observations to predict the baryonic power suppression
with the help of the baryonification model. We show that X-ray observations
point towards a relatively strong baryonic suppression signal, exceeding the
percent suppression at $k\sim 0.1-0.4$ h/Mpc and growing to a maximum of more
than 20 percent at around $k\sim 7$ h/Mpc. Such a suppression signal is in
broad agreement with OWLS, cosmo-OWLS, and BAHAMAS, but significantly stronger
than the predictions from e.g. Illustris-TNG and Horizon-AGN.

In a further step, we reduce the amount of model parameters from seven to five
to three and to one, investigating the level of agreement that can be
maintained between the simulations and the emulators. Intriguingly, we find
the agreement to remain very good for the 5 and 3 parameter models, with
deviations from simulations that stay well below 3 percent for most redshifts
and $k$-ranges. Only the one parameter model performs significantly worse,
with 3 percent agreement up to $k\sim 2$ h/Mpc and strong deviations at
smaller scales.

Finally, we study the redshift dependencies of the baryonification parameters,
that are required to reproduce the simulated power suppression signal over
cosmic time. Using our three-parameter model, we find that we need at least
one additional free parameter to describe the inherent redshift dependence of
the aforementioned hydro-dynamical simulations. We thus provide a final four-
parameter model that can reproduce simulation results at the level of a few
percent for wave-modes below $k\sim 10$ h/Mpc and redshifts below $z=2$.

The baryonic emulator from this paper is available at
https://github.com/sambit-giri/BCMemu along with the documentation.

## Acknowledgments

This work is supported by the Swiss National Science Foundation via the grant
PCEFP2_181157. All the MCMC runs were done using the emcee [65] package. For
analysis and plotting, we have used the following software packages: numpy
[66], matplotlib [67] and corner [68].

##  Appendix A Method used for reducing the number of baryonic parameters

Here we describe the procedure used to find the values of the parameters to
fix in order to reduce the number of baryonification model parameters. In this
study, we fix two parameter at a time, going from the original 7 baryonic
parameters to 5, 3, and 1. For the first step (from seven to five parameters),
we choose two out of the seven parameters and bin them within their ranges
(see Table 1) into 30 bins. After fixing these two parameters to a certain
value, we fit this five parameter model to the $\mathcal{S}(k)$ from hydro-
dynamical simulations. We sum the likelihood values of the best fitting model
of all the simulations. Note that we use $\mathrm{S}(k)$ from multiple
redshifts and we fit them independently.

Figure 12: Best fitting 5-parameter models for different choices of fixed
parameters. Each sub-panel shows a combination of fixed parameters with
$-\log\mathcal{L}$ summed over all simulations and redshifts considered. Dark-
blue (dark-red) areas correspond to 5-parameter models that provide the best
(worst) global fit to the simulations. Based on this plot, we fix $delta=7$
and $\eta=0.2$ for the 5-parameter model used in the main text. Figure 13:
Best fitting 3-parameter models for different choices of fixed parameters
(starting from the 5-parameter model selected above). Each sub-panel shows a
combination of fixed parameters with $-\log\mathcal{L}$ summed over all
simulations and redshifts. Dark-blue (dark-red) areas correspond to
3-parameter models that provide the best (worst) global fit to the
simulations. Next to fixing $\delta=7$ and $\eta=0.2$ (see construction of the
5-parameter model), we fix $\mu=1$ and $\gamma=2.5$ to obtain the 3-parameter
model used in the main text. Figure 14: Best fitting 1-parameter models for
different choices of fixed parameters (starting from the 3-parameter model
selected above). Each sub-panel shows a combination of fixed parameters with
$-\log\mathcal{L}$ summed over all simulations and redshifts. Dark-blue (dark-
red) areas correspond to the 1-parameter models that provide the best (worst)
global fit to the simulations. Next to fixing $\delta=7$, $\eta=0.2$, $\mu=1$
and $\gamma=2.5$ (see construction of the 3-parameter model), we fix
$\theta_{\rm ej}=3.5$ and $\eta_{\delta}=0.2$ to obtain the 1-parameter model
used in the main text.

For the model selection, one can use the Akaike information criterion (AIC)
and the Bayesian information criterion (BIC) that is given as,

$\displaystyle\mathrm{AIC}=2N_{p}-2\mathrm{log}\mathcal{L}\ ,~{}~{}~{}$ (A.1)
$\displaystyle\mathrm{BIC}=N_{p}\mathrm{log}(n)-2\mathrm{log}\mathcal{L}\ .$
(A.2)

In the above equations, $N_{p}$, $n$ and $\mathcal{L}$ are the number of
parameters, number of observations fitted and likelihood value respectively.
The Occam’s razor argument is built into these criteria, and models with
minimised AIC and BIC are preferred. This method is similar to hyper-parameter
optimization used in the field of machine learning to find the best parameter
value in the trained model [69]. See Refs. [61, 70] to read more about model
selection approaches.

Fig. 12 shows the negative sum of all the likelihood values of the best fit
models to the hydro-dynamical simulations (for description of the simulations,
see Sec. 3). As $N_{p}$ and $n$ is the same for the models fitted, the AIC and
BIC only depend on the log$\mathcal{L}$. In this figure, the dark-blue regions
correspond to the best five-parameter models. Based on Fig. 12, we propose to
fix $\delta$ and $\eta$ to 7 and 0.2 respectively. In Fig. 9(a), we show the
best fit models to the hydro-dynamical simulations using the five parameter
representations.

As a next step, we further reduce the numbers of parameters from 5 to 3. Fig.
13 shows the negative sum of likelihood values ($-\log\mathcal{L}$) for all
3-parameter models with fixed pairs of parameters fitted to the simulations.
Dark-blue areas correspond to models that provide the best global fits over
all simulations and redshifts considered. Based on Fig. 13, we propose to fix
$\mu$ and $\gamma$ to 1 and 2.5, respectively, ending up with the 3-parameter
model used in the main text.

The final step consists of reducing the amount of free parameters from 3 to 1.
In Fig. 14, we show the values of $-\log\mathcal{L}$ for all 1-parameter
models fitted to the simulations. Dark-blue regions again correspond to the
models that provide the best overall fit to the simulations. Based on the data
shown in Fig. 14, we fix $\theta_{\mathrm{ej}}$ to 3.5 and $\eta_{\delta}$ to
0.2 to obtain the 1-parameter model used in the main text of this paper (see
Sec. 5). Note that, contrary to the reductions from 7 to 5 and 5 to 3
parameters, the AIC and BIC values are significantly larger for the 3 to 1
parameter reduction. This suggests that the 3-parameter baryonification model
consists of the best compromise between accuracy and small number of free
parameters optimal for Bayesian inference analyses.

##  Appendix B Single parameter baryonification model

Figure 15: Fits to the simulated baryonic power suppression signal
($\mathcal{S}$) at large scales ($k\lesssim 3$ h/Mpc) for one-parameter
baryonification model. From left to right we show results from redshift 0, 1,
and 2 (which are fitted independently). Coloured lines correspond to
simulation results, dashed black lines to the emulator. We only fit up to
$k=3$ h/Mpc, beyond this range the emulator results are shown as dotted lines.

Our one-parameter baryonification model struggles to reproduce the baryonic
power suppression ($\mathcal{S}$) from simulations when fitted to scales
between $k=0.1$ h/Mpc and 12.5 h/Mpc (see Fig. 9(c) and description in Sec.
5). In this appendix, we show that $\mathcal{S}$ can be reproduced to much
better accuracy when only fitted to large scales, i.e. until $k$-modes of
$k=3$ h/Mpc.

In Fig. 15, we show the fits to the power spectra suppression for $k\leq 3$
h/Mpc. The three panels show the results at redshift 0, 1, and 2. When
focusing on large scales alone (with $k\lesssim 2$ h/Mpc), the agreement
between our model and the simulations stays within 3 percent for all
redshifts. At smaller scales (larger $k$-modes) the curves diverge from the
simulation results, leading to substantial errors of up to 10-20 percent at
$k\sim 10$ h/Mpc.

## References

  * [1] M. P. van Daalen, J. Schaye, C. M. Booth, and C. Dalla Vecchia, The effects of galaxy formation on the matter power spectrum: a challenge for precision cosmology, MNRAS 415 (Aug., 2011) 3649–3665, [arXiv:1104.1174].
  * [2] E. Semboloni, H. Hoekstra, J. Schaye, M. P. van Daalen, and I. G. McCarthy, Quantifying the effect of baryon physics on weak lensing tomography, Monthly Notices of the Royal Astronomical Society 417 (2011), no. 3 2020–2035.
  * [3] J. Harnois-Déraps, L. van Waerbeke, M. Viola, and C. Heymans, Baryons, neutrinos, feedback and weak gravitational lensing, Monthly Notices of the Royal Astronomical Society 450 (2015), no. 2 1212–1223.
  * [4] H.-J. Huang, T. Eifler, R. Mandelbaum, and S. Dodelson, Modelling baryonic physics in future weak lensing surveys, Monthly Notices of the Royal Astronomical Society 488 (2019), no. 2 1652–1678.
  * [5] N. E. Chisari, A. J. Mead, S. Joudaki, P. G. Ferreira, A. Schneider, J. Mohr, T. Tröster, D. Alonso, I. G. McCarthy, S. Martin-Alvarez, et al., Modelling baryonic feedback for survey cosmology, The Open Journal of Astrophysics 2 (2019), no. 1 4.
  * [6] J. Schaye, C. D. Vecchia, C. Booth, R. P. Wiersma, T. Theuns, M. R. Haas, S. Bertone, A. R. Duffy, I. McCarthy, and F. van de Voort, The physics driving the cosmic star formation history, Monthly Notices of the Royal Astronomical Society 402 (2010), no. 3 1536–1560.
  * [7] A. M. Le Brun, I. G. McCarthy, J. Schaye, and T. J. Ponman, Towards a realistic population of simulated galaxy groups and clusters, Monthly Notices of the Royal Astronomical Society 441 (2014), no. 2 1270–1290.
  * [8] I. G. McCarthy, A. Le Brun, J. Schaye, and G. Holder, The thermal sunyaev–zel’dovich effect power spectrum in light of planck, Monthly Notices of the Royal Astronomical Society 440 (2014), no. 4 3645–3657.
  * [9] I. G. Mccarthy, J. Schaye, S. Bird, and A. M. C. Le Brun, The BAHAMAS project: calibrated hydrodynamical simulations for large-scale structure cosmology, MNRAS 465 (2017) 2936–2965.
  * [10] I. G. Mccarthy, S. Bird, J. Schaye, J. Harnois-Deraps, A. S. Font, and L. Van Waerbeke, The BAHAMAS project: the CMB-large-scale structure tension and the roles of massive neutrinos and galaxy formation, MNRAS 476 (2018) 2999–3030.
  * [11] A. Schneider and R. Teyssier, A new method to quantify the effects of baryons on the matter power spectrum, Journal of Cosmology and Astroparticle Physics 2015 (2015), no. 12.
  * [12] A. Schneider, R. Teyssier, J. Stadel, N. E. Chisari, A. M. Brun, A. Amara, and A. Refregier, Quantifying baryon effects on the matter power spectrum and the weak lensing shear correlation, Journal of Cosmology and Astroparticle Physics 2019 (3, 2019).
  * [13] G. Aricò, R. E. Angulo, C. Hernández-Monteagudo, S. Contreras, and M. Zennaro, Simultaneous modelling of matter power spectrum and bispectrum in the presence of baryons, Monthly Notices of the Royal Astronomical Society 503 (2021), no. 3 3596–3609.
  * [14] A. Schneider, N. Stoira, A. Refregier, A. J. Weiss, M. Knabenhans, J. Stadel, and R. Teyssier, Baryonic effects for weak lensing: I. Power spectrum and covariance matrix, Journal of Cosmology and Astroparticle Physics (2019), no. 04.
  * [15] J. F. Navarro, C. S. Frenk, and S. D. M. White, The Structure of Cold Dark Matter Halos, ApJ 462 (May, 1996) 563, [astro-ph/9508025].
  * [16] E. A. Baltz, P. Marshall, and M. Oguri, Analytic models of plausible gravitational lens potentials, J. Cosmology Astropart. Phys. 2009 (Jan., 2009) 015, [arXiv:0705.0682].
  * [17] A. V. Kravtsov, A. A. Vikhlinin, and A. V. Meshcheryakov, Stellar Mass—Halo Mass Relation and Star Formation Efficiency in High-Mass Halos, Astronomy Letters 44 (Jan., 2018) 8–34, [arXiv:1401.7329].
  * [18] A. J. Mead, S. Brieden, T. Tröster, and C. Heymans, hmcode-2020 : improved modelling of non-linear cosmological power spectra with baryonic feedback, Monthly Notices of the Royal Astronomical Society 502 (2021), no. 1 1401–1422.
  * [19] M. G. Abadi, J. F. Navarro, M. Fardal, A. Babul, and M. Steinmetz, Galaxy-induced transformation of dark matter haloes, MNRAS 407 (Sept., 2010) 435–446, [arXiv:0902.2477].
  * [20] B. P. Moster, T. Naab, and S. D. M. White, Galactic star formation and accretion histories from matching galaxies to dark matter haloes, MNRAS 428 (Feb., 2013) 3121–3138, [arXiv:1205.5807].
  * [21] D. Potter, J. Stadel, and R. Teyssier, PKDGRAV3: beyond trillion particle cosmological simulations for the next era of galaxy surveys, Computational Astrophysics and Cosmology 4 (May, 2017) 2, [arXiv:1609.08621].
  * [22] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J. P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J. F. Cardoso, J. Carron, A. Challinor, H. C. Chiang, J. Chluba, L. P. L. Colombo, C. Combet, D. Contreras, B. P. Crill, F. Cuttaia, P. de Bernardis, G. de Zotti, J. Delabrouille, J. M. Delouis, E. Di Valentino, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, Y. Fantaye, M. Farhang, J. Fergusson, R. Fernandez-Cobos, F. Finelli, F. Forastieri, M. Frailis, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gruppuso, J. E. Gudmundsson, J. Hamann, W. Handley, D. Herranz, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, A. Karakci, E. Keihänen, R. Keskitalo, K. Kiiveri, J. Kim, T. S. Kisner, L. Knox, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, J. M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. L. Jeune, P. Lemos, J. Lesgourgues, F. Levrier, A. Lewis, M. Liguori, P. B. Lilje, M. Lilley, V. Lindholm, M. López-Caniego, P. M. Lubin, Y. Z. Ma, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M. A. Miville-Deschênes, D. Molinari, L. Montier, G. Morgante, A. Moss, P. Natoli, H. U. Nørgaard-Nielsen, L. Pagano, D. Paoletti, B. Partridge, G. Patanchon, H. V. Peiris, F. Perrotta, V. Pettorino, F. Piacentini, L. Polastri, G. Polenta, J. L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, E. P. S. Shellard, C. Sirignano, G. Sirri, L. D. Spencer, R. Sunyaev, A. S. Suur-Uski, J. A. Tauber, D. Tavagnacco, M. Tenti, L. Toffolatti, M. Tomasi, T. Trombetti, L. Valenziano, J. Valiviita, B. Van Tent, L. Vibert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. Wehus, M. White, S. D. M. White, A. Zacchei, and A. Zonca, Planck 2018 results. VI. Cosmological parameters, arXiv:1807.06209 (7, 2018).
  * [23] G. Aricò, R. E. Angulo, S. Contreras, L. Ondaro-Mallea, M. Pellejero-Ibañez, and M. Zennaro, The BACCO simulation project: a baryonification emulator with neural networks, Monthly Notices of the Royal Astronomical Society 506 (2021), no. 3 4070–4082.
  * [24] O. Roustant, D. Ginsbourger, and Y. Deville, Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization, Journal of Statistical Software 51 (2012), no. 1 1–55.
  * [25] M. A. Bouhlel, N. Bartoli, A. Otsmane, and J. Morlier, Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction, Structural and Multidisciplinary Optimization 53 (2016), no. 5 935–952.
  * [26] H. Wold, Soft modelling by latent variables: the non-linear iterative partial least squares (nipals) approach, Journal of Applied Probability 12 (1975), no. S1 117–142.
  * [27] M. A. Bouhlel, J. T. Hwang, N. Bartoli, R. Lafage, J. Morlier, and J. R. R. A. Martins, A python surrogate modeling framework with derivatives, Advances in Engineering Software (2019) 102662.
  * [28] R. Bellman, Dynamic programming.  Dover Publications, Mineola, N.Y, 2003.
  * [29] S. Theodoridis and K. Koutroumbas, Pattern Recognition, Fourth Edition.  Academic Press, Inc., USA, 4th ed., 2008.
  * [30] A. Schneider, A. Refregier, S. Grandis, D. Eckert, N. Stoira, T. Kacprzak, M. Knabenhans, J. Stadel, and R. Teyssier, Baryonic effects for weak lensing. Part II. Combination with X-ray data and extended cosmologies, J. Cosmology Astropart. Phys. 2020 (Apr., 2020) 020, [arXiv:1911.08494].
  * [31] V. Springel, R. Pakmor, A. Pillepich, R. Weinberger, D. Nelson, L. Hernquist, M. Vogelsberger, S. Genel, P. Torrey, F. Marinacci, et al., First results from the illustristng simulations: matter and galaxy clustering, Monthly Notices of the Royal Astronomical Society 475 (2018), no. 1 676–698.
  * [32] R. Weinberger, V. Springel, L. Hernquist, A. Pillepich, F. Marinacci, R. Pakmor, D. Nelson, S. Genel, M. Vogelsberger, J. Naiman, et al., Simulating galaxy formation with black hole driven thermal and kinetic feedback, Monthly Notices of the Royal Astronomical Society 465 (2016), no. 3 3291–3308.
  * [33] A. Pillepich, V. Springel, D. Nelson, S. Genel, J. Naiman, R. Pakmor, L. Hernquist, P. Torrey, M. Vogelsberger, R. Weinberger, et al., Simulating galaxy formation with the illustristng model, Monthly Notices of the Royal Astronomical Society 473 (2018), no. 3 4077–4106.
  * [34] A. Pillepich, D. Nelson, L. Hernquist, V. Springel, R. Pakmor, P. Torrey, R. Weinberger, S. Genel, J. P. Naiman, F. Marinacci, et al., First results from the illustristng simulations: the stellar mass content of groups and clusters of galaxies, Monthly Notices of the Royal Astronomical Society 475 (2018), no. 1 648–675.
  * [35] N. E. Chisari, M. L. Richardson, J. Devriendt, Y. Dubois, A. Schneider, A. Le Brun, R. S. Beckmann, S. Peirani, A. Slyz, and C. Pichon, The impact of baryons on the matter power spectrum from the horizon-agn cosmological hydrodynamical simulation, Monthly Notices of the Royal Astronomical Society 480 (2018), no. 3 3962–3977.
  * [36] G. Hinshaw, D. Larson, E. Komatsu, D. N. Spergel, C. L. Bennett, J. Dunkley, M. R. Nolta, M. Halpern, R. S. Hill, N. Odegard, L. Page, K. M. Smith, J. L. Weiland, B. Gold, N. Jarosik, A. Kogut, M. Limon, S. S. Meyer, G. S. Tucker, E. Wollack, and E. L. Wright, Nine-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Parameter Results, The Astrophysical Journal Supplement Series 208 (2013), no. 2 19.
  * [37] D. N. Spergel, R. Bean, O. Doré, M. Nolta, C. Bennett, J. Dunkley, G. Hinshaw, N. e. Jarosik, E. Komatsu, L. Page, et al., Three-year wilkinson microwave anisotropy probe (wmap) observations: implications for cosmology, The Astrophysical Journal Supplement Series 170 (2007), no. 2 377.
  * [38] E. Komatsu, K. M. Smith, J. Dunkley, C. L. Bennett, B. Gold, G. Hinshaw, N. Jarosik, D. Larson, M. R. Nolta, L. Page, D. N. Spergel, M. Halpern, R. S. Hill, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, G. S. Tucker, J. L. Weiland, E. Wollack, and E. L. Wright, Seven-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation, The Astrophysical Journal Supplement Series 192 (2, 2011) 18.
  * [39] Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, N. Bartolo, E. Battaner, R. Battye, K. Benabed, A. Benoît, A. Benoit-Lévy, J. P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, A. Bonaldi, L. Bonavera, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J. F. Cardoso, A. Catalano, A. Challinor, A. Chamballu, R. R. Chary, H. C. Chiang, J. Chluba, P. R. Christensen, S. Church, D. L. Clements, S. Colombi, L. P. L. Colombo, C. Combet, A. Coulais, B. P. Crill, A. Curto, F. Cuttaia, L. Danese, R. D. Davies, R. J. Davis, P. de Bernardis, A. de Rosa, G. de Zotti, J. Delabrouille, F. X. Désert, E. Di Valentino, C. Dickinson, J. M. Diego, K. Dolag, H. Dole, S. Donzelli, O. Doré, M. Douspis, A. Ducout, J. Dunkley, X. Dupac, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, M. Farhang, J. Fergusson, F. Finelli, O. Forni, M. Frailis, A. A. Fraisse, E. Franceschi, A. Frejsel, S. Galeotta, S. Galli, K. Ganga, C. Gauthier, M. Gerbino, T. Ghosh, M. Giard, Y. Giraud-Héraud, E. Giusarma, E. Gjerløw, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gregorio, A. Gruppuso, J. E. Gudmundsson, J. Hamann, F. K. Hansen, D. Hanson, D. L. Harrison, G. Helou, S. Henrot-Versillé, C. Hernández-Monteagudo, D. Herranz, S. R. Hildebrandt, E. Hivon, M. Hobson, W. A. Holmes, A. Hornstrup, W. Hovest, Z. Huang, K. M. Huffenberger, G. Hurier, A. H. Jaffe, T. R. Jaffe, W. C. Jones, M. Juvela, E. Keihänen, R. Keskitalo, T. S. Kisner, R. Kneissl, J. Knoche, L. Knox, M. Kunz, H. Kurki-Suonio, G. Lagache, A. Lähteenmäki, J. M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, J. P. Leahy, R. Leonardi, J. Lesgourgues, F. Levrier, A. Lewis, M. Liguori, P. B. Lilje, M. Linden-Vørnle, M. López-Caniego, P. M. Lubin, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marchini, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Masi, S. Matarrese, P. McGehee, P. R. Meinhold, A. Melchiorri, J. B. Melin, L. Mendes, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M. A. Miville-Deschênes, A. Moneti, L. Montier, G. Morgante, D. Mortlock, A. Moss, D. Munshi, J. A. Murphy, P. Naselsky, F. Nati, P. Natoli, C. B. Netterfield, H. U. Nørgaard-Nielsen, F. Noviello, D. Novikov, I. Novikov, C. A. Oxborrow, F. Paci, L. Pagano, F. Pajot, R. Paladini, D. Paoletti, B. Partridge, F. Pasian, G. Patanchon, T. J. Pearson, O. Perdereau, L. Perotto, F. Perrotta, V. Pettorino, F. Piacentini, M. Piat, E. Pierpaoli, D. Pietrobon, S. Plaszczynski, E. Pointecouteau, G. Polenta, L. Popa, G. W. Pratt, G. Prézeau, S. Prunet, J. L. Puget, J. P. Rachen, W. T. Reach, R. Rebolo, M. Reinecke, M. Remazeilles, C. Renault, A. Renzi, I. Ristorcelli, G. Rocha, C. Rosset, M. Rossetti, G. Roudier, B. Rouillé d’Orfeuil, M. Rowan-Robinson, J. A. Rubiño-Martín, B. Rusholme, N. Said, V. Salvatelli, L. Salvati, M. Sandri, D. Santos, M. Savelainen, G. Savini, D. Scott, M. D. Seiffert, P. Serra, E. P. S. Shellard, L. D. Spencer, M. Spinelli, V. Stolyarov, R. Stompor, R. Sudiwala, R. Sunyaev, D. Sutton, A. S. Suur-Uski, J. F. Sygnet, J. A. Tauber, L. Terenzi, L. Toffolatti, M. Tomasi, M. Tristram, T. Trombetti, M. Tucci, J. Tuovinen, M. Türler, G. Umana, L. Valenziano, J. Valiviita, F. Van Tent, P. Vielva, F. Villa, L. A. Wade, B. D. Wandelt, I. K. Wehus, M. White, S. D. M. White, A. Wilkinson, D. Yvon, A. Zacchei, and A. Zonca, Planck 2015 results. XIII. Cosmological parameters, A&A 594 (Sept., 2016) A13, [arXiv:1502.01589].
  * [40] V. Springel, E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh, Monthly Notices of the Royal Astronomical Society 401 (2010), no. 2 791–851.
  * [41] R. Teyssier, Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES, A&A 385 (Apr., 2002) 337–364, [astro-ph/0111367].
  * [42] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al., Scipy 1.0: fundamental algorithms for scientific computing in python, Nature methods 17 (2020), no. 3 261–272.
  * [43] R. L. Chambers, D. G. Steel, S. Wang, and A. Welsh, Maximum likelihood estimation for sample surveys.  CRC Press, 2012.
  * [44] R. J. Rossi, Mathematical statistics: an introduction to likelihood based inference.  John Wiley & Sons, 2018.
  * [45] G. Aricò, R. E. Angulo, C. Hernández-Monteagudo, S. Contreras, M. Zennaro, M. Pellejero-Ibañez, and Y. Rosas-Guevara, Modelling the large-scale mass density field of the universe as a function of cosmology and baryonic physics, Monthly Notices of the Royal Astronomical Society 495 (2020), no. 4 4800–4819.
  * [46] M. P. van Daalen, I. G. McCarthy, and J. Schaye, Exploring the effects of galaxy formation on matter clustering through a library of simulation power spectra, Monthly Notices of the Royal Astronomical Society 491 (6, 2019) 2424–2446.
  * [47] M. Vogelsberger, S. Genel, V. Springel, P. Torrey, D. Sijacki, D. Xu, G. Snyder, S. Bird, D. Nelson, and L. Hernquist, Properties of galaxies reproduced by a hydrodynamic simulation, Nature 509 (2014), no. 7499 177–182.
  * [48] M. Vogelsberger, S. Genel, V. Springel, P. Torrey, D. Sijacki, D. Xu, G. Snyder, D. Nelson, and L. Hernquist, Introducing the illustris project: simulating the coevolution of dark and visible matter in the universe, Monthly Notices of the Royal Astronomical Society 444 (2014), no. 2 1518–1547.
  * [49] S. Genel, M. Vogelsberger, V. Springel, D. Sijacki, D. Nelson, G. Snyder, V. Rodriguez-Gomez, P. Torrey, and L. Hernquist, Introducing the illustris project: the evolution of galaxy populations across cosmic time, Monthly Notices of the Royal Astronomical Society 445 (2014), no. 1 175–200.
  * [50] M. Haider, D. Steinhauser, M. Vogelsberger, S. Genel, V. Springel, P. Torrey, and L. Hernquist, Large-scale mass distribution in the illustris simulation, Monthly Notices of the Royal Astronomical Society 457 (2016), no. 3 3024–3035.
  * [51] A. Vikhlinin, A. Kravtsov, W. Forman, C. Jones, M. Markevitch, S. Murray, and L. Van Speybroeck, Chandra sample of nearby relaxed galaxy clusters: Mass, gas fraction, and mass-temperature relation, The Astrophysical Journal 640 (2006), no. 2 691.
  * [52] A. H. Gonzalez, S. Sivanandam, A. I. Zabludoff, and D. Zaritsky, Galaxy cluster baryon fractions revisited, The Astrophysical Journal 778 (2013), no. 1 14.
  * [53] A. J. Sanderson, E. O’Sullivan, T. J. Ponman, A. H. Gonzalez, S. Sivanandam, A. I. Zabludoff, and D. Zaritsky, The baryon budget on the galaxy group/cluster boundary, Monthly Notices of the Royal Astronomical Society 429 (2013), no. 4 3288–3304.
  * [54] L. Lovisari, T. Reiprich, and G. Schellenberger, Scaling properties of a complete x-ray selected galaxy group sample, Astronomy & Astrophysics 573 (2015) A118.
  * [55] M. Sereno and S. Ettori, Comparing masses in literature (comalit)–i. bias and scatter in weak lensing and x-ray mass estimates of clusters, Monthly Notices of the Royal Astronomical Society 450 (2015), no. 4 3633–3648.
  * [56] D. Eckert, V. Ghirardini, S. Ettori, E. Rasia, V. Biffi, E. Pointecouteau, M. Rossetti, S. Molendi, F. Vazza, F. Gastaldello, et al., Non-thermal pressure support in x-cop galaxy clusters, Astronomy & Astrophysics 621 (2019) A40.
  * [57] S. Ettori, V. Ghirardini, D. Eckert, E. Pointecouteau, F. Gastaldello, M. Sereno, M. Gaspari, S. Ghizzardi, M. Roncarelli, and M. Rossetti, Hydrostatic mass profiles in x-cop galaxy clusters, Astronomy & Astrophysics 621 (2019) A39.
  * [58] G. Hurier and R. E. Angulo, Measuring the hydrostatic mass bias in galaxy clusters by combining Sunyaev-Zel’dovich and CMB lensing data, Astronomy and Astrophysics 610 (11, 2017).
  * [59] M. Douspis, L. Salvati, and N. Aghanim, On the tension between large scale structures and cosmic microwave background, in 2nd World Summit: Exploring the Dark Side of the Universe, vol. 335, p. 037, SISSA Medialab, 2018\.
  * [60] H. Akaike, A new look at the statistical model identification, IEEE transactions on automatic control 19 (1974), no. 6 716–723.
  * [61] K. P. Burnham and D. R. Anderson, A practical information-theoretic approach, Model selection and multimodel inference 2 (2002).
  * [62] G. Schwarz, Estimating the dimension of a model, The annals of statistics (1978) 461–464.
  * [63] I. Chiu, J. Mohr, M. McDonald, S. Bocquet, S. Desai, M. Klein, H. Israel, M. Ashby, A. Stanford, B. Benson, et al., Baryon content in a sample of 91 galaxy clusters selected by the south pole telescope at 0.2< z< 1.25, Monthly Notices of the Royal Astronomical Society 478 (2018), no. 3 3072–3099.
  * [64] J. Sanders, A. Fabian, H. Russell, and S. Walker, Hydrostatic chandra x-ray analysis of spt-selected galaxy clusters–i. evolution of profiles and core properties, Monthly Notices of the Royal Astronomical Society 474 (2018), no. 1 1065–1098.
  * [65] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: the mcmc hammer, Publications of the Astronomical Society of the Pacific 125 (2013), no. 925 306.
  * [66] S. Van Der Walt, S. C. Colbert, and G. Varoquaux, The numpy array: a structure for efficient numerical computation, Computing in science & engineering 13 (2011), no. 2 22–30.
  * [67] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in science & engineering 9 (2007), no. 03 90–95.
  * [68] D. Foreman-Mackey, corner.py: Scatterplot matrices in python, The Journal of Open Source Software 1 (jun, 2016) 24.
  * [69] M. Feurer and F. Hutter, Hyperparameter optimization, in Automated machine learning, pp. 3–33.  Springer, Cham, 2019.
  * [70] P. Stoica and Y. Selen, Model-order selection: a review of information criterion rules, IEEE Signal Processing Magazine 21 (2004), no. 4 36–47.