# Shocks in the Stacked Sunyaev-Zel’dovich Profiles of Clusters II:
Measurements from SPT-SZ + Planck Compton-$y$ Map

D. Anbajagane(  ),1, 2 C. Chang,1, 2 B. Jain,3 S. Adhikari,1, 2 E. J. Baxter,4
B. A. Benson,5, 2, 1 L. E. Bleem,6, 2 S. Bocquet,7 M. S. Calzadilla,8 J. E.
Carlstrom,2, 1, 9, 6, 10C. L. Chang,6, 2, 1 R. Chown,11 T. M. Crawford,1, 2 A.
T. Crites,12, 2, 1 W. Cui,13 T. de Haan,14, 15L. Di Mascolo,16, 17, 18 M. A.
Dobbs,19, 20 W. B. Everett,21 E. M. George,22, 15 S. Grandis,7, 23 N. W.
Halverson,21, 24 G. P. Holder,25, 26 W. L. Holzapfel,15 J. D. Hrubes,27 A. T.
Lee,15, 29D. Luong-Van,27 M. A. McDonald,8 J. J. McMahon,2, 1, 9, 10 S. S.
Meyer,2, 1, 9, 10 M. Millea,15L. M. Mocanu,2, 1 J. J. Mohr,7, 23, 30 T.
Natoli,2, 1 Y. Omori,2, 1, 31, 32 S. Padin,33, 2, 1 C. Pryke,34C. L.
Reichardt,36 J. E. Ruhl,36 A. Saro16, 17, 18 K. K. Schaffer,37, 2, 10 E.
Shirokoff,2, 1, 15Z. Staniszewski,38, 36 A. A. Stark,39 J. D. Vieira,25, 26 R.
Williamson38, 2, 1

1 Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL
60637, USA
2 Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL
60637, USA
3 Department of Physics and Astronomy, Center for Particle Cosmology,
University of Pennsylvania, Philadelphia, PA 19104, USA
4 Institute for Astronomy, University of Hawai’i, 2680 Woodlawn Drive,
Honolulu, HI 96822, USA
5 Fermi National Accelerator Laboratory, MS209, P.O. Box 500, Batavia, IL
60510
6 High Energy Physics Division, Argonne National Laboratory, Argonne, IL, USA
60439
7 Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679
Munich, Germany
8 Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology, Cambridge, MA 02139, USA
9 Department of Physics, University of Chicago, Chicago, IL, USA 60637
10 Enrico Fermi Institute, University of Chicago, Chicago, IL, USA 60637
11 Department of Physics and Astronomy, McMaster University, 1280 Main St. W.,
Hamilton, ON L8S 4L8, Canada
12 Department of Astronomy & Astrophysics, University of Toronto, 50 St George
St, Toronto, ON, M5S 3H4, Canada
13 Institute for Astronomy, University of Edinburgh, Royal Observatory, EH9
3HJ Edinburgh, United Kingdom
14 High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki
305-0801, Japan
15 Department of Physics, University of California, Berkeley, CA, USA 94720
16 Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo
11, Trieste 34131, Italy
17 INAF - Osservatorio Astronomico di Trieste, via Tiepolo 11, Trieste 34131,
Italy
18 IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2,
34014 Trieste, Italy
19 Department of Physics and McGill Space Institute, McGill University,
Montreal, Quebec H3A 2T8, Canada
20 Canadian Institute for Advanced Research, CIFAR Program in Cosmology and
Gravity, Toronto, ON, M5G 1Z8, Canada
21 Center for Astrophysics and Space Astronomy, Department of Astrophysical
and Planetary Sciences, University of Colorado, Boulder, CO, 80309
22 European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748 Garching,
Germany
23 Excellence Cluster ORIGINS, Boltzmannstr. 2, 85748 Garching, Germany
24 Department of Physics, University of Colorado, Boulder, CO, 80309
25 Astronomy Department, University of Illinois at Urbana-Champaign, 1002 W.
Green Street, Urbana, IL 61801, USA
26 Department of Physics, University of Illinois Urbana-Champaign, 1110 W.
Green Street, Urbana, IL 61801, USA
27 University of Chicago, Chicago, IL, USA 60637
29 Physics Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA
94720
30 Max-Planck-Institut für extraterrestrische Physik, 85748 Garching, Germany
31 Kavli Institute for Particle Astrophysics and Cosmology, Stanford
University, 452 Lomita Mall, Stanford, CA 94305
32 Dept. of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA
94305
33 California Institute of Technology, Pasadena, CA, USA 91125
34 Department of Physics, University of Minnesota, Minneapolis, MN, USA 55455
35 School of Physics, University of Melbourne, Parkville, VIC 3010, Australia
36 Physics Department, Center for Education and Research in Cosmology and
Astrophysics, Case Western Reserve University,Cleveland, OH, USA 44106
37 Liberal Arts Department, School of the Art Institute of Chicago, Chicago,
IL, USA 60603
38 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA
91109, USA
39 Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA
Corresponding author email<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We search for the signature of cosmological shocks in stacked gas pressure
profiles of galaxy clusters using data from the South Pole Telescope (SPT).
Specifically, we stack the latest Compton-$y$ maps from the 2500 deg2 SPT-SZ
survey on the locations of clusters identified in that same dataset. The
sample contains 516 clusters with mean mass $\langle M_{\rm
200m}\rangle=10^{14.9}\,{\rm M}_{\odot}$ and redshift $\langle z\rangle=0.55$.
We analyze in parallel a set of zoom-in hydrodynamical simulations from The
Three Hundred project. The SPT-SZ data show two features: (i) a pressure
deficit at $R/R_{\rm 200m}=1.08\pm 0.09$, measured at $3.1\sigma$ significance
and not observed in the simulations, and; (ii) a sharp decrease in pressure at
$R/R_{\rm 200m}=4.58\pm 1.24$ at $2.0\sigma$ significance. The pressure
deficit is qualitatively consistent with a shock-induced thermal non-
equilibrium between electrons and ions, and the second feature is consistent
with accretion shocks seen in previous studies. We split the cluster sample by
redshift and mass, and find both features exist in all cases. There are also
no significant differences in features along and across the cluster major
axis, whose orientation roughly points towards filamentary structure. As a
consistency test, we also analyze clusters from the Planck and Atacama
Cosmology Telescope Polarimeter surveys and find quantitatively similar
features in the pressure profiles. Finally, we compare the accretion shock
radius ($R_{\rm sh,\,acc}$) with existing measurements of the splashback
radius ($R_{\rm sp}$) for SPT-SZ and constrain the lower limit of the ratio,
$R_{\rm sh,\,acc}/R_{\rm sp}>2.16\pm 0.59$.

###### keywords:

galaxies: clusters: intracluster medium – large-scale structure of Universe

††pubyear: 2021††pagerange: Shocks in the Stacked Sunyaev-Zel’dovich Profiles
of Clusters II: Measurements from SPT-SZ + Planck Compton-$y$ Map–A

##  1 Introduction

Galaxy clusters are massive structures that contain multiple components, of
which dark matter, ionized gas, and galaxies are the dominant ones. The dark
matter is essentially collisionless and responsive only to gravity, while the
ionized gas responds to hydrodynamical and electromagnetic forces in addition
to gravitational ones. Galaxies contain stars — which are collisionless like
dark matter — in addition to dark matter and multi-phase gas, and thus respond
to all three forces mentioned above. The interactions within and between the
dark matter and ionized gas, the two components that make up $>99\%$ of the
mass in the cluster, determine the cluster’s internal structure and energetics
(see Kravtsov & Borgani, 2012, for a review).

Clusters are also dynamically young, having formed recently in cosmic history
($z<2$), and are actively accreting matter from their surroundings, with this
accretion happening preferentially along directions of the filamentary large-
scale structure. Consequently, each of dark matter, ionized gas, and galaxies
can be further deconstructed into two sub-components — one belonging to the
fully-collapsed, bound structure and another to the infalling component that
originates in the large-scale structure. Naturally, the study of the galaxy
cluster outskirts (for any component) is a key part of both astrophysical and
cosmological studies as these outskirts are the transition regime between the
two sub-components, and contain an abundance of dynamical information about
clusters and their interactions with their environment (Walker et al., 2019).

This work focuses on the gaseous component of the clusters, and in particular
on the pressure profiles where sharp, shock-like features can arise from
interactions between the gas of the two sub-components. We infer these
pressure profiles via the thermal Sunyaev-Zel’dovich (tSZ) signature of
clusters (Sunyaev & Zeldovich, 1972), which arises from the inverse compton
scattering of Cosmic Microwave Background (CMB) photons off energetic
electrons in the hot intracluster medium (see Carlstrom et al., 2002;
Mroczkowski et al., 2019, for reviews). While cluster thermodynamics have
traditionally been studied using X-ray observations, the tSZ has emerged as
the more ideal probe for the cluster outskirts as the signal amplitude depends
linearly with density, whereas for X-rays this dependence is quadratic.

The study of shocks is highly relevant to cluster-based studies of both
cosmology and astrophysics given that they are a critical mechanism during
structure formation for converting gravitational potential energy into thermal
energy. Shocks can induce significant deviations in cluster pressure profiles,
and can set up thermodynamic non-equilibrium conditions that invalidate common
assumptions made in estimating hydrostatic cluster masses; these masses are a
relevant quantity for doing cosmology with cluster counts (see Allen et al.,
2011, for a review), while a clear understanding of hydrostatic equilibrium in
clusters is also necessary for certain cluster-based tests of modified gravity
(Terukina et al., 2014; Wilcox et al., 2015; Sakstein et al., 2016; Haridasu
et al., 2021). Notably, the process of shock heating generates a thermal non-
equilibrium between the electrons and ions, which can alter the expected X-ray
and tSZ emissions 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).

Shocks can also be sources for accelerating cosmic ray electrons via Diffusive
Shock Acceleration (Drury, 1983; Blandford & Eichler, 1987). 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),
and cosmic rays in general alter the total pressure support of the system.
Near the cluster core, the pressure from cosmic rays has been observationally
constrained to be subdominant to the thermal pressure (e.g., Ackermann et al.,
2014) but simulations show it can be more prominent at the outskirts (Pfrommer
et al., 2007; Vazza et al., 2012).

The location of shock features also depends closely on the mass accretion rate
of the cluster and can potentially serve as an observational proxy for the
same (Lau et al., 2015; Shi, 2016; Zhang et al., 2020, 2021). The mass
accretion rate has strong theoretical connections to key dark matter halo
properties like concentration and formation time (Wechsler et al., 2002), and
has also been shown to have significant correlations with a broad range of
halo properties (e.g., Lau et al., 2021; Anbajagane et al., 2021). However, it
has remained difficult to infer observationally.

Accurate measurements of gas profiles at the cluster outskirts — particularly
near and beyond the one-to-two-halo transition regime — improve the modelling
needed in studies of the tSZ auto-correlation (e.g., Hill & Pajer, 2013;
Horowitz & Seljak, 2017; Tanimura et al., 2021) as well as tSZ cross-
correlations with galaxy and galaxy cluster positions (e.g., Hajian et al.,
2013; Vikram et al., 2017; Hill et al., 2018; Pandey et al., 2019, 2020), with
weak lensing shears (e.g., Ma et al., 2015; Hojjati et al., 2017; Osato et
al., 2018; Osato et al., 2020; Shirasaki et al., 2020; Gatti et al., 2021;
Pandey et al., 2021), or with X-ray luminosity (Shirasaki et al., 2020); such
studies can provide strong and complementary constraints on both cosmological
and astrophysical parameters.

Certain shock features form a boundary around the gaseous halo and delineate
the cold, pristine gas of the infalling regions from the hot, thermalized gas
of the bound structure. This boundary thus marks the radius within which
galaxies are first affected thermodynamically by the cluster gas, which
consequently impacts the galaxies’ evolution (e.g., Zinger et al., 2016a) via
processes like ram-pressure stripping (Boselli et al., 2021). The cosmic rays
generated by shocks can also potentially explain the still-unconfirmed origins
of radio relics in clusters (e.g., Vazza et al., 2012; Hong et al., 2014; Ha
et al., 2018) in addition to amplifying seed magnetic fields within clusters
(see Dolag et al., 2008; Donnert et al., 2018, for reviews). The magnetic
fields also have an inherent non-thermal pressure, and so can impact the total
pressure support of a cluster and thus, the hydrostatic cluster mass
estimates, just like cosmic rays.

Even more can be learned upon combining the thermodynamic gas structure with
the distribution of dark matter and galaxies. One such combination is to
compare the shock radii with the splashback radius (e.g., Diemer & Kravtsov,
2014a; Adhikari et al., 2014; More et al., 2015; Mansfield et al., 2017; Aung
et al., 2020; Xhakaj et al., 2020; O’Neil et al., 2021; Dacunha et al., 2021),
which is a physically motivated halo boundary defined by the apocenter in the
dark matter phase space of the halo. The existence of the splashback feature
has been observationally verified 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., 2020; Shin et al., 2021), and has
been shown to play a role in galaxy formation physics (Baxter et al., 2017;
Shin et al., 2019; Adhikari et al., 2020; Dacunha et al., 2021). 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). These two features are also sensitive to different types
of mass accretion — the shock radius evolves according to smooth accretion,
which does not include accretion of subtructure, whereas the splashback radius
depends on the total accretion rate (Zhang et al., 2021) — so combining the
two could potentially constrain the amount of mass accreted via merging
substructures. Figure 1 shows a diagram of the features discussed above in
relation to more commonly used cluster radius definitions.

Figure 1: Illustration of the different cluster radii relevant to this work,
denoted by colors and labelled by text. Radii are shown on linear axes, and
are drawn to relative scales determined for the median SPT cluster mass and
redshift — $R_{\rm 500c}=0.8\,{\rm Mpc}$, $R_{\rm 200m}=1.6\,{\rm Mpc}$,
$R_{\rm sp}=1.2R_{\rm 200m}=1.9\,{\rm Mpc}$, and $R_{\rm sh,\,acc}=2.3R_{\rm
200m}=3.7\,{\rm Mpc}$. The estimate for the ratio $R_{\rm sp}/R_{\rm 200m}$
was taken from Shin et al. (2019), while $R_{\rm sh,\,acc}/R_{\rm 200m}$ was
taken from Baxter et al. (2021). The blue lines show different filaments
connected to the galaxy cluster. The accretion shock is expected to be
weak/non-existent along filamentary directions. The splashback radius of the
dark matter can be smaller than $R_{\rm 200m}$ depending on the mass accretion
rate of the cluster.

Hydrodynamical simulations show that shocks can be generated at different
radial locations via different mechanisms, and to zeroth order there are two
governing processes: (i) the accretion of gas onto the cluster, i.e. the
interaction between the “bound” and infalling components, and; (ii) the major
and minor mergers with gas clumps, galaxies, and other clusters. The accretion
of pristine cold gas — which has a low sound speed and is primarily found in
low-density regions such as cosmic voids — onto the thermalized, bound gas
subcomponent results in a shock of a high mach number and discontinuities in
the profiles of many thermodynamic quantities such as temperature, entropy,
pressure, and density. This shock is oftentimes referred to as an accretion
shock (e.g., Lau et al., 2015; Aung et al., 2020; 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
in the hot, ionized gas leads to a series of bow shocks with weak mach
numbers, which are referred to as internal shocks (Ryu et al., 2003).
Furthermore, Zhang et al. (2019, 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. This is a common process
during structure formation, and so most shocks observed in the cluster
outskirts ($R\gtrsim R_{\rm 200m}$) are expected to be MA-shocks. While the
origin of the MA-shock is rooted in merger events, the radial evolution of the
feature — once it has been generated — depends only on smooth mass accretion
(Zhang et al., 2021). Finally, we stress that all of these processes detailed
above are complex, highly aspherical, and vary significantly from cluster to
cluster.

As was noted before, the current picture of shocks has been studied
predominantly using hydrodynamical simulations. Initial studies used non-
radiative simulations that modelled gas dynamics but did not include any non-
gravitational processes such as gas cooling (Quilis et al., 1998; Miniati et
al., 2000; Ryu et al., 2003; Skillman et al., 2008; Molnar et al., 2009; Hong
et al., 2014, 2015; Schaal & Springel, 2015). More recent studies have
included the effects of gas cooling and star formation (Vazza et al., 2009;
Planelles & Quilis, 2013; Lau et al., 2015; Nelson et al., 2016; Aung et al.,
2020), and as well as the effects of feedback from supernovae and active
galactic nuclei (Kang et al., 2007; Vazza et al., 2013, 2014; Schaal et al.,
2016; Baxter et al., 2021; Planelles et al., 2021). Some work has self-
consistently modelled the evolution of cosmic-rays alongside galaxy formation
(Pfrommer et al., 2007), while a handful have also used idealized simulations
to explore the propogation of shocks and their dependence on merger events
(Pfrommer et al., 2006; Ha et al., 2018; Zhang et al., 2019, 2020, 2021).
These works use a wide variety of hydrodynamical solvers and astrophysical
model prescriptions; see Vazza et al. (2011) and Power et al. (2020) for
comparisons of different implementations.

These works are accompanied by observational studies of shocks that have
focused predominantly on small samples — often containing just one object — of
local, low-redshift clusters (Akamatsu et al., 2011; Akahori & Yoshikawa,
2012; Akamatsu et al., 2016; Basu et al., 2016; Di Mascolo et al., 2019a; Di
Mascolo et al., 2019b; 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, have normally not pushed beyond $r\gtrsim R_{\rm 500c}$
(e.g., McDonald et al., 2014; Ghirardini et al., 2017; Romero et al., 2017,
2018; Ghirardini et al., 2018), though some do exist (Planck Collaboration et
al., 2013; Sayers et al., 2013, 2016; Amodeo et al., 2021; Schaan et al.,
2021).

With the advances made by many modern CMB experiments, tSZ maps have achieved
significantly improved angular resolution and lower noise levels, and this has
happened alongside the construction of large catalogs of clusters ($N\sim
10^{3}-10^{4}$) by multiple different surveys. Together, these improvements
enable population-level analyses of the shock features in the thermodynamic
gas profiles, and especially of any features in the cluster outskirts
($R\gtrsim R_{\rm 200m}$). We present here such an analysis of $N=516$ galaxy
clusters from the SPT-SZ survey, a 2500 deg2 survey conducted with the South
Pole Telescope (SPT). This is a companion work to Baxter et al. (2021,
henceforth Paper I), and to our knowledge is the first observational
population-level study of such features in the cluster outskirts.

The key goals of this work are to (i) extract the stacked tSZ profiles of
clusters and measure deviations from theoretical expectations, such as those
deviations induced by shocks and/or by other non-equilibrium processes, (ii)
study the dependence, or lack thereof, of such deviations on cluster mass and
redshift, and their variation along the cluster major vs. minor axes, and;
(iii) compare the shock radii with the splashback radius. Additionally, our
focus on the cluster outskirts naturally provides constraints for the one-to-
two halo transition regime. Throughout our analysis, we simultaneously analyze
The Three Hundred simulations (The300, Cui et al., 2018) — which was also
employed by the study in Paper I — both to test our pipeline and to compare
simulation predictions with observations. We also compare the SPT-SZ
measurements with public data from Planck and the Atacama Cosmology Telescope
Polarimeter (ACTPol) as a consistency test.

This work is organized as follows: We describe our observational and
simulation datasets in §2, including our procedure for generating mock
catalogs. In §3, we detail both our measurement procedure for obtaining the
stacked pressure profiles as well as the formalism of the theoretical model we
compare our results to. Our main results are presented in §4, while
comparisons to external data and the splashback radii are in §5. We discuss
our findings and conclude in §6.

##  2 Data

In this work, we use data from the SPT-SZ survey to set observational
constraints, and simulated clusters from The300 to test and validate our
analysis pipelines. We also compare our SPT-SZ results to those using data
from the Planck and ACTPol surveys. We describe each of these datasets below,
including a description on how we construct a mock catalog from the
simulations to match the SPT-SZ data.

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

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

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

All three of SPT-SZ, ACTPol, and Planck infer $M_{\rm 500c}$ from the
integrated tSZ emission around each cluster, while the simulated The300
catalogs provide $M_{\rm 200c}$ which is computed directly from the particle
data. Both $M_{\rm 500c}$ and $M_{\rm 200c}$ are defined by equation 1 but
with alternative density contrasts of $\rho_{\Delta}=500\rho_{c}(z)$ and
$\rho_{\Delta}=200\rho_{c}(z)$, respectively. Here, $\rho_{c}(z)$ is the
critical density of the Universe at a given epoch. In all cases, we convert
masses of alternative definitions into $M_{\rm 200m}$ using the concentration-
mass relation from Diemer & Joyce (2019) and the publicly available routine
from the COLOSSUS111https://bdiemer.bitbucket.io/colossus/ 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., 2020).

The tSZ amplitude is reported as the dimensionless $y$ parameter,

$y\equiv\frac{k_{B}\sigma_{T}}{m_{e}c^{2}}\int n_{e}T_{e}dl,$ (2)

where $k_{B}$ is the Boltzmann constant, $\sigma_{T}$ is the Thomson cross-
section, $m_{e}c^{2}$ is the rest energy of an electron, $n_{e}$ and $T_{e}$
are the electron number density and temperature, respectively, and $l$ is the
physical line-of-sight distance. Thus $y$ represents the electron pressure
integrated along the line-of-sight. For the rest of this work, we use the
terms tSZ and $y$ interchangeably; the former is the physical process of
interest, but the latter is the actual measurement provided in the maps.

The tSZ effect corresponds to CMB photons scattering off electrons with a
thermal (i.e. Maxwellian) energy/momentum distribution. Analogous effects,
called the relativistic SZ (rSZ) and non-thermal SZ (ntSZ), correspond to non-
Maxwellian energy distributions and may leak into the measured tSZ
(Mroczkowski et al., 2019). In the rSZ effect, the presence of high-
temperature electrons ($T_{\rm e}\gtrsim 5\,{\rm keV}$) requires relativistic
corrections to the map-making procedure. However, these corrections are at
most $5\%$ and are subdominant to the features discussed in this work (Erler
et al., 2018, see Figure 1). The ntSZ effect can be generated by a cosmic ray
electron population, but is a subdominant effect within $R_{\rm 200c}$ of the
cluster, where cosmic rays make up $\lesssim 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 assume the ntSZ continues to be a
subdominant source in the outskirts, but note that the features we discuss
here are qualitatively unaffected even if the ntSZ contaminates the tSZ at the
$10\%$ level.

###  2.1 The South Pole Telescope SZ (SPT-SZ) Survey

SPT-SZ is a $2500\rm\,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 $y$-map used in our analysis was presented in Bleem et al. (2021),
has an angular resolution of $1.25^{\prime}$, 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 (multipole $\ell\lesssim 1000$). The utilized Planck data consists of
the 100, 143, 217, and 353 GHz maps from the High Frequency Instrument (HFI).
The $y$-map is constructed via the commonly used Linear Combination (LC)
algorithm (see Delabrouille & Cardoso, 2009, for a review) which is applied to
the maps of different frequencies, and here 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. (2021) 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.

The associated galaxy cluster catalog is derived from the same data used to
construct the y-map and contains 516 clusters that were first identified in
Bleem et al. (2015), and with updated redshifts and mass estimates provided in
Bocquet et al. (2019). We use the latter, updated catalog for our work. Both
the map and the cluster catalog are publicly
available.222https://lambda.gsfc.nasa.gov/product/spt/spt_prod_table.cfm Our
masses come from the M500 column and signal-to-noise ratio (SNR) from the XI
column.

We also require accurate estimates of the noise in the $y-$map, which then
translates into noise in the pressure profiles estimates, when constructing
our mock catalog from the simulations. We estimate this for SPT-SZ using the
provided half1 and half2 maps. Each was constructed using half the SPT-SZ and
Planck data, and with the same LC procedure described above. Taking the
difference of these maps nulls any signal, and results in an accurate estimate
of the non-astrophysical, instrument-based noise contribution to the $y$-map,

$\texttt{Noise}=\frac{1}{2}(\texttt{half1}-\texttt{half2}).$ (3)

We use this map in conjunction with the simulations to assess the impact of
noise on the pressure profiles measured for SPT-SZ clusters. We stress that
this noise map lacks astrophysical contaminants such as point-sources and
dust, meaning the derived noise estimates exclude astrophysical contributions,
but such contaminants are also aggressively masked in our analysis. Note also
that this noise map is only used to generate the mock catalog from
simulations; in particular, the covariance matrix used in analyzing the
observational data is built from the full maps, including all unmasked
astrophysical sources.

Figure 2: The mass-redshift plane of the cluster samples from SPT-SZ, ACTPol,
and Planck used in this work. The top and right panels show the 1D
distributions for redshift and cluster mass, respectively. The SPT-SZ and
ACTPol samples have very similar redshift distributions, with a median of
$z\approx 0.55$, while Planck is at a lower redshift. Consequently, Planck
also has more high-mass clusters. The color tones show $\log_{10}\rm 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
can be found in Table 1.

###  2.2 Atacama Cosmology Telescope Polarimeter (ACTPol)

The ACTPol survey observes at 98 and 150 GHz, and the 2014-2015 observations
cover two regions labelled as BN and D56, with an area of $1663\rm\,deg^{2}$
and $456\rm\,deg^{2}$, respectively. The combined area spans $\approx
2100\rm\,deg^{2}$ of the sky. The $y$-map was presented in ACT Data Release
4333https://lambda.gsfc.nasa.gov/product/act/act_dr4_derived_maps_info.cfm
(DR4, Madhavacheril et al., 2020), has a resolution of $1.6^{\prime}$, and
makes use of data from both ACTPol and the Planck 2015 data release; as was
the case with SPT-SZ, the former data inform small-scales and the latter, the
large-scales ($\ell\lesssim 1000$). Note that the Planck data here consist of
eight frequency channels from 30 to 545 GHz, whereas SPT-SZ used four of these
channels. Three of the four additional channels used in ACTPol contribute
minimally to the final map (Madhavacheril et al., 2020, see their Figure 5),
while the remaining one — the 70 GHz channel — has notable contributions below
$\ell\lesssim 100$. The map is made separately for each of the two regions
using an Internal Linear Combination (ILC) algorithm.

In our analysis, the map is further masked to remove point sources and dusty
regions. The ACTPol masks are continuous, not binary, and we continue with our
aggressive masking by only selecting pixels for which the mask value is 1,
meaning the pixel is uncontaminated. Note that the ACTPol maps, unlike all
other maps in this analysis, do not use the Healpy pixelation scheme and
instead work with a new scheme called
Pixell.444https://pixell.readthedocs.io/en/latest/ Also note that we perform a
combined analysis of the clusters from each region and do not treat them as
separate datasets.

While the public ACTPol $y$-map is part of DR4, the public cluster catalog is
from
DR5555https://lambda.gsfc.nasa.gov/product/act/actpol_dr5_szcluster_catalog_get.cfm
(Hilton et al., 2021), and many DR5 clusters are outside the DR4 y-map
footprint.666We define a cluster to be “outside the footprint” if a circle of
radius $R=20R_{\rm 200m}$ around the cluster does not intersect with the
footprint. The radius choice reflects the largest scales explored in our
analysis, as described in Section 3.1. So, though the full DR5 catalog
contains $N_{\rm clusters}\approx 4200$, we are only able to use $N_{\rm
clusters}\approx 1400$. As a result, the mass-redshift distribution shown in
Figure 2 differs from a similar one shown in Hilton et al. (2021, see their
Figure 18). Coincidentally, the redshift distribution of the reduced ACTPol
cluster sample used in this work is very similar to that of the SPT-SZ sample.
For our cluster mass variable, we use the M500cCal column described in Hilton
et al. (2021, see Table 1), which contains a richness-based, weak-lensing mass
calibration factor that brings the ACTPol cluster masses in better agreement
with those from SPT-SZ.

###  2.3 Planck

The Planck satellite mission is a survey of the full-sky and thus, overlaps
with both SPT-SZ and ACTPol. The $y$-map, which uses data collected up to
2015, is presented in Planck Collaboration et al. (2016a) and has an angular
resolution of $10^{\prime}$. We use the map constructed using Milca (Hurier et
al., 2013) — a Modified ILC Algorithm — on the HFI individual frequency maps
from the range 100 to 857 GHz. We apply the publicly-available masks to the
maps to remove all point sources and an additional 45% of the sky contaminated
by dust.

The cluster catalog was presented in Planck Collaboration et al. (2016b), and
has a mass/redshift distribution that significantly differs from both SPT-SZ
and ACTPol. The derived $M_{\rm 500c}$ masses (provided under column MSZ) are
also known to be biased low when compared to SPT (Planck Collaboration et al.,
2016b; Battaglia et al., 2016; Bleem et al., 2020). By comparing a subset of
$94$ clusters that overlap between Planck and SPT-SZ, we derive a ratio,
$\langle M_{\rm 500c}^{\rm SPT}/M_{\rm 500c}^{\rm Planck}\rangle=1.28\pm
0.04$. This ratio is then used to recalibrate the Planck masses. We stress
that this is only an approximate recalibration, given the mass bias is known
to depend on cluster redshift as well (see Section 5.3 in Bleem et al., 2020).
We further discuss the validity and necessity of this procedure in our results
(Section 5.1). We do not explore a more rigorous re-calibration in this work
given Planck only serves as a comparison point and is therefore not the focus
of our analysis.

###  2.4 The Three Hundred Project (The300)

The300 (Cui et al., 2018) is a set of zoom-in hydrodynamical simulations of
324 massive halos ($M_{\rm 500c}\gtrsim 10^{14.8}\,{\rm M}_{\odot}$ at $z=0$).
The simulations are run by first identifying the 324 most massive objects in
MultiDark Planck 2 (MDPL2, Klypin et al., 2016) — a $1.5\,{\rm Gpc}$ N-body,
dark matter only simulation with purely gravitational evolution — and then re-
simulating those halos with the hydrodynamics code Gadget-X (Beck et al.,
2016), a modified version of Gadget3 (last described in Springel, 2005), which
contains a full prescription of galaxy formation physics such as metal
cooling, star formation, kinetic and thermal feedback from supermassive black
holes etc. (Rasia et al., 2015). The re-simulated regions are spheres of
comoving radii $R=22\,{\rm Mpc}$ and contain multiple, smaller halos in
addition to the most massive halo that was first identified in MDPL2. Halos
and subhalos are identified with Amiga’s Halo Finder (Ahf, Knollmann & Knebe,
2009), which uses an adaptive mesh refinement grid to represent the density
field/contours and also has a binding energy criterion to remove unbound
material from halos/substructure. The mass computed in AHF follows the $M_{\rm
200c}$ definition, and so we perform the same procedure as in §2 to convert
$M_{\rm 200c}\rightarrow M_{\rm 200m}$ for our analysis.

The simulated tSZ/$y$-maps were constructed using the PyMSZ
package777https://github.com/weiguangcui/pymsz (Cui et al., 2018). At a given
redshift, there are 324 maps, each corresponding to an individual re-
simulation region, and then we have ten different maps for each region
corresponding to the redshifts, $z\in[0.1,0.2,\ldots,1.0]$; this redshift
range covers $>95\%$ that of the SPT-SZ sample. The maps are also constructed
to have the same $1.25^{\prime}$ angular resolution as SPT-SZ. When comparing
observations to simulations, we subsample The300 catalogs to replicate the
mass and redshift distributions of SPT-SZ. This is done in a hierarchial way —
for each SPT-SZ cluster, we first identify the closest discrete redshift
$z\in[0.1,0.2,\ldots,1.0]$. Next, we consider all 324 simulated clusters at
that redshift and choose the one with an $M_{\rm 200m}$ mass closest to that
of the SPT-SZ cluster. Thus, redshift matching takes precedence over mass
matching. During such subsampling and comparisons, we limit the SPT-SZ catalog
to $z<1.1$ in order to better match the available $y$-maps from The300, which
stop at $z=1$. This cut excludes only 24 out of the 516 SPT-SZ clusters.

Note that the SPT-SZ sample has 516 independent clusters whereas The300
follows the same set of 324 clusters across different redshifts. So even
though we have The300 catalogs at multiple redshifts, we only have 324
independent clusters. Thus, during the subsampling step described above, we
are forced to select the same clusters from multiple different redshifts, and
so our mass- and redshift-matched sample from The300 will not be a completely
independent set of clusters. This induced correlation will lead to an
underestimate of the cluster-to-cluster variance in the pressure profiles;
this is only a minor issue as most of the scales we study here are dominated
by noise variance instead.

To ensure we do not study radial scales that extend beyond the $R=22\,{\rm
Mpc}$ spherical volume, we limit our analysis of the simulations to $R<6R_{\rm
200m}{}$, where only 2 out of the 492 clusters from the mass- and redshift-
matched The300 sample extend slightly beyond the sphere.

##  3 Measurement and Modeling

We first describe our procedure for measuring the stacked tSZ profile and
other associated quantities in the SPT data, and then the theoretical halo
model we compare the measurements with, including how we quantify the
significance of any features in the data.

###  3.1 Measurement Procedure

We detail below our method for estimating the (i) stacked profiles, (ii)
logarithmic derivatives, (iii) bin-to-bin covariance matrix, and; (iv) the
feature locations.

Estimating stacked profiles: For each cluster in our catalog, we compute its
halo-$y$ correlation, i.e. the average $y$ profile, in 50 logarithmically
spaced radial bins that span the range $r\in[0.1,20]R_{\rm 200m}$. This is
done by measuring the average $\langle y\rangle$ in each radial bin and
subtracting the mean background value from it, where the latter is estimated
by populating the map with random points and computing profiles around them.
This entire process is performed using the fast tree-code implementation in
the software package TreeCorr (Jarvis et al., 2004). The logarithmic spacing,
as opposed to linear spacing, is an apt choice here as the signal (the tSZ
emission) drops approximately as a power law with radius, implying a lower SNR
per pixel when going from cluster core to the cluster outskirts. Thus,
including more pixels per bin towards the outskirts will help increase our SNR
per radial bin.

The profiles of the individual clusters are then stacked, with each profile
being weighted by the corresponding cluster’s detection SNR. Our final results
do not change from using alternative weighting choices (see Appendix A). Note
that for a given cluster, any radial bin that did not have a single pixel in
it — most commonly the case in the cores of high redshift clusters due to the
limited angular resolution — is ignored during the stacking. In tandem to the
stacking procedure, we perform a leave-one-out jackknife resampling,

$\langle y\rangle_{i}(r)=\frac{1}{N-1}\sum_{j\neq i}^{N}\,\,y_{j}(r),$ (4)

where $N$ is the number of measured cluster profiles, $\langle y\rangle_{i}$
is the mean profile of the sample with cluster $i$ removed, and $y_{j}$ is the
individual profile measurement from cluster $j$. The variance on the mean
profile is then given by,

$\sigma^{2}(\langle y\rangle)=\frac{N-1}{N}\sum_{j=1}^{N}\,\,\bigg{(}\langle
y\rangle_{j}-\bar{\langle y\rangle}\bigg{)}^{2},$ (5)

where $\bar{\langle y\rangle}$ is the mean of the distribution of jacknife
estimates computed in equation (4). Note that equation (5) has an additional
factor of $N-1$ compared to the traditional definition of the variance.

Estimating logarithmic derivatives: Shock features are identified by points of
steepest descent in the pressure profiles (e.g., Paper I, Aung et al., 2020),
and this corresponds to measuring minima in the logarithmic derivative.
However, derivatives are particularly susceptible to spurious noise-induced
features. To alleviate this, the stacked profiles from the jackknife sample
are all smoothed by a Gaussian of width $\sigma_{\ln r}=0.16$, which is $1.5$
times the logarithmic bin width, ${\rm d\ln}r\approx 0.11$. All the profiles
we show below have been smoothed by this scale. The smoothing step will induce
edge effects at the lower/upper radial limits of $0.1R_{\rm 200m}$ and
$20R_{\rm 200m}$ given there is no measured profile beyond those bounds. For
this reason, we only quote results for the range $0.3<R/R_{\rm 200m}<10$. The
impact of different smoothing choices, including an essentially un-smoothed
case, is discussed in Appendix A, and our results remain the same even when
using alternative choices.

We then compute the log-derivative of the smoothed mean profile, ${\rm
d}\ln\langle y\rangle/{\rm d}\ln(R/R_{\rm 200m})$, using a five-point method,

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

where $f$ is an arbitrary function of $x$, and $h={\rm d}\ln r$ is the spacing
between the sampling points. The numerical error in this derivative estimator
goes as $\mathcal{O}(h^{4})$. The uncertainty on the log-derivative is
estimated by computing equation (6) for every stacked profile in the jackknife
sample and taking the standard deviation of the resulting distribution. An
extra multiplicative factor of $\sqrt{N-1}$ is then needed to convert the
measured uncertainty to the true uncertainty, and this is entirely analogous
to the extra $N-1$ factor in equation (5).

Covariance of the log-derivative: We also need the bin-to-bin covariance
matrix, $\mathcal{C}$, of the log-derivatives when computing a detection
significance, as is discussed further below in equation (21). This covariance
is estimated using the stacked profiles of the jackknife sample,

$C_{i,j}=\frac{N-1}{N}\sum_{k=1}^{N}\,\,\bigg{(}f^{\prime}_{k,i}-\langle
f^{\prime}\rangle_{i}\bigg{)}\bigg{(}f^{\prime}_{k,j}-\langle
f^{\prime}\rangle_{j}\bigg{)},$ (7)

where $i$ and $j$ index over the different radial bins, $f^{\prime}_{k,i}$ is
the log-derivative of the $k^{\rm th}$ stacked profile in the $i^{\rm th}$
radial bin, and $\langle f^{\prime}\rangle_{i}$ is the mean log-derivative in
the $i^{\rm th}$ radial bin. For $i=j$, equation (7) reduces to equation (5),
but for the log-derivatives instead of the profiles.

Quantifying feature location: To determine the location of a given feature —
particularly, of local minima in the log-derivative — we fit cubic splines to
the log-derivative of each stacked profile in the jackknife sample and locate
the feature of interest in each. The mean and standard deviation of the
resulting distribution provides an estimate for the location of the feature
and the associated uncertainty. The $\sqrt{N-1}$ factor is needed once again
to go from the measured uncertainty to the true uncertainty.

Note that our estimates of the locations, $R/R_{\rm 200m}$, do not include the
uncertainties in the inferred $R_{\rm 200m}$ of each cluster. For the cluster
samples in this work, these uncertainties in $R_{\rm 200m}$ are $<7\%$, which
are tolerable given they increase the total uncertainty in the estimated
feature location by $<2\%$.

###  3.2 Modeling and Detection Quantification

In addition to comparing our observational results to those from simulations,
we also compare the former with commonly used theoretical models, and this
will be key in quoting a detection significance for any interesting features.
The model we employ here for the halo-$y$ correlation 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 that accounts for
contributions from nearby halos as described in Vikram et al. (2017) and later
in Pandey et al. (2019). Our two-halo term modelling is based on the linear
matter power spectrum and linear halo bias, and assumes higher-order
corrections are not required. We validate this assumption below by showing
that the linear model accurately describes the profiles of the simulated
cluster populations. Our entire modelling procedure is done using the Core
Cosmology Library (CCL) open-source python
package888https://github.com/LSSTDESC/CCL (Chisari et al., 2019).

We start by computing the 3D, halo-pressure cross-correlation function as the
sum of the one-halo and two-halo components,

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

where $\xi$ are the correlation functions, $r$ is comoving distance, and $M$
is the halo mass. We will henceforth denote the combined one-halo and two-halo
term as the “total halo model”. The one-halo term is modelled using the
pressure profile from Battaglia et al. (2012),

$P(x)=P_{\rm
500c}P_{0}\bigg{(}\frac{x}{x_{c}}\bigg{)}^{\gamma}\bigg{[}1+\bigg{(}\frac{x}{x_{c}}\bigg{)}^{\alpha}\bigg{]}^{-\beta},$
(9)

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

$P_{\rm 500c}=500\rho_{c}(z)\frac{\Omega_{b}}{\Omega_{m}}\frac{GM_{\rm
500c}}{2R_{\rm 500c}}.$ (10)

Note that while the theory includes the self-similar pressure, $P_{\rm 500c}$,
the actual pressure profile normalization is still given by the combination
$P_{\rm 500c}P_{0}$ and accounts for deviations from self-similar evolution
via the calibrated parameter $P_{0}$. The fit parameters for equation (9) are
obtained from the “Shock Heating” calibration model of Battaglia et al. (2012,
see Table 1), and these parameters have a known, calibrated scaling with both
cluster redshift, $z$, and cluster mass, $M_{\rm 500c}$.999When computing
these theoretical predictions for The300 clusters, we first convert the
provided $M_{\rm 200c}$ masses to $M_{\rm 500c}$ using the same process
described in Section 2 and then use the latter as the input to get the
parameters of the theory model.

The tSZ emission we study here is connected to the electron pressure, $P_{e}$,
whereas the Battaglia profiles are calibrated to the total gas pressure, $P$,
so we convert between the two as,

$P_{e}(r,M,z)=\frac{4-2Y}{8-5Y}P(r,M,z),$ (11)

with $Y=0.24$ being the primordial helium mass fraction. This expression now
serves as our one-halo term,

$\xi_{h,p}^{\rm one-halo}(r,M,z)=P_{e}(r,M,z).$ (12)

The two-halo term is more conveniently computed in Fourier space, so we
perform all our computations in the same and then inverse Fourier transform in
the end to get the correlation function. The two-halo term of the halo-
pressure cross-power spectrum, $P^{\rm two-halo}_{h,p}$, is written as,

$\begin{split}P^{\rm two-halo}_{h,p}(k,M,z)&=\bigg{[}b(M,z)\,P_{\rm
lin}(k,z)\,\,\times\\\
&\qquad\int_{0}^{\infty}\kern-10.00002ptdM^{\prime}\,\frac{dn}{dM^{\prime}}\,b(M^{\prime},z)\,u_{p}(k,M^{\prime},z)\,\bigg{]},\end{split}$
(13)

where $M$ is the mass of the halo we are computing the halo-pressure
correlation for, $M^{\prime}$ is the mass of a neighbouring halo contributing
to the two-halo term, $P_{\rm lin}(k,z)$ is the linear matter density power
spectrum at redshift $z$, $dn/dM^{\prime}$ is the mass function of
neighbouring halos, and $b(M,z)$ and $b(M^{\prime},z)$ 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 $u_{p}(k,M^{\prime},z)$ is the Fourier
transform of the pressure profile about the neighboring halo which, under the
assumption of spherical symmetry, is computed as,

$u_{p}(k,M^{\prime},z)=\int_{0}^{\infty}dr4\pi
r^{2}\frac{\sin(kr)}{kr}P_{e}(r,M^{\prime},z),$ (14)

where $P_{e}$ is the electron pressure profile. The halo-pressure two-point
cross-correlation is then simply the inverse Fourier transform of the cross-
power spectrum,

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

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

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

$\xi_{h,y}(r,M,z)=\frac{\sigma_{T}}{m_{e}c^{2}}\int_{-\infty}^{\infty}\frac{d\chi}{1+z}\xi_{h,p}\bigg{(}\sqrt{\chi^{2}+r^{2}},M,z\bigg{)},$
(16)

where $\sigma_{T}$ is the Thomson scattering cross-section, $m_{e}c^{2}$ is
the rest mass energy of the electron, and $\chi$ is the comoving coordinate
along the line-of-sight.

The tSZ maps we use here have finite angular resolution, which suppresses
power on small scales, and so we incorporate this resolution limit into our
theory predictions. To do so, we first calculate the angular cross-power
spectrum, using the flat sky approximation, as,

$C_{\ell}=\int\,d\theta\,2\pi\theta\,J_{0}(\ell\theta)\,\xi_{h,y}(\theta,M,z),$
(17)

where $J_{0}$ is the zeroth-order Bessel function. We then multiply $C_{\ell}$
by the Fourier-space smoothing function for SPT-SZ and then perform an
inverse-harmonic transform,

$\xi^{\rm
smooth}_{h,y}(\theta,M)=\int\frac{d\ell\ell}{2\pi}J_{0}(\ell\theta)C_{\ell}B_{\ell},$
(18)

with the smoothing function $B_{\ell}$ given as

$B_{\ell}=\exp\bigg{[}-\frac{1}{2}\ell(\ell+1)\sigma_{\rm FWHM}^{2}\bigg{]},$
(19)

where $\sigma_{\rm FWHM}=\theta_{\rm FWHM}/\sqrt{8\ln 2}$, with $\theta_{\rm
FWHM}=1.25^{\prime}$ being the full-width half-max of the Gaussian filter used
to smooth the SPT-SZ maps. The quantity $\theta_{\rm FWHM}$ will take
different values when computing theoretical predictions for the ACTPol and
Planck surveys, as their smoothing scales differ from that of SPT-SZ (see
Section 2).

To obtain our final theory curve for a given cluster sample, we compute the
smoothed total halo model, $\xi^{\rm smooth}_{h,y}$, 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.
Finally, we quantify the detection significance, which is the significance of
a deviation in the measured log-derivatives away from the theoretical model,

$\epsilon\equiv\frac{1}{\sigma}\bigg{(}\frac{{\rm d\ln}y^{\rm obs}}{{\rm
d\ln}x}-\frac{{\rm d\ln}y^{\rm th}}{{\rm d\ln}x}\bigg{)},$ (20)

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

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

$\chi^{2}=\bigg{(}\frac{{\rm d\ln}y^{\rm obs}}{{\rm d\ln}x}-\frac{{\rm
d\ln}y^{\rm th}}{{\rm d\ln}x}\bigg{)}^{T}\mathcal{C}^{-1}\bigg{(}\frac{{\rm
d\ln}y^{\rm obs}}{{\rm d\ln}x}-\frac{{\rm d\ln}y^{\rm th}}{{\rm
d\ln}x}\bigg{)},$ (21)

where $\mathcal{C}$ is the covariance matrix of the log-derivative as defined
in equation (7). We do not use all 50 radial bins for this calculation and
instead limit ourselves to 8 bins surrounding the feature of interest. Given
our choice of logarithmically space bins with $d\ln r=0.11$, the total width
of 8 bins centered on some location $r_{\star}$ is approximately $\Delta
r\approx r_{\star}$. Our results do not vary much when using anywhere from 6
to 10 bins in the calculation instead. Once $\chi^{2}$ is computed, we then
quote the total signal-to-noise, $\chi$, of a feature. We note again that our
focus on the log-derivatives as our primary detection measure for shocks is
informed by recent simulation studies (e.g., Paper I, Aung et al., 2020).

##  4 Shocks in SPT-SZ

We first present our fiducial results for the SPT-SZ $y-$profiles in Section
4.1, then study the variation in the observed features with (i) cluster mass
and redshift in Section 4.2, and; (ii) cluster major vs. minor axis in Section
4.3.

All bands show $68\%$ uncertainties estimated via the jackknife distribution
of stacked profiles. As for the detection significance, we show $\epsilon$ in
the figures but quote $\chi$ in our discussions in the text, and these
quantities are defined in equations (20) and (21), respectively. The latter is
our preferred metric as it is the significance of the entire feature across
multiple radial bins, while the former is the single-bin significance.

All constraints on feature locations and their corresponding detection
significance are provided in Table 1. The measured location of any feature is
in principle shifted further out from its “true” location due to the effects
of smoothing in the maps. However, we have estimated the magnitude of this
shift by comparing theoretical models with and without smoothing, and find
that for both SPT-SZ and ACTPol this shift is negligible in comparison to the
uncertainties. The one exception is Planck, for which the shift is
significant, and we discuss this further in Section 5.

###  4.1 Measurements of Pressure Deficit and Accretion Shock

Figure 3: The mean $y$-profile for all $516$ SPT-SZ clusters (top) and the
corrsponding log-derivative (middle). The gray dotted (dashed-dotted) line is
the one-halo (two-halo) theoretical model described in Section 3.2, while the
black dashed line representing the total halo model is a simple sum of the two
components. The bottom panel shows the detection significance, described in
equation (20). The profile shows two features — the first at $R_{\rm
sh,\,de}/R_{\rm 200m}=1.08\pm 0.09$ with $\chi=3.1\sigma$ significance, and
the second at $R_{\rm sh,\,acc}/R_{\rm 200m}=4.58\pm 1.24$ detected at
$\chi=2\sigma$ significance. The possible physical origins of each feature are
discussed in Section 4.1.

In general, the stacked $y$-profile of the SPT-SZ clusters follows the trends
from the theoretical expectations across both the one-halo and two-halo regime
(top panel, Figure 3). There are, however, two clear differences of interest,
both found prominently in the log-derivatives. We will refer to these as the
pressure deficit and accretion shock, and denote their locations as $R_{\rm
sh,\,de}$ and $R_{\rm sh,\,acc}$, respectively. Comparisons of these SPT-SZ
features with the ACTPol and Planck datasets are performed in Section 5.1.

####  4.1.1 Pressure deficit at cluster virial boundary

First, the profile shows a shock-like steepening in the form of a pressure
deficit at $R_{\rm sh,\,de}/R_{\rm 200m}=1.08\pm 0.09$ which is not present in
the theory, and the corresponding feature in the log-derivatives — a local
minimum — is measured at $3.1\sigma$ significance.

This pressure deficit can be a sign of a thermal non-equilibrium between
electrons and ions 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 — which are the primary
mechanism for converting kinetic energy to thermal energy during structure
formation — preferentially heat the ions over the electrons given the former
are more massive and constitute most of the kinetic energy. This leaves the
electrons with a lower temperature than the ions. Normally, Coulomb
interactions between the two species re-establish thermal equilibrium, but at
the cluster outskirts the particle density is low and the thermal equilibrium
time-scale exceeds the Hubble time, as is demonstrated in the above works.
Thus, the cluster outskirts will remain in significant non-equilibrium and the
temperature difference will consequently impact the tSZ emission, which is
sensitive to only the electron temperature. Rudd & Nagai (2009, see Figure 2)
showcase this effect using cluster tSZ profiles from simulations specialized
to model the electrion-ion temperature differences, while Avestruz et al.
(2015, see Figure 1) do the same but using the 3D cluster temperature profiles
from such specialized simulations. In general, this phenomenon/feature is not
resolved in most cosmological hydrodynamical simulations, including all other
simulations discussed in this work, given they a-priori assume local thermal
equilibrium between electrons and ions.

X-ray observations of local-volume galaxy clusters have shown some indication
of thermal non-equilibrium (e.g., Akamatsu et al., 2011; Akahori & Yoshikawa,
2012; Akamatsu et al., 2016), though the findings are not yet conclusive and
require more data. On the tSZ side, Planck Collaboration et al. (2013, see
Figure 6) studied the pressure profiles of 62 clusters from the early Planck
SZ data and find that their measured profiles are statistically consistent —
given the large uncertainties from the modest sample size — with many
simulations that all assumed local thermal equilibrium.

Other recent works have found deviations in the thermodynamic cluster profiles
— though none explore thermal non-equilibrium as a potential cause. Hurier et
al. (2019) found a sharp decrease in pressure at $R\sim 2R_{\rm 500c}$ for a
single cluster in the Planck data. Pratt et al. (2021) also used Planck data
and found an excess in pressure at $R\sim 2R_{\rm 500c}$ for a set of ten,
low-redshift galaxy groups. Finally, Zhu et al. (2021) found an excess in the
temperature and density profiles of the Perseus cluster at $R\approx R_{\rm
200c}$ using Suzaku X-ray data. In all three works, the deviations are found
roughly around $R\sim R_{\rm 200m}$, but the exact nature of the deviation —
specifically, whether it is an excess or deficit compared to theory — varies.

Thermal non-equilibrium is also not the only viable explanation for the
observed pressure deficit. The total halo model theory we use was calibrated
using a suite of hydrodynamical simulations (see Battaglia et al., 2012, for
more details); in principle, any physical process that was not modelled in the
simulations could cause a deviation of the observed profiles from the theory.
Two notable processes missing in all the simulations we discuss here are
magnetic fields and cosmic rays. Both processes can exert a non-thermal
pressure, can interact with the ions and electrons, and are particularly
relevant in regions containing shocks (see Dolag et al., 2008, for a review).
Turbulent gas motions also constitute a significant non-thermal pressure
component (e.g., Nelson et al., 2014; Shi & Komatsu, 2014; Shi et al., 2015),
with increasing importance towards cluster outskirts. However, they are
unlikely to be the cause of the deviations seen here as such motions, and
their resulting non-thermal pressure, are properly resolved in all
hydrodynamical simulations.

The top panel of Figure 3 shows that the theoretical one-halo term (dotted
line) accurately describes the SPT-SZ pressure profiles, including the deficit
feature. However, upon adding the two-halo term, we find the theory no longer
captures the deficit. While one could take this to mean the one-to-two halo
transition in the theory is incorrect, this is unlikely to be the case as the
profiles in the simulations — which accurately resolve this transition — can
be modelled precisely using the same theory (bottom row, Figure 4). Thus, we
do not suspect that this feature arises from a one-to-two halo modelling
issue.

####  4.1.2 Accretion shocks in cluster outskirts

Next, moving further into the cluster outskirts, the profiles have a nearly
constant $\langle y\rangle$ between $R_{\rm 200m}<r<3R_{\rm 200m}$ before
decreasing sharply. This drop is characterized by a local minimum in the log-
derivative at $R_{\rm sh,\,acc}/R_{\rm 200m}=4.58\pm 1.24$, and this is
measured at $2\sigma$ significance. The radial range of the 8 nearest bins
used in the significance computation spans $3<R/R_{\rm 200m}<6.4$. While the
detection significance is not high, the qualitative behavior is consistent
with predictions for an accretion shock from Paper I (see their Figure 5). The
location of the minimum is larger, by $\approx 2\sigma$, compared to the
findings from previous work on simulations, which locate it at $\approx
2R_{\rm 200m}$ (e.g., Paper I; Molnar et al., 2009; Lau et al., 2015; Aung et
al., 2020). However, the “start” of the feature in the SPT-SZ data — indicated
by where the log-derivative starts decreasing the second time — is at
$R/R_{\rm 200m}=2.44\pm 0.25$, whereas the same “start” in the simulation
predictions is roughly around $R/R_{\rm 200m}\approx 1.9$ (see Paper I; Aung
et al., 2020). This could imply that the observed feature starts closer to the
expected location but is much broader than in past works so the location of
the minimum in SPT-SZ is further out. For example, the evolution of the shock
location over time is known to vary significantly with differences in the mass
accretion rate of the cluster (e.g., Paper I, Aung et al., 2020; Zhang et al.,
2021). The redshift range of the cluster sample we study here is also much
wider than that of the past works.

####  4.1.3 Comparisons with simulations

Figure 4: The stacked $\langle y\rangle$ profile (left) and log-derivative
(right) for The300 simulation across four redshifts (different colors). Top
panel show the results from the fiducial maps with $1.25^{\prime}$ Gaussian
smoothing. Middle panels show results from the smoothed plus noise-added maps.
The bottom panels use the same maps as the middle panels, but the $M_{\rm
200m}-z$ distribution of the The300 cluster sample was matched to the SPT-SZ
data. We limit the SPT-SZ cluster catalog to $z<1.1$ to better match the
simulation sample (which only goes up to $z<1$), and this leave us with 492
clusters. The bands for the simulation results show our estimates for the
noise-induced measurement uncertainty, and do not represent a theoretical
uncertainty. The mass- and redshift-matched sample from The300 agrees well
with theory, including in the one-to-two halo transition regime. The300
profiles do not show any features corresponding to a pressure deficit or
accretion shock similar to those seen in the SPT-SZ data. The inset in the
bottom left panel shows the profiles for $R/R_{\rm 200m}>3$, but on linear
y-axes, where the absolute uncertainties in the mean profile can be properly
compared. See Section 4.1.3 for details. The theory curves (dashed lines) in
the bottom row are computed as described in Section 3.2.

We then perform an explicit comparison between SPT-SZ and The300 in Figure 4.
We first show the mean profiles (left) and log-derivatives (right) and
progressively modify the simulated maps with Gaussian smoothing (top row) and
SPT-SZ noise (middle row), before comparing results of the mass- and redshift-
matched sample from The300 with those of the SPT-SZ sample (bottom row). As
noted before, The300 closely follows the theoretical model and this is not
surprising given the pressures profiles of Battaglia et al. (2012) were
themselves calibrated using hydrodynamical simulations.

While we are unable to find an accretion shock feature in the mean profile of
The300, previous works have shown that different cosmological simulations,
including The300, do resolve accretion shocks (Paper I; Molnar et al., 2009;
Lau et al., 2015; Aung et al., 2020). In particular, Paper I used the same
The300 simulations as we do here and showed that the stacked $y$-profile has a
strong accretion shock feature at $R/R_{\rm 200m}\approx 2.3$ at $z=0.2$,
whereas we find none. This is because their study specifically selects relaxed
clusters — using the fraction of total cluster mass contained in substructure
as the relaxation criteria — whereas we perform no such selection as this
better represents the SPT-SZ dataset. The inclusion of unrelaxed clusters in
the analysis of The300 significantly weakens the presence of the accretion
shock feature. This could be because relaxation state affects the self-similar
scaling of the shock feature with $R_{\rm 200m}$, and deviations can cause the
feature to be washed out during stacking. Consequently, the fact that we still
see a feature in SPT-SZ indicates that accretion shocks in the data may be
stronger than those realized in The300.

Figure 4 also validates our uncertainty quantification procedure by comparing
the uncertainties in the SPT-SZ mean $y-$profile with those of the noise-added
simulations, where the noise used in the latter comes from the SPT-SZ half-
maps as described in Section 2.1. Notably, the simulations have a much larger
fractional (or log) uncertainty than the data, whereas the absolute
uncertainties of the two agree much better as shown in the inset in the bottom
left panel. This is because the fractional uncertainty is amplified by The300
mean profile being lower than the SPT-SZ mean profile. Note that the cluster-
to-cluster variance in the simulations is an underestimate of the true
variance as the mass- and redshift-matched sample is not completely
independent and contains some correlation. This is due to the same simulated
clusters being selected from more than one redshift (see Section 2.4).

###  4.2 Trends with Mass and Redshift

Figure 5: Dependence of the pressure profiles on cluster mass (left) and
redshift (right). The panels follow the same style as Figure 3. The subsamples
are not defined by a simple split in each variable; see Section 4.2 for
details. There continues to be a deviation in the log-derivatives at $R\approx
R_{\rm 200m}$, and the significance is higher/lower than the fiducial case
depending on subsample (see Table 1). All subsamples show a plateau phase
followed by a sharp drop as is characteristic of an accretion shock feature.
We do not observe any statistically significant mass/redshift evolution of
either feature.

We next focus on the variation of the two features — the pressure deficit and
the accretion shock — with cluster mass and redshift. However, the simple
approach of splitting the sample according to a median mass or redshift will
entangle the effects of mass- and redshift-dependence as the cluster mass
evolves significantly over redshift. To properly disentangle these effects, we
take a slightly more sophisticated approach which we detail below.

First, to study the temporal evolution, we label each cluster with its peak
height

$\nu_{\rm 200m}=\frac{\delta_{c}}{\sigma(M_{\rm 200m},z)},$ (22)

as it is a less redshift-dependent mass definition compared to $M_{\rm 200m}$.
Here $\sigma(M_{\rm 200m},z)$ is the root-mean square fluctuation of the
linear matter density field, at redshift $z$, smoothed on the scale $R_{\rm
200m}$, and $\delta_{c}=1.686$ is the critical overdensity for collapse. The
cluster sample is then rank-ordered according to the peak height labels and
then binned into groups of 20 clusters. In each bin, the sample is split into
either a high or low redshift group using the median cluster redshift in that
bin. The high (low) redshift groups from all peak height bins are combined to
form the final high (low) redshift subsample. Thus, each subsample has a
similar $\nu_{\rm 200m}$ distribution but different redshift ranges. We use
the same measurement procedure as Section 4.1, but now applied to the two
subsamples, to get the mean pressure profile and corresponding log-derivative.

The mass evolution is studied using similarly constructed subsamples — first
we rank-order clusters by their redshift, then separate pairs of clusters into
high/low mass groups based on $M_{\rm 200m}$, and combine all groups to get
the final high/low mass subsamples. In this case, each subsample has a similar
redshift distribution but different cluster masses. The resulting mean profile
from each subsample is also compared with a theoretical prediction, where the
latter now also accounts for the specific redshift- and mass-distribution of
the subsample. Table 1 provides the mean mass and redshift of the subsamples.

The left panels of Figure 5 show the results for the mass evolution. The low
mass sample has a lower profile normalization and this is the expected
behavior from simulations, as is evident from the good agreement between our
measurement and the theoretical model for each subsample. The location of the
pressure deficit at $R_{\rm sh,\,de}/R_{\rm 200m}\approx 1.1$ is statistically
consistent across all samples, and the significance of the feature varies
between $2\lesssim\chi\lesssim 4$ (see Table 1). For the high (low) mass
sample we find a $3.8\sigma$ ($1.8\sigma$) significance, indicating the
deficit is more strongly observed for higher mass clusters. This is likely
just reflecting the fact that higher mass objects have a higher detection SNR.
The profiles for both subsamples show a clear plateau in the pressure between
$1R_{\rm 200m}<R<3R_{\rm 200m}$ and a sharp decrease beyond that; this is
consistent with the accretion shock feature described in Paper I, as was noted
before.

For certain subsamples we do not clearly resolve a second minimum in the log-
derivative that would correspond to the accretion shock. The high mass sample
in the left panels of Figure 5 is one such example. In this case, we only show
the log-derivative up to where we have reasonable constraints, i.e. before the
error bars blow up due to noise domination, and also only quote a lower limit
for the location of the accretion shock. This limit is set by computing the
location of the first maximum in the log-derivative. Any accretion shock —
identified by the second minimum — will need to be further out from the
location of the first maximum.

The right panels of Figure 5 show our results for the redshift evolution. The
inner profile of the cluster subsamples are nearly atop one another, whereas
the outskirts deviate more but remain statistically consistent. The results
seem to suggest that at low redshift the accretion shock feature may be pushed
to slightly larger radii, but we do not quantify the significance of this
behavior given we do not observe a second minimum in the log-derivative of the
low redshift subsample and thus can only place a lower limit on the accretion
shock location. The location of the pressure deficit is consistent across
subsamples, like before, while the significance varies slightly: the high
(low) redshift sample observes the pressure deficit at $2.3\sigma$
($3.2\sigma$).

A key finding from both the analyses is that the location of pressure deficit
is statistically consistent across all subsamples. It is constantly found at
$R_{\rm sh,\,de}/R_{\rm 200m}\approx 1.1$, where the log-derivative is always
steeper than in the theory. The characteristic plateauing plus subsequent
drop-off of the $y$-profiles, indicative of an accretion shock, also exists in
all subsamples.

###  4.3 Connection to Filaments

Figure 6: The stacked pressure profiles of SPT-SZ sample (left) and the mass-
and redshift-matched The300 (right) sample, computed in quadrants that are
along/across the cluster major axis (denoted as “Cl” in the legend), which we
take to be aligned with cosmic filaments. Panels adopt the same format as
Figure 3. The profiles are consistent with being stretched along the major
axis and squeezed along the minor axis, which leads to a shallower and steeper
log-derivative, respectively. The accretion shock along the minor axis is
potentially closer to the cluster than the same feature along the major axis,
and a similar feature appears to show up in the simulated cluster sample as
well. The simulations are cut off at $R<6R_{\rm 200m}$ due to the limited
scales available in the resimulated regions, and these excluded scales are
shaded out.

The morphology of the accretion shock is known to depend on the filamentary
structure connected to the galaxy cluster (e.g., Molnar et al., 2009; Lau et
al., 2015; Zhang et al., 2020; Aung et al., 2020; Zhang et al., 2021). The
amount of matter accreted onto the cluster is not isotropic and is
preferentially higher along the direction of filaments. However, the gas
within filaments is already preshocked and is unlikely cold enough to generate
an accretion shock (e.g., Zinger et al., 2016b). Thus, the accretion shock
feature is expected to be anisotropic — it is a weak-to-nonexistent feature in
directions directly towards filaments. However, the accretion shock boundary
is also elliptical and has been shown to align with the filamentary large-
scale structure (Aung et al., 2020, see Figure 1).

To test this in the data, we assume the cluster major axis — as determined by
the $y$-map image of the cluster — is preferentially aligned in the direction
with more filamentary structure. So we split this $y$-map image into
quadrants, with two quadrants containing the major axis and the other two
containing the minor axis. By reorienting and stacking the cluster quadrants
appropriately, we compute the mean pressure profile along and across the
cluster major axis. Note that the quadrants across the cluster minor axis will
still contain filamentary structures, but will contain fewer structures than
the quadrants along the major axis. The orientation of the cluster in the map
is obtained by fitting a 2D Gaussian to the $y$-map using the AstroPy library
(Astropy Collaboration et al., 2013, 2018). We only fit to pixels within
$R/R_{\rm 200m}<0.5$, and find that extending the aperture further leads to
systematic biases due to the noise-domination within individual pixels beyond
those scales.

The left panels of Figure 6 show the results of such an analysis on the SPT-SZ
data. The profile mostly follows our expectations — along the longer, major
axis, it is essentially stretched out, while along the minor axis it is
compressed. So at fixed radii, the pressure measured along the major axis
would be higher than that along the minor axis. The stretching naturally
results in a shallower log-derivative, while the squeezing results in a
steeper one. We also see features resembling an accretion shock — i.e. the
plateauing followed by a sharp drop off — but they are not very prominent. The
feature along the cluster minor axis is closer in than a similar feature along
the cluster major axis, which is consistent with results from Aung et al.
(2020). However, note that these are low significance features. We also stress
that the theoretical model plotted in Figure 6 does not account for the
ellipticity of the halos, and so should not be used to quantify the detection
significance of an observed feature. For this reason we quote neither
$\epsilon$ or $\chi$ for this analysis.

The right panels of Figure 6 show the results from the mass- and redshift-
matched cluster sample from The300. We find similar “stretching” and
“compressing” behaviors in the one-halo term, as expected of clusters that are
elliptical, not spherical, in nature. The log-derivatives follow the total
halo model theory, though there are some deviations in The300 at $r>2R_{\rm
200m}$ in the form of a second minimum. These could be the accretion shock,
but are far too weak to make any strong claims. They are however located
around $R/R_{\rm 200m}\approx 2.3$, which is where Paper I found their
accretion shock for the relaxed sample of The300 clusters at $z=0.2$. If these
two features are indeed accretion shocks, they follow the SPT-SZ data
behavior, in that the feature measured along the minor axis is closer to the
cluster than that measured along the major axis. The fact that we find some
features here, in contrast to finding no clear features in the azimuthally-
averaged profiles, also indicates that splitting profiles across cluster major
and minor axes may be a better way of searching for accretion shock features.

##  5 Comparison with External Data

Figure 7: Same format as Figure 3, but for SPT-SZ and ACTPol (left), and for
Planck and SPT-SZ + ACTPol (right). SPT-SZ and ACTPol show a consistent
pressure deficit at $R_{\rm sh,\,de}/R_{\rm 200m}\approx 1.1$, with
$3.1\sigma$ and $3.6\sigma$ significance each. Combining both SPT-SZ and
ACTPol data leads to a $4.6\sigma$ significance detection of the pressure
deficit, and a $2.7\sigma$ detection of the accretion shock. Planck has a
pressure deficit further out at $R_{\rm sh,\,de}/R_{\rm 200m}=1.62$, and this
is due to its order-of-magnitude larger, $10^{\prime}$ smoothing scale. The
unsmoothed theory prediction for Planck is shown in the dotted blue line.

We have thus far focused on the gas pressure profile in the SPT-SZ dataset. In
this section, we first repeat our pressure profile analysis using cluster
samples and corresponding $y$-maps from two other CMB experiments — the public
ACTPol and Planck surveys — in Section 5.1. Next, we also compare our SPT-SZ
results with existing splashback radius measurements — made for the same SPT-
SZ cluster catalog — in Section 5.2. Similar to before, the location and
significance of the detected features are summarized in Table 1.

###  5.1 Consistency with Alternative Datasets

Here we repeat the analysis presented in Figure 3, but for three additional
cluster samples: (i) ACTPol, and (ii) a combination of SPT-SZ + ACTPol, and;
(iii) Planck. In each case, the measurements for clusters from a particular
survey are made using only the $y$-map released by that particular survey. The
descriptions of these maps, and the data used to make them, can be found in
Section 2.

Comparisons of SPT-SZ with ACTPol and Planck provide an important validation
step of our findings as the latter two datasets contain a nearly independent
set of clusters from SPT-SZ, and their $y$-maps were made using different
techniques applied to different data. Note that while the large scales
($\ell\lesssim 1000$) in both SPT-SZ and ACTPol maps are primarily informed by
the same Planck data, the features we study are localized in real-space, and
thus their corresponding harmonic-space features will span a broad range of
scales instead of being localized to a range in $\ell$. Additionally, the
average angular distance of the pressure deficit feature from the cluster
center corresponds to $\ell\sim 2000$, where the SPT-SZ and ACTPol data are
more informative than Planck. The same for the accretion shock feature
corresponds to $\ell\sim 700$.

Figure 7 shows the results of this comparison analysis, alongside theoretical
predictions constructed in a similar fashion to previous plots. The left
column compares SPT-SZ and ACTPol, where we find statistical consistency on
both the depth and location of the pressure deficit at $R_{\rm sh,\,de}/R_{\rm
200m}\approx 1.1$ (see Table 1 for specific constraints). This provides strong
evidence that the feature is of physical origin and also provides a robust
consistency test of the feature’s radial location. The log-derivatives of
ACTPol show some semblance of a second minimum at $4-5R_{\rm 200m}$.
Furthermore, the $y$-profiles from both SPT-SZ and ACTPol showcase a
plateauing phase where the pressure is nearly constant over a wide radii
range. As previously noted, this feature is consistent with the simulation
predictions for an accretion shock as shown in Paper I.

Given that the log-derivatives of the SPT-SZ and ACTPol samples are
statistically consistent, we combine the two cluster samples and redo our
analysis. The ACTPol and SPT-SZ cluster measurements are still made using only
their respective $y-$maps, and the profile data vectors are simply combined to
get $N=N_{\rm SPT}+N_{\rm ACT}$ profiles. Note that the SPT-SZ and ACTPol
smoothing scales are slightly different ($1.25^{\prime}$ for SPT-SZ and
$1.6^{\prime}$ for ACTPol) but we have confirmed this difference does not
impact our results for the combined sample. The purple lines on the right
columns of Figure 7 show the results of the combined SPT-SZ + ACTPol sample,
where a key finding is the shock-like feature, found at $R_{\rm
sh,\,de}/R_{\rm 200m}=1.09\pm 0.06$, is now measured at $4.6\sigma$
significance. The significance of the accretion shock, which is located at
$R_{\rm sh,\,acc}/R_{\rm 200m}=4.17\pm 0.32$, also increases slightly from
$2\sigma\rightarrow 2.7\sigma$. The large reduction of uncertainties in
$R_{\rm sh,\,acc}/R_{\rm 200m}$, from $\pm 1.24$ for SPT-SZ to $\pm 0.32$ for
SPT-SZ + ACTPol, is mostly due to the minor change in the shape of the
accretion shock feature between the two cases.

Figure 7 also shows results from the Planck sample, which are significantly
different from SPT-SZ and ACTPol primarily due to both the different cluster
redshift and mass distributions (see Figure 2), and also the nearly order-of-
magnitude larger smoothing scale difference ($\approx 1^{\prime}$ vs.
$10^{\prime}$). We still see a pressure deficit in Planck, with its location
now pushed out to $R_{\rm sh,\,de}/R_{\rm 200m}=1.62\pm 0.14$ due to the
smoothing; the dotted blue line in the panel shows the theoretical prediction
when no smoothing is included and here the minimum is significantly closer to
the cluster core. We estimate the expected shift due to smoothing by measuring
the location of the first minimum in the smoothed and unsmoothed theory
curves, and find the ratio of locations is $0.5$. This implies the true,
corrected location of the measured pressure deficit is closer to $R_{\rm
sh,\,de}^{\rm corr}/R_{\rm 200m}=0.5\times 1.62=0.81$. This pressure deficit
is found at $3.0\sigma$ significance. The uncertainties in both the profile
and log-derivative are significantly smaller for Planck and this is due to the
larger smoothing in the Planck $y$-maps.

We do not see any indication of an accretion shock feature for Planck in
neither the profile nor the log-derivative, and so do not set any lower limits
on its location. We suspect this is caused by smoothing, as our tests with the
simulations have shown that smoothing at $10^{\prime}$ significantly
suppresses features that were seen at the $1.25^{\prime}$ smoothing level.
Note that there are also minor, nearly-constant vertical offsets between the
theory and the observations for Planck. This arises because the masses in the
Planck 2015 cluster catalog are biased low, as was discussed in Section 2.3.
We have only approximately corrected for this with our re-calibration step,
whereas a full re-analysis of these masses would result in better agreement
between observations and theory. For these reasons mentioned above, we do not
explore a detailed comparison of Planck with SPT-SZ and ACTPol and instead
consider the Planck results as a more qualitative comparison point.

###  5.2 Connecting Splashback and Shock Features

Dataset | $R_{\rm sh,\,de}/R_{\rm 200m}$ | $R_{\rm sh,\,acc}/R_{\rm 200m}$ | $\frac{{\rm d\ln}y}{{\rm d\ln}R}(\frac{R_{\rm sh,\,de}}{R_{\rm 200m}})$ | $\frac{{\rm d\ln}y}{{\rm d\ln}R}(\frac{R_{\rm sh,\,acc}}{R_{\rm 200m}})$ | $\chi_{\rm sh,\,de}$ | $\chi_{\rm sh,\,acc}$ | $\langle\log_{10}M_{\rm 200m}\rangle[\,{\rm M}_{\odot}]$ | $\langle z\rangle$
---|---|---|---|---|---|---|---|---
SPT-SZ | $1.08\pm 0.09$ | $4.58\pm 1.24$ | $-3.95\pm 0.50$ | $-2.77\pm 1.65$ | $3.1$ | $2.0$ | $14.9$ | $0.57$
ACTPol | $1.09\pm 0.08$ | $>2.97\pm 0.21$ | $-3.74\pm 0.32$ | — | $3.6$ | — | $14.7$ | $0.54$
Planck | $1.62\pm 0.14$ | — | $-3.04\pm 0.14$ | — | $3.0$ | — | $15.0$ | $0.25$
SPT-SZ + ACTPol | $1.09\pm 0.06$ | $4.17\pm 0.32$ | $-3.80\pm 0.27$ | $-2.40\pm 0.82$ | $4.6$ | $2.7$ | $14.8$ | $0.55$
SPT-SZ (DES-matched) | $0.93\pm 0.07$ | $>2.56\pm 0.41$ | $-3.88\pm 0.51$ | — | $3.0$ | — | $14.9$ | $0.49$
SPT-SZ (High M) | $1.10\pm 0.10$ | $>2.43\pm 0.36$ | $-4.11\pm 0.67$ | — | $3.8$ | — | $14.9$ | $0.57$
SPT-SZ (Low M) | $1.07\pm 0.13$ | $4.03\pm 0.38$ | $-3.69\pm 0.77$ | $-3.89\pm 2.48$ | $1.8$ | $1.3$ | $14.8$ | $0.57$
SPT-SZ (High z) | $1.04\pm 0.10$ | $3.80\pm 0.52$ | $-3.77\pm 0.68$ | $-3.34\pm 1.96$ | $2.3$ | $2.0$ | $14.8$ | $0.68$
SPT-SZ (Low z) | $1.13\pm 0.12$ | $>2.57\pm 0.15$ | $-4.15\pm 0.74$ | — | $3.2$ | — | $14.7$ | $0.46$
SPT-SZ (Major axis) | $1.09\pm 0.44$ | $5.44\pm 0.56$ | $-3.11\pm 0.31$ | $-3.44\pm 1.65$ | — | — | $14.9$ | $0.57$
SPT-SZ (Minor axis) | $1.09\pm 0.06$ | $3.66\pm 0.38$ | $-4.53\pm 0.91$ | $-2.49\pm 1.15$ | — | — | $14.9$ | $0.57$
Table 1: Table containing the constraints presented in this work. All
uncertainties are $\pm 1\sigma$ estimates. From left to right the columns
show: (i) the sample name, (ii - iii) location of the pressure deficit and
accretion shock, (iv - v) the value of the log-derivative at the location of
the two features, (vi - vii) detection significance of the two features,
extracted using equation (21), (viii - ix) the mean log-mass and redshift of
the cluster sample. When we do not observe a clear second minimum in the log-
derivative, we quote a lower limit for $R_{\rm sh,\,acc}$ by quantifying the
location of the first maximum. We also do not quote the log-derivative or the
detection significance in this case. For Planck, we cannot set lower limits as
we do not observe any accretion shock features. The $R_{\rm sh,\,de}$ estimate
for Planck is also strongly impacted by smoothing, making it significantly
different from SPT-SZ and ACTPol, and a smoothing-corrected value is provided
in the text. We also do not quantify the significance of features in the
major/minor axis analyses of SPT-SZ as we do not have an appropriate theory
model to use as a baseline. Figure 8: The log-derivative of the mean pressure
profile alongside that of the mean galaxy number density profile from Shin et
al. (2019). Both measurements are made on the same subset of SPT-SZ clusters,
which is different from the data used for Figure 3. See Section 5.2 for
details on sample selection. The two profiles have different shapes at all
radii, evidenced by the different log-derivatives. The pressure deficit is
located inside the splashback radius, but the locations differ by just over
$1\sigma$. The accretion shock was not measurable due to increased noise
levels from using fewer clusters in this analysis. The location of the
pressure profile features are shown in Table 1, and a lower limit for the
ratio of accretion shock radius and splashback radius is given in Equation
(23).

Simple, self-similar evolution models of structure formation predict that the
locations of the accretion shock and splashback feature coincide for a gas
adiabatic index $\gamma=5/3$ (Bertschinger, 1985). Assuming different gas
adiabatic indices in the range $\gamma\in[4/3,3]$, and also allowing the mass
accretion rate to vary, can lead to nearly factor of 2 differences in the
location of the accretion shock in self-similar solutions (Shi, 2016).
However, more complex, 3D hydrodynamical simulations have all found that the
accretion shock feature is consistently located well beyond the splashback
radius (Paper I; Lau et al., 2015; Aung et al., 2020; Zhang et al., 2020,
2021), and this is related to the MA-shocks we discussed earlier (Zhang et
al., 2020).

Here, we compare features in the gas pressure profiles presented in this work
with the dark matter splashback feature previously measured for the SPT-SZ
sample by Shin et al. (2019, S19). In S19, the splashback feature is measured
by correlating SPT-SZ galaxy clusters with a galaxy sample from the Dark
Energy Survey (DES, Flaugher, 2005) and constructing galaxy density profiles
that serve as a proxy for the shape of the dark matter radial density
distribution. The minimum in the log-derivative of that profile provides the
location of the splashback features; in this sense the splashback and shock
measurements are highly analogous.

S19 were unable to use all 516 SPT-SZ clusters as they required clusters to be
in both the SPT-SZ and DES footprints, and to be in the redshift range of the
DES galaxy sample, $0.25<z<0.7$. In order to perform a fair comparison between
works, we use the same selections as S19 on our cluster sample and redo the
stacking and log-derivative measurements for this new subsample. Specifically,
we only use clusters that are in the redshift range $0.25<z<0.7$ and that are
contained in the DES footprint. This selection lowers our sample from
$N=516\rightarrow 256$.

Figure 8 compares log-derivatives of the pressure profiles and of the galaxy
number density profiles of SPT-SZ clusters. We change our upper radial limit
to $R/R_{\rm 200m}=10\rightarrow 5$ as the factor of 2 reduction in sample
size, due to selection cuts, makes measurements of the mean pressure profile
beyond $R/R_{\rm 200m}=5$ noise-dominated. The log-derivatives, which quantify
the shape of the profile, are broadly similar within $R<R_{\rm 200m}$, and
deviate significantly at $R>R_{\rm 200m}$. This is consistent with inner
profiles being determined to zeroth order by simple gravitational evolution,
while at the outskirts the presence of shocks causes significant differences
in the behavior of infalling dark matter and gas. The pressure deficit is now
found at $R_{\rm sh,\,de}/R_{\rm 200m}=0.93\pm 0.07$ and is statistically
consistent with the location of the splashback feature, which is at $R_{\rm
sp}/R_{\rm 200m}=1.22\pm 0.25$ as quoted in S19. While $R_{\rm sh,\,de}$ in
the fiducial sample is further out than the DES-matched sample described here,
the difference is a little over $1\sigma$. Moreover, the DES-matched sample
has a significantly different redshift distribution — it excludes 148 high-
redshift clusters ($z>0.7$) and 61 low-redshift ones ($z<0.25$) — which could
cause the minor difference in $R_{\rm sh,\,de}$.

We also do not find a visible second minimum in the log-derivatives that could
indicate the presence of an accretion shock feature. However, we can still set
a lower limit on the radial location of such a feature by quantifying the
location of the first maximum, like we did in Section 4.2. We find $R_{\rm
sh,\,acc}/R_{\rm 200m}>2.56\pm 0.41$ for this lower limit. Using this value,
we can also set a lower limit on the ratio between the accretion shock radius
and the splashback radius, which is

$R_{\rm sh,\,acc}/R_{\rm sp}>2.12\pm 0.59.$ (23)

This ratio is statistically consistent with expectations from simulations,
which have been in the range $1.5\lesssim R_{\rm sh,\,acc}/R_{\rm sp}\lesssim
2.5$ (Paper I; Aung et al., 2020; Zhang et al., 2021). If we assume the true
location of the accretion shock in the subsample is roughly the same as its
location in the fiducial sample (see Section 4.1.2 or Table 4), then the ratio
is closer to $R_{\rm sh,\,acc}/R_{\rm sp}\approx 4$.

##  6 Summary and Discussion

The outskirts of galaxy clusters are where bound material within the cluster
interacts with infalling material that is accreted from the filamentary large-
scale structure. Naturally, these outskirts are dynamically active and contain
many instances of shocks. The study of such shocks is of particular interest
given the impact of shocks on a wide array of cluster science; from their
relevance to cosmology, via their connection to cluster hydrostatic mass
estimates, cluster dynamical states, and modelling of tSZ auto/cross-
correlations, to their relevance for astrophysics, via their influence on
galaxy formation, radio relics, magnetic fields, and cosmic ray production.
These cluster outskirts have only recently become observable due to the
availability of large samples of clusters and of $y-$maps with improved
angular resolution and noise levels; both are necessary advances for making
meaningful measurements in the noise-dominated regimes. In this work, we use
516 galaxy clusters from the SPT-SZ survey, and search for features in the
stacked gas pressure profiles between $0.3<R/R_{\rm 200m}<10$. Our key
quantitative results are all provided in Table 1, and our main findings are
summarized below:

  * •

There is a pressure deficit at $R_{\rm sh,\,de}/R_{\rm 200m}=1.08\pm 0.09$
detected at $3.1\sigma$ (Figure 3), and is consequently missing in both
simulations and analytic predictions (Figure 4). Past works show that such a
feature can be expected from a shock-driven thermal non-equilibrium between
electrons and ions, and also suggest that the lack of a similar feature in the
simulation-calibrated theory and The300 simulations is due to the simulations’
assumption that electrons and ions are in thermal equilibrium at all times.

  * •

The mean pressure profiles undergo a plateau phase where the pressure is
constant between $1\lesssim R/R_{\rm 200m}\lesssim 3$ and then drops sharply.
The corresponding minimum in the log-derivative is located at $R_{\rm
sh,\,acc}/R_{\rm 200m}=4.58\pm 1.24$, and measured at a low ($2\sigma$)
significance, but the qualitative features are consistent with predictions of
an accretion shock from Paper I.

  * •

Both the deficit and the plateau features are present when splitting the
sample based on mass or redshift, and do not show any statistically
significant evolution with either quantity. The pressure deficit is found both
along and across the cluster major axis, and does not show any significant
differences either. There is a mild indication that the boundary shell formed
by the accretion shock is elliptical and points in the same direction as the
cluster major axis (and thus, in the direction of filamentary structure).

  * •

Comparisons with similarly measured pressure profiles from ACTPol and Planck
show pressure deficits that are quantitatively consistent with the deficit
seen in SPT-SZ. The ACTPol pressure deficit in particular has the same
location and depth as the SPT-SZ feature. ACTPol also has qualitative features
indicative of an accretion shock, but we see no clear second minimum in the
log-derivatives.

  * •

Comparing the pressure profiles with the splashback radius estimate of Shin et
al. (2019), we find the location of the pressure deficit is statistically
consistent with the splashback radius, $R_{\rm sp}$. The additional selections
needed for this comparison deteriorate the SNR of our measurements, and so we
are unable to resolve the accretion shock. However, we place a lower limit on
the accretion shock location, which translates to a lower limit on the ratio
of shock radius to splashback radius. We find $R_{\rm sh,\,acc}/R_{\rm
sp}>2.12\pm 0.59$, which is consistent with previous simulation studies.

Our work shows that the pressure deficit feature can already be measured with
strong significance using existing data. Future studies can probe its
existence in other SZ datasets to perform further consistency checks at higher
precision and/or, more interestingly, probe the physical origin of the deficit
as well. In this work, we have discussed how the deficit could arise from a
shock-induced thermal non-equilibrium between ions and electrons. However,
other physical processes may also resolve this difference. More detailed
analysis of this feature — on both observational and simulation fronts — will
be necessary to better understand the thermodynamic structure of the cluster
outskirts and update our models/assumptions appropriately.

Measurements of the accretion shock feature, on the other hand, are currently
still dominated by the instrument noise. We have partially worked around this
issue by stacking $N\sim 10^{3}$ clusters, which significantly lowers the
noise amplitude, but this still results in only a weak measurement of the
feature. Ongoing and future SZ surveys — like SPT-3G (Benson et al., 2014),
Advanced ACT (Henderson et al., 2016), Simons Observatory (Ade et al., 2019),
and CMB-S4 (Abazajian et al., 2019) — will both have a higher sensitivity and
a larger sample of clusters. Both factors will greatly improve the ability to
make precise measurements of the stacked mean pressure profile, leading to
more precise measurements of the accretion shock. It will also make it
possible to better study the dependence of the accretion shock on redshift,
mass, and orientation. Studying the dependence on orientation will also
provide insight into how shocks respond to the dynamics of mass accretion,
which differs significantly between directions toward filaments and toward
voids.

Another consideration — particularly regarding the accretion shock — is the
relaxation state of the cluster. Studies of shocks will benefit from
preferentially selecting systems where the shock feature evolves more self-
similarly with cluster radii, i.e. relaxed systems. A potential avenue is to
make SZ measurements around X-ray selected clusters, given the latter are
preferentially more relaxed. Such a technique is viable with the current
eROSITA All-sky X-ray mission (Merloni et al., 2012) which will detect
$N\approx 10^{5}$ galaxy groups and clusters, and of a median redshift,
$z=0.35$ (Pillepich et al., 2012). Selecting relaxed clusters could also be
done using the SZ map morphology alone (e.g., Capalbo et al., 2021). Any such
selection would allow for an observational measurement that more closely
mimics the study of Paper I and improves the strength of the observed
accretion shock feature in stacked profiles.

## Acknowledgements

We thank Tara DaCunha and Abigail Lee for their initial exploratory analyses
of tSZ profiles in the Planck dataset. We thank Greg Bryan, Damiano Caprioli,
Mark Devlin, Andrey Kravtsov, Congyao Zhang, and Irina Zhuravleva for helpful
discussions on gas physics in clusters, Shivam Pandey for discussions on the
theoretical modelling of the tSZ profiles, Tae-Hyeon Shin for providing us
with the SPT-SZ log-derivative curves from S19, and Maya Mallaby-Kay for
kindly helping us navigate the Pixell software library. We also thank the
anonymous referee for raising useful points that added to the discussion
presented here. Finally, we are grateful to the SPT-SZ, ACTPol, and Planck
collaborations for making their data products publicly available to the
community. This project strongly benefited from the free-flow of information
across surveys and collaborations. We additionally thank Colin Hill and Mathew
Madhavacheril for their generous support in including the ACT data in this
study.

DA is supported by the National Science Foundation Graduate Research
Fellowship under Grant No. DGE 1746045. CC is supported by the Henry Luce
Foundation and DOE grant DE-SC0021949. BJ is supported in part by NASA ATP
Grant No. NNH17ZDA001N and DOE Grant No. DE-SC0007901. WC is supported by the
European Research Council under grant number 670193 and by the STFC AGP Grant
ST/V000594/1. He further acknowledges the science research grants from the
China Manned Space Project with NO. CMS-CSST-2021-A01 and CMS-CSST-2021-B01.
LDM is supported by the ERC-StG ‘ClustersXCosmo’ grant agreement 716762. CR
acknowledges support from the Australian Research Council Discovery Projects
scheme (DP200101068). AS is supported by the ERC-StG ‘ClustersXCosmo’ grant
agreement 716762, by the FARE-MIUR grant ’ClustersXEuclid’ R165SBKTMA, and by
INFN InDark Grant. The South Pole Telescope program is supported by the
National Science Foundation (NSF) through award OPP-1852617

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

## Data Availability

All data used in our analyses of SPT-SZ, ACTPol, and Planck are publicly
available at the repositories linked to in this paper. Data products from
The300 simulations are not hosted on a public repository, but the interested
reader is encouraged to reach out to The300 Collaboration for data access.

The code used to generate the theoretical tSZ profile of a halo, including
both one-halo and two-halo contributions, is made available at
https://github.com/DhayaaAnbajagane/tSZ_Profiles. This repository also
contains the profiles and log-derivatives shown in the figures of this work.

## References

  * Abazajian et al. (2019) Abazajian K., et al., 2019, arXiv e-prints, p. arXiv:1907.04473
  * Ackermann et al. (2014) Ackermann M., et al., 2014, ApJ, 787, 18
  * Ade et al. (2019) Ade P., et al., 2019, J. Cosmology Astropart. Phys., 2019, 056
  * Adhikari et al. (2014) Adhikari S., Dalal N., Chamberlain R. T., 2014, J. Cosmology Astropart. Phys., 2014, 019
  * Adhikari et al. (2020) Adhikari S., et al., 2020, arXiv e-prints, p. arXiv:2008.11663
  * Akahori & Yoshikawa (2010) Akahori T., Yoshikawa K., 2010, PASJ, 62, 335
  * Akahori & Yoshikawa (2012) Akahori T., Yoshikawa K., 2012, PASJ, 64, 12
  * Akamatsu et al. (2011) Akamatsu H., Hoshino A., Ishisaki Y., Ohashi T., Sato K., Takei Y., Ota N., 2011, PASJ, 63, S1019
  * Akamatsu et al. (2016) Akamatsu H., et al., 2016, A&A, 593, L7
  * Allen et al. (2011) Allen S. W., Evrard A. E., Mantz A. B., 2011, ARA&A, 49, 409
  * Amodeo et al. (2021) Amodeo S., et al., 2021, Phys. Rev. D, 103, 063514
  * Anbajagane et al. (2021) Anbajagane D., Evrard A. E., Farahi A., 2021, arXiv e-prints, p. arXiv:2109.02713
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Aung et al. (2020) Aung H., Nagai D., Lau E. T., 2020, arXiv e-prints, p. arXiv:2012.00977
  * Avestruz et al. (2015) Avestruz C., Nagai D., Lau E. T., Nelson K., 2015, ApJ, 808, 176
  * Basu et al. (2016) Basu K., Sommer M., Erler J., Eckert D., Vazza F., Magnelli B., Bertoldi F., Tozzi P., 2016, ApJ, 829, L23
  * Battaglia et al. (2012) Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., 2012, ApJ, 758, 75
  * Battaglia et al. (2016) Battaglia N., et al., 2016, J. Cosmology Astropart. Phys., 2016, 013
  * Baxter et al. (2017) Baxter E., et al., 2017, ApJ, 841, 18
  * Baxter et al. (2021) Baxter E. J., Adhikari S., Vega-Ferrero J., Cui W., Chang C., Jain B., Knebe A., 2021, arXiv e-prints, p. arXiv:2101.04179
  * Beck et al. (2016) Beck A. M., et al., 2016, MNRAS, 455, 2110
  * Benson et al. (2014) Benson B. A., et al., 2014. (arXiv:1407.2973)
  * Bertschinger (1985) Bertschinger E., 1985, ApJS, 58, 39
  * Blandford & Eichler (1987) Blandford R., Eichler D., 1987, Phys. Rep., 154, 1
  * Bleem et al. (2015) Bleem L. E., et al., 2015, ApJS, 216, 27
  * Bleem et al. (2020) Bleem L. E., et al., 2020, ApJS, 247, 25
  * Bleem et al. (2021) Bleem L. E., et al., 2021, arXiv e-prints, p. arXiv:2102.05033
  * Bocquet et al. (2019) Bocquet S., et al., 2019, ApJ, 878, 55
  * Boselli et al. (2021) Boselli A., Fossati M., Sun M., 2021, arXiv e-prints, p. arXiv:2109.13614
  * Brunetti & Jones (2014) Brunetti G., Jones T. W., 2014, International Journal of Modern Physics D, 23, 1430007
  * Capalbo et al. (2021) Capalbo V., De Petris M., De Luca F., Cui W., Yepes G., Knebe A., Rasia E., 2021, MNRAS, 503, 6155
  * Carlstrom et al. (2002) Carlstrom J. E., Holder G. P., Reese E. D., 2002, ARA&A, 40, 643
  * Carlstrom et al. (2011) Carlstrom J. E., et al., 2011, PASP, 123, 568
  * Chang et al. (2018) Chang C., et al., 2018, ApJ, 864, 83
  * Child et al. (2018) Child H. L., Habib S., Heitmann K., Frontiere N., Finkel H., Pope A., Morozov V., 2018, ApJ, 859, 55
  * Chisari et al. (2019) Chisari N. E., et al., 2019, ApJS, 242, 2
  * Cui et al. (2018) Cui W., et al., 2018, MNRAS, 480, 2898
  * Dacunha et al. (2021) Dacunha T., Belyakov M., Adhikari S., Shin T.-h., Goldstein S., Jain B., 2021, arXiv e-prints, p. arXiv:2111.06499
  * Delabrouille & Cardoso (2009) Delabrouille J., Cardoso J. F., 2009, Diffuse Source Separation in CMB Observations. pp 159–205, doi:10.1007/978-3-540-44767-2_6
  * Di Mascolo et al. (2019a) Di Mascolo L., Churazov E., Mroczkowski T., 2019a, MNRAS, 487, 4037
  * Di Mascolo et al. (2019b) Di Mascolo L., et al., 2019b, A&A, 628, A100
  * Diemer (2018) Diemer B., 2018, ApJS, 239, 35
  * Diemer & Joyce (2019) Diemer B., Joyce M., 2019, ApJ, 871, 168
  * Diemer & Kravtsov (2014a) Diemer B., Kravtsov A. V., 2014a, ApJ, 789, 1
  * Diemer & Kravtsov (2014b) Diemer B., Kravtsov A. V., 2014b, ApJ, 789, 1
  * Dolag et al. (2008) Dolag K., Bykov A. M., Diaferio A., 2008, Space Sci. Rev., 134, 311
  * Donnert et al. (2018) Donnert J., Vazza F., Brüggen M., ZuHone J., 2018, Space Sci. Rev., 214, 122
  * Drury (1983) Drury L. O., 1983, Reports on Progress in Physics, 46, 973
  * Erler et al. (2018) Erler J., Basu K., Chluba J., Bertoldi F., 2018, MNRAS, 476, 3360
  * Ettori & Fabian (1998) Ettori S., Fabian A. C., 1998, MNRAS, 293, L33
  * Flaugher (2005) Flaugher B., 2005, International Journal of Modern Physics A, 20, 3121
  * Fox & Loeb (1997) Fox D. C., Loeb A., 1997, ApJ, 491, 459
  * Gatti et al. (2021) Gatti M., et al., 2021, arXiv e-prints, p. arXiv:2108.01600
  * Ghirardini et al. (2017) Ghirardini V., Ettori S., Amodeo S., Capasso R., Sereno M., 2017, A&A, 604, A100
  * Ghirardini et al. (2018) Ghirardini V., Ettori S., Eckert D., Molendi S., Gastaldello F., Pointecouteau E., Hurier G., Bourdin H., 2018, A&A, 614, A7
  * Ha et al. (2018) Ha J.-H., Ryu D., Kang H., 2018, ApJ, 857, 26
  * Hajian et al. (2013) Hajian A., Battaglia N., Spergel D. N., Bond J. R., Pfrommer C., Sievers J. L., 2013, J. Cosmology Astropart. Phys., 2013, 064
  * Haridasu et al. (2021) Haridasu B. S., Karmakar P., De Petris M., Cardone V. F., Maoli R., 2021, arXiv e-prints, p. arXiv:2111.01101
  * Henderson et al. (2016) Henderson S. W., et al., 2016, Journal of Low Temperature Physics, 184, 772
  * Hill & Pajer (2013) Hill J. C., Pajer E., 2013, Phys. Rev. D, 88, 063526
  * Hill et al. (2018) Hill J. C., Baxter E. J., Lidz A., Greco J. P., Jain B., 2018, Phys. Rev. D, 97, 083501
  * Hilton et al. (2021) Hilton M., et al., 2021, ApJS, 253, 3
  * Hojjati et al. (2017) Hojjati A., et al., 2017, MNRAS, 471, 1565
  * Hong et al. (2014) Hong S. E., Ryu D., Kang H., Cen R., 2014, ApJ, 785, 133
  * Hong et al. (2015) Hong S. E., Kang H., Ryu D., 2015, ApJ, 812, 49
  * Horowitz & Seljak (2017) Horowitz B., Seljak U., 2017, MNRAS, 469, 394
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Hurier et al. (2013) Hurier G., Macías-Pérez J. F., Hildebrandt S., 2013, A&A, 558, A118
  * Hurier et al. (2019) Hurier G., Adam R., Keshet U., 2019, A&A, 622, A136
  * Ishiyama et al. (2020) Ishiyama T., et al., 2020, arXiv e-prints, p. arXiv:2007.14720
  * Jarvis et al. (2004) Jarvis M., Bernstein G., Jain B., 2004, MNRAS, 352, 338
  * Kang et al. (2007) Kang H., Ryu D., Cen R., Ostriker J. P., 2007, ApJ, 669, 729
  * Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
  * Knollmann & Knebe (2009) Knollmann S. R., Knebe A., 2009, ApJS, 182, 608
  * Kravtsov & Borgani (2012) Kravtsov A. V., Borgani S., 2012, ARA&A, 50, 353
  * Lau et al. (2015) Lau E. T., Nagai D., Avestruz C., Nelson K., Vikhlinin A., 2015, ApJ, 806, 68
  * Lau et al. (2021) Lau E. T., Hearin A. P., Nagai D., Cappelluti N., 2021, MNRAS, 500, 1029
  * Ma et al. (2015) Ma Y.-Z., Van Waerbeke L., Hinshaw G., Hojjati A., Scott D., Zuntz J., 2015, J. Cosmology Astropart. Phys., 2015, 046
  * Madhavacheril et al. (2020) Madhavacheril M. S., et al., 2020, Phys. Rev. D, 102, 023534
  * Mansfield et al. (2017) Mansfield P., Kravtsov A. V., Diemer B., 2017, ApJ, 841, 34
  * McDonald et al. (2014) McDonald M., et al., 2014, ApJ, 794, 67
  * McKinney (2011) McKinney W., 2011, Python for High Performance and Scientific Computing, 14
  * Merloni et al. (2012) Merloni A., et al., 2012, arXiv e-prints, p. arXiv:1209.3114
  * Miniati et al. (2000) Miniati F., Ryu D., Kang H., Jones T. W., Cen R., Ostriker J. P., 2000, ApJ, 542, 608
  * Miniati et al. (2001) Miniati F., Jones T. W., Ryu D., Kang H., 2001, in International Cosmic Ray Conference. pp 2062–2065
  * Molnar et al. (2009) Molnar S. M., Hearn N., Haiman Z., Bryan G., Evrard A. E., Lake G., 2009, ApJ, 696, 1640
  * More et al. (2015) More S., Diemer B., Kravtsov A. V., 2015, ApJ, 810, 36
  * More et al. (2016) More S., et al., 2016, ApJ, 825, 39
  * Mroczkowski et al. (2019) Mroczkowski T., et al., 2019, Space Sci. Rev., 215, 17
  * Murata et al. (2020) Murata R., Sunayama T., Oguri M., More S., Nishizawa A. J., Nishimichi T., Osato K., 2020, PASJ, 72, 64
  * Nelson et al. (2014) Nelson K., Lau E. T., Nagai D., 2014, ApJ, 792, 25
  * Nelson et al. (2016) Nelson D., Genel S., Pillepich A., Vogelsberger M., Springel V., Hernquist L., 2016, MNRAS, 460, 2881
  * O’Neil et al. (2021) O’Neil S., Barnes D. J., Vogelsberger M., Diemer B., 2021, MNRAS, 504, 4649
  * Osato et al. (2018) Osato K., Flender S., Nagai D., Shirasaki M., Yoshida N., 2018, MNRAS, 475, 532
  * Osato et al. (2020) Osato K., Shirasaki M., Miyatake H., Nagai D., Yoshida N., Oguri M., Takahashi R., 2020, MNRAS, 492, 4780
  * Pandey et al. (2019) Pandey S., et al., 2019, Phys. Rev. D, 100, 063519
  * Pandey et al. (2020) Pandey S., Baxter E. J., Hill J. C., 2020, Phys. Rev. D, 101, 043525
  * Pandey et al. (2021) Pandey S., et al., 2021, arXiv e-prints, p. arXiv:2108.01601
  * Pfrommer et al. (2006) Pfrommer C., Springel V., Enßlin T. A., Jubelgas M., 2006, MNRAS, 367, 113
  * Pfrommer et al. (2007) Pfrommer C., Enßlin T. A., Springel V., Jubelgas M., Dolag K., 2007, MNRAS, 378, 385
  * Pillepich et al. (2012) Pillepich A., Porciani C., Reiprich T. H., 2012, MNRAS, 422, 44
  * Planck Collaboration et al. (2013) Planck Collaboration et al., 2013, A&A, 550, A131
  * Planck Collaboration et al. (2016a) Planck Collaboration et al., 2016a, A&A, 594, A22
  * Planck Collaboration et al. (2016b) Planck Collaboration et al., 2016b, A&A, 594, A27
  * Planelles & Quilis (2013) Planelles S., Quilis V., 2013, MNRAS, 428, 1643
  * Planelles et al. (2021) Planelles S., et al., 2021, arXiv e-prints, p. arXiv:2108.09670
  * Power et al. (2020) Power C., et al., 2020, MNRAS, 491, 3923
  * Pratt et al. (2021) Pratt C. T., Qu Z., Bregman J. N., 2021, arXiv e-prints, p. arXiv:2105.01123
  * Quilis et al. (1998) Quilis V., M J., Ibáñez ., Sáez D., 1998, ApJ, 502, 518
  * Rasia et al. (2015) Rasia E., et al., 2015, ApJ, 813, L17
  * Romero et al. (2017) Romero C. E., et al., 2017, ApJ, 838, 86
  * Romero et al. (2018) Romero C., et al., 2018, A&A, 612, A39
  * Rudd & Nagai (2009) Rudd D. H., Nagai D., 2009, ApJ, 701, L16
  * Ryu et al. (2003) Ryu D., Kang H., Hallman E., Jones T. W., 2003, ApJ, 593, 599
  * Sakstein et al. (2016) Sakstein J., Wilcox H., Bacon D., Koyama K., Nichol R. C., 2016, J. Cosmology Astropart. Phys., 2016, 019
  * Sayers et al. (2013) Sayers J., et al., 2013, ApJ, 768, 177
  * Sayers et al. (2016) Sayers J., et al., 2016, ApJ, 832, 26
  * Schaal & Springel (2015) Schaal K., Springel V., 2015, MNRAS, 446, 3992
  * Schaal et al. (2016) Schaal K., et al., 2016, MNRAS, 461, 4441
  * Schaan et al. (2021) Schaan E., et al., 2021, Phys. Rev. D, 103, 063513
  * Shi (2016) Shi X., 2016, MNRAS, 461, 1804
  * Shi & Komatsu (2014) Shi X., Komatsu E., 2014, MNRAS, 442, 521
  * Shi et al. (2015) Shi X., Komatsu E., Nelson K., Nagai D., 2015, MNRAS, 448, 1020
  * Shin et al. (2019) Shin T., et al., 2019, MNRAS, 487, 2900
  * Shin et al. (2021) Shin T., et al., 2021, arXiv e-prints, p. arXiv:2105.05914
  * Shirasaki et al. (2020) Shirasaki M., Lau E. T., Nagai D., 2020, MNRAS, 491, 235
  * Skillman et al. (2008) Skillman S. W., O’Shea B. W., Hallman E. J., Burns J. O., Norman M. L., 2008, ApJ, 689, 1063
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Sunyaev & Zeldovich (1972) Sunyaev R. A., Zeldovich Y. B., 1972, Comments on Astrophysics and Space Physics, 4, 173
  * Tanimura et al. (2021) Tanimura H., Douspis M., Aghanim N., Salvati L., 2021, arXiv e-prints, p. arXiv:2110.08880
  * Terukina et al. (2014) Terukina A., Lombriser L., Yamamoto K., Bacon D., Koyama K., Nichol R. C., 2014, J. Cosmology Astropart. Phys., 2014, 013
  * Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
  * Tinker et al. (2010) Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878
  * Van der Walt et al. (2011) Van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science and Engineering, 13, 22
  * Vazza et al. (2009) Vazza F., Brunetti G., Gheller C., 2009, MNRAS, 395, 1333
  * Vazza et al. (2011) Vazza F., Dolag K., Ryu D., Brunetti G., Gheller C., Kang H., Pfrommer C., 2011, MNRAS, 418, 960
  * Vazza et al. (2012) Vazza F., Brüggen M., van Weeren R., Bonafede A., Dolag K., Brunetti G., 2012, MNRAS, 421, 1868
  * Vazza et al. (2013) Vazza F., Brüggen M., Gheller C., 2013, MNRAS, 428, 2366
  * Vazza et al. (2014) Vazza F., Gheller C., Brüggen M., 2014, MNRAS, 439, 2662
  * Vikram et al. (2017) Vikram V., Lidz A., Jain B., 2017, MNRAS, 467, 2315
  * Vink et al. (2015) Vink J., Broersen S., Bykov A., Gabici S., 2015, A&A, 579, A13
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Walker et al. (2019) Walker S., et al., 2019, Space Sci. Rev., 215, 7
  * Wechsler et al. (2002) Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52
  * Wilcox et al. (2015) Wilcox H., et al., 2015, MNRAS, 452, 1171
  * Wong & Sarazin (2009) Wong K.-W., Sarazin C. L., 2009, ApJ, 707, 1141
  * Xhakaj et al. (2020) Xhakaj E., Diemer B., Leauthaud A., Wasserman A., Huang S., Luo Y., Adhikari S., Singh S., 2020, MNRAS, 499, 3534
  * Zhang et al. (2019) Zhang C., Churazov E., Forman W. R., Lyskova N., 2019, MNRAS, 488, 5259
  * Zhang et al. (2020) Zhang C., Churazov E., Dolag K., Forman W. R., Zhuravleva I., 2020, MNRAS, 494, 4539
  * Zhang et al. (2021) Zhang C., Zhuravleva I., Kravtsov A., Churazov E., 2021, MNRAS, 506, 839
  * Zhu et al. (2021) Zhu Z., et al., 2021, in Galaxy Cluster Formation II (GCF 2021) - Virtual Workshop. p. 56 (arXiv:2107.01663), doi:10.5281/zenodo.5013757
  * Zinger et al. (2016a) Zinger E., Dekel A., Kravtsov A. V., Nagai D., 2016a, arXiv e-prints, p. arXiv:1610.02644
  * Zinger et al. (2016b) Zinger E., Dekel A., Birnboim Y., Kravtsov A., Nagai D., 2016b, MNRAS, 461, 412
  * Zürcher & More (2019) Zürcher D., More S., 2019, ApJ, 874, 184

##  Appendix A Tests of Robustness

Figure 9: The response of our SPT-SZ results to changes in the analysis
pipeline. From left to right we test: (i) smoothing kernels, (ii) weighting
schemes, and (iii) bin widths. For each panel, we denote our fiducial choice
as “Fid”, and also show the theory curves, which are the same in every column.
The results change slightly when increasing our Gaussian smoothing scale by a
factor $2$, while the rest of the method variations have no relevant effects.

We have made certain methodological choices in our analysis, and here we test
the sensitivity of our results to these choices. In Figure 9, we show the
response of the SPT-SZ mean $y-$profile and log-derivative (previously shown
in Figure 3) to variations in the (i) smoothing kernels, (ii) weighting
schemes, and; (iii) bin widths. We detail our findings below.

Smoothing: In general, all smoothing kernel choices still resolve the shock
features at the same detection significance. When using a narrower Gaussian
with width $\sigma_{\rm\ln R}=0.08<{\rm d\ln}R\approx 0.11$ —which is
effectively an unsmoothed case — the log-derivative is noisier, particularly
further out from the cluster center. Conversely, using a wider Gaussian with
$\sigma_{\rm\ln R}=0.32$ results in oversmoothing which both changes the
profile significantly and dampens the deviations from the theory model as
showcased in the fiducial SPT-SZ result (Figure 3). The oversmoothing also
shifts the profiles and log-derivative to larger scales. Using a Savitsky-
Golay filter — the original choice of Paper I — still resolves the shock. Here
we used a window length of 7 and a polynomial order of 3 for the input
parameters of the Savitsky-Golay filter.

Weights: Our results are quite insensitive to the choices of weights used when
stacking profiles. Using $w=\texttt{SNR}$ leads to the highest detection
significance for both features, but the improvement is miniscule. Even in the
unweighted case, we can resolve both shocks at close to the fiducial detection
significance. Note that for SPT-SZ, $\texttt{SNR}\in[4.5,42]$, and so
$\sqrt{\texttt{SNR}}\in[2,6.5]$ which is closer to the unweighted scenario
given the smaller variation in weights. The mean profile beyond $R/R_{\rm
200m}>5$ does show some differences but these are not statistically
significant. We do not consider the case of $w=\texttt{SNR}^{2}$ as the wide
dynamic range of the resulting weights leads to a small fraction of clusters
(10%) dominating the final result.

Binning: Finally, changing the bin widths by factors of 2 has no impact on the
stacked profile and log-derivative. Note that we continue to use a Gaussian
smoothing step with $\sigma_{\rm\ln R}=0.16$ in all three cases.