Cosmological shocks around galaxy clusters: A coherent investigation with DES, SPT & ACT

We search for signatures of cosmological shocks in gas pressure profiles of galaxy clusters using the cluster catalogs 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 < 𝑀 200m / M ⊙ < 10 15 . 5 and 0 . 1 < 𝑧 < 2, and the total sky coverage of the maps is ≈ 15 , 000 deg 2 . We find a clear pressure deficit at 𝑅 / 𝑅 200m ≈ 1 . 1 in SZ profiles around both ACT and SPT clusters, estimated at 6 𝜎 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 behavior. 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 catalogs — 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.


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 onto a halo.The gravitational infall velocity of the cold gas will generically exceed the sound speed of the gas, especially for infall around massive halos, and this results in a high Mach number shock ( ∼ 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 et al. 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. 2020Zhang et al. , 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 et al. 2022a;Shin & Diemer 2023).However, it has remained difficult to infer observationally.
This process of shock heating generates a thermal nonequilibrium 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 largescale 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 et al. 2017;Hill et al. 2018;Pandey et al. 2019Pandey et al. , 2020;;Sánchez et al. 2023), with weak-lensing shears (e.g., Ma et al. 2015;Hojjati et al. 2017;Osato et al. 2018Osato et al. , 2020;;Shirasaki et al. 2020;Gatti et al. 2022;Pandey et al. 2022), or with X-ray luminosity (Shirasaki et al. 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 onto 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 behaviors.This phenomenology has been extensively studied in simulations over the past many decades.The first studies used nonradiative 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. 2014Hong et al. , 2015;;Schaal & Springel 2015).These were then followed by studies using simulations that include gas cooling and star formation (Vazza et al. 2009;Planelles & Quilis 2013;Lau et al. 2015;Nelson et al. 2016;Aung et al. 2021), and also include the effects of feedback from supernovae and active galactic nuclei (Kang et al. 2007;Vazza et al. 2013Vazza et al. , 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 et al. 2018;Zhang et al. 2019bZhang et al. , 2020Zhang et al. , 2021)).A number of works have also theoretically estimated the potential signal-to-noise of shocks from various surveys/instruments (e.g., Kocsis et al. 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 gaswhich has a low sound speed and is primarily found in low-density regions such as cosmic voids -onto the thermalized, bound gas component results in a shock of a high Mach number ( ∼ 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 clusteris oftentimes referred to as an accretion shock (e.g., Lau et al. 2015;Aung et al. 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. (2019bZhang et al. ( , 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 ≲ / 200m ≲ 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 et al. 2014;Diemer & Kravtsov 2014).This feature has been found in various datasets (Baxter et al. 2017;Chang et al. 2018;Shin et al. 2019;Adhikari et al. 2021;Shin et al. 2021) and its connection to cosmological shocks has been explored via both analytic calculations and simulations (Shi 2016;Aung et al. 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 et al. 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 et al. 2019a,b;Hurier et al. 2019;Pratt et al. 2021;Zhu et al. 2021).More general studies of gas thermodynamic profiles, without a specific focus on shocks, do not push beyond  ≳  500c (e.g., McDonald et al. 2014;Ghirardini et al. 2017;Romero et al. 2017Romero et al. , 2018;;Ghirardini et al. 2018), though some do exist (Planck Collaboration et al. 2013;Sayers et al. 2013Sayers et al. , 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 102 − 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 catalogs, and between measurements from different SZ maps (ACT and SPT), (iii) measure cluster profile outskirts for lower-mass clusters,  200m < 10 14.5 M ⊙ , (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 catalog from the Dark Energy Survey (DES) Year 3 dataset, and an SZ map from ACT Data Release (DR) 6.Both datasets 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 catalog, whereas A22 were limited to using a subset (≈ 25%) of the catalog that overlapped with the ACT DR4 map.
We organize this work as follows: in §2 we describe the survey datasets used in this work and in §3 our choices for the profile measurement procedure and the theoretical modelling.Our results on shocks are shown in §4 and their connections to other large-scale structure features are explored in §5.We conclude in §6.

DATA
We use data from three wide-field 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 -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 Section 2.2 and Section 2.3).The former choice is because the 10 ′ resolution of the Planck SZ map (which is an order-of-magnitude larger than the 1 ′ resolution of SPT and ACT) is a limiting factor in detecting shock features.The Planck cluster catalog (Planck Collaboration et al. 2016) also has significant overlap with the SPT and ACT catalogs used in this work; 45% of the 1093 Planck clusters are found within either the ACT or SPT footprints.
The clusters in our samples are labeled by their spherical overdensity mass,  200m , which is defined as, with  Δ = 200  (), where   () is the mean matter density of the Universe at a given epoch.The associated radius is denoted as  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  500c from the integrated tSZ emission around each cluster, while DES infers  500c from the cluster richness, where richness is the probabilistic number of satellite galaxies in the cluster.We then convert the  500c estimate into  200c and  200m using the concentration-mass relation from Diemer & Joyce (2019) and the publicly available routine from the COLOSSUS1 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.The tSZ amplitude is reported as the dimensionless  parameter, where   is the Boltzmann constant,   is the Thomson crosssection,    2 is the rest energy of an electron,   and   are the electron number density and temperature, respectively, and  is the physical line-of-sight distance.Thus  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 scatteringoff 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 ( e ≳ 5 keV) requires relativistic corrections to the procedure for making the SZ maps.These corrections, however, are ≲ 5% (Erler et al. 2018, see Figure 1) and are subdominant to the amplitudes of the features discussed in our work. 2The ntSZ effect can be generated by a cosmic ray electron population, but is a subdominant effect within  200c of the cluster, where cosmic rays make up ≲ 1% 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% level.

The Dark Energy Survey (DES) Year 3
DES Y3 is a 5000 deg 2 photometric survey of the southern sky in five bands ( ).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", , which is analogous to the number of red galaxies in the cluster.redMaPPer assigns each galaxy  a probability that it is a satellite of galaxy cluster .The richness of cluster  is then the sum of these probabilities.
This richness is used alongside a richness-mass relationwhich 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 richnessmass 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 observable-mass relation -for clusters with  ≥ 20.The observable-mass relation was in turn calibrated with targeted weaklensing measurements.Note that the catalogs we use have objects of lower richness ( ≈ 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 signal-to-noise ratio 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, /Δ, 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 catalog (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. 3Then, we also use the magnitude-limited lens galaxy catalog, Maglim (Porredon et al. 2021), to infer the density field in the DES footprint, from which we can estimate a cluster orientation based on largescale structure.This follows the methods of Lokken et al. (2022), and is discussed further in Section 5.1.Both datasets are part of the publicly available DES Y3 data release. 4that this is a factor two difference in an effect that contributes < 5% to the total signal. 3We 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. 4https://des.ncsa.illinois.edu/releases/y3a2

The South Pole Telescope (SPT) SZ Survey
SPT-SZ is a 2500 deg 2 survey of the southern sky at 95, 150, and 220 GHz, and was conducted using the South Pole Telescope (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 ′ , and is made using data from both SPT-SZ and the Planck 2015 data release; the former provides lower-noise measurements of the small scales, whereas the latter does the same for larger scales (multipoles ℓ ≲ 1000).The Planck data consists of the 100, 143, 217, and 353 GHz maps from the High Frequency Instrument (HFI).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 linear combination 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% 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 catalog from this data contains 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 catalog for our work, where the mass is estimated via a joint modeling of SZ, X-ray, and weak lensing measurements.Both the map and the cluster catalog are publicly available. 5Our masses come from the M500 column and signal-tonoise ratio (SNR) from the XI column.This dataset is the exact same as the SPT-SZ data used in A22.

Atacama Cosmology Telescope (ACT) data releases 4, 5, and 6
The ACT data covers 90, 150, and 220 GHz frequencies, and the maps from data release (DR) 6 cover ≈ 13,000 deg 2 of the sky (after applying the relevant masks; see discussion below).The SZ map (Coulton et al. 2023) has a resolution of 1.6 ′ , and makes use of data from both ACT and the Planck NPIPE data release (Planck Collaboration et al. 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 (NILC) 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 Pixell6 , 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.16).The color tones of the points show log 10 SNR, the signal-to-noise ratio of each cluster detection, with lighter colors indicating a higher SNR.The mean redshift and mass of the different samples are listed in Table 1.
We also use the ≈ 4200 clusters from ACT DR5 7 catalog (Hilton et al. 2021), which covers the same area as the ACT DR6 map.Note that only the subset of the ACT DR5 catalog, that corresponded to the 2000 deg 2 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 catalog 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.

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

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 7 https://lambda.gsfc.nasa.gov/product/act/actpol_dr5_szcluster_catalog_info.html 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 ⟨⟩ profile in 50 logarithmically spaced radial bins in the range  ∈ [0.1, 20]  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 catalog, 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. 8reviously, 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, 3D space but appear close in projected, 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  > 20 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 lineof-sight separation is not relevant.Thus, selecting clusters where the line-of-sight separation is greater than 20 200m isolates such interlopers.The choice of 20 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 (ΛCDM) cosmology with Ω m = 0.3 and ℎ = 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 signalto-noise (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 -th jackknife sample of the stacked profile can be written as where   is the SNR per cluster used in the weighted average,   () is 1 if the datapoint for radius  in cluster  is unmasked and 0 otherwise,  cl is the total number of clusters.In this notation, ⟨⟩  is the mean profile of the sample with cluster  removed, and   is the individual profile measurement from cluster .The variance on the mean profile is then given by, where ⟨ ⟩ is the mean of the distribution of jackknife estimates computed in Equation (3).Note that Equation ( 5) has an additional factor of  −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 et al. 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  ln  = 0.16, which is 1.5 times the logarithmic bin width, Δ ln  ≈ 0.11.All profiles are smoothed by this scale, and we present results only for the range 0.3 < / 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, where  is an arbitrary function of , and ℎ = Δ ln  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 √  − 1 is applied to convert the measured uncertainty to the unbiased uncertainty, and this is analogous to the extra  − 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, where  and  index over the different radial bins,  ′ , is the logderivative of the mean profile in the  th bin for the  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 Figure 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 √  − 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  200m (as determined from the catalogs) is 15%, and so the uncertainty in  200m is around 5%.This is tolerable as it increases the total uncertainty in the estimated feature location by < 2%.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.

Modeling 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- 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 halos as described in Vikram et al. (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(Cui et al. , 2022)).The entire model is implemented in the Core Cosmology Library (CCL) open-source python package9 (Chisari et al. 2019) and is public. 10e begin by representing the 3D, halo-pressure crosscorrelation function as a composition of the one-halo and two-halo components, where  are the correlation functions,  is comoving distance, and  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), where  0 ,   , , , and  are the fit parameters calibrated from hydrodynamical simulations,  = / 200c is the distance in units of cluster radius, and  200c is the thermal pressure expectation from self-similar evolution, Equation ( 10) accounts for deviations from self-similar evolution via the calibrated mass and redshift dependencies of the parameters  0 ,   , and .The model also includes the effects of non-thermal pressure support within halos -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 with both cluster redshift, , and cluster mass,  200c .
The calibration matches simulations within < 10% in the one-halo regime (Battaglia et al. 2012, see their Figure 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 smallscales 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,   , whereas the profiles of Battaglia et al. (2012) are calibrated to the total gas pressure, .We convert between them as with  = 0.24 being the primordial helium mass fraction.This provides our one-halo term, 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 crosspower spectrum,  two−halo ℎ,  , is written as, where  is the mass of the halo we are computing the halo-pressure correlation for,  ′ is the mass of a neighbouring halo contributing to the two-halo term,  lin (, ) is the linear, matter density power spectrum at redshift , / ′ is the mass function of neighbouring halos, and (, ) and ( ′ , ) are the linear bias factors for the target halo and neighboring halos, respectively.The mass function model comes from Tinker et al. (2008) and the linear halo bias model from Tinker et al. (2010).The term   (,  ′ , ) is the Fourier transform of the pressure profile about the neighboring halo which, under the assumption of spherical symmetry, is computed as, where   is the electron pressure profile.The halo-pressure twopoint cross-correlation is obtained as the inverse Fourier transform of the cross-power spectrum, (, , ). ( The terms shown in equations ( 13) and ( 16) can be combined according to Equation ( 9) to get the total halo model,  ℎ,  .We have thus far described the real-space 3D pressure, whereas the Compton- parameter is the integrated (or projected) pressure along the line of sight.The halo- correlation is therefore obtained by a projection integral, where   is the Thomson scattering cross-section,    2 is the rest mass energy of the electron, and  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, where  0 is the zeroth-order Bessel function.We then multiply  ℓ by the Fourier-space smoothing function for the given survey of interest and then perform an inverse-harmonic transform, with the smoothing function  ℓ given as where  FWHM =  FWHM / √ 8 ln 2, with  FWHM = 1.25 ′ ( FWHM = 1.6 ′ ) 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,  smooth ℎ, , for each individual cluster in our catalog, 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 miscentering 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 where  is the uncertainty in the log-derivative measurement.The quantity  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, where C −1 is the inverse of covariance matrix for the log-derivative, accounting for the Hartlap factor (Hartlap et al. 2007) as, where  jk are the number of jackknife samples (more than 500 for almost all samples),  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 Δ log 10  = 0.1 of the location of the feature.Once the  2 is computed, we quote the total signal-to-noise of a feature, as following the definition of Secco et al. (2022b, see their Equation C15), with  dof = 5 as mentioned above.This definition of signalto-noise improves on that used in A22 as it is more robust to noise fluctuations and binning choices.

Miscentering model for optically selected clusters
An additional component to our theoretical model, in comparison to that of A22, is the impact of cluster miscentering.For SZ-selected clusters, the offset between the cluster center and the true center (called "miscentering") is negligible when compared to the radial scale of features we study, which are ∼  200m .When using optically selected clusters, however, the optically determined center can be significantly offset from the center of the gas distribution (Sehgal et al. 2013;Zhang et al. 2019a;Bleem et al. 2020). 11The impact of miscentering in the profile is to transfer power from small scales to large scales.The total observed profile, with miscentering, can be modelled as, where  miscen is the fraction of miscentered objects and  true ( miscen ) is the profile of correctly centered (miscentered) clusters.For a given miscentering offset,  mis , the average miscentered profile is, and the total model is obtained by marginalizing over the distribution of possible offsets, Following previous works (e.g., Baxter et al. 2017;Chang et al. 2018;Shin et al. 2019Shin et al. , 2021)), we assume the offsets follow a Rayleigh distribution, 11 SZ-selected clusters also incur a noise-induced miscentering effect, with a scale of  miscen = √︃  2 500c +  2 beam /SNR.For  ≳  200m , the miscentering scale is at/below the bin width and is negligible as our features of interest span multiple bins.The average SZ-selected cluster ( 200m ≈ 10 14.8 M ⊙ and  ≈ 0.6) has  miscen = 0.3 ′ , while the same for the average optically selected cluster ( 200m ≈ 10 14.6 M ⊙ and  ≈ 0.4) is factors of 5 to 10 larger (Zhang et al. 2019a;Bleem et al. 2020).
where  is the cluster richness.The free parameters of this model are  miscen and  miscen which set the fraction of miscentered objects, and the amplitude of the miscentering offset, respectively.The impact of miscentering -and the choice of the parameter values -for DES cluster profile model is discussed in Section 4.1 and further in Appendix A.

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 halos, in Section 4.3.We will use the format CATALOG x MAP as a shorthand reference for measurements for a given cluster catalog using a given SZ map (e.g., SPTxSPT, DESxACT).All bands show 68% uncertainties estimated via jackknife resampling of the profiles.As for the detection significance, we show  in the figures but quote  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  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 behaviors in the one-halo and two-halo regimes.Shocks occur at the transition between the bound halo component (one-halo term) and the surrounding largescale structure (two-halo term), so studying shock-induced features also requires studying these regimes.Thus, some of our discussions below will include behaviors of the one-halo and two-halo terms, as changes in these terms affect the overall shape of the halo profile.

Measurements from fiducial cluster samples
In Figure 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 analyzed 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 (/ 200m ≲ 0.5) and also in the far outskirts (/ 200m ≳ 5), but has significant deviations at / 200m ≈ 1, and potentially also at / 200m ≈ 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  21).The theoretical prediction (dashed lines), which is a sum of one-halo and two-halo contributions (each shown as gray dotted and dashed-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  200m > 10 14.5 M ⊙ .For the SZ-selected samples, the derivative is lower at  200m than the theory curve, consistent with A22.The behavior for optically selected clusters is less clear due to potential inaccuracies in the theoretical model, such as the miscentering 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 gray 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 miscentering effects included; including the miscentering (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 Figure D1.
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 Figure 2) use simulations specialized to model the electron-ion temperature differences and show that this effect causes a deficit in the cluster tSZ profiles 12 , while Avestruz et al.
(2015, see their Figure 1) do the same but focus on the 3D cluster temperature profiles.This pressure deficit feature would not be present in most cosmological hydrodynamical simulations as they a priori assume local thermal equilibrium between electrons and ions.
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,  ≈  200m , and are instead limited to much smaller radii where the higher number densities allow the ion and electrons to quickly achieve temperature equilibrium.
We will henceforth refer to the pressure deficit as a shock feature and denote its location the shock radius,  sh .
As Figure 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.Given the new, more robust definition of signal-to-noise 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 compared to the estimate of 3.1 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 / 200m = 1.09 ± 0.08 and / 200m = 1.16 ± 0.04, respectively.The minima in the log-derivatives are consistent as well, with 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 artifact 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 Figure 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 catalogs which have  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 Figure 2 also present the quantity , 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,  takes a maximum value at / 200m ≈ 1, corresponding to the pressure deficit.In optically selected clusters, which we will discuss below, the maximum values of  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 miscentering 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 / 200m ≈ 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 datasets.2013) finds that the 3D pressure profiles have a deficit, relative to the theoretical predictions of Battaglia et al. (2012), for  ≳  500c ( ≳ 0.3 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  ≈  200c = 0.5 200m using Suzaku X-ray data.Hou et al. (2023) study the radio emission around galaxy clusters and find a signal at  = 2.5 500c ≈  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  ≈  200m , consistent with the shock radius  sh .
The right panels of Figure 2 show, for the first time, the outskirts of SZ profiles for optically selected clusters.We have placed a mass cut of  200m > 10 14.5 M ⊙ (where  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 Figure 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 ≈ 30% lower normalization than those of the SZ-selected catalogs, 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  200m , which is inbetween the one-halo and two-halo regime, the profile for the DESxACT measurement has a minimum logderivative ( d ln d ln = −3.1 ± 0.15) that is more negative than that of the DESxSPT measurement ( d ln d ln = −2.7±0.2), with a significance of 1.6.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 Figure 3 below that the SPT and ACT maps provide statistically indistinguishable results across the full range of scales considered in this work.
Looking at the DESxACT and the ACTxACT results, we see the location of the log-derivative minima is consistent at 0.2, while the depth of minima is deeper in ACTxACT at 2.The comparison of DESxSPT and SPTxSPT is similar, where the location of the logderivative minima is consistent while the depth deviates at 2.4.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-analyze 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 Figure 2 is a qualitatively good match to the data across the whole range of presented scales.The prediction for the DESxSPT and DESxACT measurements closely overlap one another.This model includes the miscentering effects described in Section 3.2.1, using values of  miscen = 0.9 and  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 onehalo regime, near the cluster core.The preferred values for  miscen and  miscen are both near the 3 − 4 upper limit of the miscentering 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  miscen is within 1 of the estimate from Bleem et al. (2020, see their Table 6), which is based on a SPT-DES matched cluster sample.Figure 2 shows the theory matches the data better (in the 1-halo regime) when we include this miscentering effect, and the dotted lines show the theory without such effects.
In Appendix A, we discuss how the profiles and log-derivatives depend on miscentering parameters.We emphasize that in our work we only focus on results from optically selected clusters that are insensitive to the choice of miscentering 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) SZselected 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 miscentering, contamination, or mass estimation errors) causes the feature to be suppressed.In Section 4.2 below, 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 The two measurements are consistent across the whole range of scales, with  2 / dof = 1.1 and  = 0.14, validating that the datasets and map-making procedures of the two surveys are consistent in both high and low signalto-noise regimes.The SPT and ACT datasets are independently calibrated and mapped, and the statistical consistency in the measurements above is determined at the ≈ 1% level given the precise measurements in the high signal-to-noise regime.
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, ⟨ 200m ⟩ = 10 14.85 M ⊙ , 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 Figure 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  = 2.We verify in Appendix E that our results are unchanged if we use SZ maps that minimize/deproject the CIB signal.

Sensitivity to map-making and cluster selection
The results in Figure 2 show a difference between the profiles of SZ-selected and optically selected clusters -the former sees a clear pressure deficit at / 200m ≈ 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 Figure 3, we test an aspect of the former effect, namely noise-based SZ-selection effects. 13These 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 ( = 669 clusters), and then by measuring the subsample's average SZ profile using either the SPT map or the ACT map.We find consistency ( 2 / dof = 1.1 with  = 0.14) in both the profiles and the log-derivatives of the two measurements.While this implies 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. 14 It validates the maps' consistency in both the high signal-to-noise regime at the location of massive clusters, as well as in the noisedominated, low-signal-to-noise 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. 15We first remove all ACT clus- 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 Figure 4. 14 A similar analysis using all SPT clusters in both footprints finds  2 / dof = 1.06 with  = 0.36.However, the profile measurement uncertainties are broader as the SPT cluster sample size is half that of ACT. 15 The zeroth-order effect of SZ and optical selection on the cluster sample is in its mass and redshift distributions (see Figure 1).The reweighting accounts for these selection effects, and thus any further differences in the ters with redshifts/masses outside the ranges of the DES sample.We therefore use all ACT clusters within 0.1 <  < 0.8 and 13.7 < log 10  200m < 15.35 to create the subsample used in this analysis, which has  = 3392 clusters.The reweighting is then done by computing the weighted counts of clusters in a 2D grid of  200m and , 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 where   is a delta function -with values of 0 or 1 -that denotes whether cluster  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  200m −  grid cell it is associated with.We have checked that our results do not change if we use a 20x20 grid instead.The final weight of the DES cluster is, which uses the original signal-to-noise 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.
Figure 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 / 200m ≈ 1 than the theoretical predictionand 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).
Figure 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 < / 200m < 10, the two profiles are consistent( 2 / dof = 1.4 with  = 0.1).Under this reweighting, the weighted mean mass of the DES sample increases from ⟨ 200m ⟩ = 10 14.6 M ⊙ → 10 14.85 M ⊙ , which is a fractional change of 80%, while the mean redshift is left unchanged at ⟨⟩ = 0.46 (see Table 4).The depth of the minima is now consistent across the two samples, whereas it was inconsistent at the 2 level for the fiducial ACT and DES cluster samples (Figure 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 Figure 3 and Figure 4 imply that -for a mass and redshift range corresponding to clusters in SZ surveys (see Figure 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 SPTxSPT and ACTxACT measurements, as reweighted profiles can be attributed to selection effects beyond those on the cluster samples' mass and redshift distributions.
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 (Figure 2) could possibly be attributed to the clusters in the lower mass end of the sample.We explore the behavior of such systems further in the following section.

Towards galaxy groups
The profiles of massive clusters have many observational constraints, especially near the cluster core (/ 200m ≲ 0.5, see for example, McDonald et al. 2014;Ghirardini et al. 2017;Romero et al. 2017Romero et al. , 2018;;Ghirardini et al. 2018), and the further outskirts have only been recently explored observationally (Planck Collaboration et al. 2013;Sayers et al. 2013Sayers et al. , 2016;;Amodeo et al. 2021;Schaan et al. 2021;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 (Figure 2) was not identified in previous cross-correlation works but was easily identified in A22 by measuring individual profiles.SZselected halo samples are ideal for studying massive, cluster-scale halos but are not viable for probing lower masses.Here, we use the DES redMaPPer sample to obtain a catalog of lower mass objects ( 200m ≳ 10 13.8 M ⊙ ) and measure their SZ profiles across a wide range of scales.
In Figure 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 Figure 5, we see the pressure deficit exists for all mass bins, at 3.8, 2.5, and 2.9 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 Figure 5, however, we find that splitting by mass -which is directly proportional to signalto-noise -does not notably change the significance of the deficit.In Appendix B, we also redo this test by splitting directly on SNR instead of  200m , and find consistent results.
The other two columns in Figure 5 (left and middle) show the SZ profile for three mass bins of optically selected clusters.The mass range  200m > 10 14.5 M ⊙ corresponds to  > 30, while the range 10 14 M ⊙ <  200m < 10 14.5 M ⊙ corresponds to 15 <  < 30, and finally, 10 13.8 M ⊙ <  200m < 10 14 M ⊙ corresponds to 10 <  < 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 Figure 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

R/R 200m
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 ( 200m < 10 14.5 M ⊙ ) 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 Figure 2.
lower mass objects the log-derivatives are closer to zero in the onehalo term (/ 200m ≲ 1) and similar to the high mass bin results for the two-halo term (/ 200m ≳ 4).The log-derivative minima in each mass bin are found at similar radii of / 200m ≈ 1.1 (see Table 1).
The theoretical model also significantly deviates from the measurements in these two lower mass bins.For  200m ∈ [10 13.8 , 10 14 ] M ⊙ , the deviation is a factor of ≈ 5 in the halo core.In the intermediate mass bin,  200m ∈ [10 14 , 10 14.5 ] M ⊙ , it is a factor of ≈ 2 and is also consistent with previous analyses in this intermediate mass range.Saro et al. (2017, see their Table 1 and  Figures 5/6) found a factor of ≈ 2 difference when measuring the integrated SZ effect around clusters from the DES Science Verification data, while Planck Collaboration et al. (2011, see their Figure 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 behavior 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 halos, 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  ∝  ℎ () ∝  0.5 for the halo bias model of Tinker et al. (2010) at high halo masses, whereas the one-halo term scales as  ∝  5/3 .The deviations in Figure 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  ∝  5/3 200m -would be 50% to 150% in  200m .Myles et al. (2021) show the richness bias due to contamination is ≈ 20% for clusters of 5 <  < 20 (see their Section 4.3), and also that richness depends on halo mass as  ∝  1.0 (see their Section 4.5).This implies the mass bias is 20% × 1.0 = 20%, lower than the required values of 50% to 150% 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 ≈ 30% changes in the final mass estimate (Costanzi et al. 2021, see their Equations 16 and 17).
The second effect,  −  relation deviations, are deviations in the pressure profile model for lower mass halos.This work uses the model of Battaglia et al. (2012), and while it is accurate for higher mass halos (e.g., see Figure 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 behaviors 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 Figure 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 SZmass scaling relation of the selected cluster sample.The effect has been detailed in the analytical work of Evrard et al. (2014, see their Figure 4 for an example).The scaling relation for the optically selected sample is now written as where − 1 is the slope of the halo mass function at a chosen mass scale,   is the slope of the richness-mass relation, and cov(. ..) 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 ⟨ln | ln ⟩ can be higher (lower) than ⟨ln | ln  200m (ln )⟩, 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 ≲  ≲ 0.5 at 95% 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 Figure 5) while that of the stellar mass and richness are correlated (Anbajagane et al. 2020, see their Figure 7).A negative correlation/covariance suppresses the SZ signal of richness-selected clusters, which could cause the observed suppression.For conservative values of − 1 = 1.5,   = 1.5,  ln  = 0.3,  ln  = 0.8,  = −0.6, 16we find the bias is at most 30%.Thus, this effect cannot be the main cause of the behaviors found in Figure 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.

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 16  1 is the slope at a pivot mass of  200m ∼ 10 14 M ⊙ , computed using the halo mass function of Tinker et al. (2008)

Inter. Minor
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.
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 behaviors 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 behaviors of these two regimes.

Connections to filaments
We have discussed previously that cosmological shocks form from the accretion of collisional matter onto 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 et al. 2021, see their Figure 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 shock feature further into the cluster core and/or completely destroy it (Zhang et al. 2020, see their Figure 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 subsection of the image.The geometry of this split is shown in Figure 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 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 catalog, and from a 2D Gaussian fit to the SZ image cutout 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.3 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.
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  < 0.5 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  < 0.3 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 catalog, and also using the largescale 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, , where  is the overdensity field of the galaxy number density and   are the projected coordinates.The Hessian is then diagonalized to find the orientation of the major axis.In this work,  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 half-max (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.Figure 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.
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 (/ 200m ≳ 5).In the far outskirts, / 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, / 200m ≈ 1, the pressure profile has a steeper derivative along the minor axis.The profiles also show a plateauing feature at / 200m ≈ 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 et al. 2021, see their Figure 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 ACTxACT and SPTxSPT results in Figure 2. We have verified that all shock behaviors 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). 17 The consistency of the anisotropic profiles on small-scales is not an artifact of angular resolution limits.For the average cluster, the scale  = 0.5 200m is a factor of ≈ 2.5 larger than the FWHM of the smoothing for the maps.Furthermore, as we will show below, in Figure 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 other orientation estimates are sensitive to, and the direction towards large-scale structure can change noticeably across different scales (Lokken et al. 2022, see their Figure 16).We find that orienting by the BCG shape impacts only the one-halo regime.The observed striation in Figure 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 behavior 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 / 200m ≲ 0.8.
In summary, we observe potential behaviors of the logderivative 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 et al. (2021), who showed the accretion shock is elliptical and pointed along the largescale 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. 17The 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.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% 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  sp / sh = 1.17 ± 0.20, meaning the two projected radii are within 0.9 of each other.The dashed line is the prediction of the SZ profile log-derivative for this cluster sample.

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 et al. 2014;More et al. 2015;Mansfield et al. 2017;Aung et al. 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 et al. 2019).A22 performed the first comparison of the splashback and pressure deficit features using the SPT-SZ dataset.In this work, we supplement this result with a complementary analysis using the ACT data.A subset of the ACT cluster catalog 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 catalog selections as Shin et al. (2021), taking all ACT DR5 clusters within the DES Y3 footprint and with 0.1 <  < 0.7.We then measure the log-derivative of the average SZ profile for this subsample, and present the result in Figure 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 Figure 8, the minimum in the SZ log-derivative, corresponding to the pressure deficit and which we denote as the shock radius  sh , coincides with the splashback radius,  sp .In particular, we find the ratio to be  sp / sh = 1.17 ± 0.20. (33) The estimate for  sh used above is presented in Table 1 under the name "ACT x ACT ( sp comparison)".Our results above show the splashback radius and shock radius are consistent within 0.9.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 et al. 2021;Zhang et al. 2021;Baxter et al. 2021), this is for the merger-accelerated accretion shock which is expected to be at  >  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.

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 datasets -the Dark Energy Survey Year 3, the South Pole Telescope SZ survey, and the Atacama Cosmology Telescope 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: • Consistent with A22, there is a pressure deficit at / 200m ≈ 1.1 detected at 2. 7𝜎 and 6.1𝜎 in SPT and ACT,respectively (Figure 2).This feature is consistent with a shock-driven thermal nonequilibrium between electrons and ions.We do not quote a detection significance for DES clusters given uncertainties in the theoretical modelling (see Section 4.1).
• The SZ maps from SPT and ACT are consistent in both high and low signal-to-noise regimes (Figure 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 < / 200m < 10.
• We construct ACT and DES subsamples with similar mass and redshift distributions and find their mean SZ profiles to be consistent (Figure 4).This implies that for clusters of a higher mass, ⟨ 200m ⟩ = 10 14.85 M ⊙ , 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 ACTxACT and SPTxSPT measurements.
• For optically selected clusters of lower masses,  200m < 10 14.5 M ⊙ , the radial location of the log-derivative minima are consistent at / 200m = 1.1 while the depth becomes shallower with decreasing mass (Figure 5).
• The SZ profiles measured around group-scale halos also differ significantly from the theoretical model in the one-halo regime, and are consistent with the model for the two-halo regime (Figure 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.
• 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 (Figure 7).
• The location of the pressure deficit,  sh , is consistent with the splashback radius measured with galaxy number density profiles in Shin et al. (2021).The ratio is  sp / sh = 1.17 ± 0.20, and this consistency between shock and splashback radii further signifies the variety of dynamical processes happening at  ≈  200m .
Our work, through the use of multiple independent datasets, 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(Zhang et al. , 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 catalog and/or lower-noise SZ maps, or perform selection cuts on the current catalogs (particulary related to cluster relaxation) that can maximize the signal-to-noise of this feature.
Moving forward there are still additional ongoing and future surveys/datasets that could be used for this work -such as SPT- From left to right the columns show: (i) the sample name, denoted as "cluster catalog 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 miscentering 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.
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 datasets can also provide maps with a depth comparable to the upcoming CMB-S4 experiment, and the corresponding cluster catalogs -such as the SPT Megadeep catalog (Kornoelje et.al, in prep.)-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 dataset.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 et al. 2022b).Techniques have also been developed to extract such profile correlations in a data-driven manner, with minimal assumptions, such as Gaussian processes (Farahi et al. 2021) and local linear regression (Farahi et al. 2022a).
Thus, there are many synergistic opportunities for crosscorrelating the different types of datasets -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.
Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.

APPENDIX A: IMPACT OF MISCENTERING ON PROFILES
As we discussed previously, our theoretical model for the SZ profiles of optically selected clusters depends on the miscentering model parameters assumed for the cluster sample.In Figure A1 we show the model for the DES clusters' SZ profile as we vary the amplitude of miscentering,  miscen , and the fraction of miscentered clusters,  miscen , defined in Section 3.2.1.We do not have an external, calibrated constraint for the miscentering 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 Figure 2. In practice, we do this by making predictions in a 5x5 grid of parameter values, and find  miscen = 0.9 and  miscen = 0.4 to be the best combination.As was noted before, both values are near the 3 − 4 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  miscen is within 1 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  miscen ≈ 1.5.Given these potential limitations of the implemented miscentering 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 miscentering 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 miscentering as the deviations are significantly larger than those from miscentering effects alone.
Figure A1 shows that for the high mass sample (left panels), the variation in  miscen changes the location of the logderivative minimum from  min = 0.7 200m → 1.2 200m as we vary  miscen = 0.1 → 1.2.The minimum value of the log-derivative goes from −2.5 → −3.0 as we vary  miscen = 0.1 → 0.8.In particular, the result for  miscen = 0.8 and  miscen = 0.9 appears to replicate a shock-esque feature at  ≈  200m .However, this is not an indication that shock features can be explained by miscentering.For SZ-selected clusters, the miscentering is much smaller than than the values being considered in Figure A1.While the predicted profile for large miscentering 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 miscentering is the cause of the deficit feature, and is instead evidence of the miscentering model's ability to capture steep drops in the pressure profile.The average SZ profiles of SZ-selected clusters measured on their respective SZ maps, with the samples split by their signal-to-noise ratio.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.
We then fit the model in Equation (C1) to the measured profiles and obtain constraints on the three parameters, ,  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, −10 <  < 0 0.8 <  < 2.5 (C2) 0.05 <  < 0.9.
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  < 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  > 0.6 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  2 for the log-likelihood: where  is the covariance matrix of the profile,  obs is the measured profile, and  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 catalogs, each measured on their corresponding SZ maps.The parameters corresponding to this fit are shown in Figure C1, while the fits and the data are compared in Figure 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 Figure 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 across the two samples.We also list the  2 of the fit from the original halo model and the modified halo model.For both datasets, the  2 is noticeably improved, and we can also see this visually in the fits of Figure 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  = 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  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 Figure 2 for these samples is consistent.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
Figure D1 presents the correlation matrix of the ACTxACT logderivative 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  2 shown in Table 1.The top panel shows the log-derivative measurement, now in discrete points corresponding to the binning, corresponding to the ACTxACT measurement of Figure 2.

APPENDIX E: IMPACT OF COSMIC INFRARED BACKGROUND CONTAMINATION
We explicitly test the impact of the cosmic infrared background (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); Coulton et al. (2023) for details on the deprojection procedure of SPT and ACT, respectively.In Figure 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 multiple CIB-deprojected maps, and we use the fiducial one (Coulton et al. 2023, see their Section C.3 and Equation 18).

Figure 1 .
Figure 1.The mass-redshift plane of the cluster samples from SPT, ACT, and DES used in this work.The Planck catalog 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  = 5000.The 1D distributions are estimated using the full samples.The SPT and ACT samples have similar redshift distributions, with a median of  ≈ 0.55, while DES Y3 is limited to 0.1 <  < 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 Equation16).The color tones of the points show log 10 SNR, the signal-to-noise ratio of each cluster detection, with lighter colors indicating a higher SNR.The mean redshift and mass of the different samples are listed in Table1.

Figure 2 .
Figure2.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 gray dotted and dashed-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  200m > 10 14.5 M ⊙ .For the SZ-selected samples, the derivative is lower at  200m than the theory curve, consistent with A22.The behavior for optically selected clusters is less clear due to potential inaccuracies in the theoretical model, such as the miscentering 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 Table1.The gray 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 miscentering effects included; including the miscentering (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 FigureD1.
Hurier et al. (2019) see a sharp decrease in pressure at  = 3 500c ≈  200m for a single cluster in the Planck data.Pratt et al. (2021) also use Planck data and find an excess in pressure at  = 2 500c ≈ 0.7 200m for a set of ten, low-redshift galaxy groups.The analysis of Planck Collaboration et al. (

Figure 3 .
Figure3.The average SZ profile of an ACT cluster subsample ( = 669 clusters) measured using either the ACT map or SPT map.The subsample is defined as all clusters whose centers lie in both the ACT and SPT footprints.The two measurements are consistent across the whole range of scales, with  2 / dof = 1.1 and  = 0.14, validating that the datasets and map-making procedures of the two surveys are consistent in both high and low signalto-noise regimes.The SPT and ACT datasets are independently calibrated and mapped, and the statistical consistency in the measurements above is determined at the ≈ 1% level given the precise measurements in the high signal-to-noise regime.

Figure 4 .
Figure 4.The average SZ profiles of two different cluster samples: an SZselected 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 <  < 0.8, and then reweight the DES clusters to match the ACT subsample's  200m −  distribution.See Section 4.2 for details.The two profiles are consistent ( 2 / dof = 1.1 with  = 0.14), suggesting that in this mass/redshift range there are no SZ or optical selection effects that generate/suppress the pressure deficit feature.

Figure 8 .
Figure8.The log-derivative of the average SZ profile around ACT clusters with 0.15 <  < 0.7 and whose centers are in the DES footprint.We perform this selection so as to use the same ACT DR5 subsample asShin 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% 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  sp / sh = 1.17 ± 0.20, meaning the two projected radii are within 0.9 of each other.The dashed line is the prediction of the SZ profile log-derivative for this cluster sample.
Figure B1.The average SZ profiles of SZ-selected clusters measured on their respective SZ maps, with the samples split by their signal-to-noise ratio.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 Table1.

Figure C1 .
Figure C1.The constraints on the parameters of Equation (C1) obtained from the SPT or ACT data.The values are listed in TableC1.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 Figure2for these samples is consistent.

Figure C2 .
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.

Table 1 .
A summary of the numerical results presented in this work.All uncertainties are ±1 estimates.