# Determining the Baryon Impact on the Matter Power Spectrum with Galaxy
Clusters

Sebastian Grandis,1 Giovanni Aricò,2 Aurel Schneider,2 Laila Linke1
1Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr.
25/8, 6020 Innsbruck, Austria
2Institute for Computational Science, University of Zurich,
Winterthurerstrasse 190, 8057 Zurich, Switzerland
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The redistribution of baryonic matter in massive halos through processes like
active galactic nuclei feedback and star formation leads to a suppression of
the matter power spectrum on small scales. This redistribution can be measured
empirically via the gas and stellar mass fractions in galaxy clusters, and
leaves imprints on their electron density profiles. We constrain two semi-
analytical baryon correction models with a compilation of recent Bayesian
population studies of galaxy groups and clusters sampling a mass range above
$\sim 3\times 10^{13}$ $M_{\odot}$, and with cluster gas density profiles
derived from deep, high-resolution X-ray observations. We are able to fit all
the considered observational data, but highlight some anomalies in the
observations. The constraints allow us to place precise, physically informed
priors on the matter power spectrum suppression. At a scale of $k=1h$ Mpc-1 we
find a suppression of $0.042^{+0.012}_{-0.014}$ ($0.049^{+0.016}_{-0.012}$),
while at $k=3h$ Mpc-1 we find $0.184^{+0.026}_{-0.031}$
($0.179^{+0.018}_{-0.020}$), depending on the model used. We also predict at
97.5 percent credibility, that at scales $k<0.37h$ Mpc-1 baryon feedback
impacts the matter power less than $1\%$. This puts into question if baryon
feedback is the driving factor for the discrepancy between cosmic shear and
primary CMB results. We independently confirm results on this suppression from
small-scale cosmic shear studies, while we exclude some hydro-dynamical
simulations with too strong and too weak baryonic feedback. Our empirical
prediction of the power spectrum suppression shows that studies of galaxy
groups and clusters will be instrumental in unlocking the cosmological
constraining power of future cosmic shear experiments like Euclid and Rubin-
LSST.

###### keywords:

large-scale structure of Universe – galaxies: clusters: general – methods:
data analysis

††pubyear: 2015††pagerange: Determining the Baryon Impact on the Matter Power
Spectrum with Galaxy Clusters–C

##  1 Introduction

Cosmological inference on the matter distribution of the Universe from future
surveys like Euclid111http://www.euclid-ec.org/, LSST222https://www.lsst.org/,
or Roman333https://roman.gsfc.nasa.gov/ will be limited by our knowledge of
non-gravitational effects on the matter power spectrum (for reviews, see
Chisari et al., 2019; Eckert et al., 2021, and reference therein).
Specifically, on physical scales below $\sim 10$ Mpc, active galactic nuclei
are known to lead to a significant redistribution of a fraction of the
baryonic matter, while on even smaller scales star formation allows another
fraction of baryons to condensate into galaxies inhabiting the halo centers.
Both these redistribution mechanisms lead to a gravitational back reaction on
the collisionless dark matter, further altering the matter distribution. These
effects are collectively referred to as baryon feedback, and are degenerate
with interesting cosmological signals like neutrino masses, dark energy
equation of state modifications, modified gravity signatures or non cold dark
matter components (e.g. Harnois-Déraps et al., 2015; Schneider et al., 2020).

While hydro-dynamical simulations are in principle able to accurately predict
baryon feedback effects on cosmology, practical concerns like run-time and
resolution limits require the tuning of several recipes that attempt the
modelling of baryon feedback by summarizing sub resolution processes (Schaye
et al., 2015b; McCarthy et al., 2017; Pillepich et al., 2018). While providing
physically self-consistent solutions, hydro-dynamical simulations are thus
only able to present point estimates for individual, best guess feedback
models. They lack the predictive power and flexibility to be employed in
cosmological inference tasks, as they probe only dozens of feedback models
(Mead et al., 2015; Chisari et al., 2018; van Daalen et al., 2020), though
noticeable progress has been made by the recent CAMELS (Villaescusa-Navarro et
al., 2021), MillenniumTNG (Pakmor et al., 2022)444https://www.mtng-
project.org/, ANTILLES (Salcido et al., 2023), and
FLAMINGO555https://flamingo.strw.leidenuniv.nl/ (Schaye et al., 2023; Kugel et
al., 2023) projects.

Filling this gap are so called Baryon Correction Models (hereafter BCM,
originally proposed by Schneider & Teyssier, 2015). Starting from assumptions
on the fraction and distribution of baryonic matter in and around halos, these
models alter gravity-only simulations in a semi-analytical way, to mimic the
impact of baryon feedback. Recently, Aricò et al. (2020); Aricò et al. (2021a)
and Giri & Schneider (2021) have shown that their respective BCMs are able
simultaneously to fit the stellar and gas components of halos over several
orders of magnitude in mass, as well as the power and bispectrum in different
hydro-dynamical simulations with a handful of physically motivated parameters.
The resulting gas and stellar fractions as well as suppression of the matter
power and bispectrum are easy to compute. This provides a significant benefit
compared to feedback prescriptions that only parameterize the matter power
spectrum suppression, as for instance presented by Mead et al. (2021), and
offer no predictive power on other quantities.

First attempts to constrain these BCMs have been undertaken with a variety of
data sets. Schneider et al. (2022) used Kilo Degree Survey cosmic shear,
Atacama Cosmology Telescope kinematic Sunyaev-Zeldovich profiles and a
compilation of galaxy cluster and group measurements to simultaneously fit for
the cosmological parameters and the BCM. On the cluster and group side that
work was limited by the use of individual hydrostatic mass estimators from
compilations of X-ray observed objects. Such compilations suffer from a lack
of selection effects and mass calibration modelling, as they assume that the
compiled list of clusters and groups is a fair sample of the underlying halo
population, and that their halo masses can be calibrated by specifying on
single parameter, the hydrostatic mass bias. Chen et al. (2022b) and Aricò et
al. (2023) have derived first constraints on BCMs using Dark Energy Survey
year 3 cosmic shear on small scale, and on all scales, respectively. Both
works find weak constraints on one of the seven parameters of the BCM.

Parallel to these studies, major progress has been achieved in the
observational study of galaxy clusters, which inhabit massive halos. Bayesian
Population models simultaneously describe cluster selection and mass
calibration, preferentially via the weak gravitational lensing (WL) signal
these massive objects impress on background galaxies (Bocquet et al., 2015;
Mantz et al., 2015; Sereno & Ettori, 2017; Dietrich et al., 2019; Chiu et al.,
2022). The robustness of these modelling techniques is demonstrated by the
fact that the number counts of clusters are an independent, competitive
cosmological probe (Mantz et al., 2015; Bocquet et al., 2019; Chiu et al.,
2023). This indicates that the mass distribution of these samples can be
reliably reconstructed in the aforementioned Bayesian frameworks.
Additionally, multi-wavelength studies of samples of galaxy clusters measure
the gas and stellar mass in massive halos, and its mass trend (see for
instance Mantz et al., 2016; Chiu et al., 2018; Akino et al., 2022; Chiu et
al., 2022). Crucially, these studies forego the use of hydrostatic masses in
favor of the more accurate WL signal. This lifts the uncertainties on sample
selection and mass calibration that plague compilation of X-ray observed
clusters, like the one recently presented by Kugel et al. (2023).

In this work, we combine recent weak lensing informed gas and stellar mass
fraction measurements of massive halos (Mantz et al., 2016; Chiu et al., 2018;
Akino et al., 2022; Chiu et al., 2022) to inform the BCM and predict the
matter power spectrum suppression due to baryon feedback. The introduction of
weak lensing mass calibration allows us to avoid using hydrostatic masses,
significantly improving our mass accuracy. We also include measurements of the
gas profile from deep, high-resolution X-ray observations (Ghirardini et al.,
2019) to narrow down the shape of the gas density profiles of massive halos.

Halo masses in this work are reported as spherical over density masses
$M_{\Delta\rm{c}}$, with $\Delta=500,200$. This means that they are defined
via the radius $R_{\Delta\rm{c}}$ enclosing a sphere of average density
$\Delta$ times the critical density of the Universe at that redshift, i.e.
$M_{\Delta\rm{c}}=\frac{4\pi}{3}\Delta\rho_{\text{crit}}(z)R_{\Delta\rm{c}}^{3}$.
Throughout this work, we use a flat $\Lambda$CDM cosmology with parameters
$\Omega_{\text{M}}=0.315$, $\sigma_{8}=0.83$, $n_{\text{S}}=0.96$ and $h=0.67$
as reference cosmology.

##  2 Data

Figure 1: Approximate mass distributions of the galaxy cluster and group samples used. These distributions are not used in our inference and are presented for visualisation purposes only. Our compilation of studies samples the mass range $M_{500\text{c}}\gtrapprox 3\times 10^{13}$ $M_{\odot}$. Sources: eFEDS-HSC Chiu et al. (2022, Tab. C1); SPT Bocquet et al. (2019, with $\xi>6.8$); HSC-XXL redshifts from Umetsu et al. (2020, Tab. 2) together the mass distribution parameters from Akino et al. (2022); WtG Mantz et al. (2016, Tab. 2 and Eq. 1). Table 1:  Baryonic fractions and mass trends of weak lensing calibrated cluster and group population studies compiled by this work. | comp. |  $M_{\text{piv}}$ [$M_{\odot}$] | fraction | slope
---|---|---|---|---
SPT | ICM | 4.80e+14 | $0.12\pm 0.01$ | $0.33\pm 0.09$
SPT | stars | 4.80e+14 | $0.008\pm 0.001$ | $-0.20\pm 0.12$
eFEDS | ICM | 2.00e+14 | $0.05\pm 0.01$ | $0.19\pm 0.11$
HSC-XXL | ICM | 1.00e+14 | $0.08\pm 0.01$ | $0.23\pm 0.12$
HSC-XXL | stars | 1.00e+14 | $0.021\pm 0.003$ | $-0.20\pm 0.11$
WtG | ICM | 1.00e+15 | $0.12\pm 0.00$ |

###  2.1 Weak lensing calibrated mass fractions

In this work we use measurements of the relation between the gas / stellar
mass and the halo mass in galaxy clusters to constrain Baryon Correction
Models. We use the following compilation of stellar and gas mass fractions
extracted from Bayesian Population models, summarized also in Table 1:

  * •

The analysis of 139 ROSAT selected clusters, followed up with Chandra imaging
analysed by Mantz et al. (2016), 26 of which have dedicated weak lensing
measurements. This analysis, called Weighing the Giants (WtG), reports a hot
gas fraction $f_{\text{ICM}}=0.125\pm 0.005$ at pivot mass
$M_{\text{piv}}=10^{15}$ $M_{\odot}$666Note that all these observational
results measure baryonic fractions within $r_{\text{500c}}$, and w.r.t. halo
masses $M_{\text{500c}}$.. The mass information comes from a combination of
cluster abundance measurements and weak lensing. We discard the slope
measurement of this work, as it has implausibly tight error bars that would
otherwise dominate our fitting.

  * •

The study of 91 South Pole Telescope (SPT) selected clusters, whose Chandra,
DES, Wise and Spitzer follow-up has been studied by Chiu et al. (2018),
reporting a hot gas fraction $f_{\text{ICM}}=0.119\pm 0.013$ and a mass trend
of the gas mass $B_{\text{ICM}}=1.33\pm 0.09$ at $M_{\text{piv}}=4.8\times
10^{14}$ $M_{\odot}$, as well as stellar mass fraction $10^{3}f_{\star}=8.3\pm
0.6$ and mass trend of the stellar mass $B_{\star}=0.80\pm 0.13$. The mass
information in this analysis comes from fitting the cluster abundance
measurements (de Haan et al., 2016). The resulting mass calibration has been
confirmed via weak lensing (Dietrich et al., 2019; Stern et al., 2019; Bocquet
et al., 2019; Schrabback et al., 2021) and dynamical analysis (Capasso et al.,
2019).

  * •

The analysis of 136 galaxy clusters and groups by Akino et al. (2022), which
we call HSC-XXL, reporting $f_{\text{ICM}}=0.075\pm 0.008$ and
$B_{\text{ICM}}=1.23\pm 0.12$ at $M_{\text{piv}}=10^{14}$ $M_{\odot}$, as well
as $10^{3}f_{\star}=21.2\pm 2.8$ and $B_{\star}=0.80\pm 0.11$. These objects
were selected in the XMM-Newton XXL survey, and have HSC-SSP weak lensing
measurements, XMM-Newton gas mass measurements and stellar mass measurements
for HSC and SDSS photometry.

  * •

The study by Chiu et al. (2022) of 434 clusters and groups selected in the
eROSITA Final Equatorial Depth Survey (eFEDS), 313 of which have HSC-SSP weak
lensing measurements, reporting $f_{\text{ICM}}=0.0540\pm 0.0065$ and
$B_{\text{ICM}}=1.19\pm 0.11$ at $M_{\text{piv}}=2\times 10^{14}$ $M_{\odot}$.

The selection of these samples is performed via their extended Bremsstrahlung
emission in the X-rays (HSC-XXL, eFEDS, WtG), or via the Sunyaev-Zel’dovich
effect (SZe, see Carlstrom et al. (2002) for a review) in (sub-)millimeter
observation (SPT), i.e via inverse Compton scattering signatures in cosmic
microwave background observations, and confirmed in optical wavelength. These
methods provide high purity samples of massive halos, albeit with different
redshift trends in the limiting mass. The information on the mass scale of
these samples comes from weak lensing measurements in dedicated deep
photometric follow-up (WtG) or deep photometric survey data (HSC-XXL, eFEDS),
and in some cases from abundance matching to the theoretical halo mass
function (SPT, WtG). Their gas mass is measured in high angular resolution
X-ray follow-ups (SPT, WtG) or X-ray survey data (HSC-XXL, eFEDS). The stellar
masses are determined with optical and near-infrared observations of the
cluster and group member galaxies. We refer the interested reader to the
respective papers for more information.

Of interest to this work is that the different samples have minimal to no
overlap in actual objects. We can thus treat the different measurements as
mutually independent. Also, their mass distributions probe significantly
different mass ranges, as seen in Fig. 1. These mass distributions are derived
outcomes of the Bayesian Population analyses described below (cf. section
3.1), and are not directly used in this analysis. We instead employ only the
reconstructed observable halo mass relations.

Table 2:  Mean electron density profile $n_{e}$, and its error $\delta n_{e}$, as a function of radius, from the X-COP sample at median redshift $z_{\text{med}}=0.064$ and median hydro static mass $M^{\text{hydro}}_{\text{med},500c}=3.79\times 10^{14}\,h^{-1}M_{\odot}$, adapted from Ghirardini et al. (2019, Tab. 2). $r$ [Mpc$h^{-1}$] |  $n_{e}$ [$10^{-3}$ cm-3] |  $\delta n_{e}$ [$10^{-3}$ cm-3]
---|---|---
0.034 | 6.159 | 0.605
0.085 | 3.371 | 0.380
0.145 | 2.240 | 0.196
0.221 | 1.395 | 0.085
0.327 | 0.886 | 0.031
0.502 | 0.445 | 0.011
0.791 | 0.179 | 0.006
1.339 | 0.053 | 0.004

###  2.2 Electron density profiles

We complement the measurements of the relation between gas / stellar mass and
halo mass for large populations of clusters with information on the gas
density profile of 13 high mass, low redshift clusters observed with deep and
high resolution dedicated X-ray imaging, carried out by the XMM-Newton Cluster
Outskirts Project (X-COP, Eckert et al., 2017). Analysis of their X-ray
surface brightness and spectral information by Ghirardini et al. (2019)
resulted in a measurement of the hydro-static mass of these objects, with
median hydro static mass $M^{\text{hydro}}_{\text{med},500c}=3.79\times
10^{14}\,h^{-1}M_{\odot}$, as well as their electron density profile
$n_{e}\left(r/r^{\text{hydro}}_{500c}\right)$. We utilize the results from the
piece wise linear interpolation of the electron density profiles given in
Ghirardini et al. (2019, Tab. 2) and the median redshift of the X-COP sample,
$z_{\text{med}}=0.064$, to generate a measurement of the electron density
profile $n_{e}$ in units of cm-3 at 8 radii $r$ in units of Mpc$h^{-1}$. The
resulting data vector is reported in Tab. 2. The resulting electron density
profile is shown in Fig. 3 as black points. Note that we scaled the electron
density profile by the radius $r$ to reduce the dynamical range in the plot.
Many more high resolution X-ray measurements of the electron density of
clusters would, in principle, be available (e.g. Croston et al., 2008;
McDonald et al., 2013; Sanders et al., 2018; Bulbul et al., 2019). Adapting
their results to our fitting procedure would, however, far exceed the scope of
this work.

##  3 Method

We shall describe in the following how the multi-wavelength observations of
galaxy groups and clusters are analysed to extract the gas and stellar mass
fractions as a function of halo mass. We then describe how we use the BCM to
fit these fractions, as well as the electron density profile of halos.
Finally, we outline our fitting procedure and the adopted priors.

###  3.1 Bayesian Population Models

The galaxy cluster and group studies used in this work, and listed in Section
2, infer the relation between the gas / stellar mass $M_{\text{ICM},\star}$
enclosed in $R_{500\text{c}}$ and the halo mass $M_{500\text{c}}$ via Bayesian
Population models. These models are underpinned by a statistical prescription
that describes the generation of the cluster catalog data starting from the
parent distribution in mass and redshift space, $P(M,z)$ (often the halo mass
function times the cosmological volume). A stochastical mapping from halo mass
$M$ and redshift $z$ to intrinsic (noise free) observables $\mathcal{O}$ is
then applied, $P(\mathcal{O}|M,z)$. Motivated by empirical evidence and
simulation results (Angulo et al., 2012), this mapping is modelled as a
multivariate log-normal distribution in the observables. The means are power
laws in mass and redshift, called _mass observable relations_. The simple
power law relation is justified by the fact that clusters deviate only to
second order from self similar behaviour characteristic of gravity and
adiabatic thermodynamics777Mulroy et al. (2019, sec. 3) presents an excellent
review of the original derivation by Kaiser (1986).. The multivariate scatter
around these mean relations is an expression of the heterogeneity of the
cluster population at a given mass and redshift (see, for instance, Farahi et
al., 2019). Measurement noise in the observables is modelled via a
stochastical mapping between intrinsic (noise-free) observables and measured
observables, $P(\hat{\mathcal{O}}|\mathcal{O},z)$. This last mapping is
determined directly from the noise properties of the observations, or via
image simulations.

The distribution of observed cluster properties $P(\hat{\mathcal{O}},z)$ is
then obtained by marginalising of mass and intrinsic properties

$P(\hat{\mathcal{O}},z)=\int\text{d}M\,P(M,z)\int\text{d}\mathcal{O}P(\mathcal{O}|M,z)P(\hat{\mathcal{O}}|\mathcal{O},z).$
(1)

This distribution is normalized to account for the selection criteria imposed
in the sample’s selection. The likelihood of a sample of clusters
$\big{\\{}\hat{\mathcal{O}}_{i},z_{i}\big{\\}}$ results from evaluating the
probability of this sample given the constructed distribution

$\ln\mathcal{L}=\sum_{i}\ln P(\hat{\mathcal{O}}_{i},z_{i}),$ (2)

as in Bayesian inference, the likelihood is simply the probability of the data
(here the catalog) given the model (the distribution of observed cluster
properties). This analysis approach ensures that selection biases and mass
calibration uncertainties are correctly accounted for.

In the analyses we use for this work, the measured observables include the gas
mass $\hat{M}_{\text{ICM}}$, the weak lensing signal in the form of a shear
profile, or a best fit mass, and, where applicable, the stellar mass
$\hat{M}_{\star}$, among other observables. For further details on each
cluster study we refer the reader to the relative papers and references
therein, noting that many details of the implementation, as well as the
individual notations used to describe the Bayesian population model might
differ significantly between the different works. Our notations follow most
closely (Chiu et al., 2022, sec. 5).

Among the model parameters of the population likelihood, there are the
amplitude of the gas / stellar mass – halo mass relation
$A_{\text{ICM},\star}$, and the mass slope $B_{\text{ICM},\star}$ of that
relation. They define scaling relations that follow power laws in mass around
a fixed pivot mass $M_{\text{piv}}$, chosen close to the median mass of the
sample, reading

$\ln\Bigg{(}\frac{M_{\text{ICM},\star}}{M_{\text{ICM},\star}^{\text{piv}}}\Bigg{)}=A_{\text{ICM},\star}+B_{\text{ICM},\star}\ln\Big{(}\frac{M_{500\text{c}}}{M^{\text{piv}}}\Big{)}.$
(3)

Here, we marginalize over any redshift dependencies, that are usually also
studied, as they have all empirically been shown to be consistent with zero.
Crucially, the statistically and systematic uncertainties incurring in the gas
/ stellar content and total mass measurement process, as well as in the
selection of the cluster sample, are accounted for and distilled into
posteriors on the scaling relation parameters
$(A_{\text{ICM},\star},B_{\text{ICM},\star})$. We therefore base our analysis
on these summary statistics, together with their reported uncertainties.

###  3.2 Baryon Correction Model

The Baryon Correction Model proposed by Schneider & Teyssier (2015) provides a
semi-analytical framework to modify N-body, gravity-only simulations
accounting for the effects of gas, stars and feedback. The framework is based
on small, radial shifts of simulation particles around halo centres. In
particular, the gravity-only profile of a given halo is mapped to the sum of a
dark matter, a gas and a stellar matter profile,

$\rho(r)\mapsto\rho_{\text{DM}}(r)+\rho_{\text{ICM}}(r)+\rho_{\star}(r).$ (4)

The model parameters determine the amplitude and shape of these profiles for
each input halo mass. For each set of model parameters, the particles in the
gravity-only simulation are displaced in order to obtain the new,
_baryonified_ halo profiles. Quasi-adiabatic relaxation of the corrected
profile is used to account for the gravitational back-reaction of the
displaced baryons. In this work, we use two BCM models: the "S19 model"
proposed by Schneider & Teyssier (2015); Schneider et al. (2019); Giri &
Schneider (2021), and the bacco-model proposed by Aricò et al. (2020); Aricò
et al. (2021a). In spirit, the two models are similar, and we shall only
quickly outline their differences in the following.

####  3.2.1 S19 model

The BCM from Schneider et al. (2019, hereafter S19) and Giri & Schneider
(2021) provides a direct parametrisation of the density profiles for the gas
and the stellar components. For the gas component, a power-law profile with an
additional core in the centre and steep truncation in the outskirts is
assumed. The functional form is given by

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

with a mass dependent inner slope

$\beta(M_{\text{vir}})=\frac{3\left(M_{\text{vir}}/M_{\text{c}}\right)^{\mu}}{1+\left(M_{\text{vir}}/M_{\text{c}}\right)^{\mu}},$
(6)

where $M_{\text{vir}}$ and $r_{\text{vir}}$ are the virial mass and the virial
radius, and $\Omega_{\text{b,M}}$ are the cosmic baryon/matter abundance.

The stellar profile is parametrised as a central stellar component ("cga") and
a satellite component ("sat"), the latter following the dark matter profile
due to its collisionless nature. The fractions of the two components follow a
power law in mass

$f_{\text{cga,sat}}=0.055\left(\frac{2\times
10^{11}M_{\odot}/h}{M_{\text{vir}}}\right)^{\eta_{\text{cga,sat}}},$ (7)

with slopes $\eta_{\text{sat}}=\eta$ and $\eta_{\text{cga}}=\eta+\delta\eta$.
Upon specifying the baryonic components, the collisionless profiles of the
dark matter and satellite galaxies are subjected to adiabatic relaxation. This
leaves us with 7 free parameters for the S19-model:
$(\log_{10}M_{\text{c}}h/M_{\odot},\mu,\theta_{\text{ej}},\gamma,\delta)$ for
the gas profile and $(\eta,\delta\eta)$ for the stellar component.

Integration of the relative profiles to $r_{500\text{c}}$ allows us to compute
the total mass at that over-density $M_{500\text{c}}$, as well as the enclosed
gas and stellar masses $M_{\text{ICM},\star}$, from which we can readily
derive the stellar and gas fractions.

####  3.2.2 bacco model

The bacco simulation projects (Contreras et al., 2020; Angulo et al., 2021;
Zennaro et al., 2023) proposed a BCM that follows a slightly different
approach to define baryonic profiles. For each gravity-only halo of mass
$M_{200\text{c}}$, a fraction of its mass $f_{\text{ICM}}(M_{200\text{c}})$ is
redistributed into an empirically motivated gas distribution

$\rho_{\text{ICM}}(r)\propto\begin{cases}\left(1+\frac{r}{r_{\text{inn}}}\right)^{-\beta_{\text{inn}}}\left(1+\left(\frac{r}{r_{\text{out}}}\right)^{2}\right)^{-2},&\text{
if }r<r_{\text{out}}\\\ \rho_{\text{NFW}}(r),&\text{otherwise,}\end{cases}$
(8)

with an inner scale $r_{\text{inn}}=\theta_{\text{inn}}r_{\text{200c}}$, the
inner slope of the gas profiles
$\beta_{\text{inn}}=3-\left(\frac{M_{\text{inn}}}{M_{\text{200c}}}\right)^{\mu_{\text{inn}}}$,
with $\mu_{\text{inn}}=0.31$, and the outer scale
$r_{\text{out}}=\theta_{\text{out}}^{-1}r_{\text{200c}}$888This definition
deviates from the notation in previous works (Aricò et al., 2020; Aricò et
al., 2021a, b). It reflects the actual implementation in the code at the time
of writing. Changing $\theta_{\text{out}}^{-1}\mapsto\theta_{\text{out}}$ is
planned for future releases.. $\rho_{\text{NFW}}(r)$ is the dark matter
density profile proposed by Navarro et al. (1996, hereafter NFW) with the
concentration mass relation fitted on the respective simulation. The
normalisations of the two components are adjusted such that the transition is
continuous and the integrated profiles attain the gas fraction
$f_{\text{ICM}}(M_{200\text{c}})$ at the corresponding radius.

Another fraction $f_{\star}(M_{200\text{c}})$ is assumed to be stellar matter
and distributed between a component mimicking the central galaxy, and another
representing the satellite galaxies. This provides the prescription for the
stellar matter density profile $\rho_{\star}(r)$. Thus, this model defines the
fractions at $r_{200\text{c}}$ first, then defines the shapes of baryonic
profiles.

Aricò et al. (2020); Aricò et al. (2021a) refined this model by allowing for a
late time re-accreted gas component, and an ejected gas component.
Specifically, they model the halo gas fraction within $r_{200\text{c}}$ as

$f_{\text{ICM,
200c}}=\Big{(}\frac{\Omega_{\text{b}}}{\Omega_{\text{M}}}-f_{\star}(M_{200\text{c}})\Big{)}\Bigg{(}1+\Big{(}\frac{M_{\text{c}}}{M_{200\text{c}}}\Big{)}^{\beta}\Bigg{)}^{-1},$
(9)

where $M_{\text{c}}$ and $\beta$ are free parameters of the bacco-model. Note
that both models define a parameter $M_{\text{c}}$, which however is not the
same between the models. The stellar fraction is modulated by the parameter
$M_{1}$, which is the characteristic halo mass with a stellar mass fraction
0.023. Fitting for $M_{1}$ is thus equivalent to finding the halo mass, whose
stellar mass fraction is 0.023. The stellar-halo mass relation is modelled
with the sub-halo abundance matching proposed by Behroozi et al. (2013), and
best-fitting parameters found by Kravtsov et al. (2018). The free parameters
of the gas profile are the characteristic radii $\theta_{\text{inn}}$ and
$\theta_{\text{out}}$, as well as $M_{\text{inn}}$, which models the inner
slope of the gas profile. Finally, the parameter $\eta$ describes the position
of the ejected gas cut-off. This gives a total of seven free parameters for
the bacco model.

Contrary to the S19-model, the bacco model only displaces particles within a
given radius $r_{\text{p}}$ of halos for computational ease. Deviating from
previous works, which fixed $r_{\text{p}}=r_{\rm 200c}$, we set
$r_{\text{p}}=r_{\text{out}}$, and discuss the resulting differences in
Appendix B. We justify this choice with the argument that beyond
$r_{\text{out}}$ the hot gas is assumed to perfectly trace the dark matter.
The only component which deviates from the NFW profile beyond $r_{\text{out}}$
is the ejected gas, which is going to be traced with particles displaced from
inside $r_{\text{out}}$. We note that previous works always assumed
$r_{\text{out}}\leq r_{\rm 200c}$, thus this condition was always satisfied in
previous work.

####  3.2.3 Link to scaling relations

The BCM is directly linked to the parameters of the scaling relation via a
Taylor expansion of the gas / stellar mass fraction enclosed in
$R_{500\text{c}}$ around the pivot mass, yielding

$f_{\text{ICM},\star}^{500\text{c}}\Big{|}_{M_{500\text{c}}=M_{\text{piv}}}=e^{A_{\text{ICM},\star}}\frac{M_{\text{ICM},\star}^{\text{piv}}}{M^{\text{piv}}},\text{
and}$ (10)

$\frac{\text{d}\ln f_{\text{ICM},\star}^{500c}}{\text{d}\ln
M_{500\text{c}}}\Big{|}_{M_{500\text{c}}=M_{\text{piv}}}=B_{\text{ICM},\star}-1.$
(11)

This expression provides the link between the baryonification models and the
scaling relation studies performed with Bayesian Population models. This
constitutes a crucial methodological advancement with respect to the use of
compilations of individual clusters that do not account for mass accuracy and
selection modelling. WL calibrated Bayesian population models of cluster
samples are the tool of choice to account for selection effects and mass
accuracy. Compared to un-binned compilations of heterogeneously selected
clusters, like the one used by Kugel et al. (2023), they are able to
accurately reconstruct the underlying true halo mass distribution of the
cluster samples, as well as its relation to observable quantities like the gas
mass. This methodological improvement is corroborated by the fact that WL
calibrated cluster population models are able to derive competitive
cosmological constraints from the abundance clusters (Mantz et al., 2015;
Bocquet et al., 2019; Chiu et al., 2023).

Practically speaking, we also need to convert between the different mass
definitions used in X-ray studies of clusters ($M_{500\text{c}}$) and the BCM
models ($M_{200\text{c}}$). In the S19-model the baryon fractions at
$r_{500\text{c}}$ are directly predicted by integration of the respective
profiles. For the bacco-model the conversion from $M_{500\text{c}}$ to
$M_{200\text{c}}$ is performed using the concentration mass relation by
Ragagnin et al. (2021) at its pivot cosmology, while the enclosed ICM mass is
corrected using the gas distribution of the BCM, as further described in
Appendix A.

####  3.2.4 Predicting electron density profiles

The ICM is, to good approximation, a fully ionized gas, such that the electron
density profiles $n_{\text{e}}(r)$ can readily be estimated from the gas
density as

$n_{\text{e}}(r)=\frac{\rho_{\text{ICM}}(r)}{\mu_{\text{e}}m_{\text{p}}},$
(12)

with the mean molecular weight of electrons $\mu_{\mathrm{e}}=1.17$ (following
Bulbul et al., 2019), and the proton mass $m_{\text{p}}$. We can thus readily
transform the ICM profiles assumed by the BCM into predictions for the
electron density profiles.

One operational note is that studies of X-ray electron density profiles of
massive clusters have to date not been carried out together with accurate WL
mass calibration, as described for the gas and stellar fractions. We are
therefore forced to rely on the hydro static mass reported by Ghirardini et
al. (2019), and transform this to a halo mass using the hydrostatic mass bias
$b_{\text{HS}}$, with prior
$p(b_{\text{HS}})=\mathcal{N}(b_{\text{HS}}|0.26,0.07^{2})$ (Hurier & Angulo,
2018).

###  3.3 Fitting Procedure

Table 3:  Prior choices of this analysis. All priors are chosen uniform within the bounds that we report, except for the hydrostatic mass bias, where we use a Gaussian prior and indicate the mean and variance of the distribution. Note that both BCMs use the parameter $\log_{10}M_{\text{c}}h/M_{\odot}$, but it has different meanings in the two models, as described in section 3.2.  AS-model | bacco
---|---
$\log_{10}M_{\text{c}}h/M_{\odot}$ | $(11,\,15)$ | $\log_{10}M_{\text{c}}h/M_{\odot}$ | $(12,\,15)$
$\mu$ | $(0,\,2)$ | $\log_{10}\beta$ | $(-1,\,0.7)$
$\theta_{\text{ej}}$ | $(2,\,8)$ | $\log_{10}M_{1}h/M_{\odot}$ | $(11.5,\,15)$
$\gamma$ | $(1,\,4)$ | $\log_{10}\eta$ | $(-0.7,\,0.7)$
$\delta$ | $(3,\,11)$ | $\log_{10}\theta_{\text{out}}$ | $(-1,\,0.5)$
$\eta$ | $(0.05,\,0.4)$ | $\log_{10}M_{\text{inn}}h/M_{\odot}$ | $(5,\,15)$
$\delta\eta$ | $(0.05,\,0.4)$ | $\log_{10}\theta_{\text{inn}}$ | $(-2,\,0)$
cosmological parameters
$\Omega_{\text{M}}=0.315$ | $h=0.67$ | $\Omega_{\text{b}}\in(0.0473,0.0535)$
hydro static mass bias
$b_{\text{HS}}\sim\mathcal{N}(0.26,0.07^{2})$

Figure 2: Top: Prediction for the gas (magenta) and stellar (orange) mass
fraction together with the fitted data points from the cluster analyses
considered. We show the 1 sigma region as filled and represent also the 2
sigma region with faded lines. Bottom: Prediction for the mass slope of the
gas (magenta) and stellar (orange) mass fraction together with the fitted data
points from the cluster analyses considered. The data point for the slope a
measured by clsuter population studies, and are statistically independent of
the baryon fraction measurements. The model for the slope is the log-
derivative w.r.t. to log-mass of model predictions for the fractions. In grey
the slope measurement we discard due to its implausibly tight error bar. The
left panel shows the bacco fit, while the right shows the S19-model fit. Due
to anomalies in the data, the BCMs are only able to capture the qualitative
trends in the data.

We sample a Gaussian likelihood for the measurements reported in Table 1 and
Table 2. The mean values of the individual data point are given by the model
prediction of the BCMs, and the respective variances are given by the square
of the measurement uncertainties reported above (cf. section 2, Table 1, and
Table 2). The measurements are treated as independent, as justified in section
2. The model prediction for the gas and stellar fraction from the S19-model is
computed by integrating the respective profiles to $R_{\text{500c}}$, while
the prediction of the bacco model relies on a combination of functions
provided by the package bacco999https://bacco.dipc.org/ (Aricò et al., 2021b)
and the necessary conversions described in section 3.2. Predictions for the
electron density profiles directly result from appropriate re-scaling of the
gas density profile (cf. Section. 3.2.4)

We employ the Monte Carlo Markov Chain sampler
emcee101010https://emcee.readthedocs.io/en/stable/ (Goodman & Weare, 2010) to
sample the posteriors. We sample all 7 parameters of the BCMs within flat
priors (for the exact ranges, see Table 3). The cosmic baryon density
$\Omega_{\text{b}}$ is sampled in all cases within the prior ranges
$(0.0473,0.0535)$ to emulate a residual uncertainty in the cosmic baryon
fraction, which impacts the gas mass fraction predictions and the amplitude of
the predicted electron density profiles.

The cosmological parameters are kept fixed at $\Omega_{\text{M}}=0.315$, and
$h=0.67$. The choice of cosmic baryon density $\Omega_{\text{b}}$ reflects the
range of cosmic baryon fractions sampled by Schneider et al. (2022),
$0.15<f_{\text{b}}<0.17$. The Hubble constant is kept fixed as our analysis
is, to first order, independent of the Hubble parameter. The BCMs require
input masses in units of $M_{\odot}$/$h$ by design. Observations only have
access to the angular scale of clusters or groups. As such, even if not
explicitly reported, observational masses are effectively in units of
$M_{\odot}$/$h$. Similarly, the electron density profile is effectively
measured $h$-scaled units, matching the density profile predictions. To first
order, this analysis thus depends on the main cosmological parameters only via
the cosmic baryon fraction. In future work, especially when used as priors for
cosmological analyses, we plan to take account also of the residual
cosmological dependence.

##  4 Results

In this section, we first discuss the posteriors of our analysis, both with
regard to the model parameters themselves, as well as the posterior predictive
distributions of the data. We then derive the posterior predictive
distributions for the matter power spectrum suppression due to baryonic
effects.

###  4.1 Parameter Posteriors and Goodness of Fit

Figure 3: Fit to the electron density profile (black points), re-scaled by the
radius to reduce the dynamical range of the plot. Color coded the different
BCM fits. Computing the chi-squared of the median posterior predictions and
the data, we find that both models provide a good fit to the data.

When constructing a posterior sample, we discard the burn-in phase of our
sampling and visually inspect the trace plots to ensure convergence,
especially the trace plot of the total log-likelihood. For sake of brevity,
the 1- and 2-d marginal posteriors are reported in Appendix, Fig. 6 and Fig.
7.

For the bacco model, we find $\log_{10}M_{\text{c}}h/M_{\odot}=13.82\pm 0.36$,
$\log_{10}\beta=-0.29\pm 0.11$, and $\log_{10}M_{1}h/M_{\odot}=12.00\pm 0.06$.
For the parameters defining the scales of the gas profile, we only find upper
limits: $\log_{10}\theta_{\text{inn}}<-0.65$ and
$\log_{10}\theta_{\text{out}}<-0.11$ (at the 97.5th percentile). The position
of the expelled gas $\eta$, and the characteristic mass of the inner gas
profile slope $M_{\text{inn}}$ remain unconstrained. As seen in the 1- and 2-d
marginal contours shown in Fig. 6, we find a major degeneracy between the
characteristic mass of the inner gas profile slope $M_{\text{inn}}$ and
corresponding scale $\log_{10}\theta_{\text{inn}}$.

For the S19 model, we find $\log_{10}M_{\text{c}}h/M_{\odot}=14.53\pm 0.20$,
and $\mu=0.54\pm 0.10$. We get upper limits on $\theta_{\text{ej}}<6.72$, and
$\gamma<1.42$ (at the 2.5th percentile); and a lower limit on $\delta>5.37$
(at the 97.5th percentile). The mass slopes of the satellite stellar mass
$\eta$, and the difference between the satellite and central stellar mass
slopes $\delta\eta$ are strongly degenerate with each other (see Fig. 7). The
best constrained parameter combination is
$\eta\left(\delta\eta/0.24\right)^{0.19}=0.239\pm 0.009$.

In Fig. 2 we present the model prediction for the ICM (magenta) and stellar
(orange) mass fraction (upper panels) as a function of mass over plotting the
data points used for the fitting. Similarly, in the lowers panel of Fig. 2 we
show the model prediction for the logarithmic mass slope of the ICM and
stellar mass fraction (same color scheme as upper panel), together with the
slope measurements. The left images show the best fit bacco model for the
baryon fractions and their mass slopes. The right plot shows the baryon
fraction and mass slope predictions for the S19-model. Both models that we
investigate are able to capture the main trends, while the overall quality of
the fit is hampered by anomalies in the data.

A closer inspection reveals a few anomalies. In grey we plot the discarded
slope measurement by Mantz et al. (2016), which has implausibly small error
bars, as can be seen in relation to the other analyses. Furthermore, the
eFEDS-HSC ICM fraction measurement seems to be several $\sigma$ low with
respect to the predicted range of gas fractions at that mass, independently on
the BCM. We discuss this below (cf. Section 5). The model predicts very
accurately the logarithmic slope of the ICM mass fraction. The logarithmic
mass slope of the stellar mass fraction from the bacco model is almost 2 sigma
lower than our measurements. Conversely, the AS model seems to better predict
the stellar mass fraction slopes, but at the cost of under predicting the
stellar fraction of the HSC-XXL analysis. We discuss these minor anomalies
below (cf. Section 5).

The fit of both models fares better on the electron density data, shown in
Fig. 3. Here the median posterior predictive derived from the bacco model
attains a chi-squared of $\chi^{2}=3.8$ on the 8 data points of the electron
density profile, while effectively constraining 3 model parameters and having
a tight degeneracy on another pair of parameters. For the S19 model, the
median posterior predictive has a chi-squared of $\chi^{2}=4.8$ on the 8 data
points of the electron density profile, while effectively constraining 2 model
parameters, placing two limits on other parameters, and having a tight
degeneracy on another pair of parameters. The chi-squared values of the median
model predictions are in line with the expectation from the effective number
of degrees of freedom. This means that the BCM are able to fit the electron
density data well.

###  4.2 Impact on the matter power spectrum

Given our posterior of the BCM model parameters, we can predict the
suppression of the matter power spectrum induced by the evacuation of baryons
from massive halos. Note that our upper limit on the outer characteristic
scale of the bacco model ($\theta_{\text{out}}<-0.11$) falls outside of the
range of the bacco emulator ($0<\theta_{\text{out}}<0.5$, Aricò et al.,
2021b). We therefore recompute the matter power spectrum suppression for this
model for $\mathcal{O}(1000)$ points drawn from our posterior. For the
S19-model, we use the BCemu (Giri & Schneider, 2021) to compute the baryon
suppression of the matter power spectrum. The resulting posterior predictive
distributions are shown in Fig. 4, with the full line representing the 2.5th,
50th, and 97.5th percentile at each wave-number $k$, while the filled regions
encompass the range between the 16th and 84th percentile. In light green we
show the suppression for the bacco-model, while the orange shows the
suppression for the S19-model. Via the BCM models, our data is able to predict
matter power spectrum suppressions on all scales of interest with a noteworthy
precision.

We find excellent agreement between the matter power suppression predictions
from the two models on scales larger than $k<6h$ Mpc-1. Both models quickly
converge towards no suppression at scales larger than $k<0.5h$ Mpc-1.
Specifically, for $k<0.37h$ Mpc-1 ($k<0.29h$ Mpc-1) 97.5 percent of our
predictive posterior shows as suppression weaker than 1% for the bacco (S19)
model, indicating that baryon feedback processes do not impact the matter
power spectrum at larger scales.

The two BCMs also show the same strong suppression trend at intermediate
scales, $0.5h$ Mpc-1 $<k<6h$ Mpc-1. The predicted suppression is still rather
weak at $k=1h$ Mpc-1: $0.042^{+0.012}_{-0.014}$ for bacco, and
$0.049^{+0.016}_{-0.012}$ for the S19-model. The suppression grows sizeably at
$k=3h$ Mpc-1, predicting $0.184^{+0.026}_{-0.031}$ for bacco, and
$0.179^{+0.018}_{-0.020}$ for the S19-model. At smaller scales, the bacco
prediction has the characteristic up-turn at larger scales than the S19
prediction. The prediction of the two models thus diverges. At $k=10h$ Mpc-1
we find $0.185^{+0.023}_{-0.017}$ for bacco, and $0.271^{+0.023}_{-0.033}$ for
the S19-model. This discrepancy at small scales is to be expected given our
lack of data for $M_{\text{500c}}<10^{14}M_{\odot}$. Aricò et al. (2021a)
highlight that halos in the mass range
$10^{14}h^{-1}M_{\odot}<M_{\text{200c}}<10^{15}h^{-1}M_{\odot}$, contribute
most to the matter power spectrum suppression above scales of $k\lesssim 3h$
Mpc-1, while halos in the mass range of
$10^{13}h^{-1}M_{\odot}<M_{\text{200c}}<10^{14}h^{-1}M_{\odot}$ contribute the
most to scales $k\gtrsim 3h$ Mpc-1. We highlight this by shading the area in
Fig. 4 that we do not directly constrain with our data. Discrepancies in the
stellar mass fraction results might also play a role in the divergence of the
two model predictions, as discussed below.

Figure 4: Posterior predictive distribution of the suppression in the matter
power spectrum, that is the ratio of the matter power spectrum with baryonic
effects, and the gravity only matter power spectrum, $P_{\rm
with\,baryons}/P_{\rm gravity\,only}$, obtained with the S19 and bacco models
(orange and green shaded area, respectively) informed with the galaxy cluster
data collected in this work. The full lines represent the 2.5th, 50th, and
97.5th percentile at each wave-number $k$, while the filled regions encompass
the range between the 16th and 84th percentile. Greyed out the scales impacted
by extrapolating our data to lower halo masses. Upper left: Comparison with
S19-model prediction from KiDS 1000 cosmic shear, kSZ measurements and
individual cluster with hydrostatic masses (Schneider et al., 2022, hatched
blue,), which find a two sigma stronger suppression. Upper right: Comparison
with the bacco prediction from DES Y3 cosmic shear (Aricò et al., 2023,
hatched red), which is in good statistical agreement. Lower left: in cyan the
suppression needed to reconcile DES Y3 cosmic shear and Planck primary CMB
constraints (a.k.a. S8 tension, hatched in cyan, Preston et al., 2023).
Predictions based on galaxy clusters can not explain the large scale
suppression needed to solve the S8 tension. Lower right: the prediction of
several state-of-the-art hydrodynamical simulations, according to the legend.
Our data suggest a intermediate to strong feedback scenario.

##  5 Discussion

In the following, we will discuss several aspects of this work.

##### Weak lensing masses

This work provides a clear improvement compared to prior attempts to use the
gas mass – halo mass relation of galaxy clusters and groups to constrain the
baryonic suppression of the matter power spectrum (Schneider et al., 2019;
Debackere et al., 2020; Giri & Schneider, 2021), which use hydrostatic mass
estimates to interpret bayron fractions. Those works showed that the limiting
factor in using hydrostatic masses was the weakly constrained hydrostatic mass
bias. Also recent empirical calibrations of hydro-dynamical simulations are
limited by the use of hydrostatic masses (Kugel et al., 2023). This
uncertainty is reduced by using weak lensing mass calibration for the baryon
fractions. Furthermore, the cluster and group measurements in those works did
not account for selection effects via Bayesian population models, while the
analyses used in this work do. In future work, we plan to also use weak
lensing mass calibration and model selection effects also for the exploitation
of electron density profiles.

##### Comparison to previous results

Previous work by Giri & Schneider (2021); Schneider et al. (2022) fitted the
S19 model with various data sets. While neither of the two works resulted in a
constraint on any of the BCM model parameters, their posteriors suggest that
$\log_{10}M_{\text{c}}<14$, which is in slight tension with our measurement.
Both works also indicate that $\mu\sim 0.4-0.5$, though at low significance.
Our measurement $\mu=0.54\pm 0.10$ is in agreement with that. Also our
constraint that $\eta=0.239\pm 0.009$ at $\delta\eta=0.24$ coincides with the
posterior distributions they plot. Aricò et al. (2023) recently analysed the
Dark Energy survey year 3 cosmic shear auto-correlation with the bacco model.
They find a value of $\log_{10}M_{\text{c}}h/M_{\odot}=14.38^{+0.60}_{-0.56}$
– in agreement with our measurement. We therefore confirm the best fit small
scale baryonic effects found by Aricò et al. (2023) when analyzing cluster gas
density profiles, and gas and stellar fractions. Our results are also in
excellent agreement with the analysis by Chen et al. (2022b) on cosmic shear
data only, though in that work, 6 of the 7 parameters where fixed to the best
fit values from the BAHAMAS simulation. Compared to those works, which only
constrain one BCM parameter, we are able to measure 3 out of 7 parameters
tightly, while placing upper limits on the other 2. In Fig. 4, upper panels,
we directly compare our matter power spectrum suppression predictions with the
ones from the bacco fit to DES Y3 cosmic shear Aricò et al. (2023) (bacco, DES
Y3 WL in hatched red), and from the fit of the S19-model to KiDS cosmic shear,
kSZ measurements and a compilation of gas and stellar fractions (Schneider et
al., 2022). While the former is in excellent agreement with our results, the
prediction by Schneider et al. (2022) is 2-sigma lower on all scales. If this
is due to our inclusion of the X-ray profiles, or their use of kSZ data, is
left to future investigation.

##### S8 tension

Current cosmological constraints from wide photometric surveys suggest that
the amplitude of fluctuation in the low redshift Universe is lower than
predicted by extrapolating cosmic Microwave background experiments (see
Huterer, 2022, for a review). The recent CMB lensing analysis of ACT (Qu et
al., 2023), which is in agreement with primary CMB predictions, indicates that
this tension might not be present at large scales and redshift larger than
$\sim 1$. As pointed out by Amon & Efstathiou (2022); Preston et al. (2023),
the solution would likely have to be non-linear, and at low redshift. To
elaborate on this, we compare the matter power spectrum suppression derived by
redistributing baryons in and around halos, with the suppression needed to
reconcile DES Y3 cosmic shear and Planck primary CMB constraints, derived by
Preston et al. (2023), whose 1-sigma region is shown in cyan in Fig. 4, lower
left panel. That work, as well as a comparable analysis with KiDS cosmic shear
(Amon & Efstathiou, 2022), find that a solution to the "S8 tension" would
require strong suppression of large scale, $0.2h$ Mpc${}^{-1}<k<1h$ Mpc-1.
Given our predictions, it seems unlikely the baryon feedback can provide such
a suppression by redistribution baryons in and around halos. Using
$k\approx\pi/R$, solving the S8 tension would require a suppression at scales
of from $3h^{-1}$ Mpc to $15h^{-1}$ Mpc, far outside even the outermost
regions of halos. Following the interpretation by Amon & Efstathiou (2022),
this would hint at yet unknown physical effects impacting the dark matter at
(mildly) non linear scales. This also confirms that in the context of cosmic
shear, appropriate scale cuts can make that experiments impervious to baryonic
effects by only considering large scales (Krause et al., 2021; Aricò et al.,
2023; Dark Energy Survey & Kilo-Degree Survey Collaboration et al., 2023).
Previous constraints on BCM model confirm that baryon feedback only plays a
sub-ordinate role in solving the S8 tension. Indeed, a re-analysis of the Kilo
Degree Survey 1000 cosmic shear data by Schneider et al. (2022) using the S19
model, shows that the baryon feedback can not resolve, and only marginally
alleviate the discrepancy between cosmic shear and CMB.

##### Hydro-dynamical simulations

We compare in Fig. 4, lower right panel, our predicted matter power
suppression to several point estimates from state-of-the-art hydro-dynamical
simulations: Magneticum (Martinet et al., 2021, and references therein),
TNG300 (Pillepich et al., 2018; Springel et al., 2018), BAHAMAS (McCarthy et
al., 2017; McCarthy et al., 2018), EAGLE (Schaye et al., 2015a; Crain et al.,
2015; McAlpine et al., 2016), Illustris (Vogelsberger et al., 2014), and OWLS
(Schaye et al., 2010). Both very strong feedback models (like the Illustris
simulation), and very weak feedback models (like TNG300 and EAGLE) are
effectively outside of the combined range of our constraints, even when
considering the variation introduced by the difference between the
baryonification models. Our posteriors qualitatively follow most closely the
trend of the OWLS simulations. Our work further supports the strong feedback
scenario previously empirically suggested by Eckert et al. (2016); Giri &
Schneider (2021).

##### Selection effects

The poor fit the BCM provides to the gas fraction of eFEDS-HSC (see Fig. 2,
upper panel) could be indicative of some of the challenges faced by Bayesian
Population studies of galaxy clusters and groups. Objects of the same mass
with higher gas mass fraction would plausibly also be more X-ray luminous as
objects who have retained, at the same mass, less baryons (for simulational
evidence, see Ragagnin et al., 2022). Popesso et al. (2023) have shown, that
X-ray selected groups are more concentrated, both in X-rays as well as in
total mass. They also preferentially live in nodes of cosmic web and have
redder central galaxies than groups detected only via friends-of-friends
algorithms in spectroscopic galaxy data. In a Bayesian Population model, this
can be accounted for by a correlation between the intrinsic scatters of
different observables (Angulo et al., 2012; Mantz et al., 2015; Bocquet et
al., 2015; Farahi et al., 2019; Grandis et al., 2021a). While Akino et al.
(2022) did not explicitly fit for a possible correlation between selection
observable scatter and gas mass scatter, Chiu et al. (2022) find that such
correlation is consistent with zero, albeit with large error bars. Both these
analyses will be soon superseded by the publication of the first eROSITA all
sky survey (Predehl et al., 2021). In that context, multi wavelength cross
check to validate the selection modelling, such as proposed by Grandis et al.
(2020, 2021a), should be performed. We plan to revisit the predictions of the
matter power spectrum suppression once low mass halo selection in better
understood.

##### Internal consistency

Further anomalies can be found in our compilation of data sets. For instance,
if one was to take the two measurements of the stellar mass fraction as nodes
for a linear interpolation, one would find, after accounting for errors, a
stellar fraction slope of $B_{\star}-1=-0.58\pm 0.10$. This is more than
2$\sigma$ (and less than 3$\sigma$) steeper than the measured slopes
$B_{\star}-1=-0.20\pm 0.11$, and $B_{\star}-1=-0.20\pm 0.12$, from HSC-XXL,
and SPT respectively. This explains why neither of our models is able to
simultaneously fit the stellar mass fractions and their mass slopes. Such
inconsistencies could be avoided by fitting larger future cluster and groups
samples directly with the prescriptions from baryon correction models. The
relative stiffness of the BCM models also results in rather low stellar
fractions, especially in the S19 model. This is likely the reason for the
different behaviours of the BCM predictions at small scales. Less stellar
material would imply more baryons as hot gas. In light of the measured gas
fraction, it would seem that we do not find this excess hot gas in the halos.
To reconcile this, we would need strong feedback. In summary, in our inference
set-up, under-predicted stellar fractions lead to more feedback.

##### Mass measurement accuracy

Another limiting factor for the accuracy of gas and stellar mass fraction
measurements for galaxy clusters and groups is the accuracy of the weak
lensing mass estimation. Many of the observational systematics in shape and
photometric redshift measurement, as well as background selection will be
alleviated by the increased data quality of future data provided by surveys
like Euclid and LSST. Grandis et al. (2021b) however showed that also the
theoretical knowledge of the total matter distribution of halos can lead to
systematic uncertainties no smaller than 2% (see also Debackere et al., 2021).
This systematics floor results from the comparison between Magneticum and
TNG300. Given that the weak feedback scenario in the latter is excluded by
this work, the actual uncertainty on baryon feedback might well be smaller
than the reported 2$\%$. In this context data constrained BCM models, like the
one we present here, can be used to marginalise over the residual theoretical
uncertainty in baryon feedback effects on the matter profiles of galaxy
clusters.

##### Modelling improvements

At a qualitative level the performance of the BCMs in fitting our data
satisfies us, so we see no immediate source of concern. It is also noticeable
that the high precision of the used data can be effectively propagated by the
BCMs to high precision matter power spectrum suppression predictions.
Concerning the predictive accuracy of the BCMs, Schneider et al. (2019); Giri
& Schneider (2021); Schneider et al. (2022) explicitly show that when fitting
the baryon fraction measured in different hydro-dynamical simulations, the
S19-model is able to correctly predict the matter power spectrum suppression,
and vice versa. While the S19 model is also able to fit gas density profiles
derived from X-ray observations, it has not been shown explicitly that fits to
the electron density profiles and baryonic fraction correctly predict the
matter power spectrum. For the bacco-model, Aricò et al. (2021a) have shown
that is able to jointly fit the matter power spectrum and bi-spectrum in
hydrodynamical simulations, and thereby reproduce baryon fractions. Separate
fits are able to describe the gas and stellar density profiles. The predictive
accuracy on the matter power spectrum of fits to baryon fraction and electron
density profiles remains to be tested. Similarly, future work should consider
the redshift evolution of baryon feedback, both in the collection of data, as
well as the modelling.

##### Other data sets

While cluster baryon fractions are clearly a potent way to empirically
constrain baryon feedback, the only put constraints of a rather narrow range
of scales, given the limited mass range they probe. For scales $k\gtrapprox
1h$ Mpc-1, the X-ray surface brightness profile of clusters provide further
empirical constraint, as shown in this work. We only used a small fraction of
the recorded data. It will surely be fruitful to ingest more X-ray data,
especially at lower halo masses, into the BCM fitting procedure, an effort we
hope this work provides the motivation for. At even smaller scales, the
relation between central galaxy stellar mass and host halo mass for low mass
groups will be crucial. Also here, ample observational information exists,
though, in the authors opinion, rarely in an easy to access format. For both
of these observational channels, the challenge is to interface existing and
ongoing observations with the BCM. It remains to be tested if data covering a
wider mass range would be well fitted by the BCM, or if the current
extrapolation to lower masses leads to artificially precise, but inaccurate
prediction. On scales larger than $k<2h$ Mpc-1, the hot gas can be traced with
the Sunyeav-Zeldovich effect (SZe), as recently shown by Pandey et al. (2019);
Chen et al. (2022a); Gatti et al. (2022); Tröster et al. (2022) for thermal
SZe. Schneider et al. (2022); Chaves-Montero et al. (2021) showed constraints
from kinematic SZe, which requires less modelling assumptions, but is also
lower signal to noise than thermal SZe. The dispersion of localized fast radio
bursts also probes the cosmic baryon fraction (Macquart et al., 2020), and
future larger samples of localised fast radio bursts will be able to probe the
baryon distribution (Nicola et al., 2022). Another approach is to directly
exploit the high signal to noise of small scale galaxy and cosmic shear
lensing to constrain the matter power spectrum suppression (Yoon & Jee, 2021;
Huang et al., 2021; Amon & Efstathiou, 2022; Chen et al., 2022b; Preston et
al., 2023; Aricò et al., 2023), with the disadvantage that it would not be
clear, if such a suppression signal is physically due to baryon feedback, or
some other physically interesting effect on non-linear scales. The tightest
and most robust constraints will naturally result from the combination of all
these observations.

##  6 Conclusions

In this work, we use gas and stellar mass fraction measurements derived by
several weak lensing informed multi wavelength galaxy cluster studies (Mantz
et al., 2016; Chiu et al., 2018, 2022; Akino et al., 2022), as well as the
electron density profile derived by deep, high resolution X-ray observations
(Ghirardini et al., 2019), to empirically inform the Baryon Correction Model
proposed by Aricò et al. (2020); Aricò et al. (2021a) and by Schneider &
Teyssier (2015); Schneider et al. (2019). Despite some anomalies in the
observed data, we find that the data points used are qualitatively described
by the proposed models, once the free parameters of these models are
constrained by Bayesian inference. Using our constraints on the BCM
parameters, we predict the suppression of the non linear matter power spectrum
with sub-percent precision. We find that at scales smaller than $k<0.37h$
Mpc-1 the baryon suppression is $<1\%$ at 97.5 percent credibility. At $k=1h$
Mpc-1 we find a suppression of $0.042^{+0.012}_{-0.014}$ for bacco, and
$0.049^{+0.016}_{-0.012}$ for the S19-model, while at $k=3h$ Mpc-1, we find
$0.184^{+0.026}_{-0.031}$ for bacco, and $0.179^{+0.018}_{-0.020}$ for the
S19-model. The predictions diverge at small scales for $k>5h$ Mpc-1 most
likely due to data anomalies driving different model fits, and our limited
range in halo mass leading to extrapolation inaccuracies.

Given our predictions, we can exclude a series on simulations that feature
both too strong, or too weak feedback effects. Our predicted matter power
suppression seems to be reproduced most closely by the hydro-dynamical
simulation results from OWLS (Schaye et al., 2010), while other simulations
deviated from our posterior predictive on small or large scales. Constraints
on the baryonic effects derived from small scale cosmic shear by Chen et al.
(2022b); Aricò et al. (2023) are empirically confirmed by our fit to stellar
and gas mass fractions, as well as electron density profiles. Comparing our
results with the matter power spectrum suppression needed to solve the S8
tension (Amon & Efstathiou, 2022; Preston et al., 2023), we find that baryon
feedback can likely not provide the necessary large scale suppression, and
therefore is unlikely to solve the S8 tension.

In summary, we demonstrate that current measurements of galaxy cluster gas and
stellar fraction, and cluster gas profiles can strongly constrain the baryon
correction to the matter power spectrum. The cluster data can thus act as
ancillary data products, calibrating astrophysical uncertainties in cosmic
shear experiments. This opens up the possibility of using small scale cosmic
shear measurement to constrain deviation from the standard cosmological model,
while controlling the baryonic effects via external data, if other systematics
are similar well constrained. The improved understanding of baryon feedback
will also reduce theoretical uncertainties in weak lensing calibrated cluster
number counts.

## Acknowledgements

The authors thank Sebastian Bocquet, Inon Chiu, Tim Schrabback, Daisuke Nagai,
Vittorio Ghirardini, Dominique Eckert, Alexandra Amon and Raul Angulo, as well
as the eROSITA Cluster Working Group for the useful comments provided at
different stages of this work. Marginal contour plots of high-dimensional
samples are visualized with pyGTC (Bocquet & Carter, 2016).

## Data Availability

All data will be made available upon reasonable request to the authors.

## References

  * Akino et al. (2022) Akino D., et al., 2022, PASJ, 74, 175
  * Amon & Efstathiou (2022) Amon A., Efstathiou G., 2022, MNRAS, 516, 5355
  * Angulo et al. (2012) Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
  * Angulo et al. (2021) Angulo R. E., Zennaro M., Contreras S., Aricò G., Pellejero-Ibañez M., Stücker J., 2021, MNRAS, 507, 5869
  * Aricò et al. (2020) Aricò G., Angulo R. E., Hernández-Monteagudo C., Contreras S., Zennaro M., Pellejero-Ibañez M., Rosas-Guevara Y., 2020, MNRAS, 495, 4800
  * Aricò et al. (2021a) Aricò G., Angulo R. E., Hernández-Monteagudo C., Contreras S., Zennaro M., 2021a, MNRAS, 503, 3596
  * Aricò et al. (2021b) Aricò G., Angulo R. E., Contreras S., Ondaro-Mallea L., Pellejero-Ibañez M., Zennaro M., 2021b, MNRAS, 506, 4070
  * Aricò et al. (2023) Aricò G., Angulo R. E., Zennaro M., Contreras S., Chen A., Hernández-Monteagudo C., 2023, arXiv e-prints, p. arXiv:2303.05537
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57
  * Bocquet & Carter (2016) Bocquet S., Carter F. W., 2016, The Journal of Open Source Software, 1
  * Bocquet et al. (2015) Bocquet S., et al., 2015, ApJ, 799, 214
  * Bocquet et al. (2019) Bocquet S., et al., 2019, ApJ, 878, 55
  * Bulbul et al. (2019) Bulbul E., et al., 2019, ApJ, 871, 50
  * Capasso et al. (2019) Capasso R., et al., 2019, MNRAS, 482, 1043
  * Carlstrom et al. (2002) Carlstrom J. E., Holder G. P., Reese E. D., 2002, ARA&A, 40, 643
  * Chaves-Montero et al. (2021) Chaves-Montero J., Hernández-Monteagudo C., Angulo R. E., Emberson J. D., 2021, MNRAS, 503, 1798
  * Chen et al. (2022a) Chen Z., Zhang P., Yang X., 2022a, arXiv e-prints, p. arXiv:2201.12591
  * Chen et al. (2022b) Chen A., et al., 2022b, arXiv e-prints, p. arXiv:2206.08591
  * Chisari et al. (2018) Chisari N. E., et al., 2018, MNRAS, 480, 3962
  * Chisari et al. (2019) Chisari N. E., et al., 2019, The Open Journal of Astrophysics, 2, 4
  * Chiu et al. (2018) Chiu I., et al., 2018, MNRAS, 478, 3072
  * Chiu et al. (2022) Chiu I. N., et al., 2022, A&A, 661, A11
  * Chiu et al. (2023) Chiu I. N., Klein M., Mohr J., Bocquet S., 2023, MNRAS, 522, 1601
  * Contreras et al. (2020) Contreras S., Angulo R. E., Zennaro M., Aricò G., Pellejero-Ibañez M., 2020, MNRAS, 499, 4905
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * Croston et al. (2008) Croston J. H., et al., 2008, A&A, 487, 431
  * Dark Energy Survey & Kilo-Degree Survey Collaboration et al. (2023) Dark Energy Survey & Kilo-Degree Survey Collaboration et al., 2023, arXiv e-prints, p. arXiv:2305.17173
  * Debackere et al. (2020) Debackere S. N. B., Schaye J., Hoekstra H., 2020, MNRAS, 492, 2285
  * Debackere et al. (2021) Debackere S. N. B., Schaye J., Hoekstra H., 2021, MNRAS, 505, 593
  * Dietrich et al. (2019) Dietrich J. P., et al., 2019, MNRAS, 483, 2871
  * Eckert et al. (2016) Eckert D., et al., 2016, A&A, 592, A12
  * Eckert et al. (2017) Eckert D., Ettori S., Pointecouteau E., Molendi S., Paltani S., Tchernin C., 2017, Astronomische Nachrichten, 338, 293
  * Eckert et al. (2021) Eckert D., Gaspari M., Gastaldello F., Le Brun A. M. C., O’Sullivan E., 2021, Universe, 7, 142
  * Ettori et al. (2010) Ettori S., Gastaldello F., Leccardi A., Molendi S., Rossetti M., Buote D., Meneghetti M., 2010, A&A, 524, A68
  * Farahi et al. (2019) Farahi A., et al., 2019, Nature Communications, 10, 2504
  * Gatti et al. (2022) Gatti M., et al., 2022, Phys. Rev. D, 105, 123525
  * Ghirardini et al. (2019) Ghirardini V., et al., 2019, A&A, 621, A41
  * Giri & Schneider (2021) Giri S. K., Schneider A., 2021, J. Cosmology Astropart. Phys., 2021, 046
  * Goodman & Weare (2010) Goodman J., Weare J., 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Grandis et al. (2020) Grandis S., et al., 2020, MNRAS, 498, 771
  * Grandis et al. (2021a) Grandis S., et al., 2021a, MNRAS, 504, 1253
  * Grandis et al. (2021b) Grandis S., Bocquet S., Mohr J. J., Klein M., Dolag K., 2021b, MNRAS, 507, 5671
  * Harnois-Déraps et al. (2015) Harnois-Déraps J., van Waerbeke L., Viola M., Heymans C., 2015, MNRAS, 450, 1212
  * Huang et al. (2021) Huang H.-J., et al., 2021, MNRAS, 502, 6010
  * Hurier & Angulo (2018) Hurier G., Angulo R. E., 2018, A&A, 610, L4
  * Huterer (2022) Huterer D., 2022, arXiv e-prints, p. arXiv:2212.05003
  * Kaiser (1986) Kaiser N., 1986, MNRAS, 222, 323
  * Krause et al. (2021) Krause E., et al., 2021, arXiv e-prints, p. arXiv:2105.13548
  * Kravtsov et al. (2018) Kravtsov A. V., Vikhlinin A. A., Meshcheryakov A. V., 2018, Astronomy Letters, 44, 8
  * Kugel et al. (2023) Kugel R., et al., 2023, arXiv e-prints, p. arXiv:2306.05492
  * Macquart et al. (2020) Macquart J. P., et al., 2020, Nature, 581, 391
  * Mantz et al. (2015) Mantz A. B., et al., 2015, MNRAS, 446, 2205
  * Mantz et al. (2016) Mantz A. B., et al., 2016, MNRAS, 463, 3582
  * Martinet et al. (2021) Martinet N., Castro T., Harnois-Déraps J., Jullo E., Giocoli C., Dolag K., 2021, A&A, 648, A115
  * McAlpine et al. (2016) McAlpine S., et al., 2016, Astronomy and Computing, 15, 72
  * McCarthy et al. (2017) McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465, 2936
  * McCarthy et al. (2018) McCarthy I. G., Bird S., Schaye J., Harnois-Deraps J., Font A. S., van Waerbeke L., 2018, MNRAS, 476, 2999
  * McDonald et al. (2013) McDonald M., et al., 2013, ApJ, 774, 23
  * Mead et al. (2015) Mead A. J., Peacock J. A., Heymans C., Joudaki S., Heavens A. F., 2015, MNRAS, 454, 1958
  * Mead et al. (2021) Mead A. J., Brieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401
  * Mulroy et al. (2019) Mulroy S. L., et al., 2019, MNRAS, 484, 60
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Nicola et al. (2022) Nicola A., et al., 2022, J. Cosmology Astropart. Phys., 2022, 046
  * Pakmor et al. (2022) Pakmor R., et al., 2022, arXiv e-prints, p. arXiv:2210.10060
  * Pandey et al. (2019) Pandey S., et al., 2019, Phys. Rev. D, 100, 063519
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Popesso et al. (2023) Popesso P., et al., 2023, arXiv e-prints, p. arXiv:2302.08405
  * Predehl et al. (2021) Predehl P., et al., 2021, A&A, 647, A1
  * Preston et al. (2023) Preston C., Amon A., Efstathiou G., 2023, arXiv e-prints, p. arXiv:2305.09827
  * Qu et al. (2023) Qu F. J., et al., 2023, arXiv e-prints, p. arXiv:2304.05202
  * Ragagnin et al. (2021) Ragagnin A., Saro A., Singh P., Dolag K., 2021, MNRAS, 500, 5056
  * Ragagnin et al. (2022) Ragagnin A., Andreon S., Puddu E., 2022, A&A, 666, A22
  * Salcido et al. (2023) Salcido J., McCarthy I. G., Kwan J., Upadhye A., Font A. S., 2023, MNRAS, 523, 2247
  * Sanders et al. (2018) Sanders J. S., Fabian A. C., Russell H. R., Walker S. A., 2018, MNRAS, 474, 1065
  * Schaye et al. (2010) Schaye J., et al., 2010, MNRAS, 402, 1536
  * Schaye et al. (2015a) Schaye J., et al., 2015a, MNRAS, 446, 521
  * Schaye et al. (2015b) Schaye J., et al., 2015b, MNRAS, 446, 521
  * Schaye et al. (2023) Schaye J., et al., 2023, arXiv e-prints, p. arXiv:2306.04024
  * Schneider & Teyssier (2015) Schneider A., Teyssier R., 2015, J. Cosmology Astropart. Phys., 2015, 049
  * Schneider et al. (2019) Schneider A., Teyssier R., Stadel J., Chisari N. E., Le Brun A. M. C., Amara A., Refregier A., 2019, J. Cosmology Astropart. Phys., 2019, 020
  * Schneider et al. (2020) Schneider A., et al., 2020, J. Cosmology Astropart. Phys., 2020, 020
  * Schneider et al. (2022) Schneider A., Giri S. K., Amodeo S., Refregier A., 2022, MNRAS, 514, 3802
  * Schrabback et al. (2021) Schrabback T., et al., 2021, MNRAS, 505, 3923
  * Sereno & Ettori (2017) Sereno M., Ettori S., 2017, MNRAS, 468, 3322
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Stern et al. (2019) Stern C., et al., 2019, MNRAS, 485, 69
  * Tröster et al. (2022) Tröster T., et al., 2022, A&A, 660, A27
  * Umetsu et al. (2020) Umetsu K., et al., 2020, ApJ, 890, 148
  * Villaescusa-Navarro et al. (2021) Villaescusa-Navarro F., et al., 2021, Astrophys. J., 915, 71
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * Yoon & Jee (2021) Yoon M., Jee M. J., 2021, ApJ, 908, 13
  * Zennaro et al. (2023) Zennaro M., Angulo R. E., Pellejero-Ibáñez M., Stücker J., Contreras S., Aricò G., 2023, MNRAS,
  * de Haan et al. (2016) de Haan T., et al., 2016, ApJ, 832, 95
  * van Daalen et al. (2020) van Daalen M. P., McCarthy I. G., Schaye J., 2020, MNRAS, 491, 2424

##  Appendix A Conversion of gas fraction between over densities

The gas mass fraction w.r.t. a given over density $\Delta$c is given by the
ratio of the gas mass $M_{\text{ICM}}(<r_{\Delta\text{c}})$ enclosed within
the radius $r_{\Delta\text{c}}$ corresponding to that over density, and the
total mass $M_{\Delta\text{c}}$ w.r.t. that over density, i.e.

$f^{\Delta\text{c}}_{\text{ICM}}=\frac{M_{\text{ICM}}(<r_{\Delta\text{c}})}{M_{\Delta\text{c}}}.$
(13)

Thus, the ratio between the gas mass fractions at two different over densities
$\Delta$c and $\Delta^{\prime}$c, can be readily expressed by the ratios of
gas masses and total masses, as

$\frac{f^{\Delta\text{c}}_{\text{ICM}}}{f^{\Delta^{\prime}\text{c}}_{\text{ICM}}}=\frac{M_{\text{ICM}}(<r_{\Delta\text{c}})}{M_{\text{ICM}}(<r_{\Delta^{\prime}\text{c}})}\frac{M_{\Delta^{\prime}\text{c}}}{M_{\Delta\text{c}}}.$
(14)

The second term is computed via the customary transformations between over
density masses at different over densities (see for instance Ettori et al.,
2010, App. A), assuming the concentration mass relation by Ragagnin et al.
(2021) at its pivot cosmology.

The ratio between the gas masses can be estimated from the bound gas profile
$\rho_{\text{bg}}(r)$ as

$\frac{M_{\text{ICM}}(<r_{\Delta\text{c}})}{M_{\text{ICM}}(<r_{\Delta^{\prime}\text{c}})}=\int_{0}^{r_{\Delta\text{c}}}r^{2}\text{d}r\,\rho_{\text{ICM}}(r)\Big{/}\int_{0}^{r_{\Delta^{\prime}\text{c}}}r^{2}\text{d}r\,\rho_{\text{ICM}}(r),$
(15)

which is a function of the BCM parameters.

We assume that the stellar fraction stays constant, as the brightest central
galaxy is included in both $r_{500\text{c}}$ and $r_{500\text{c}}$, and the
rest of the profile follows the dark matter profile.

##  Appendix B Particle Radius in bacco

Figure 5: Matter power suppression predicted in the bacco model when
displacing all particles initially within $r_{\text{200c}}$ (cyan), or
$r_{\text{out}}$ (light green), evaluated on our posterior. The former model
causes larger particle displacement, which leads to higher suppression at
large scales. Vice versa, the latter model features shorter particle
displacement and thus smaller suppression at large scales.

One of the main differences between S19 and the bacco model is that the former
displaces formally all the particles in the simulation (although the
displacement tends to zero at large distances from halos), whereas the latter
displaces only particles in haloes, inside a radius that we shall call the
particle radius, $r_{\text{p}}$. The BCM analytical halo density profiles are
opportunely truncated at $r_{\text{p}}$, so that a long-range displacement
causes particles to be ejected from the haloes and trace the analytical
ejected gas profiles (see Aricò et al., 2020, for more details). In Aricò et
al. (2020), and following works, it was fixed $r_{\text{p}}=r_{\text{200c}}$.
This choice was motivated by the fact that the parametrisation and priors used
were such that, beyond $r_{\text{200c}}$, the ejected gas was the only
component with a shape different from NFW. However, in this work, we utilize
X-ray density profile to constrain $r_{\text{out}}>r_{\text{200c}}$, meaning
that the hot gas does not follow the dark matter beyond $r_{\text{200c}}$ (we
remind the reader that in bacco the hot gas profile follows the NFW after
$r_{\text{out}}$). For this work, we decided then to set
$r_{\text{p}}=r_{\text{out}}$, to enforce the particle radius to be at the
exact position of the change in slope of the hot gas to NFW. The choice of
$r_{\text{p}}$ leads to a qualitatively different – though statistically
consistent – matter power suppression, as shown in Fig 5. Specifically, we
find that setting $r_{\text{p}}=r_{\text{out}}$ the suppression in the matter
power spectrum is reduced at large scales, and enhanced at small scales, with
respect to $r_{\text{p}}=r_{\text{200c}}$. We note that, while the free
parameters of bacco act on the shape and normalisation of the analytical
density profiles, the particle radius $r_{\text{p}}$ impacts directly the
displacement fields applied to the particles. Smaller $r_{\text{p}}$ drive
larger displacements fields, and therefore more suppression at large scales.
e.g. in Aricò et al. (2021b). Vice versa, larger $r_{\text{p}}$ cause on
average smaller displacements in the particles, and therefore the suppression
in the matter power spectrum is more similar to S19. A priori, it is not clear
which type of displacement leads to more realistic results. More detailed
comparisons against hydrodynamical simulations, as well as consistency checks
between different observables such as cosmic shear and X-ray data, can inform
the model on the choice of $r_{\text{p}}$.

##  Appendix C Marginal Plots of the Posteriors

In Fig. 6 and Fig. 7 we show the 1-d and 2- marginalized posterior samples
resulting from the fit of the bacco, and S19 model, respectively. The
interpretation of these plots is provided in section 4.1.

Figure 6: 1-d and 2-d marginalized posterior samples resulting from the fit of
the bacco model. The parameters $\log_{10}M_{\text{c}}h/M_{\odot}$,
$\log_{10}\beta$ and $\log_{10}M_{1}h/M_{\odot}$ are well constrained (for the
sake of a compact representation, we omit the mass units in the plot labels).
We place upper limits on $\log_{10}\theta_{\text{out}}$ and
$\log_{10}\theta_{\text{inn}}$, while the other parameters remain
unconstrained, with the exception of the hydrostatic mass bias, on which we
placed a prior. Figure 7: 1-d and 2-d marginalized posterior samples resulting
from the fit of the S+19-model. The parameters
$\log_{10}M_{\text{c}}h/M_{\odot}$, and $\mu$ are well constrained (for the
sake of a compact representation, we omit the mass units in the plot labels).
We place upper limits on $\theta_{\text{ej}}$ and $\gamma$, and a lower limit
on $\delta$. A specific combination of the parameters $\eta$ and $\delta\eta$
is also well constrained.