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

# Cosmological shocks around galaxy clusters: a coherent investigation with DES, SPT, and ACT

D. Anbajagane (துப்)¹,²* C. Chang¹,² E. J. Baxter,³ S. Charney,⁴ M. Lokken⁵,⁶,⁷ M. Aguena¹,⁸

S. Allam,⁹ O. Alves,¹⁰ A. Amon,¹¹,¹² R. An,¹³ F. Andrade-Oliveira,¹⁰ D. Bacon,¹⁴ N. Battaglia,¹⁵

K. Bechtol,¹⁶ M. R. Becker,¹⁷ B. A. Benson,¹,²,¹⁸ G. M. Bernstein,⁴ L. Bleem,²,¹⁹ S. Bocquet,²⁰

J. R. Bond,²¹ D. Brooks,²² A. Carnero Rosell,⁸,²³,²⁴ M. Carrasco Kind,²⁵,²⁶ R. Chen,²⁷ A. Choi,²⁸

M. Costanzi,²⁹,³⁰,³¹ T. M. Crawford,¹,² M. Crocce,³²,³³ L. N. da Costa,⁸ M. E. S. Pereira,³⁴ T. M. Davis,³⁵

J. De Vicente,³⁶ S. Desai,³⁷ M. J. Devlin,⁴ H. T. Diehl,⁹ P. Doel,²² C. Doux,⁴,³⁸ A. Drlica-Wagner,¹,²,⁹

J. Elvin-Poole,³⁹ I. Ferrero,⁴⁰ A. Ferté,⁴¹ B. Flaugher,⁹ P. Fosalba,³²,³³ D. Friedel,²⁶ J. Frieman,²,⁹

J. García-Bellido,⁴² M. Gatti,⁴ G. Giannini,⁴³ S. Grandis,⁴⁴ D. Gruen,⁴⁵ R. A. Gruendl,²⁵,²⁶ G. Gutierrez,⁹

I. Harrison,⁴⁶ J. C. Hill,⁴⁷ M. Hilton,⁴⁸,⁴⁹ S. R. Hinton,³⁵ D. L. Hollowood,⁵⁰ K. Honscheid,⁵¹,⁵² B. Jain,⁴

D. J. James,⁵³ M. Jarvis,⁴ K. Kuehn,⁵⁴,⁵⁵ M. Lin,⁴ N. MacCrann,⁵⁶ J. L. Marshall,⁵⁷ J. McCullough,⁵⁸

J. J. McMahon,¹,²,⁵⁹,⁶⁰ J. Mena-Fernández,³⁶ F. Menanteau,²⁵,²⁶ R. Miquel,⁴³,⁶¹ K. Moodley,⁴⁹,⁶²

T. Mroczkowski,⁶³ J. Myles,⁴¹,⁵⁸,⁶⁴ S. Naess,⁶⁵ A. Navarro-Alsina,⁶⁶ R. L. C. Ogando,⁶⁷ L. A. Page,⁶⁸

A. Palmese,⁶⁹ S. Pandey,⁴ B. Patridge,⁷⁰ A. Pieres,⁸,⁶⁷ A. A. Plazas Malagón,⁴¹,⁵⁸ A. Porredon,⁷¹,⁷²

J. Prat,¹,² C. Reichardt,⁷³ K. Reil,⁴¹ M. Rodriguez-Monroy,³⁶,⁷⁴ R. P. Rollins,⁷⁵ A. K. Romer,⁷⁶

E. S. Rykoff,⁴¹,⁵⁸ E. Sanchez,³⁶ C. Sánchez,⁴ D. Sanchez Cid,³⁶ E. Schaan,⁴¹,⁷⁷ M. Schubnell,¹⁰

L. F. Secco,² I. Sevilla-Noarbe,³⁶ E. Sheldon,⁷⁸ T. Shin,⁷⁹ C. Sifón,⁸⁰ M. Smith,⁸¹ S. T. Staggs,⁶⁸

E. Suchyta,⁸² M. E. C. Swanson,⁸³ G. Tarle,¹⁰ C. To,⁵² M. A. Troxel,²⁷ I. Tutusaus,⁸⁴ E. M. Vavagiakis,⁸⁵

N. Weaverdyck,¹⁰,⁸⁶ J. Weller,⁸⁷,⁸⁸ P. Wiseman,⁸¹ E. J. Wollack⁸⁹ and B. Yanny⁹

*Affiliations are listed at the end of the paper*

Accepted 2023 November 29. Received 2023 November 28; in original form 2023 September 30

## ABSTRACT

We search for signatures of cosmological shocks in gas pressure profiles of galaxy clusters using the cluster catalogues from three surveys: the Dark Energy Survey (DES) Year 3, the South Pole Telescope (SPT) SZ survey, and the Atacama Cosmology Telescope (ACT) data releases 4, 5, and 6, and using thermal Sunyaev–Zeldovich (SZ) maps from SPT and ACT. The combined cluster sample contains around $10^5$ clusters with mass and redshift ranges $10^{13.7} < M_{200m}/M_{\odot} < 10^{15.5}$ and $0.1 < z < 2$, and the total sky coverage of the maps is $\approx 15\,000\,\text{deg}^2$. We find a clear pressure deficit at $R/R_{200m} \approx 1.1$ in SZ profiles around both ACT and SPT clusters, estimated at $6\sigma$ significance, which is qualitatively consistent with a shock-induced thermal non-equilibrium between electrons and ions. The feature is not as clearly determined in profiles around DES clusters. We verify that measurements using SPT or ACT maps are consistent across all scales, including in the deficit feature. The SZ profiles of optically selected and SZ-selected clusters are also consistent for higher mass clusters. Those of less massive, optically selected clusters are suppressed on small scales by factors of 2–5 compared to predictions, and we discuss possible interpretations of this behaviour. An oriented stacking of clusters – where the orientation is inferred from the SZ image, the brightest cluster galaxy, or the surrounding large-scale structure measured using galaxy catalogues – shows the normalization of the one-halo and two-halo terms vary with orientation. Finally, the location of the pressure deficit feature is statistically consistent with existing estimates of the splashback radius.

**Key words:** galaxies: clusters: intracluster medium – large-scale structure of Universe.

<footnote>* E-mail: dhayaa@uchicago.edu</footnote>

<footnote>© 2023 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

Cosmological shocks are violent, high-energy phenomena that are a natural consequence of cosmic structure formation, and form in the far outskirts of massive, collapsed objects like galaxy clusters. They impact astrophysical processes like cosmic ray production and galaxy evolution, and are generated when colder gas is accreted on to a halo. The gravitational infall velocity of the cold gas will generically exceed the sound speed of the gas, especially for infall around massive haloes, and this results in a high Mach number shock ($M \sim 100$, e.g. Molnar et al. 2009).

The presence of such shocks impacts a wide array of astrophysical processes. These shocks are a natural thermodynamic boundary around the cluster, at the interface between the cluster-dominated gas component and the surrounding large-scale structure. They thereby also set the boundary within which the cluster has a thermodynamic impact on objects, such as galaxy quenching via ram-pressure stripping (e.g. Zinger et al. 2016; Boselli, Fossati & Sun 2022). Shocks are sites for accelerating cosmic ray electrons via Diffusive Shock Acceleration (Drury 1983; Blandford & Eichler 1987), and such accelerated cosmic ray electrons form a non-thermal tail in the energy distribution of the electron population (Miniati et al. 2001; Ryu et al. 2003; Brunetti & Jones 2014). The radial location of shock features also depends on the mass accretion rate of the cluster and can potentially serve as an observational proxy for the same (Lau et al. 2015; Shi 2016; Zhang et al. 2020, 2021). The mass accretion rate has strong theoretical connections to key dark matter halo properties such as the concentration and formation time (Wechsler et al. 2002), and has significant correlations with a wider range of halo properties (e.g. Lau et al. 2021; Anbajagane, Evrard & Farahi 2022a; Shin & Diemer 2023). However, it has remained difficult to infer observationally.

This process of shock heating generates a thermal non-equilibrium between the electrons and ions, which can alter the expected thermodynamic profiles and will consequently need to be considered in analyses that include these cluster outskirts (Fox & Loeb 1997; Ettori & Fabian 1998; Wong & Sarazin 2009; Rudd & Nagai 2009; Akahori & Yoshikawa 2010; Avestruz et al. 2015; Vink et al. 2015). Specifically, shocks preferentially heat ions over electrons given the mass difference of the two species, and at the low-number densities of the cluster outskirts, the two species may not interact often enough to equilibrate. This will lead to a deficit in the measured SZ profiles – which traces the *electron*, not ion, temperature – near a shock, and such a deficit has been observed previously with SPT data (Anbajagane et al. 2022c). In addition, an accurate model of these cluster outskirts – particularly near the transition regime between the bound component and the large-scale structure – will be beneficial for studies of the large-scale gas pressure fields (e.g. Hill & Pajer 2013; Horowitz & Seljak 2017; Tanimura et al. 2022) as well as cross-correlations of the gas pressure with galaxy and galaxy cluster positions (e.g. Hajian et al. 2013; Vikram, Lidz & Jain 2017; Hill et al. 2018; Pandey et al. 2019; Pandey, Baxter & Hill 2020; Sánchez et al. 2023), with weak-lensing shears (e.g. Ma et al. 2015; Hojjati et al. 2017; Osato et al. 2018, 2020; Shirasaki, Lau & Nagai 2020; Gatti et al. 2022; Pandey et al. 2022), or with X-ray luminosity (Shirasaki, Lau & Nagai 2020); these kinds of studies are positioned to provide strong and complementary constraints on astrophysical, as well as cosmological processes. The model will also be beneficial for understanding the impact from the gas dynamics of the outskirts on the weak lensing signal (via the impact of gas dynamics on the total matter field) – this impact is a significant limitation in extracting cosmological information from the lensing signal (e.g. Gatti et al. 2020; Krause et al. 2021; Secco et al. 2022a; Amon et al. 2022;

Anbajagane et al. 2023a, b) – and for subsequently modelling the impact via a halo-model approach (e.g. Schneider et al. 2019; Chen et al. 2023).

While a wide variety of physical processes are influenced by the presence of shocks, the cosmological shocks are themselves simple, as their formation has two basic requirements: a matter component that is collisional and thus behaves hydrodynamically ('gas'), and an influx of this collisional matter on to a halo via gravitational infall. However, both hydrodynamics and gravitational infall are highly asymmetric processes with complicated geometries, and so, in practice, these shocks have a rich phenomenology with intricate, subtle behaviours.

This phenomenology has been extensively studied in simulations over the past many decades. The first studies used non-radiative simulations with gas dynamics but no astrophysical processes (Quilis et al. 1998; Miniati et al. 2000; Ryu et al. 2003; Skillman et al. 2008; Molnar et al. 2009; Hong et al. 2014; Hong, Kang & Ryu 2015; Schaal & Springel 2015). These were then followed by studies using simulations that include gas cooling and star formation (Vazza, Brunetti & Gheller 2009; Planelles & Quilis 2013; Lau et al. 2015; Nelson et al. 2016; Aung, Nagai & Lau 2021), and also include the effects of feedback from supernovae and active galactic nuclei (Kang et al. 2007; Vazza, Brüggen & Gheller 2013; Vazza, Gheller & Brüggen 2014; Schaal et al. 2016; Baxter et al. 2021; Planelles et al. 2021; Baxter et al. 2023; Sayers et al. 2023). Some works have also opted to model the evolution of cosmic rays – which are generated at the shocks – alongside galaxy formation (Pfrommer et al. 2007), while others employ idealized simulations to understand the propagation of shocks and their dependence on different merger events (Pfrommer et al. 2006; Ha, Ryu & Kang 2018; Zhang et al. 2019b, 2020, 2021). A number of works have also theoretically estimated the potential signal-to-noise ratio (SNR) of shocks from various surveys/instruments (e.g. Kocsis, Haiman & Frei 2005; Baxter et al. 2021).

In the current picture, cosmological shocks form at different radial locations around the galaxy cluster depending on the mechanism that generates them. The accretion of pristine cold gas – which has a low sound speed and is primarily found in low-density regions such as cosmic voids – on to the thermalized, bound gas component results in a shock of a high Mach number ($M \sim 100$) and discontinuities in the profiles of many thermodynamic quantities such as temperature, entropy, pressure, and density. This shock – approximately located near the virial radius of the cluster – is oftentimes referred to as an accretion shock (e.g. Lau et al. 2015; Aung, Nagai & Lau 2021; Baxter et al. 2021) or an external shock (Ryu et al. 2003), and has a theoretical foundation that goes back many decades (Bertschinger 1985). Closer to the cluster core, the supersonic infall of galaxies and gas clumps into the hot, ionized gas leads to a series of bow shocks with weak Mach numbers, that are referred to as internal shocks (Ryu et al. 2003). Furthermore, Zhang et al. (2019b, 2020) found that these bow shocks detach from the infalling substructure, leading to a runaway merger shock that then collides with the accretion shock. This generates a new shock, named the *Merger-accelerated Accretion Shock* or MA shock, that is both further out and longer lived than the original accretion shock. The infall of substructure is a common process during structure formation, and so most shocks observed in the cluster outskirts are expected to be MA shocks and can have radial locations between $1 \le R/R_{200m} \le 2.5$ depending on the accretion history of the cluster (Zhang et al. 2021). These structures, given their origin in the large-scale accretion of matter, are connected to other features in the cluster outskirts such as the splashback radius (Adhikari, Dalal & Chamberlain 2014; Diemer & Kravtsov 2014).
This feature has been found in various data sets (Baxter et al. 2017; Chang et al. 2018; Shin et al. 2019, 2021; Adhikari et al. 2021) and its connection to cosmological shocks has been explored via both analytic calculations and simulations (Shi 2016; Aung, Nagai & Lau 2021; Baxter et al. 2021; Zhang et al. 2021).

While many simulation-based studies exist on the formation and evolution of these shocks, there are only a few observational studies of these features. A key observable for studying these shocks is the cluster gas pressure profiles, measured via the thermal Sunyaev-Zel'dovich (SZ) signature of clusters (Sunyaev & Zeldovich 1972). The SZ effect is the inverse Compton scattering of cosmic microwave background (CMB) photons off energetic electrons in the hot intra-cluster medium (see Carlstrom, Holder & Reese 2002; Mroczkowski et al. 2019 for reviews). While cluster thermodynamic properties have traditionally been studied using X-ray observations, the SZ effect has emerged as the more ideal probe for the cluster outskirts as its signal amplitude depends linearly with density, whereas for X-rays this dependence is quadratic. Many of the existing observational works – using either X-ray or SZ – do not explicitly focus on shocks and most are limited to small, often single, cluster samples at lower redshifts (Akamatsu et al. 2011; Akahori & Yoshikawa 2012; Akamatsu et al. 2016; Basu et al. 2016; Di Mascolo, Churazov & Mroczkowski 2019a; Di Mascolo et al. 2019b; Hurier, Adam & Keshet 2019; Pratt, Qu & Bregman 2021; Zhu et al. 2021). More general studies of gas thermodynamic profiles, without a specific focus on shocks, do not push beyond $r \gtrsim R_{500c}$ (e.g. McDonald et al. 2014; Ghirardini et al. 2017; Romero et al. 2017, 2018; Ghirardini et al. 2018), though some do exist (Planck Collaboration V 2013; Sayers et al. 2013, 2016; Amodeo et al. 2021; Schaan et al. 2021; Melin & Pratt 2023; Lyskova et al. 2023).

Anbajagane et al. (2022c), henceforth A22, performed the first analysis of the cluster outskirts with a large statistical sample of $10^2-10^3$ clusters, and found evidence of a pressure deficit at the cluster virial radius. This work is a follow-up on A22, and our goals are to (i) to strengthen the evidence for the pressure deficit with additional, sensitive SZ data, (ii) compare the SZ profiles and their pressure deficit feature, between SZ-selected and optically selected cluster catalogues, and between measurements from different SZ maps (ACT and SPT), (iii) measure cluster profile outskirts for lower mass clusters, $M_{200m} < 10^{14.5} M_{\odot}$, (iv) extract anisotropic features of the profile outskirts, using SZ image shapes, the brightest cluster galaxy (BCG) shapes, or the large-scale density field, and finally (v) compare the location of detected features with other physical cluster radii, namely the splashback radius. We achieve all of the above by expanding our study to include additional surveys: an optically selected cluster catalogue from the Dark Energy Survey (DES) Year 3 data set, and an SZ map from ACT Data Release (DR) 6. Both data sets were not used in the work of A22. The availability of the ACT DR6 map also allows us to now use the full ACT DR5 cluster catalogue, whereas A22 were limited to using a subset ($\approx 25$ per cent) of the catalogue that overlapped with the ACT DR4 map.

We organize this work as follows. In Section 2, we describe the survey data sets used in this work and in Section 3 our choices for the profile measurement procedure and the theoretical modelling. Our results on shocks are shown in Section 4 and their connections to other large-scale structure features are explored in Section 5. We conclude in Section 6.

## 2 DATA

We use data from three wide-field surveys – the DES Year 3, the South Pole Telescope (SPT) SZ survey, and the Atacama Cosmology

Telescope (ACT) DRs 4, 5, and 6 – to constrain the cluster pressure profile on large scales. In contrast to A22, we do not consider profiles from the Planck SZ map, though Planck data are used in the construction of the ACT and SPT maps (described below in Sections 2.2 and 2.3). The former choice is because the 10 arcmin resolution of the Planck SZ map (which is an order of magnitude larger than the 1 arcmin resolution of SPT and ACT) is a limiting factor in detecting shock features. The Planck cluster catalogue (Planck Collaboration XXVII 2016) also has significant overlap with the SPT and ACT catalogues used in this work; 45 per cent of the 1093 Planck clusters are found within either the ACT or SPT footprints.

The clusters in our samples are labelled by their spherical overdensity mass, $M_{200m}$, which is defined as

$$M_{\Delta} = \rho_{\Delta} \frac{4\pi}{3} R_{\Delta}^{3}, \qquad (1)$$

with $\rho_{\Delta} = 200\rho_{m}(z)$, where $\rho_{m}(z)$ is the mean matter density of the Universe at a given epoch. The associated radius is denoted as $R_{200m}$. Features at the cluster outskirts, such as shocks, follow a more self-similar evolution when normalized by this radius definition (Diemer & Kravtsov 2014; Lau et al. 2015).

Both SPT and ACT infer $M_{500c}$ from the integrated tSZ emission around each cluster, while DES infers $M_{500c}$ from the cluster richness, where richness is the probabilistic number of satellite galaxies in the cluster. We then convert the $M_{500c}$ estimate into $M_{200c}$ and $M_{200m}$ using the concentration-mass relation from Diemer & Joyce (2019) and the publicly available routine from the COLOSSUS<fnref n="1" /> open-source PYTHON package (Diemer 2018). We find our results are insensitive to assuming other choices for the concentration-mass relation (e.g. Child et al. 2018; Ishiyama et al. 2021). The impact of baryons on this relation is also negligible at these halo masses and so is not considered here (e.g. Beltz-Mohrmann & Berlind 2021; Anbajagane, Evrard & Farahi 2022a; Shao, Anbajagane & Chang 2023; Shao & Anbajagane 2023). Both $M_{500c}$ and $M_{200c}$ are defined by equation (1) but with alternative density contrasts of $\rho_{\Delta} = 500\rho_{c}(z)$ and $\rho_{\Delta} = 200\rho_{c}(z)$, respectively. Here, $\rho_{c}(z)$ is the critical density of the Universe at a given epoch. The mass and redshift distributions of the different cluster samples are shown in Fig. 1.

The tSZ amplitude is reported as the dimensionless y parameter,

$$y \equiv \frac{k_B \sigma_T}{m_e c^2} \int n_e T_e dl, \qquad (2)$$

where $k_B$ is the Boltzmann constant, $\sigma_T$ is the Thomson cross-section, $m_e c^2$ is the rest energy of an electron, $n_e$ and $T_e$ are the electron number density and temperature, respectively, and $l$ is the physical line-of-sight distance. Thus, y represents the electron pressure integrated along the line of sight.

The tSZ effect corresponds to CMB photons scattering off electrons with a thermal (i.e. Maxwellian) energy/momentum distribution. There exist similar effects, called the relativistic SZ (rSZ) and non-thermal SZ (ntSZ), which correspond to photons scattering off electrons with non-Maxwellian energy distributions, and may leak into the measured tSZ signal (Mroczkowski et al. 2019). In the rSZ effect, the presence of high-temperature electrons ($T_e \gtrsim 5$ keV) requires relativistic corrections to the procedure for making the SZ maps. These corrections, however, are $\lesssim 5$ per cent (see fig. 1 of Erler et al. 2018) and are subdominant to the amplitudes of the features discussed in our work.<fnref n="2" /> The ntSZ effect can be generated

<footnote n="1">¹https://bdiemer.bitbucket.io/colossus/</footnote>

<footnote n="2">²The work of Lee et al. (2022) shows the rSZ effect in simulations scales self-similarly as $\propto M^{2/3}$, or alternatively $\propto Y^{2/5}$, and so given our cluster</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The mass–redshift plane of the cluster samples from SPT, ACT, and DES used in this work. The *Planck* catalogue is shown in grey for reference. The top and right panels show the 1D distributions for redshift and cluster mass, respectively. For visibility, we only plot a randomly chosen subset of DES clusters, with *N* = 5000. The 1D distributions are estimated using the full samples. The SPT and ACT samples have similar redshift distributions, with a median of *z* ≈ 0.55, while DES Y3 is limited to 0.1 < *z* < 0.8. DES also extends to much lower masses across all redshifts, where the masses are computed using the mass–richness relation of Costanzi et al. (2021, see their equation 16). The colour tones of the points show log₁₀ SNR, the SNR of each cluster detection, with lighter colours indicating a higher SNR. The mean redshift and mass of the different samples are listed in Table 1.</figcaption></figure>

by a cosmic ray electron population, but is a subdominant effect within R<sub>200c</sub> of the cluster, where cosmic rays make-up ≲ 1 per cent of the total pressure (Ackermann et al. 2014). Beyond this radius, the cosmic ray energy fraction is not well constrained. For this work, we follow A22 in assuming the ntSZ continues to be subdominant in the outskirts, and point out that the features we discuss are unaffected even if the ntSZ contaminates the tSZ at the 10 per cent level.

## 2.1 The DES Year 3

DES Y3 is a 5000 deg² photometric survey of the southern sky in five bands (*grizY*). Galaxy clusters are identified using the REDMAPPER algorithm (Rykoff et al. 2014), which identifies clusters from overdensities of red-sequence galaxies. Each cluster is assigned a ‘richness’, $\lambda$, which is analogous to the number of red galaxies in the cluster. REDMAPPER assigns each galaxy *i* a probability that it is a satellite of galaxy cluster *j*. The richness of cluster *j* is then the sum of these probabilities.

This richness is used alongside a richness–mass relation – which can be calibrated using various methods such as galaxy lensing (McClintock et al. 2019), CMB lensing (Baxter et al. 2018), cross-correlations of probes (To et al. 2021), galaxy velocity dispersion (Farahi et al. 2016; Anbajagane et al. 2022b), etc. – to obtain a mass estimate for each cluster. In this work, we use the richness–mass relation from Costanzi et al. (2021, see their equation 16), which is calibrated using a combination of optical and SZ cluster measurements – namely, the DES cluster number counts and the SPT

sample spans across an order-of-magnitude in mass, the rSZ effect would change at most by a factor of two across our cluster sample. Note, however, that this is a factor of 2 difference in an effect that contributes < 5 per cent to the total signal.

observable-mass relation – for clusters with $\lambda \ge 20$. The observable–mass relation was in turn calibrated with targeted weak-lensing measurements. Note that the catalogues we use have objects of lower richness ($\lambda \approx 10$) and thus the inferred mass of these objects could be biased given we must extrapolate the scaling relation of Costanzi et al. (2021) to this regime. There are no well-calibrated richness–mass relations in this regime, and thus extrapolation is a necessity. In Section 4.3, we discuss the impact of such mass biases in our analysis.

We also use a cluster SNR as a weight when averaging the profiles across the sample (see Section 3.1). For DES, this signal-to-noise is taken to be the ratio of the richness over the richness uncertainty, $\lambda/\Delta\lambda$, where richness and the uncertainty are taken from the REDMAPPER columns LAMBDA_CHISQ and LAMBDA_CHISQ_E, respectively.

Finally, we also use two different galaxy samples to enable oriented stacking of the cluster profiles. First, we use the DES Y3 source galaxy shape catalogue (Gatti et al. 2021) – where the shapes were measured using the METACALIBRATION code (Sheldon & Huff 2017) – to obtain the orientation of the BCG of each cluster.<fnref n="3" /> Then, we also use the magnitude-limited lens galaxy catalogue, MAGLIM (Porredon et al. 2021), to infer the density field in the DES footprint, from which we can estimate a cluster orientation based on large-scale structure. This follows the methods of Lokken et al. (2022), and is discussed further in Section 5.1. Both data sets are part of the publicly available DES Y3 DR.<fnref n="4" />

## 2.2 The SPT SZ Survey

SPT-SZ is a 2500 deg² survey of the southern sky at 95, 150, and 220 GHz, and was conducted using the SPT (Carlstrom et al. 2011). The SZ map used in our analysis was presented in Bleem et al. (2022), has an angular resolution of 1.25 arcmin, and is made using data from both SPT-SZ and the *Planck* 2015 DR; the former provides lower noise measurements of the small scales, whereas the latter does the same for larger scales (multipoles $\ell \lesssim 1000$). The *Planck* data consists of the 100, 143, 217, and 353 GHz maps from the high-frequency instrument. The SZ map is constructed with the Linear Combination (LC) algorithm (see Delabrouille & Cardoso 2009 for a review), applied to the maps of different frequencies. The weights of the LC are chosen so as to minimize the total variance in the output map. The weights are also modified to reduce contamination from the cosmic infrared background [CIB; see section 3.5 in Bleem et al. (2022) for more details]. In our analysis, the map is further masked to remove point sources as well as the top 5 per cent of map regions most dominated by galactic dust. This is done using the binary masks provided in Bleem et al. (2022, see point 4 in their Appendix A).

The galaxy cluster catalogue from this data contain 516 clusters that were first identified in Bleem et al. (2015), and were assigned updated redshifts and mass estimates in Bocquet et al. (2019). We use the latter, updated catalogue for our work, where the mass is estimated via a joint modelling of SZ, X-ray, and weak lensing measurements. Both the map and the cluster catalogue are publicly available.<fnref n="5" /> Our masses come from the M500 column and SNR from the XI column. This data set is the exact same as the SPT-SZ data used in A22.

<footnote n="3">³We have verified that using alternative shape measurements, such as those from the single object fitting procedure (Sevilla-Noarbe et al. 2021), results in similar orientations for the galaxies.</footnote>

<footnote n="4">⁴https://des.ncsa.illinois.edu/releases/y3a2</footnote>

<footnote n="5">⁵https://lambda.gsfc.nasa.gov/product/spt/spt.prod.table.cfm</footnote>
## 2.3 ACT DRs 4, 5, and 6

The ACT data covers 90, 150, and 220 GHz frequencies, and the maps from DR 6 cover ≈ 13,000 deg² of the sky (after applying the relevant masks; see discussion below). The SZ map (Coulton et al. 2023) has a resolution of 1.6 arcmin and makes use of data from both ACT and the Planck NPIPE DR (Planck Collaboration LVII 2020); as was the case with SPT, the former data inform small-scales and the latter, the large-scales (ℓ ≲ 1000). Note that the Planck data here consist of eight frequency channels from 30 to 545 GHz, whereas the map from Bleem et al. (2022) used four of these channels. The map is made using a Needlet Internal Linear Combination algorithm.

In our analysis, the map is further masked to remove point sources and dusty regions. The ACT DR6 mask is an apodized, continuous mask, not a binary one, and we continue with our aggressive masking by only selecting pixels for which the mask value is 1, meaning the impact of point sources and dust is negligible in this pixel. Note that this map does not use the HEALPIX pixelation scheme implemented in HEALPY and instead uses the Plate Carrée scheme implemented in PIXELL,<fnref n="6" /> a package optimized to work with partial sky maps in the flat-sky approximation. We use the ACT DR6 map in its native scheme and do not convert it to a HEALPIX format.

We also use the ≈4200 clusters from ACT DR5⁷ catalogue (Hilton et al. 2021), which covers the same area as the ACT DR6 map. Note that only the subset of the ACT DR5 catalogue, that corresponded to the 2000 deg² area of the ACT DR4 map was used in A22. The redshift distribution of the ACT DR5 cluster sample is similar to that of the SPT-SZ sample. As was the case in A22, the cluster masses come from the M500cCal column described in Hilton et al. (2021, see their table 1), which contains a weak lensing mass calibration factor. While other lensing-based calibrations also exist for the ACT data (e.g. Robertson et al. 2023), we use the fiducial calibration included in the catalogue of Hilton et al. (2021). The SPT and ACT masses are similar (e.g. Hilton et al. 2021, see their section 5.1), with the agreement at a level adequate for astrophysical analyses.

# 3 MEASUREMENT AND MODELLING

We first describe our procedure for measuring the stacked SZ profile in Section 3.1, and then in Section 3.2 the theoretical halo model we compare the measurements with, including how we quantify the significance of any features in the data.

## 3.1 Measurement procedure

Our measurement procedure closely follows that described in A22, with some notable changes. We reproduce the main aspects of the measurement here for completeness but also point readers to A22 for a more detailed discussion on some elements of the procedure. Overall, the measurement procedure can be broken into four steps: (i) stacked profiles, (ii) logarithmic derivatives, (iii) bin-to-bin covariance matrix, and (iv) feature locations.

**Estimating stacked profiles:** For each cluster, we compute the $\langle y \rangle$ profile in 50 logarithmically spaced radial bins in the range $r \in [0.1, 20] R_{200m}$. We convert between angular and physical scales using the angular diameter distance estimated at the redshift of each cluster. The profile also has a mean background value subtracted from it. Previously, this background was estimated by measuring

the average profile around uniform random points across the whole map. This method was adequate for maps with mostly homogeneous survey properties, but can cause biases for maps with inhomogeneous survey properties, such as ACT DR6 where some regions of the sky are observed to significantly higher depth than other regions. We have thus updated our background subtraction procedure to capture this inhomogeneity. We take the region spanned by the cluster catalogue, and split it into different 'tiles' based on HEALPIX pixelization of NSIDE = 4. We have verified that our results below are robust if we instead use NSIDE = 8 or NSIDE = 16. We continue using NSIDE = 4 for our analysis given it is computationally cheaper. Once we tile the maps, we estimate the background separately in each tile by measuring profiles around all random points in the chosen tile. During background subtraction for a given cluster, we choose the background profile of the tile closest to that cluster.<fnref n="8" /> Previously, all clusters had a common background profile subtracted from them, whereas now the subtracted profile varies across the sky.

In A22, we did not consider the contamination in a cluster's measured profile due to interloper clusters in the foreground/background. Interlopers are distant in physical, three-dimensional (3D) space but appear close in projected, two-dimensional (2D) space. We have explicitly checked this effect – by masking out all potential interlopers when measuring the profiles of a given cluster – and found it does not impact the features we discuss in this work. In our test, an interloper is defined as any cluster whose line-of-sight distance from the target cluster is $R > 20R_{200m}$. An object with a large line-of-sight separation from a given cluster is not part of the latter's local large-scale environment but can appear so in projected 2D space where the line-of-sight separation is not relevant. Thus, selecting clusters where the line-of-sight separation is greater than $20R_{200m}$ isolates such interlopers. The choice of $20R_{200m}$ is because that is the largest radius we measure the profiles to. We convert the cluster redshift to physical distance assuming a fiducial Lambda cold dark matter cosmology with $\Omega_m = 0.3$ and $h = 0.7$, and use the distances to identify the interlopers. Photometric redshift uncertainties and cluster line-of-sight peculiar velocities will affect the accuracy of the distance estimate. Even so, this test is useful as an approximate check of the interlopers' impact. For our main analysis below, we do not perform any interloper masking as we have confirmed it is a negligible effect.

The profiles of the individual clusters are then stacked, with each profile being weighted by the corresponding cluster's SNR. Performing a standard average/stack with no weights does not change the result (see appendix A in A22). Note that for a given cluster, any radial bin that did not have any pixels in it – most commonly the case in the cores of high redshift clusters due to the limited angular resolution – is masked, and thus ignored, during the stacking. The uncertainty of the stacked profile is obtained through a leave-one-out jackknife resampling. The ith jackknife sample of the stacked profile can be written as

$$ \langle y \rangle_i(r) = \frac{1}{W_i(r)} \sum_{j \neq i}^{N_{\text{el}}} y_j(r) w_j \delta_j(r), \quad (3) $$

$$ W_i(r) = \sum_{j \neq i}^{N_{\text{el}}} w_j \delta_j(r), \quad (4) $$

<footnote n="8">⁸Alternatively, one could also produce a catalogue of random points that sample the sky in a manner consistent with the cluster catalogue of a given survey, and this can be produced by using maps of multiple survey properties. We have pursued our inhomogeneous background subtraction method as it can be performed without requiring this additional data product.</footnote>

<footnote n="6">⁶https://pixell.readthedocs.io/en/latest/</footnote>

<footnote n="7">⁷https://lambda.gsfc.nasa.gov/product/act/actpol_dr5_szcluster catalogue.info.html</footnote>
where $w_j$ is the SNR per cluster used in the weighted average, $\delta_j(r)$ is 1 if the data point for radius $r$ in cluster $j$ is unmasked and 0 otherwise, $N_{cl}$ is the total number of clusters. In this notation, $\langle y \rangle_i$ is the mean profile of the sample with cluster $i$ removed, and $y_j$ is the individual profile measurement from cluster $j$. The variance on the mean profile is then given by

$$ \sigma^2(r) = \frac{N(r) - 1}{N(r)} \sum_{j=1}^{N_{cl}} \left( \langle y \rangle_j(r) - \langle \bar{y} \rangle(r) \right)^2 \delta_j(r), \quad (5) $$

$$ N(r) = \sum_{j=1}^{N_{cl}} \delta_j(r), \quad (6) $$

where $\langle \bar{y} \rangle$ is the mean of the distribution of jackknife estimates computed in equation (3). Note that equation (5) has an additional factor of $N-1$ compared to the traditional definition of the variance, as required when using a jackknife estimator for the variance.

**Estimating logarithmic derivatives:** Shocks are generally characterized by sharp changes in thermodynamic quantities, and have been identified in some previous works as the point of steepest descent in the pressure profiles (e.g. Aung, Nagai & Lau 2021; Baxter et al. 2021). This corresponds to measuring minima in the logarithmic derivative. Derivatives, however, are affected by noise and we alleviate this by smoothing the stacked profiles with a Gaussian of width $\sigma_{lnr} = 0.16$, which is 1.5 times the logarithmic bin width, $\Delta \ln r \approx 0.11$. All profiles are smoothed by this scale, and we present results only for the range $0.3 < R/R_{200m} < 10$ which does not contain any edge effects due to the smoothing. A22 (see their Appendix A) have already shown that smoothing choices have negligible impact on the final results.

The log-derivative of the smoothed mean profile is computed using a five-point method,

$$ \frac{df}{dx} = \frac{-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)}{12h}, \quad (7) $$

where $f$ is an arbitrary function of $x$, and $h = \Delta \ln r$ is the spacing between the sampling points. We estimate the uncertainty on the log-derivative by computing equation (7) for every jackknifed mean profile and taking the standard deviation of the resulting distribution. An extra multiplicative factor of $\sqrt{N-1}$ is applied to convert the measured uncertainty to the unbiased uncertainty, and this is analogous to the extra $N-1$ factor used in the variance estimator, as shown in equation (5).

**Covariance of the log-derivative:** To compute a detection significance for any feature, we require the bin-to-bin covariance matrix, $C$, of the measured mean log-derivative, as is discussed further below in equation (22). This covariance is estimated using a jackknife sampling of the profiles,

$$ C_{ij} = \frac{N(r) - 1}{N(r)} \sum_{k=1}^{N_{cl}} \left( f'_{k,i} - \langle f' \rangle_i \right) \left( f'_{k,j} - \langle f' \rangle_j \right) \delta_j, \quad (8) $$

where $i$ and $j$ index over the different radial bins, $f'_{k,i}$ is the log-derivative of the mean profile in the $i$th bin for the $k$th jackknifed sample. All quantities in the sum are implicit functions of radius, and we have suppressed the notation for brevity. The correlation matrix is shown in Fig. D1.

**Quantifying feature location:** We are interested in the location of a given feature – particularly, of local minima in the log-derivative – and this is estimated by fitting cubic splines to the log-derivative of each mean profile in the jackknifed sample and then locating the feature of interest in each profile. The mean and standard deviation of the resulting distribution provide estimates of the location of the

feature and the associated uncertainty. Given our use of the jackknife method to estimate the uncertainty, the $\sqrt{N-1}$ factor is needed once again to convert from the measured uncertainty to the unbiased uncertainty. For the SZ-selected samples, the median uncertainty in $M_{200m}$ (as determined from the catalogues) is 15 per cent, and so the uncertainty in $R_{200m}$ is around 5 per cent. This is tolerable as it increases the total uncertainty in the estimated feature location by < 2 per cent. Note that the uncertainty in the feature location comes from variations in the shape of the profile. This depends both on the raw signal-to-noise of the measurement and on the intrinsic shape of the profiles. Thus, profiles that appear noisy can still have precise feature locations if the shape of the profile has less variation.

## 3.2 Modelling and detection quantification

As was done in A22, we look for features in the profile outskirts by comparing the measurements with theoretical predictions. The model we employ here for the halo-$y$ correlation follows that used in A22 with some changes that we highlight.

The model consists of two components: a one-halo term given by the projected version of the pressure profile from Battaglia et al. (2012), who calibrated the profiles using hydrodynamical simulations, and a two-halo term which accounts for contributions from nearby haloes as described in Vikram, Lidz & Jain (2017) and later in Pandey et al. (2019). The two-halo term prediction uses a linear matter power spectrum and linear halo bias, and assumes higher order corrections are not required. We have validated this assumption in A22 checking the model matches the two-halo term of profiles from THE THREE HUNDRED simulations (Cui et al. 2018, 2022). The entire model is implemented in the CORE COSMOLOGY LIBRARY (CCL) open-source PYTHON package<fnref n="9" /> (Chisari et al. 2019) and is public.<fnref n="10" />

We begin by representing the 3D, halo-pressure cross-correlation function as a composition of the one-halo and two-halo components,

$$ \xi_{h,p}(r, M, z) = \xi_{h,p}^{\text{one-halo}}(r, M, z) + \xi_{h,p}^{\text{two-halo}}(r, M, z), \quad (9) $$

where $\xi$ are the correlation functions, $r$ is comoving distance, and $M$ is the halo mass. We denote the combined one-halo and two-halo term as the ‘total halo model’. The one-halo term is obtained via the pressure profile of Battaglia et al. (2012),

$$ P(x) = P_{200c} P_0 \left( \frac{x}{x_c} \right)^\gamma \left[ 1 + \left( \frac{x}{x_c} \right)^{\alpha\gamma} \right]^{-\beta}, \quad (10) $$

where $P_0$, $x_c$, $\alpha$, $\beta$, and $\gamma$ are the fit parameters calibrated from hydrodynamical simulations, $x = r/R_{200c}$ is the distance in units of cluster radius, and $P_{200c}$ is the thermal pressure expectation from self-similar evolution,

$$ P_{200c} = 200 \rho_c(z) \frac{\Omega_b G M_{200c}}{\Omega_m 2 R_{200c}}. \quad (11) $$

Equation (10) accounts for deviations from self-similar evolution via the calibrated mass and redshift dependencies of the parameters $P_0$, $x_c$, and $\beta$. The model also includes the effects of non-thermal pressure support within haloes – which is generated by the incomplete thermalization of gas – as it is calibrated on simulations that include this phenomenon. The fit parameters for equation (10) are obtained from the ‘200 AGN’ calibration model of Battaglia et al. (2012, see Table 1), and these parameters have a known, calibrated scaling

<footnote n="9">⁹https://github.com/LSSTDESC/CCL</footnote>

<footnote n="10">¹⁰https://github.com/DhayaaAnbajagane/tSZ_Profiles</footnote>
with both cluster redshift, z, and cluster mass, M$_{200c}$. The calibration matches simulations within < 10 per cent in the one-halo regime (Battaglia et al. 2012, see their fig. 2 and section 4.2). While A22 used the ‘500 SH’ model, the ‘200 AGN’ model opted for here provides a better fit to the measured profiles on small-scales and is the model choice for other works that we compare to below (e.g. in Section 4.3). The pressure deficit we discuss below is observed regardless of the model chosen to be the comparison point.

The tSZ emission is connected to the electron pressure, P_e, whereas the profiles of Battaglia et al. (2012) are calibrated to the total gas pressure, P. We convert between them as

$$P_e(r, M, z) = \frac{4 - 2Y}{8 - 5Y} P(r, M, z), \qquad (12)$$

with Y = 0.24 being the primordial helium mass fraction. This provides our one-halo term,

$$\xi_{h,p}^{\text{one-halo}}(r, M, z) = P_e(r, M, z). \qquad (13)$$

It is more convenient to compute the two-halo term in Fourier space, so our computations are done in the same. We inverse Fourier transform the model in the end to obtain the required real-space correlation function. The two-halo term of the halo-pressure cross-power spectrum, $P_{h,p}^{\text{two-halo}}$, is written as

$$P_{h,p}^{\text{two-halo}}(k, M, z) = \left[ b(M, z) P_{\text{lin}}(k, z) \times \int_{0}^{\infty} dM' \frac{dn}{dM'} b(M', z) u_{p}(k, M', z) \right], \quad (14)$$

where M is the mass of the halo we are computing the halo-pressure correlation for, M' is the mass of a neighbouring halo contributing to the two-halo term, $P_{\text{lin}}(k, z)$ is the linear, matter density power spectrum at redshift z, $dn/dM'$ is the mass function of neighbouring haloes, and $b(M, z)$ and $b(M', z)$ are the linear bias factors for the target halo and neighbouring haloes, respectively. The mass function model comes from Tinker et al. (2008) and the linear halo bias model from Tinker et al. (2010). The term $u_p(k, M', z)$ is the Fourier transform of the pressure profile about the neighbouring halo which, under the assumption of spherical symmetry, is computed as

$$u_p(k, M', z) = \int_0^\infty dr 4\pi r^2 \frac{\sin(kr)}{kr} P_e(r, M', z), \qquad (15)$$

where $P_e$ is the electron pressure profile. The halo-pressure two-point cross-correlation is obtained as the inverse Fourier transform of the cross-power spectrum,

$$\xi_{h,p}^{\text{two-halo}}(r, M, z) = \int_{0}^{\infty} \frac{dk}{2\pi^2} k^2 \frac{\sin(kr)}{kr} P_{h,p}^{\text{two-halo}}(k, M, z). \qquad (16)$$

The terms shown in equations (13) and (16) can be combined according to equation (9) to get the total halo model, $\xi_{h,p}$.

We have thus far described the real-space 3D pressure, whereas the Compton-y parameter is the integrated (or projected) pressure along the line of sight. The halo-y correlation is therefore obtained by a projection integral,

$$\xi_{h,y}(r, M, z) = \frac{\sigma_T}{m_e c^2} \int_{-\infty}^{\infty} \frac{d\chi}{1+z} \xi_{h,p} \left( \sqrt{\chi^2 + r^2}, M, z \right), \qquad (17)$$

where $\sigma_T$ is the Thomson scattering cross-section, $m_e c^2$ is the rest mass energy of the electron, and $\chi$ is the comoving coordinate along the line of sight.

All SZ maps have a finite angular resolution, where the resolution limitation suppresses power on small scales. We incorporate this into

our model by smoothing the prediction. We first calculate the angular cross-power spectrum, using the flat sky approximation, as

$$C_{\ell} = \int d\theta 2\pi\theta J_0(\ell\theta) \xi_{h,y}(\theta, M, z), \qquad (18)$$

where $J_0$ is the zeroth-order Bessel function. We then multiply $C_\ell$ by the Fourier-space smoothing function for the given survey of interest and then perform an inverse-harmonic transform,

$$\xi_{h,y}^{\text{smooth}}(\theta, M) = \int \frac{d\ell\ell}{2\pi} J_0(\ell\theta) C_{\ell} B_{\ell}, \qquad (19)$$

with the smoothing function $B_\ell$ given as

$$B_{\ell} = \exp \left[ -\frac{1}{2} \ell(\ell + 1) \sigma_{\text{FWHM}}^2 \right], \qquad (20)$$

where $\sigma_{\text{FWHM}} = \theta_{\text{FWHM}}/\sqrt{8\ln 2}$, with $\theta_{\text{FWHM}} = 1.25$ arcmin ($\theta_{\text{FWHM}} = 1.6$ arcmin) being the full-width half-max of the Gaussian filter used to smooth the SPT (ACT) maps.

Our final theory curve for a given cluster sample is obtained as follows: we compute the smoothed total halo model, $\xi_{h,y}^{\text{smooth}}$, for each individual cluster in our catalogue, and then perform a weighted stack identical to that done on the data, i.e. where the weights are the SNR of the observed clusters. The only inputs to this model are the cluster mass, redshift, and SNR (which is used as a weight). Thus, the theoretical curves shown below are true predictions and are not model fits made on the profile measurements. The one exception is the model for DES clusters, which includes a miscentring component (described in Section 3.2.1) which does have free parameters that we vary. The approach to fixing those parameters is described in that same section. We generally only discuss results for DES clusters that do not require a theoretical model.

Finally, we estimate the significance of any deviation between the measured log-derivatives and the theoretical model as

$$\epsilon \equiv \frac{1}{\sigma} \left( \frac{d \ln y^{\text{obs}}}{d \ln x} - \frac{d \ln y^{\text{th}}}{d \ln x} \right), \qquad (21)$$

where $\sigma$ is the uncertainty in the log-derivative measurement. The quantity $\epsilon$ is the number of sigma by which the log-derivative in the data differs from that of the theory.

We also measure a standard chi-squared significance for the feature of interest as a whole,

$$\chi^2 = \left( \frac{d \ln y^{\text{obs}}}{d \ln x} - \frac{d \ln y^{\text{th}}}{d \ln x} \right)^T C^{-1} \left( \frac{d \ln y^{\text{obs}}}{d \ln x} - \frac{d \ln y^{\text{th}}}{d \ln x} \right), \qquad (22)$$

where $C^{-1}$ is the inverse of covariance matrix for the log-derivative, accounting for the Hartlap factor (Hartlap, Simon & Schneider 2007) as

$$C^{-1} \rightarrow \frac{N_{jk} - N_{\text{bin}} - 2}{N_{jk} - 1} C^{-1}. \qquad (23)$$

where $N_{jk}$ are the number of jackknife samples (more than 500 for almost all samples), $N_{\text{bins}} = 5$ are the number of bins used to estimate the significance of a particular feature (i.e. the pressure deficit). The rescaling accounts for the bias due to limited realizations being used to numerically estimate the covariance matrix. The covariance $C$ is defined in equation (8). As mentioned above, we do not use all 50 radial bins for this calculation and instead limit ourselves to all bins whose radii are within $\Delta\log_{10}r = 0.1$ of the location of the feature. Once the $\chi^2$ is computed, we quote the total signal-to-noise of a feature, as

$$\chi_{\text{sh}} = \sqrt{\chi^2 - N_{\text{dof}}} \qquad (24)$$
following the definition of Secco et al. (2022b, see their Equation C15), with $N_{\text{dof}} = 5$ as mentioned above. This definition of signal to noise improves on that used in A22 as it is more robust to noise fluctuations and binning choices.

### 3.2.1 Miscentring model for optically selected clusters

An additional component to our theoretical model, in comparison to that of A22, is the impact of cluster miscentring. For SZ-selected clusters, the offset between the cluster centre and the true centre (called ‘miscentring’) is negligible when compared to the radial scale of features we study, which are ~$R_{200m}$. When using optically selected clusters, however, the optically determined centre can be significantly offset from the centre of the gas distribution (Sehgal et al. 2013; Zhang et al. 2019a; Bleem et al. 2020).<fnref n="11" /> The impact of miscentring in the profile is to transfer power from small scales to large scales. The total observed profile, with miscentring, can be modelled as

$$y(R) = (1 - f_{\text{miscen}})y^{\text{true}}(R) + f_{\text{miscen}}y^{\text{miscen}}(R), \quad (25)$$

where $f_{\text{miscen}}$ is the fraction of miscentred objects and $y^{\text{true}}(y^{\text{miscen}})$ is the profile of correctly centred (miscentred) clusters. For a given miscentring offset, $R_{\text{mis}}$, the average miscentred profile is

$$y^{\text{miscen}}(R | R_{\text{mis}}) = \int_{0}^{2\pi} d\theta y^{\text{true}}\left(\sqrt{R^2 + R_{\text{mis}}^2 + 2 \cos \theta R R_{\text{mis}}^2}\right), \quad (26)$$

and the total model is obtained by marginalizing over the distribution of possible offsets,

$$y^{\text{miscen}}(R) = \int dR P(R_{\text{mis}}) y^{\text{miscen}}(R | R_{\text{mis}}). \quad (27)$$

Following previous works (e.g. Baxter et al. 2017; Chang et al. 2018; Shin et al. 2019, 2021), we assume the offsets follow a Rayleigh distribution,

$$P(R_{\text{mis}}) = \frac{R_{\text{mis}}}{\sigma_{R}^{2}} \exp \left[ - \frac{R_{\text{mis}}^{2}}{2\sigma_{R}^{2}} \right], \quad (28)$$

$$\sigma_R = \tau_{\text{miscen}} \left( \frac{\lambda}{100} \right)^{0.2} \text{Mpc.} \quad (29)$$

where $\lambda$ is the cluster richness. The free parameters of this model are $f_{\text{miscen}}$ and $\tau_{\text{miscen}}$ which set the fraction of miscentred objects, and the amplitude of the miscentring offset, respectively. The impact of miscentring – and the choice of the parameter values – for DES cluster profile model is discussed in Section 4.1 and further in Appendix A.

## 4 SHOCKS IN GALAXY CLUSTERS

We first present our main results in Section 4.1 using the cluster samples of the different surveys, then study the variation of the profiles (i) with cluster selection and choice of SZ map in Section 4.2, and; (ii) with halo mass, towards group-scale haloes, in Section 4.3. We will use the format CATALOG x MAP as a shorthand reference

for measurements for a given cluster catalogue using a given SZ map (e.g. SPT x SPT, DES x ACT).

All bands show 68 per cent uncertainties estimated via jackknife resampling of the profiles. As for the detection significance, we show $\epsilon$ in the figures but quote $\chi_{\text{sh}}$ in our discussions in the text as the total signal-to-noise of a feature. These are defined in equations (21) and (22), respectively. The latter is the combined significance of the feature across multiple radial bins, while the former is the single-bin significance and is useful for identifying the radial range of a signal.

Constraints on feature locations and their corresponding detection significance are provided in Table 1. In general, the measured location of the feature is expected to be offset from the true location due to the impact of beam smoothing in the SZ maps. However, we have verified previously, using simulations, that this difference is negligible for the SPT and ACT resolution level (A22). Note that, for the average cluster in our samples, the scale of $R_{200m}$ is a factor of ≈5 larger than the full-width half-max of the smoothing scale in these maps.

While the specific focus of this work is on finding pressure deficits and other shock-induced features in the SZ profile outskirts, this focus also requires we discuss profile behaviours in the one-halo and two-halo regimes. Shocks occur at the transition between the bound halo component (one-halo term) and the surrounding large-scale structure (two-halo term), so studying shock-induced features also requires studying these regimes. Thus, some of our discussions below will include behaviours of the one-halo and two-halo terms, as changes in these terms affect the overall shape of the halo profile.

### 4.1 Measurements from fiducial cluster samples

In Fig. 2, we present the average SZ profiles of different cluster samples measured using different SZ maps. The SPT result is from the exact same data as A22, but analysed using the slightly updated measurement pipeline described in Section 3.1. As was the case in A22, the theoretical prediction matches the measurements in the cluster core ($R/R_{200m} \lesssim 0.5$) and also in the far outskirts ($R/R_{200m} \gtrsim 5$), but has significant deviations at $R/R_{200m} \approx 1$, and potentially also at $R/R_{200m} \approx 3$. These two deviations were denoted a pressure deficit and accretion shock, respectively, in A22 and we use the same nomenclature here.

This pressure deficit was discussed in A22 as a possible sign of thermal non-equilibrium between electrons and ions, where the non-equilibrium is generically caused by shock heating (Fox & Loeb 1997; Ettori & Fabian 1998; Wong & Sarazin 2009; Rudd & Nagai 2009; Akahori & Yoshikawa 2010; Avestruz et al. 2015; Vink et al. 2015). Shocks are the primary mechanism for converting kinetic energy to thermal energy during structure formation. They preferentially heat the ions over the electrons given the former are more massive. Thus, shock-heated plasma has colder electrons than protons, and the low density of particles in the cluster outskirts implies these two particle species never equilibrate. Rudd & Nagai (2009, see their fig. 2) use simulations specialized to model the electron-ion temperature differences and show that this effect causes a deficit in the cluster tSZ profiles,<fnref n="12" /> while Avestruz et al. (2015, see their fig. 1) do the same but focus on the 3D cluster temperature profiles. This pressure deficit feature would not be present in most

<footnote n="11">¹¹SZ-selected clusters also incur a noise-induced miscentring effect, with a scale of $\theta_{\text{miscen}} = \sqrt{\theta_{500c}^2 + \theta_{\text{beam}}^2}/\text{tt SNR}$. For $R \gtrsim R_{200m}$, the miscentring scale is at/below the bin width and is negligible as our features of interest span multiple bins. The average SZ-selected cluster ($M_{200m} \approx 10^{14.8} M_{\odot}$ and $z \approx 0.6$) has $\theta_{\text{miscen}} = 0.3'$, while the same for the average optically selected cluster ($M_{200m} \approx 10^{14.6} M_{\odot}$ and $z \approx 0.4$) is factors of 5–0 larger (Zhang et al. 2019a; Bleem et al. 2020).</footnote>

<footnote n="12">¹²Such a deficit should also be present in electron temperature profiles measured through X-ray data. However, our current X-ray observations do not extend to such large radii, $R \approx R_{200m}$, and are instead limited to much smaller radii where the higher number densities allow the ion and electrons to quickly achieve temperature equilibrium.</footnote>
<table><caption>**Table 1.** A summary of the numerical results presented in this work. All uncertainties are ±1σ estimates. From left to right the columns show: (i) the sample name, denoted as 'cluster catalogue source' x 'SZ map source', (ii) location of the pressure deficit, (iii) the value of the log-derivative at the location, (iv) detection significance of the feature, extracted using equations (22) and (24), (v–vi) the weighted mean of the log-mass and the redshift of the cluster sample (using cluster SNR as weights), (vii) the number of clusters in the sample, and (viii) the Figure in this work containing the profile corresponding to the result. We do not quote a detection significance for the optically selected clusters given the dependence of this significance on the assumed miscentring model (see Section 4.1 and Appendix A for details). The uncertainties are estimated via jackknife resampling (see Section 3.1) and do not include systematic uncertainties.</caption><thead><tr><th>Data set</th><th>R<sub>sh</sub>/R<sub>200m</sub></th><th>d ln y / d ln R (R<sub>sh</sub> / R<sub>200m</sub>)</th><th>χ<sub>sh</sub></th><th>⟨log<sub>10</sub> M<sub>200m</sub>⟩ [M<sub>☉</sub>]</th><th>⟨z⟩</th><th>N<sub>cl</sub></th><th>Figure</th></tr></thead><tbody><tr><td>SPT x SPT</td><td>1.09 ± 0.08</td><td>-3.98 ± 0.48</td><td>2.7σ</td><td>14.94</td><td>0.57</td><td>503</td><td>2</td></tr><tr><td>ACT x ACT</td><td>1.16 ± 0.04</td><td>-3.53 ± 0.12</td><td>6.1σ</td><td>14.84</td><td>0.55</td><td>4045</td><td>2</td></tr><tr><td>DES x SPT</td><td>0.95 ± 0.09</td><td>-2.71 ± 0.2</td><td>-</td><td>14.63</td><td>0.44</td><td>1990</td><td>2</td></tr><tr><td>DES x ACT</td><td>1.14 ± 0.07</td><td>-3.11 ± 0.15</td><td>-</td><td>14.61</td><td>0.44</td><td>4340</td><td>2</td></tr><tr><td>ACT x SPT (overlap)</td><td>1.15 ± 0.07</td><td>-3.79 ± 0.49</td><td>2.6σ</td><td>14.83</td><td>0.58</td><td>669</td><td>3</td></tr><tr><td>ACT x ACT (overlap)</td><td>1.15 ± 0.05</td><td>-3.77 ± 0.44</td><td>2.7σ</td><td>14.83</td><td>0.58</td><td>669</td><td>3</td></tr><tr><td>ACT x ACT (sel. effect)</td><td>1.11 ± 0.06</td><td>-3.26 ± 0.18</td><td>4.5σ</td><td>14.84</td><td>0.46</td><td>3297</td><td>4</td></tr><tr><td>DES x ACT (sel. effect)</td><td>1.14 ± 0.03</td><td>-3.17 ± 0.17</td><td>-</td><td>14.83</td><td>0.46</td><td>4034</td><td>4</td></tr><tr><td>DES x SPT (high M)</td><td>0.95 ± 0.09</td><td>-2.71 ± 0.2</td><td>-</td><td>14.63</td><td>0.44</td><td>1990</td><td>5</td></tr><tr><td>DES x SPT (med M)</td><td>0.83 ± 0.19</td><td>-1.87 ± 0.16</td><td>-</td><td>14.19</td><td>0.51</td><td>20712</td><td>5</td></tr><tr><td>DES x SPT (low M)</td><td>1.44 ± 0.1</td><td>-2.76 ± 0.88</td><td>-</td><td>13.91</td><td>0.55</td><td>20973</td><td>5</td></tr><tr><td>DES x ACT (high M)</td><td>1.14 ± 0.07</td><td>-3.11 ± 0.15</td><td>-</td><td>14.61</td><td>0.44</td><td>4340</td><td>5</td></tr><tr><td>DES x ACT (med M)</td><td>1.10 ± 0.12</td><td>-2.23 ± 0.14</td><td>-</td><td>14.21</td><td>0.50</td><td>45851</td><td>5</td></tr><tr><td>DES x ACT (low M)</td><td>1.27 ± 0.45</td><td>-1.82 ± 0.36</td><td>-</td><td>13.88</td><td>0.55</td><td>47426</td><td>5</td></tr><tr><td>ACT x ACT (high M)</td><td>1.13 ± 0.08</td><td>-3.39 ± 0.19</td><td>3.8σ</td><td>14.97</td><td>0.51</td><td>1635</td><td>5</td></tr><tr><td>ACT x ACT (med M)</td><td>1.15 ± 0.08</td><td>-3.56 ± 0.33</td><td>2.5σ</td><td>14.75</td><td>0.57</td><td>1183</td><td>5</td></tr><tr><td>ACT x ACT (low M)</td><td>1.23 ± 0.15</td><td>-4.15 ± 0.68</td><td>2.9σ</td><td>14.64</td><td>0.62</td><td>1217</td><td>5</td></tr><tr><td>ACT x ACT (R<sub>sp</sub> comparison)</td><td>1.00 ± 0.17</td><td>-3.19 ± 0.20</td><td>4.3σ</td><td>14.86</td><td>0.45</td><td>1138</td><td>8</td></tr><tr><td>SPT x SPT (high SNR)</td><td>1.11 ± 0.04</td><td>-4.27 ± 0.66</td><td>2.4σ</td><td>15.02</td><td>0.56</td><td>259</td><td>B1</td></tr><tr><td>SPT x SPT (low SNR)</td><td>0.97 ± 0.15</td><td>-3.68 ± 0.94</td><td>1.4σ</td><td>14.81</td><td>0.58</td><td>272</td><td>B1</td></tr><tr><td>ACT x ACT (high SNR)</td><td>1.18 ± 0.08</td><td>-3.41 ± 0.23</td><td>2.5σ</td><td>14.97</td><td>0.55</td><td>1401</td><td>B1</td></tr><tr><td>ACT x ACT (med SNR)</td><td>1.19 ± 0.05</td><td>-3.65 ± 0.38</td><td>2.5σ</td><td>14.75</td><td>0.57</td><td>1394</td><td>B1</td></tr><tr><td>ACT x ACT (low SNR)</td><td>1.13 ± 0.06</td><td>-4.01 ± 0.68</td><td>4.0σ</td><td>14.69</td><td>0.53</td><td>1400</td><td>B1</td></tr></tbody></table>

cosmological hydrodynamical simulations as they a priori assume local thermal equilibrium between electrons and ions. We will henceforth refer to the pressure deficit as a shock feature and denote its location the shock radius, $R_{sh}$.

As Fig. 2 and Table 1 show the ACT DR6 data strengthen the evidence for a pressure deficit feature near the cluster virial radius. This is the same feature first noted in A22 with SPT-SZ data and with ACT DR5 clusters measured on the ACT DR4 map. We estimate the significance of the feature in the ACT data at $6.1\sigma$. Given the new, more robust definition of SNR in equation (24) and the switch from the ‘Shock heating’ model of Battaglia et al. (2012) to the ‘200 AGN’ model, the estimated significance of the feature in SPT-SZ is $2.7\sigma$ compared to the estimate of $3.1\sigma$ from A22. We have verified that our pipeline reproduces the SPT-SZ result of the previous work if we revert back to the previous signal-to-noise definition and model choice.

The deficit in both SPT and ACT is found at consistent radial locations, with $R/R_{200m} = 1.09 \pm 0.08$ and $R/R_{200m} = 1.16 \pm 0.04$, respectively. The minima in the log-derivatives are consistent as well, with $\frac{d\ln y}{d\ln R} = -4 \pm 0.5$ and $\frac{d\ln y}{d\ln R} = -3.5 \pm 0.1$, respectively. These estimates are detailed further in Table 1. The similarity of the deficit seen in SPT and ACT suggests the feature is physical and not an artefact introduced in either the map-making or the cluster-finding procedures in each survey. We have also independently verified the consistency of these features using a complementary fitting method, described in Appendix C. In A22, we validated that the theoretical model used in this work matches cosmological hydrodynamical

simulations (see their fig. 4). In specific, we used the THE300 suite which simulates a sizable number of massive clusters, and provides a sample relevant for SZ-selected cluster catalogues which have $M_{200m} > 10^{14.5} M_{☉}$. Thus, any differences between the measurements and the theoretical profiles can be equivalently interpreted as differences between the measurements and simulations.

The bottom panels of Fig. 2 also present the quantity $\epsilon$, defined in equation (21), which is the bin-by-bin deviation between the measured and predicted log-derivatives, normalized by the measurement uncertainty. In SZ-selected clusters, $\epsilon$ takes a maximum value at $R/R_{200m} \approx 1$, corresponding to the pressure deficit. In optically selected clusters, which we will discuss below, the maximum values of $\epsilon$ are at smaller scales. This is because the measurement is much more precise on these scales so small deviations between the data and theory – such as those caused by imperfections in the miscentring model – can have large statistical significance.

We do not discuss the potential accretion shock features in detail as these are currently still low-significance features dominated by noise, as was the case in A22. We simply note it is intriguing that the log-derivatives of the SPT and ACT profile measurements both have a maximum at $R/R_{200m} \approx 3$, followed by a sharp drop. The maximum corresponds to a plateauing phase in the profiles, which is a feature of the accretion shock as presented in Baxter et al. (2021). More detailed work is required to robustly verify this feature as arising from the presence of a shock.

Other studies also find features in the cluster outskirts using a variety of different data sets. Hurier, Adam & Keshet (2019) see a
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The average SZ profiles of different cluster samples, and measured using different SZ maps (top), their associated log-derivative (middle), and the difference between the log-derivatives of the data and model (bottom) as defined in equation (21). The theoretical prediction (dashed lines), which is a sum of one-halo and two-halo contributions (each shown as grey dotted and dash-dotted lines, respectively, in the left panel for the SPT predictions alone), is described in Section 3.2. The left panels show results for SZ-selected clusters from SPT and ACT, while the right is optically selected clusters in DES with a mass cut $M_{200m} > 10^{14.5}$ M⊙. For the SZ-selected samples, the derivative is lower at $R_{200m}$ than the theory curve, consistent with A22. The behaviour for optically selected clusters is less clear due to potential inaccuracies in the theoretical model, such as the miscentring model. All profile measurements have a two-halo component, seen most prominently at large radii, that is consistent with the model. Estimates for the location and depth of the first log-derivative minimum ('pressure deficit') in each measurement are shown in Table 1. The grey band in the left panels demarcates the range of radii used to quantify the significance of the pressure deficit as shown in the table. The dotted lines in the right panel are the theory models without any miscentring effects included; including the miscentring (dashed line) changes/improves the model. The two dashed lines in the right panels overlap with one another. The correlation matrix of the log-derivative is shown in Fig. D1.</figcaption></figure>

sharp decrease in pressure at $R = 3R_{500c} \approx R_{200m}$ for a single cluster in the Planck data. Pratt, Qu & Bregman (2021) also use Planck data and find an excess in pressure at $R = 2R_{500c} \approx 0.7R_{200m}$ for a set of ten, low-redshift galaxy groups. The analysis of Planck Collaboration et al. (2013) finds that the 3D pressure profiles have a deficit, relative to the theoretical predictions of Battaglia et al. (2012), for $R \gtrsim R_{500c}$ ($R \gtrsim 0.3R_{200m}$) while being a good match for scales below that radius. Zhu et al. (2021) find an excess in the temperature and density profiles of the Perseus cluster at $R \approx R_{200c} = 0.5R_{200m}$ using Suzaku X-ray data. Hou, Hallinan & Keshet (2023) study the radio emission around galaxy clusters and find a signal at $R = 2.5R_{500c} \approx R_{200m}$. They interpret this as the presence of a non-thermal electron population and find that the corresponding electron energy distribution is consistent with one generated by strong shocks. In all works, the deviations are found around $R \approx R_{200m}$, consistent with the shock radius $R_{sh}$.

The right panels of Fig. 2 show, for the first time, the outskirts of SZ profiles for optically selected clusters. We have placed a mass cut of $M_{200m} > 10^{14.5}$ M⊙ (where $M_{200m}$ is the mass inferred from the cluster richness, see Section 2.1) on this sample as this reduces

the impact of systematic effects (such as projection, contamination, etc.); this cut is also consistent with the minimum mass of the SZ-selected samples (see Fig. 1). We discuss the results of lower mass objects, which are removed by this cut, in Section 4.3. The SZ profiles of DES clusters have a $\approx 30$ per cent lower normalization than those of the SZ-selected catalogues, and this is due to the differences in mass distributions and the mean mass of the samples (see Table 1). The normalization of the theoretical model (dashed lines) also decreases a similar amount if we input the DES cluster mass/redshift distribution rather than the SPT or ACT ones. At $R_{200m}$, which is inbetween the one-halo and two-halo regime, the profile for the DESxACT measurement has a minimum log-derivative ($\frac{d\ln v}{d\ln R} = -3.1 \pm 0.15$) that is more negative than that of the DESxSPT measurement ($\frac{d\ln v}{d\ln R} = -2.7 \pm 0.2$), with a significance of $1.6\sigma$. The two results use different cluster subsamples, defined as all DES clusters within the ACT/SPT footprint. We interpret this difference as a statistical variation and do not examine it further. We verify in Fig. 3 that the SPT and ACT maps provide statistically indistinguishable results across the full range of scales considered in this work.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The average SZ profile of an ACT cluster subsample ($N = 669$ clusters) measured using either the ACT map or SPT map. The subsample is defined as all clusters whose centres lie in both the ACT and SPT footprints. The two measurements are consistent across the whole range of scales, with $\chi^2/N_{\text{dof}} = 1.1$ and $p = 0.14$, validating that the data sets and map-making procedures of the two surveys are consistent in both high and low signal-to-noise regimes. The SPT and ACT data sets are independently calibrated and mapped, and the statistical consistency in the measurements above is determined at the $\approx 1$ per cent level given the precise measurements in the high signal-to-noise regime.</figcaption></figure>

Looking at the DES x ACT and the ACT x ACT results, we see the location of the log-derivative minima is consistent at $0.2\sigma$, while the depth of minima is deeper in ACT x ACT at $2\sigma$. The comparison of DES x SPT and SPT x SPT is similar, where the location of the log-derivative minima is consistent while the depth deviates at $2.4\sigma$. The mass and redshift distributions of the DES cluster sample are notably different from those of ACT and SPT, which could lead to differences in this depth. In Section 4.2, we re-analyse the ACT and DES data after accounting for such mass/redshift differences, and find that the depth becomes consistent across the two measurements.

The model (dashed line) for the DES-related results in Fig. 2 is a qualitatively good match to the data across the whole range of presented scales. The prediction for the DES x SPT and DES x ACT measurements closely overlap one another. This model includes the miscentring effects described in Section 3.2.1, using values of $\tau_{\text{miscen}} = 0.9$ and $f_{\text{miscen}} = 0.4$. These values were chosen after exploring a sparsely sampled 2D grid of parameter values and picking the parameters that provided the visually best fit to the one-halo regime, near the cluster core. The preferred values for $\tau_{\text{miscen}}$ and $f_{\text{miscen}}$ are both near the $3-4\sigma$ upper limit of the miscentring parameter constraints of Zhang et al. (2019a, see their Chandra-DES constraints in Table 1) for the DES Y1 cluster sample. However, the value of $\tau_{\text{miscen}}$ is within $1\sigma$ of the estimate from Bleem et al. (2020, see their Table 6), which is based on a SPT-DES matched cluster sample. Fig. 2 shows the theory matches the data better (in the 1-halo regime) when we include this miscentring effect, and the dotted lines show the theory without such effects.

In Appendix A, we discuss how the profiles and log-derivatives depend on miscentring parameters. We emphasize that in our work we only focus on results from optically selected clusters that are insensitive to the choice of miscentring model and parameters. For

example, Table 1 does not quote any detection significance of a pressure deficit for DES clusters. However, we still measure and quote the location and depth of the log-derivative minimum for the DES cluster profiles as it does not depend on an assumed theoretical model.

Our results show that the SZ-selected clusters have a clear pressure deficit while such a deficit is not seen as clearly in optically selected clusters. In general, this difference could occur if (i) SZ-selected clusters have a selection effect that preferentially picks out objects with such features, (ii) an aspect of the richness-SZ-mass correlations makes optically selected clusters suppress the deficit feature, and (iii) systematic effect(s) in optically selected clusters (e.g. the miscentring, contamination, or mass estimation errors) causes the feature to be suppressed. In Section 4.2, we verify that the first two possibilities are not the cause for the difference between the results of SZ-selected and optically selected clusters. The third possibility – the systematic effects in optically selected clusters – is an intricate issue spanning many different parts of the cluster detection/processing pipeline, and so we do not explore this direction as it is beyond the scope of our work. However, in Section 4.2, we will show that limiting the DES clusters to higher masses, $\langle M_{200m} \rangle = 10^{14.85} M_{\odot}$, results in the profile measurement showing a deficit that is consistent with those of the SPT and ACT clusters. This in turn implies that the three effects mentioned above have negligible impact on the measurements if we use optically selected clusters that are limited to higher masses than those of the fiducial sample used in Fig. 2.

One SZ-related systematic effect is the CIB, which is sourced by dusty, star-forming galaxies. DES clusters are selected on richness (i.e. galaxy counts) and preferentially contain clusters with more satellite galaxies compared to an SZ-selected sample. Thus, the amplitude of the infrared signal for such a sample could be higher. However, SZ-selected samples probe higher redshifts than optically selected clusters, which are closer to the peak of cosmic star-formation at $z = 2$. We verify in Appendix E that our results are unchanged if we use SZ maps that minimize/deproject the CIB signal.

## 4.2 Sensitivity to map-making and cluster selection

The results in Fig. 2 show a difference between the profiles of SZ-selected and optically selected clusters – the former sees a clear pressure deficit at $R/R_{200m} \approx 1$, while the latter either sees a less significant feature or no feature at all – and this could be caused by SZ-selection preferentially picking out clusters with such deficit-like features, or by the optical selection effects preferentially missing such clusters (in this case, due to a correlation this feature may have with cluster richness).

In Fig. 3, we test an aspect of the former effect, namely noise-based SZ-selection effects.<fnref n="13" /> These effects correspond to the fact that the clusters are identified in the same (noisy) maps used to measure their SZ profiles. We test the impact of this effect by taking all ACT clusters that fall into the intersection of the ACT and SPT footprints ($N = 669$ clusters), and then by measuring the subsample’s average SZ profile using either the SPT map or the ACT map. We find consistency $(\chi^2/N_{\text{dof}} = 1.1$ with $p = 0.14)$ in both the profiles and the log-derivatives of the two measurements. While this implies

<footnote n="13">¹³We consider this a systematics-based selection effect, in contrast to physical selection effects such as, for example, SZ-selected clusters being preferentially more/less dynamically active compared to mass-selected clusters. An aspect of these physical SZ-selection effects is tested in Fig. 4.</footnote>
that noise-based SZ-selection effects are not the cause of the pressure deficit feature, the agreement is also a check on the data and map-making procedures of the SPT and ACT surveys.<fnref n="14" /> It validates the maps' consistency in both the high SNR regime at the location of massive clusters, as well as in the noise-dominated, low SNR regime of the cluster outskirts.

Next, we test the impact of optical selection on this deficit feature by comparing profiles around ACT and DES clusters that are reweighted to have the same mass/redshift distribution. The reweighting is done to minimize any differences in the measured average SZ profiles due to differences in just the mass/redshift distribution of the samples.<fnref n="15" /> We first remove all ACT clusters with redshifts/masses outside the ranges of the DES sample. We therefore use all ACT clusters within 0.1 < z < 0.8 and 13.7 < log₁₀ M₂₀₀m < 15.35 to create the subsample used in this analysis, which has N = 3392 clusters. The reweighting is then done by computing the weighted counts of clusters in a 2D grid of M₂₀₀m and z, and then using the ratio of ACT counts to DES counts. The weight used in the weighted counts is the signal-to-noise per cluster, consistent with the rest of our analysis. The exact expression of the re-weighting is

$$w(M_{200m}, z) = \frac{\sum_{i=1}^{N_{\text{cl,ACT}}} \delta_i \text{SNR}_i^{\text{ACT}}}{\sum_{i=1}^{N_{\text{cl,DES}}} \delta_i \text{SNR}_i^{\text{DES}}}, \quad (30)$$

where $\delta_i$ is a delta function – with values of 0 or 1 – that denotes whether cluster i falls into a given mass and redshift bin. We compute the weights in a 10-by-10 grid and assign each DES cluster a new weight based on the M₂₀₀m – z grid cell it is associated with. We have checked that our results do not change if we use a 20 x 20 grid instead. The final weight of the DES cluster is

$$w^{\text{mod}}(M_{200m}, z) = \text{SNR} \times w(M_{200m}, z), \quad (31)$$

which uses the original SNR weights of our analysis alongside the mass/redshift-based reweighting of equation (30). We have tested that our results, shown below, are unchanged if we exclude all ACT clusters in the DES footprint, where this exclusion would remove any overlap between the cluster samples.

Fig. 4 shows the average SZ profile around the ACT subsample and the reweighted DES sample. The two profiles are consistent with one another. The ACT subsample shows a clear pressure deficit in the log-derivatives – evidenced by the measured profile dropping more steeply at R/R₂₀₀m ≈ 1 than the theoretical prediction – and the DES measurement matches this feature. This consistency is partially expected as the DES reweighting increases the contribution of the most massive clusters to the average SZ profile, and any systematic effects on the pressure deficit measurement could be less prominent in this mass regime. However, it is still a valuable check as even in the high-mass regime, optically selected clusters have shown differences in their total matter density profiles that were due to optical selection effects (Baxter et al. 2017; Chang et al. 2018; Shin et al. 2019).

Fig. 4 provides evidence that at high mass, the optical selection does not result in biased SZ profiles for the one-halo and two-halo regimes. Across the range 0.5 < R/R₂₀₀m < 10, the two profiles are

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The average SZ profiles of two different cluster samples: an SZ-selected one from ACT, and an optically selected one from DES. In both cases, the samples are modified from the original distribution. We use all ACT clusters with 0.1 < z < 0.8, and then reweight the DES clusters to match the ACT subsample's M₂₀₀m – z distribution (see Section 4.2 for details). The two profiles are consistent ($\chi^2/N_{\text{dof}} = 1.1$ with $p = 0.14$), suggesting that in this mass/redshift range there are no SZ or optical selection effects that generate/suppress the pressure deficit feature.</figcaption></figure>

consistent ($\chi^2/N_{\text{dof}} = 1.4$ with $p = 0.1$). Under this reweighting, the weighted mean mass of the DES sample increases from $\langle M_{200m} \rangle = 10^{14.6} M_{\odot} \rightarrow 10^{14.85} M_{\odot}$, which is a fractional change of 80 per cent, while the mean redshift is left unchanged at $\langle z \rangle = 0.46$ (see Table 1). The depth of the minima is now consistent across the two samples, whereas it was inconsistent at the 2$\sigma$ level for the fiducial ACT and DES cluster samples (Fig. 2). Given that agreement between the samples is recovered after accounting for their mass/redshift differences, we infer that the earlier disagreement was due to these differences.

The results of Figs 3 and 4 imply that – for a mass and redshift range corresponding to clusters in SZ surveys (see Fig. 1) – the SZ or optical selection has negligible impact on the measured pressure deficit. This adds to the robustness of the deficit features found in the SPT x SPT and ACT x ACT measurements, as clusters identified with a completely different type of data (i.e. optical images) still show a pressure deficit. These results also show that the DES sample exhibits a clear deficit (given its agreement with the ACT measurement) when limited to higher masses, implying that the shallower log-derivative depth found in our fiducial measurement (Fig. 2) could possibly be attributed to the clusters in the lower mass end of the sample. We explore the behaviour of such systems further in the following section.

## 4.3 Towards galaxy groups

The profiles of massive clusters have many observational constraints, especially near the cluster core ($R/R_{200m} \lesssim 0.5$, see for example, McDonald et al. 2014; Ghirardini et al. 2017; Romero et al. 2017, 2018; Ghirardini et al. 2018), and the further outskirts have only been recently explored observationally (Planck Collaboration V 2013; Sayers et al. 2013, 2016; Amodeo et al. 2021; Schaan et al. 2021;

<footnote n="14">¹⁴A similar analysis using all SPT clusters in both footprints finds $\chi^2/N_{\text{dof}} = 1.06$ with $p = 0.36$. However, the profile measurement uncertainties are broader as the SPT cluster sample size is half that of ACT.</footnote>

<footnote n="15">¹⁵The zeroth-order effect of SZ and optical selection on the cluster sample is in its mass and redshift distributions (see Fig. 1). The reweighting accounts for these selection effects, and thus any further differences in the reweighted profiles can be attributed to selection effects beyond those on the cluster samples' mass and redshift distributions.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The average SZ profiles, binned according to inferred halo mass, for different combinations of cluster samples and SZ maps. The dashed lines are the theoretical prediction described in Section 3.2. Note that the mass ranges for the ACT sample (right) are significantly narrower than those of the DES sample (left, middle). The measured profiles of lower mass clusters ($M_{200m} < 10^{14.5} M_{\odot}$) deviate significantly from the theoretical predictions, but the two-halo term is still consistent between data and theory. This two-halo term is prominent at large radii, as shown in Fig. 2.</figcaption></figure>

Melin & Pratt 2023; Lyskova et al. 2023). Less massive objects – galaxy group scales and down to Milky Way scales – have not been studied on a profile level, and the presence/absence of any features in the outskirts is relatively unknown. Previous works have studied the cross-correlation function of the tSZ field with galaxy counts (Hill et al. 2018; Amodeo et al. 2021; Schaan et al. 2021; Sánchez et al. 2023), which is an observable that is sensitive to halo profiles but cannot always distinguish features in the profiles. For example, the pressure deficit in the cluster outskirts (Fig. 2) was not identified in previous cross-correlation works but was easily identified in A22 by measuring individual profiles. SZ-selected halo samples are ideal for studying massive, cluster-scale haloes but are not viable for probing lower masses. Here, we use the DES REDMAPPER sample to obtain a catalogue of lower mass objects ($M_{200m} \gtrsim 10^{13.8} M_{\odot}$) and measure their SZ profiles across a wide range of scales.

In Fig. 5, we show the average SZ profile for three mass bins of DES clusters, measured on both the SPT map (left) and the ACT map (middle), and also the profiles for the ACT cluster sample measured on the ACT map (right). The mass bins of the latter differs significantly from those of the former two. Focusing first on the ACT cluster results in the rightmost column of Fig. 5, we see the pressure deficit exists for all mass bins, at $3.8\sigma$, $2.5\sigma$, and $2.9\sigma$ from highest to lowest mass bins. The three minima from the log-derivative measurements are all statistically consistent with each other. This result also serves as an additional validation check – if the pressure deficit in SZ-selected samples is caused by noise-based selection effects, then its amplitude will be higher for clusters detected at the low signal-to-noise regime, which is right near the cluster detection threshold. In Fig. 5, however, we find that splitting by mass – which is directly proportional to SNR – does not notably change the significance of the deficit. In Appendix B, we also redo this test by splitting directly on SNR instead of $M_{200m}$, and find consistent results.

The other two columns in Fig. 5 (left and middle) show the SZ profile for three mass bins of optically selected clusters. The

mass range $M_{200m} > 10^{14.5} M_{\odot}$ corresponds to $\lambda > 30$, while the range $10^{14} M_{\odot} < M_{200m} < 10^{14.5} M_{\odot}$ corresponds to $15 < \lambda < 30$, and finally, $10^{13.8} M_{\odot} < M_{200m} < 10^{14} M_{\odot}$ corresponds to $10 < \lambda < 15$. These masses are not exact translations of the richness but rather approximate conversions for interpreting the discussions to follow. The highest mass bin (purple) is the same result as Fig. 2 and shows good, qualitative agreement between the measurements and the theoretical predictions. When comparing the measurements of lower mass bins to those of the highest mass bin, we see that for lower mass objects the log-derivatives are closer to zero in the one-halo term ($R/R_{200m} \lesssim 1$) and similar to the high mass bin results for the two-halo term ($R/R_{200m} \gtrsim 4$). The log-derivative minima in each mass bin are found at similar radii of $R/R_{200m} \approx 1.1$ (see Table 1).

The theoretical model also significantly deviates from the measurements in these two lower mass bins. For $M_{200m} \in [10^{13.8}, 10^{14}] M_{\odot}$, the deviation is a factor of $\approx 5$ in the halo core. In the intermediate mass bin, $M_{200m} \in [10^{14}, 10^{14.5}] M_{\odot}$, it is a factor of $\approx 2$ and is also consistent with previous analyses in this intermediate mass range. Saro et al. (2017, see their table 1 and figs 5/6) found a factor of $\approx 2$ difference when measuring the integrated SZ effect around clusters from the DES Science Verification data, while Planck Collaboration XII (2011, see their fig. 2) finds similar suppression in the SZ-richness scaling relation. These differences generically point to some inaccuracy in the theoretical model.

The discussion of the mass trends thus far focuses on the behaviour of the log-derivative minima, rather than of the 'pressure deficit'. The latter is defined as significant deviations between measurement and model in the shape of the profile. However, for the lower mass bins, the model has inaccuracies as noted above, which limit our ability to identify such a feature. There are a few known reasons why such inaccuracies could occur: (i) the contamination of the cluster sample at low masses, e.g. two or more low-mass clusters are projected together on the sky and are observed as one large cluster, causing a mass estimation bias (ii) inaccuracy in the utilized pressure profile model for lower mass haloes, and; (iii) significant correlations in
the richness and SZ scatter at fixed halo mass which, in tandem with the optical selection effect at low richness, could become an important effect. We briefly discuss each to check if it can explain the deviations and thereby provide an avenue to correct the existing model prediction.

**The first, contamination of the sample**, causes an overestimate of the cluster mass (and thus, the SZ profile) compared to the truth. This overestimate is more significant in the one-halo regime than the two-halo regime as the latter's mass-dependence is weaker. The two-halo term scales as $y \propto b_h(M) \propto M^{0.5}$ for the halo bias model of Tinker et al. (2010) at high halo masses, whereas the one-halo term scales as $y \propto M^{5/3}$. The deviations in Fig. 5 are roughly factors of 2–5 in the SZ signal, and suggest the corresponding maximum bias in the mass – assuming a self-similar scaling of $y \propto M_{200m}^{5/3}$ – would be 50 per cent to 150 per cent in $M_{200m}$. Myles et al. (2021) show the richness bias due to contamination is $\approx 20$ per cent for clusters of $5 < \lambda < 20$ (see their Section 4.3), and also that richness depends on halo mass as $M \propto \lambda^{1.0}$ (see their section 4.5). This implies the mass bias is $20$ per cent $\times 1.0 = 20$ per cent, lower than the required values of $50$ per cent to $150$ per cent denoted above, and provides evidence that contamination from projection cannot be the dominant cause of the suppression. Similarly, variations in the assumed projection model of the mass–richness relation show $\approx 30$ per cent changes in the final mass estimate (Costanzi et al. 2021, see their equations 16 and 17).

**The second effect, Y – M relation deviations**, are deviations in the pressure profile model for lower mass haloes. This work uses the model of Battaglia et al. (2012), and while it is accurate for higher mass haloes (e.g. see Fig. 2), observational analyses find a preference for deviations from this model at lower halo masses (Hill et al. 2018; Pandey et al. 2022). Such deviations can arise from differences between the assumed galaxy formation process in the simulations, and the relevant processes in the data. In particular, these above works suggest the SZ signal for the lower mass bins we consider here is suppressed by factors of 3–4 and that the suppression grows stronger with decreases in halo mass. Both these behaviours are consistent with our findings. However, the uncertainties on the inferred suppression are not precise enough to confirm that this effect is the dominant cause of the deviations in Fig. 5.

**The third effect, correlations in the richness and SZ scatter at fixed mass**, is relevant as our work involves the simultaneous use of cluster mass, SZ, and richness; we select clusters using richness, infer a halo mass from this richness, and then use the inferred halo mass to predict the SZ profile. The correlations between the three properties require non-trivial corrections to the model for the SZ–mass scaling relation of the selected cluster sample. The effect has been detailed in the analytical work of Evrard et al. (2014, see their fig. 4 for an example). The scaling relation for the optically selected sample is now written as

$$ \langle \ln y | \ln \lambda \rangle = \langle \ln y | \ln M_{200m} (\ln \lambda) \rangle + \beta_1 \alpha_{\lambda} \times \operatorname{cov}(\ln y, \ln \lambda | \ln M_{200m}(\lambda)), \quad (32) $$

where $-\beta_1$ is the slope of the halo mass function at a chosen mass scale, $\alpha_\lambda$ is the slope of the richness–mass relation, and $\operatorname{cov}(\dots)$ is the covariance of the SZ and richness scatter. A general form of this expression can be found in Evrard et al. (2014, see their Equation 6). Inspecting equation (32) shows $\langle \ln y | \ln \lambda \rangle$ can be higher (lower) than $\langle \ln y | \ln M_{200m} (\ln \lambda) \rangle$, for a positive (negative) sign of correlation in the SZ–richness scatter at fixed mass. Farahi et al. (2019) observationally constrain this correlation coefficient to be $-0.5 \lesssim r \lesssim 0.5$ at 95 per cent confidence, and their results indicate

the correlation of gas-based and stellar-based cluster observables is negative (see their table 2). Cosmological simulations also show the correlation is negative – the scatter of gas mass and stellar mass are anti-correlated (Farahi et al. 2018, see their fig. 5) while that of the stellar mass and richness are correlated (Anbajagane et al. 2020, see their fig. 7). A negative correlation/covariance suppresses the SZ signal of richness-selected clusters, which could cause the observed suppression. For conservative values of $-\beta_1 = 1.5$, $\alpha_\lambda = 1.5$, $\sigma_{\ln y} = 0.3$, $\sigma_{\ln \lambda} = 0.8$, $r = -0.6$,<fnref n="16" /> we find the bias is at most 30 per cent. Thus, this effect cannot be the main cause of the behaviours found in Fig. 5.

Our discussions and estimates above indicate the deviations between measurement and model are unlikely to be explained by just one of these effects. Thus, the model we use for SZ profiles cannot be easily corrected to match our measurements in the one-halo regime of lower mass clusters. Furthermore, the latter two effects we discuss – deviations in the Y–M relation and the correlated richness and SZ scatter – are also functions of radius that are not well-known and would be required to accurately correct our model. Previous works (including all works cited above) have only discussed these effects for volume-integrated quantities, rather than for radial profiles. Accurate predictions for these profiles, however, are necessary to study a pressure deficit (i.e. shock-induced deviations between the data and model). Given this limitation, our main results of this section focus on the raw log-derivative measurements (rather than inferring a pressure deficit from them by comparing to theory), which have a clear striation with mass in the one-halo regime and weak-to-no striation in the two-halo regime. The minima of these derivatives are located at similar radii for all three mass bins.

# 5 CONNECTIONS TO STRUCTURE FORMATION FEATURES

Having explored the average SZ profiles using different combinations of cluster samples and SZ maps, we now connect these profiles to broader features from structure formation. First, we detail the connection to cosmic filaments in Section 5.1 via oriented stacking of the profiles. Then in Section 5.2, we compare the pressure deficit seen in the SZ profile to the splashback feature observed in the galaxy number density profile measured around clusters.

As mentioned prior, discussing shock-induced features also requires discussing behaviours in the one-halo and two-halo regimes as shocks occur at the transition between the two. Therefore, some of our discussions below include the behaviours of these two regimes.

## 5.1 Connections to filaments

We have discussed previously that cosmological shocks form from the accretion of collisional matter on to bound objects, and the accreted matter originates primarily from cosmic filaments. Simulations suggest that the shock boundary generally follows the same ellipticity/orientation as the cluster’s, which in turn is informed by the filaments’ topology around the cluster (Aung, Nagai & Lau 2021, see their fig. 1). However, along the specific line of sight connecting the cluster core and the filament, the accretion rate of cold gas (cold relative to the hot gas bound in the cluster) is highest and can push the

<footnote n="16">¹⁶$\beta_1$ is the slope at a pivot mass of $M_{200m} \sim 10^{14} M_\odot$, computed using the halo mass function of Tinker et al. (2008), $\alpha_\lambda$, $\sigma_{\ln y}$, and $\sigma_{\ln \lambda}$, are chosen to be larger than constraints from Costanzi et al. (2021, see their table 4) and $r$ is set by the 95 per cent bound from Farahi et al. (2019).</footnote>
shock feature further into the cluster core and/or completely destroy it (Zhang et al. 2020, see their fig. 6).

In our analysis, we use various orientation measures, each probing a different range of scales, as estimates of the orientation of the nearby filamentary structure around the cluster. We then split the 2D SZ image of each cluster into three equal-area sections – according to how close a section is to the major axis of the orientation – and compute the profile using pixels within each sub-section of the image. The geometry of this split is shown in Fig. 6. In A22, we split the cluster into two equal areas, corresponding to the major and minor axis. In this work, we add a third area that probes the intermediate region. Through this, we can more easily distinguish coherent trends across the orientations from any noise fluctuations. This increase in subsections is made possible by our larger cluster sample and thus, greater statistical constraining power.

We now have multiple choices for determining the orientation of the cluster. In A22, we fit a 2D Gaussian to the SZ image and determined the cluster orientation accordingly. However, it is problematic to measure the orientation using the same data used to measure the profiles, as this can lead to a noise bias. For example, A22 limited their fits to $R < 0.5R_{200m}$ as at larger radii the noise in the maps biased the shape measurements and the ensuing oriented profile measurements. In this work, we further alleviate this issue by only using pixels within $R < 0.3R_{200m}$ as we do not use or show the profiles in this radial range. This can only partially, not totally, alleviate the noise bias, as the noise in the SZ map can be correlated on large scales due to the presence of the CMB and CIB contaminants.

To make measurements that do not have such biases, we also leverage the optical survey data to obtain two completely independent estimates of the cluster orientation. In particular, we orient the clusters using the shape measurements of the brightest cluster galaxy (BCG) from the DES Y3 shape catalogue, and also using the large-scale density field estimated from the distribution of DES Y3 galaxy positions. The BCG of each cluster is identified with REDMAPPER, and its shape is measured using the METACALIBRATION estimator. To estimate the orientation of the large-scale density field, we compute the Hessian of the smoothed, projected galaxy overdensity field. This is obtained using the methods of Lokken et al. (2022). As a brief description, this Hessian is a matrix of second derivatives with respect to the 2D projected coordinates, $H_{ij} = \frac{\partial^2 \delta}{\partial x_i \partial x_j}$, where $\delta$ is the overdensity field of the galaxy number density and $x_i$ are the projected coordinates. The Hessian is then diagonalized to find the orientation of the major axis. In this work, $\delta$ is given by the galaxy positions of the DES Y3 MAGLIM sample (Porredon et al. 2021) and is smoothed with a Gaussian filter with full width at half-maximum (FWHM) of 20 Mpc. In practice, this produces orientations similar to top-hat smoothing of radius 15 Mpc. On such scales, the shape measurement is dominated by the surrounding large-scale structure (i.e. filaments) and is not impacted by the cluster's own shape. More details on the method can be found in section 3 of Lokken et al. (2022), and the choices used for this analysis are identical to those of that work.

Given two of the three orientation estimates come from the optical data, we focus the analysis of this section on DES clusters. For simplicity, we only show measurements made on the ACT map but note that those of the SPT map are qualitatively similar. Fig. 7 shows the average SZ profiles of DES clusters measured in the three sections, where the orientation is obtained from each of the three methods listed above: the density field's Hessian, the BCG, or the SZ image. We will discuss the results of each orientation method separately.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** A diagram of how we split the cluster into three regions based on the angle away from the major axis. The lightest region falls along the major axis, the darkest along the minor axis, and we also add an intermediate region that is at 45 deg to both axes. Having three regions allows us to more clearly and robustly identify trends as we move from major to minor axis.</figcaption></figure>

**First, the LSS orientation.** The one-halo term of the profile is consistent across all three sections. This is expected as this method measures orientations of the density field on scales of ≈ 15 Mpc, which is well into the two-halo regime of the cluster ($R/R_{200m} \gtrsim 5$). In the far outskirts, $R/R_{200m} > 4$, the measured two-halo term shows a clear striation, where the amplitude grows in the direction of more structure (i.e. the major axis). In the transition regime, $R/R_{200m} \approx 1$, the pressure profile has a steeper derivative along the minor axis. The profiles also show a plateauing feature at $R/R_{200m} \approx 3-4$, where this plateau is found at larger radii along the major axis than the minor axis. This plateau could indicate a shock as has been shown in previous simulation work (Baxter et al. 2021). If this is indeed a shock feature, its dependence on orientation would be consistent with previous work showing the shock boundary is elliptical with the major axis aligned towards the surrounding LSS from which matter is accreted (Aung, Nagai & Lau 2021, see their fig. 1). The prevalence of this feature in all three data subsets suggests it is physical, and also adds some validity to the second minimum seen in the angle-averaged ACT x ACT and SPT x SPT results in Fig. 2. We have verified that all shock behaviours discussed above are also found when the LSS orientations are computed with a different DES Y3 galaxy sample, REDMAGIC (Rodríguez-Monroy et al. 2022).<fnref n="17" /> The consistency of the anisotropic profiles on small-scales is not an artifact of angular resolution limits. For the average cluster, the scale $R = 0.5R_{200m}$ is a factor of ≈2.5 larger than the FWHM of the smoothing for the maps. Furthermore, as we will show below, in Fig. 7, orienting with the SZ image does show clear striations on these small-scales.

**Second, the BCG orientation.** There is now a small striation in the one-halo term, where the profile along the major axis (solid line) has a slightly higher amplitude than that along the minor axis (dotted line). The BCG in massive clusters has a size of roughly ~ 100 kpc, which is a much smaller physical scale than those the

<footnote n="17">¹⁷The galaxies in REDMAGIC are more conservatively selected than MAGLIM. This enables better photometric redshift precision but at the cost of a smaller sample; REDMAGIC has approximately one-third of the galaxy counts of MAGLIM.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Oriented stacking of DES clusters on the ACT SZ map using three different methods to obtain the orientation of each cluster: the large-scale density field estimated using DES MAGLIM galaxies, the brightest cluster galaxy (BCG) shapes from the DES Y3 shape catalogue, and from a 2D Gaussian fit to the SZ image cut-out of each cluster. The LSS (BCG) orientation primarily impacts the two-halo (one-halo) regime of the profile. The SZ orientation (which is measured within $0.3R_{200m}$) causes a much larger difference in the one halo term given we measure the shapes and the profiles on the same map. The BCG shape probes the orientation on scales of the order 100 kpc, while the SZ shape probes 0.5 Mpc to 1 Mpc scales, and the LSS-based method probes ≈ 15 Mpc scales.</figcaption></figure>

other orientation estimates are sensitive to, and the direction towards large-scale structure can change noticeably across different scales (see their fig. 16 in Lokken et al. 2022). We find that orienting by the BCG shape impacts only the one-halo regime. The observed striation in Fig. 7 is the expected consequence of an elliptical cluster profile. This can be seen by taking a circular pressure profile and stretching/squeezing it to make it elliptical. At a fixed physical radius, the profile value along the minor axis will be lower (since the profile has been squeezed radially) compared to the profile value along the major axis (where the profile has been stretched radially). We do not see any clear trends in the two-halo term nor in the transition regime with the log-derivative minimum.

**Finally, the SZ orientation.** This is the technique used previously in A22. There is a significant striation in the one-halo term that is suppressed as we move to the two-halo term. This behaviour is expected as the orientation is measured on the same image used to measure the profiles. Thus the striation in the one-halo term is stronger than when using the BCG or LSS-based orientations. Note, however, that the orientation was measured using data at smaller radial scales than the lower radial limit of the profiles shown here. Near the one-to-two halo transition regime of these profiles, the log-derivative minimum along the cluster major axis (solid line) is steeper than that along the minor axis (dotted line). While this is more statistically significant than the striations seen in the LSS and BCG orientation cases, it is still not significant enough to consider a definite detection. Also, note that though the amplitude of the one-halo term varies between minor axis to major axis, the actual shape of the profile – as seen in the log-derivatives – is consistent in all three directions, up to a radius of $R/R_{200m} \le 0.8$.

In summary, we observe potential behaviours of the log-derivative minima as we shift from major axis to minor axis: when orienting by the large-scale density field, the minimum along the minor axis is steeper than that along the major axis. We also see a potential sign of a shock at much larger radii; namely, the plateauing phase of the

profiles. If this corresponds to an accretion shock, it implies that such oriented stacking could be a more optimal way to detect such features. This is consistent with Aung, Nagai & Lau (2021), who showed the accretion shock is elliptical and pointed along the large-scale structure, and also consistent with Baxter et al. (2021), who found the shock signal is more prominent in the azimuthally averaged profiles of relaxed clusters, as such clusters are predominantly spherical.

## 5.2 Connections to splashback radius

The process of matter accretion can also cause/impact distinct features in other halo profiles, and not just the pressure profile we study here. The splashback radius, which is one such feature, is a physically motivated halo boundary defined by the apocenter in the dark matter phase space of the halo (e.g. Diemer & Kravtsov 2014; Adhikari, Dalal & Chamberlain 2014; More, Diemer & Kravtsov 2015; Mansfield, Kravtsov & Diemer 2017; Aung, Nagai & Lau 2021; Xhakaj et al. 2020; O’Neil et al. 2021; Dacunha et al. 2022). The existence of the splashback feature has been observed by various analyses (More et al. 2016; Baxter et al. 2017; Chang et al. 2018; Shin et al. 2019; Zürcher & More 2019; Murata et al. 2020; Adhikari et al. 2021; Shin et al. 2021), where it is identified as a minimum in the log-derivative of the lensing or galaxy number density profile – similar to how the pressure deficit is a minimum in the log-derivative of the pressure profile – and has been shown to play a role in galaxy formation physics (Baxter et al. 2017; Shin et al. 2019; Adhikari et al. 2021; Dacunha et al. 2022). The ratio of the shock radius and splashback radius, alongside appropriate theoretical models (e.g. Shi 2016), can provide observational constraints on both the adiabatic index of the gas and the mass accretion rate of the cluster (e.g. Hurier, Adam & Keshet 2019).

A22 performed the first comparison of the splashback and pressure deficit features using the SPT-SZ data set. In this work, we supplement this result with a complementary analysis using the ACT
data. A subset of the ACT cluster catalogue has already been used to identify the splashback radius, using galaxy number density profiles and weak lensing profiles measured with DES Y3 galaxies (Shin et al. 2021). Here, we perform the same ACT cluster catalogue selections as Shin et al. (2021), taking all ACT DR5 clusters within the DES Y3 footprint and with $0.1 < z < 0.7$. We then measure the log-derivative of the average SZ profile for this subsample, and present the result in Fig. 8. We also overplot the constraint from Shin et al. (2021) for the splashback radius. Note that this radius was obtained by taking their fits to the observed 2D galaxy number density profile, and using the pipeline from this work to compute the minima of the log-derivative. Thus, we consistently compare the shock and splashback features in 2D, projected profiles. We do not show the weak-lensing result from Shin et al. (2021) but it was shown in that work to be consistent with the galaxy splashback radius.

In Fig. 8, the minimum in the SZ log-derivative, corresponding to the pressure deficit and which we denote as the shock radius $R_{sh}$, coincides with the splashback radius, $R_{sp}$. In particular, we find the ratio to be

$$R_{sp}/R_{sh} = 1.17 \pm 0.20. \qquad (33)$$

The estimate for $R_{sh}$ used above is presented in Table 1 under the name 'ACT x ACT ($R_{sp}$ comparison)'. Our results above show the splashback radius and shock radius are consistent within $0.9\sigma$. These results match those of A22, who found the shock radius in SPT-SZ clusters was also consistent with the splashback radius of that same cluster sample as measured in Shin et al. (2019). While some theoretical works quote a ratio of a shock radius to splashback radius (Molnar et al. 2009; Shi 2016; Aung, Nagai & Lau 2021; Baxter et al. 2021; Zhang et al. 2021), this is for the merger-accelerated accretion shock which is expected to be at $R > R_{200m}$ (Zhang et al. 2019b), and is distinct from the pressure deficit we discuss here. If the deficit is indeed generated by a shock, it would be linked to a typical accretion shock formed via the accretion of gas onto the halo; this is the shock we discussed in the introduction of this work and forms around the virial radius, similar to the splashback feature. Other merger-related shocks within the halo (such as bow shocks from infalling galaxies and gas clumps) can collide with this accretion shock and form a merger-accelerated accretion shock (Zhang et al. 2019b) at larger radii than the original accretion shock.

## 6 SUMMARY AND DISCUSSION

The outskirts of galaxy clusters are where the collapsed halo component interacts most dynamically with the surrounding large-scale structure. A striking feature of this dynamic environment is shocks. The formation and evolution of these shocks have a rich and interesting phenomenology; they form due to the interplay between gravitational infall and hydrodynamical forces, and impact a wide array of cluster astrophysical processes once formed. In this work, we advance on previous studies and use nearly $10^5$ clusters across three data sets – the Dark Energy Survey Year 3, the SPT SZ survey, and the ACT DR4, DR5, and DR6 – to search for shock-generated features in the average pressure profiles of different cluster samples, as measured in different SZ maps. Our key findings are summarized below:

(i) Consistent with A22, there is a pressure deficit at $R/R_{200m} \approx 1.1$ detected at $2.7\sigma$ and $6.1\sigma$ in SPT and ACT, respectively (Fig. 2). This feature is consistent with a shock-driven thermal non-equilibrium between electrons and ions. We do not quote a detection significance

<figure><img src="image_9.png" /><figcaption>**Figure 8.** The log-derivative of the average SZ profile around ACT clusters with $0.15 < z < 0.7$ and whose centres are in the DES footprint. We perform this selection so as to use the same ACT DR5 subsample as Shin et al. (2021), who measured the splashback radius around this sample from the galaxy number density profile (using DES Y3 galaxies) of these clusters. We show their 68 per cent bounds for the projected splashback radius as the vertical blue band. The minimum corresponding to the pressure deficit coincides with the splashback radius. The ratio $R_{sp}/R_{sh} = 1.17 \pm 0.20$, meaning the two projected radii are within $0.9\sigma$ of each other. The dashed line is the prediction of the SZ profile log-derivative for this cluster sample.</figcaption></figure>

for DES clusters given uncertainties in the theoretical modelling (see Section 4.1).

(ii) The SZ maps from SPT and ACT are consistent in both high and low SNR regimes (Fig. 3). For a subset of clusters that lie within both SPT and ACT footprints, we measure the mean SZ profiles using either the SPT map or ACT map and show the profiles are consistent across the entire radial range of our analysis, $0.3 < R/R_{200m} < 10$.

(iii) We construct ACT and DES subsamples with similar mass and redshift distributions and find their mean SZ profiles to be consistent (Fig. 4). This implies that for clusters of a higher mass, $(M_{200m}) = 10^{14.85} M_{\odot}$, the SZ and optical selection effects do not amplify/suppress the deficit feature, and this adds to the robustness of the pressure deficit found in the ACT x ACT and SPT x SPT measurements.

(iv) For optically selected clusters of lower masses, $M_{200m} < 10^{14.5} M_{\odot}$, the radial location of the log-derivative minima are consistent at $R/R_{200m} = 1.1$ while the depth becomes shallower with decreasing mass (Fig. 5).

(v) The SZ profiles measured around group-scale haloes also differ significantly from the theoretical model in the one-halo regime, and are consistent with the model for the two-halo regime (Fig. 5). We discuss three potential causes for this: (i) mass estimation biases, (ii) deviations from the model of Battaglia et al. (2012), (iii) a non-zero correlation in the richness and SZ scatter at fixed mass. All three are more prominent for low-mass clusters, and we find none provide a clear explanation of the observations.

(vi) We perform an oriented stacking of the clusters – with the orientation determined by (i) the large-scale density field comprised of the surrounding structure, (ii) the brightest cluster galaxy, and (iii) a 2D Gaussian fit to the SZ image – and split the profiles into three regions closest-to-furthest from the major axis. When using the LSS orientations, the two-halo term amplitude increases towards the major axis, while the log-derivative depth is steeper along the minor axis (Fig. 7).

(vii) The location of the pressure deficit, $R_{sh}$, is consistent with the splashback radius measured with galaxy number density profiles
in Shin et al. (2021). The ratio is $R_{sp}/R_{sh} = 1.17 \pm 0.20$, and this consistency between shock and splashback radii further signifies the variety of dynamical processes happening at $R \approx R_{200m}$.

Our work, through the use of multiple independent data sets, shows the robustness of the pressure deficit feature in the outskirts of galaxy clusters. While we have discussed this feature as arising from the temperature difference between ions and electrons induced by shock heating, other physical processes could potentially cause this difference. The best way to identify the source of the feature is to obtain the electron number density and electron temperature profiles around these clusters. However, this is quite challenging given the deficit is in the outskirts of the cluster, and it is not possible for X-ray observations – which are the primary way to measure these profiles – to probe these regions for a large enough sample of clusters. Instead, it may be more possible to use high-fidelity X-ray observations of nearby individual clusters to look for such shocks in a small sample of low-redshift clusters.

Compared to A22, we have focused less on the accretion shock feature in this work. While the ACT data shows some potential signs of a feature consistent with SPT data, the amplitude of the signal – and thus the significance of the feature – is low. This is not particularly surprising in that accretion shocks are highly irregular, in both their radial location around the clusters as well as in their geometry (Zhang et al. 2020, 2021). In fact, the simulation-based work of Baxter et al. (2021) found the signal was clearly seen only when selecting relaxed clusters alone. To further pursue a detection of this feature, we can either redo our analyses with the release of a larger cluster catalogue and/or lower noise SZ maps, or perform selection cuts on the current catalogues (particularly related to cluster relaxation) that can maximize the SNR of this feature.

Moving forward there are still additional ongoing and future surveys/data sets that could be used for this work – such as SPT-3G (Benson et al. 2014), Simons Observatory (Ade et al. 2019), and CMB-S4 (Abazajian et al. 2019) – that will all either have higher sensitivity and/or a larger sample of clusters across a broader range in mass and redshift. This would allow the study of the pressure deficit across redshift and mass. Finding clear trends may shed some light on the physical origin of the features. Combining existing data sets can also provide maps with a depth comparable to the upcoming CMB-S4 experiment, and the corresponding cluster catalogues – such as the SPT Megadeep catalogue Kornoelje et al. (in preparation) – will be particularly relevant for comparing the profile outskirts of low-mass SZ-selected and optically selected clusters.

From the optical survey side, we have the cluster samples observed in the Kilo-degree survey (Maturi et al. 2019) and the Dark energy spectroscopic instrument legacy imaging survey (Zou et al. 2021) using richness selection techniques like in DES (but with different algorithms), and samples observed in Hyper Suprime-Cam using a weak-lensing mass selection (Miyazaki et al. 2018; Chen et al. 2020). The Hyper Suprime-Cam sample in particular accesses much higher redshifts than the DES data set. Recently, the sample of X-ray selected clusters has also grown considerably, in part due to the eROSITA All-sky X-ray mission (Liu et al. 2022). While X-ray samples have significantly lower redshift than the SZ and optical samples, they allow the pursuit of unique science cases – X-ray clusters are bimodal in whether or not they have a cool core, and measuring the SZ profile outskirts around the two different types of clusters could shed light on the interplay between the physics of the outskirts and that of the cluster core. Opportunities also exist for studying the correlations between profiles, and these can have strong astrophysical signatures (e.g. Farahi, Nagai & Anbajagane

2022b). Techniques have also been developed to extract such profile correlations in a data-driven manner, with minimal assumptions, such as Gaussian processes (Farahi, Nagai & Chen 2021) and local linear regression (Farahi, Anbajagane & Evrard 2022a).

Thus, there are many synergistic opportunities for cross-correlating the different types of data sets – both ongoing and upcoming – and each combination will allow us to access different science cases regarding the physics of these cluster outskirts. The use of three independent, wide-field surveys in this work – all analyzed under a common, coherent framework – has given us the ability to easily cross-check and validate the signatures we see, and in general, be less sensitive to both known, and unknown, systematic effects. The ability to perform such tests and explorations will only grow, as we move into the age of even larger surveys with higher overlap and greater synergies.

## ACKNOWLEDGEMENTS

DA is supported by National Science Foundation (NSF) grant no. 2108168. CC is supported by the Henry Luce Foundation and DOE grant DE-SC0021949. MHi acknowledges support from the National Research Foundation of South Africa (grant no. 137975). KM acknowledges support from the National Research Foundation of South Africa. CS acknowledges support from the Agencia Nacional de Investigación y Desarrollo (ANID) through FONDECYT grant no. 11191125 and BASAL project FB210003. JCH acknowledges support from NSF grant AST-2108536, NASA grants 21-ATP21-0129 and 22-ADAP22-0145, the Sloan Foundation, and the Simons Foundation.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. NSF, 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's 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 Cerro Tololo Inter-American Observatory at NSF's 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 NSF.

The DES data management system is supported by the NSF under Grant Numbers 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 programme 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.

The SPT programme is supported by the NSF through the grant no. OPP-1852617. Partial support is also provided by the Kavli Institute of Cosmological Physics at the University of Chicago.

Support for ACT was through the U.S. NSF through awards AST-0408698, AST-0965625, and AST-1440226 for the ACT project, as well as awards PHY-0355328, PHY-0855887 and PHY-1214379. Funding was also provided by Princeton University, the University of Pennsylvania, and a Canada Foundation for Innovation (CFI) award to UBC. ACT operated in the Parque Astronómico Atacama in northern Chile under the auspices of the Agencia Nacional de Investigación y Desarrollo (ANID). The development of multichroic detectors and lenses was supported by NASA grants NNX13AE56G and NNX14AB58G. Detector research at NIST was supported by the NIST Innovations in Measurement Science programme. Computing for ACT was performed using the Princeton Research Computing resources at Princeton University, the National Energy Research Scientific Computing Center (NERSC), and the Niagara supercomputer at the SciNet HPC Consortium. SciNet is funded by the CFI under the auspices of Compute Canada, the Government of Ontario, the Ontario Research Fund-Research Excellence, and the University of Toronto. We thank the Republic of Chile for hosting ACT in the northern Atacama, and the local indigenous Licanantay communities whom we follow in observing and learning from the night sky.

All analysis in this work was enabled greatly by the following software: PANDAS (McKinney 2011), NUMPY (Van der Walt, Colbert & Varoquaux 2011), SciPy (Virtanen et al. 2020), and MATPLOTLIB (Hunter 2007). We have also used the Astrophysics Data Service (ADS) and arXiv preprint repository extensively during this project and the writing of the paper.

## DATA AVAILABILITY

All SPT data used in our analyses are publicly available at the repositories linked to in this paper. The cluster catalogue from ACT DR5 is publicly available, alongside the maps from which these catalogues are constructed. The SZ map of ACT DR6, and the raw DR6 maps used to construct it, will be made public shortly. The DR4 maps used to construct this SZ map are already available. The DES Year 3 shape catalogue is also available, while the cluster catalogue is not yet public. The links to the online portals hosting the publicly available catalogues can be found under the relevant data subsection in Section 2.

The code used to generate the theoretical tSZ profile of a halo, including both one-halo and two-halo contributions, is made available at https://github.com/DhayaaAnbajagane/tSZ.Profiles.

## REFERENCES

Abazajian K. et al., 2019, preprint (arXiv:1907.04473)

Ackermann M. et al., 2014, ApJ, 787, 18

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

Adhikari S., Dalal N., Chamberlain R. T., 2014, J. Cosmol. Astropart. Phys., 2014, 019

Adhikari S. et al., 2021, ApJ, 923, 37

Akahori T., Yoshikawa K., 2010, PASJ, 62, 335

Akahori T., Yoshikawa K., 2012, PASJ, 64, 12

Akamatsu H., Hoshino A., Ishisaki Y., Ohashi T., Sato K., Takei Y., Ota N., 2011, PASJ, 63, S1019

Akamatsu H. et al., 2016, A&A, 593, L7

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

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

Anbajagane D., Evrard A. E., Farahi A., Barnes D. J., Dolag K., McCarthy I. G., Nelson D., Pillepich A., 2020, MNRAS, 495, 686

Anbajagane D., Evrard A. E., Farahi A., 2022a, MNRAS, 509, 3441

Anbajagane D. et al., 2022b, MNRAS, 510, 2980

Anbajagane D. et al., 2022c, MNRAS, 514, 1645

Anbajagane D., Chang C., Lee H., Gatti M., 2023a, preprint (arXiv:2310.02349)

Anbajagane D. et al., 2023b, MNRAS, 526, 5530

Aung H., Nagai D., Lau E. T., 2021, MNRAS, 508, 2071

Avestruz C., Nagai D., Lau E. T., Nelson K., 2015, ApJ, 808, 176

Basu K., Sommer M., Erler J., Eckert D., Vazza F., Magnelli B., Bertoldi F., Tozzi P., 2016, ApJ, 829, L23

Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., 2012, ApJ, 758, 75

Baxter E. et al., 2017, ApJ, 841, 18

Baxter E. J. et al., 2018, MNRAS, 476, 2674

Baxter E. J., Adhikari S., Vega-Ferrero J., Cui W., Chang C., Jain B., Knebe A., 2021, MNRAS, 508, 1777

Baxter E. J., Pandey S., Adhikari S., Cui W., Shin T.-h., Li Q., Rasia E., 2023, preprint (arXiv:2304.08731)

Beltz-Mohrmann G. D., Berlind A. A., 2021, ApJ, 921, 112

Benson B. A. et al., 2014, preprint (arXiv:1407.2973)

Bertschinger E., 1985, ApJS, 58, 39

Blandford R., Eichler D., 1987, Phys. Rep., 154, 1

Bleem L. E. et al., 2015, ApJS, 216, 27

Bleem L. E. et al., 2020, ApJS, 247, 25

Bleem L. E. et al., 2022, ApJS, 258, 36

Bocquet S. et al., 2019, ApJ, 878, 55

Boselli A., Fossati M., Sun M., 2022, A&A Rev., 30, 3

Brunetti G., Jones T. W., 2014, Int. J. Mod. Phys. D, 23, 1430007

Carlstrom J. E., Holder G. P., Reese E. D., 2002, ARA&A, 40, 643

Carlstrom J. E. et al., 2011, PASP, 123, 568

Chang C. et al., 2018, ApJ, 864, 83

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

Chen K.-F., Oguri M., Lin Y.-T., Miyazaki S., 2020, ApJ, 891, 139

Child H. L., Habib S., Heitmann K., Frontiere N., Finkel H., Pope A., Morozov V., 2018, ApJ, 859, 55

Chisari N. E. et al., 2019, ApJS, 242, 2

Costanzi M. et al., 2021, Phys. Rev. D, 103, 043522

Coulton W. R. et al., 2023, preprint (arXiv:2307.01258)

Cui W. et al., 2018, MNRAS, 480, 2898

Cui W. et al., 2022, MNRAS, 514, 977

Dacunha T., Belyakov M., Adhikari S., Shin T.-h., Goldstein S., Jain B., 2022, MNRAS, 512, 4378

Delabrouille J., Cardoso J. F., 2009, Diffuse Source Separation in CMB Observations. p. 159

Di Mascolo L., Churazov E., Mroczkowski T., 2019a, MNRAS, 487, 4037

Di Mascolo L. et al., 2019b, A&A, 628, A100

Diemer B., 2018, ApJS, 239, 35
Diemer B., Joyce M., 2019, ApJ, 871, 168

Diemer B., Kravtsov A. V., 2014, ApJ, 789, 1

Drury L. O., 1983, Rep. Prog. Phys., 46, 973

Erler J., Basu K., Chluba J., Bertoldi F., 2018, MNRAS, 476, 3360

Ettori S., Fabian A. C., 1998, MNRAS, 293, L33

Evrard A. E., Arnault P., Huterer D., Farahi A., 2014, MNRAS, 441, 3562

Farahi A., Evrard A. E., Rozo E., Rykoff E. S., Wechsler R. H., 2016, MNRAS, 460, 3900

Farahi A., Evrard A. E., McCarthy I., Barnes D. J., Kay S. T., 2018, MNRAS, 478, 2618

Farahi A. et al., 2019, Nat. Commun., 10, 2504

Farahi A., Nagai D., Chen Y., 2021, AJ, 161, 30

Farahi A., Anbajagane D., Evrard A. E., 2022a, ApJ, 931, 166

Farahi A., Nagai D., Anbajagane D., 2022b, ApJ, 933, 48

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASF, 125, 306

Fox D. C., Loeb A., 1997, ApJ, 491, 459

Gatti M. et al., 2020, MNRAS, 498, 4060

Gatti M. et al., 2021, MNRAS, 504, 4312

Gatti M. et al., 2022, Phys. Rev. D, 105, 123525

Ghirardini V., Ettori S., Amodeo S., Capasso R., Sereno M., 2017, A&A, 604, A100

Ghirardini V., Ettori S., Eckert D., Molendi S., Gastaldello F., Pointecouteau E., Hurier G., Bourdin H., 2018, A&A, 614, A7

Hajian A., Battaglia N., Spergel D. N., Bond J. R., Pfrommer C., Sievers J. L., 2013, J. Cosmol. Astropart. Phys., 2013, 064

Ha J.-H., Ryu D., Kang H., 2018, ApJ, 857, 26

Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399

Hill J. C., Pajer E., 2013, Phys. Rev. D, 88, 063526

Hill J. C., Baxter E. J., Lidz A., Greco J. P., Jain B., 2018, Phys. Rev. D, 97, 083501

Hilton M. et al., 2021, ApJS, 253, 3

Hojjati A. et al., 2017, MNRAS, 471, 1565

Hong S. E., Ryu D., Kang H., Cen R., 2014, ApJ, 785, 133

Hong S. E., Kang H., Ryu D., 2015, ApJ, 812, 49

Horowitz B., Seljak U., 2017, MNRAS, 469, 394

Hou K.-C., Hallinan G., Keshet U., 2023, MNRAS, 521, 5786

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Hurier G., Adam R., Keshet U., 2019, A&A, 622, A136

Ishiyama T. et al., 2021, MNRAS, 506, 4210

Kang H., Ryu D., Cen R., Ostriker J. P., 2007, ApJ, 669, 729

Kocsis B., Haiman Z., Frei Z., 2005, ApJ, 623, 632

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

Lau E. T., Nagai D., Avestruz C., Nelson K., Vikhlinin A., 2015, ApJ, 806, 68

Lau E. T., Hearin A. P., Nagai D., Cappelluti N., 2021, MNRAS, 500, 1029

Lee E. et al., 2022, MNRAS, 517, 5303

Liu A. et al., 2022, A&A, 661, A2

Lokken M. et al., 2022, ApJ, 933, 134

Lyskova N., Churazov E., Khabibullin I. I., Burenin R., Starobinsky A. A., Sunyaev R., 2023, preprint (arXiv:2305.07080)

Ma Y.-Z., Van Waerbeke L., Hinshaw G., Hojjati A., Scott D., Zuntz J., 2015, J. Cosmol. Astropart. Phys., 2015, 046

Mansfield P., Kravtsov A. V., Diemer B., 2017, ApJ, 841, 34

Maturi M., Bellagamba F., Radovich M., Roncarelli M., Sereno M., Moscardini L., Bardelli S., Puddu E., 2019, MNRAS, 485, 498

McClintock T. et al., 2019, MNRAS, 482, 1352

McDonald M. et al., 2014, ApJ, 794, 67

McKinney W., 2011, Python for High Performance and Scientific Computing. p. 14

Melin J. B., Pratt G. W., 2023, preprint (arXiv:2304.09041)

Miniati F., Ryu D., Kang H., Jones T. W., Cen R., Ostriker J. P., 2000, ApJ, 542, 608

Miniati F., Jones T. W., Ryu D., Kang H., 2001, in Simon M., Lorenz E., Pohl M., eds, 27th Internat. Cosmic Ray Conf., Nonthermal Radiation from Groups and Clusters of Galaxies: Numerical Simulation Perspective. Hamburg, p. 2062

Miyazaki S. et al., 2018, PASJ, 70, S27

Molnar S. M., Hearn N., Haiman Z., Bryan G., Evrard A. E., Lake G., 2009, ApJ, 696, 1640

More S., Diemer B., Kravtsov A. V., 2015, ApJ, 810, 36

More S. et al., 2016, ApJ, 825, 39

Mroczkowski T. et al., 2019, Space Sci. Rev., 215, 17

Murata R., Sunayama T., Oguri M., More S., Nishizawa A. J., Nishimichi T., Osato K., 2020, PASJ, 72, 64

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

Nelson D., Genel S., Pillepich A., Vogelsberger M., Springel V., Hernquist L., 2016, MNRAS, 460, 2881

O'Neil S., Barnes D. J., Vogelsberger M., Diemer B., 2021, MNRAS, 504, 4649

Osato K., Flender S., Nagai D., Shirasaki M., Yoshida N., 2018, MNRAS, 475, 532

Osato K., Shirasaki M., Miyatake H., Nagai D., Yoshida N., Oguri M., Takahashi R., 2020, MNRAS, 492, 4780

Pandey S. et al., 2019, Phys. Rev. D, 100, 063519

Pandey S., Baxter E. J., Hill J. C., 2020, Phys. Rev. D, 101, 043525

Pandey S. et al., 2022, Phys. Rev. D, 105, 123526

Pfrommer C., Springel V., Enßlin T. A., Jubelgas M., 2006, MNRAS, 367, 113

Pfrommer C., Enßlin T. A., Springel V., Jubelgas M., Dolag K., 2007, MNRAS, 378, 385

Planck Collaboration XII, 2011, A&A, 536, A12

Planck Collaboration V, 2013, A&A, 550, A131

Planck Collaboration XXVII, 2016, A&A, 594, A27

Planck Collaboration LVII, 2020, A&A, 643, A42

Planelles S., Quilis V., 2013, MNRAS, 428, 1643

Planelles S. et al., 2021, MNRAS, 507, 5703

Porredon A. et al., 2021, Phys. Rev. D, 103, 043503

Pratt C. T., Qu Z., Bregman J. N., 2021, ApJ, 920, 104

Quilis V. M. J., Ibáñez, Sáez D., 1998, ApJ, 502, 518

Robertson N. C. et al., 2023, preprint (arXiv:2304.10219)

Rodríguez-Monroy M. et al., 2022, MNRAS, 511, 2665

Romero C. E. et al., 2017, ApJ, 838, 86

Romero C. et al., 2018, A&A, 612, A39

Rudd D. H., Nagai D., 2009, ApJ, 701, L16

Rykoff E. S. et al., 2014, ApJ, 785, 104

Ryu D., Kang H., Hallman E., Jones T. W., 2003, ApJ, 593, 599

Sánchez J. et al., 2023, MNRAS, 522, 3163

Saro A. et al., 2017, MNRAS, 468, 3347

Sayers J. et al., 2013, ApJ, 768, 177

Sayers J. et al., 2016, ApJ, 832, 26

Sayers J., Mantz A. B., Rasia E., Allen S. W., Cui W., Golwala S. R., Morris R. G., Wan J. T., 2023, ApJ, 944, 221

Schaal K., Springel V., 2015, MNRAS, 446, 3992

Schaal K. et al., 2016, MNRAS, 461, 4441

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

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

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

Secco L. F. et al., 2022b, Phys. Rev. D, 105, 103537

Sehgal N. et al., 2013, ApJ, 767, 38

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

Shao M., Anbajagane D., 2023, preprint (arXiv:2311.03491)

Shao M. J., Anbajagane D., Chang C., 2023, MNRAS, 523, 3258

Sheldon E. S., Huff E. M., 2017, ApJ, 841, 24

Shi X., 2016, MNRAS, 461, 1804

Shin T. et al., 2019, MNRAS, 487, 2900

Shin T. et al., 2021, MNRAS, 507, 5758

Shin T.-h, Diemer B., 2023, MNRAS, 521, 5570

Shirasaki M., Lau E. T., Nagai D., 2020, MNRAS, 491, 235

Skillman S. W., O'Shea B. W., Hallman E. J., Burns J. O., Norman M. L., 2008, ApJ, 689, 1063

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

Tanimura H., Douspis M., Aghanim N., Salvati L., 2022, MNRAS, 509, 300
Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709

Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878

To C. et al., 2021, Phys. Rev. Lett., 126, 141301

Van der Walt S., Colbert S. C., Varoquaux G., 2011, Comput. Sci. Eng., 13, 22

Vazza F., Brunetti G., Gheller C., 2009, MNRAS, 395, 1333

Vazza F., Brüggen M., Gheller C., 2013, MNRAS, 428, 2366

Vazza F., Gheller C., Brüggen M., 2014, MNRAS, 439, 2662

Vikram V., Lidz A., Jain B., 2017, MNRAS, 467, 2315

Vink J., Broersen S., Bykov A., Gabici S., 2015, A&A, 579, A13

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52

Wong K.-W., Sarazin C. L., 2009, ApJ, 707, 1141

Xhakaj E., Diemer B., Leauthaud A., Wasserman A., Huang S., Luo Y., Adhikari S., Singh S., 2020, MNRAS, 499, 3534

Zhang Y. et al., 2019a, MNRAS, 487, 2578

Zhang C., Churazov E., Forman W. R., Lyskova N., 2019b, MNRAS, 488, 5259

Zhang C., Churazov E., Dolag K., Forman W. R., Zhuravleva I., 2020, MNRAS, 494, 4539

Zhang C., Zhuravleva I., Kravtsov A., Churazov E., 2021, MNRAS, 506, 839

Zhu Z. et al., 2021, A&A, 652, A147

Zinger E., Dekel A., Kravtsov A. V., Nagai D., 2016, preprint (arXiv:1610.02644)

Zou H. et al., 2021, ApJS, 253, 56

Zürcher D., More S., 2019, ApJ, 874, 184

## APPENDIX A: IMPACT OF MISCENTRING ON PROFILES

As we discussed previously, our theoretical model for the SZ profiles of *optically selected* clusters depends on the miscentring model parameters assumed for the cluster sample. In Fig. A1, we show the model for the DES clusters’ SZ profile as we vary the amplitude of miscentring, $\tau_{miscen}$, and the fraction of miscentred clusters, $f_{miscen}$, defined in Section 3.2.1. We do not have an external, calibrated constraint for the miscentring effect in this specific DES cluster

sample. Thus, as an alternative, we vary the parameter values until the theory visually matches the data for the one-halo term as shown in Fig. 2. In practice, we do this by making predictions in a 5x5 grid of parameter values, and find $\tau_{miscen} = 0.9$ and $f_{miscen} = 0.4$ to be the best combination. As was noted before, both values are near the $3-4\sigma$ upper limit of constraints on the DES Y1 cluster sample, depending on the parameter (Zhang et al. 2019a, see their Chandra-DES constraints in Table 1), while the value of $\tau_{miscen}$ is within $1\sigma$ of the estimate from Bleem et al. (2020, see their table 6), which is based on a SPT-DES matched cluster sample. It is also generally consistent with the work of Sehgal et al. (2013), who find the offsets in individual clusters seen in ACT have upper limits of 1.5 Mpc, which corresponds to $\tau_{miscen} \approx 1.5$.

Given these potential limitations of the implemented miscentring effects in our work, we focus our analysis of optically selected cluster on results that do not require accurate theoretical estimates of pressure profiles for these clusters. We specifically avoid quoting a detection significance of shock features in these clusters given the uncertainty in the miscentring model parameters of the theoretical model. Our results in Section 4.3, which *does* compare theory and data for low mass DES clusters and finds large deviations, are insensitive to miscentring as the deviations are significantly larger than those from miscentring effects alone.

Fig. A1 shows that for the high-mass sample (left panels), the variation in $\tau_{miscen}$ changes the location of the log-derivative minimum from $r_{min} = 0.7R_{200m} \rightarrow 1.2R_{200m}$ as we vary $\tau_{miscen} = 0.1 \rightarrow 1.2$. The minimum value of the log-derivative goes from $-2.5 \rightarrow -3.0$ as we vary $f_{miscen} = 0.1 \rightarrow 0.8$. In particular, the result for $f_{miscen} = 0.8$ and $\tau_{miscen} = 0.9$ appears to replicate a shock-esque feature at $R \approx R_{200m}$. However, this is not an indication that shock features can be explained by miscentring. For SZ-selected clusters, the miscentring is much smaller than the values being considered in Fig. A1. While the predicted profile for large miscentring values forms a deficit-like feature, the agreement in the one-halo term is significantly degraded as a result. Thus, this is not evidence that miscentring is the cause of the deficit feature, and is instead evidence of the miscentring model’s ability to capture steep drops in the pressure profile.
<figure><img src="image_10.png" /></figure>

**Figure A1.** The DESxSPT theory predictions – both mean profile (top) and log-derivative (bottom) – when varying the parameters of the miscentring model. We vary $\tau_{\text{miscen}}$, which is the miscentring length scale (top two panels), and miscentring fraction $f_{\text{miscen}}$ (bottom two panels). The different columns show the model for different mass ranges. Miscentring transfers power from small scales ($R < 0.3R_{200m}$) to large scales ($R \sim R_{200m}$). This is distinct from the impacts of contamination, which will suppress the one-halo term over all scales, but similar to feedback which will push gas from small scales out to large scales.

## APPENDIX B: DEPENDENCE ON CLUSTER SNR

The analysis in Section 4.2 implies the pressure deficit is not formed due to SZ selection effects. A characteristic of a feature driven by noise-effects is its amplitude grows near the limit of the selection threshold. For SZ-selected clusters, this is the SNR threshold, which

is SNR > 4.5 for ACT and SNR > 4 for SPT. In Fig. B1, we show the average SZ profiles of cluster subsamples split by their SNR. The pressure deficit feature exists in all SNR bins, and the values for the location and log-derivative of the deficit (listed in Table 1) are consistent within <1$\sigma$. This adds further to the evidence that the feature is not formed from SZ noise-based selection effects.
<figure><img src="image_11.png" /><figcaption>**Figure B1.** The average SZ profiles of SZ-selected clusters measured on their respective SZ maps, with the samples split by their SNR. All subsamples show pressure deficits, confirming that there is no SNR dependence on the deficit feature. The location and depth of the log-derivative minima are listed in Table 1.</figcaption></figure>

## APPENDIX C: ANALYTIC FIT FOR PRESSURE DEFICIT

In Fig. 2, the halo model is a good match to the measured profiles in both the one-halo and two-halo regimes, but the pressure deficit feature in the transition region is only seen in the measurements. We

<figure><img src="image_12.png" /><figcaption>**Figure C1.** The constraints on the parameters of equation (C1) obtained from the SPT or ACT data. The values are listed in Table C1. The constraints from SPT and ACT are consistent with each other within <0.5σ, which independently validates our statement that the pressure deficit features shown in Fig. 2 for these samples is consistent.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure C2.** The unsmoothed average SZ profile of the ACT sample (top) and SPT sample (bottom). Overlaid is the original theory of the total halo model in black dashed line, and the modified halo model, described in equation (C1), in the solid black line. The modified model provides a better fit to the data, particularly to the pressure deficit term. Note that the original halo model is a prediction and not a fit, whereas the modified halo theory fits the additional Gaussian component while keeping the original halo theory component fixed.</figcaption></figure>
<table><caption>**Table C1.** The best-fitting parameters of the modified halo model (A, μ, σ) for both SPT and ACT data. The $\chi^2_{orig}$ and $\chi^2_{mod}$ columns are the chi-squared for the original and modified halo models. The modified model is significantly better in both cases.</caption><thead><tr><th>Data set</th><th>A</th><th>μ</th><th>σ</th><th>χ<sup>2</sup><sub>mod</sub></th><th>χ<sup>2</sup><sub>orig</sub></th></tr></thead><tbody><tr><td>SPT x SPT</td><td>-0.47<sup>+0.21</sup><sub>-0.21</sub></td><td>1.39<sup>+0.36</sup><sub>-0.15</sub></td><td>0.31<sup>+0.18</sup><sub>-0.10</sub></td><td>23.41</td><td>32.86</td></tr><tr><td>ACT x ACT</td><td>-0.37<sup>+0.10</sup><sub>-0.10</sub></td><td>1.49<sup>+0.13</sup><sub>-0.10</sub></td><td>0.29<sup>+0.06</sup><sub>-0.05</sub></td><td>65.97</td><td>86.43</td></tr></tbody></table>

implement a simple modification to the existing halo model theory to match this effect. Our modification multiplies the original theory by a Gaussian,

$$ \xi_{h,p}^{\text{mod}}(r, A, \mu, \sigma) = \xi_{h,p} \times \left( 1 + A \frac{\mathcal{M}(r, \mu, \sigma)}{\mathcal{M}(\mu, \mu, \sigma)} \right), \quad \text{(C1)} $$

where A is the amplitude of the Gaussian, μ is the mean/location in units of $R/R_{200m}$, σ is the width of the Gaussian, r is the comoving distance bins mentioned previously in Section 3.2, and $\xi_{h,p}$ is the halo-pressure correlation computed in that same section, accounting for beam-smoothing effects. Equation (C1) shows that we additionally normalize the Gaussian feature, $\mathcal{M}(r, \mu, \sigma)$, by another Gaussian evaluated at the mean, $\mathcal{M}(\mu, \mu, \sigma)$, and this ensures that only the parameter A controls the amplitude of the Gaussian.<fnref n="18" />

We then fit the model in equation (C1) to the measured profiles and obtain constraints on the three parameters, A, μ, and σ. The fit is done using the Markov Chain Monte Carlo (MCMC) technique as implemented in the EMCEE package (Foreman-Mackey et al. 2013). We use a set of mostly uninformative priors for all parameters,

$$ \begin{aligned} & -10 < A < 0, \\ & 0.8 < \mu < 2.5, \\ & 0.05 < \sigma < 0.9. \end{aligned} \quad \text{(C2)} $$

The bounds of μ, and σ are chosen to prevent the fitting of either random noise fluctuations in the profiles or differences between data and theory near the cluster core. We enforce $A < 0$ as we are fitting a deficit feature and so the Gaussian must suppress (and not amplify) the pressure in the existing theory prediction. The fit is performed using only $R > 0.6R_{200m}$, and this is done primarily to prevent the MCMC from focusing on any deviations between data and theory in the cluster core. Limiting our analysis to this radial range also helps obtain a numerically stable covariance matrix for use in the MCMC. The fitting is performed by minimizing the $\chi^2$ for the log-likelihood:

$$ \chi^2 = (y^{\text{obs}} - y^{\text{th}}) C^{-1} (y^{\text{obs}} - y^{\text{th}}), \quad \text{(C3)} $$

where C is the covariance matrix of the profile, $y^{\text{obs}}$ is the measured profile, and $y^{\text{th}}$ is the modified halo model. The MCMC is run on

the raw, unsmoothed profiles using 300 walkers and 40 000 steps per walker.

We show the results for the SPT and ACT cluster catalogues, each measured on their corresponding SZ maps. The parameters corresponding to this fit are shown in Fig. C1, while the fits and the data are compared in Fig. C2. The latter shows that the modified halo model is a good fit to the data, and better than the original total halo model, in the one-to-two halo transition regime with the pressure deficit feature. The constraints from the ACT and SPT samples shown in Fig. C1 are consistent with each other. Table C1 lists the amplitude, size, and location of the pressure deficit, and we see the values are consistent within $<0.5\sigma$ across the two samples. We also list the $\chi^2$ of the fit from the original halo model and the modified halo model. For both data sets, the $\chi^2$ is noticeably improved, and we can also see this visually in the fits of Fig. C2. These fits can be used as a simple technique to include SPT/ACT-like pressure deficits in an existing halo model.

In our main analysis, we have not performed any fits, which has primarily been due to the lack of a model for the pressure deficit. In this section, we now show we do have such a model. However, it is not used in our main analysis as we have not yet studied its robustness and validated it against any potential biases. For example, we have already found that it is fairly easy for this model to fit features other than the pressure deficit – like fluctuations on small scales – and the priors must be slightly hand-tuned to make the model focus on the deficit. In our case, the bounds on σ and μ were tuned so as to avoid such issues when fitting the two mean SZ profiles we present here. It is unlikely these priors can be used generically for all measured mean profiles without running into fit failures or prior boundary effects. Thus, while the fits we describe here are a useful phenomenological model, they have not been validated at the same rigor as our current pipeline (which was tested extensively in A22) and so we have continued with the original pipeline for the main analysis in this work.

In the future, one could also use this technique – namely, the posteriors of the model parameters – as an alternative estimator of the detection significance for the pressure deficit, where $A = 0$ would denote no detection of the deficit. While the results of Table C1 already provide the relevant numbers for this work, we do not quote this detection significance as we have yet to validate our profile-fitting technique adequately. Thus, the main purpose and result of this section remain the fits that enable a simple, data-driven replication of the deficit feature in a halo model prediction.

## APPENDIX D: CORRELATION MATRIX

Fig. D1 presents the correlation matrix of the ACT x ACT log-derivative measurements. It is a typical diagonal matrix, with some correlations in adjacent bins due to the smoothing procedure (see Section 3.1). The white box highlights the bins used to estimate the $\chi^2$ shown in Table 1. The top panel shows the log-derivative measurement, now in discrete points corresponding to the binning, corresponding to the ACT x ACT measurement of Fig. 2.

<footnote n="18">¹⁸Without this renormalization, the parameter σ would also control the amplitude in addition to A, given the Gaussian goes as $N \propto 1/\sqrt{2\pi\sigma^2}$. This then causes degeneracies in the parameter space that lead to problems in the Bayesian inference.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure D1.** The top panel shows the log-derivative for the ACT x ACT measurement (presented in Fig. 2), now presented as discrete points corresponding to the measurements within each bin. The dashed line is the theoretical prediction, and the grey band shows the range of scales used to quantify the detection of a pressure deficit. The bottom panel shows the correlation matrix of the log-derivative. The matrices for the other measurements have the same structure. The white box demarcates the bins used in computing the $\chi^2$ of the pressure deficit, as listed in Table 1.</figcaption></figure>

## APPENDIX E: IMPACT OF CIB CONTAMINATION

We explicitly test the impact of the CIB on our SZ profile measurements by comparing the fiducial maps with those where the CIB signal is deprojected/minimized in the final maps. See Bleem et al. (2022) and Coulton et al. (2023) for details on the deprojection procedure of SPT and ACT, respectively. In Fig. E1, we compare measurements of the log-derivative made on these two sets of maps for both SPT and ACT. The measurements are statistically consistent across the maps with and without CIB deprojection. The fiducial SPT map already removes a significant fraction of the CIB, as discussed in Bleem et al. (2022, see their section 3.5). The ACT data contain

<figure><img src="image_15.png" /><figcaption>**Figure E1.** The log-derivatives of the SPTxSPT and ACTxACT data sets, measured using the fiducial maps (colored lines) and maps with the CIB deprojected (black lines). The impact of the CIB on this measurement is negligible as the two versions are statistically consistent across the range of scales.</figcaption></figure>

multiple CIB-deprojected maps, and we use the fiducial one (Coulton et al. 2023, see their section C.3 and equation 18).

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

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

<footnote n="3">³Institute for Astronomy, University of Hawaii'i, 2680 Woodlawn Drive, Honolulu, HI 96822, USA</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="5">⁵David A. Dunlap Department of Astronomy & Astrophysics, University of Toronto, 50 St George Str, Toronto, ON M5S 3H4, Canada</footnote>

<footnote n="6">⁶Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St George Str, Toronto, ON M5S 3H4, Canada</footnote>

<footnote n="7">⁷Dunlap Institute of Astronomy & Astrophysics, 50 St George Str, Toronto, ON M5S 3H4, Canada</footnote>

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

<footnote n="9">⁹Fermi National Accelerator Laboratory, PO Box 500, Batavia, IL 60510, USA</footnote>

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

<footnote n="11">¹¹Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>
<footnote n="12">¹²Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="13">¹³Department of Physics and Astronomy, University of Southern California, Los Angeles, CA 90089, USA</footnote>

<footnote n="14">¹⁴Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK</footnote>

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

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

<footnote n="17">¹⁷Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA</footnote>

<footnote n="18">¹⁸Fermi National Accelerator Laboratory, MS209, PO Box 500, Batavia, IL 60510, USA</footnote>

<footnote n="19">¹⁹High Energy Physics Division, Argonne National Laboratory, Argonne, IL 60439, USA</footnote>

<footnote n="20">²⁰Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81679 Munich, Germany</footnote>

<footnote n="21">²¹Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, ON, M5S 3H8, Canada</footnote>

<footnote n="22">²²Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

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

<footnote n="24">²⁴Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain</footnote>

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

<footnote n="26">²⁶Center for Astrophysical Surveys, National Center for Supercomputing Applications, 1205 West Clark Str, Urbana, IL 61801, USA</footnote>

<footnote n="27">²⁷Department of Physics, Duke University Durham, NC 27708, USA</footnote>

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

<footnote n="29">²⁹Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11, I-34131 Trieste, Italy</footnote>

<footnote n="30">³⁰INAF – Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="31">³¹Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy</footnote>

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

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

<footnote n="34">³⁴Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029 Hamburg, Germany</footnote>

<footnote n="35">³⁵School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia</footnote>

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

<footnote n="37">³⁷Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India</footnote>

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

<footnote n="39">³⁹Department of Physics and Astronomy, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada</footnote>

<footnote n="40">⁴⁰Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029 Blindern, NO-0315 Oslo, Norway</footnote>

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

<footnote n="42">⁴²Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

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

<footnote n="44">⁴⁴Universität Innsbruck, Institute für Astro- und Teilchenphysik, Technik-strasse 25/8, 6020 Innsbruck, Austria</footnote>

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

<footnote n="46">⁴⁶School of Physics and Astronomy, Cardiff University, CF24 3AA, UK</footnote>

<footnote n="47">⁴⁷Department of Physics, Columbia University, 538 West 120th Street, New York, NY 10027, USA</footnote>

<footnote n="48">⁴⁸Wits Centre for Astrophysics, School of Physics, University of the Witwatersrand, Private Bag 3, 2050, Johannesburg, South Africa</footnote>

<footnote n="49">⁴⁹Astrophysics Research Centre, School of Mathematics, Statistics, and Computer Science, University of KwaZulu-Natal, Westville Campus, Durban 4041, South Africa</footnote>

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

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

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

<footnote n="53">⁵³Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="54">⁵⁴Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia</footnote>

<footnote n="55">⁵⁵Lowell Observatory, 1400 Mars Hill Rd, Flagstaff, AZ 86001, USA</footnote>

<footnote n="56">⁵⁶Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, UK</footnote>

<footnote n="57">⁵⁷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="58">⁵⁸Kavli Institute for Particle Astrophysics & Cosmology, PO Box 2450, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="59">⁵⁹Department of Physics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="60">⁶⁰Enrico Fermi Institute, University of Chicago, Chicago, IL 60637, USA</footnote>

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

<footnote n="62">⁶²Astrophysics Research Centre, University of KwaZulu-Natal, Westville Campus, Durban 4041, South Africa</footnote>

<footnote n="63">⁶³European Southern Observatory, Karl-Schwarzschild-Str. 2, Garching D-85748, Germany</footnote>

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

<footnote n="65">⁶⁵Institute of Theoretical Astrophysics, University of Oslo, 0371 Oslo, Norway</footnote>

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

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

<footnote n="68">⁶⁸Joseph Henry Laboratories of Physics, Jadwin Hall, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="69">⁶⁹Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania 15312, USA</footnote>

<footnote n="70">⁷⁰Department of Physics and Astronomy, Haverford College, Haverford, PA 19041, USA</footnote>

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

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

<footnote n="73">⁷³School of Physics, University of Melbourne, Parkville, VIC 3010, Australia</footnote>

<footnote n="74">⁷⁴Laboratoire de physique des 2 infinis Irène Joliot-Curie, CNRS Université Paris-Saclay, Bât. 100, Faculté des sciences, F-91405 Orsay Cedex, France</footnote>

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

<footnote n="76">⁷⁶Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote n="77">⁷⁷Kavli Institute for Particle Astrophysics and Cosmology and Department of Physics, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="78">⁷⁸Brookhaven National Laboratory, Bldg 510, Upton, NY 11973, USA</footnote>

<footnote n="79">⁷⁹Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA</footnote>

<footnote n="80">⁸⁰Instituto de Física, Pontificia Universidad Católica de Valparaíso, Casilla 4059, Valparaíso, Chile</footnote>

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

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

<footnote n="83">⁸³National Center for Supercomputing Applications, 1205 West Clark Str, Urbana, IL 61801, USA</footnote>
This paper has been typeset from a TEX/LATEX file prepared by the author.

<footnote n="84">⁸⁴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="85">⁸⁵Department of Physics, Cornell University, Ithaca, NY 14853, USA</footnote>

<footnote n="86">⁸⁶Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA</footnote>

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

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

<footnote n="89">⁸⁹NASA/Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>