ROYAL ASTRONOMICAL SOCIETY

MNRAS 534, 655–682 (2024)

https://doi.org/10.1093/mnras/stae2100

Advance Access publication 2024 September 09

# Weak lensing combined with the kinetic Sunyaev–Zel’dovich effect: a study of baryonic feedback

L. Bigwood,¹,²* A. Amon,¹,²,³* A. Schneider,⁰,⁴ J. Salcido,⁰,⁵ I. G. McCarthy,⁰,⁵ C. Preston,⁰,¹

D. Sanchez,⁶ D. Sijacki,¹,² E. Schaan,⁷,⁸ S. Ferraro,⁹,¹⁰ N. Battaglia,¹¹ A. Chen,¹²,¹³ S. Dodelson,¹⁴,¹⁵

A. Roodman,¹⁶,¹⁷ A. Pieres,¹⁸,¹⁹ A. Ferté,¹⁷ A. Alarcon,²⁰,²¹ A. Drlica-Wagner,²²,²³,²⁴ A. Choi,²⁵

A. Navarro-Alsina,²⁶ A. Campos,¹⁴,¹⁵ A. J. Ross,²⁷ A. Carnero Rosell,¹⁸,²⁸ B. Yin,¹⁴ B. Yanny,²³

C. Sánchez,²⁹ C. Chang,²²,²⁴ C. Davis,¹⁶ C. Doux,²⁹,³⁰ D. Gruen,³¹ E. S. Rykoff,¹⁶,¹⁷ E. M. Huff,³²

E. Sheldon,³³ F. Tarsitano,³⁴ F. Andrade-Oliveira,¹² G. M. Bernstein,²⁹ G. Giannini,²⁴,³⁵ H. T. Diehl,²³

H. Huang,³⁶,³⁷ I. Harrison,³⁸ I. Sevilla-Noarbe,³⁹ I. Tutusaus,⁴⁰ J. Elvin-Poole,⁴¹ J. McCullough,¹⁶

J. Zuntz,⁴² J. Blazek,⁴³ J. DeRose,⁴⁴ J. Cordero,⁴⁵ J. Prat,²²,⁴⁶ J. Myles,³ K. Eckert,²⁹ K. Bechtol,⁴⁷

K. Herner,²³ L. F. Secco,²⁴ M. Gatti,²⁹ M. Raveri,⁴⁸ M. Carrasco Kind,⁴⁹,⁵⁰ M. R. Becker,²⁰ M. A. Troxel,⁵¹

M. Jarvis,²⁹ N. MacCrann,⁵² O. Friedrich,⁵³ O. Alves,¹² P.-F. Leget,¹⁶ R. Chen,⁵¹ R. P. Rollins,⁴⁵

R. H. Wechsler,¹⁶,¹⁷,⁵⁴ R. A. Gruendl,⁴⁹,⁵⁰ R. Cawthon,⁵⁵ S. Allam,²³ S. L. Bridle,⁴⁵ S. Pandey,²⁹

S. Everett,³² T. Shin,⁵⁶ W. G. Hartley,⁵⁷ X. Fang,³⁶,⁵⁸ Y. Zhang,⁵⁹ M. Aguena,¹⁸ J. Annis,²³ D. Bacon,⁶⁰

E. Bertin,⁶¹,⁶² S. Bocquet,³¹ D. Brooks,⁶³ J. Carretero,³⁵ F. J. Castander,²¹,⁶⁴ L. N. da Costa,¹⁸

M. E. S. Pereira,⁶⁵ J. De Vicente,³⁹ S. Desai,⁶⁶ P. Doel,⁶³ I. Ferrero,⁶⁷ B. Flaugher,²³ J. Frieman,²³,²⁴

J. García-Bellido,⁶⁸ E. Gaztanaga,²¹,⁶⁴,⁶⁰ G. Gutierrez,²³ S. R. Hinton,⁶⁹ D. L. Hollowood,⁷⁰

K. Honscheid,²⁷,⁷¹ D. Huterer,¹² D. J. James,⁷² K. Kuehn,⁷³,⁷⁴ O. Lahav,⁶³ S. Lee,³² J. L. Marshall,⁷⁵

J. Mena-Fernández,⁷⁶ R. Miquel,³⁵,⁷⁷ J. Muir,⁷⁸ M. Paterno,²³ A. A. Plazas Malagón,¹⁶,¹⁷ A. Porredon,⁷⁹

A. K. Romer,⁸⁰ S. Samuroff,⁴³ E. Sanchez,³⁹ D. Sanchez Cid,³⁹ M. Smith,⁸¹ M. Soares-Santos,¹²,⁸²

E. Suchyta,⁸³ M. E. C. Swanson,⁴⁹ G. Tarle,¹² C. To,²⁷ N. Weaverdyck,⁴⁴,⁵⁸ J. Weller,⁸⁴,⁸⁵ P. Wiseman⁸¹ and M. Yamamoto⁵¹

*Affiliations are listed at the end of the paper*

Accepted 2024 September 5. Received 2024 September 5; in original form 2024 April 15

## ABSTRACT

Extracting precise cosmology from weak lensing surveys requires modelling the non-linear matter power spectrum, which is suppressed at small scales due to baryonic feedback processes. However, hydrodynamical galaxy formation simulations make widely varying predictions for the amplitude and extent of this effect. We use measurements of Dark Energy Survey Year 3 weak lensing (WL) and Atacama Cosmology Telescope DR5 kinematic Sunyaev–Zel’dovich (kSZ) to jointly constrain cosmological and astrophysical baryonic feedback parameters using a flexible analytical model, ‘baryonification’. First, using WL only, we compare the $S_8$ constraints using baryonification to a simulation-calibrated halo model, a simulation-based emulator model, and the approach of discarding WL measurements on small angular scales. We find that model flexibility can shift the value of $S_8$ and degrade the uncertainty. The kSZ provides additional constraints on the astrophysical parameters, with the joint WL + kSZ analysis constraining $S_8 = 0.823_{-0.020}^{+0.019}$. We measure the suppression of the non-linear matter power spectrum using WL + kSZ and constrain a mean feedback scenario that is more extreme than the predictions from most hydrodynamical simulations. We constrain the baryon fractions and the gas mass fractions and find them to be generally lower than inferred from X-ray observations and simulation predictions. We conclude that the WL + kSZ measurements provide a new and complementary benchmark for building a coherent picture of the impact of gas around galaxies across observations.

**Key words:** gravitational lensing: weak – large-scale structure of Universe – cosmology: observations.

<footnote>* E-mails: lmb224@cam.ac.uk (LB); alexandra.amon@princeton.edu (AA)</footnote>

<footnote>© 2024 The Author(s).</footnote>

<footnote>Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.</footnote>
# 1 INTRODUCTION

The standard cosmological model, lambda-cold dark matter (ΛCDM), has been very successful when tested against observations of the cosmic microwave background (CMB; Bennett et al. 2003; Planck Collaboration VI 2020a), the lensing of the CMB at intermediate redshifts (Planck Collaboration VIII 2020b; Madhavacheril et al. 2024) and low-redshift observations of the expansion history as probed by baryon acoustic oscillations (Eisenstein et al. 2005). Measurements of weak galaxy lensing provide a strong test of ΛCDM at relatively small scales in the low-redshift Universe. In order to extract unbiased cosmological constraints from weak lensing (WL), accurate modelling of the non-linear matter distribution at $k > 0.1\ h\ Mpc^{-1}$ is crucial. This requires understanding both the non-linear dark matter evolution due to gravity to percent-level accuracy, as well as the impact of the baryons. In particular, a number of physical processes associated with baryons redistribute gas and impact the non-linear matter power spectrum by up to ~30 per cent (see e.g. Chisari et al. 2019; van Daalen, McCarthy & Schaye 2020). This ‘baryonic feedback’ encapsulates the energetic effect of active galactic nuclei (AGNs) heating gas and ejecting it to the outskirts of groups and clusters, as well as the likely subdominant effects of stellar winds, supernovae feedback, and gas cooling (van Daalen et al. 2011). At present, astrophysical model uncertainties, such as those due to baryonic feedback, have been shown to limit the precision of WL surveys (Amon et al. 2022; DES & KiDS Collaborations 2023). Therefore, to extract maximal information from WL data demands improved modelling of baryonic effects and their impact on the matter distribution. Moreover, it has been proposed that the ‘$S_8$ tension’ – the finding that WL constraints on the clustering amplitude parameter, $S_8 = \sigma_8(\Omega_m/0.3)^{0.5}$, are lower than predictions from the CMB – could be explained by a suppression of the non-linear matter power spectrum (Amon & Efstathiou 2022; Preston, Amon & Efstathiou 2023). This could either be caused by a more extreme baryonic feedback effects than hydrodynamical simulations predict or extensions to the standard model of cosmology. In order to isolate a departure from the standard cosmological model, baryonic effects must be better understood.

Powerful AGN feedback is believed to have the ability to eject baryons beyond the virial radius of galaxies, redistributing the gas to the outskirts of galaxy groups and clusters (McCarthy et al. 2011; Dubois et al. 2016; Henden et al. 2018; Springel et al. 2018). Indeed, studies of hydrodynamical simulations have demonstrated that AGN feedback alters the total matter distribution relative to dark-matter-only simulations, and that it causes a suppression of the power spectrum at scales $0.1 \lesssim k \lesssim 10\ h\ Mpc^{-1}$, whereas increased star formation can enhance power on the smallest scales (see Chisari et al. 2019, for a review). These simulations reproduce many of the observed properties of galaxies, including optical properties, galaxy group/cluster profiles, scaling relations, and Sunyaev–Zel’dovich counts. However, despite these successes, the scale, amplitude, and redshift dependence of the larger scale power suppression remain largely uncertain, with significant variation between simulations. These differences are direct outcomes of the ‘subgrid’ modelling of astrophysical processes, which take place below the resolution scale of the simulation.

Specifically, subgrid models are required to follow the formation, growth, and feedback of black holes, as well as gas cooling, metal enrichment, star formation, and associated stellar feedback. The AGN feedback may operate in either kinetic or thermal modes (generally associated with the radio and quasar modes, respectively), or alternate between the two depending on the black hole accretion

rate in a ‘two-mode’ feedback scenario (see Sijacki et al. 2007). The feedback model normally has a number of associated ill-constrained parameters encoding the feedback strength, such as the efficiency of thermal/kinetic coupling, the black hole accretion rate, and, in some models, the minimum heating temperature of gas cells before a feedback event occurs. While physical arguments can be used to narrow the plausible range of some of these parameters, this is normally not sufficiently constraining for precision cosmology purposes. Thus, it is often the case that the parameters are calibrated against key observables.

X-ray measurements of the hot gas fractions of groups and clusters within the virial radius are most widely used to benchmark the simulations, along with galaxy stellar mass function, star formation history, and galaxy sizes (McCarthy et al. 2017; Henden et al. 2018; Nelson et al. 2019; Kugel et al. 2023; Schaye et al. 2023). Even with identical subgrid physics, the box size and resolution of the simulation can also have a non-negligible impact on the matter distribution (van Daalen et al. 2011; Pandey et al. 2023). Indeed, there is a large parameter space of feedback prescriptions, modelling choices and simulation properties that result in significant variation in the suppression of the matter power spectrum (Henden et al. 2018; Davé et al. 2019; Hernández-Aguayo et al. 2023; Salcido et al. 2023; Schaye et al. 2023).

Recent WL analyses have devised various approaches to mitigate the impact of baryonic feedback on cosmological constraints. The DES Y3 (Dark Energy Survey Year 3) cosmic shear analysis opted to discard measurements from the analysis on angular scales that are impacted by baryonic effects from the analysis (Krause et al. 2021; Amon et al. 2023; Secco, Samuroff et al. 2022). Alternatively, baryon feedback has been modelled using a halo model approach (Asgari et al. 2021; Li et al. 2023), and using a halo model that is calibrated to a hydrodynamical simulations (Mead et al. 2021; DES & KiDS Collaborations 2023). More recently, Salcido et al. (2023) have developed an emulator trained using a suite of hydrodynamical simulations with varied feedback efficiencies.

Instead of relying on the hydrodynamical simulations, the baryonification (bfc) model is another approach which shifts particle outputs in gravity-only simulations to attain modified halo profiles, modelling the re-arrangement of baryon material caused by feedback effects (Schneider & Teyssier 2015). This approach has been used in Schneider et al. (2022), Chen et al. (2023), and Aricò et al. (2023). Other approaches include a principal component analysis (Huang et al. 2019).

An alternative approach to using models that are informed by hydrodynamical simulations (and therefore indirectly benchmarked against X-ray data) is to jointly analyse WL data with observations of the gas content in and around galaxy groups and clusters. This has been done to improve cosmological constraints by reducing the model space of the nuisance parameters and exploiting the complementary dependence of the cross-correlated data sets on cosmological parameters (Tröster et al. 2020). The suppression of the matter power spectrum has also been constrained using this approach (Schneider et al. 2022).

A highly complementary observable to X-ray is the kinetic Sunyaev–Zel’dovich (kSZ) effect, caused by the Thomson scattering of the CMB photons by free electrons moving with bulk motion in groups and clusters of galaxies relative to the CMB rest frame (Sunyaev & Zeldovich 1972, 1980). This causes a shift in the CMB temperature while preserving the blackbody spectrum. If the bulk line-of-sight velocity is known, the kSZ effect directly measures the free electron number density, independent of temperature (Hand et al. 2012; Planck Collaboration XXXVII 2016; Soergel et al. 2016;
Schaan et al. 2021). The kSZ effect is well suited to probe low density and low temperature environments like the outskirts of galaxies and clusters, whereas X-ray measurements are more sensitive to the inner regions (Amodeo et al. 2021).

The goals of this work are threefold:

I. We test the performance of four baryon feedback mitigation strategies for analysing mock and DES Y3 WL data: the DES Y3 scale cut approach (Krause et al. 2021), a halo model approach (Mead et al. 2021) calibrated to a hydrodynamical simulation, an emulator built using a suite of hydrodynamical simulations (Salcido et al. 2023) and the bfc model (Schneider & Teyssier 2015; Schneider et al. 2019).

II. As an alternative to simulation-driven models, we use the most flexible model, bfc, to jointly analyse the lensing data with the Atacama Cosmology Telescope (ACT) kSZ measurements for improved constraints on the baryonic feedback parameters and therefore the cosmological parameters.

III. We constrain astrophysical observables for the first time in a joint WL and kSZ analysis, providing a new avenue to benchmark the hydrodynamical simulations.

The paper is structured as follows. Section 2 describes the DES Y3 cosmic shear and ACT DR5 kSZ data sets used in this analysis. Section 3 outlines the modelling of the cosmic WL and kSZ measurements, including the four baryon models.

In Section 4, we summarize the findings of a mock analysis. We compare the cosmological constraints when analysing the shear data with different baryon mitigating models and model complexities in Section 5. In Section 6, we present our constraints on both cosmological and baryonic parameters obtained from a joint WL + kSZ analysis. Finally, in Section 7, we consider our constraints on the observables that simulations benchmark against. We summarize key findings and conclude in Section 8.

## 2 DATA

### 2.1 Dark Energy Survey cosmic shear

The DES has completed six years of photometric observations in the grizY bands, using the 4-m Blanco telescope located at the Cerro Tololo Inter-American Observatory. The survey spans ~5000 deg² in the Southern Hemisphere.

For this analysis, we use data taken during the survey's first three years of operation (DES Y3), between 2013 and 2016 (Sevilla-Noarbe et al. 2021). The DES Y3 footprint covers 4143 deg² with a number density of 5.59 galaxies arcmin⁻² to a depth of i ~ 23.5. The METACALIBRATION WL catalogue has over 100 million galaxies that have passed a raft of validation tests (Gatti et al. 2022). The sample has been divided into four redshift bins and the calibrated redshift distributions and associated uncertainty are defined in Myles et al. (2021). Remaining biases in the shape measurement and redshift distributions, primarily due to blending, are calibrated using image simulations, and the associated corrections for each redshift bin are reported in MacCrann et al. (2022). The DES Y3 cosmic shear tomographic two-point correlation functions, ξ_±, are measured in 20 angular logarithmic bins spanning 2.5–250.0 arcmin (Amon et al. 2022; Secco et al. 2022).

### 2.2 Atacama Cosmology Telescope kinetic Sunyaev–Zel'dovich

The ACT is a 6-m millimeter waveband telescope, observing the CMB. Since first light in 2007, it has had three generations of

receivers, the most recent of which is the polarization-sensitive Advanced ACTPol (AdvACT), which extended the frequency coverage to five bands spanning 28–230 GHz. The fifth data release (herinafter DR5) co-adds maps collected from 2008 to 2018 covering approximately 18 000 deg² and utilizes data from all three generations of receiver (Naess et al. 2020).

This work uses the kSZ measurements presented in Schaan et al. (2021). These are stacked measurements of the ACT DR5 and Planck CMB temperature maps at 98 GHz (hereinafter called f90 for consistency with Schaan et al. 2021) and 150 GHz (f150) with the reconstructed velocities of the spectroscopic Baryon Oscillation Spectroscopic Survey (BOSS) CMASS galaxy catalogue. The kSZ signal is detected at a 6.5σ significance. The galaxy sample spans the redshift range 0.4 < z < 0.7 with a median redshift of z = 0.55. It corresponds to a selection of relatively massive galaxies with a mean stellar mass of log₁₀(M<sub>star</sub>/M☉) ≈ 11.3 and an assumed mean halo mass of ~10¹³ M☉ (Amodeo et al. 2021; Schaan et al. 2021), though the latter is quite uncertain. Given its importance to the modelling of the stacked kSZ sample, we will discuss the mean halo mass of this sample later in the work.

### 2.3 Hydrodynamical simulations

Throughout the paper, we compare our constraints on the matter power spectrum to predictions from a range of hydrodynamical simulations: the FLAMINGO (1 Gpc)³ box with baryonic particle masses of 10⁹ M☉ (Schaye et al. 2023, solid line); BAHAMAS (400 Mpc h⁻¹)³ box with Θ<sub>AGN</sub> = 7.8 and 1024³ dark matter and baryonic particles (McCarthy et al. 2017, dashed-dotted line); SIMBA (100 Mpc h⁻¹)³ box with 1024³ gas elements (Davé et al. 2019, dotted line); MillenniumTNG 740 Mpc³ with baryonic mass resolution of 3.1 × 10⁷ M☉ (Pakmor et al. 2023, dashed line); and FABLE (100 Mpc h⁻¹)³ box with 1280³ dark matter particles and 1280³ baryonic resolution elements (double-dotted-dashed line, Henden et al. 2018, Bigwood et al. in preparation). These simulations not only span a range of box sizes and resolutions, but also feedback implementations, hydrodynamical schemes, and calibration strategies. Fig. 1 shows the suppression of the non-linear matter power spectrum due to baryonic effects, P(k)/P<sub>DMonly</sub>, as predicted by each simulation, demonstrating the spread in the predicted amplitude and scale dependence. Several additional hydro simulations are also used to calibrate the baryonic feedback models tested throughout this work. We therefore also plot the prediction from the OWLS-AGN (100 Mpc h⁻¹)³ box with 512³ dark matter and baryonic particles, respectively (Schaye et al. 2010; van Daalen et al. 2011, dark green solid line), the span in the suppression predicted by BAHAMAS when Θ<sub>AGN</sub> is varied within the range 7.6–8.0 (see Section 3.2.2), and range modification measured by the ANTILLES suite of 400 simulations (Salcido et al. 2023, light green region), each with box size (100 Mpc h⁻¹)³ and 256³ dark matter and baryonic particles, respectively. ANTILLES spans significantly more extreme feedback scenarios than the other simulations considered, such that baryonic effects impact the power spectrum with greater amplitude and at smaller k scales.

## 3 MODELLING

### 3.1 Cosmic shear signal

The shear two-point correlation functions, ξ_±(θ), for a given angular separation, θ, computed for redshift bins i, j, can be related to the 3D non-linear matter power spectrum. First, it can be expressed as a decomposition of the 2D convergence power spectrum C<sub>k</sub><sup>i,j</sup>(ℓ) at an
<figure><img src="image_1.png" /><figcaption>**Figure 1.** The suppression of the matter power spectrum due to baryonic effects predicted by hydrodynamical simulations. We show OWLS-AGN (Schaye et al. 2010; van Daalen et al. 2011, dark green solid line), which informed the scale cuts baryonic feedback mitigation approach. We also plot the matter power spectrum suppression measured in the 400 simulations of the ANTILLES suite, which calibrated the SP(k) emulator (Salcido et al. 2023, light green region). The range in the suppression predicted by the BAHAMAS suite spanning $Θ_{AGN} = 7.6-8.0$ from which HM20 was calibrated (McCarthy et al. 2017, dark blue region). We further plot FLAMINGO (Schaye et al. 2023, black solid line); BAHAMAS (McCarthy et al. 2017, black dashed-dotted line); SIMBA (Davé et al. 2019, black dotted line); MillenniumTNG (Pakmor et al. 2023, black dashed line); and FABLE (double-dotted-dashed line, Henden et al. 2018; Bigwood et al. in preparation). Finally, we show in red the three baryonic feedback scenarios we test in our mock analysis (Section 4.2): BAHAMAS $T_{heat} = 8.0$ (dotted red line), cosmo-OWLS $T_{heat} = 8.5$ (dashed red line), and $A_{mod} = 0.858$ (solid red line, Preston et al. 2023).</figcaption></figure>

angular wave number, $ℓ$, as

$$ \xi_{\pm}(\theta) = \sum_{\ell} \frac{2\ell + 1}{4\pi} G_{\ell}^{\pm}(\cos\theta) [C_{\kappa,EE}^{ij}(\ell) \pm C_{\kappa,BB}^{ij}(\ell)]. \quad (1) $$

We note that weak gravitational lensing does not produce B modes. However, we show the more general expression here, as a B-mode contribution from intrinsic alignments (IA) is possible. The functions $G_{\ell}^{\pm}$ are computed from the Legendre polynomials following Stebbins (1996).

Under the Limber and flat-sky approximations (Limber 1953; LoVerde & Afshordi 2008), we can relate $C_{\kappa}^{ij}(\ell)$ to the 3D non-linear matter power spectrum, P, via

$$ C_{\kappa}^{ij}(\ell) = \int_{0}^{\chi_{H}} d\chi \frac{W_{i}(\chi)W_{j}(\chi)}{\chi^{2}} P\left(k = \frac{\ell + 0.5}{\chi(z)}, z\right). \quad (2) $$

Here, we assume a spatially flat Universe, $\chi$ is the comoving angular diameter distance and $\chi_H$ is the horizon distance. $W_i(\chi)$ are the lensing efficiency kernels, given by:

$$ W_i(\chi) = \frac{3H_0^2 \Omega_m}{2c^2} \frac{\chi}{a(\chi)} \int_{\chi}^{\chi_H} d\chi' n_i(\chi') \frac{\chi' - \chi}{\chi'}, \quad (3) $$

where $a(\chi)$ is the scale factor at comoving distance $\chi$, and $n_i(\chi)d\chi$ is the effective number density of galaxies in $d\chi$, normalized to unity.

In this analysis, the linear matter power spectrum is calculated using CAMB (Lewis, Challinor & Lasenby 2000), and the non-linear

correction is determined using HMCODE2020 (Mead et al. 2021). We refer the reader to DES & KiDS Collaborations (2023) for a detailed comparison of cosmic shear analysis choices, which we use as a guide to formulate the baseline model in this work. We assume three neutrino species with two massless states and one massive state with a mass fixed at the minimum mass allowed by oscillation experiments, $m_v = 0.06\text{eV}$ (Patrignani 2016).

The IA of galaxies with their local environment also contributes to the shear correlation function and must be modelled. We chose to do this using the non-linear linear-alignment model (NLA), which describes the linear tidal alignment of galaxies with the density field (Hirata & Seljak 2004), with a non-linear correction to the linear matter power spectrum (Bridle & King 2007). This approach requires two additional free nuisance parameters: $A_{IA}$, modulating the amplitude of the IA model (see equations 3–5 in Bridle & King 2007 for the NLA IA power spectra, $C_{GI}$ and $C_{II}$) and a redshift-dependence parameter, using a power law with $[(1+z)/(1+0.62)]^{\eta_{IA}}$. The choice of modelling IA with NLA, and the dependence of our cosmological results on this model choice, are further discussed in Section 6.1 and Appendix C.

We model the uncertainty in mean redshift and the shear calibration for each redshift bin $i$ as the free parameters $\Delta z^i$ and $m^i$, respectively, following Amon et al. (2022) and Secco et al. (2022), and preserve values of the uncertainty determined by Myles et al. (2021) and MacCrann et al. (2022). In the cases where the small-scale measurements are analysed, we refer the reader to Chen et al. (2023) for validation that higher order cosmic shear modelling corrections remain subdominant.

## 3.2 Modelling baryonic feedback for cosmic shear

Strategies have been devised to account for baryonic effects when analysing WL measurements in order to extract unbiased cosmological information. In this section, we outline the four approaches we investigate in this work: (1) restrict angular scales, (2) a halo model approach, (3) a hydrodynamical simulation-based emulator, and (4) an analytical N-body simulation model. In Table 1, we summarize the free parameters of each baryon feedback model and prior choices.

### 3.2.1 Restricting angular scales

The DES Y3 cosmic shear analysis mitigates the impact of baryonic effects by discarding the measurements at small angular scales (Krause et al. 2021; Amon et al. 2023; Secco et al. 2022) and analysing the data assuming a dark-matter-only model. In this work, we adopt their ‘ΛCDM Optimized’ scale cuts, which were designed to minimize the bias on ΛCDM cosmological parameters due to unmodelled baryonic effects to be less than 0.3σ$_{2D}$ in the $\Omega_m - S_8$ parameter space, for a joint lensing and clustering analysis. For a ΛCDM analysis of cosmic shear alone, this corresponds to up to a 0.14σ$_{2D}$ potential bias. Note that to determine the angular scales to be used, the OWLS-AGN hydrodynamical simulation (Schaye et al. 2010; van Daalen et al. 2011) was chosen as a representative feedback scenario (shown as the green line, Fig. 1). Synthetic cosmic shear data vectors were contaminated according to

$$ P_b(k, z) = \frac{P_{\text{hydro}}(k, z)}{P_{\text{DMO}}(k, z)} P(k, z), \quad (4) $$

where $P_{\text{hydro}}(k, z)$ and $P_{\text{DMO}}(k, z)$ are the full hydrodynamical and dark-matter-only matter power spectra from the OWLS-AGN suite and $P(k, z)$ is the non-linear matter power spectrum.
<table><caption>**Table 1.** Qualitative descriptions of the parameters associated with the HM20, SP(*k*), and BCEmu methods, which each model the impact of baryonic feedback on the non-linear power spectrum. We also show the priors utilized on each parameter in this analysis, and the conservative prior alternative if applicable. U[] brackets indicate flat uniform priors within the range shown. N() brakets indicate Gaussian priors described by their mean and 1σ width. The BCEmu1 model has log<sub>10</sub> M<sub>c</sub> as the only free parameter, BCEmu3 varies log<sub>10</sub> M<sub>c</sub>, θ<sub>ej</sub> and η<sub>cga</sub>, and all seven parameters are free in BCEmu7. In the case of using the reduced model complexity of BCEmu1 or BCEmu3 models, we also show the values that the BCEmu parameters are fixed to. Note that the M<sub>200</sub> mass parameter is only included in the BCEmu joint WL + kSZ analysis and the choice of the fixed value is adopted from Schaan et al. (2021).</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th><th>Wide prior</th><th>Fixed value</th></tr></thead><tbody><tr><td>Halo model: HM20<br>Θ<sub>AGN</sub> log<sub>10</sub>(T<sub>AGN</sub>/K)</td><td>Subgrid heating parameter calibrated to the BAHAMAS simulations designed to modulate the amplitude and shape of the 'one-halo' term in the halo model</td><td>U[7.6, 8.0]</td><td>U[7.3, 9.0]</td><td>—</td></tr><tr><td>Simulation-based: SP(k)<br>α</td><td>Power-law normalization for the f<sub>b</sub>-M<sub>halo</sub> relation (equation 7)</td><td>N(4.16, 0.07)</td><td>U[2.85, 4.50]</td><td>—</td></tr><tr><td>β</td><td>Power-law slope for the f<sub>b</sub>-M<sub>halo</sub> relation (equation 7)</td><td>N(1.20, 0.05)</td><td>U[0.95, 1.60]</td><td>—</td></tr><tr><td>γ</td><td>Redshift dependence of the power-law normalization of the f<sub>b</sub>-M<sub>halo</sub> relation (equation 7)</td><td>N(0.39, 0.09)</td><td>U[0.12, 0.85]</td><td>—</td></tr><tr><td>Baryonification: BCEmu<br>log<sub>10</sub> M<sub>c</sub></td><td>The characteristic mass scale at which the slope of the gas profile becomes shallower than −3 (equations 9 and 10)</td><td>U[11.0, 15.0]</td><td>—</td><td>—</td></tr><tr><td>θ<sub>ej</sub></td><td>Specifies the maximum radius of gas ejection relative to the virial radius</td><td>U[2.0, 8.0]</td><td>—</td><td>3.5</td></tr><tr><td>η<sub>δ</sub></td><td>Related to the stellar fraction of the central galaxy: η<sub>δ</sub> = η<sub>cga</sub> − η (equation 11)</td><td>U[0.05, 0.4]</td><td>—</td><td>0.20</td></tr><tr><td>μ</td><td>Defines how fast the slope of the gas profile becomes shallower towards small halo masses (equations 9 and 10)</td><td>U[0.0, 2.0]</td><td>—</td><td>1.0</td></tr><tr><td>γ</td><td>Exponent in gas profile parametrization (equation 9)</td><td>U[1.0, 4.0]</td><td>—</td><td>2.5</td></tr><tr><td>δ</td><td>Exponent in the gas profile parametrization (equation 9)</td><td>U[3.0, 11.0]</td><td>—</td><td>7.0</td></tr><tr><td>η</td><td>Specifies the total stellar fraction within a halo (equation 11)</td><td>U[0.05, 0.4]</td><td>—</td><td>0.20</td></tr><tr><td>M<sub>h,200</sub> (M<sub>☉</sub>)</td><td>Halo mass of the kSZ sample, used in the joint analysis only</td><td>U[5 × 10<sup>12</sup>, 7 × 10<sup>13</sup>]</td><td>—</td><td>3 × 10<sup>13</sup></td></tr></tbody></table>

A benefit of this approach is that it is agnostic to the exact shape and physics of the matter power spectrum suppression, once the feedback in the Universe is lower in amplitude and scale extent than the simulation chosen (in this case, OWLS-AGN). However, this approach misses the opportunity to extract high-signal-to-noise information about the underlying cosmological model and the astrophysical effects.

### 3.2.2 Halo-model approach: HM20

HMCODE2020, hereinafter HM20, models the non-linear power spectrum and includes a free parameter to modulate the amount of baryonic feedback, Θ<sub>AGN</sub> = log<sub>10</sub>(T<sub>AGN</sub>/K) (Mead et al. 2021). This parameter scales the halo concentration and the stellar and gas content, leading to a modification in the overall amplitude and shape of the 'one-halo' term in the halo model. The model is calibrated to fit the power spectrum 'response' (the matter-matter power spectrum divided by the same measurement in an equivalent dark-matter-only box) of the BAHAMAS hydrodynamical simulations (McCarthy et al. 2017; van Daalen et al. 2020) in the range Θ<sub>AGN</sub> = [7.6 - 8.0] (blue shaded region, Fig. 1). We note that T<sub>AGN</sub> is related to ΔT<sub>heat</sub>, which is the BAHAMAS subgrid heating parameter, where an AGN feedback event will only occur after the black hole has stored sufficient energy to heat a fixed number of gas particles by ΔT<sub>heat</sub>.

We define two prior ranges for this case. The first spans the range of Θ<sub>AGN</sub> values that bracket the BAHAMAS Θ<sub>AGN</sub> = 7.6 - 8.0 simulations that the model was calibrated against. The 'wide prior' chosen here to be Θ<sub>AGN</sub> = 7.3 - 9.0 extends beyond the calibration range to span more extreme scenarios and allow for a dark-matter-only case.

This approach has been shown to be accurate at the level of < 2.5 per cent to k < 10 h Mpc<sup>-1</sup> (Mead et al. 2021) when fitting

simulated power spectra at z < 2 with a range of cosmologies thus allowing all measured angular scales of the DES Y3 data set to be utilized (in this case, to 2.5 arcmin). A downside of this model is that it relies on the accuracy of the specific feedback implementation and predicted power suppression of a particular simulation and may not be flexible enough to capture the true scenario.

### 3.2.3 Hydrodynamical simulation emulator: SP(k)

SP(*k*) is a flexible empirical model trained on the ANTILLES suite of 400 cosmological hydrodynamical simulations (plotted as the green lines in Fig. 1, Salcido et al. 2023). The model predicts the power spectrum suppression given the baryon fraction-halo mass relation of galaxy groups and clusters as input, building upon previous insights from van Daalen et al. (2020). The ANTILLES suite spans a range of feedback scenarios, allowing the emulator to achieve a ~2 per cent level accuracy to describe baryonic effects at scales of up to k < 10 h Mpc<sup>-1</sup> and redshifts up to z = 3.

In particular, SP(*k*) casts the suppression in terms of the baryon fraction at the optimal mass, M&#770;<sub>k</sub>, defined as the halo mass that maximizes the strength of the correlation between the suppression of the total matter power spectrum and the total baryon fraction of haloes of different mass. It uses an exponential plateau function to model for the fractional impact of baryons as,

$$P_{\text{hydro}}(k)/P_{\text{DM}}(k) = \lambda(k, z) - [\lambda(k, z) - \mu(k, z)] \times \exp[-\nu(k, z)\tilde{f}_{\text{b}}], \quad (5)$$

where $\tilde{f}_{\text{b}}$ is the baryon fraction at the optimal halo mass normalized by the universal baryon fraction, that is:

$$\tilde{f}_{\text{b}} = f_{\text{b}}(\hat{M}_{k,\text{SO}}(k, z)) / (\Omega_{\text{b}} / \Omega_{\text{m}}), \quad (6)$$
and $\hat{M}_{k,SO}$, $\lambda(k, z)$, $\mu(k, z)$, and $\nu(k, z)$ are functions with best-fitting parameters given in Salcido et al. (2023).

For the mass range that can be relatively well probed in current X-ray and Sunyaev-Zel'dovich effect observations ($10^{13} \lesssim M_{200}$ [M$_{\odot}$] $\lesssim 10^{15}$), the total baryon fraction of haloes can be roughly approximated by a power law with constant slope (e.g. Mulroy et al. 2019; Akino et al. 2022). Salcido et al. (2023) find that a modified version of the functional form,<fnref n="1" /> presented in Akino et al. (2022) provides a reasonable agreement with simulations up to redshift $z=1$,

$$f_b/(\Omega_b/\Omega_m) = \left(\frac{e^\alpha}{100}\right) \left(\frac{M_{500c}}{10^{14} M_{\odot}}\right)^{\beta-1} \left(\frac{E(z)}{E(0.3)}\right)^\gamma, \quad (7)$$

where $\alpha$ sets the power-law normalization, $\beta$ sets power-law slope, $\gamma$ provides the redshift dependence, $E(z)$ is the dimensionless Hubble parameter, and $f_b$ is the baryon fraction measured within $R_{500}$. We use this function to facilitate the marginalization over the uncertainties in the observed baryon fractions, introducing $\alpha$, $\beta$, and $\gamma$ as free parameters in our WL analysis.

SP(k) effectively depends only on a single physically meaningful parameter, that is, the baryon fraction. The benefit of this approach is that observational constraints on the baryon fraction could be used to inform the priors used in cosmological analysis. For our study, we use two different sets of priors to marginalize over the uncertainties in the observed baryon fractions: wide (conservative) priors consistent with the range of feedback models probed by the ANTILLES simulations used to calibrate the SP(k) model (Salcido et al. 2023), and ‘observational’ priors that encompass current observational constraints on the baryon fraction-halo mass relation from Akino et al. (2022). Table 1 reports the priors on $\alpha$, $\beta$, and $\gamma$ for the two choices.

As with the scale cuts and HM20 model, we caveat that the SP(k) method is limited by the specific feedback implementation used in the hydrodynamical simulation it was trained on, in this case the ANTILLES suite. While ANTILLES is currently the largest and widest suite of hydro simulations in terms of feedback variations, it still may not cover all possible baryonic responses (see e.g. Appendix B4).

### 3.2.4 Analytical N-body simulation model: Baryonification

Bfc is a method for including the effects of baryonic feedback in dark-matter-only N-body simulations based on perturbative shifts of particles that mimic the effects of feedback at cosmological scales (Schneider & Teyssier 2015; Schneider et al. 2019; Aricò et al. 2020). The particles are shifted in order to obtain modified halo profiles that include the presence of gas and stars which are shaped by feedback effects. We provide a summary of the method including some important aspects of the parametrization and refer the reader to Schneider et al. (2019) and Giri & Schneider (2021) for a more detailed explanation.

In practice, the bfc method relies on a modification of profiles via spherically symmetrical particle shifts around the centres of N-body haloes with Navarro-Frenk-White (NFW)-like profiles, $\rho_{nfw}$, following

$$\rho_{nfw}(r) \rightarrow \rho_{bfc}(r) = \rho_{clm}(r) + \rho_{gas}(r) + \rho_{cga}(r). \quad (8)$$

The final baryonified profiles, $\rho_{bfc}$, consist of three components: the collisionless matter, gas, and central galaxy. The collisionless

matter ($\rho_{clm}$) profile is dominated by dark matter but also contains satellite galaxies and halo stars. Its shape is modified with respect to the original NFW shape via adiabatic contraction and expansion (Teyssier et al. 2011). The central galaxy profile, $\rho_{cga}$, is parametrized as a power law with an exponential cut-off. This component affects only the innermost part of the halo, rather than cosmological scales.

At cosmological scales, baryonic effects are primarily caused by feedback-induced changes of the gas distribution around haloes, described by the gas profile, $\rho_{gas}$. Motivated by X-ray observations (Eckert et al. 2016; Schneider et al. 2019), these effects are parametrized with five model parameters as

$$\rho_{gas}(r) \propto \frac{\Omega_b / \Omega_m - f_{star}(M)}{\left[1 + \left(\frac{r}{r_{core}}\right)^{\beta(M)}\right] \left[1 + \left(\frac{r}{r_{ej}}\right)^{\gamma}\right]^{\frac{\delta - \beta(M)}{\gamma}}}, \quad (9)$$

which consists of a cored power-law profile with a truncation at the ejection radius, $r_{ej} = \theta_{ej} R_{200}$, where $\theta_{ej}$ is a dimensionless-free model parameter.<fnref n="2" /> The shape of the truncation beyond the ejection radius is controlled by the $\gamma$ and $\delta$ parameters, where the former defines the abruptness and the latter defines the slope of the function. The core of the profile is fixed to $r_{core} = \theta_{core} R_{200}$ with $\theta_{core} = 0.1$. The slope of the power law $\beta(M)$ is a function that varies with halo mass,<fnref n="3" /> and is parametrized as

$$\beta(M) = \frac{3(M/M_c)^\mu}{1 + (M/M_c)^\mu}. \quad (10)$$

The halo mass dependence of $\beta$ accounts for the fact that AGN feedback is more efficient in removing gas around galaxy groups while large clusters tend to keep most of their gas inside the virial radius. The free model parameters $M_c$ and $\mu$ thereby define the scale and the abruptness of the transition when the slope of the profile goes from 3 to 0 for decreasing halo masses.

The total fraction of stars, $f_{star}$, and the fraction of stars that belong to the central galaxy, $f_{cga}$, indirectly affect the available gas that can be pushed out by feedback processes. They are parametrized as

$$f_i(M) = 0.055 \left( \frac{M}{M_s} \right)^{\eta_i} \quad (11)$$

with $i = \{\text{star, cga}\}$ and where $M_s = 2.5 \times 10^{11} M_{\odot} h^{-1}$. The power law is constructed to match the Moster relation (Moster, Naab & White 2013), which is reasonably well understood for galaxy groups and clusters. Note that following Giri & Schneider (2021), we redefine the parameters as $\eta = \eta_{star}$ and $\delta_i = \eta_{cga} - \eta_{star}$. These two additional bfc parameters, together with the five gas parameters, are summarized in Table 1.

An efficient cosmological analysis will marginalize over a minimum number of baryonic feedback parameters. How many parameters are required is an open research question that ultimately depends on the unknown baryonic feedback realized in nature. For now, we adhere to the requirement that a given parametrization needs to be able to fit the matter power spectrum suppression predicted by a range of hydrodynamical simulations. Following Giri & Schneider (2021), in this work we consider the models BCEmu7, BCEmu3, and BCEmu1, referring to the number of model parameters varied in the analysis. While in BCEmu7 all parameters introduced above are kept free, BCEmu3 only allows $\log_{10} M_c$, $\theta_{ej}$, and $\eta_\delta$ to vary

<footnote n="2">²Note that the gas profile is normalized so that an integral over $r^2/2\pi^2 \times \rho$ gives the total halo mass.</footnote>

<footnote n="3">³Here, we define the halo mass as the mass enclosed within a radius, centred on the group or cluster, within which the mean density is 200 times the critical density of the Universe. The notation $M_{200}$ is also used throughout the paper.</footnote>

<footnote n="1">¹See https://github.com/jemme07/pyspk.</footnote>
and BCEmu1 only uses $log_{10} M_c$ as a free parameter. The fixed parameters in the BCEmu3 and BCEmu1 models are listed in Table 1 and have been selected so that they provide the best fit to a variety of hydrodynamical simulations (see Giri & Schneider 2021, for more information about the method). Both BCEmu7 and BCEmu3 are able to reproduce the baryonic suppression or the power spectrum predicted by a variety of hydrodynamical simulations to better than one percent. The BCEmu1 model, on the other hand, shows deviations of order five per cent, hinting towards the possibility that one parameter is generally insufficient to describe the variety of existing results from simulations (Giri & Schneider 2021).

A key feature of the bfc model is that it is based on physically motivated profiles around halo centres. The model is not restricted to the power spectrum but can also be used to obtain the 3D baryonified density field and therefore many corresponding summary statistics. In the following we will take advantage of the connection between power spectrum and halo profiles to obtain simultaneous predictions for the cosmic shear and the stacked kSZ signal. Compared to other approaches, the bfc model BCEmu7 is not calibrated to specific hydrodynamical simulations. It rather depends on the parameters of empirical density profiles that are broadly motivated by observations. The model therefore provides a independent check with very different modelling choices and systematics compared to the subgrid modelling in hydrodynamical simulations. We note that despite the model's valuable flexibility, it can result in potentially non-physical scenarios. Other limitations include the fact that the model parameters are currently assumed to have no redshift dependence, as well as the fact that the gas profiles do not separate a hot and cold gas component.

## 3.3 Kinetic Sunyaev-Zel'dovich signal

The kSZ measurements can be used to constrain the gas density. This effect arises from the bulk motion of the ionized gas in and around galaxies, galaxy groups, and clusters, which imparts a Doppler shift on CMB photons. It preserves the blackbody frequency spectrum of the CMB and instead the thermodynamic temperature fluctuates as

$$ \frac{\Delta T_{\text{kSZ}}}{T_{\text{CMB}}} = \frac{\sigma_{\text{T}}}{c} \int_{\text{los}} e^{-\tau} n_{\text{e}} v_{\text{p}} dl, \qquad (12) $$

where $T_{\text{CMB}}$ is the present-day temperature of the CMB, $\sigma_{\text{T}}$ is the Thomson cross-section, c is the speed of light, $n_e$ is the free-electron physical number density, $v_p$ is the peculiar velocity, and $\tau$ is the optical depth due to Thomson scattering along the line-of-sight distance, dl. Following Schaan et al. (2021) and Amodeo et al. (2021), for the redshift range of the kSZ measurements used in this work, the mean optical depth is observed to be below the percent level (Planck Collaboration XXXVII 2016) and the CMASS galaxy groups are optically thin, therefore we can assume that the integral $e^{-\tau} \approx 1$. Furthermore, as the measurements are stacked, the velocity field is independently estimated from the large-scale distribution of galaxies via a reconstruction method, thereby eliminating the dependence on the velocity. We therefore follow Schaan et al. (2021) and Amodeo et al. (2021) in simplifying equation (12), so that the resulting shift in the CMB temperature can be approximated as;

$$ \frac{\Delta T_{\text{kSZ}}}{T_{\text{CMB}}} = \tau_{\text{gal}}(\theta) \frac{v_r}{c}, \qquad (13) $$

where $\tau_{\text{gal}}$ refers to the contribution of the optical depth to Thomson scattering of the galaxy group considered and $v_r$ is the root mean square (RMS) of the peculiar velocities projected along the line

of sight. For the median redshift of the CMASS sample in the linear approximation, $z = 0.55$, the RMS of the peculiar velocities projected along the line of sight is $v_r = 1.06 \times 10^{-3} c$ (Schaan et al. 2021). The uncertainty on the velocity reconstruction is estimated to be less than a few percent, which we ignore given the statistical precision (Schaan et al. 2021; Ried Guachalla et al. 2024; Hadzhiyska et al. 2024).

To model the ACT kSZ measurements it is necessary to convolve equation (13) with the beam profiles utilized at the f90 and f150 frequencies. We follow Schaan et al. (2021) and approximate the beam using a Gaussian with FWHM = 2.1 arcmin for the former band and FWHM = 1.3 arcmin for the latter. Furthermore, to minimize noise due to degree-scale CMB fluctuations, compensated aperture photometry filters were also used on the observations. We therefore apply the same filter function that was used in the analysis of the data. The smoothing function is +1 between $\theta < \theta_d$, -1 between $\theta_d < \theta < \sqrt{2}\theta_d$ and 0 otherwise, where $\theta_d$ is the aperture radius centred around each galaxy (Schaan et al. 2021, equation 9).

To calculate $\tau_{\text{gal}}$ measured within a disk of radius $\theta$ centred on the group or cluster, we assume spherical symmetry and integrate the electron number density, $n_e$, over the line of sight as

$$ \tau_{\text{gal}}(\theta) = 2\sigma_{\text{T}} \int_{0}^{\infty} n_{\text{e}}(\sqrt{l^2 + d_{\text{A}}(z)^2\theta^2}) dl, \qquad (14) $$

where $d_{\text{A}}$ is the angular diameter distance. We assume a fully ionized medium with primordial abundances to describe the electron density in terms of the gas density as

$$ n_{\text{e}}(r) = \frac{(X_{\text{H}} + 1) \rho_{\text{gas}}(r)}{2 m_{\text{amu}}}, \qquad (15) $$

with $X_{\text{H}} = 0.76$ being the hydrogen mass fraction and $m_{\text{amu}}$ the atomic mass unit. We evaluate equation (14) at the median redshift of the CMASS sample, $z = 0.55$. Future work will test the validity of this assumption, that is, that the kSZ signal does not have significant redshift evolution in the CMASS sample, and that our assumption is representative over the full redshift range we are sensitive to.

In order to calibrate the model for the kSZ profile, the mean halo mass of the CMASS galaxy sample is needed. Note that because the integrated kSZ signal scales with the gas mass and this quantity approximately tracks the halo mass, it is important that the theoretical predictions are for a sample with the same mean halo mass as the CMASS sample. Given the significant scatter in the literature, we choose to include an additional model parameter in the analysis, $M_{\text{h,200}}$, corresponding to the mean $M_{200}$ of the CMASS sample, with a prior range provided in Table 1. The justification for this prior choice, and an investigation of its impact are given in Appendix B3.

## 3.4 Modelling kSZ with baryonification

The bfc method provides a model for the gas density (equation 9) and describes the gas content as measured by kSZ (equations 13–15). We select this model to jointly analyse the kSZ measurements with the lensing data given that it is agnostic to any choice of hydrodynamical simulation.

In Fig. 2, we aim to build a better intuition on the BCEmu model by exploring the impact of each parameter on the suppression of the matter power spectrum, $P(k)/P_{\text{DMonly}}(k)$ (top row) and the kSZ radial
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Top: the impact of varying parameters of the bfc model, BCEmu, on the ratio of the total matter power spectrum compared to a dark-matter-only power spectrum, $P(k)/P_{DMonly}(k)$. Each panel varies one bfc parameter at a time within the prior bounds (reported in Table 2), whilst keeping the remaining six parameters at their fiducial fixed value, corresponding to the fit to a range of hydrodynamical simulations (Giri & Schneider 2021). For reference, various predictions for the suppression of the matter power spectrum from simulations are overplotted in black: FLAMINGO (Schaye et al. 2023, solid line); BAHAMAS (McCarthy et al. 2017, dashed-dotted line); SIMBA (Davé et al. 2019, dotted line); MillenniumTNG (Pakmor et al. 2023, dashed line); and FABLE (double-dotted-dashed line, Henden et al. 2018, Bigwood et al. in preparation). Bottom: the stacked kSZ temperature profile at 98 GHz as a function of angular radius, $\theta$, centred on the group or cluster (bottom) for the same bfc parameters as above. The ACT CMASS measurements at 98 GHz are shown as the black data points in the bottom panels and the model profiles are convolved with the f90 beam profile for comparison.</figcaption></figure>

temperature profile $T_{\text{kSZ}}(\theta)$ (bottom row). At a fixed cosmology,<fnref n="4" /> we test the dependence on each of the seven BCEmu parameters from left to right within the defined prior range, whilst fixing the remaining six model parameters at their fiducial fixed value (see Table 1).

The $M_c$ parameter controls the mass or proportion of groups and clusters that have had feedback-induced gas removal and has been previously identified as the most important in the model (Schneider & Teyssier 2015; Giri & Schneider 2021). Indeed, we see that within its prior range, this parameter modulates the amplitude, slope and extent of the suppression of the matter power spectrum at scales $k \gtrsim 0.1$ and the amplitude and slope of the kSZ profile. For clusters of mass greater than $M_c$, the slope of the gas profile (equation 10) tends to $\beta = 3$ and the gas profile approaches a truncated NFW profile. However, for groups of mass smaller than $M_c$, the slope of the power law decreases and the gas profile more closely resembles one that has experienced AGN feedback and had gas ejected from the halo. Hence, a larger value of $M_c$ results in a greater proportion of groups and clusters having gas profiles that mimic the effect of baryonic feedback, leading to a greater suppression of the matter power spectrum and simultaneously, a decrease in the amount of gas within an aperture centred on the galaxy, which gives a smaller kSZ signal.

In the reduced complexity BCEmu1 model, the remaining parameters are kept fixed, although it is clear that their value choice can have a significant impact on the predictions for the matter power spectrum and the kSZ signal. We observe that even at fixed $M_c$, a larger value of $\theta_{ej}$ increases the radius that gas is ‘ejected to’ in the gas profile and effectively causes matter to be redistributed on smaller $k$-scales. This corresponds to a decrease in the amount of gas as measured by the amplitude of the kSZ signal. Similarly, we find that increasing the exponent in the gas profile parametrization $\gamma$ results in a greater suppression of the power spectrum and lower kSZ

amplitude, although the impact is on smaller scales ($k \gtrsim 1$) than the effect of $\log_{10} M_c$ and $\theta_{ej}$. The remaining gas parameters, $\delta$ and $\mu$, have the reverse effect, that is, decreasing them results in an enhanced power spectrum suppression and a lower kSZ signal.

Decreasing the values of stellar parameters, $\eta$ and $\eta_\delta$, makes a less steep stellar-halo mass relation, so that a larger amount of stars condense out of the gas and form galaxies. This causes a boost in the matter power spectrum at large $k$, as we observe in Fig. 2. The kSZ radial temperature profile is only dependent on the gas density profile of groups and clusters (equation 15) and so not directly impacted by the stellar parameters $\eta$ and $\eta_\delta$. However, $\eta$ can indirectly have a small impact on the kSZ signal as it alters the number of stars in groups and clusters, and therefore the reservoir of available gas.

We note two limitations in our modelling of the kSZ signal. While we vary the halo mass of the kSZ sample in our analysis, we assume that the mass dependence of the model is valid beyond the range probed by the kSZ measurement and at varying redshift. The kSZ measurements span a mass range $M_{200} \approx [0.5-7] \times 10^{13} M_\odot$, which is similar to the halo mass range that cosmic shear is most sensitive to ($M_{200} \approx 10^{13.5} M_\odot$). Nevertheless, cosmic shear still has contributions from higher and lower mass haloes (e.g. Salcido et al. 2023).

Future kSZ measurements that span multiple mass and redshift bins will provide a better understanding of the dependence of the kSZ signal on galaxy properties and the suitability of this assumption. Another issue that we do not address in the present study is the role of centrals versus satellites in the observed stacked kSZ profiles. Our theoretical predictions correspond to central galaxies that are assumed to be perfectly centred within their host haloes and we may reasonably expect some degree of bias (with respect to theoretical predictions) to be introduced by the inclusion of satellites. Without forward modelling the BOSS CMASS selection function, which is beyond the scope of this work, it is difficult to predict the magnitude or sign of this effect. On the one hand, satellites will obviously be mis-centred with respect to their host haloes and one may expect this to lead to a lower kSZ signal. On the other hand, a stellar mass-based

<footnote n="4">⁴We chose the parameters $h = 0.742$, $\Omega_m = 0.255$, and $f_b = \Omega_b/\Omega_m = 0.166$.</footnote>
<table><caption>**Table 2.** Summary of cosmological, observational, and astrophysical parameters and priors used in the analysis. In the case of flat priors, the prior is bound to the range indicated in the 'value' column, while Gaussian priors are described by their mean and $1\sigma$ width.</caption><thead><tr><th>Parameter</th><th>Type</th><th>Value</th></tr></thead><tbody><tr><td><strong>Cosmological</strong></td><td></td><td></td></tr><tr><td>Ω<sub>m</sub>, total matter density</td><td>Flat</td><td>[0.1, 0.9]</td></tr><tr><td>Ω<sub>b</sub>, baryon density</td><td>Flat</td><td>[0.03, 0.07]</td></tr><tr><td>10<sup>-9</sup> A<sub>s</sub>, scalar spectrum amplitude</td><td>Flat</td><td>[0.5, 5.0]</td></tr><tr><td><em>h</em>, Hubble parameter</td><td>Flat</td><td>[0.55, 0.91]</td></tr><tr><td>n<sub>s</sub>, spectral index</td><td>Flat</td><td>[0.87,1.07]</td></tr><tr><td>Ω<sub>ν</sub> h<sup>2</sup>, neutrino mass density</td><td>Flat</td><td>0.06</td></tr><tr><td><strong>Observational</strong></td><td></td><td></td></tr><tr><td>Δz<sup>1</sup>, source redshift 1</td><td>Gaussian</td><td>(0.0, 0.018)</td></tr><tr><td>Δz<sup>2</sup>, source redshift 2</td><td>Gaussian</td><td>(0.0, 0.015)</td></tr><tr><td>Δz<sup>3</sup>, source redshift 3</td><td>Gaussian</td><td>(0.0, 0.011)</td></tr><tr><td>Δz<sup>4</sup>, source redshift 4</td><td>Gaussian</td><td>(0.0, 0.017)</td></tr><tr><td>m<sup>1</sup>, shear calibration 1</td><td>Gaussian</td><td>(-0.006, 0.009)</td></tr><tr><td>m<sup>2</sup>, shear calibration 2</td><td>Gaussian</td><td>(-0.020, 0.008)</td></tr><tr><td>m<sup>3</sup>, shear calibration 3</td><td>Gaussian</td><td>(-0.024, 0.008)</td></tr><tr><td>m<sup>4</sup>, shear calibration 4</td><td>Gaussian</td><td>(-0.037, 0.008)</td></tr><tr><td><strong>Intrinsic alignment</strong></td><td></td><td></td></tr><tr><td>a<sub>1</sub>, tidal alignment amplitude</td><td>Flat</td><td>[-5, 5]</td></tr><tr><td>η<sub>1</sub>, tidal alignment redshift index</td><td>Flat</td><td>[-5, 5]</td></tr></tbody></table>

selection implies that the satellites will typically be in hosts that are more massive than a host which has a central of similar stellar mass. This will tend to boost the kSZ signal. For the present study we neglect these uncertainties, leaving their careful consideration for a future study.

# 4 MODEL PIPELINE AND VALIDATION

In this section, we briefly describe the set up of the cosmological inference pipeline (Section 4.1). We validate the robustness of this pipeline using each of the four baryon models with a mock analysis, described in Appendix A. Here, we briefly motivate the choices made in the construction of the mock data (Section 4.2) and summarize the findings (Section 4.3).

## 4.1 Inference pipeline

To analyse the cosmic shear data, we build upon the public DES Y3 cosmological inference pipeline. We utilize one parameter for redshift and shear calibration per tomographic bin, with prior ranges set to those used in DES Y3. Without feedback models, we have 15 parameters that we marginalize over in the analysis. Table 2 summarizes the cosmological, observational, and astrophysical priors used in this work.

Parameters are estimated via nesting sampling using MULTINEST<fnref n="5" /> (Feroz, Hobson & Bridges 2009) within the COSMOSIS framework (Zuntz et al. 2015), with the sampler settings listed in Appendix D. We note however that DES & KiDS Collaborations (2023) demonstrate that MULTINEST can underestimate the 68 per cent confidence levels for $S_8$, at the level of ~10 per cent – 20 per cent, while POLYCHORD (Handley, Hobson & Lasenby 2015) is more accurate. In agreement with this finding, Appendix D reports that sampling with MULTINEST instead of POLYCHORD in a WL-only analysis with BCEmu7 results in a 18 per cent smaller 68 per cent confidence level

for $S_8$, and a 9 per cent smaller 68 per cent confidence level for $S_8$ in a joint WL and kSZ analysis with BCEmu7. We follow Amon et al. (2022) and Secco et al. (2022) when reporting the parameter constraints and quote the mean of the 1D marginal distribution, along with the 68 per cent confidence limit, which defines the area around the peak of the posterior within which 68 per cent of the probability lies.<fnref n="6" />

## 4.2 Mock data

In order to assess the robustness of our lensing inference pipeline, non-linear power spectrum model and baryon feedback models, we perform analyses using synthetic data and test the ability to recover unbiased cosmology. The mock data were created using the best-fitting cosmological parameters obtained from the DES Y3 joint lensing and clustering analysis (DES Collaboration 2022).<fnref n="7" />

We create two dark-matter-only mock data vectors, using different models for the non-linear matter power spectrum. The first uses HM20, the same model used to analyse the data throughout this work. This mock is important for testing that the analysis pipeline can accurately recover cosmological parameters before considering feedback effects. The second mock uses the EuclidEmulator2 (EE2), which has been shown to be accurate to 1 per cent for $k \le 10 h \text{ Mpc}^{-1}$ (Euclid Collaboration: Knabenhans et al. 2019, 2021; Adamek et al. 2023). This mock is used to ensure that HM20 provides a sufficiently accurate description for the non-linear matter power spectrum when compared to EE2, for the full angular scale range of the DES Y3 cosmic shear data.

To test the four baryon mitigation strategies, we chose three baryon feedback scenarios to build the mock data and their predictions for the suppression of the matter power spectrum are shown as the red lines in Fig. 1. These choices include extreme scenarios as they are designed to test the flexibility and limits of the baryon models, rather than an attempt to chose the most accurate prescription. First, we consider the upper limit of the BAHAMAS hydrodynamic simulation suite, with $T_{\text{heat}} = 8.0$ (McCarthy et al. 2017). Next, we consider a mock universe with a more extreme power spectrum prescribed by cosmo-OWLS $T_{\text{heat}} = 8.5$ (Le Brun et al. 2014), although this simulation does not replicate the local gas fractions in groups and clusters. As it is possible that hydrodynamical simulations do not capture the complexities of feedback, we want to test the ability of the baryon models to accurately capture a scenario that modulates the matter power spectrum with a different shape than that in typical simulations. We consider a mock with a suppression described by the $A_{\text{mod}}$ parameter (Amon & Efstathiou 2022; Preston et al. 2023), where $A_{\text{mod}} \approx 0.858$ is that required to reconcile the DES Y3 cosmic shear with Planck\&CDM cosmology.

## 4.3 Mock results

The results of these mock tests are detailed in Appendix A. Here, we summarize the findings, which are shown in Fig. A2:

(i) When analysing an EE2-generated dark-matter-only mock with the HM20 dark-matter-only model, $S_8$ is overestimated by ~0.4σ. (This is reduced to ~0.2σ when analysing restricted angular scale

<footnote n="6">⁶Formally, these are credible intervals; however, we choose to use the term 'confidence interval' in this paper to retain consistency with the language used in Amon et al. (2022) and Secco et al. (2022).</footnote>

<footnote n="7">⁷That is, with $S_8 = 0.7805$, $\Omega_m = 0.3380$ and $\sigma_8 = 0.7353$ (see Appendix A for more detail).</footnote>

<footnote n="5">⁵MULTINEST: https://github.com/farhanferoz/MultiNest</footnote>
measurements.) Furthermore, we find that $\Omega_m$ is underestimated by ~0.9$\sigma$. However, we identify projection effects in this parameter of ~0.5$\sigma$ by analysing an HM20 mock with an HM20 model. In this work, we focus on the $S_8$ parameter, and we note that further testing is needed to assess the reliability of the $\Omega_m$ constraints.

(ii) The DES Y3 'ACDM-optimized' scale cuts underestimate $S_8$ in all three mock baryonic feedback scenarios. This is as expected, as these scale cuts were defined with the OWLS-AGN scenario, which predicts less power suppression than BAHAMAS 8.0, cosmo-OWLS 8.5, and $A_{\text{mod}}$. While the scale cuts, by design, remove the sensitivity of the analysis to the impact of baryon feedback, this method's success relies on the true feedback scenario to be less extreme than the simulation used to define the cuts.

(iii) HM20, as used with their fiducial BAHAMAS-based prior, underestimates $S_8$ by ~0.7$\sigma$ for the BAHAMAS 8.0 mock, and by more than 1$\sigma$ when analysing a mock with a more extreme feedback scenario. Using a wide prior alleviates this, with the model recovering the true cosmology within ~0.5$\sigma$ for all mock scenarios, although with a cost of almost a factor of two in the precision of the $S_8$ constraint.

(iv) The SP(k) model, used with both a wide and X-ray prior, can recover the input cosmology to within ~0.2$\sigma$ for a BAHAMAS 8.0 mock. With the more extreme cosmo-OWLS 8.5 and $A_{\text{mod}}$ mocks, this model underestimates the value of $S_8$ by up to ~0.5$\sigma$. This is as expected, as both cosmo-OWLS 8.5 and $A_{\text{mod}}$ are outside the expected baryon fractions as compared to observations by Akino et al. (2022).

(v) When allowing only one parameter to vary in BCEmu emulator, BCEmu1, we recover $S_8$ to within ~0.2$\sigma$ for the BAHAMAS 8.0 mock, and ~0.5$\sigma$ for cosmo-OWLS 8.5. When we use the more flexible BCEmu7, all mock scenarios recover the true cosmology within ~0.5$\sigma$, and the error bar on $S_8$ is up to 1.5 times wider.

Overall, we find that when more restrictive modelling choices are used, we tend to underestimate $S_8$. The bias is worsened when restrictive choices are used to analyse mocks with the more extreme baryonic feedback scenarios, that is, cosmo-OWLS 8.5 and $A_{\text{mod}}$. We note that despite the greater accuracy of using more conservative priors or marginalizing over a greater number of baryonic feedback parameters, it is at the expense of the precision. We find that the uncertainty on the constraint of $S_8$ can degrade by up to a factor of two when switching to more flexible modelling choices.

# 5 RESULTS: ASSESSING MODELS FOR BARYONIC EFFECTS

The results of the WL-only DES Y3 analysis are divided into three sections. In Section 5.1, we present the headline cosmological constraints using the four baryon feedback strategies outlined in Section 3.2. In Section 5.2, we show the constraints on the suppression of the power spectrum. Finally, in Section 5.3, we explore the dependence of our results on the model complexity and prior choices within each strategy.

## 5.1 Cosmological parameter constraints

The 1D marginalized constraints obtained for $S_8$ are summarized in Fig. 3 for analyses using all model variants. Here, we compare those from the DES Y3 'ACDM optimized' scale cuts, HM20: BAHAMAS, SP(k): Xray and BCEmu7. For these approaches, the mean marginal values of $S_8$ are found with 68 per cent confidence

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Summary of the 1D marginalized constraints on $S_8$ from analysing the DES Y3 cosmic shear data with different baryonic feedback mitigation strategies. The mean of the $S_8$ marginalized posterior is indicated by the symbol and 68 per cent confidence levels are shown as horizontal bars. The primary result from the joint analysis of WL + kSZ using BCEmu7 is represented as the purple shaded region. We compare to the Planck TTTEE result presented in Efstathiou & Gratton (2021).</figcaption></figure>

levels to be

$$ 
\begin{align} 
\text{BCEmu7: } S_8 &= 0.818^{+0.017}_{-0.024} \nonumber \\ 
\text{Spk: Xray: } S_8 &= 0.806^{+0.015}_{-0.013} \tag{16} \\ 
\text{HM20: BAHAMAS: } S_8 &= 0.811^{+0.013}_{-0.014} \nonumber \\ 
\text{Scale cuts: } S_8 &= 0.805^{+0.018}_{-0.018}. \nonumber 
\end{align} 
$$

For reference, we show the result when all angular scales are analysed without any model for baryonic effects and the Planck TTTEE.⁸ACDM result (orange, Efstathiou & Gratton 2021). The 2D marginalized posteriors for $S_8$, $\Omega_m$, and $\sigma_8$ using DES Y3 WL data are plotted in the left panel of Fig. 4, also showing no baryon mitigation (dotted grey), DES Y3 'ACDM optimized' scale cuts (dashed green), HM20: BAHAMAS (navy), SP(k): Xray (solid green), BCEmu7 (pink), and Planck TTTEE (orange). Table 3 lists the mean constraints on $S_8$ and $\Omega_m$ and quantifies the goodness of fit to the data of each modelling variant by quoting the minimum and reduced $\chi^2$, $\chi^2_{\text{min}}$, and $\chi^2_{\text{red}}$. We find that each baryon feedback analysis variant demonstrates a suitable fit to the measurements.

When analysing all angular scales without modelling baryonic feedback, we attain a low value of $S_8 = 0.794^{+0.013}_{-0.011}$. This is up to 1$\sigma$ lower than constraints attained with modelled baryonic effects, highlighting the importance of mitigating feedback to avoid biased cosmology. When accounting for baryonic effects, we find, in agreement with the mock analysis for the models tested, that the measured $S_8$ is consistent at the level of 0.6$\sigma$ (~2 per cent).<fnref n="9" /> However, the error bar on the $S_8$ constraint varies by a factor of 1.5. In more detail, we see that HM20: BAHAMAS and SP(k): Xray

<footnote n="8">⁸TTTEE refers to the high multipole likelihood attained from combining the temperature power spectra (TT), temperature-polarization E-mode cross spectra (TE) and polarization E-mode power spectra (EE).</footnote>

<footnote n="9">⁹Note that throughout the following sections we quantify the shift in the measured value of $S_8$ by two analyses using the metric $\Delta S_8 / [(\sigma_{S_8}^1)^2 + (\sigma_{S_8}^2)^2]^{1/2}$, where $\Delta S_8$ is the difference between the respective mean values and $\sigma_{S_8}^1, \sigma_{S_8}^2$ are the 1$\sigma$ errors on $S_8$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Left: the marginalized posteriors for $\Omega_m$, $\sigma_8$, and $S_8$ using the DES Y3 cosmic shear data and different baryon feedback models. We compare the DES Y3 optimized scale cut approach (green dashed) to scenarios where all angular scales of the DES Y3 lensing data are modelled: HM20, using their recommended BAHAMAS-based prior (HM20: BAHAMAS, navy), the seven-parameter bfc model (BCEmu7, pink), and the SP(k) emulator, using their 'X-ray observational prior' (SPk: Xray, solid green). For reference, we show the case where all scales of the DES data vector are used but baryons are not modelled (grey dotted) and the Planck TTTEE likelihood (orange, Efstathiou & Gratton 2021). The inner and outer contours show the 68 per cent and 95 per cent confidence levels, respectively. Right: the corresponding constraints on the suppression of the total matter power spectrum compared to a dark-matter-only scenario, $P(k)/P_{DMonly}(k)$, at $z=0$ using the DES Y3 cosmic shear and again, HM20 (HM20: BAHAMAS, navy), BCEmu (BCEmu7, pink), and SP(k) (SPk: Xray, green). The solid lines show the mean suppression and the shaded regions indicate the 68 per cent confidence levels. For reference, various predictions from hydrodynamical simulations are overplotted in black: FLAMINGO (Schaye et al. 2023, solid line); BAHAMAS (McCarthy et al. 2017, dashed-dotted line); SIMBA (Davé et al. 2019, dotted line); MillenniumTNG (Pakmor et al. 2023, dashed line); and FABLE (double-dotted-dashed line, Henden et al. 2018; Bigwood et al. in preparation).</figcaption></figure>

<table><caption>**Table 3.** The constraints on $S_8$, $\Omega_m$ attained for each method of baryonic feedback mitigation. We report the mean value of each parameter with errors given by the 68 per cent confidence levels. We also demonstrate the quality of the fit by reporting $\chi^2_{min}$, the minimum value of $\chi^2$, for each analysis variant. We report $\chi^2_{red} = \chi^2_{min}/N_{dof}$, with $N_{dof}$ being the number of degrees of freedom $N_{dof} = N_{dp} - N_{param}$, where $N_{dp}$ and $N_{param}$ are the number of data points and model parameters utilized in the analysis. With the exception of the scale cuts method, angular scales down to 2.5 arcmin of the DES Y3 data were used.</caption><thead><tr><th>Model</th><th>S<sub>8,mean</sub></th><th>Ω<sub>m,mean</sub></th><th>χ<sup>2</sup><sub>min</sub></th><th>N<sub>dp</sub></th><th>N<sub>param</sub></th><th>N<sub>dof</sub></th><th>χ<sup>2</sup><sub>red</sub></th></tr></thead><tbody><tr><td>No baryon mitigation</td><td>0.794<sup>+0.013</sup><sub>-0.011</sub></td><td>0.278<sup>+0.033</sup><sub>-0.044</sub></td><td>418.37</td><td>400</td><td>15</td><td>385</td><td>1.09</td></tr><tr><td>Scale cuts</td><td>0.805<sup>+0.018</sup><sub>-0.018</sub></td><td>0.281<sup>+0.035</sup><sub>-0.051</sub></td><td>284.86</td><td>273</td><td>15</td><td>258</td><td>1.10</td></tr><tr><td>HM20 Θ<sub>AGN</sub> = 7.6 – 8.0</td><td>0.811<sup>+0.013</sup><sub>-0.014</sub></td><td>0.261<sup>+0.026</sup><sub>-0.034</sub></td><td>415.25</td><td>400</td><td>16</td><td>384</td><td>1.08</td></tr><tr><td>HM20 Θ<sub>AGN</sub> = 7.3 – 9.0</td><td>0.822<sup>+0.017</sup><sub>-0.026</sub></td><td>0.252<sup>+0.027</sup><sub>-0.039</sub></td><td>415.82</td><td>400</td><td>16</td><td>384</td><td>1.08</td></tr><tr><td>BCEmu1</td><td>0.804<sup>+0.016</sup><sub>-0.017</sub></td><td>0.274<sup>+0.033</sup><sub>-0.042</sub></td><td>414.49</td><td>400</td><td>16</td><td>384</td><td>1.08</td></tr><tr><td>BCEmu3</td><td>0.814<sup>+0.019</sup><sub>-0.021</sub></td><td>0.261<sup>+0.029</sup><sub>-0.044</sub></td><td>414.97</td><td>400</td><td>18</td><td>382</td><td>1.09</td></tr><tr><td>BCEmu7 WL</td><td>0.818<sup>+0.017</sup><sub>-0.024</sub></td><td>0.255<sup>+0.027</sup><sub>-0.038</sub></td><td>414.21</td><td>400</td><td>22</td><td>378</td><td>1.10</td></tr><tr><td>BCEmu7 WL + kSZ</td><td>0.823<sup>+0.019</sup><sub>-0.020</sub></td><td>0.250<sup>+0.025</sup><sub>-0.036</sub></td><td>439.33</td><td>418</td><td>23</td><td>395</td><td>1.11</td></tr><tr><td>BCEmu7, Σm<sub>v</sub> : [0.06, 0.6]</td><td>0.813<sup>+0.019</sup><sub>-0.023</sub></td><td>0.269<sup>+0.028</sup><sub>-0.043</sub></td><td>414.64</td><td>400</td><td>23</td><td>377</td><td>1.10</td></tr><tr><td>BCEmu7, TATT</td><td>0.802<sup>+0.028</sup><sub>-0.024</sub></td><td>0.239<sup>+0.021</sup><sub>-0.041</sub></td><td>408.54</td><td>400</td><td>25</td><td>375</td><td>1.09</td></tr><tr><td>SP(k) conservative prior</td><td>0.817<sup>+0.015</sup><sub>-0.019</sub></td><td>0.255<sup>+0.025</sup><sub>-0.040</sub></td><td>415.44</td><td>400</td><td>18</td><td>382</td><td>1.09</td></tr><tr><td>SP(k) Xray et al. (2022) prior</td><td>0.806<sup>+0.015</sup><sub>-0.013</sub></td><td>0.261<sup>+0.025</sup><sub>-0.036</sub></td><td>415.02</td><td>400</td><td>18</td><td>382</td><td>1.09</td></tr></tbody></table>

give the tightest constraints on the $S_8$ parameter, which are in excellent agreement with each other. This is expected, as these models are calibrated on hydrodynamical simulations informed by X-ray constraints. The BCEmu7 analysis gives the highest value of $S_8$. As in the case of the mock analysis, this supports our findings that

restrictive modelling choices for baryonic feedback leads to lower value of $S_8$ when compared to more flexible models. The flexibility of BCEmu7 comes at a cost, as the uncertainty on $S_8$ is a factor of 1.5 larger than that in the HM20: BAHAMAS analysis. This is to be expected given the degeneracy between the extremity of feedback and
$S_8$; greater flexibility in the modelling of baryons inevitably results in a larger error bar on the $S_8$ constraint. The analysis using the DES Y3 scale cuts does not suffer as substantial a loss in constraining power as the BCEmu7 case, but results in the lowest value for $S_8$.$^{10}$

## 5.2 Power suppression constraints

Baryon feedback processes modify the gravitational evolution of the cosmic density field and suppress the matter power spectrum compared to a dark-matter-only scenario on non-linear scales, as seen in hydrodynamic simulations. This effect has been previously observed by analyses of WL data using variations of the bfc model (Schneider et al. 2022; Aricò et al. 2023; Chen et al. 2023). In this section, we constrain the amplitude and scale dependence of the suppression of the matter power spectrum due to baryonic effects, $P(k)/P_{DMonly}(k)$, using the DES Y3 cosmic shear. For the first time, we show the model dependence of the constraints by considering the model complexity of the bfc model and the comparison to the SP(k) and HM20 models. For each baryonic feedback model, we record the power spectrum suppression at each step in the chain.

In the right panel of Fig. 4 we plot the mean suppression and the 68 per cent confidence levels inferred from analyses with our three baseline models: BCEmu7, HM20: BAHAMAS ($\Theta_{AGN} = 7.6 - 8.0$), and SP(k): Xray. We find that the suppression inferred by the three models are consistent within the 68 per cent confidence limits up to $k \approx 3 \ h\ Mpc^{-1}$. However, we find that BCEmu7 allows more extreme suppression of the power spectrum at all non-linear scales. There are substantial differences in the size of the uncertainties, correlated with the flexibility of the model. BCEmu7, the most flexible model, has the largest uncertainty. SP(k): Xray provides the tightest constraints on the power spectrum suppression and constrains a less extreme feedback scenario in terms of the amplitude and the scale extent of the suppression.

We compare our constraints to five hydrodynamical simulations: FLAMINGO, (Schaye et al. 2023, solid line); BAHAMAS (McCarthy et al. 2017, dashed-dotted line); SIMBA (Davé et al. 2019, dotted line); MillenniumTNG (Pakmor et al. 2023, dashed line); and FABLE (double-dotted-dashed line, Henden et al. 2018, Bigwood et al. in preparation). HM20 predicts a feedback strength that encompasses BAHAMAS $T_{\text{heat}} = 7.8$, which is unsurprising, given that the model is calibrated to span these simulations, but it is notable that the mean constraint is more extreme on all scales. Interestingly, all three models find the mean suppression to be more extreme than FLAMINGO on scales $k \sim 0.2-4 \ h^{-1}\ Mpc$, at the level of $1.1\sigma$ (BCEmu7), $3.0\sigma$ (SP(k): Xray), and $1.4\sigma$ (HM20: BAHAMAS) at $k = 2 \ h^{-1}\ Mpc$.

The power suppression constraints are broadly consistent with those from previous WL analyses (Schneider et al. 2022) and slightly more extreme than the constraints of Aricò et al. (2023), Chen et al. (2023), García-García et al. (2024), and Terasawa et al. (2024). Here, we make note of some details. Owing to the enhanced statistical power of the DES Y3 data over that of the Kilo-Degree Survey, we

find substantially improved constraints from the WL-only analysis compared to Schneider et al. (2022), even though we include an additional parameter for the IA model and eight additional nuisance parameters to account for the uncertainty in the shear and redshift calibration. Furthermore, we use a different bfc model to that of Chen et al. (2023) and Aricò et al. (2023) and we do not impose X-ray priors on any of the baryon parameters. One way in which this model is different from BCEmu is that only particles within the virial radius, $R_{200}$, are displaced (see Grandis et al. 2024, for a more detailed discussion of the model comparison). One implication of the model differences is that it is not straightforward to compare constraints on the $M_c$ parameter. We note however that in our choice to use BCEmu over the bfc emulator used in Chen et al. (2023) and Aricò et al. (2023), BACCOemu, the posterior on the $M_c$ parameter is not limited by the upper value of the prior. This may allow our analysis to constrain more extreme feedback scenarios than analyses with BACCOemu, and therefore constrain a greater matter power spectrum suppression than Aricò et al. (2023).

## 5.3 Impact of model complexity

The mock analysis revealed that marginalizing over a greater number of baryonic nuisance parameters, or utilizing wider priors on these parameters, generally improved the accuracy of the cosmological constraints. However, we saw that this was at the expense of inflated errors on the cosmological parameters, which is clearly suboptimal for an effective cosmological analysis. In this section, we explore the impact of altering the complexity and prior choices of each baryon feedback model on the measured cosmological and baryonic constraints when analysing the DES Y3 $\xi_{\pm}$ measurements.

The upper left panel of Fig. 5 shows the impact of the BAHAMAS-informed prior on the HM20 feedback parameter $\Theta_{AGN}$ on the marginalized $S_8$ posterior. In the light blue constraint, we extend the prior range for $\Theta_{AGN}$ outside of the calibration range to encompass more extreme feedback scenarios, as well as a dark-matter-only scenario. The parameters are degenerate and opening up the $\Theta_{AGN}$ prior leads to long tails that extend to higher values of $S_8$, such that the mean constraint is > $0.5\sigma$ higher. This illustrates how high values of $S_8$ are disfavoured by the restricted prior on the baryonic feedback model, and suggests that WL data may favour a higher value of $\Theta_{AGN}$ than the BAHAMAS simulations span. (Although we note that the HM20 model was only calibrated within the BAHAMAS range, so the mapping between the power spectrum suppression and the baryon fraction outside this range is uncertain.) Similarly, for the case shown in the right-hand panel of Fig. 5, when the X-ray informed prior is lifted, the data constrains higher values of $S_8$ by ~ $0.5\sigma$ (the posteriors on the SP(k) parameters are shown in Appendix B1). These shifts are consistent with those we determine in the mock analysis (Section 4.2, Appendix A). In Section 7, we discuss the implications of these results on our understanding of the gas models and observations.

We test the impact of limiting the bfc model complexity to the one (BCEmu1) and three-parameter (BCEmu3) case, compared to the fiducial BCEmu7. The central panel of Fig. 5 shows the marginalized posteriors on $\Omega_m$, $S_8$, and the bfc parameter $\log_{10} M_c$. The posteriors on all of the baryonic feedback parameters are shown in Appendix B2. BCEmu7 and BCEmu3 produce comparable constraints on cosmological and feedback parameters, generating a shift in $S_8$ of only ~ $0.2\sigma$. We also do not see any significant improvement of the precision on these constraints when marginalizing over four less baryonic parameters. However, switching to BCEmu1 from BCEmu7 results in a substantially lower values of $S_8$ by > $0.5\sigma$, consistent

<footnote>$^{10}$Note that the DES Y3 ‘ΛCDM optimized’ cosmic shear analysis (Amon et al. 2023; Secco et al. 2022) obtains a lower value of $S_8$ than that obtained here ($S_8 = 0.772_{-0.018}^{+0.018}$). Based on the study of the impact of analysis choices in (DES & KiDS Collaborations 2023), we attribute the difference in our results primarily to the use of HM20 to model the dark matter non-linear matter power spectrum, which was shown to be more accurate than HALOFIT, as well as the IA model that we chose for this analysis, and the choice to fix the neutrino mass in the analysis. In Appendix C, we investigate the impact of these choices further.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Top: the degeneracy of the marginalized $S_8$ posterior with varying the baryon model complexity or prior choices on baryon model parameters. For each panel, the inner and outer contours show the 68 per cent and 95 per cent confidence levels, respectively. Left: $\text{HM20} - S_8$ and $\Theta_{\text{AGN}}$ attained when using a prior bracketing the BAHAMAS simulations (dark blue, $\theta_{\text{AGN}} = 7.6 - 8.0$) and a less-informative prior choice, (light blue, $\theta_{\text{AGN}} = 7.3 - 9.0$). Centre: $\text{BCEMu} - S_8$ and $\log_{10}(M_c)$ attained for the full seven parameter (pink, BCEmu7), three parameter (purple, BCemu3), and one parameter (blue, BCEmu1) models. Right: $\text{SP}(k) - S_8$ and $\beta$ attained for the X-ray informed model (dark green) and a less-informative prior (light green). Bottom: the corresponding constraints on the suppression of the total matter power spectrum compared to a dark-matter-only scenario, $P(k)/P_{\text{DMonly}}(k)$, at $z=0$ for each of the models. In general, when allowing for a more flexible model, the constraints indicate more extreme suppression of power, although with degraded constraining power. For each panel, the solid lines show the mean suppression predicted, and the shaded regions the 68 per cent confidence levels. For reference, various predictions for the suppression of the matter power spectrum from simulations are overplotted in black: FLAMINGO (Schaye et al. 2023, solid line); BAHAMAS (McCarthy et al. 2017, dashed-dotted line); SIMBA (Davé et al. 2019, dotted line); MillenniumTNG (Pakmor et al. 2023, dashed line); and FABLE (double-dotted-dashed line, Henden et al. 2018; Bigwood et al. in preparation).</figcaption></figure>

with the mock analysis. This reduction in the flexibility of the model forces $\log_{10} M_c$ to a larger value, since the full extremity of feedback has to be captured by only the one parameter. We attribute this shift in the cosmology and feedback parameters to the values at which the remaining six ‘under-the-hood’ bfc parameters are fixed to in the BCEmu1 model, quoted in Table 1.<fnref n="11" /> Given that the analysis using BCEmu1 gives a larger value of $\log_{10} M_c$ than that attained using BCEmu7, and that for example, the posterior on $\theta_{\text{ej}}$ is toward higher values than the fixed 3.5, this implies that the values for the simulation-informed fixed parameters of BCEmu1 represent a less-extreme feedback scenario than those constrained by an analysis of DES cosmic shear.

The lower panels of Fig. 5 illustrate how the corresponding constraints on the power spectrum suppression are sensitive to the restrictiveness of the choices within each model. Analogously to the upper panels, the left and rightmost panels test switching to the

more conservative prior choices of the HM20 and SP(k) models, that is, $\Theta_{\text{AGN}} = 7.3 - 9.0$ and a prior spanning the feedback landscape of the ANTILLES suite, respectively. The central panel shows the impact of restricting the model complexity of BCEmu, by allowing only one (BCEmu1) out of the full seven (BCEmu7) bfc parameters to vary in the analysis. As before, we compare the data constraints to predictions from hydrodynamical simulations.

For each model, when we opt for the more flexible modelling choices, the mean constraint on the power spectrum suppression tends to more extreme scenarios at all non-linear $k$ with respect to their more restrictive counterpart. In particular, each of the flexible models is most consistent with more extreme feedback scenarios (e.g. SIMBA). This could suggest that the restrictive baryonic modelling choices do not have the flexibility to capture the full extremity of feedback that the data prefer and therefore the higher values of $S_8$. It is clear, however, that the WL data cannot place strong constraints on feedback on its own. The use of the more conservative model choices comes at the expense of reduced constraining power. For example, SP(k) displays over a factor of two increase in the size of the 68 per cent confidence level at all scales when switching from the restrictive to conservative modelling choice. A complementary

<footnote n="11">¹¹The fixed parameters were determined by fitting to the baryonic suppression of the matter power spectrum of a number of hydrodynamical simulations, then summing the likelihoods to find the best-fitting parameters to all of the simulations (Giri & Schneider 2021).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Left: the marginalized posteriors for $\Omega_m$ and $S_8$ attained modelling baryonic feedback with the BCEmu7 model using the DES Y3 cosmic shear data set only (pink), or a combined analysis of DES Y3 cosmic shear and ACT DR5 kSZ measurements (purple). The inner and outer contours show the 68 per cent and 95 per cent confidence levels, respectively. We compare to the CMB $\Lambda$CDM constraint measured by Planck TTTEE E likelihood (Efstathiou & Gratton 2021). Right: the constraints on the suppression of the total matter power spectrum compared to a dark-matter-only scenario, $P(k)/P_{DMonly}(k)$, at $z=0$ when modelling baryonic feedback with the BCEmu7 model using the DES Y3 cosmic shear data set only (pink), or a combined analysis of DES Y3 cosmic shear and ACT DR5 kSZ measurements (purple). The solid lines show the mean suppression and the shaded regions indicate the 68 per cent confidence levels. For reference, various predictions for the suppression of the matter power spectrum from simulations are overplotted in black: FLAMINGO (Schaye et al. 2023, solid line); BAHAMAS (McCarthy et al. 2017, dashed-dotted line); SIMBA (Davé et al. 2019, dotted line); MillenniumTNG (Pakmor et al. 2023, dashed line); and FABLE (double-dotted-dashed line, Henden et al. 2018; Bigwood et al. in preparation).</figcaption></figure>

avenue to simulation-based models that still maximizes cosmological constraints is to use the flexible model framework and jointly analyse the lensing with observations of the gas.

# 6 RESULTS: JOINT WEAK LENSING + KINETIC SZ

In the previous section, we demonstrated that WL $S_8$ constraints are degenerate with the amount of baryonic feedback. A complementary approach to simulation-informed baryon models is to use a flexible model to jointly analyse cosmic shear with probes of the gas distribution in order to better constrain the model parameters. In this section, we report the results of a joint analysis of the DES Y3 cosmic shear and ACT kSZ measurements, using the BCEmu7 baryon mitigation model, described in Section 3.3.

## 6.1 Cosmological parameter constraints

The left panel of Fig. 6 compares the marginalized posteriors on $S_8$ and $\Omega_m$ attained from the WL-only analysis using BCEmu7 (pink) to those obtained from a joint analysis with kSZ (purple). The Planck TTTEE E $\Lambda$CDM posteriors are shown for reference (Efstathiou & Gratton 2021). The mean marginal value of $S_8$ for the joint analysis is found with 68 per cent credible levels to be

$$ \text{Lensing + kSZ } S_8 = 0.823^{+0.019}_{-0.020} \qquad (17) $$

which corresponds to a ~ 0.2$\sigma$ shift towards higher values with respect to the result of the WL-only BCEmu7 analysis of $S_8 = 0.818^{+0.017}_{-0.024}$.

Amon & Efstathiou (2022) and Preston et al. (2023) have proposed that the $S_8$ tension could be resolved if the non-linear matter power spectrum is suppressed more strongly than is currently assumed in WL analyses, either due to unmodelled baryonic feedback effects or non-standard dark matter. For baryonic feedback to be the source, its effects on the matter power spectrum would be more extreme than is currently predicted by the hydrodynamical simulations. The impact of the kSZ on the $S_8$ constraint to reduce the uncertainty and shift the value toward Planck is minor given the low signal-to-noise ratio of the kSZ measurements. However, it motivates us to investigate the impact of the kSZ on the $P(k)$ constraints and discuss the possibility of a more extreme suppression.

We find that incorporating a joint analysis with kSZ results in a significant improvement in the constraint on $\log_{10} M_c$, reducing the uncertainty by a factor of ~3 with respect to the WL-only analysis. The joint kSZ and WL data prefer larger values of $\theta_{ej}$ and $\gamma$, and lower values of $\mu$ (see Fig. B2 and Table B2), suggesting that gas is ejected to larger radii, redistributing matter on larger scales (Fig. 2). This supports the idea of a more extreme feedback scenario, resulting in higher values of $S_8$. The improvement on the WL constraint on $S_8$ with the inclusion of kSZ is ~ 10 per cent. Although this is modest, it is clear that the parameter space of the bfc model is better constrained, even in this case of a kSZ measurement with signal-to-noise ratio of ~7.

We compare the best-fitting models to the measured DES Y3 cosmic shear two-point correlation functions, $\xi_{\pm}$ in Fig. 7. The dark green line indicates the best fit for our reanalysis of DES Y3 using their scale cuts, which are indicated by the green shaded region. The purple line shows that for the BCEmu joint analysis of all angular
<figure><img src="image_7.png" /><figcaption>**Figure 7.** The DES Y3 cosmic shear two-point correlation function measurements, $\xi_{\pm}$, as a function of angular separation, $\theta$, (black data points, upper panels) from Amon et al. (2022) and Secco et al. (2022) for each pair of redshift bins, indicated by the panel label. The measurements are scaled by $\theta$ for visual aid. The error bars are calculated as the square root of the diagonal of the analytic covariance matrix. We show the best-fitting $\Lambda$CDM theoretical predictions to the large angular scale measurements (indicated by the shaded regions) without a baryon model (green line) and the best fit obtained from jointly analysing all angular scales of the DES shear measurements with ACT kSZ data using the BCEmu7 model (purple line). The residuals between the measurement and best-fitting model predictions $\delta\xi_{\pm}/\xi_{\pm}$ are shown in the lower panels.</figcaption></figure>

scales of the lensing measurement and the kSZ. The lower panels highlight the fractional residuals between the measurements and the model, $(\xi_{\pm} - \xi_{\pm}^{\text{model}})/\xi_{\pm}^{\text{model}}$, following the same colour scheme. While the fits are indistinguishable at large scales, at small scales, particularly for $\xi_{-}$, the predictions differ and the best-fitting line for the joint analysis has a lower amplitude. Both of these model choices provide a good fit to the data, although their $S_8$ values differ by $\sim 1\sigma$ and their non-linear matter spectrum predictions differ substantially. This highlights the degeneracy between a low-$S_8$ cosmology and a higher-$S_8$ cosmology with baryonic effects modelled on non-linear scales. Fig. 8 displays the joint WL + kSZ BCEmu7 constraints on the stacked kSZ radial temperature profile at 98 and 150 GHz (purple). We verify that like the DES data vector, the WL + kSZ model provides a good fit to the data. The $\chi^2_{\text{red}}$ values reported in Table 3 further show that the best-fitting models attained from the shear and WL + kSZ analyses are an equally good fit to the data sets.

It is important to note that the constraints on $S_8$ that we obtain are dependent on the choice of NLA as the IA model. Appendix C tests the impact of using the Tidal Alignment and Tidal Torquing ('TATT'; Blazek et al. 2019) superspace IA model. In a WL-only analysis with BCEmu7 and using TATT, we find a value of $S_8 \sim 0.5\sigma$ lower than that we obtain in our fiducial analysis with NLA. This shift is found to be consistent across the various baryon models that we consider,

as well as with previous findings (e.g. Secco et al. 2022; DES & KiDS Collaborations 2023). Therefore, when considering the results of this work in the context of the $S_8$ tension, it is important to bear in mind the existing uncertainty in IA modelling and the shifts in cosmological parameters that can occur as a result.

## 6.2 Power suppression constraints

We investigate how the constraint on the suppression of the matter power spectrum changes with the addition of the kSZ data. The right panel of Fig. 6 shows the mean constraint and 68 per cent confidence level for the BCEmu7 cosmic shear analysis (pink), as shown previously, compared to the joint WL and kSZ analysis (purple). The joint analysis results in a mean suppression that is more extreme at all displayed $k$-scales. For example, at $k = 2\,h\,\text{Mpc}^{-1}$, the $1\sigma$ bounds on the power suppression from the WL-only analysis range from 5 per cent to 20 per cent, and with the inclusion of the kSZ, the suppression ranges from 10 per cent to 25 per cent. Consistent with Schneider et al. (2022), we find that the inclusion of kSZ data in the analysis improves the constraint from WL only to favour more extreme scenarios. Compared to the hydrodynamical simulations, the mean suppression from the joint analysis is more extreme than all of MillenniumTNG, BAHAMAS, FLAMINGO, FABLE, and SIMBA.
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Stacked ACT DR5 kSZ temperature profile measurements (black points) as a function of angular radius, $\theta$, at 98 and 150 GHz (Schaan et al. 2021) and the best fit from the joint analysis with the DES Y3 cosmic shear using the BCEmu7 model (purple). The model profiles are convolved with the f90 and f150 beam profile for comparison to the 98 and 150 GHz data, respectively.</figcaption></figure>

at $k = 2h \, \text{Mpc}^{-1}$ by more than ~2.5$\sigma$, ~1.4$\sigma$, ~1.9$\sigma$, ~2.4$\sigma$, and ~0.9$\sigma$ respectively. Of the simulations, SIMBA shows the best agreement with the WL and kSZ constraint. We note, however, that for BAHAMAS and FLAMINGO we compare to only their fiducial feedback variants here, but both suites have more extreme feedback variations which are in better agreement with our measurements.

Finally, we consider our findings in the context of Amon & Efstathiou (2022) and Preston et al. (2023), which proposed that the $S_8$ tension could be resolved if a more extreme baryon feedback scenario than that predicted by the state-of-the-art hydrodynamical situations existed. These works analysed the KiDS and DES WL data assuming the Planck cosmology on linear scales, and modulating the non-linear power spectrum suppression via the $A_{\text{mod}}$ parameter. Preston et al. (2023) therefore quantified the small-scale suppression required to resolve the suppression, which corresponds to $A_{\text{mod}} \approx 0.858$. In this scenario, the matter power suppression is suppressed enough to reconcile the difference in $S_8$ between DES Y3 cosmic shear and the Planck $\Lambda$CDM model. The joint WL + kSZ constraints on the suppression of the matter spectrum that we obtain are consistent with the $A_{\text{mod}} = 0.858$ prediction.

## 7 DISCUSSION: CONSISTENCY OF X-RAY AND KINETIC SZ DATA

Our WL + kSZ approach is complementary to existing efforts toward the goal of a complete model of baryonic feedback that is consistent with a wide range of observables. X-ray measurements of cluster gas mass fractions are the primary observable used to calibrate or benchmark many hydrodynamical simulations (e.g. McCarthy et al. 2017; Henden et al. 2018; Kugel et al. 2023). Observations of the thermal Sunyaev-Zel'dovich (tSZ) power spectrum and tSZ flux-halo

mass relation have also been used to assess the simulations' realism (Henden et al. 2018; McCarthy et al. 2023; Pakmor et al. 2023; Schaye et al. 2023). In general, while tSZ and X-ray measurements are more sensitive to the inner regions of galaxy groups and clusters, the kSZ effect is well suited to probe the outskirts of haloes, through its sensitivity to low-density and low-temperature environments (Schaan et al. 2021). Furthermore, X-ray and tSZ observations are typically derived from massive cluster haloes ($M_{500} > 10^{14}$ M$_\odot$), while our kSZ measurements represent haloes of mass $M_{200} \sim 10^{13}$, closer to the mass range that WL is most sensitive to. Since the kSZ effect probes the gas in haloes of a different mass regime and on different scales to that which is currently used to calibrate feedback effects in simulations, it may allow new insights to be gained.

Our WL constraints on the matter power spectrum suppression point to a feedback scenario that is more extreme than most simulations predict. This observation holds in all three flexible model scenarios tested (Fig. 5). The addition of the kSZ data pushes the mean constraint towards an even more disruptive feedback scenario: at $k = 2h \, \text{Mpc}^{-1}$, the fiducial FLAMINGO simulation is disfavoured at ~2$\sigma$ and MillenniumTNG at ~2.5$\sigma$. Beyond the comparison to the simulations, our findings point to a more disruptive feedback scenario than inferred from the predicted $P(k)$ suppression using X-ray gas and stellar fraction observations (Grandis et al. 2024), as well as tSZ-mass relation of clusters (To et al. 2024), which could point to interesting differences between the X-ray view of baryon feedback, compared to that from WL and kSZ. We note that a strong feedback scenario was also found using an analysis of the cross-correlation of diffuse X-ray and WL (Ferreira et al. 2024), and hints of a stronger scenario with the cross-correlation of tSZ with galaxies (Pandey et al. 2023).

### 7.1 Mass dependence of the total baryon fraction

As discussed in Section 3.2.3, the baryon fraction, $f_b$, measured in haloes of mass $M_{500} \approx 10^{14} M_\odot$ can be related directly to the matter power spectrum suppression in a manner which is robust to a number of baryonic feedback prescriptions (van Daalen et al. 2020). In this section, we scrutinize the WL constraints on the halo mass dependence of the total baryon fraction in haloes with respect to the cosmic baryon fraction, $f_b/(\Omega_b/\Omega_m)$. For a given set of parameters, we can use the BCEmu model to compute the baryon fraction through the summation of the integrated gas and stellar profiles (see equations 9 and 11) out to $R_{500}$. For both the WL-only and WL + kSZ BCEmu7 analyses, we compute the $f_b/(\Omega_b/\Omega_m)-M_{500}$ relation for the cosmology and baryonic feedback parameters sampled at each step of the chain.

Fig. 9 shows the constraints for the WL-only analysis in pink and the WL + kSZ in purple, with the mean relation indicated by a solid line and 68 per cent confidence level as a shaded region. Measurements of halo masses are generally reliant on the assumption of hydrostatic equilibrium and hence may suffer a bias due to non-thermal pressure support in the halo (see e.g. Rasia et al. 2006). However, estimates have also been obtained from WL data which do not require this assumption (Hoekstra et al. 2015; Eckert et al. 2016; Mulroy et al. 2019; Akino et al. 2022). As a result we compare to the $1\sigma$ region of the X-ray Hyper Suprime-Cam (HSC)-XXL constraint (Akino et al. 2022) in hatched black. The X-ray data have a dependence on cosmology through $E(z)$, hence we scale the data to the mean cosmology constrained by the WL + kSZ analysis. We also show the baryon fractions measured in BAHAMAS (McCarthy et al. 2017), FLAMINGO (Schaye et al. 2023), MillenniumTNG (Pakmor
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Left: the total baryon fraction, $f_b/\Omega_b/\Omega_m$, as a function of halo mass, $M_{500}$, where $\Omega_b/\Omega_m$ is the cosmic baryon fraction and $f_b$ is the fraction of mass in baryons to the total halo mass in groups and clusters. We plot the constraint attained from analysing the DES Y3 cosmic shear data with BCEmu7 (pink) in addition to the constraint attained from the joint WL + kSZ analysis (purple). Solid lines show the mean baryon fraction halo mass relation, and the shaded regions enclose the $1\sigma$ constraints. Right: the fraction of mass in gas to the total halo mass in groups and clusters, $M_{\text{gas,500}}/M_{500}$, as a function of halo mass, $M_{500}$. We plot the constraint attained from analysing the DES Y3 cosmic shear data with BCEmu7 (pink) in addition to the constraint attained from the joint WL + kSZ analysis (purple). Solid lines show the mean gas fraction halo mass relation, and the shaded regions enclose the $1\sigma$ constraints. In both panels, we plot the X-ray constraints from HSC-XXL $1\sigma$ (Akino et al. 2022) as the black hatched region, scaled to the mean cosmology obtained from the WL + kSZ analysis. For reference, we also plot the measurements from BAHAMAS (McCarthy et al. 2017, dashed-dotted line); FLAMINGO (Schaye et al. 2023, solid line); MillenniumTNG (Pakmor et al. 2023, dashed line); and FABLE (double-dotted-dashed line, Henden et al. 2018; Bigwood et al. in preparation).</figcaption></figure>

et al. 2023), and FABLE (Henden et al. 2018, Bigwood et al. in preparation).

Overall, the BCEmu7 WL-only and WL + kSZ constraints are in good agreement. The addition of the kSZ reduces the uncertainty by a factor of ~1.5 at $M_{500} \approx 10^{13} M_\odot$ and prefers lower baryon fractions for all masses. For large groups with $M_{500} \gtrsim 5 \times 10^{13} M_\odot$, we find good consistency between the 68 per cent confidence levels on $f_b/(\Omega_b/\Omega_m)$ as predicted by BCEmu7 WL-only and WL + kSZ analyses compared to the X-ray data. However, for lower mass groups with masses $M_{500} \approx 10^{13} M_\odot$, the WL + kSZ data prefers a total baryon fraction that is lower than the X-ray data by ~1.4$\sigma$. With the exception of BAHAMAS which is in best agreement, all of the simulations predict higher values of $f_b/(\Omega_b/\Omega_m)$ for all masses compared to the mean of the WL + kSZ constraint. At $M_{500} \sim 4 \times 10^{13} M_\odot$, FLAMINGO predicts ~1.5$\sigma$ higher values, MillenniumTNG ~3.1$\sigma$ and FABLE ~1.4$\sigma$.

In Appendix B4, we compare the relationship between the baryon fraction and power suppression that our WL + kSZ constrains to that proposed by van Daalen et al. (2020), which is a good fit to many of the hydrodynamical simulations. We note that the BCEmu model does not impose a prior on this relationship and therefore provides a route to place constraints using data. Given that our constraints differ from the simulation-based relationship, we can conclude that either there are unaccounted for systematics in the WL and kSZ data, the BCEmu model allows non-physical scenarios, or simulations do not currently capture the full possible range of feedback effects, thereby overestimating the relationship (see e.g. Debackere, Schaye & Hoekstra 2020). Future work is needed to further understand this relationship using data.

The SP(k) baryon model also provides a direct mapping from the matter power spectrum suppression to the baryon fraction. In Appendix B1, we show the mean and 68 per cent confidence levels for the baryon fraction attained from the WL-only SP(k) analyses. Similar to the BCEmu case, there is good agreement of the lensing

analysis with SP(k):wide and X-ray data at high masses $M_{500} \sim 1 \times 10^{14} M_\odot$. In agreement with the findings with BCEmu, we find that lensing prefers slightly lower baryon fractions for haloes $M_{500} \sim 10^{13} M_\odot$, with the SP(k): wide analysis lying lower by ~1.3$\sigma$. This highlights that the lensing data alone, when analysed with flexible modelling choices, prefers a lower baryon content in haloes to that measured by X-ray observations and that predicted by simulations.

## 7.2 Mass dependence of the gas mass fraction

The majority of baryonic mass in galaxy groups and clusters exists as diffuse gas and measurements of the fraction of the halo mass in gas, $M_{\text{gas}}/M_{500}$, are also sensitive to the matter power spectrum suppression (Schneider & Teyssier 2015; van Daalen et al. 2020; Aricò et al. 2023). The hot intracluster medium is X-ray luminous, hence studies of the X-ray emissivity of groups and clusters allow measurements of the gas mass to be derived.

In this section, we compare X-ray derived measurements of the $M_{\text{gas}}/M_{500}-M_{500}$ relation to those constrained by the BCEmu7 WL-only and joint WL + kSZ analyses in this work. As with the baryon fractions, we calculate the $M_{\text{gas}}/M_{500}-M_{500}$ relation at each step in the chain by integrating equation (9) to $r_{500}$. The right panel of Fig. 9 shows the mean and 68 per cent confidence levels on the gas fraction-halo mass relation for both the WL-only (pink) and WL + kSZ chains (purple). We compare to the $1\sigma$ region of the X-ray HSC-XXL constraint (Akino et al. 2022) in hatched black (also scaled to the mean cosmology obtained in the WL + kSZ analysis) and the baryon fractions measured in BAHAMAS (McCarthy et al. 2017), FLAMINGO (Schaye et al. 2023), MillenniumTNG (Pakmor et al. 2023), and FABLE (Henden et al. 2018, Bigwood et al. in preparation).

In consistency with the baryon fraction result, we find that the BCEmu7 WL-only and WL + kSZ constraints are in good agreement with each other and the X-ray data for massive groups
of $M_{500} \gtrsim 5 \times 10^{13} M_{\odot}$. For groups of lower masses we find that the joint WL + kSZ analysis significantly improves the constraining power on the gas fractions, with the 68 per cent confidence levels on $M_{gas}/M_{500}$ at $M_{500} \approx 10^{13} M_{\odot}$ shrinking by a factor of $\sim 2$ between the WL-only and WL + kSZ analyses. At $M_{500} \approx 10^{13} M_{\odot}$, the WL + kSZ data constrain a total baryon fraction that is lower than the X-ray data by $\sim 1.6\sigma$. The fiducial BAHAMAS simulation is in best agreement with the WL + kSZ constraint, with the remaining simulations predicting higher values of $M_{gas,500}/M_{500}$ at all halo masses. At $M_{500} \sim 4 \times 10^{13} M_{\odot}$, the fiducial FLAMINGO simulation predicts $\sim 2.3\sigma$ higher values, MillenniumTNG $\sim 7.2\sigma$ and FABLE $\sim 2.1\sigma$.

The baryon and gas fraction constraints by the WL + kSZ data imply that lower mass groups are expelling a greater amount of baryonic matter due to feedback than predicted by X-ray measurements. When combined with the findings of Section 6.2 that WL + kSZ predicts a greater mean matter power spectrum than all of the simulations, we build a consistent picture: that X-ray observations constrain a weaker feedback scenario than that preferred by WL + kSZ, subject to the uncertainties we have discussed.

## 8 SUMMARY AND CONCLUSIONS

WL measurements at small angular scales are statistically powerful. Not only do they offer a window to test $\Lambda$CDM in the non-linear regime, but also to constrain baryonic feedback effects and benchmark hydrodynamical simulations. Extracting accurate cosmological information relies on accurate modelling of physical processes associated with baryons which can re-distribute matter on small scales. Otherwise, unmodelled baryonic effects can bias the cosmological constraints from lensing analyses (Amon & Efstathiou 2022; Preston et al. 2023) or severely limit their precision (Amon et al. 2022; DES & KiDS Collaborations 2023). The aim of this work has been twofold. First, we test four baryon model approaches, comparing their ability to constrain cosmological parameters and the suppression of the matter power spectrum. Secondly, we perform a joint analysis of DES cosmic shear with ACT kSZ measurements. We demonstrate that a combined WL and kSZ analysis provides an exciting opportunity to not only improve constraints on cosmological parameters, but also on astrophysical effects. The main results of this study are:

(i) We perform the most extensive mock baryon model comparison for cosmic shear to date, in order to assess the robustness of the models to recover the underlying cosmology in a DES Y3-like analysis. We consider three different mock 'feedback scenarios' to test four baryon strategies; a halo model approach, a simulation-based emulator, an analytical $N$-body simulation model, and the approach of discarding small angular scales. In general, using restrictive modelling choices which do not capture the input matter power spectrum suppression can underestimate the recovery of $S_8$. Given the spread in the true suppression of the matter power spectrum as predicted by simulations and the lack of observational constraints on this quantity, model flexibility is crucial to ensure accurate cosmology, but it comes at a cost: the uncertainty on $S_8$ can degrade by up to a factor of $\sim 2$ with different model choices.

(ii) We analyse the DES Y3 cosmic shear data with the three baryon feedback models, with the scale cut approach, and with no feedback model. We find that each baryon model provides a good fit to the DES Y3 data, but the measured value of $S_8$ varies by $\sim 0.5\sigma - 2\sigma$ and the errorbar by a factor of $\sim 1.5$.

(iii) We output the posterior for the suppression of the matter power spectrum using the three models, with and without their informative priors. For all three of the models, without their informative priors, the mean suppression of the matter power spectrum constrained using DES Y3 WL is more extreme than the prediction from the hydrodynamical simulations considered.

(iv) We jointly analyse the DES Y3 cosmic shear with kSZ measurements from ACT DR5 and CMASS, in order to constrain the BCAemu7 model parameters. We find a slightly higher value of $S_8 = 0.823_{-0.020}^{+0.019}$, compared to the value attained by the WL-only analysis. If instead we analyse the DES Y3 cosmic shear jointly with X-ray baryon fraction constraints using SP(k), we find a lower value of $S_8 = 0.806_{-0.013}^{+0.015}$.

(v) The kSZ significantly improves the constraint on the suppression of the matter power spectrum from WL. The joint WL+kSZ predicts a more extreme suppression of the matter power spectrum than the WL scenario, with a mean constraint predicting a greater suppression than the fiducial BAHAMAS, MillenniumTNG, fiducial FLAMINGO, and FABLE simulations. At $k \gtrsim 1\hbar\,\mathrm{Mpc}^{-1}$, only SIMBA falls within the $1\sigma$ bound.

(vi) We constrain the baryon fraction-halo mass and the gas fraction-halo mass relations using WL + kSZ data. Both the baryon fraction and gas fraction is consistent with that from Akino et al. (2022) X-ray data within the 68 per cent confidence level for groups $M_{500} \approx 10^{14} M_{\odot}$. However for lower mass groups $M_{500} \approx 10^{13} M_{\odot}$ the baryon fraction lies $\sim 1.4\sigma$ lower than the X-ray data and the gas fraction lies $\sim 1.6\sigma$ lower.

(vii) Our constraints on the matter power spectrum suppression, baryon fractions and gas fractions all point towards a tension between the feedback of groups and clusters predicted by WL + kSZ and X-rays, or X-ray-calibrated models.

The next generation of shear surveys, such as Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST Dark Energy Science Collaboration 2012), Euclid (Laureijs et al. 2011), and the Nancy Grace Roman Space Telescope (Spergel et al. 2015) will deliver unprecedented statistical power to test the cosmological model in the non-linear regime ($k > 0.1\hbar\,\mathrm{Mpc}^{-1}$) using WL (Foreman, Becker & Wechsler 2016). This work highlights the importance of determining an accurate model for baryonic effects, which is flexible enough to not bias cosmological constraints. It will be crucial to test models of feedback with thorough mock analyses. In order for the model uncertainty in WL analyses to not limit the statistical power of the survey, it is crucial to either incorporate external probes of the gas content to constrain the additional baryon parameters, or to reduce the uncertainty in simulation-based priors. The latter requires a consistent picture for baryonic feedback effects, supported by a range of observations.

We demonstrate that joint analyses of gas measurements with WL data not only improve cosmological constraints, but provide valuable constraints on astrophysical feedback models and benchmark hydrodynamical simulations. We find a consistent picture that could imply that the WL and kSZ data is in tension with the X-ray measurements, and as a result, the predictions from simulations calibrated to X-ray data. This is not particularly surprising as X-ray measurements are generally sensitive to the hot gas content in the inner regions of clusters, compared to the outer regions and lower mass haloes that kSZ measurements probe.

Looking ahead, kSZ measurements as a function of mass and redshift will provide a handle for improved baryonic feedback models. Spectroscopic galaxy surveys, such as the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016) and the
Prime Focus Spectrograph (Sugai et al. 2012), will greatly increase the sample size of galaxy catalogues, in combination with the state-of-the-art CMB observations, for example, from Simons Observatory (Ade et al. 2019).

## ACKNOWLEDGEMENTS

LB is supported by a Science and Technology Facilities Council studentship. AA was supported by a Kavli Fellowship for part of the duration of this work. JS and IGM were supported by an European Research Council (ERC) Consolidator Grant under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 769130).

The author's contributions are as follows. LB and AA performed the analysis and wrote the bulk of the manuscript. The additional core authors of AA, JS, and IGM contributed BCEmu and SP(k) code and continuous insight. CP and DS assisted in the early stages, including the mock analysis. DSi provided useful discussions and detailed comments on the manuscript. ES, SF, and NB provided the ACT kSZ data set. The paper has gone through internal review by the DES collaboration, with A. Chen and S. Dodelson serving as excellent internal reviewers. The remaining authors have made contributions to the DES instrumentation, data collection, data processing, and calibration, in addition to other code and validation tests.

This project has used public archival data from the DES. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, NSF NOIRLab, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at NSF Cerro Tololo Inter-American Observatory at NSF NOIRLab (NOIRLab Prop. ID 2012B-0001; PI: J. Frieman), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under grant nos AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MICINN under grants ESP2017-89838, PGC2018-094773, PGC2018-102021, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) do e-Universo (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under contract no. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.

## DATA AVAILABILITY

The data underlying this article may be shared on reasonable request to the corresponding authors.

## REFERENCES

Adamek J. et al., 2023, J. Cosmol. Astropart. Phys., 2023, 035

Ade P. et al., 2019, J. Cosmol. Astropart. Phys., 2019, 056

Akino D. et al., 2022, PASJ, 74, 175

Amodeo S. et al., 2021, Phys. Rev. D, 103, 063514

Amon A., Efstathiou G., 2022, MNRAS, 516, 5355

Amon A. et al., 2022, Phys. Rev. D, 105, 023514

Amon A. et al., 2023, MNRAS, 518, 477

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ò G., Angulo R. E., Zennaro M., Contreras S., Chen A., Hernández-Monteagudo C., 2023, A&A, 678, A109

Asgari M. et al., 2021, A&A, 645, A104

Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143

Bennett C. L. et al., 2003, ApJS, 148, 1

Blazek J. A., MacCrann N., Troxel M. A., Fang X., 2019, Phys. Rev. D, 100, 103506

Bridle S., King L., 2007, New J. Phys., 9, 444

Chen A. et al., 2023, MNRAS, 518, 5340

Chisari N. E. et al., 2019, Open J. Astrophys., 2, 4

DES, KiDS Collaborations, 2023, Open J. Astrophys., 6, 53

DES Collaboration, 2022, Phys. Rev. D, 105, 023520

DESI Collaboration, 2016, preprint (arXiv:1611.00036)

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827

Debackere S. N. B., Schaye J., Hoekstra H., 2020, MNRAS, 492, 2285

Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948

Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359

Eckert D. et al., 2016, A&A, 592, A12

Efstathiou G., Gratton S., 2021, Open J. Astrophys., 4, 8

Eisenstein D. J. et al., 2005, ApJ, 633, 560

Euclid Collaboration: Knabenhans M. et al., 2019, MNRAS, 484, 5509

Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601

Ferreira T., Alonso D., Garcia-Garcia C., Chisari N. E., 2024, Phys. Rev. Lett., 133, 051001

Foreman S., Becker M. R., Wechsler R. H., 2016, MNRAS, 463, 3326

García-García C., Zennaro M., Aricò G., Alonso D., Angulo R. E., 2024, J. Cosmol. Astropart. Phys., 2024, 024

Gatti M., et al., 2022, Phys. Rev. D, 106, 083509
Giri S. K., Schneider A., 2021, J. Cosmol. Astropart. Phys., 2021, 046

Grandis S., Aricò G., Schneider A., Linke L., 2024, MNRAS, 528, 4379

Hadzhiyska B., Ferraro S., Ried Guachalla B., Schaan E., 2024, Phys. Rev. D, 109, 103534

Hand N. et al., 2012, Phys. Rev. Lett., 109, 041101

Handley W. J., Hobson M. P., Lasenby A. N., 2015, MNRAS, 453, 4384

Henden N. A., Puchwein E., Shen S., Sijacki D., 2018, MNRAS, 479, 5385

Hernández-Aguayo C. et al., 2023, MNRAS, 524, 2556

Hirata C. M., Seljak U., 2004, Phys. Rev. D, 70, 063526

Hoekstra H., Herbonnet R., Muzzin A., Babul A., Mahdavi A., Viola M., Cacciato M., 2015, MNRAS, 449, 685

Huang H.-J., Eifler T., Mandelbaum R., Dodelson S., 2019, MNRAS, 488, 1652

Knabenhans M. et al., 2021, MNRAS, 505, 2840

Krause E. et al., 2021, preprint (arXiv:2105.13548)

Kravtsov A. V., Vikhlinin A. A., Meshcheryakov A. V., 2018, Astron. Lett., 44, 8

Kugel R. et al., 2023, MNRAS, 526, 6103

LSST Dark Energy Science Collaboration, 2012, preprint (arXiv:1211.0310)

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Le Brun A. M. C., McCarthy I. G., Schaye J., Ponman T. J., 2014, MNRAS, 441, 1270

Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473

Li X. et al., 2023, Phys. Rev. D, 108, 123518

Limber D. N., 1953, ApJ, 117, 134

LoVerde M., Afshordi N., 2008, Phys. Rev. D, 78, 123506

MacCrann N. et al., 2022, MNRAS, 509, 3371

Madhavacheril M. S. et al., 2024, ApJ, 962, 113

McCarthy I. G., Schaye J., Bower R. G., Ponman T. J., Booth C. M., Dalla Vecchia C., Springel V., 2011, MNRAS, 412, 1965

McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465, 2936

McCarthy I. G. et al., 2023, MNRAS, 526, 5494

Mead A. J., Brieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401

Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121

Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822

Mulroy S. L. et al., 2019, MNRAS, 484, 60

Myles J., et al., 2021, MNRAS, 505, 4249

Naess S. et al., 2020, J. Cosmol. Astropart. Phys., 2020, 046

Nelson D. et al., 2019, Comput. Astrophys. Cosmol., 6, 2

Pakmor R. et al., 2023, MNRAS, 524, 2539

Pandey S. et al., 2023, MNRAS, 525, 1779

Patrignani C., 2016, Chin. Phys. C, 40, 100001

Planck Collaboration XXXVII, 2016, A&A, 586, A140

Planck Collaboration VI, 2020a, A&A, 641, A6

Planck Collaboration VIII, 2020b, A&A, 641, A8

Preston C., Amon A., Efstathiou G., 2023, MNRAS, 525, 5554

Rasia E. et al., 2006, MNRAS, 369, 2013

Ried Guachalla B., Schaan E., Hadzhiyska B., Ferraro S., 2024, Phys. Rev. D, 109, 103533

Salcido J., McCarthy I. G., Kwan J., Upadhye A., Font A. S., 2023, MNRAS, 523, 2247

Schaan E. et al., 2021, Phys. Rev. D, 103, 063513

Schaye J. et al., 2010, MNRAS, 402, 1536

Schaye J. et al., 2023, MNRAS, 526, 4978

Schneider A., Teyssier R., 2015, J. Cosmol. Astropart. Phys., 2015, 049

Schneider A., Teyssier R., Stadel J., Chisari N. E., Le Brun A. M. C., Amara A., Refregier A., 2019, J. Cosmol. Astropart. Phys., 2019, 020

Schneider A., Giri S. K., Amodeo S., Refregier A., 2022, MNRAS, 514, 3802

Secco L. F., Samuroff S. et al., 2022, Phys. Rev. D, 105, 023515

Sevilla-Noarbe I. et al., 2021, ApJS, 254, 24

Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877

Soergel B. et al., 2016, MNRAS, 461, 3172

Sonnenfeld A., Wang W., Bahcall N., 2019, A&A, 622, A30

Spergel D. et al., 2015, preprint (arXiv:1503.03757)

Springel V. et al., 2018, MNRAS, 475, 676

Stebbins A., 1996, preprint (arXiv:astro-ph/9609149)

Sugai H. et al., 2012, Proc. SPIE, 8446, 84460Y

Sunyaev R. A., Zeldovich Y. B., 1972, Comment. Astrophys. Space Phys., 4, 173

Sunyaev R. A., Zeldovich Y. B., 1980, MNRAS, 190, 413

Terasawa R. et al., 2024, preprint (arXiv:2403.20323)

Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2011, MNRAS, 414, 195

To C.-H., Pandey S., Krause E., Dalal N., Anbajagane D., Weinberg D. H., 2024, J. Cosmol. Astropart. Phys., 2024, 037

Tröster T. et al., 2020, A&A, 633, L10

van Daalen M. P., Schaye J., Booth C. M., Dalla Vecchia C., 2011, MNRAS, 415, 3649

van Daalen M. P., McCarthy I. G., Schaye J., 2020, MNRAS, 491, 2424

White M. et al., 2011, ApJ, 728, 126

Zuntz J. et al., 2015, Astron. Comput., 12, 45

## APPENDIX A: MOCK ANALYSIS

In this section, we provide more detail on the mock analysis presented in Section 4.2. We create the synthetic cosmic shear data using the best-fitting cosmological parameters attained from the DES Y3 joint lensing and clustering analysis (DES Collaboration 2022).<fnref n="12" /> As described in Section 4.2, the mocks were created with HM20 as the dark-matter-only non-linear correction to the matter power spectrum, unless otherwise stated. We create data vectors with three variants for the impact on the non-linear power spectrum of baryonic feedback; we use the power suppression predicted by BAHAMAS 8.0, cosmo-OWLS 8.5, and $A_{\text{mod}} = 0.858$. In this section, we present a more detailed summary of the dark matter and baryonic feedback model validation tests ran using synthetic data vectors. We refer to Table A1 which presents the MULTINEST mean values of $S_8$ attained for the full mock suite, reporting the relative shifts from the input cosmology. The results are also summarized in Fig. A2, where we additionally plot the mean and best-fitting $\Omega_m$ and $\sigma_8$ for each mock test.

We validate the choice to model the dark-matter-only power spectrum with HM20 using dark-matter-only mocks created with HM20 and EE2. Here, we summarize the findings of mock tests 1–4.

(i) This test analyses a mock with the same model choices used to create it and is useful for identifying projection effects. We find these to be present for $\Omega_m$ and $\sigma_8$, which are under- and overestimated, respectively, by ~0.5σ.

(ii) When modelling dark matter with HM20 in a HM20 generated mock, allowing the neutrino mass to vary with the prior $\sum m_\nu$: [0.06,6] eV improves the recovery of $\Omega_m$ and $\sigma_8$. It however decreases the accuracy in the recovery of $S_8$, resulting in an underestimation of 0.5σ.

(iii) Modelling the non-linear dark-matter-only power spectrum with HM20 in an EE2 generated mock overestimates $S_8$ by ~0.4σ, which is reduced marginally to ~0.2σ when scale cuts are applied. We note that DES & KiDS Collaborations (2023) find a smaller bias in the mean of ~0.1σ when scale cuts are used. In this mock analysis, there are two setup differences that could explain this: we

<footnote n="12">¹²$\Omega_m = 0.3380$, $\Omega_b = 0.0450$, $10^{-9}A_s = 1.8418$, $h = 0.615$, $n_s = 0.949$, $S_8 = 0.7805$, and $\sigma_8 = 0.7353$.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure A1.** Marginalized posteriors for $S_8$ and the baryon model parameter using each mock for analyses using the HM20 (upper), BCEmu (middle), and SP(k) (lower) models. The mocks are labelled in the format 'model used to create the mock: model/analysis choices used to analyse the mock'. The inner and outer contours show the 68 per cent and 95 per cent confidence levels, respectively and the dashed line indicates the input cosmology.</figcaption></figure>

use a higher value of $\Omega_m$ to create the mock and we marginalize over $A_s$, instead of $S_8$.

Next we validate each of our four baryon models (scale cuts, HM20 with free $\Theta_{AGN}$, BCEmu, and SP(k)) and their respective prior and analysis variants. We test the modelling choices on three 'baryonic feedback scenarios' of increasing 'extremity' in terms of their impact on the matter power spectrum suppression: BAHAMAS 8.0, cosmo-OWLS 8.5, and $A_{\text{mod}} = 0.858$. Here, we summarize the findings of the mock tests 5–29 and Fig. A1. Note that in the figure, mocks are labelled in the format 'model used to create the mock: model/analysis choices used to analyse the mock'. For example, 'BAHAMAS: HM20 BAHAMAS' denotes a mock created with BAHAMAS-like baryon feedback scenario and analysed with HM20 with the BAHAMAS $\Theta_{AGN} = 7.6 - 8.0$ prior.

(i) Optimized scale cuts underestimates $S_8$ by $0.5\sigma - 1.2\sigma$ in the three feedback scenarios we test, with the tension worsening with more extreme feedback. As noted in Section 4.2, this is as expected given our choice of feedback scenarios used to build the synthetic data. The OWLS-AGN scenario, which was used to decide which angular scales of the measurement would be discarded in the DES Y3 analysis predicts a less extreme feedback scenario than the three scenarios used here.

(ii) The left panel of Fig. A1 compares prior choices of the HM20 model with $\Theta_{AGN}$ as a free parameter. We find that HM20 with the BAHAMAS calibrated prior choice of $\Theta_{AGN} = 7.6 - 8.0$ underestimates $S_8$ by $\sim 0.7\sigma - 1.2\sigma$, with the tension worsening with a more extreme feedback scenario. Using the wide prior choice of $\Theta_{AGN} = 7.3 - 9.0$ overestimates $S_8$ in BAHAMAS 8.0 and cosmo-OWLS 8.5 by $\sim 0.5\sigma$, but still underestimates $S_8$ in an $A_{\text{mod}} = 0.858$-like feedback scenario by $\sim 0.2\sigma$.

(iii) The right panel of Fig. A1 shows that SP(k) with the X-ray informed prior demonstrates the best recovery of $S_8$ in a BAHAMAS 8.0-like feedback scenario of all the baryonic feedback modelling choices. In more extreme feedback scenarios, it however underestimates $S_8$ by $\sim 0.3\sigma - 1.6\sigma$. Using the wide prior choice improves the recovery of $S_8$ to be within $\sim 0.4\sigma$ for the three feedback scenarios.

(iv) The central panel of Fig. A1 demonstrates that BCEmu7 overestimates $S_8$ in both BAHAMAS 8.0 and cosmo-OWLS 8.5 baryonic feedback scenarios by $\sim 0.4\sigma$, but underestimates $S_8$ in an $A_{\text{mod}} = 0.858$-like feedback scenario by $\sim 0.2\sigma$. BCEmu1 underestimates $S_8$ by $0.3\sigma - 1.2\sigma$, again with the tension worsening when testing more extreme feedback scenarios.

(v) In the mock tests which overestimate $S_8$, freeing the neutrino mass with the prior $\sum m_\nu$: [0.06,6] eV appears to improve recovery of the input cosmology. The shifts we see in $S_8$ are however consistent with the $\sim 0.2\sigma$ shift to lower values we find when testing the dark-matter-only modelling of HM20 with and without free neutrinos, that is, tests 1 and 2.

Ultimately we find that choosing the restrictive analysis variant of a model tends to result in an underestimation of $S_8$. This motivates the use of more flexible modelling choices in cosmic shear analyses, since these tend to result in improved recovery of $S_8$, however this is at the consequence of a loss in precision.
<table><caption>**Table A1.** Summary of the mock tests we performed to validate the modelling of the non-linear matter power spectrum. Mock tests 1–3 validate the non-linear correction to the matter power spectrum due to dark matter only. The remaining tests validate the baryon feedback mitigation strategies we test in this work; scale cuts (4–7), HM20 (8–14), SP(*k*) (15–20), and BCEmu (21–29). We test analysing the mocks with various prior and analysis choices for each of the four baryon model approaches. If not otherwise specified, the mocks were created with HM20 as the dark-matter-only non-linear correction to the matter power spectrum, and three ‘baryonic feedback scenarios’; BAHAMAS, cosmo-OWLS, and $A_{\text{mod}} = 0.858$. The ‘mock’ column therefore states the model/baryonic feedback scenario used to create the mock and the ‘model’ column labels the model/analysis choices used to analyse the mock’. For example, a mock of ‘EuclidEmu’ and model of HM20-DM’ denotes a mock created with Euclid Emulator for the non-linear correction to the matter power spectrum, and analysed with HM20. We report the $S_8$ constraints and 68 per cent confidence level using the mean-marginal approach $S_8$, with $\Delta S_8$ quantifying the offset from the true $S_8 = 0.7805$, as a fraction of the $1\sigma$ error.</caption><thead><tr><th>No.</th><th>Mock</th><th>Model</th><th>S<sub>8</sub></th><th>ΔS<sub>8</sub></th></tr></thead><tbody><tr><td>1</td><td>HM20-DM</td><td>HM20-DM</td><td>0.781<sup>+0.012</sup><sub>0.011</sub></td><td>0.044σ</td></tr><tr><td>2</td><td>HM20-DM</td><td>HM20-DM ∑ m<sub>ν</sub>: [0.06,6] eV</td><td>0.775<sup>+0.012</sup><sub>0.011</sub></td><td>-0.461σ</td></tr><tr><td>3</td><td>EuclidEmu-DM</td><td>HM20-DM</td><td>0.785<sup>+0.012</sup><sub>0.011</sub></td><td>0.403σ</td></tr><tr><td>4</td><td>EuclidEmu-DM</td><td>HM20-DM + Scale cuts</td><td>0.784<sup>+0.017</sup><sub>0.017</sub></td><td>0.239σ</td></tr><tr><td>5</td><td>BAHAMAS 8.0</td><td>Scale cuts</td><td>0.771<sup>+0.015</sup><sub>0.018</sub></td><td>-0.566σ</td></tr><tr><td>6</td><td>cosmo-OWLS 8.5</td><td>Scale cuts</td><td>0.763<sup>+0.013</sup><sub>0.017</sub></td><td>-1.191σ</td></tr><tr><td>7</td><td>A<sub>mod</sub> = 0.858</td><td>Scale cuts</td><td>0.762<sup>+0.016</sup><sub>0.018</sub></td><td>-1.114σ</td></tr><tr><td>8</td><td>BAHAMAS 8.0</td><td>HM20 BAHAMAS</td><td>0.772<sup>+0.013</sup><sub>0.012</sub></td><td>-0.675σ</td></tr><tr><td>9</td><td>cosmo-OWLS 8.5</td><td>HM20 BAHAMAS</td><td>0.765<sup>+0.013</sup><sub>0.014</sub></td><td>-1.182σ</td></tr><tr><td>10</td><td>A<sub>mod</sub> = 0.858</td><td>HM20 BAHAMAS</td><td>0.765<sup>+0.012</sup><sub>0.013</sub></td><td>-1.270σ</td></tr><tr><td>11</td><td>BAHAMAS 8.0</td><td>HM20 wide</td><td>0.792<sup>+0.021</sup><sub>0.023</sub></td><td>0.509σ</td></tr><tr><td>12</td><td>BAHAMAS 8.0</td><td>HM20 wide ∑ m<sub>ν</sub>: [0.06,6] eV</td><td>0.788<sup>+0.017</sup><sub>0.023</sub></td><td>0.358σ</td></tr><tr><td>13</td><td>cosmo-OWLS 8.5</td><td>HM20 wide</td><td>0.792<sup>+0.021</sup><sub>0.022</sub></td><td>0.537σ</td></tr><tr><td>14</td><td>A<sub>mod</sub> = 0.858</td><td>HM20 wide</td><td>0.777<sup>+0.019</sup><sub>0.024</sub></td><td>-0.180σ</td></tr><tr><td>15</td><td>BAHAMAS 8.0</td><td>SP(*k*): X-ray</td><td>0.781<sup>+0.017</sup><sub>0.018</sub></td><td>0.017σ</td></tr><tr><td>16</td><td>cosmo-OWLS 8.5</td><td>SP(*k*): X-ray</td><td>0.775<sup>+0.016</sup><sub>0.018</sub></td><td>-0.347σ</td></tr><tr><td>17</td><td>A<sub>mod</sub> = 0.858</td><td>SP(*k*): X-ray</td><td>0.761<sup>+0.012</sup><sub>0.012</sub></td><td>-1.643σ</td></tr><tr><td>18</td><td>BAHAMAS 8.0</td><td>SP(*k*) wide</td><td>0.780<sup>+0.015</sup><sub>0.016</sub></td><td>-0.010σ</td></tr><tr><td>19</td><td>cosmo-OWLS 8.5</td><td>SP(*k*) wide</td><td>0.774<sup>+0.015</sup><sub>0.017</sub></td><td>-0.424σ</td></tr><tr><td>20</td><td>A<sub>mod</sub> = 0.858</td><td>SP(*k*) wide</td><td>0.773<sup>+0.015</sup><sub>0.016</sub></td><td>-0.479σ</td></tr><tr><td>21</td><td>BAHAMAS 8.0</td><td>BCEmu1</td><td>0.776<sup>+0.017</sup><sub>0.018</sub></td><td>-0.248σ</td></tr><tr><td>22</td><td>cosmo-OWLS 8.5</td><td>BCEmu1</td><td>0.773<sup>+0.015</sup><sub>0.017</sub></td><td>-0.489σ</td></tr><tr><td>23</td><td>A<sub>mod</sub> = 0.858</td><td>BCEmu1</td><td>0.759<sup>+0.016</sup><sub>0.019</sub></td><td>-1.239σ</td></tr><tr><td>24</td><td>BAHAMAS 8.0</td><td>BCEmu7</td><td>0.792<sup>+0.020</sup><sub>0.025</sub></td><td>0.510σ</td></tr><tr><td>25</td><td>BAHAMAS 8.0</td><td>BCEmu7 ∑ m<sub>ν</sub>: [0.06,6] eV</td><td>0.785<sup>+0.020</sup><sub>0.025</sub></td><td>0.202σ</td></tr><tr><td>26</td><td>BAHAMAS 8.0 + EuclidEmu-DM</td><td>BCEmu7</td><td>0.791<sup>+0.020</sup><sub>0.025</sub></td><td>0.466σ</td></tr><tr><td>27</td><td>BAHAMAS 8.0 + EuclidEmu-DM</td><td>BCEmu7 ∑ m<sub>ν</sub>: [0.06,6] eV</td><td>0.787<sup>+0.019</sup><sub>0.027</sub></td><td>0.287σ</td></tr><tr><td>28</td><td>cosmo-OWLS 8.5</td><td>BCEmu7</td><td>0.791<sup>+0.022</sup><sub>0.026</sub></td><td>0.439σ</td></tr><tr><td>29</td><td>A<sub>mod</sub> = 0.858</td><td>BCEmu7</td><td>0.777<sup>+0.018</sup><sub>0.024</sub></td><td>-0.165σ</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Figure A2.** Summary of the 1D marginalized constraints on $S_8$, $\Omega_m$, and $\sigma_8$ attained by our mock analysis with respect to the input cosmology. We plot the mean and 68 per cent confidence levels as listed in Table A1, with crosses showing the MULTINEST best-fitting result. The mocks were created with DES-Y3 covariance, with the input cosmology indicated by the vertical yellow line. The top panel shows the results of validating the non-linear correction to the matter power spectrum due to dark matter only. The mocks are labelled in the format 'model used to create the mock: model/analysis choices used to analyse the mock'. For example, 'EuclidEmu: HM20-DM' denotes a mock created with Euclid Emulator for the non-linear correction to the matter power spectrum, and analysed with HM20. The remaining panels validate the four methods of mitigating for baryonic feedback we test in this work; optimized scale cuts, HM20, BCEmu, and SP(k). If not otherwise specified, the mocks were created with HM20 as the dark-matter-only non-linear correction to the matter power spectrum, and three 'baryonic feedback scenarios': BAHAMAS, cosmo-OWLS, and $A_{\text{mod}} = 0.858$. We test analysing the mocks with various prior and analysis choices for each of the four baryon model approaches. Each mock was analysed with HM20 as the dark-matter-only non-linear correction to the matter power spectrum.</figcaption></figure>
## APPENDIX B: BARYON FEEDBACK
CONSTRAINTS

### B1 SP(k)

The mean and 68 per cent confidence levels on the SP(k) parameters constrained by the WL-only SP(k) analyses with the wide and X-ray informed prior choices are reported in Table B1. The corresponding marginalized posteriors for these analyses are shown in the upper panel of Fig. B1. Since SP(k) directly maps between the baryon fraction in groups and clusters, and the matter power spectrum suppression, we also show the constraints on the $f_b/\Omega_b/\Omega_m - M_{500}$ relation predicted when analysing the DES cosmic shear with the two SP(k) analysis variants in the lower panel of Fig. B1.

### B2 Baryonification

The mean and 68 per cent confidence levels on the bfc parameters constrained by the WL-only BCEmu1, BCEmu3 and BCEmu7 analyses, in addition to the WL + kSZ BCEmu7 joint analysis (for which the average halo mass of the kSZ sample $M_{h,200}$ is also reported) are given in Table B2. The corresponding marginalized posteriors for these analyses are shown in Fig. B2.

### B3 Prior choice for the halo mass of the kSZ sample

Given the significant scatter in the literature, we choose to include an additional model parameter in the analysis, $M_{h,200}$, corresponding to the mean $M_{200}$ of the CMASS sample, with a prior range provided in Table 1. For example, the stacked stellar mass-halo mass relation of Sonnenfeld, Wang & Bahcall (2019) of CMASS galaxies derived from HSC galaxy-galaxy lensing measurements, when combined with the stellar mass distribution of CMASS, implies a mean $M_{200} \approx 0.5 \times 10^{13} M_\odot$, whereas the abundance matching methods

<table><caption>**Table B1.** The constraints on the SP(k) parameters attained by the WL-only SP(k) analyses with both the wide and X-ray informed prior choices. We report the mean and 68 per cent confidence levels.</caption><thead><tr><th>Parameter</th><th>Wide prior</th><th>X-ray prior</th></tr></thead><tbody><tr><td>α</td><td>3.830<sup>+0.409</sup><sub>-0.363</sub></td><td>4.153<sup>+0.061</sup><sub>-0.058</sub></td></tr><tr><td>β</td><td>1.174<sup>+0.104</sup><sub>-0.158</sub></td><td>1.206<sup>+0.043</sup><sub>-0.042</sub></td></tr><tr><td>γ</td><td>0.519<sup>+0.265</sup><sub>-0.157</sub></td><td>0.388<sup>+0.077</sup><sub>-0.076</sub></td></tr></tbody></table>

<table><caption>**Table B2.** The constraints on the bfc parameters attained by the WL-only BCEmu1, BCEmu3 and BCEmu7 analyses, as well as the BCEmu7 WL + kSZ BCEmu7 analysis. For the WL + kSZ analysis we also report the constraint on the average halo mass of the kSZ sample ($10^{13} M_\odot$). We report the mean and 68 per cent confidence levels.</caption><thead><tr><th>Parameters</th><th>BCEmu1</th><th>BCEmu3</th><th>BCEmu7</th><th>WL + kSZ</th></tr></thead><tbody><tr><td>log<sub>10</sub> M<sub>c</sub></td><td>13.42<sup>+0.60</sup><sub>-0.54</sub></td><td>13.13<sup>+0.79</sup><sub>-0.63</sub></td><td>13.06<sup>+1.01</sup><sub>-0.87</sub></td><td>13.22<sup>+0.42</sup><sub>-0.29</sub></td></tr><tr><td>θ<sub>ej</sub></td><td>—</td><td>5.31<sup>+2.17</sup><sub>-1.30</sub></td><td>5.14<sup>+1.75</sup><sub>-1.74</sub></td><td>5.15<sup>+1.46</sup><sub>-1.48</sub></td></tr><tr><td>η<sub>δ</sub></td><td>—</td><td>0.22<sup>+0.09</sup><sub>-0.12</sub></td><td>0.23<sup>+0.10</sup><sub>-0.10</sub></td><td>0.22<sup>+0.10</sup><sub>-0.11</sub></td></tr><tr><td>μ</td><td>—</td><td>—</td><td>1.00<sup>+0.50</sup><sub>-0.61</sub></td><td>0.68<sup>+0.25</sup><sub>-0.58</sub></td></tr><tr><td>γ</td><td>—</td><td>—</td><td>2.57<sup>+1.00</sup><sub>-0.71</sub></td><td>2.66<sup>+0.98</sup><sub>-0.67</sub></td></tr><tr><td>δ</td><td>—</td><td>—</td><td>6.83<sup>+2.09</sup><sub>-2.43</sub></td><td>7.35<sup>+2.16</sup><sub>-2.10</sub></td></tr><tr><td>η</td><td>—</td><td>—</td><td>0.22<sup>+0.10</sup><sub>-0.11</sub></td><td>0.21<sup>+0.08</sup><sub>-0.13</sub></td></tr><tr><td>M<sub>h,200</sub></td><td>—</td><td>—</td><td>—</td><td>4.10<sup>+1.55</sup><sub>-1.64</sub></td></tr></tbody></table>

<figure><img src="image_12.png" /></figure>

SPk: wide

SPk: Xray

<figure><img src="image_13.png" /><figcaption>**Figure B1.** SP(k) marginalized posteriors when analysing the DES Y3 cosmic shear data with the wide prior on SP(k) (light green) and X-ray informed prior (dark green). We show the 68 per cent and 95 per cent confidence levels of the $S_8$ and the three model parameters (upper panel) and the mean total baryon fraction, $f_b/\Omega_b/\Omega_m$, and $1\sigma$ uncertainty as a function of halo mass, $M_{500}$ (lower). The X-ray prior is derived from HSC-XXL $1\sigma$ constraints (black hatched, Akino et al. 2022), in this plot scaled to the mean cosmology obtained from the SPk: X-ray analysis.</figcaption></figure>

of Kravtsov, Vikhlinin & Meshcheryakov (2018), Moster, Naab & White (2018), and Behroozi et al. (2019) imply mean $M_{200}$ values of $2.7 \times 10^{13}$, $4.4 \times 10^{13}$, and $6.6 \times 10^{13} M_\odot$, respectively. A halo occupation distribution (HOD)-based analysis of the clustering of BOSS CMASS galaxies by White et al. (2011) found a mean $M_{200}$ of $3.6 \times 10^{13} M_\odot$. Given this large study-to-study variance, we choose a flat prior range of $[0.5-7] \times 10^{13} M_\odot$ and marginalize over this parameter.

The combined WL + kSZ analysis constrains $M_{h,200} = 4.098^{+1.548}_{-1.639} \times 10^{13} M_\odot$, indicating we are prior constrained at the $2\sigma$ level. For comparison to the remainder of this work which
<figure><img src="image_14.png" /><figcaption>**Figure B2.** The marginalized posteriors for $S_8$, the seven bfc parameters, and the average halo mass of the kSZ sample $M_{h,200}$ attained by the WL-only BCEmu1 (red solid line), BCEmu3 (light purple dotted line), and BCEmu7 (pink dashed line) analyses, as well as the WL + kSZ BCEmu7 (dark purple line) analysis. The inner and outer contours show the 68 per cent and 95 per cent confidence levels, respectively. The grey dashed lines show the values the BCEmu parameters are fixed to in the case of using BCEmu1 and BCEmu3, and the mean halo mass of the CMASS sample reported by Schaan et al. (2021).</figcaption></figure>

generally quotes halo masses in $M_{500}$, this corresponds to $M_{h,500} = 3.01 \times 10^{13} M_{\odot}$, assuming an NFW profile and a concentration-mass relation from Dutton & Macciò (2014).

In this appendix, we explore the impact on our cosmology and baryon model parameter constraints of choosing a wider prior $M_{h,200}$: [0.8, 30] $\times 10^{13} M_{\odot}$. We also consider a fixed $M_{h,200}$ analysis using the mean mass of the CMASS sample as determined by Schaan et al. (2021) and used in Amodeo et al. (2021): $M_{h,200} = 3 \times 10^{13} M_{\odot}$. Fig. B3 shows the marginalized posteriors for $\Omega_m$, $S_8$, $\log_{10}(M_c)$, $\theta_{ej}$, and $M_{h,200}$ in a WL + kSZ analysis with the different prior choices on $M_{h,200}$. The halo mass, $M_{h,200}$, is correlated with

$\log_{10}(M_c)$ and in the case of a wider prior, both $M_{h,200}$ and $\log_{10}(M_c)$ prefer higher values. However, $S_8$ is relatively stable to the halo mass prior, shifting by $0.3\sigma$.

## B4 The relationship between baryon fraction and matter power spectrum suppression

It has been demonstrated that in hydrodynamical simulations, the mean baryon fraction in haloes of $M \sim 10^{14} M_{\odot}$ can be predictive for the suppression of the matter power spectrum due to baryonic feedback effects, robust to a number of feedback prescriptions (van
<figure><img src="image_15.png" /><figcaption>**Figure B3.** Upper: the stacked kSZ temperature profile at 98 GHz as a function of angular radius, $\theta$, centred on the group or cluster (bottom) when varying the mean halo mass of CMASS galaxy sample modelled, $M_{h,200}$, within the limits of the fiducial prior choice $M_{h,200}: [5 \times 10^{12}, 7 \times 10^{13}] M_{\odot}$. The ACT CMASS measurements at 98 GHz are shown as the black data points in the bottom panels and the model profiles are convolved with the f90 beam profile for comparison. Lower: the marginalized posteriors for $S_8$, $\log_{10}(M_c)$, and $M_{h,200}$ in a WL + kSZ analysis with different prior choices on $M_{h,200}$. We show the prior used in the fiducial analysis $M_{h,200}: [5 \times 10^{12}, 7 \times 10^{13}] M_{\odot}$ (purple), as well as a wide prior $M_{h,200}: [8 \times 10^{12}, 3 \times 10^{14}] M_{\odot}$ (light blue) and fixed $M_{h,200} = 3 \times 10^{13} M_{\odot}$ (dark blue, also shown as the dashed grey line). The inner and outer contours show the 68 per cent and 95 per cent confidence levels, respectively.</figcaption></figure>

Daalen et al. 2020). This relationship can be described by an empirical fitting function relating the mean baryon fraction, $f_b/\Omega_b/\Omega_m$, measured within $R_{500}$ for haloes of mass $M_{500} = 10^{14} M_{\odot}$ to the suppression of the matter power spectrum $P(k)/P_{\text{DMonly}}(k)$ at $k = 1h\,\text{Mpc}^{-1}$ (equation 5 in their work). The best-fitting relation, shown as the solid grey line in Fig. B4, was fit to the cosmo-OWLS (Le Brun

<figure><img src="image_16.png" /><figcaption>**Figure B4.** The relation between the matter power spectrum suppression $P(k)/P_{\text{DMonly}}(k)$ at $k = 1h\,\text{Mpc}^{-1}$ and the mean baryon fraction, $f_b/\Omega_b/\Omega_m$, measured within $R_{500}$ for haloes of mass $M_{500} = 10^{14} M_{\odot}$. We plot the empirical best-fitting relation of van Daalen et al. (2020) as the grey solid line, with the 1 per cent accuracy on the relation's ability to predict $P(k)/P_{\text{DMonly}}(k)$ shown as the shaded grey region. We plot the constraints obtained in our analysis in the $P(k)/P_{\text{DMonly}}(k)-f_b/\Omega_b/\Omega_m$ plane. We show the WL-only (pink starred) and WL + kSZ (pruple starred) analyses with BCEmu7, and the WL-only results with SPk, showing both the wide (light green starred) and X-ray (dark green starred) prior choices. We also plot the simulations FLAMINGO (Schaye et al. 2023); BAHAMAS (McCarthy et al. 2017); SIMBA (Davé et al. 2019); MillenniumTNG (Pakmor et al. 2023); and FABLE (Henden et al. 2018, Bigwood et al. in preparation) as the black data points.</figcaption></figure>

et al. 2014) and BAHAMAS (McCarthy et al. 2017) simulations and is accurate to 1 per cent for the simulations they test, shown as the grey shaded region in Fig. B4. In this section, we discuss where constraints of our analysis lie in the $P(k)/P_{\text{DMonly}}(k)-f_b/\Omega_b/\Omega_m$ plane with respect to the van Daalen et al. (2020) relation. In Fig. B4, we plot the result for the WL-only and WL + kSZ analyses with BCEmu7 baryon modelling, as well as the WL-only results with SP(k)'s wide and X-ray informed prior choices. For comparison, we also show the simulations discussed throughout this work: FLAMINGO (Schaye et al. 2023); BAHAMAS (McCarthy et al. 2017); SIMBA (Davé et al. 2019); MillenniumTNG (Pakmor et al. 2023); and FABLE (Henden et al. 2018, Bigwood et al. in preparation).

This van Daalen et al. (2020) relationship is at the core of the the SP(k) model (Salcido et al. 2023). Given that the ANTILLES simulations follow this relationship, the mean baryon fraction is used to calibrate their emulator. Naturally, both of our WL-only analyses with SP(k) fall on the van Daalen et al. (2020) relation to within 1 per cent. However, the BCEmu model does not enforce the van Daalen et al. (2020) relationship between the baryon fraction and matter power suppression, and allows for greater flexibility in the impact of feedback on the matter distribution, including scenarios which, according to the simulations, could be deemed unphysical. Both the WL-only and WL + kSZ analyses with BCEmu give a mean constraint below the van Daalen et al. (2020) relation. That is, they allow for a more extreme matter power spectrum suppression for their predicted mean baryon fraction predicted. In the future, it will be important to test this relationship with a range of observations, including those probe lower halo masses (see also Pandey et al. 2023).
# APPENDIX C: INTRINSIC ALIGNMENT MODEL

In this appendix, we consider the impact of the IA model choice on our results. Throughout this work, we use the NLA IA model, which is a subspace of the TATT model. A WL-only analysis using BCEmu7 with TATT constrains $S_8 = 0.802_{-0.024}^{+0.028}$, lying ~ 0.5σ lower than that constrained using NLA. This shift in $S_8$ between the two IA models is consistent when using BCEmu1 instead of BCEmu7, and similar to that found in previous work when using HM20 (DES & KiDS Collaborations 2023) and scale cuts (Secco et al. 2022; Amon et al. 2023). As in the literature, both IA model choices give comparable $\chi^2_{\text{red}}$ values (see Table 3). In future work, it is important to determine the more accurate IA model. Here, we use this test as validation that our baryon model choice and the results of our baryon model comparison are independent of the choice of IA model.

# APPENDIX D: SAMPLING ALGORITHM CHOICE

In this appendix, we compare the cosmological parameter estimates attained using the MULTINEST and POLYCHORD samplers. Throughout the work, we use the MULTINEST settings $n_{\text{live}} = 500$, efficiency = 0.3, tolerance = 0.1, constant efficiency = False, and max. iterations = 50 000 for computing efficiency. For POLYCHORD, we use $n_{\text{live}} = 500$, tolerance = 0.01, $n_{\text{repeats}} = 60$, and fast fraction = 0.1. Considering the BCEmu7 WL-only analysis, we find that sampling with POLYCHORD estimates a mean value of $S_8 = 0.817_{-0.025}^{+0.025}$, consistent with that attained by MULTINEST of $S_8 = 0.818_{-0.024}^{+0.017}$. However, in agreement with the findings of DES & KiDS Collaborations (2023), we find that the 68 per cent confidence level for $S_8$ attained using MULTINEST is 18 per cent smaller than that estimated with POLYCHORD. Similarly we find a WL + kSZ analysis with POLYCHORD returns $S_8 = 0.821_{-0.023}^{+0.020}$, consistent with the $S_8 = 0.823_{-0.020}^{+0.019}$ using MULTINEST but with a 9 per cent larger confidence region.

<footnote n="1">¹Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="2">²Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="3">³Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA</footnote>

<footnote n="4">⁴Institute for Computational Science, University of Zurich, Winterthur-erstrasse 190, CH-8057 Zurich, Switzerland</footnote>

<footnote n="5">⁵Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK</footnote>

<footnote n="6">⁶CIEMAT, Centre for Energy, Environmental and Technological Research, Av. Complutense, 40, Moncloa - Aravaca, E-28040 Madrid, Spain</footnote>

<footnote n="7">⁷Kavli Institute for Particle Astrophysics and Cosmology, 382 Via Pueblo Mall, Stanford, CA 94305-4060, USA</footnote>

<footnote n="8">⁸SLAC National Accelerator Laboratory, 2575 Sand, Hill Road, Menlo Park, CA 94025, USA</footnote>

<footnote n="9">⁹Physics Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA</footnote>

<footnote n="10">¹⁰Berkeley Center for Cosmological Physics, Department of Physics, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="11">¹¹Department of Astronomy, Cornell University, Ithaca, NY 14853, USA</footnote>

<footnote n="12">¹²Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="13">¹³Kavli Institute for the Physics and Mathematics of the Universe (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="14">¹⁴Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15312, USA</footnote>

<footnote n="15">¹⁵NSF AI Planning Institute for Physics of the Future, Carnegie Mellon University, Pittsburgh, PA 15213, USA</footnote>

<footnote n="16">¹⁶Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, PO Box 2450, Stanford, CA 94305, USA</footnote>

<footnote n="17">¹⁷SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA</footnote>

<footnote n="18">¹⁸Laboratório Interinstitucional de e-Astronomia - LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ-20921-400, Brazil</footnote>

<footnote n="19">¹⁹Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ-20921-400, Brazil</footnote>

<footnote n="20">²⁰Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA</footnote>

<footnote n="21">²¹Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain</footnote>

<footnote n="22">²²Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="23">²³Fermi National Accelerator Laboratory, PO Box 500, Batavia, IL 60510, USA</footnote>

<footnote n="24">²⁴Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="25">²⁵NASA Goddard Space Flight Center, 8800 Greenbelt Rd, Greenbelt, MD 20771, USA</footnote>

<footnote n="26">²⁶Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, Campinas, SP 13083-859, Brazil</footnote>

<footnote n="27">²⁷Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="28">²⁸Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="29">²⁹Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="30">³⁰CNRS, LPSC-IN2P3, Université Grenoble Alpes, F-38000 Grenoble, France</footnote>

<footnote n="31">³¹University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81679 Munich, Germany</footnote>

<footnote n="32">³²Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA</footnote>

<footnote n="33">³³Brookhaven National Laboratory, Bldg 510, Upton, NY 11973, USA</footnote>

<footnote n="34">³⁴Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093 Zurich, Switzerland</footnote>

<footnote n="35">³⁵Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona), Spain</footnote>

<footnote n="36">³⁶Department of Astronomy/Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA</footnote>

<footnote n="37">³⁷Department of Physics, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="38">³⁸School of Physics and Astronomy, Cardiff University, Wales, CF24 3AA, UK</footnote>

<footnote n="39">³⁹Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain</footnote>

<footnote n="40">⁴⁰Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, F-31400 Toulouse, France</footnote>

<footnote n="41">⁴¹Department of Physics and Astronomy, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada</footnote>

<footnote n="42">⁴²Institute for Astronomy, University of Edinburgh, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="43">⁴³Department of Physics, Northeastern University, Boston, MA 02115, USA</footnote>

<footnote n="44">⁴⁴Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA</footnote>

<footnote n="45">⁴⁵Jodrell Bank Center for Astrophysics, School of Physics and Astronomy, University of Manchester, Oxford Road, Manchester, M13 9PL, UK</footnote>

<footnote n="46">⁴⁶Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfvénsväg 12, SE-10691 Stockholm, Sweden</footnote>

<footnote n="47">⁴⁷Physics Department, University of Wisconsin-Madison, 2320 Chamberlin Hall, 1150, University Avenue Madison, WI 53706-1390, USA</footnote>

<footnote n="48">⁴⁸Department of Physics, University of Genova and INFN, Via Dodecaneso 33, I-16146 Genova, Italy</footnote>

<footnote n="49">⁴⁹Center for Astrophysical Surveys, National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA</footnote>

<footnote n="50">⁵⁰Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA</footnote>

<footnote n="51">⁵¹Department of Physics, Duke University, Durham, NC 27708, USA</footnote>

<footnote n="52">⁵²Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, UK</footnote>
<footnote n="53">⁵³Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="54">⁵⁴Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA</footnote>

<footnote n="55">⁵⁵Physics Department, William Jewell College, Liberty, MO 64068, USA</footnote>

<footnote n="56">⁵⁶Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA</footnote>

<footnote n="57">⁵⁷Department of Astronomy, University of Geneva, ch. d'Écogia 16, CH-1290 Versoix, Switzerland</footnote>

<footnote n="58">⁵⁸Department of Astronomy, University of California, Berkeley, 501 Campbell Hall, Berkeley, CA 94720, USA</footnote>

<footnote n="59">⁵⁹Cerro Tololo Inter-American Observatory, NSFs National Optical-Infrared Astronomy Research Laboratory, Casilla 603, La Serena, Chile</footnote>

<footnote n="60">⁶⁰Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK</footnote>

<footnote n="61">⁶¹CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014 Paris, France</footnote>

<footnote n="62">⁶²UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, Sorbonne Universités, F-75014 Paris, France</footnote>

<footnote n="63">⁶³Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="64">⁶⁴Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain</footnote>

<footnote n="65">⁶⁵Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029 Hamburg, Germany</footnote>

<footnote n="66">⁶⁶Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India</footnote>

<footnote n="67">⁶⁷Institute of Theoretical Astrophysics, University of Oslo. PO Box 1029 Blindern, NO-0315 Oslo, Norway</footnote>

<footnote n="68">⁶⁸Instituto de Fisica Teorica UAM/CSIC, Universidad Autonoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="69">⁶⁹School of Mathematics and Physics, University of Queensland, Brisbane QLD 4072, Australia</footnote>

<footnote n="70">⁷⁰Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA</footnote>

<footnote n="71">⁷¹Department of Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="72">⁷²Center for Astrophysics | Harvard and Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="73">⁷³Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia</footnote>

<footnote n="74">⁷⁴Lowell Observatory, 1400 Mars Hill Rd, Flagstaff, AZ 86001, USA</footnote>

<footnote n="75">⁷⁵George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA</footnote>

<footnote n="76">⁷⁶LPSC Grenoble - 53, Avenue des Martyrs, F-38026 Grenoble, France</footnote>

<footnote n="77">⁷⁷Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain</footnote>

<footnote n="78">⁷⁸Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada</footnote>

<footnote n="79">⁷⁹Faculty of Physics and Astronomy, Astronomical Institute, German Centre for Cosmological Lensing, Ruhr University Bochum, D-44780 Bochum, Germany</footnote>

<footnote n="80">⁸⁰Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote n="81">⁸¹School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK</footnote>

<footnote n="82">⁸²Physik-Institut, University of Zürich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland</footnote>

<footnote n="83">⁸³Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA</footnote>

<footnote n="84">⁸⁴Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, D-85748 Garching, Germany</footnote>

<footnote n="85">⁸⁵Fakultät für Physik, Ludwig-Maximilians Universität München, Universitäts-Sternwarte, Scheinerstr. 1, D-81679 München, Germany</footnote>

<footnote>This paper has been typeset from a TeX/Latex file prepared by the author.</footnote>