# Inferring the impact of feedback on the matter distribution using the
Sunyaev Zel’dovich effect: Insights from CAMELS simulations and ACT+DES data

Shivam Pandey,1,2 Kai Lehman,3,4 Eric J. Baxter,3 Yueying Ni,5 Daniel Angles-
Alcazar,6,7 Shy Genel,7,1 Francisco Villaescusa-Navarro,7 Ana Maria Delgado,5
Simeon Bird,8 Tiziana di Matteo,9 Rupert Croft9, …10,
1 Department of Physics, Columbia University, New York, NY, USA 10027
2 Department of Physics and Astronomy, University of Pennsylvania,
Philadelphia, PA 19104, USA
3Institute for Astronomy, University of Hawai‘i, 2680 Woodlawn Drive,
Honolulu, HI 96822, USA
4Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-
Universität, Scheinerstr. 1, 81679 München, Germany
5Center for Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA 02138, US
6Department of Physics, University of Connecticut, 196 Auditorium Road,
U-3046, Storrs, CT, 06269, USA
7Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue,
New York, NY, 10010, USA
8Department of Physics & Astronomy, University of California, Riverside, 900
University Ave., Riverside, CA 92521, US
9McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University, Pittsburgh, PA 15213
10Affil

###### Abstract

Feedback from active galactic nuclei and stellar processes changes the matter
distribution on small scales, leading to significant systematic uncertainty in
weak lensing constraints on cosmology. We investigate how the observable
properties of group-scale halos can constrain feedback’s impact on the matter
distribution using the CAMELS simulations. Extending the results of previous
work to smaller halo masses and higher $k$, we find that the baryon fraction
in halos contains significant information about the impact of feedback on the
matter power spectrum. We explore how the thermal Sunyaev Zel’dovich (tSZ)
signal from group-scale halos contains similar information. Using recent Dark
Energy Survey (DES) weak lensing and Atacama Cosmology Telescope (ACT) tSZ
cross-correlation measurements and models trained on CAMELS, we obtain 10%
constraints on feedback effects at $k\sim 5\,h/{\rm Mpc}$. We show that with
catalogs from future surveys, such as halo positions from Dark Energy
Spectroscopic Instrument and tSZ from Simons Observatory, we will be able to
constrain baryonic effects on the power spectrum to $O(<1\%)$ at $k=1\,h/{\rm
Mpc}$ and $O(3\%)$ at $k=5\,h/{\rm Mpc}$. Finally, we investigate the impact
of feedback on the matter bispectrum, finding that tSZ observables are highly
informative in this case.

###### keywords:

††pagerange: Inferring the impact of feedback on the matter distribution using
the Sunyaev Zel’dovich effect: Insights from CAMELS simulations and ACT+DES
data–15

##  1 Introduction

The statistics of the matter distribution on scales $k\gtrsim 0.1\,h{\rm
Mpc}^{-1}$ are tightly constrained by current weak lensing surveys (e.g.
Asgari:2021; DESY3cosmo). However, modeling the matter distribution on small
scales $k\gtrsim 1\,h{\rm Mpc}^{-1}$ to extract cosmological information is
complicated by the effects of baryonic feedback (Rudd:2008). Energetic output
from active galactic nuclei (AGN) and stellar processes (e.g. winds and
supernovae) directly impacts the distribution of gas on small scales, thereby
changing the total matter distribution (e.g. Chisari:2019).111Changes to the
gas distribution can also gravitationally influence the dark matter
distribution, further modifying the total matter distribution. The coupling
between these processes and the large-scale gas distribution is challenging to
model theoretically and in simulations because of the large dynamic range
involved, from the scales of individual stars to the scales of galaxy
clusters. While it is generally agreed that feedback leads to a suppression of
the matter power spectrum on scales $0.1\,h{\rm Mpc}^{-1}\lesssim k\lesssim
20\,h{\rm Mpc}^{-1}$, the amplitude of this suppression remains uncertain by
tens of percent (vanDaalen:2020; Villaescusa-Navarro:2021:ApJ:) (see also Fig.
1). This systematic uncertainty limits constraints on cosmological parameters
from current weak lensing surveys (e.g. DESY3cosmo; Asgari:2021). For future
surveys, such as the Vera Rubin Observatory LSST
(TheLSSTDarkEnergyScienceCollaboration:2018:arXiv:) and Euclid
(EuclidCollaboration:2020:A&A:), the problem will become even more severe
given expected increases in statistical precision. In order to reduce the
uncertainties associated with feedback, we would like to identify observable
quantities that carry information about the impact of feedback on the matter
power distribution and develop approaches to extract this information (e.g.
Nicola:2022:JCAP:).

Recently, vanDaalen:2020 showed that the halo baryon fraction, $f_{b}$, in
halos with $M\sim 10^{14}\,M_{\odot}$ carries significant information about
suppression of the matter power spectrum caused by baryonic feedback. Notably,
they found that the relation between $f_{b}$ and matter power suppression was
robust to changing feedback prescription. Note that $f_{b}$ as defined by
vanDaalen:2020 counts baryons in both the intracluster medium as well as those
in stars. The connection between $f_{b}$ and feedback is expected, since one
of the main drivers of feedback’s impact on the matter distribution is the
ejection of gas from halos by AGN. Therefore, when feedback is strong, halos
will be depleted of baryons and $f_{b}$ will be lower. The conversion of
baryons into stars — which will not significantly impact the matter power
spectrum on large scales — does not impact $f_{b}$, since $f_{b}$ includes
baryons in stars as well as the ICM. vanDaalen:2020 specifically consider the
measurement of $f_{b}$ in halos with $6\times 10^{13}M_{\odot}\lesssim
M_{500c}\lesssim 10^{14}\,M_{\odot}$. In much more massive halos, the energy
output of AGN is small compared to the binding energy of the halo, preventing
gas from being expelled. In smaller halos, vanDaalen:2020 find that the
correlation between power spectrum suppression and $f_{b}$ is less clear.

Although $f_{b}$ carries information about feedback, it is somewhat unclear
how one would measure $f_{b}$ in practice. Observables such as the kinematic
Sunyaev Zel’dovich (kSZ) effect can be used to constrain the gas density;
combined with some estimate of stellar mass, $f_{b}$ could then be inferred.
However, measuring the kSZ is challenging, and current measurements have low
signal-to-noise (Hand:2012; Hill:2016; Soergel:2016). Moreover, vanDaalen:2020
consider a relatively limited range of feedback prescriptions. It is unclear
whether a broader range of feedback models could lead to a greater spread in
the relationship between $f_{b}$ and baryonic effects on the power spectrum.
In any case, it is worthwhile to consider other potential observational probes
of feedback.

Another potentially powerful probe of baryonic feedback is the thermal SZ
(tSZ) effect. The tSZ effect is caused by inverse Compton scattering of CMB
photons with a population of electrons at high temperature. This scattering
process leads to a spectral distortion in the CMB that can be reconstructed
from multi-frequency CMB observations. The amplitude of this distortion is
sensitive to the line-of-sight integral of the electron pressure. Since
feedback changes the distribution and thermodynamics of the gas, it stands to
reason that it could impact the tSZ signal. Indeed, several works using both
data (e.g Pandey:2019; Pandey:2022; Gatti:2022) and simulations (e.g.
Scannapieco:2008; Bhattacharya:2008; Wadekar:2022) have shown that the tSZ
signal from low-mass (group scale) halos is sensitive to feedback. Excitingly,
the sensitivity of tSZ measurements is expected to increase dramatically in
the near future due to high-sensitivity CMB measurements from e.g. SPT-3G
(Benson:2014:SPIE:), Advanced ACTPol (Henderson:2016:JLTP:), Simons
Observatory (Ade:2019:JCAP:), and CMB Stage 4 (CMBS4).

The goal of this work is to investigate what information the tSZ signals from
low-mass halos contain about the impact of feedback on the small-scale matter
distribution. The tSZ signal, which we denote with the Compton $y$ parameter,
carries different information from $f_{b}$. For one, $y$ is sensitive only to
the gas and not to stellar mass. Moreover, $y$ carries sensitivity to both the
gas density and pressure, unlike $f_{b}$ which depends only on the gas
density.222Of course, sensitivity to gas temperature does not necessarily mean
that the tSZ is a more useful probe of feedback. The $y$ signal is also be
easier to measure than $f_{b}$, since it can be estimated simply by cross-
correlating halos with a tSZ map. The signal-to-noise of such cross-
correlation measurements is already high with current data, on the order of
10s of $\sigma$ (Vikram:2017; Pandey:2019; Pandey:2022; Sanchez:2022).

In this paper, we investigate the information content of the tSZ signal using
the CAMELS simulations. As we describe in more detail in §2, CAMELS is a suite
of many hydrodynamical simulations run across a range of different feedback
prescriptions and different cosmological parameters. The relatively small
volume of the CAMELS simulations ($(25/h)^{3}\,{\rm Mpc^{3}}$) means that we
are somewhat limited in the halo masses and scales that we can probe. We
therefore view our analysis as an exploratory work that investigates the
information content of low-mass halos for constraining feedback and how to
extract this information; more accurate results over a wider range of halo
mass and $k$ may be obtained in the future using the same methods applied to
larger volume simulations.

By training statistical models on the CAMELS simulations, we explore what
information about feedback exists in tSZ observables, and how robust this
information is to changes in subgrid prescriptions. We consider three very
different prescriptions for feedback based on the SIMBA (Dave:2019:MNRAS:),
Illustris-TNG (Pillepich:2018:MNRAS:) and Astrid (Bird:2022:MNRAS:;
Ni:2022:MNRAS:) models across a wide range of possible parameter values
(including variations in cosmology). The flexibility of the statistical models
we employ means that it is possible to uncover more complex relationships
between e.g. $f_{b}$, $y$, and the baryonic suppression of the power spectrum
than considered in vanDaalen:2020.

Finally, we apply our trained statistical models to recently published
measurements of the $y$ signal in low-mass halos. In particular, we consider
the inferred values of $Y$ from the lensing-tSZ correlation analysis of
Atacama Cosmology Telescope (ACT) and Dark Energy Survey (DES)
Madhavacheril:2020:PhRvD:; Amon:2022:PhRvD:; Secco:2022:PhRvD:b data presented
in Gatti:2022 and Pandey:2022. In addition to providing interesting
constraints on the impact of feedback, these results highlight the potential
of future similar analyses with e.g. Dark Energy Spectroscopic Experiment
(DESI; DESI) and CMB Stage 4 (CMBS4).

Two recent works — Moser:2021 and Wadekar:2022 — have used the CAMELS
simulations to explore the information content of the tSZ signal for
constraining feedback. These works focus on the ability of tSZ observations to
constrain the parameters of subgrid feedback models in hydrodynamical
simulations. Here, in contrast, we attempt to connect the observable
quantities directly to the impact of feedback on the matter power spectrum.
Additionally, unlike some of the results presented in Moser:2021 and
Wadekar:2022, we consider the full parameter space explored by the CAMELS
simulations rather than the small variations around a fiducial point that are
relevant to calculation of the Fisher matrix. Finally, we only focus on the
intra-halo gas profile of the halos in the mass range captured by the CAMELS
simulations (c.f. Moser:2021). We do not expect the inter-halo gas pressure to
be captured by the small boxes used here as it may be sensitive to higher halo
masses (Pandey:2020).

Nonlinear evolution of the matter distribution induces non-Gaussianity, and
hence there is additional information to be recovered beyond the power
spectrum. Recent measurements detect higher-order matter correlations at
cosmological scales at $O(10\sigma)$(Secco:2022:PhRvD:b; Gatti:2022:PhRvD:),
and the significance of these measurements is expected to rapidly increase
with up-coming surveys (Pyne:2021:MNRAS:). Jointly analyzing two-point and
three-point correlations of the matter field can help with self-calibration of
systematic parameters and improve cosmological constraints. As described in
Foreman:2020:MNRAS:, the matter bispectrum is expected to be impacted by
baryonic physics at $O(10\%)$ over the scales of interest. With these
considerations in mind, we also investigate whether the SZ observations carry
information about the impact of baryonic feedback on the matter bispectrum.

The plan of the paper is as follows. In §2 we discuss the CAMELS simulation
and the data products that we use in this work. In §3, we present the results
of our explorations with the CAMELS simulations, focusing on the information
content of the tSZ signal for inferring the amount of matter power spectrum
suppression. In §4, we apply our analysis to the DES and ACT measurements. We
summarize our results and conclude in §5.

##  2 CAMELS simulations and observables

###  2.1 Overview of CAMELS simulations

Simulation | Type/Code |  | Astrophysical parameters varied
---
& its meaning
IllustrisTNG |  | Magneto-hydrodynamic/
---
AREPO
|  $A_{\rm SN1}$: (Energy of Galactic winds)/SFR
---
$A_{\rm SN2}$: Speed of galactic winds
$A_{\rm AGN1}$: Energy/(BH accretion rate)
$A_{\rm AGN2}$: Jet ejection speed or burstiness
SIMBA | Hydrodynamic/GIZMO |  |  $A_{\rm SN1}$ : Mass loading of galactic winds
---
$A_{\rm SN2}$ : Speed of galactic winds
$A_{\rm AGN1}$ : Momentum flux in QSO-mode of feedback
$A_{\rm AGN2}$ : Jet speed in kinetic mode of feedback
Astrid | Hydrodynamic/pSPH |  |  $A_{\rm SN1}$: (Energy of Galactic winds)/SFR
---
$A_{\rm SN2}$: Speed of galactic winds
$A_{\rm AGN1}$: Energy/(BH accretion rate)
$A_{\rm AGN2}$: Thermal feedback efficiency
Table 1: Summary of three varieties of simulations used in this analysis. In
addition to the four astrophysical parameters mentioned, all simulations vary
two cosmological parameters, $\Omega_{\rm m}$ and $\sigma_{8}$.

We investigate the use of SZ signals for constraining the impact of feedback
on the matter distribution using approximately 3000 cosmological simulations
run by the CAMELS collaboration (Villaescusa-Navarro:2021:ApJ:). One half of
these are gravity-only N-body simulations and the other half are
hydrodynamical simulations with matching initial conditions. The simulations
are run using three different hydrodynamical sub-grid codes, Illustris-TNG
(Pillepich:2018:MNRAS:), SIMBA (Dave:2019:MNRAS:) and Astrid
(Bird:2022:MNRAS:; Ni:2022:MNRAS:). As detailed in Villaescusa-
Navarro:2021:ApJ:, for each sub-grid implementation six parameters are varied:
two cosmological parameters ($\Omega_{m}$ and $\sigma_{8}$) and four dealing
with baryonic astrophysics. Of these, two deal with supernovae feedback
($A_{\rm SN1}$ and $A_{\rm SN2}$) and two deal with AGN feedback ($A_{\rm
AGN1}$ and $A_{\rm AGN2}$). The meanings of these parameters for each subgrid
model are summarized in Table 1.

Note that the astrophysical parameters might have somewhat different physical
meanings for different subgrid prescriptions and there is usually a complex
interplay between them regarding their impact on the properties of galaxy and
gas. For example, the parameter $A_{\rm SN1}$ approximately corresponds to the
pre-factor for the overall energy output in galactic wind feedback per-unit
star-formation in both TNG and Astrid simulations. However, in the SIMBA
simulations it corresponds to the wind mass outflow rate per unit star-
formation. Similarly, the $A_{\rm AGN2}$ parameter controls the burstiness and
the temperature of the heated gas during the AGN bursts in the TNG
simulations. Whereas in the SIMBA suite, it corresponds to the speed of the
continuously-driven AGN jets. However, in the Astrid suite, it corresponds to
the thermal feedback efficiency. As we describe in § 3.2, this can result in
counter-intuitive impact on the matter power spectrum in the Astrid
simulation, relative to TNG and SIMBA.

For each of the sub-grid physics prescriptions, three varieties of simulations
are provided. These include 27 sims for which the parameters are fixed and
initial conditions are varied (cosmic variance, or CV, set), 66 sims varying
only one parameter at a time (1P set) and 1000 sims varying parameters in a
six dimensional latin hyper-cube (LH set). We use the CV sims to estimate the
variance expected in the matter power suppression due to stochasticity (see
Fig. 1). We use the 1P sims to understand how the matter suppression responds
to variation in each parameter individually. Finally we use the full LH set to
effectively marginalize over the full parameter space varying all six
parameters. We use publicly available power spectrum and bispectrum
measurements for these simulation boxes.333https://www.camel-
simulations.org/data Where unavailable, we calculate the power spectrum and
bispectrum, using the publicly available code
Pylians.444https://github.com/franciscovillaescusa/Pylians3

###  2.2 Baryonic effects on the power spectrum and bispectrum

Figure 1: Far left: Baryonic suppression of the matter power spectrum, $\Delta
P$, in the CAMELS simulationsat . The dark-blue, red and orange shaded region
correspond to $1\sigma$ error estimated with the fiducial values CV suite of
the feedback parametersTNG, SIMBA and Astrid respectively. Each curve The
light-blue region corresponds to a different prescription for the subgrid
physics as indicated in $1\sigma$ error estimated with the legendLH suite of
TNG, showing a significantly larger spread. Middle and right panels: the
impact of baryonic physics on the matter bispectrum suppression for the same
set of simulations for equilaterla and squeezed triangle configurations
respectively.

The left panel of Fig. 1 shows the measurement of the power spectrum
suppression caused by baryonic effects in the Illustris, SIMBA, and Astrid
simulations for their fiducial feedback settings. To compute the suppression,
we use the matter power spectra and bispectra of the hydrodynamical
simulations (hydro) and the dark-matter only (DMO) simulations generated at
varying initial conditions (ICs). The power spectrum and bispectrum for all
the simulations are provided by the CAMELS collaboration and are publicly
available. For each of the 27 unique IC runs, we calculate the ratios $\Delta
P/P_{\rm DMO}=(P_{\rm hydro}-P_{\rm DMO})/P_{\rm DMO}$ and $\Delta B/B_{\rm
DMO}=(B_{\rm hydro}-B_{\rm DMO})/B_{\rm DMO}$. As the hydro-dynamical and the
N-body simulations are run with same initial conditions, the ratios $\Delta
P/P_{\rm DMO}$ and $\Delta B/B_{\rm DMO}$ are independent of sample variance.

It is clear that the amplitude of suppression of the small-scale matter power
spectrum can be significant: suppression on the order of tens of percent is
reached for all three simulations. It is also clear that the impact is
significantly different between the three simulations. Even for the
simulations in closest agreement (Illustris-TNG and Astrid), the measurements
of $\Delta P/P_{\rm DMO}$ disagree by a factor of five at $k=5\,h/{\rm Mpc}$.
The width of the curves in Fig. 1 represents the standard deviation measured
across the cosmic variance simulations, which all have the same parameter
values but different initial conditions. For bispectrum, we show both the
equilateral and squeezed triangle configurations with cosine of angle between
long sides fixed to $\mu=0.9$. Interestingly, the spread in $\Delta P/P_{\rm
DMO}$ and $\Delta B/B_{\rm DMO}$ increases with increasing $k$ over the range
$0.1\,h/{\rm Mpc}\lesssim k\lesssim 10\,h/{\rm Mpc}$. This increase is driven
by stochasticity arising from baryonic feedback. The middle and right panels
show the impact of feedback on the bispectrum for the equilateral and squeezed
triangle configurations, respectively.

Throughout this work, we will focus on the regime $0.3\,h/{\rm
Mpc}<k<10\,h/{\rm Mpc}$. Larger scales are not well constrained by modes are
not present in the $(25{\rm Mpc}/h)^{3}$ CAMELS simulations, and in any case,
the impact of feedback on large scales is typically small. Much smaller
scales, on the other hand, are difficult to model even in the absence of
baryonic feedback (Schneider:2016:JCAP:). In Appendix. A we show how the
matter power suppression changes when changing the resolution of the volume of
the simulation boxes. When comparing with the original IllustrisTNG boxes, we
find that while the box sizes do not change the measured power suppression
significantly, the resolution of the boxes has a non-negligible impact. This
is expected since the physical effect of feedback mechanisms depend on the
resolution of the simulations. Note that the error-bars presented in Fig. 1
will also depend on the resolution and size of the simulation box as well as
on the feedback parameter values assumed. We defer a detailed study of the
covariance dependence on the simulation properties to a future study.

###  2.3 Measuring gas profiles around halos

We use 3D grids of various fields (e.g. gas density and pressure) made
available by the CAMELS team to extract the profiles of these fields around
dark matter halos. The grids are generated with resolution of 0.05 Mpc/$h$.
Following vanDaalen:2020, we define $f_{b}$ as $(M_{\rm gas}+M_{\rm
stars})/M_{\rm total}$, where $M_{\rm gas}$, $M_{\rm stars}$ and $M_{\rm
total}$ are the mass in gas, stars and all the components, respectively. The
gas mass is computed by integrating the gas number density profile around each
halo. We typically measure $f_{b}$ within the spherical overdensity radius
$r_{\rm 500c}$.555We define spherical overdensity radius ($r_{\Delta c}$,
where $\Delta=200/500$$\Delta=200,500$) and overdensity mass ($M_{\Delta c}$)
such that the mean density within $r_{\Delta}$ is $\Delta$ times the critical
density $\rho_{\rm crit}$, $M_{\Delta}=\Delta\frac{4}{3}\pi
r^{3}_{\Delta}\rho_{\rm crit}$.

The SZ effect is sensitive to the electron pressure. We compute the electron
pressure profiles, $P_{e}$, using $P_{e}=2(X_{\rm H}+1)/(5X_{\rm H}+3)P_{\rm
th}$, where $P_{\rm th}$ is the total thermal pressure, and $X_{\rm H}=0.76$
is the primordial hydrogen fraction. Given the electron pressure profile, we
measure the integrated SZ signal within $r_{\rm 500c}$ as:

$Y_{\rm 500c}=\frac{\sigma_{\rm T}}{m_{e}c^{2}}\int_{0}^{r_{\rm 500c}}4\pi
r^{2}\,P_{e}(r)\,dr,$ (1)

where, $\sigma_{\rm T}$ is the Thomson scattering cross-section, $m_{e}$ is
the electron mass and $c$ is the speed of light.

We normalize the SZ observables by the self-similar expectation
(Battaglia:2012:ApJ:a),666Note that we use spherical overdensity mass
corresponding to $\Delta=500$ and hence adjust the coefficients accordingly,
while keeping other approximations used in their derivations as the same.

$Y^{\rm
SS}=131.7h^{-1}_{70}\,\bigg{(}\frac{M_{500c}}{10^{15}h^{-1}_{70}M_{\odot}}\bigg{)}^{5/3}\frac{\Omega_{b}}{0.043}\frac{0.25}{\Omega_{m}}\,{\rm
kpc^{2}},$ (2)

where, $M_{200c}$ is mass inside the $r_{200c}$ and $h_{70}=h/0.7$. This
calculation, which scales as $M^{5/3}$, assumes hydrostatic equilibrium and
that the baryon fraction is equal to cosmic baryonic fraction. Hence
deviations from this self-similar scaling provide a probe of the effects of
baryonic feedback. Our final SZ observable is defined as $Y_{500c}/Y^{\rm
SS}$. On the other hand, that the amplitude of the pressure profile
approximately scales as $M^{2/3}$. Therefore, when considering the pressure
profile as the observable, we factor out a $M^{2/3}$ scaling.

##  3 Results I: Simulations

###  3.1 Inferring feedback parameters from $f_{b}$ and $y$

We first consider how the halo $Y$ signal can be used to constrain the
parameters describing the subgrid physics models. This question has been
previously investigated using the CAMELS simulations by Moser:2021 and
Wadekar:2022. The rest of our analysis will focus on constraining changes to
the power spectrum and bispectrum, and our intention here is mainly to provide
a basis of comparison for those results. Similar to Wadekar:2022, we treat the
mean $\log Y$ value in two mass bins ($10^{12}<M(M_{\odot}/h)<5\times 10^{12}$
and $5\times 10^{12}<M(M_{\odot}/h)<10^{14}$) as our observable, which we
refer to as $\vec{d}$. Here and throughout our investigations with CAMELS we
ignore the contributions of measurement uncertainty since our intention is
mainly to assess the information content of the SZ signals. We therefore use
the CV simulations to determine the covariance, $\mathbf{C}$, of the
$\vec{d}$. Note that the level of cosmic variance will depend on the volume
probed, and can be quite large for the CAMELS simulations. The Fisher matrix,
$F_{ij}$ is then given by

$F_{ij}=\frac{\partial\vec{d}^{T}}{\partial\theta_{i}}\mathbf{C}^{-1}\frac{\partial\vec{d}}{\partial\theta_{i}},$
(3)

where $\theta_{i}$ refers to the $i$th parameter value.

Calculation of the derivatives $\partial\vec{d}/\partial\theta_{i}$ is
complicated by the large amount of stochasticity between the CAMELS
simulations. To perform the derivative calculation, we use a spline radial
basis function interpolation method based on Moser:2022; Cromer:2022. We show
an example of the derivative calculation in Appendix E. We additionally assume
a prior of $\sigma(\ln A)=1$ on a parameter $p$ of $\sigma(\ln p)=1$ on the
feedback parameters and $\sigma(A)=1$ $\sigma(p)=1$ on the cosmological
parameters.

The parameter constraints corresponding to our calculated Fisher matrix are
shown in Fig. 2. We show results only for $\Omega_{m}$, $A_{\rm SN1}$ and
$A_{\rm AGN2}$, but additionally marginalize over $\sigma_{8}$, $A_{\rm SN2}$
and $A_{\rm AGN1}$. The degeneracy directions seen in our results are
consistent with those in Wadekar:2022. We we find a weaker constraint on
$A_{\rm AGN2}$, likely owing to the large sample variance contribution to our
calculation.

It is clear from Fig. 2 that the marginalized constraints on the feedback
parameters are weak. If information about $\Omega_{m}$ is not used, we
effectively have no information about the feedback parameters. Even when
$\Omega_{m}$ is fixed, the constraints on the feedback parameters are not very
precise. This finding is consistent with Wadekar:2022, for which measurement
uncertainty was the main source of variance rather than sample variance. Part
of the reason for the poor constraints is the degeneracy between the AGN and
SN parameters. Degeneracies between the impacts of feedback parameters and
cosmology on $Y$, as well as the potentially complex relation between feedback
parameters and the changing matter distribution motivate us to consider
instead direct inference of changes to the statistics of the matter
distribution from the $Y$ observables. However, note that the conclusions here
will depend upon the simulation volume, which would change the covariance and
will capture the effects like super sample covariance.

Figure 2: Forecast parameter constraints on the feedback parameters when $\log
Y$ in two halo mass bins is treated as the observable. We assume that the only
contribution to the variance of this observable is sample variance coming from
the finite volume of the CAMELS simulations.

###  3.2 $f_{b}$ and $y$ as probes of baryonic effects on the matter power
spectrum

As discussed above, vanDaalen:2020 observed a tight correlation between
suppression of the matter power spectrum and the baryon fraction, $f_{b}$, in
halos with $6\times 10^{13}M_{\odot}\lesssim M_{500c}\lesssim
10^{14}\,M_{\odot}$. That relation was found to hold regardless of the details
of the feedback implementation, suggesting that by measuring $f_{b}$ in high-
mass halos, one could robustly infer the impact of baryonic feedback on the
power spectrum. We begin by investigating the connection between matter power
spectrum suppression and $f_{b}$ in low-mass, $M\sim 10^{13}\,M_{\odot}$,
halos. We also consider a wider range of feedback models than vanDaalen:2020,
including the SIMBA and Astrid models.

Figure 3: We show the relation between matter power suppression at $k=2h/{\rm
Mpc}$ and baryon fraction of halos in mass range
$10^{13}<M\,(M_{\odot}/h)<10^{14}$ in the SIMBA simulation suite. In each of
six panels, the points are colored corresponding to the parameter value given
in the associated colorbar.

Fig. 3 shows the impact of cosmological and feedback parameters on the
relationship between the power spectrum suppression ($\Delta P/P_{\rm DMO}$)
and the ratio $Y_{\rm 500c}/Y^{\rm SS}$ for the SIMBA simulations. Each point
corresponds to a single simulation, taking the average over all halos with
$10^{13}<M(M_{\odot}/h)<10^{14}$ when computing $Y_{\rm 500c}/Y^{\rm SS}$. We
observe that the largest $\Delta P/P_{\rm DMO}$ occurs when $A_{\rm AGN2}$ is
large. This is caused by powerful AGN feedback ejecting gas from halos,
leading to a significant reduction in the matter power spectrum, as described
by e.g. vanDaalen:2020. For SIMBA, the parameter $A_{\rm AGN2}$ controls the
velocity of the ejected gas, with higher velocities (i.e. higher $A_{\rm
AGN2}$) leading to more ejected gas. On the other hand, when $A_{\rm SN2}$ is
large, $\Delta P/P_{\rm DMO}$ is small. This is because supernovae prevent the
gas from accreting onto the supermassive blackhole, reducing the strength of
efficient supernovae feedback prevents the formation of massive galaxies which
host AGN and hences reduces the strength of the AGN feedback. The parameter
$A_{\rm AGN1}$, on the other hand, controls the radiative quasar mode of
feedback, which has slower gas outflows and thus a smaller impact on the
matter distribution.

It is also clear from Fig. 3 that increasing $\Omega_{\rm m}$ reduces $|\Delta
P/P_{\rm DMO}|$, relatively independently of the other parameters. By
increasing $\Omega_{m}$, the ratio $\Omega_{b}/\Omega_{m}$ decreases, meaning
that halos of a given mass have fewer baryons, and the impact of feedback is
therefore reduced. We propose a very simple toy model for this effect in §3.3.

The impact of $\sigma_{8}$ in Fig. 3 is less clear. For halos in the mass
range shown, we find that increasing $\sigma_{8}$ leads to a roughly monotonic
decrease in $f_{b}$, presumably because higher $\sigma_{8}$ means that there
are more halos amongst which the same amount of baryons must be distributed.
This effect would not occur for cluster-scale halos because these are rare and
large enough to gravitationally dominate their local environments, giving them
$f_{b}\sim\Omega_{b}/\Omega_{m}$, regardless of $\sigma_{8}$. In any case, no
clear trend with $\sigma_{8}$ is seen in Fig. 3 because $\sigma_{8}$ does not
correlate strongly with $\Delta P/P_{\rm DMO}$.

Fig. 4 shows the relationship between $\Delta P/P_{\rm DMO}$ and $f_{b}$ or
$Y_{500}$ in different halo mass bins and for different feedback models,
colored by the value of $A_{AGN2}$. As in Fig. 3, each point represents an
average over all halos in the indicated mass range for a particular CAMELS
simulation (i.e. at fixed values of cosmological and feedback parameters).
Note that the meaning of $A_{\rm AGN2}$ is not exactly the same across the
different feedback models, as noted in §2. For TNG and SIMBA we expect
increasing $A_{\rm AGN2}$ to lead to stronger AGN feedback driving more gas
out of the halos, leading to more power suppression. However, for Astrid,
increasing $A_{\rm AGN2}$ parameter would more strongly regulate and suppress
the black hole growth in the box. This drastically reduces the number of high
mass black holes and hence effectively reducing the amount of jet feedback
that can push the gas out of the halos, leading to less matter power
suppression. Therefore, in Fig. 4, we redefine the $A_{\rm AGN2}$ parameter
for Astrid to be $1/A_{\rm AGN2}$ when plotting.

For the highest mass bin ($M>5\times 10^{13}M_{\odot}/h$, rightmost column)
our results are in agreement with vanDaalen:2020: we find that there is a
fairly tight relation between between $f_{b}/(\Omega_{b}/\Omega_{m})$ and the
matter power suppression. This relation is roughly consistent across different
feedback subgrid models, although the different models appear to populate
different parts of this relation. Moreover, varying $A_{AGN2}$ appears to move
points along this relation, rather than broadening the relation. This is in
contrast to $\Omega_{m}$, which as shown in Fig. 3, tends to move simulations
in the direction orthogonal to the narrow $f_{b}$-$\Delta P/P_{\rm DMO}$
locus. For this reason, and given current constraints on $\Omega_{m}$, we
restrict this plot to simulations with $0.2<\Omega_{m}<0.4$. The dashed curves
shown in Fig. 3 correspond to the toy model discussed in §3.3.

At low halo mass, the relation between $f_{b}/(\Omega_{b}/\Omega_{m})$ and
$\Delta P/P_{\rm DMO}$ appears similar to the high mass bin, although it is
somewhat flatter at high $f_{b}$, and somewhat steeper at low $f_{b}$. Again
the results are fairly consistent cross the different feedback prescriptions,
although points with high $f_{b}/(\Omega_{b}/\Omega_{m})$ are largely absent
for SIMBA. This is largely because the feedback mechanisms are highly
efficient in SIMBA, driving the gas out of their parent halos.

The relationships between $Y$ and $\Delta P/P$ appear quite similar to those
between $\Delta P/P$ and $f_{b}$$f_{b}/(\Omega_{b}/\Omega_{m})$. This is not
too surprising because $Y$ is sensitive to the gas density, which dominates
$f_{b}$$f_{b}/(\Omega_{b}/\Omega_{m})$. However, $Y$ is also sensitive to the
gas temperature. Our results suggest that variations in the gas temperature
are not significantly impacting the $Y$-$\Delta P/P$ relation. These results
suggest the possibility of using the tSZ signal rather than $f_{b}$
$f_{b}/(\Omega_{b}/\Omega_{m})$ to infer the impact of feedback on the matter
distribution. This will be the focus of the remainder of the paper.

Figure 4: Impact of baryonic physics on the matter power spectrum at
$k=2h/{\rm Mpc}$ for the Illustris, SIMBA and Astrid simulations (top, middle,
and bottom rows). Each point corresponds to an average across halos in the
indicated mass ranges in a different CAMELS simulation. We restrict the figure
to simulations that have $0.2<\Omega_{\rm m}<0.4$. The dashed curves
illustrate the behavior of the model described in §3.3 in the regime that the
radius to which gas is ejected by AGN is larger than the halo radius, and
larger than $2\pi/k$. Note that for Astrid simulations (marked by asterisk),
we take inverse of $A_{\rm AGN2}$ parameter while plotting as described in the
main text.

Fig. 5 shows the same quantities as Fig. 4, but now for a fixed halo mass
range ($10^{13}<M/(M_{\odot}/h)<10^{14}$), fixed subgrid prescription
(Astrid), and varying values of $k$. We find roughly similar results when
using the different subgrid physics prescriptions. At low $k$, we find that
there is a regime at high $Y$ $f_{b}/(\Omega_{b}/\Omega_{m})$ for which
$\Delta P/P_{\rm DMO}$ changes negligibly. It is only when $Y$
$f_{b}/(\Omega_{b}/\Omega_{m})$ becomes very low that $\Delta P/P_{\rm DMO}$
begins to change. On the other hand, at high $k$, there is a near-linear
relation between $Y$ $f_{b}/(\Omega_{b}/\Omega_{m})$ and $\Delta P/P$.

Figure 5: Similar to Fig. 4, but for different values of $k$. For simplicity,
we show only the Astrid simulations for halos in the mass range
$10^{13}<M(M_{\odot}/h)<10^{14}$. The dashed curves illustrate the behavior of
the model described in §3.3 in the regime that the radius to which gas is
ejected by AGN is larger than the halo radius, and larger than $2\pi/k$. As
expected, this model performs best in the limit of high $k$ and large halo
mass.

###  3.3 A toy model for power suppression

We now describe a simple model for the effects of feedback on the relation
between $f_{b}$ or $Y$ and $\Delta P/P$ that explains some of the features
seen in Figs. 3, 4 and 5. Following expectations from the literature, we
assume in this model that it is ejection of gas from halos by AGN feedback
that is responsible for changes to the matter power spectrum. SN feedback, on
the other hand, prevents gas from accreting onto the SMBH, and therefore
reduces the impact of AGN feedback. This scenario is consistent with the fact
that at high SN feedback, we see that $\Delta P/P_{\rm DMO}$ goes to zero
(second panel from the bottom in Fig. 3).

We identify three relevant scales: (1) the halo radius, $R_{h}$, (2) the
distance to which gas is ejected by the AGN, $R_{\rm ej}$, and (3) the scale
at which the power spectrum is measured, $2\pi/k$. If $R_{\rm ej}<2\pi/k$,
then there will be no impact on $\Delta P$ at $k$: this corresponds to a
rearrangement of the matter distribution on scales below where we measure the
power spectrum. If, on the other hand, $R_{\rm ej}<R_{h}$, then there will be
no impact on $f_{b}$ or $Y$, since the gas is not ejected out of the halo.
Thus, we can consider four different regimes:

  * •

Regime 1: $R_{\rm ej}<R_{h}$ and $R_{\rm ej}<2\pi/k$. In this regime, changes
to the feedback parameters have no impact on $f_{b}$ or $\Delta P$.

  * •

Regime 2: $R_{\rm ej}>R_{h}$ and $R_{\rm ej}<2\pi/k$. In this regime, changes
to the feedback parameters result in movement along the $f_{b}$ or $Y$ axis
without changing $\Delta P$. Gas is being removed from the halo, but the
resultant changes to the matter distribution are below the scale at which we
measure the power spectrum. Note that Regime 2 cannot occur when
$R_{h}>2\pi/k$ (i.e. high-mass halos at large $k$).

  * •

Regime 3: $R_{\rm ej}>R_{h}$ and $R_{\rm ej}>2\pi/k$. In this regime, changing
the feedback amplitude directly changes the amount of gas ejected from halos
as well as $\Delta P/P_{\rm DMO}$.

  * •

Regime 4: $R_{\rm ej}<R_{h}$ and $R_{\rm ej}>2\pi/k$. In this regime, gas is
not ejected out of the halo, so $f_{b}$ and $Y$ should not change. In
principle, the redistribution of gas within the halo could lead to changes in
$\Delta P/P_{\rm DMO}$. However, as we discuss below, this does not seem to
happen in practice.

Let us now consider the behavior of $\Delta P/P_{\rm DMO}$ and $f_{b}$ or $Y$
as the feedback parameters are varied in Regime 3. A halo of mass $M$ is
associated with an overdensity $\delta_{m}$ in the absence of feedback, which
is changed to $\delta^{\prime}_{m}$ due to ejection of baryons as a result of
feedback. In Regime 3, some amount of gas, $M_{\rm ej}$, is completely removed
from the halo. This changes the size of the overdensity associated with the
halo to

$\displaystyle\frac{\delta_{m}^{\prime}}{\delta_{m}}$ $\displaystyle=$
$\displaystyle 1-\frac{M_{\rm ej}}{M}.$ (4)

The change to the power spectrum is then

$\displaystyle\frac{\Delta P}{P_{\rm DMO}}$ $\displaystyle\sim$
$\displaystyle\left(\frac{\delta_{m}^{\prime}}{\delta_{m}}\right)^{2}-1\approx-2\frac{M_{\rm
ej}}{M},$ (5)

where we have assumed that $M_{\rm ej}$ is small compared to $M$. We have
ignored the $k$ dependence here, but in Regime 3, the ejection radius is
larger than the scale of interest, so the calculated $\Delta P/P_{\rm DMO}$
should apply across a range of $k$ in this regime.

The ejected gas mass can be related to the gas mass in the absence of
feedback. We write the gas mass in the absence of feedback as
$f_{c}(\Omega_{b}/\Omega_{m})M$, where $f_{c}$ encapsulates non-feedback
processes that result the halo having less than the cosmic baryon fraction. We
then have

$\displaystyle M_{\rm ej}$ $\displaystyle=$ $\displaystyle
f_{c}(\Omega_{b}/\Omega_{m})M-f_{b}M-M_{0},$ (6)

where $M_{0}$ is the mass that has been removed from the gaseous halo, but
that does not change the power spectrum, e.g. the conversion of gas to stars.
Substituting into Eq. 5, we have

$\displaystyle\frac{\Delta P}{P_{\rm
DMO}}=-2\frac{f_{c}\Omega_{b}}{\Omega_{m}}\left(1-\frac{f_{b}\Omega_{m}}{f_{c}\Omega_{b}}-\frac{\Omega_{m}M_{0}}{f_{c}\Omega_{b}M}\right).$
(7)

In other words, for Regime 3, we find a linear relation between $\Delta
P/P_{\rm DMO}$ and $f_{b}\Omega_{m}/\Omega_{b}$. For high mass halos, we
should have $f_{c}\approx 1$ and $M_{0}/M\approx 0$. In this limit, the
relationship between $f_{b}$ and $\Delta P/P$ becomes

$\displaystyle\frac{\Delta P}{P_{\rm
DMO}}=-2\frac{\Omega_{b}}{\Omega_{m}}\left(1-\frac{f_{b}\Omega_{m}}{\Omega_{b}}\right),$
(8)

which is linear between $(\Delta P/P_{\rm
DMO},f_{b}\Omega_{m}/\Omega_{b})=(-2\Omega_{b}/\Omega_{m},0)$ and $(\Delta
P/P_{\rm DMO},f_{b}\Omega_{m}/\Omega_{b})=(0,1)$. We show this relation as the
dashed line in the $f_{b}$ columns of Figs. 4.

We can repeat the above argument for $Y$. Unlike the case with $f_{b}$,
processes other than the removal of gas may reduce $Y$; these include, e.g.,
changes to the gas temperature in the absence of AGN feedback, or nonthermal
pressure support. We account for these with a term $Y_{0}$, defined such that
when $M_{\rm ej}=M_{0}=0$, we have
$Y+Y_{0}=f_{c}(\Omega_{b}/\Omega_{m})MT/\alpha$, where we have assumed
constant gas temperature, $T$, and $\alpha$ is a dimensionful constant of
proportionality. We then have

$\displaystyle\frac{\alpha(Y+Y_{0})}{T}=f_{c}(\Omega_{b}/\Omega_{m})M-M_{\rm
ej}-M_{0}.$ (9)

Substituting the above equation into Eq. 5 we have

$\displaystyle\frac{\Delta P}{P_{\rm DMO}}$ $\displaystyle=$
$\displaystyle-2\frac{f_{c}\Omega_{b}}{\Omega_{m}}\left(1-\frac{\alpha(Y+Y_{0})\Omega_{m}}{f_{c}TM\Omega_{b}}-\frac{\Omega_{m}M_{0}}{f_{c}\Omega_{b}M}\right).$
(10)

Following Eq. 2, we define the self-similar value of $Y$, $Y^{\rm SS}$, via

$\displaystyle\alpha Y^{\rm SS}/T=(\Omega_{b}/\Omega_{m})M,$ (11)

leading to

$\displaystyle\frac{\Delta P}{P_{\rm DMO}}$ $\displaystyle=$
$\displaystyle-2\frac{f_{c}\Omega_{b}}{\Omega_{m}}\left(1-\frac{(Y+Y_{0})}{f_{c}Y^{\rm
SS}}-\frac{\Omega_{m}M_{0}}{f_{c}\Omega_{b}M}\right).$ (12)

Again taking the limit that $f_{c}\approx 1$ and $M_{0}/M\approx 0$, we have

$\displaystyle\frac{\Delta P}{P_{\rm DMO}}$ $\displaystyle=$
$\displaystyle-2\frac{\Omega_{b}}{\Omega_{m}}\left(1-\frac{(Y+Y_{0})}{Y^{\rm
SS}}\right).$ (13)

Thus, we see that in Regime 3, the relation between $Y/Y^{SS}$ and $\Delta
P/P_{\rm DMO}$ is linear. The $Y/Y^{SS}$ columns of Figs. 4 show this
relationship, assuming $Y_{0}=0$.

In summary, we interpret the results of Figs. 4 and 5 in the following way.
Starting at low feedback amplitude, we are initially in Regime 1. In this
regime, the simulations cluster around $f_{b}f_{c}\Omega_{m}/\Omega_{b}\approx
1$ (or $Y\approx Y_{0}$) and $\Delta P/P\approx 0$ since changing the feedback
parameters in this regime does not impact $f_{b}$ or $\Delta P/P$. For high
mass halos, we have $f_{c}\approx 1$ and $Y_{0}\approx 0$ (although SIMBA
appears to have $Y_{0}>0$, even at high mass); for low mass halos, $f_{c}<1$
and $Y_{0}>0$. As we increase the AGN feedback amplitude, the behavior is
different depending on halo mass and $k$:

  * •

For low halo masses or low $k$, increasing the AGN feedback amplitude leads
the simulations into Regime 2. Increasing the feedback amplitude in this
regime moves points to lower $Y/Y^{\rm SS}$ (or $f_{b}\Omega_{m}/\Omega_{b}$)
without significantly impacting $\Delta P/P_{\rm DMO}$. Eventually, when the
feedback amplitude is sufficiently strong, these halos enter Regime 3, and we
see a roughly linear decline in $\Delta P/P_{\rm DMO}$ with decreasing
$Y/Y^{\rm SS}$ (or $f_{b}\Omega_{m}/\Omega_{b}$), as discussed above.

  * •

For high mass halos and high $k$, we never enter Regime 2 since it is not
possible to have $R_{\rm ej}>R_{h}$ and $R_{\rm ej}<2\pi/k$ when $R_{h}$ is
very large. In this case, we eventually enter Regime 3, leading to a linear
trend of decreasing $\Delta P/P_{\rm DMO}$ with decreasing $Y/Y^{\rm SS}$ or
$f_{b}\Omega_{m}/\Omega_{b}$, as predicted by the above discussion. This
behavior is especially clear in Fig. 5: at high $k$, the trend closely follows
the predicted linear relation. At low $k$, on the other hand, we see a more
prominent Regime 2 region. The transition between these two regimes is
expected to occur when $k\sim 2\pi/R_{h}$, which is roughly $5\,h^{-1}{\rm
Mpc}$ for the halo mass regime shown in the figure. This expectation is
roughly confirmed in the figure.

Interestingly, we never see Regime 4 behavior: when the halo mass is large and
$k$ is large, we do not see rapid changes in $\Delta P/P$ with little change
to $f_{b}$ and $Y$. This could be because this regime corresponds to movement
of the gas entirely within the halo. If the gas has time to re-equilibrate, it
makes sense that we would see little change to $\Delta P/P$.

###  3.4 Predicting the power spectrum suppression from the halo observables

While the toy model described above roughly captures the trends between $Y$
(or $f_{b}$) and $\Delta P/P_{\rm DMO}$, it of course does not capture all of
the physics associated with feedback. It is also clear that there is
significant scatter in the relationships between observable quantities and
$\Delta P$. It is possible that this scatter is reduced in some higher
dimensional space that includes more observables. To address both of these
issues, we now train statistical models to learn the relationships between
observable quantities and $\Delta P/P_{\rm DMO}$. We will focus on results
obtained with random forest regression (Breiman2001). We have also tried using
neural networks to infer these relationships, but have not found any
significant improvement with respect to the random forest results, presumably
because the space is low-dimensional (i.e. we consider at most about five
observable quantities at a time). We leave a detailed comparison with other
decision tree based approaches, such as gradient boosted tree
(Friedman_boosted_tree:01) to a future study.

Figure 6: We show the results of the random forest (RF) regressor at
predicting $\Delta P/P$ for variations in the test and train sample using the
different subgrid physics models (TNG, SIMBA, Astrid) on LH suite of
simulations. The RF model is trained on $f_{b}$ for halos with $5\times
10^{12}<M(M_{\odot}/h)<10^{14}$ as well as $\Omega_{m}$. The errorbars
correspond to the 16-th and 84-th percentile error on the recovered
probability density function (PDF) of the power suppression on the test set of
LH simulations, where as the marker corresponds to the peak of the PDF. The
gray band corresponds to the expected $1\sigma$ error from the CV simulations
on the power suppression. We find that we can predict the power suppression
robustly when the training and test simulation codes are the same. At high
$k$, training on two simulations and predicting on a third works well. At low
$k$, however, the RF model becomes biased when testing on a third simulation.
For the rest of the paper, we present results trained on all three
simulations.

We train a random forest model to go from observable quantities (e.g.
$f_{b}/(\Omega_{b}/\Omega_{m})$ and $Y_{500}/Y^{SS}$) to a prediction for
$\Delta P/P_{\rm DMO}$ at multiple $k$ values. The random forest model uses
100 trees with a ${\rm max}_{\rm depth}=10$.777We use a publicly available
code: https://scikit-
learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html.
We also verified that our conclusions are robust to changing the settings of
the random forest. In this section we analyze the halos in the mass bin
$5\times 10^{12}<M_{\rm halo}(M_{\odot}/h)<10^{14}$ but we also show the
results for halos with lower masses in Appendix C. We also consider supplying
the values of $\Omega_{\rm m}$ as input to the random forest, since it can be
constrained precisely through other observations (e.g. CMB observations), and
as we showed in §3.2, the cosmological parameters can impact the
observables.888One might worry that using cosmological information to
constrain $\Delta P/P_{\rm DMO}$ defeats the whole purpose of constraining
$\Delta P/P_{\rm DMO}$ in order to improve cosmological constraints. However,
observations, such as those of CMB primary anisotropies, already provide
precise constraints on the matter density without using information in the
small-scale matter distribution.

Ultimately, we are interested in making predictions for $\Delta P/P_{\rm DMO}$
using observable quantities. However, the sample variance in the CAMELS
simulations limits the precision with which we can measure $\Delta P/P_{\rm
DMO}$. It is not possible to predict $\Delta P/P_{\rm DMO}$ to better than
this precision. We will therefore normalize the uncertainties in the RF
predictions by the cosmic variance error. In order to obtain the uncertainty
in the prediction, we randomly split the data into 70% training and 30% test
set. After training the RF regressor using the training set and a given
observable, we make compute the 16th and 84th percentile of the distribution
of prediction errors evaluated on the test set. This constitutes our
assessment of prediction uncertainty.

Figure 7:  Similar to Fig. 6 but when training the RF model on different
observables from all three simulations (TNG, SIMBA and Astrid) to predict
$\Delta P/P$ of a random subset of the the three simulations not used in
training. We find that jointly training on the deviation of the integrated SZ
profile from the self-similar expectation, $Y_{500c}/Y^{\rm SS}$ and
$\Omega_{m}$ results in inference of power suppression that is comparable to
cosmic variance errors, with small improvements when additionally adding the
baryon fraction ($f_{b}$) of halos in the above mass range.  Figure 8:  Same
as Fig. 7 but when using the full pressure profile and electron number density
profiles instead of the integrated quantities. We again find that with
pressure profile and $\Omega_{m}$ information we can recover robust and
precise constraints on the matter power suppression.

Fig. 6 shows the accuracy of the RF predictions for $\Delta P/P$ when trained
on $f_{b}$ and $\Omega_{m}$, normalized to the sample variance error in
$\Delta P/P$. As we will show later in this section, this combination of
inputs results in precise constraints on the matter power suppression. We show
the results of training and testing on a single simulation suite, and also the
results of training/testing across different simulation suites. It is clear
that when training and testing on the same simulation suite, the RF learns a
model that comes close to the best possible uncertainty on $\Delta P/P_{\rm
DMO}$ (i.e. cosmic variance). When training on one or two simulation suites
and testing another, however, the predictions show bias at low $k$. This
suggests that the model learned from one simulation does not generalize very
well to another in this regime. This result is somewhat different from the
findings of vanDaalen:2020, where it was found that the relationship between
$f_{b}$ and $\Delta P/P_{\rm DMO}$ did generalize to different simulations.
This difference may result from the fact that we are considering a wider range
of feedback prescriptions than in vanDaalen:2020, as well as considering
significant variations in cosmological parameters.

Fig. 6 also shows the results of testing and training on all three simulations
(black points with errorbars). Encouragingly, we find that in this case, the
predictions are of comparable accuracy to those obtained from training and
predicting on the same simulation suite. This suggests that there is a general
relationship across all feedback models that can be learned to go from
$\Omega_{m}$ and $f_{b}$ to $\Delta P/P_{\rm DMO}$. Henceforth, we will show
results trained on all simulation suites and tested on all simulations suites.
Of course, this result does not imply that our results will generalize to some
completely different feedback prescription.

In Fig. 7 we show the results of training the random forest on different
combinations of $f_{b}$, $Y_{500}$ and $\Omega_{m}$. Consistent with the
findings of vanDaalen:2020, we find that $f_{b}/(\Omega_{b}/\Omega_{m})$
results in robust constraints on the matter power suppression (blue points
with errors). These constraints come close to the cosmic variance limit across
a wide range of $k$.

We additionally find that providing $f_{b}$ and $\Omega_{m}$ as separate
inputs to the RF improves on the combination $f_{b}/(\Omega_{b}/\Omega_{m})$,
yielding smaller variance in the predicted $\Delta P/P_{\rm DMO}$, with the
most dramatic largest improvement at small scales. This is not surprising
given the predictions of our simple model, for which it is clear that $\Delta
P/P$ can be impacted by both $\Omega_{m}$ and $f_{b}/(\Omega_{b}/\Omega_{b})$
independently. Similarly, it is clear from Fig. 3 that changing $\Omega_{m}$
changes the relationship between $\Delta P/P$ and the halo gas-derived
quantities (like $Y$ and $f_{b}$).

We next consider a model trained on $Y_{500c}/Y^{SS}$ (orange points in Fig.
7). This model yields reasonable predictions for $\Delta P/P$, although not
quite as good as the model trained on $f_{b}/(\Omega_{b}/\Omega_{m})$. The
$Y/Y^{SS}$ model yields somewhat larger errorbars, and the distribution of
$\Delta P/P$ predictions is highly asymmetric. When we train the RF model
jointly on $Y_{500c}/Y^{SS}$ and $\Omega_{m}$ (green points), we find that the
predictions improve considerably, particularly at high $k$. In this case, the
predictions are typically symmetric around the true $\Delta P/P$, have smaller
uncertainty compared to the model trained on $f_{b}/(\Omega_{b}/\Omega_{m})$,
and comparable uncertainty to the model trained on
$\\{f_{b}/(\Omega_{b}/\Omega_{m})$,$\Omega_{m}\\}$. We thus conclude that when
combined with matter density information, $Y/Y^{\rm SS}$ provides a powerful
probe of baryonic effects on the matter power spectrum.

Above we have considered the integrated tSZ signal from halos, $Y_{500c}$.
Measurements in data, however, can potentially probe the tSZ profiles rather
than only the integrated tSZ signal (although the instrumental resolution may
limit the extent to which this is possible). In Fig. 8 we consider RF models
trained on the full profiles instead of just the integrated quantities. The
electron pressure and number density profiles are measured in eight
logarithmically spaced bins between $0.1<r/r_{200c}<1$. We find that while the
ratio $P_{e}(r)/P^{\rm SS}$ results in robust constraints, jointly providing
the information on $\Omega_{m}$ makes them more precise. Similar to the
integrated profile case, we find that additionally providing the electron
density profile information only marginally improves the constraints. We also
show the results when jointly using the measured pressure profiles for both
the low and high mass halos to infer the matter power suppression. We find
that this leads to only a marginal improvements in the constraints. This would
suggest that the deviation of the thermal pressure from the expected self-
similar relation already captures the strength of feedback adequately with
minimal improvements with adding information from lower mass halos.

Note that we have input the 3D pressure and electron density profiles in this
case. Even though observed SZ maps are projected quantities, we can infer the
3D pressure profiles from the model used to analyze the projected
correlations.

Figure 9: Same as Fig. 7, but for bispectrum suppression for equilateral
triangle configurations at different scales. We find that having pressure
profile information results in unbiased constraints here as well.

###  3.5 Predicting the bispectrum suppression with $f_{b}$ and electron
pressure

In Fig. 9, we test our methodology on bi-spectrum suppressionbispectrum
suppression, $\Delta B(k)/B(k)$. Similar to the matter power spectrum, we
train and test our model on a combination of the three simulations. We train
and test on equilateral triangle bi-spectrum bispectrum configurations with
different scales $k$. We again see that information about the electron
pressure and $\Omega_{m}$ results in a precise and unbiased constraints on the
bispectrum suppression. The constraints improve as we go to the small scales.
In Appendix D we show similar methodology applied to squeezed bi-spectrum
bispectrum configurations.

However several caveats are important to consider here. The bispectrum is more
sensitive to the higher mass ($M>5\times 10^{13}M_{\odot}/h$) halos, which are
missing from the CAMELS simulations. This results in the sensitivity of the
bispectrum suppression to also be biased in these small boxes. We also note
that there additionally can be more important dependencies on the simulation
resolution which are beyond the scope of this study. However, we expect the
qualitative methodology applied here to remain valid with the set of larger
simulations. Finally, there would be some degeneracy between the power
spectrum suppression and the bispectrum suppression as they both stem from
same underlying physics which we defer to a future study.

##  4 Results II: ACTxDES measurements and forecast

Figure 10: Constraints on the matter power suppression obtained from the
inferred $Y_{\rm 500c}/Y^{\rm SS}$ (and fixing $\Omega_{\rm m}=0.3$) from the
shear-$y$ correlations obtained from DESxACT analysis (Pandey:2022). We also
show the expected improvements from future halo-$y$ correlations from DESIxSO
using the constraints in (Pandey:2020). We compare these to the inferred
constraints obtained using cosmic shear (Chen:2022:MNRAS:) and additionally
including X-ray and kSZ data (Schneider:2022:MNRAS:). We also compare with the
results from larger simulations; OWLS (Schaye:2010:MNRAS:), BAHAMAS
(McCarthy:2017:MNRAS:) and TNG-300 (Springel:2018:MNRAS:).

Our analysis above has resulted in a statistical model (random forest) that
predicts $\Delta P/P_{\rm DMO}$ (and the corresponding uncertainties) given
values of $Y_{500c}$ for low-mass halos and $\Omega_{m}$. This model is robust
to significant variations in the feedback prescription, at least across the
SIMBA, Illustris and Astrid models. We now apply this model to constraints on
$Y_{500c}$ coming from the cross-correlation of galaxy lensing shear with tSZ
maps.

Gatti:2022 and Pandey:2022 measured cross-correlations of DES galaxy lensing
with Compton $y$ maps from a combination of Advanced-ACT
(Madhavacheril:2020:PhRvD:) and Planck data (PlanckCollaboration:2016:A&A:)
over an area of 400 deg. sq. sq. deg. They analyze these cross-correlations
using a halo model framework, where the pressure profile in halos is
parameterized using a generalized Navarro-Frenk-White profile
Navarro:1996:ApJ:; Battaglia:2012:ApJ:b. This pressure profile is described
using four free parameters, allowing for scaling with mass, redshift and
distance from halo center. A tomographic analysis of this shear-y correlation
constrains these parameter and hence the pressure profile of halos in a wide
range of masses. The constraints on these parameterized profiles can be
translated directly into constraints on $Y_{500c}$ for halos in the mass range
that we have used to infer the constraints on matter power suppression from
the trained random forest model as described in previous section. Note that
the shear-$y$ correlation has sensitivity across the range relevant to our
trained model ($M>5\times 10^{12}M_{\odot}/h$), but that the sensitivity is
reduced towards the low end of this range.

Fig. 10 shows the results of feeding the inferred $Y_{\rm 500c}$ constraints
from Pandey:2022 into our random forest model to infer the impact of baryonic
feedback on the matter power spectrum (black points with errorbars). Note that
in this inference we fix the matter density parameter, $\Omega_{m}=0.3$, same
value as used by the CAMELS CV simulations as we use these to estimate the
halo mass function. The shear-tSZ analysis correlation provides constraints on
the parameters of the 3D pressure profiles of halos is a wide range of masses.
We profile of halos and its evolution with mass and redshift. We first use
these parameter constraints to first generate 400 samples of the inferred 3D
profiles of the halos at $z=0$ in 10 logarithmic mass bins in range
$12.7<\log_{10}(M)<14$. Then we perform the volume integral of these profiles
to infer the $Y_{\rm 500c}(M,z)$ . (see Eq. 1). Afterwards, we generate the
stacked normalized integrated pressure for each sample $j$ as: by integrating
over the halo masses as:

$\bigg{\langle}\frac{Y_{\rm 500c}}{Y^{\rm
SS}}\bigg{\rangle}^{j}=\frac{1}{\bar{n}^{j}}\int
dM\bigg{(}\frac{dn}{dM}\bigg{)}^{j}_{\rm
CAMELS}{\color[rgb]{1,0,0}\sout{\frac{Y_{\rm 500c}(M)}{Y^{\rm
SS}}}}{\color[rgb]{0,0,1}\uwave{\frac{Y^{j}_{\rm 500c}(M)}{Y^{\rm SS}}}}$ (14)

where, $\bar{n}^{j}=\int dM(dn/dM)^{j}_{\rm CAMELS}$ and $(dn/dM)^{j}_{\rm
CAMELS}$ is a randomly chosen halo mass function from the CV set of boxes of
TNG, SIMBA or Astrid. This allows for incorporating the uncertainties and the
impact impact (and its corresponding uncertainties) of having a small box size
on the halo mass function. Note that due to small box size, there is deficit
of high mass halos and hence the functional form differs from other fitting
functions in literature, e.g. Tinker:2008:ApJ:. Thereafter, we feed-in these
stacked pressure profiles to the random forest regressor, jointly trained on
TNG, SIMBA and Astrid. For each of the input pressure profile, we recover the
value of matter power suppression, $\Delta P/P_{\rm DMO}$. Finally, in Fig.
10, we plot the mean and the 18th and 64th percentile of the recovered $\Delta
P/P_{\rm DMO}$ distribution . from the 400 samples. We note that our inference
of uncertainties is robust to the number of samples considered.

In the same figure, we also show the constraints from Chen:2022:MNRAS: and
Schneider:2022:MNRAS: obtained using the analysis of complementary datasets.
Chen:2022:MNRAS: analyze the small scale cosmic shear measurements from DES
Year-3 data release using a baryon correction model. Note that in this
analysis, they only use a limited range of cosmologies, particularly
restricting to high $\sigma_{8}$ range due to the emulator calibration.
Moreover they also impose cosmology constraints from the large scale analysis
of the DES data. Note that unlike the procedure presented here, their modeling
and constraints are sensitive to the priors on $\sigma_{8}$. Therefore, their
constraints might be optimistic in this case. Schneider:2022:MNRAS: analyze
the X-ray data (as presented in Giri:2021:JCAP:) and kSZ data from ACT and
SDSS (Schaan:2021:PhRvD:) and the cosmic shear measurement from kids
(Asgari:2021), using another version of baryon correction model. A joint
analysis from these complementary dataset leads to crucial degeneracy breaking
in the parameters. It would be interesting to include the tSZ observations
presented here in the same framework as it can potentially make the
constraints more precise.

Several caveats about our analysis with data are in order. First, the lensing-
SZ correlation is most sensitive to halos in the mass range of $M_{\rm
halo}\geq 10^{13}M_{\odot}/h$. However, our RF model operates on halos with
mass in the range of $5\times 10^{12}\geq M_{\rm halo}\leq
10^{14}M_{\odot}/h$, with limited volume of the simulations restricting the
number of halos above $10^{13}M_{\odot}/h$. We have attempted to account for
this selection effect by using the halo mass function from the CV sims of the
CAMELS simulations when calculating the stacked profile. However, using a
larger volume simulation suite would be a better alternative (also see
discussion in Appendix A). Moreover, the CAMELS simulation suite fix the
$\Omega_{b}$ to a fiducial value. There might be some non-trivial effects on
the inferences when varying that parameter, as that would impact the
distribution of baryons, especially to low mass halos and its interplay with
the changing baryonic feedback.

In order to shift the sensitivity of the data correlations to lower halo
masses, it would be preferable to analyze the galaxy-SZ and halo-SZ
correlations. In Pandey:2020 we forecast the constraints on the inferred 3D
pressure profile from the future halo-SZ correlations using DESI and CMB-S4 SZ
maps for a wide range of halo masses. In Fig. 10 we also show the expected
constraints on the matter power suppression using the halo-SZ correlations
from halos in $M_{h}>5\times 10^{12}M_{\odot}/h$. We again follow the same
methodology as described above to create a stacked normalized integrated
pressure (see Eq. 14). Moreover, we also fix $\Omega=0.3$ to predict the
matter power suppression. Note that we shift the mean value of $\Delta
P/P_{\rm DMO}$ to the recovered value from BAHAMAS high-AGN simulations
(McCarthy:2017:MNRAS:). As we can see in Fig. 10, we can expect to obtain
significantly more precise constraints from these future observations.

##  5 ConclusionConclusions

We have shown that the tSZ signals from low-mass halos contain significant
information about the impacts of baryonic feedback on the small-scale matter
distribution. Using models trained on hydrodynamical simulations with a wide
range of feedback implementations, we demonstrate that information about
baryonic effects on the power spectrum and bispectrum can be robustly
extracted. By applying these same models to measurements with ACT and DES, we
have shown that current tSZ measurements already constrain the impact of
feedback on the matter distribution. Our results suggest that using
simulations to learn the relationship between halo gas observables and
baryonic effects on the matter distribution is a promising way forward for
constraining these effects with data.

Our main findings from our explorations with the CAMELS simulations are the
following:

  * •

In agreement with vanDaalen:2020, we find that baryon fraction in halos
correlates with the power spectrum suppression. We find that the correlation
is especially robust in small scales.

  * •

We find that there can be significant scatter in relationship between baryon
fraction and power spectrum suppression at low halo mass, and that the
relationship varies to some degree with feedback implementation. However, the
bulk trends appear to be consistent regardless of feedback implementation.

  * •

We propose a simple model that qualitatively (and in some cases
quantitatively) captures the broad features in the relationships between
feedback, $\Delta P/P_{\rm DMO}$ (at different values of $k$), and halo gas-
related observables like $f_{b}$ and $Y$ (at different halo masses).

  * •

Despite significant scatter in the relations between $Y$ and $\Delta P/P$ at
low halo mass, we find that simple random forest models yield tight and robust
constraints on $\Delta P/P$ given information about $Y$ in low-mass halos and
$\Omega_{m}$.

  * •

Using the pressure profile instead of just the integrated $Y_{\rm 500c}$
signal provides additional information about $\Delta P/P_{\rm DMO}$, leading
to 20-50% improvements when not using any cosmological information. When
additionally providing the $\Omega_{m}$ information, the improvements in
constraints on power or bi-spectrum bispectrum suppression are modest when
using the full pressure profile relative to the integrated quantities.

  * •

The pressure profiles and baryon fractions also carry information about
baryonic effects on the bi-spectrumbispectrum.

Our main findings from our analysis of constraints from the DESxACT shear-$y$
correlation analysis are

  * •

The measured electron pressure profile from the data analysis of tSZ and LSS
can then be used infer the matter power suppression. We infer competitive
constraints on this suppression using the measurements from shear-y
correlation from DES-ACT data.

  * •

We also show that the constraints would improve significantly in the future,
particularly using the halo catalog from DESI and tSZ map from S4.

With data from future galaxy and CMB surveys, we expect constraints on the tSZ
signal from halos across a wide mass and redshift range to improve
significantly. These improvements will come from both the galaxy side (e.g.
halos detected over larger areas of the sky down and out to higher redshifts)
and the CMB side (more sensitive tSZ maps over larger areas of the sky). Our
forecast for DESI and CMB Stage 4 in Fig. 10 suggests that very tight
constraints can be obtained on the impact of baryonic feedback on the matter
power spectrum. By combining these results with weak lensing constraints on
the small-scale matter distribution, we expect to be able to extract
significantly more cosmological information.

##  Appendix A Impact of limited volume of CAMELS simulations

In order to analyze the impact of varying box sizes and resolution on the
matter power suppression, we use the Illustris-TNG simulations as presented in
Springel:2018:MNRAS:. Particularly we use their boxes with side lengths of 210
Mpc/$h$, 75 Mpc/$h$ and 35 Mpc/$h$ (which they refer to as TNG-300, TNG-100
and TNG-50 as it corresponds to side length in the units of Mpc). We then make
the comparison to 25 Mpc/$h$ TNG boxes run from CAMELS. We use the CV set of
simulations and use them to infer the expected variance due to stochasticity
induced by changing initial conditions. Note that the hydrodynamical model is
identical between CAMELS CV runs and the bigger TNG boxes. In Fig. 11, we show
the power suppression for these boxes, including the runs at varying
resolution. We find that while changing box sizes gives relatively robust
values of power suppression, changing resolution can have non-negligible
impact. However, all the TNG boxes are consistent at 2-3$\sigma$ level
relative to the CAMELS boxes.

Figure 11: Comparison of the suppression of matter power in the CAMELS-TNG
simulation and simulations using same sub-grid prescription but larger box
sizes. We also show 1$\sigma$ and 2$\sigma$ uncertainty due to cosmic
variance. In the top panel we change the TNG box sizes, while preserving the
resolution, where as in the bottom panel we preserve the TNG box size while
changing the resolution.

##  Appendix B Robustness of results to different train sims

In Fig. 12, we test the impact of changing the simulations used to train the
random forest regressor. We then use these different trained models to infer
the constraints on the matter power suppression from the same stacked $\langle
Y_{500c}/Y^{\rm SS}\rangle$ as described in § 4. We see that our inferred
constraints remain consistent when changing the simulations.

Figure 12: In this figure we change the simulations used to train the RF when
inferring the power suppression from the data measurements.

##  Appendix C Test with lower halo masses

In Fig. 13, we show the constraints on the power suppression obtained by
analyzing the observables obtained from halos with lower masses, $1\times
10^{12}<M(M_{\odot}/h)<5\times 10^{12}$. We see that remarkably, even these
lower halo masses provide unbiased constraints on the matter power suppression
with robust inference especially at smaller scales. However, when compared to
the results descibed in § 3.2, we obtain less precise constraints. This is
expected as lower halos with lower masses are more susceptible to
environmental effects which induces a larger scatter in the relation between
their observables (such as $f_{b}$ or $Y_{500c}$) and their halo masses
governing feedback processes.

Figure 13: Same as Fig. 7 and Fig. 8, but obtained on lower halo masses,
$1\times 10^{12}<M(M_{\odot}/h)<5\times 10^{12}$. We find that having pressure
profile information results in unbiased constraints here as well, albeit with
a larger errorbars

##  Appendix D Test with other bispectrum configurations

In Fig. 14, we show the constraints obtained on the suppression of the
squeezed bi-spectrum bispectrum configurations. We fix the the angle between
the long sides of the triangle to correspond to $\mu=0.9$. We again find
robust inference on the bi-spectrum bispectrum suppression when using either
the integrated pressure profile or full radial pressure profile.

Figure 14: Same as Fig. 9, but for squeezed triangle configurations
($\mu=0.9$).

##  Appendix E Example of Emulation

We present an example of a constructed emulator according to § 3.1 for
$A_{AGN1}$ in Fig. 15. This shows how we aim to conservatively compute a
derivative without overestimating due to stochasticity.

Figure 15: The constructed emulator and resulting derivative for the
$A_{AGN1}$ parameter in the mass bin $10^{12}<M(M_{\odot}/h)<5\times 10^{12}$.