# Surface Brightness Fluctuations in Two SPT Clusters: A Pilot Study

Charles E. Romero<fnref n="1" />, Massimo Gaspari<fnref n="2" />, Gerrit Schellenberger<fnref n="1" />, Bradford A. Benson<fnref n="3" />,<fnref n="4" />,<fnref n="5" />, Lindsey E. Bleem<fnref n="5" />,<fnref n="6" />, Esra Bulbul<fnref n="7" />, Matthias Klein<fnref n="8" />, Ralph Kraft<fnref n="1" />, Paul Nulsen<fnref n="1" />,<fnref n="9" />, Christian L. Reichardt<fnref n="10" />, Laura Salvati<fnref n="11" />, Taweewat Somboonpanyakul<fnref n="12" />,<fnref n="13" />, and Yuanyuan Su<fnref n="14" />

## Abstract

Studies of surface brightness fluctuations in the intracluster medium present an indirect probe of turbulent properties such as the turbulent velocities, injection scales, and the slope of the power spectrum of fluctuations toward smaller scales. With the advancement of Sunyaev-Zel'dovich (SZ) studies and surveys relative to X-ray observations, we seek to investigate surface brightness fluctuations in a sample of South Pole Telescope (SPT)-SZ clusters which also have archival XMM-Newton data. Here we present a pilot study of two typical clusters in that sample: SPT-CLJ0232-4421 and SPT-CLJ0638-5358. We infer injection scales larger than 500 kpc in both clusters and Mach numbers ≈ 0.5 in SPT-CLJ0232-4421 and Mach numbers ≈ 0.6-1.6 in SPT-CLJ0638-5358, which has a known shock. We find hydrostatic bias values for $M_{500}$ less than 0.2 for SPT-CLJ0232-4421 and less than 0.1 for SPT-CLJ0638-5358. These results show the importance to assess quantitative values via a detailed multiwavelength approach and suggest that the drivers of turbulence may occur at quite large scales.

*Unified Astronomy Thesaurus concepts: Galaxy clusters (584); Intracluster medium (858)*

## 1. Introduction

Surface brightness fluctuations of the intracluster medium (ICM) provide a window to probe the turbulent properties of the ICM. The surface brightness fluctuations of the ICM in the X-ray and millimeter (via the Sunyaev-Zel'dovich (SZ) effect; Sunyaev & Zel'dovich 1972) regimes can be deprojected to thermodynamic quantities because the ICM is optically thin. The SZ signal, parameterized in terms of Compton y, is proportional to the electron pressure, $P_e$ along the line of sight:

$$y = \frac{\sigma_T}{m_e c^2} \int P_e dz, \quad (1)$$

where $\sigma_T$ is the Thomson cross section, $m_e$ is the electron mass, c is the speed of light, and z is taken to be the axis along the line of sight. The X-ray surface brightness is proportional to the emission integral (EI):

$$\text{EI} \equiv \int n_p n_e dz, \quad (2)$$

where $n_p$ is the proton density and $n_e$ is the electron density. This value is similar to the classically defined emission

measure (EM);<fnref n="15" /> one can consider the EI to be the EM per unit area. The proportionality is modulated by the cooling function<fnref n="16" /> defined for a given energy band b, $\Lambda_b(Z, T_g)$, which is a function of metallicity Z and gas temperature $T_g$. Noting that $n_p$ will have some proportionality to $n_e$ depending on Z, we see that EM $\propto n_e^2$. When measuring the X-ray surface brightness of hot clusters (such as those in our sample) in a soft-energy band (e.g., 0.5-2 keV), the cooling function (and thus the surface brightness) is relatively temperature insensitive (e.g., Sarazin 1988).

Due to these proportionalities, the surface brightness of the SZ and X-ray (in soft bands) can be deprojected to determine pressure (e.g., Khatri & Gaspari 2016) and density (e.g., Churazov et al. 2012) fluctuations, respectively. This deprojection can be done in Fourier space (Peacock 1999) such that a power spectrum of the surface brightness fluctuations can be transformed into a power spectrum of the respective thermodynamic quantity governing the surface brightness signal. The power spectrum of thermodynamic quantities is particularly informative about the physics of the

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA; charles.romero@cfa.harvard.edu</footnote>

<footnote n="2">²Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544-1001, USA</footnote>

<footnote n="3">³Fermi National Accelerator Laboratory, MS209, P.O. Box 500, Batavia, IL 60510, USA</footnote>

<footnote n="4">⁴Department of Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue, Chicago IL 60637, USA</footnote>

<footnote n="5">⁵Kavli Institute for Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA</footnote>

<footnote n="6">⁶High Energy Physics Division, Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA</footnote>

<footnote n="7">⁷Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse 1, 85748 Garching, Germany</footnote>

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

<footnote n="9">⁹ICRAR, University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia</footnote>

<footnote n="10">¹⁰School of Physics, University of Melbourne, Parkville, VIC 3010, Australia</footnote>

<footnote n="11">¹¹Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatial, F-91405 Orsay, France</footnote>

<footnote n="12">¹²Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, 452 Lomital Mall, Stanford, CA 94305, USA</footnote>

<footnote n="13">¹³Department of Physics, Faculty of Science, Chulalongkorn University 254 Phayathai Road, Pathumwan, Bangkok 10330, Thailand</footnote>

<footnote n="14">¹⁴Department of Physics and Astronomy, University of Kentucky, 505 Rose Street, Lexington, KY 40506, USA</footnote>

<footnote>Received 2023 November 10; revised 2024 January 31; accepted 2024 February 13; published 2024 July 17</footnote>

<footnote n="15">¹⁵ EM is the integral of the expression in Equation (2) over a volume, $dV$, rather than $dz$.</footnote>

<footnote n="16">¹⁶ Formally the cooling function is bolometric; i.e., it encapsulates the emission over all frequencies.</footnote>
ICM plasma as it will depend on the injection scale, turbulent motions, and transport properties such as viscosity.

In practice, probing the viscosity (turbulent cascade) is very difficult due to the required dynamic range, both in angular scales and amplitudes. That is, these measurements generally require high resolution and high sensitivity (e.g., Churazov et al. 2012). In contrast, estimating the turbulent velocity is less stringent, since in a stratified atmosphere we expect a linear relation between thermodynamic fluctuations and turbulent velocities, v, (i.e., Mach number, v/c_s, that quantifies gas motion relative to the speed of sound, c_s; Gaspari & Churazov 2013; Gaspari et al. 2014; Zhuravleva et al. 2014, 2023; Simonte et al. 2022), where this relation is dominated by the peak of the amplitude spectra. Thus, assuming an injection scale, Mach numbers can be discerned over a sample of clusters via Fourier analysis of surface brightness fluctuations (e.g., Eckert et al. 2017). Similarly, Hofmann et al. (2016) binned regions by photon count and obtained deprojected thermodynamic profiles and took the scatter in those profiles to reflect the fluctuations and then infer Mach numbers.

As a complementary application of power spectra of ICM fluctuations, we can investigate the hydrostatic mass bias within clusters. From simulations and observations, many studies find average hydrostatic mass biases between 10% to 30% (e.g., Mahdavi et al. 2013; Applegate et al. 2014; Martino et al. 2014; Hoekstra et al. 2015; Mantz et al. 2016; Hurier & Angulo 2018; Siegel et al. 2018), while a hydrostatic mass bias of over 40% would be needed to resolve the tension between Planck cosmology and Planck cluster number counts (Planck Collaboration et al. 2016a; see also Pratt et al. 2019 for a recent review of cluster mass estimates and their cosmological impact). The majority of recent observational constraints on hydrostatic mass bias values compare some hydrostatic mass estimate (from X-ray or SZ data) and a total mass estimate from weak lensing, though one can also estimate the mass bias from other methods; e.g., calculating total mass via caustic methods (e.g., Maughan et al. 2016) or relating an observed gas fraction to an assume gas fraction (e.g., Allen et al. 2008; Eckert et al. 2019; Wicker et al. 2023).

In order to infer a hydrostatic mass bias from thermodynamic fluctuations, we seek to relate the fluctuations to the nonthermal pressure support. Given that we expect the nonthermal pressure support to be predominantly from turbulent motions, we can relate the Mach numbers, which we infer from the spectra, to the nonthermal pressure support as:

$$ \frac{P_{\text{NT}}}{P_{\text{th}}} = (\gamma/3) M_{\text{3D}}^2, \qquad (3) $$

where P<sub>th</sub> is the thermal pressure, γ is the adiabatic index, and M<sub>3D</sub> is the Mach number of the turbulent velocity in 3D (e.g., Lau et al. 2009; Khatri & Gaspari 2016). When calculating the canonical hydrostatic mass one needs the slope of the thermal pressure (with respect to radius); similarly, one needs the slope of the nonthermal pressure to infer a hydrostatic mass bias. In practice we calculate the logarithmic slope of the turbulent Mach number (see Section 6) in the calculation of a hydrostatic mass bias. As such, we must determine Mach numbers which can be associated with at least two distinct radii.

In this study, we build on our previous theoretical and observational investigations (e.g., Gaspari & Churazov 2013; Gaspari et al. 2014; Khatri & Gaspari 2016; Romero et al. 2023) and seek to leverage key scaling relations between the relative thermodynamic fluctuations and the turbulent

kinematics in the ICM via the Fourier power spectra (e.g., Mach numbers, injection scales, and spectral slopes) when applied to a cluster sample. In particular, we seek to apply these scaling relations to both SZ and X-ray data where both data sets have overlap in angular (physical) scales accessed.

We discuss our full and pilot sample selection in Section 2. Together with Gaspari et al. (2014), Khatri & Gaspari (2016), and Romero et al. (2023), this work serves as our last pilot SZ/ X-ray spectral study to complete the foundational methodology to be applied to the full sample. As investigations of surface brightness fluctuations for the ICM require the relative fluctuations, i.e., δI/Î, or the residual surface brightness divided by some model, we discuss how we obtain Î in Sections 3 and 4 for XMM-Newton and the South Pole Telescope (SPT), respectively. Section 5 introduces the δI/Î images and discusses the methodology to obtain 2D amplitude spectra and subsequent deprojection to obtain 3D amplitude spectra. Results and inferences are presented in Section 6 and we reflect on the impact of methodological choices on our results in Section 7. We provide concluding remarks in Section 8.

We adopt a concordance cosmology with H<sub>0</sub> = 70 km s<sup>-1</sup> Mpc<sup>-1</sup>, Ω<sub>M</sub> = 0.3, and Ω<sub>Λ</sub> = 0.7. Unless otherwise stated all uncertainties are reported as the standard deviation (for distributions taken to be symmetric) or the distance from the median to the 16th and 84th percentiles (when allowing for asymmetric distributions).

## 2. Sample Selection

We seek to analyze a sample of clusters across a wide redshift range and which have been observed both via the SZ effect and in the X-rays. The SPT-SZ Survey (Bleem et al. 2015) identified 516 galaxy clusters in 2500 deg<sup>2</sup> which span redshifts 0 < z < 1.8 and with masses of M<sub>500</sub> &ge; 3 &times; 10<sup>14</sup> M<sub>☉</sub>. We opt to select SPT clusters which have been observed with XMM-Newton as photon counts are likely to be a limiting factor. Such a selection has already been studied with respect to thermodynamic profiles as determined from X-ray data in Bulbul et al. (2019) and thus provides an excellent sample which complements our study of fluctuations. Our sample selection is thus that used in Bulbul et al. (2019), though additional archival XMM-Newton data will augment the sample size from the 59 used in Bulbul et al. (2019).

Of this sample, we choose two clusters which are well resolved by SPT to pilot our analysis of surface brightness fluctuations. We visually inspected a subsample of clusters for which R<sub>500</sub> is several times the FWHM of the SPT beam (FWHM = 1.25). We wish to have a dynamically relaxed and a disturbed cluster, which we crudely estimate from the SPT Compton γ images. Based on this visual inspection, we choose SPT-CLJ0232-4421 as a dynamically relaxed cluster and SPT-CLJ0638-5358 as our dynamically disturbed cluster. Figure 1 presents adaptively smoothed X-ray images with SPT contours overlaid. We list M<sub>500</sub> and corresponding R<sub>500</sub> as reported in the SPT-SZ catalog (Bleem et al. 2015) in Table 1; the angular scale of R<sub>500</sub> is computed with our assumed cosmology.

In SPT-CLJ0232-4421 Parekh et al. (2021) find substructure about 1' southwest of the cluster core using XMM-Newton and Chandra data and posit that while the cool core remains relatively undisturbed (and which Hudson et al. 2010 classify as a weak cool core) this substructure could be related to the detected radio halo (Kale et al. 2019). More recently MeerKAT has revealed candidate relics (Kale et al. 2022): an eastern relic
<figure><img src="image_2.png" /><figcaption>**Figure 1.** X-ray (XMM) images of SPT-CLJ0232-4421 (left) and SPT-CLJ0638-5358 (right) with SZ (SPT) contours in orange. The X-ray images are adaptively smoothed; the SZ contours begin at $3\sigma$ with $2\sigma$ intervals. The physical scale of $R_{500}$ is captioned above a (pink) bar showing the angular scale of $R_{500}$. The dashed cyan circles enclose $R_{500}$ and the dashed magenta circles separate the inner and outer regions (rings; see Section 5). The red curve in the image of SPT-CLJ0638-5358 indicates the shock location as found in Botteon et al. (2018).</figcaption></figure>

<table><caption>**Table 1**
Fundamental Properties of the Clusters in Our Pilot Sample</caption><thead><tr><th>Cluster</th><th>Alternative Names</th><th>z<sup>a</sup></th><th>M<sub>500</sub><sup>a</sup><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>R<sub>500</sub><sup>a</sup><br>(Mpc)</th><th>θ<sub>500</sub><br>(arcmin)</th><th>kT<sub>ce</sub><sup>a</sup><br>(keV)</th><th>kT<sub>ce</sub><sup>b</sup><br>(keV)</th></tr></thead><tbody><tr><td rowspan="2">SPT-CLJ0232-4421</td><td>RXCJ0232.2-4420</td><td rowspan="2">0.28</td><td rowspan="2">9.45</td><td rowspan="2">1.36</td><td rowspan="2">5.32</td><td rowspan="2">7.19<sup>+0.46</sup><sub>-0.50</sub></td><td rowspan="2">10.06 ± 2.31</td></tr><tr><td>SPT-CLJ0232-4420</td></tr><tr><td>SPT-CLJ0638-5358</td><td>Abell S0592</td><td>0.23</td><td>9.42</td><td>1.38</td><td>6.26</td><td>8.44<sup>+0.86</sup><sub>-0.48</sub></td><td>9.45 ± 0.95</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Values are from Bulbul et al. (2019), who use XMM-Newton data.

<sup>b</sup> Values are from Mantz et al. (2010), who use Chandra data.

which lies within $R_{500}$ and a southern relic which lies at 1.9 Mpc (~$R_{200}$) from the cluster center. However, Kale et al. (2022) suspect that their origin is not from a cluster merger.

SPT-CLJ0638-5358, also listed as RXC J0638.7-5358 and Abell-S 592, had been noted as having a high ICM temperature, mass $M_{500} \approx 10^{15} M_{☉}$, and a noncircular shape in data from the Atacama Cosmology Telescope (e.g., Hincks et al. 2010; Mantz et al. 2010; Menanteau et al. 2010). More recently Botteon et al. (2018) have found a shock to the southwest of the cluster core (indicated with a red curve in Figure 1) using Chandra data. They derive a Mach number of $M = 1.72^{+0.15}_{-0.12}$ from the surface brightness (i.e., density jump) and they derive a consistent Mach value ($M = 1.6^{+0.54}_{-0.42}$) from the temperature jump. Despite searching for a counterpart shock to the northeast, they do not find such a counterpart. Botteon et al. (2018) also identify two cool cores (visible in their temperature map of SPT-CLJ0638-5358) and note the low entropies in those cool cores despite the clear dynamical disturbance present in the cluster.2

### 3. XMM-Newton Data Analysis

There are two XMM-Newton observations (ObsIDs) of SPT-CLJ0232-4421: 0042340301 and 0827350201; these are

<table><caption>**Table 2**
Overview of Imaging XMM-Newton Observations of SPT-CLJ0232-4421 and SPT-CLJ0638-5358</caption><thead><tr><th></th><th>SPT-CLJ0638-5358</th><th colspan="2">SPT-CLJ0232-4421</th></tr><tr><th></th><th>0650860101</th><th>0042340301</th><th>0827350201</th></tr></thead><tbody><tr><th scope="row">Cluster:</th><td></td><td></td><td></td></tr><tr><th scope="row">ObsID</th><td></td><td></td><td></td></tr><tr><th scope="row">Date</th><td>2010 May 22</td><td>2002 Jul 11</td><td>2021 Jan 1</td></tr><tr><th scope="row">Exposure (ks)</th><td>46.8</td><td>13.4</td><td>27.8</td></tr><tr><th scope="row">Clean exposure (ks)</th><td>MOS1: 24.6<br>MOS2: 31.7<br>PN: 7.6</td><td>MOS1: 11.6<br>MOS2: 12.1<br>PN: 6.8</td><td>MOS1: 24.6<br>MOS2: 25.7<br>PN: 18.0</td></tr><tr><th scope="row">PN mode</th><td>FF</td><td>EFF</td><td>FF</td></tr><tr><th scope="row">PI</th><td>J. Hughes</td><td>H. Böhringer</td><td>M. Arnaud & S. Ettori</td></tr></tbody></table>

**Note.** FF refers to the “full frame” mode.

summarized in Table 2. There is a single XMM-Newton observation of SPT-CLJ0638-5358 with ObsID of 0650860101. Our data processing and analysis are nearly identical to those used in Romero et al. (2023) where we use
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Surface brightness profiles for SPT-CLJ0232-4421 (left) and SPT-CLJ0638-5358 (right) and the respective single β-model fits and residuals. For brevity, the profiles presented here are the combination of profiles from each CCD and each ObsID. In practice the fits are done to each individual CCD and ObsID. Residuals are shown both in terms of χ and perhaps more relevant to the fluctuations, as δS/δ̃. The vertical red line marks R$_{500}$.</figcaption></figure>

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

Heasoft v6.28 and SAS 19.0 and the Extended Source Analysis Software (ESAS) data reduction package (Snowden et al. 2008) to produce event files and eventually images for the three EPIC detectors: MOS1, MOS2, and PN.

### 3.1. Image Creation

We choose to extract images in the [0.4–1.25] keV and [2.0–5.0] keV bands thus avoiding the energies where fluorescent lines are present in the EPIC cameras. Images and vignetted exposures are extracted for each detector over the entire detector area while masking point sources (see Section 3.2 for point-source identification) via the tasks `mos-spectra` or `pn-spectra`.

### 3.2. Point-source Exclusion

We generate a list of point sources using `cheese` and initially mask the sources with circles of radius 25". If the signal-to-noise ratio (S/N) of the point source is above 50 we modify the masking radius to 35". Finally, we perform a manual inspection to identify any remaining point sources or any necessary modifications to the masking circles. In the case of SPT-CLJ0232-4421 a very bright source at (R.A., decl.) = (38.°158650, -44.°363767) requires masking out to a radius of 80".

### 3.3. Profile Fitting

We use the Python package `pyproffit` (Eckert et al. 2017) to determine centers, ellipticity, and extract profiles for each of our images. As we have low-energy (LE; [0.4–1.25] keV) and high-energy (HE; [2.0–5.0] keV) images, per camera and per ObsID, we have 12 images in total for SPT-CLJ0232-4421 and six images in total for SPT-CLJ0638-5358. We adopt the average centroid (across all respective images) as our X-ray center, where `pyproffit` accounts for masked regions, point sources, and chip gaps, where we consider the gaps to include low-exposure pixels (Romero et al. 2023).

The best-fit profile parameters and uncertainties are determined via `emcee` (Foreman-Mackey et al. 2013) separately for each detector, each energy band, and each

observation. We adopt a circular β-model with the form:

$$S(r) = S_0[(1 + (r/r_c)^2)^{-3\beta+0.5} + B, \quad (4)$$

where r is the radius, $r_c$ is the “core” (scaling) radius, $S_0$ is the surface brightness normalization, and B is the background. We find that a uniform background component is sufficient and that to constrain the background we should fit (from r = 0) out to at least 10'. The β-model leaves evident residuals in Figure 2, where we have combined the profiles of each CCD across ObsIDs for brevity; the parameter results of these abbreviated profiles are also listed in Table 4. To reiterate, we fit the β-models to each CCD and each ObsID separately. In all cases we find that aside from the innermost radii the β-model does quite well. More specifically, when considering the context of developing a method to be applied to a sample of clusters with varying depth, we wish to adopt a model which should do well across the sample. Finally, the actual impact of the residuals in the center is likely small given the area it subtends relative to the regions used for power spectra extraction (see Sections 5 and 7). We also explore the impact of ellipticity in Section 7.1.

## 4. SPT Image Analysis

The SPT-SZ (Bleem et al. 2015) products are available in spherical coordinates (HEALpix format) or a flat sky (Sanson-Flamsteed projection).<fnref n="17" /> As the angular size of $R_{500}$ for our clusters is less than 10', we expect that the flat-sky projection is adequate and allows for the application of the same analysis tools for the X-ray and SZ data. There are also several flavors of Compton y maps available: minimum variance, CMB-Nulled, and CMB-CIB-Nulled (Bleem et al. 2022). When confined to the vicinity of a cluster, the cosmic microwave background (CMB) and cosmic infrared background (CIB) are subdominant on the scales probed and so we opt for the minimum variance maps, which provide the best resolution and signal significance.

<footnote n="17">¹⁷ https://lambda.gsfc.nasa.gov/product/spt/spt_sz_comp Maps_get.html</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The power spectra of the SPT-SZ fields for the minimum variance Compton y maps. Comments on the departure(s) from a flat spectrum are given in the text. The fields ra3h30dec-42.5 and ra6h30dec-55 contain SPT-CLJ0232-4421 and SPT-CLJ0638-5358, respectively. The dotted blue line is at 0.8 arcmin⁻¹, corresponding to the reported FWHM of 1/25. The black dashed curve shows an arbitrarily normalized power spectrum for a Gaussian with a FWHM of 1/25.</figcaption></figure>

## 4.1. Characterizing SPT Noise

The SPT Compton y maps are provided in fields (over 10° on a side) of 19 contiguous subpatches (e.g., Bleem et al. 2015). The noise in each field is expected to be homogeneous. To characterize the noise (on scales of 20' and smaller), we identify 100 regions within a field that are devoid of point sources or any other evident signal. This is achieved by avoiding regions which are masked.<fnref n="18" />

The noise in the maps is assumed to be stationary; as such, the power spectrum of the 100 blank regions identified above serve to characterize the noise. To ensure further that we avoid signal in our noise characterization, we calculate the power spectrum on the blank regions of difference map (half 1–half 2). In particular, we average over the 100 regions to calculate the average power spectrum shown in Figure 3. We identify three features in the noise power spectrum that warrant comment: (1) at $k \approx 10^{-1}$ arcmin⁻¹ the steepening of the slope is due to the transition of weights between SPT and Planck (with Planck’s weight increasing toward lower values of $k$; Bleem et al. 2022), where the weight is especially dominated by the 217 GHz band in Planck data; (2) just below $k = 1$ arcmin⁻¹ the power drops quickly (toward higher $k$); and (3) a similar drop occurs again around $k \approx 2$ arcmin⁻¹. These are features of the projection from HEALpix to flat sky (Sanson–Flamsteed projection) that are inconsequential in our analyses.

## 4.2. β-profile Fitting

To fit the cluster profile we extract cutouts out to $2 R_{500}$, i.e., the sides are 4 $R_{500}$ (rounding up for selecting pixels). Although no point sources are found within $R_{500}$ for the two clusters investigated, we do see a point source beyond $R_{500}$ in the cutout of SPT-CLJ0232-4421. We thus pass the cutout images through the Daophot utility within astropy. A point-source map is then constructed from the locations and peaks as

found with the Daophot utility. We assume a circular Gaussian with a FWHM of 1/25 for each point source and subtract any sources found from the maps.

To fit a surface brightness profile to our point-source-subtracted maps, we employ emcee (Foreman-Mackey et al. 2013) and forward model a circular β-model of the cluster. The forward modeling includes the convolution with the beam (circular Gaussian with a FWHM of 1/25) and transfer function (provided as a function of $ℓ^{19}$). The impact of the point source (if not subtracted) on the profile fit is negligible and so we are not concerned with a more detailed treatment of point sources (e.g., such as simultaneous fitting with the cluster).

We perform the fitting in image space. We take our principal results to be those which assume the centroid from XMM-Newton data and a Fourier (power spectrum) characterization of the noise.<fnref n="20" /> We check both assumptions with fits where both the centroid is allowed to vary (see Table 3) and fits where we assume the noise is uniformly white, which we denote by rms (see Figure 4). The fitted parameters of both (abbreviated) X-ray and SZ β-models are presented in Table 4. The profiles do not change much with the change of centroid nor with noise characterization, though we do note a tendency for the uniform white noise to prefer a steeper slope near $R_{500}$.

Residual images of both clusters are shown in Figure 5 and we see that both clusters exhibit substructure. Additionally, Figure 5 shows the regions that are used in the spectral analysis (Section 5) and for SPT-CLJ0638-5358 we have also marked the location of a known shock (Botteon et al. 2018) in red and possible shocks as seen in the SPT residuals in dark green.

# 5. Power Spectrum Analysis

As indicated in Section 1, our first step in calculating spectra of surface brightness fluctuations is to obtain fractional residual images. Sections 3 and 4 have discussed how we obtain surface brightness images and fit and subtract a model to derive residual images. For each residual image, we need only divide by the respective fitted ICM surface brightness model (i.e., excluding a fitted background component if present); these fractional residuals are presented in Figure 6.

To calculate power spectra we employ the Δ-variance method presented in Arévalo et al. (2012), which is designed to provide low-resolution spectra for arbitrary data masking. As noted in Arévalo et al. (2012), Zhuravleva et al. (2015), and Romero et al. (2023), there are biases associated with this method. To the extent that one knows the underlying spectra (or can at least make reasonable estimates of them), these biases can be corrected. We consider power spectra of the form $P(k) \propto e^{-k_c/k}k^{-\alpha}$, with a low cutoff wavenumber $k_c = 1/(5 R_{500})$ such that we effectively consider a single power law within our sampled range. In particular, we correct for the bias when $\alpha = 3$ and we consider the spread of biases over the range ($2 < \alpha < 4$) as a source of systematic uncertainty that we fold into our reported uncertainties as in Romero et al. (2023). We privilege the bias for $\alpha = 3$ because this will be the slope of the power spectrum at the peak of the 3D amplitude spectrum (by definition; see Equation (6)) and Section 6 clarifies our interest in the peak of the amplitude power spectrum). Appendix C investigates biases associated with large scales due to masking.

<footnote n="18">¹⁸ Masks are provided at https://lambda.gsfc.nasa.gov/data/suborbital/SPT/bleem_2021/masks.zip.</footnote>

<footnote n="19">¹⁹ https://lambda.gsfc.nasa.gov/data/suborbital/SPT/bleem_2021/sptsz_trough_filter_1d.txt</footnote>

<footnote n="20">²⁰ Again, this choice comes from the expectation that this will provide good consistency across our full sample.</footnote>
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Profiles (data and models) for SPT-CLJ0232-4421 (left) and SPT-CLJ0638-5358 (right). The model profile includes the effect of beam convolution. In the case of SPT-CLJ0232-4421, the positive Compton y values well beyond R<sub>500</sub> are due to primary CMB fluctuations as this feature is not present in the CMB-nulled map.</figcaption></figure>

<table><caption>**Table 3**
Centers of the Two Clusters Taken to Be the Centroid as Found with *pyproffit* (XMM<sub>p</sub>), Taken to Be the Peak as Found with *pyproffit* (XMM<sub>r</sub>), as Indicated in the SPT-SZ Catalog (SPT), and as Found When Allowing the Center to Be a Parameter in the β-model Fit (SZ)</caption><thead><tr><th rowspan="2"></th><th colspan="3">SPT-CLJ0232-4432</th><th colspan="3">SPT-CLJ0638-5531</th></tr><tr><th>R.A. (deg)</th><th>Decl. (deg)</th><th>Δ(′′)</th><th>R.A. (deg)</th><th>Decl. (deg)</th><th>Δ(′′)</th></tr></thead><tbody><tr><td>XMM<sub>r</sub></td><td>38.072751</td><td>-44.348505</td><td>...</td><td>99.699212</td><td>-53.976286</td><td>...</td></tr><tr><td>XMM<sub>p</sub></td><td>38.078971</td><td>-44.346976</td><td>16.9</td><td>99.702212</td><td>-53.974323</td><td>9.5</td></tr><tr><td>SPT</td><td>38.070109</td><td>-44.354129</td><td>21.4</td><td>99.697766</td><td>-53.974866</td><td>6.0</td></tr><tr><td>SZ</td><td>38.069225</td><td>-44.349572</td><td>9.4</td><td>99.697749</td><td>-53.976731</td><td>3.5</td></tr></tbody></table>

<table><caption>**Table 4**
β-model Fit Parameters for the SZ (SPT) and X-Ray (XMM) Data Sets, Where I<sub>0</sub> Is in Units of Compton y for the SZ fits and in Units of photon counts s<sup>-1</sup> arcmin<sup>-2</sup> for the X-Ray Fits. In the SZ case, the Two Subsets (rms and Fourier) Refer to the Same Data but Two Different Characterizations of Noise When Fitting</caption><thead><tr><th>Cluster</th><th>Data Set</th><th>Subset</th><th>I<sub>0</sub></th><th>θ<sub>c</sub></th><th>β</th></tr></thead><tbody><tr><td rowspan="4">SPT-CLJ0232-4432</td><td rowspan="2">SPT</td><td>rms</td><td>1.34<sup>+0.04</sup><sub>-0.04</sub> × 10<sup>-4</sup></td><td>2.49<sup>+0.19</sup><sub>-0.17</sub></td><td>1.11<sup>+0.13</sup><sub>-0.10</sub></td></tr><tr><td>Fourier</td><td>1.39<sup>+0.02</sup><sub>-0.02</sub> × 10<sup>-4</sup></td><td>1.97<sup>+0.09</sup><sub>-0.09</sub></td><td>0.76<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td rowspan="2">XMM</td><td>0.4–1.25 keV</td><td>1.66<sup>+0.03</sup><sub>-0.03</sub></td><td>0.38<sup>+0.01</sup><sub>-0.01</sub></td><td>0.560<sup>+0.003</sup><sub>-0.003</sub></td></tr><tr><td>2–5 keV</td><td>0.87<sup>+0.01</sup><sub>-0.01</sub></td><td>0.42<sup>+0.01</sup><sub>-0.01</sub></td><td>0.572<sup>+0.003</sup><sub>-0.003</sub></td></tr><tr><td rowspan="4">SPT-CLJ0638-5358</td><td rowspan="2">SPT</td><td>rms</td><td>1.69<sup>+0.16</sup><sub>-0.15</sub> × 10<sup>-4</sup></td><td>1.29<sup>+0.21</sup><sub>-0.18</sub></td><td>0.44<sup>+0.08</sup><sub>-0.07</sub></td></tr><tr><td>Fourier</td><td>1.73<sup>+0.04</sup><sub>-0.04</sub> × 10<sup>-4</sup></td><td>1.21<sup>+0.05</sup><sub>-0.05</sub></td><td>0.41<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td rowspan="2">XMM</td><td>0.4–1.25 keV</td><td>1.21<sup>+0.03</sup><sub>-0.02</sub></td><td>0.68<sup>+0.02</sup><sub>-0.02</sub></td><td>0.67<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>2–5 keV</td><td>0.76<sup>+0.01</sup><sub>-0.01</sub></td><td>0.73<sup>+0.02</sup><sub>-0.02</sub></td><td>0.69<sup>+0.01</sup><sub>-0.01</sub></td></tr></tbody></table>

Here again we apply corrections for the inferred biases for underlying power spectra of $\alpha=3$. Independent of the $\Delta$-variance method, our images are the convolution of each instrument's point-spread function (PSF), which we correct for by dividing by the normalized power spectrum of the appropriate PSF (where we consider the CCD position and use an *ellbeta* PSF model to determine this for the XMM data).

We wish to sample scales between (roughly) the FWHM of the PSF of each instrument and $R_{500}$ and we space apart the angular nodes where the power spectra are sampled by roughly a factor of 2. These constraints lead us to sample the SPT data

over four nodes and the XMM data over seven nodes, of which we keep the four smallest nodes (largest scales) to match those used for the SPT data. As mentioned in Section 1, if we hope to constrain the hydrostatic mass bias we will need an estimate of the radial slope of the rms turbulent velocity (Mach number) and thus we desire at least two regions: an inner circle and outer annulus (or annuli). If we are to sample scales of $R_{500}$, then we expect that our inner circle should have a radius of at least $r = R_{500}/2$. Conversely, we do not want any region (annulus) to be thinner than the largest beamwidth (taken as the FWHM) of the instruments and so we set a minimum width of annuli to 1/25. There are additional considerations with respect
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left: the residual (Compton y) image of SPT-CLJ0232-4421. Right: the residual (Compton y) image of SPT-CLJ0638-5358 with the known shock front location (Botteon et al. 2018) in red and plausible shock fronts from the SPT data in dark green. In both clusters, R<sub>500</sub> is shown as the black dashed circle. The division between Ring 1 and Ring 2 used in Section 5 is marked by the dashed magenta circle. The rms within the fields for both clusters is roughly 1.8 × 10<sup>-5</sup> in Compton y.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Left: fractional residuals for SPT-CLJ0232-4421. Right: fractional residuals for SPT-CLJ0638-5358. The top panels are from XMM-Newton data (δS/S), specifically images from the PN CCD for our LE band (0.4–1.25 keV for a single ObsID). The bottom panels are from SPT data (δy/y). The known and candidate shock fronts in SPT-CLJ0638-5358 presented in Figure 5 are marked in green on the X-ray fractional residual.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Left: the 2D amplitude spectra for SPT-CLJ0232-4421 color coded by energy band (legend reports in units of eV). The solid-color filled regions show the 1σ range for Ring 1 and the hatched and stippled regions show the 1σ range for Ring 2. Right: the 2D amplitude spectra for SPT-CLJ0638-5358 with the same scheme for the energy bands and rings.</figcaption></figure>

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

to the division of regions that are discussed and investigated in Section 7, which still do not offer a clear choice. Left with some ambiguity, we opt to divide our rings at $R_{500}/\phi$, where we take $\phi = (1 + \sqrt{5})/2$, the golden ratio. We thus define Ring 1 to be a circle of radius $3.3$ for SPT-CLJ0232-4421 and a circle of radius $3.9$ for SPT-CLJ0638-5358 and Ring 2 to be the annulus from the edge of Ring 1 to each cluster's respective $R_{500}$ (5.3 and 6.3).

As elsewhere, we report our results in the form of amplitude spectra; given 2D ($P_{2D}$) and 3D ($P_{3D}$) power spectra, amplitude spectra are given as:

$$A_{2D}(k) = [(2\pi)k^2 P_{2D}]^{1/2}, \qquad (5)$$

$$A_{3D}(k) = [(4\pi)k^3 P_{3D}]^{1/2}, \qquad (6)$$

where $k = 1/\lambda$ is the wavenumber which corresponds to (angular or physical) scales of length $\lambda$.

## 5.1. Application to XMM Data

As with our profile fits, power spectra are measured on each fractional image (i.e., per CCD, per ObsID, and per energy band) and mask the same pixels as masked in the profile fits. We do not measure the spectra on a single image (e.g., as presented in Figure 6), but for each CCD, energy band, and ObsID, we take 100 surface brightness models from the respective Markov Chain Monte Carlo chain (of the surface brightness profiles, long after burn-in). For each set of model parameters, a single Poisson noise realization is generated, from which we can create a noise realization of the fractional residual in addition to the (data) fractional residual for that model. The spectra of these two images are measured and their difference is our debiased spectrum. The debiased spectra also have the same treatment to account for faint point sources as in Romero et al. (2023).

We take the expected values and associated uncertainties to be the mean and standard deviation, respectively, of the debiased spectra for each CCD and ObsID. We combine the spectra for each energy band across CCDs and ObsIDs as their weighted averages. Figure 7 shows mild (more so visually than

statistically) differences in the amplitude spectra of each energy band for both clusters and both rings.

## 5.2. Application to SPT Data

We take the power spectrum of each SPT cluster to be that derived from the fractional residuals of its best-fit profile. In order to debias the spectrum, we make fractional residuals on each cluster's respective half maps and take the cross spectrum of the two half maps for each cluster. The uncertainty is taken to be the standard deviation of the cross spectra of the half-map blank regions (from the respective fields). In case of any leakage of instrumental noise, we subtract the mean across these noise cross spectra (the effect of this is negligible). Spectra are not reported for Ring 2 as none are statistically significant (Figure 8).

## 5.3. Deprojection to 3D

We wish to deproject our 2D fluctuation spectra to 3D fluctuations. To do this we use a common formalism (e.g., Peacock 1999; Churazov et al. 2012; Khatri & Gaspari 2016):

$$P_{2D}(k_{\theta}) = \int P_{3D}(\sqrt{k_{\theta}^2 + k_z^2}) |\tilde{W}(k_z)|^2 dk_z, \qquad (7)$$

where $z$ is the axis along the line of sight, $k_{\theta}^2 = k_x^2 + k_y^2$ is in the plane of the sky, and $|\tilde{W}(k_z)|^2$ is the 1D power spectrum of the window function, which normalizes the distribution of the relevant (unperturbed) 3D signal generation to the (unperturbed) 2D surface brightness. In the case of Compton y, the units of our SZ surface brightness maps, we already have Equation (1). For our X-ray surface brightness, $S$, maps in units of photon counts $s^{-1} \text{arcmin}^{-2}$, we can write:

$$S = \int \epsilon dz, \qquad (8)$$

where $\epsilon$ encapsulates the EI and band-averaged cooling function, $\Lambda_b$, presented in Section 1. We denote the smooth (i.e., model) 3D distributions of $P_e$ and $\epsilon$ as $\bar{P}_e$ and $\bar{\epsilon}$, respectively, which when integrated along the line of sight, produce $\bar{y}$ and $\bar{S}$, the 2D (circular, unperturbed) surface brightness models. The SZ and
<figure><img src="image_12.png" /><figcaption>Figure 8. Amplitude spectrum of δy/y for Ring 1 in SPT-CLJ0232-4421 (left) and SPT-CLJ0638-5358 (right).</figcaption></figure>

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

X-ray window functions are, respectively:

$$W_{SZ}(\theta, z) \equiv \frac{\sigma_T}{m_e c^2} \frac{\bar{P}_e(\theta, z)}{\bar{y}(\theta)}, \qquad (9)$$

$$W_X(\theta, z) \equiv \frac{\bar{\bar{\epsilon}}(\theta, z)}{\bar{S}(\theta)}. \qquad (10)$$

Equation (7) can be approximated as:

$$P_{2D}(k_{\theta}) \approx P_{3D}(k) \int |\tilde{W}(k_z)|^2 dk_z, \qquad (11)$$

(e.g., Churazov et al. 2012; Khatri & Gaspari 2016), where $k = \sqrt{k_{\theta}^2 + k_z^2}$. We adopt the notation used in Khatri & Gaspari (2016), Romero et al. (2023) and define:

$$N(\theta) \equiv \int |\tilde{W}(k_z)|^2 dk_z. \qquad (12)$$

This approximation is quite robust (e.g., Khatri & Gaspari 2016; Romero et al. 2023) for most scales probed; we verify the validity of this approximation for our clusters in Appendix B. Another concern with deprojection is the evolution of the window function within a region (Zhuravleva et al. 2015), which we consider in Appendix C where the result is that we weight a sampling of $N(\theta)$ by the annular area to derive an effective value for a given ring, $N_{eff}$.

Once we have $N_{eff}$, from the best-fit profiles of each image (again, for XMM data this is considered per CCD, per energy band, and per ObsID), we calculate the weighted average of deprojected spectra from individual deprojected spectra as before with the 2D spectra. We further combine both energy bands considering the overall agreement. If we simply divide $P_{2D,X-ray}$ for a given ring by the respective $N_{eff}$, then we have only calculated the power, $P_{3D}$ of the underlying emmisivity fluctuations. Dividing by a further factor of 4 (e.g., Churazov et al. 2012) yields the density fluctuations.<fnref n="21" /> The 3D spectra ($P_{3D}$) for SPT deprojected from their 2D counterparts are shown in Figure 8 and directly yield pressure fluctuations. We present both SZ and X-ray 3D amplitude spectra in Figure 9 with dotted gray lines indicating a Kolmogorov spectrum

($\alpha = 11/3$ for the convention $P_{3D}(k) \propto k^{-\alpha}$; thus the plotted amplitude spectra have a logarithmic slope of $-1/3$).

## 6. Results and Inferences from the Spectra

We are able to constrain confidently the spectra for both clusters at the largest scales we probe, and we do capture the peaks of the amplitude spectra. We recover injection scales ($\lambda_{peak}$) of several hundreds of kiloparsecs (see Table 5), which is largely within expectations (see Gaspari et al. 2014; Khatri & Gaspari 2016; Eckert et al. 2017). The values of the amplitude spectra at scales smaller than the inferred injection scale are generally smaller, with the exception of the amplitude spectrum of SPT-CLJ0232-4421 at $k = 1.6 \text{ arcmin}^{-1}$. Unfortunately the uncertainties at scales smaller than our inferred injection scales are often large and thus we cannot exclude with high confidence the possibility of slightly smaller injection scales. Similarly, in both clusters and both rings, our uncertainties at scales of $1'$ or smaller do not allow us to constrain the slope of any turbulent cascade. Even so, we have drawn a guide line with an arbitrary normalization and a spectral slope of $\alpha = 11/3$ (logarithmic slope of $-1/3$ for the 3D amplitude spectra), i.e., a Kolmogorov spectrum, in Figure 9.

There are two principle avenues to derive turbulent velocities from power spectra: use either the peak of the amplitude spectrum (e.g., Gaspari & Churazov 2013; Khatri & Gaspari 2016; Eckert et al. 2017) or some measure of the dispersion of fluctuations in resolution elements of numerical simulations (e.g., Zhuravleva et al. 2012; Mohapatra et al. 2021; Simonte et al. 2022). The (total) variance can be calculated from the integrated power spectrum:

$$\sigma_{3D}^2 = \int P_{3D}(k) 4\pi k^2 dk. \qquad (13)$$

We report the linear standard deviation ($\sigma_{3D}$) and equivalent (natural) logarithmic standard deviation ($\sigma_{ln}$) in Table 5 along with the peak of the amplitude spectra of density fluctuations. Though the linear and logarithmic standard deviations are nearly identical, we bring our distribution values into compliance with the width in Zhuravleva et al. (2023), which is equivalent of a FWHM of a $\log_{10}$-normal distribution ($\delta\xi/\xi$

<footnote n="21">²¹ This holds especially well for low photon energies in hot, massive clusters where the photon counts depend negligibly on gas temperature.</footnote>
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 9.** Amplitude spectra of $\delta P/P$ and $\delta\rho/\rho$ for SPT-CLJ0232-4421 (left) and SPT-CLJ0638-5358 (right). The dotted lines have a spectral index, $\alpha = 11/3$, which equates to a slope of $-1/3$ in the amplitude spectra. Arrows indicate $3\sigma$ upper limits. The constraints on the amplitude spectra are not adequate to constrain the cascade (slope) of fluctuations at $k > k_{\text{peak}}$. The roughly flat spectra (plateau) of Ring 1 for both clusters is likely indicative of multiple injection scales (see Section 6).</figcaption></figure>

<table><caption>**Table 5**
Key Values from the Amplitude and Power Spectra</caption><thead><tr><th>Cluster</th><th>Data Set</th><th>Ring</th><th>A<sub>3D</sub>(k<sub>peak</sub>)</th><th>σ<sub>3D</sub></th><th>σ<sub>ln</sub></th><th>k<sub>peak</sub> (arcmin<sup>-1</sup>)</th><th>λ<sub>peak</sub> (kpc)</th></tr></thead><tbody><tr><td rowspan="2">SPT-CLJ0232-4421</td><td>X-ray</td><td>Ring 1</td><td>0.16 ± 0.03</td><td>0.21 ± 0.17</td><td>0.20 ± 0.11</td><td>0.49</td><td>516</td></tr><tr><td>X-ray</td><td>Ring 2</td><td>0.14 ± 0.06</td><td>0.11 ± 0.06</td><td>0.11 ± 0.05</td><td>0.30</td><td>837</td></tr><tr><td rowspan="3">SPT-CLJ0638-5358</td><td>X-ray</td><td>Ring 1</td><td>0.17 ± 0.05</td><td>0.17 ± 0.04</td><td>0.17 ± 0.04</td><td>0.27</td><td>807</td></tr><tr><td>SZ</td><td>Ring 1</td><td>0.19 ± 0.10</td><td>0.15 ± 0.88</td><td>0.15 ± 0.31</td><td>0.16</td><td>1380</td></tr><tr><td>X-ray</td><td>Ring 2</td><td>0.44 ± 0.20</td><td>0.28 ± 0.11</td><td>0.27 ± 0.11</td><td>0.27</td><td>807</td></tr></tbody></table>

in their notation) and note that this width is equal to $1.02 \sigma_{\text{ln}}$. The integrated spectra are calculated up until $P_{3D} < 0$.

Unfortunately much of the recovered pressure fluctuations are not statistically significant; however, the largest scale probed in SPT-CLJ0638-5358 is significant and so we report spectral properties of pressure fluctuations for SPT-CLJ0638-5358 in Table 5. In all cases we require that any selected peak be at least $2\sigma$ significance, that is, we select the largest amplitude that meets this criterion.

To derive Mach numbers from the peak amplitude spectra we use as reference the linear relation first found via high-resolution hydrodynamical simulations by Gaspari & Churazov (2013):

$$ \begin{aligned} \mathcal{M}_{3D,\text{peak}} &\approx 4 A_{\rho}(k_{\text{peak}}) \left( \frac{l_{\text{inj}}}{500 \text{ kpc}} \right)^{\alpha_H} \\ &\approx 2.4 A_{P}(k_{\text{peak}}) \left( \frac{l_{\text{inj}}}{500 \text{ kpc}} \right)^{\alpha_H}, \end{aligned} \quad (14) $$

where $l_{\text{inj}} = 1/k_{\text{peak}}$ (here in units of kpc) and $0.2 < \alpha_H < 0.3$ is the slight dependence on injection scale. We also test a relation from large-scale cosmological simulations (with 80 clusters; Zhuravleva et al. 2023); noting that $\mathcal{M}_{3D} = \sqrt{3} \mathcal{M}_{1D}$, they find the following radially averaged ($r < 2 R_{500}$) scaling with the total standard deviation:

$$ \begin{aligned} \mathcal{M}_{3D,\text{int}} &\approx \sqrt{3} * \langle \eta_{\rho} \rangle * 1.02 * \sigma_{\text{ln},\rho} \\ &\approx \sqrt{3} * \langle \eta_{P} \rangle * 1.02 * \sigma_{\text{ln},P}, \end{aligned} \quad (15) $$

where $\eta_{\rho}$ and $\eta_{P}$ are functions of radius, cluster geometry (spherical or ellipsoidal), and dynamical state (relaxed, in-between, or unrelaxed). With respect to the choice of $\eta$, we thus utilize the value for the in-between state for SPT-CLJ0232-4421 and the unrelaxed state for SPT-CLJ0638-5358. For the spherical case, this gives the radially averaged values of $\langle \eta_{\rho} \rangle = 1.2$ for both the in-between state and unrelaxed state, $\langle \eta_{P} \rangle = 1.4$ for the in-between state, and $\langle \eta_{P} \rangle = 1.5$ for the unrelaxed state. Mach numbers corresponding to each estimation method are reported in Table 6; a combined value is also calculated as the average of the values derived from the two methods. The uncertainty on this combined value takes the spread (difference) between the two methods as a systematic uncertainty. Further quantities derived from $\mathcal{M}_{3D,\text{comb}}$ propagate the combined uncertainties (systematic and statistical uncertainties are summed).

By including nonthermal pressure support due to turbulence in the standard hydrostatic equilibrium (HE) equation and the relation of nonthermal pressure, thermal pressure, and Mach number in Equation (3) one can derive the hydrostatic bias in terms of the Mach number, logarithmic Mach slope, and logarithmic pressure slope. We briefly reiterate this here and use the notation and convention$^{22}$ from Khatri & Gaspari (2016) for the hydrostatic bias, $b_M \equiv M_x/M_{\text{tot}} - 1$, derived in

<footnote>$^{22}$ N.B. that the sign convention of $b_M$ is opposite that of the conventional hydrostatic bias, $b$.</footnote>
<table><caption>**Table 6**
Derived Mach numbers</caption><thead><tr><th>Cluster</th><th>Data Set</th><th>Ring</th><th>M<sub>3D,peak</sub></th><th>M<sub>3D,int</sub></th><th>M<sub>3D,comb</sub></th></tr></thead><tbody><tr><td rowspan="2">SPT-CLJ0232-4421</td><td>X-ray</td><td>Ring 1</td><td>0.63 ± 0.12</td><td>0.43 ± 0.23</td><td>0.53 ± 0.10 ± 0.11</td></tr><tr><td>X-ray</td><td>Ring 2</td><td>0.45 ± 0.18</td><td>0.24 ± 0.10</td><td>0.34 ± 0.11 ± 0.09</td></tr><tr><td rowspan="3">SPT-CLJ0638-5358</td><td>X-ray</td><td>Ring 1</td><td>0.60 ± 0.17</td><td>0.35 ± 0.08</td><td>0.48 ± 0.12 ± 0.07</td></tr><tr><td>SZ</td><td>Ring 1</td><td>0.59 ± 0.32</td><td>0.40 ± 0.83</td><td>0.50 ± 0.09 ± 0.30</td></tr><tr><td>X-ray</td><td>Ring 2</td><td>1.58 ± 0.72</td><td>0.57 ± 0.23</td><td>1.07 ± 0.51 ± 0.22</td></tr></tbody></table>

<table><caption>**Table 7**
Derived Hydrostatic Biases via two Mach Scaling Relations, as Well as a Third Estimation from the Combination of the Mach Estimations</caption><thead><tr><th></th><th></th><th>-b<sub>M,peak</sub></th><th>-b<sub>M,int</sub></th><th>-b<sub>M,comb</sub></th></tr></thead><tbody><tr><td rowspan="2">SPT-CLJ0232-4421</td><td>Ring 1</td><td>0.34 ± 0.20</td><td>0.26 ± 0.26</td><td>0.30 ± 0.27</td></tr><tr><td>Ring 2</td><td>0.13 ± 0.10</td><td>0.05 ± 0.04</td><td>0.08 ± 0.09</td></tr><tr><td rowspan="2">SPT-CLJ0638-5358</td><td>Ring 1</td><td>-0.42 ± 0.65</td><td>-0.01 ± 0.09</td><td>-0.16 ± 0.38</td></tr><tr><td>Ring 2</td><td>0.02 ± 0.73</td><td>0.08 ± 0.10</td><td>0.10 ± 0.44</td></tr></tbody></table>

this manner:

$$b_M = \frac{-\gamma M_{3D}^2}{3} \frac{d \ln P_{NT}}{d \ln P_{th}} \left( 1 + \frac{\gamma M_{3D}^2}{3} \frac{d \ln P_{NT}}{d \ln P_{th}} \right)^{-1}, \quad (16)$$

where $\gamma$ is the adiabatic index, taken to be $5/3$ for the ICM. We can reformulate $(d \ln P_{NT})/(d \ln P_{th})$ in terms of variables more directly calculated from our SZ surface brightness fits and derived Mach numbers (Table 6):

$$\frac{d \ln P_{NT}}{d \ln P_{th}} = \frac{d \ln P_{NT}/d \ln r}{d \ln P_{th}/d \ln r} = 1 + 2 \frac{d \ln M_{3D}/d \ln r}{d \ln P_{th}/d \ln r}. \quad (17)$$

Taking the $M_{3D,peak}$ values, we derive hydrostatic biases, which we report in Table 7; the logarithmic slopes of the Mach numbers and pressure profiles are reported in Appendix E. As we only have two rings, we can only estimate a single power law for $M_{3D}(r)$ across all radii of concern. For the logarithmic pressure slope we calculate the average slope within a ring based on our deprojected SZ $\beta$-model.

## 6.1. Interpretation

We selected our clusters in part based on a visual inspection of the SZ and X-ray images to assess dynamical states; we note that while there is correlation of the cluster morphology with the dynamics of the ICM (e.g., Battaglia et al. 2012; Zhuravleva et al. 2023), cluster morphology is not strictly indicative of the dynamics of the ICM. We note that turbulence itself does not require a clear visual disturbance (such as an elongation or surface brightness edge, as would emerge for a shock or cold front). Even turbulence driven inside-out via AGN feeding or feedback can occur in a quasi-spherical manner over several cycles, superposed on large-scale cosmological chaotic motions (Gaspari et al. 2013, 2020; Lau et al. 2017; Wittor & Gaspari 2020, 2023). Moreover, shocks or other injection mechanisms may exist but be masked via several observational effects, e.g., due to projection along the line of sight.

The Mach numbers for both clusters span $0.5 < M_{3D} < 1.6$, with both clusters having $M_{3D} \approx 0.6$ in Ring 1. It is perhaps surprising that the Mach numbers are similar in Ring 1 for the two clusters given that the known shock in SPT-

CLJ0638-5358 (Botteon et al. 2018) occurs within 2 arcmin of the cluster center, thus in Ring 1. However, neither the length of the detected X-ray edge associated with the shock front nor its radial extent are large. That is, the effect of the shock, which is already subdominant (see below and Appendix D), is dampened by the filling factor of its area relative to the area of the region (Ring 1). Even so, taking the X-ray Mach values, $M_{3D} = 0.64 \pm 0.03$ is not much larger than the derived value for SPT-CLJ0232-4421 in Ring 1 of $M_{3D} = 0.51 \pm 0.06$. Again, we come to the somewhat ambiguous state of SPT-CLJ0232-4421: many X-ray metrics indicate that it is relaxed (e.g., Lovisari et al. 2017), but the substructure found in Parekh et al. (2021) is evident in our data too and undoubtedly contribute to our derived Mach number. In Ring 2, both clusters show substructure (prominent residuals) in the SPT data, but distinct features are not readily visible in the X-ray (fractional) residuals. That said, the fluctuations appear to have a larger amplitude to the southwest and northeast of the cluster and an investigation of fluctuations separating directions in the cluster (see Appendix D) echoes this.

Our Mach numbers are generally larger than those found in other studies, e.g., Gaspari & Churazov (2013; $M_{3D} \approx 0.45$), Zhuravleva et al. (2015; $M_{3D} \approx 0.2$), Hitomi Collaboration et al. (2016; $M_{3D} \approx 0.3$), Hofmann et al. (2016; $M_{3D} \approx 0.3$), Eckert et al. (2017; $M_{3D} \approx 0.3$), Zhang et al. (2023; $M_{3D} \approx 0.2$), and Dupourqué et al. (2023; $M_{3D} \approx 0.2$). As noted in Hofmann et al. (2016), Mach numbers may be expected to grow with cluster-centric radius and with the scales being probed. That is, at larger cluster-centric radii we expect larger nonthermal pressure support (e.g., Battaglia et al. 2012) and thus larger Mach numbers. Similarly, the Mach number inferred at scales below the injection scale will be less than the Mach number inferred at the injection scale, which tends toward scales of a few hundred kiloparsecs (as most of the studies quoted above; e.g., Zhuravleva et al. 2015; Hofmann et al. 2016). Our analysis probes scales out to $0.62 R_{500} \approx 0.85$ Mpc for the inner ring and scales up to $R_{500}$, where we recover injection scales $l_{inj} \ge 500$ kpc; thus, on both accounts we might expect larger Mach numbers than found in many other studies. Even if one considers that due to poor spectral constraints, the injection scales may be smaller than inferred, one is still left with substantial fluctuations at our
inferred injection scales. While not many studies probe our large scale, a similar multiwavelength analysis by Khatri & Gaspari (2016) finds $M_{3D} \simeq 0.8$ at megaparsec scales. Another aspect when comparing Mach numbers across the literature is that there is a significant scatter in the Mach scaling relation (Gaspari & Churazov 2013; Gaspari et al. 2014; Zhuravleva et al. 2014; Mohapatra et al. 2021; Simonte et al. 2022; Zhuravleva et al. 2023). While nearly all of the literature results agree on a linear conversion between ICM fluctuations and turbulent velocities, the different numerical and reduction techniques can introduce differences in the normalization. However, the inferred Mach numbers as determined via the peak and via the integrated spectra are generally in agreement, though the inferred values from integrated spectra tend to be slightly less than those from the peak (Table 6). A similar trend was seen in Zwicky 3146 (Romero et al. 2023), although the pressure spectrum in their Ring 1 yields a larger Mach number from the integrated relation than the peak relation. Thus it appears that scaling relations derived from high-resolution hydrodynamical simulations (Gaspari & Churazov 2013; Gaspari et al. 2014) are in agreement with those derived from large-scale cosmological simulations (Zhuravleva et al. 2014, 2023). An outlier to this converged trend is the scaling from Simonte et al. (2022) cosmological runs, who find a shallower sample relation (by roughly $2\times$). This can be understood by the different ad hoc filtering of “clumps”/substructures (not driven by turbulence) which must be adopted in large-scale runs with necessarily coarse resolution. For example, in Zhuravleva et al. (2014), the clump filtering is minimal with a 3.5 rms cut of the right tail in the density distribution, implying that 99.98% of the distribution is retained. In Simonte et al. (2022), the filtering is instead aggressive, with a cut above ~1.6 rms. They report that the related density fluctuation threshold cut is $\sigma_{\rho,3D}^2 = 0.22$, implying that the potentially included Mach numbers remain below 0.5, hence the retrieved low scaling normalization and small pressure support in their study, also compared with our findings.

We know that our recovered spectrum of fluctuations cannot be interpreted as solely due to turbulence, especially in the case of SPT-CLJ0638-5358, which has a known shock (Botteon et al. 2018) and for which we see suggestive hints of edges in the SZ (SPT) data (see Section 4). Nonetheless, the observed X-ray surface brightness fluctuations from the known shock are subdominant (see Appendix D) though the known shock may contribute to flattening the amplitude spectrum at smaller scales than $l_{inj} \approx 800$ kpc, but otherwise does not significantly impact the results we present (i.e., Mach numbers or hydrostatic mass biases). Were we to probe the turbulent cascade (with deeper data), it would be necessary to mask the shock to constrain the fluctuations due to turbulent motion properly.

We find that the hydrostatic bias for $M_{500}$ (those for Ring 2) derived for SPT-CLJ0232-4421 is well within the expected range (e.g., Hurier & Angulo 2018), while the hydrostatic bias (for Ring 2, i.e., $M_{500}$) for SPT-CLJ0638-5358 is nearly zero, which is plausible but rare for clusters. It is generally expected that a hydrostatic mass will underestimate the mass because it fails to account for nonthermal pressure. However, this perspective may overlook the conditions in which a hydrostatic mass may overestimate the mass (which are likely to be conditions in which pressure equilibrium itself is not entirely accurate). In particular, studies such as those by Wik et al. (2008) and Krause et al. (2012) show that both the central

Compton y and integrated Compton Y parameter (proxies for mass) will increase rapidly with the generation of merger shocks. In the case of a boosted SZ signal, we can expect that a hydrostatic mass derived either from a scaling relation (e.g., $Y-M$) or through an explicit calculation of thermal pressure support will have boosted Y and $P_e$ values and would overestimate the true mass. By extension, inferences from X-ray observations should be similarly affected: if $P_e$ is being boosted, under roughly adiabatic conditions, then both $n_e$ and $T_e$ will be boosted.

Numerical simulations which have estimated the hydrostatic mass bias have found average bias values of $0.1 \lesssim b \lesssim 0.3$ (e.g., Nagai et al. 2007; Lau et al. 2009; Rasia et al. 2012; Biffi et al. 2016; Ruppin et al. 2019; Gianfagna et al. 2023; Jennings & Davé 2023). However, they all find that the hydrostatic mass bias will have a distribution (and will also depend on the radius or density contrast at which the mass is being measured), and nearly all of these studies find that for an individual cluster, the bias for $M_{500}$ can be close to zero, and in some cases even negative. Given that they all find (positive) nonthermal pressure support, the small biases cannot be due to a lack of nonthermal pressure support. Moreover, Ruppin et al. (2019) classify clusters into relaxed and disturbed and find that both have a median hydrostatic mass bias $b=0.29$ while the standard deviation of the hydrostatic mass bias is narrow for the relaxed sample and broad for the disturbed sample. In other words, simulations appear to corroborate the case for sufficient departures from equilibrium resulting in low (and even negative) hydrostatic mass biases.

## 6.2. Other Hydrostatic Mass Estimates for Our Clusters

We present here a variety of alternative calculations of the hydrostatic mass bias leaning heavily on mass estimates in the literature which have been compiled in Appendix F. Our additional methods of calculating a hydrostatic bias are (1) employing the gas fraction, $f_{gas}$, (2a) comparing hydrostatic masses (from X-ray data) to total masses as calculated by weak lensing, and (2b) comparing hydrostatic masses (from SZ data) to calibrated total masses; i.e., some scaling of the SZ data with either a lensing calibration (e.g., Hilton et al. 2021; H21 hereafter) or via abundance matching (e.g., Bleem et al. 2015; de Haan et al. 2016; Bocquet et al. 2019; Bulbul et al. 2019; B15, H16, B19, and Bu19 hereafter, respectively).

In the case of deriving a hydrostatic mass bias from $f_{gas}$, we employ the relation presented in Allen et al. (2008) and Wicker et al. (2023), which gives a theoretical expectation for the gas fraction:

$$f_{\text{gas,Theory}} = K \frac{\Upsilon(M, z)}{1-b} A(z) \left(\frac{\Omega_b}{\Omega_m}\right) \left(\frac{D_A^{\text{ref}}(z)}{D_A(z)}\right)^{3/2} - f_*, \quad (18)$$

where K is an instrumental calibration factor, $\Upsilon(M, z)$ is the baryon depletion factor, b is the hydrostatic mass bias ($b=1-M_{HE}/M_{tot}$), $A(z)$ is an angular correction, $\Omega_b/\Omega_m$ is the universal baryon fraction, $D_A$ is the angular diameter distance, and $f_*$ is the stellar mass fraction. As we are not comparing different cosmologies (and the works from which we derive our masses assume the same cosmology), the terms A (z) and $D_A^{\text{ref}}(z)/D_A(z)$ can be taken as unity. We follow Wicker et al. (2023) and take $K=1 \pm 0.1$, $\Omega_b/\Omega_m=0.156 \pm 0.03$ (Planck Collaboration et al. 2020), $\Upsilon(M, z)=\Upsilon_0=0.85 \pm 0.03$
<figure><img src="image_16.png" /><figcaption>**Figure 10.** Hydrostatic mass biases calculated via various methods for SPT-CLJ0232-4421 (top) and SPT-CLJ0638-5358 (bottom). We remind the reader that the sign convention for the hydrostatic mass biases calculated via our fluctuation analyses, $b_M$ is opposite that of the more common $b = 1 - M_{HE}/M_{tot}$. Biases are grouped by method: the top group are those derived from this work and estimation from spectra of surface brightness fluctuations and the second grouping is derived from the gas fraction. The third grouping shows biases when comparing HE mass estimates to total (weak lensing (WL)) mass estimates from Fox et al. (2022). The final grouping shows HE mass estimates relative to cosmological abundance-matched total masses. Citation abbreviations (indicated in the text) reflect corresponding references for HE and total (or abundance-matched) masses in the third and fourth groupings, or $M_{gas}$ and $M_{HE}$, respectively, for biases from $f_{gas}$. The mass estimate sub-scripted with $A_{10} \times R_{23}$ is an HE mass estimate described below and in Appendix F.</figcaption></figure>

(Planelles et al. 2013), and $f_* = 0.015 \pm 0.005$ (Eckert et al. 2019). Mass biases derived from $f_{gas}$ in Figure 10 make use of $M_{gas}$ and $M_{HE}$ either from a single paper (if both gas mass and HE mass are reported in that paper), or the gas mass in one paper, cited first, and the HE mass from the second paper.

We choose this method to derive a hydrostatic mass bias from $f_{gas}$, as opposed to the method employed in Eckert et al. (2019), because it requires no estimate of the nonthermal pressure and its derivative. Here again, assumptions would need to be made (assume a nonthermal pressure profile as given in Nelson et al. 2014) or we would need to incorporate values from our fluctuation analysis, meaning the resulting hydrostatic mass bias would not be independent of our fluctuation analysis.

The remaining method(s) of hydrostatic mass bias is simply using its definition ($b = 1 - M_{HE}/M_{tot}$) and varying what we

take as a hydrostatic mass estimate and what we take as a total mass estimate. In one set, we take predominantly X-ray-derived masses and compare to mass estimates from weak lensing (Fox et al. 2022). In this category, we include an estimate of $M_{500}$ by fixing the pressure profile shape to the universal pressure profile (UPP) shape found in Arnaud et al. (2010) and fitting the profile to the SPT data, allowing only $M_{500}$ as a free parameter, which artificially constrains the uncertainty (see Appendix F). Finally, we also compare mass estimates from SZ surveys which are derived, in basis, from the assumption of HE to mass estimates from SZ surveys where the mass is estimated from either weak-lensing calibration or abundance matching. To avoid exploring all possible permutations, Figure 10 shows only the lowest and highest biases one can calculate from comparable hydrostatic and total mass estimates for an assumed cosmology. This is discussed further in Appendix F. Remaining citation abbreviations are Z08 and M10, which refer to Zhang et al. (2008) and Mantz et al. (2010), respectively.

Figure 10 then compares hydrostatic mass biases as calculated by our fluctuation analyses ($-b_M$) to the other methods laid out in this section. We see the expected clustering of hydrostatic mass biases from SZ + HE to SZ abundance matched between $0.25 < b < 0.55$. That is, the mass biases when the total mass is matched to be consistent with a cosmology tend toward higher bias values ($b$), such as that found in Planck Collaboration et al. (2016a; P16 hereafter; $1 - b = 0.58 \pm 0.04$). By comparison, the other methods find bias values much closer to zero (for these two clusters), with a scattering of positive and negative values. Thus, beyond the statistical agreement of our fluctuation-derived hydrostatic mass bias values (in part due to large uncertainties, especially for SPT-CLJ0638-5358), the values themselves are easily within the range of values one can derive from mass estimates for individual clusters in the literature. Thus, a more significant comparison will arise when comparing the distribution of hydrostatic mass biases derived from fluctuations over the full SPT and XMM sample to the distribution from other methods of estimating the hydrostatic mass bias.

## 7. Impact of Choices on the Results

Every step in analyzing surface brightness fluctuations (and those of their corresponding deprojected thermodynamic quantities) presents choices. We have addressed major algorithmic choices in the corresponding analysis sections. Here we address choices (alternatively, assumptions) about the underlying morphology of the ICM and distribution of its fluctuations. In this section we consider the impact of (1) the assumed geometry for the surface brightness profiles and (2) which regions we use in calculating the spectra of our fluctuations.

### 7.1. Choice of Profile Fitting Parameters

For our surface brightness profiles we adopted circular $\beta$-models for both the SZ and X-ray images where we fixed the centers to the X-ray centroids. Each of these choices ([1] choosing to fit a $\beta$-model to each image, [2] adopting circular symmetry, and [3] using the X-ray centroid) have been justified earlier in this article. Given that the $\beta$-model appears quite sufficient for both data sets, we do not explore the range of other parameterizations available in the literature. However, we
<figure><img src="image_17.png" /><figcaption>**Figure 11.** Axial ratio and orientation vs. radius for SPT-CLJ0232-4421 (left) and SPT-CLJ0638-5358 (right) as a function of fitting radius about the cluster centroid. The lines in all subplots are as indicated in the legend of the upper left plot. Neither the axis ratio nor rotation angle are constant.</figcaption></figure>

do explore the impact of our center choice and elliptical geometry on the resultant X-ray power spectra.

For the center choice we consider the X-ray peak as an alternative to the centroid. Regarding an elliptical model, we have multiple choices as there is no singular choice of ellipse as illustrated in Figure 11. The axial ratios and rotation angles are calculated via pyproffit<fnref n="23" /> within a circle of radius *r* about the centroid for each image. We vary *r* in steps of 0.5 out to R₅₀₀.

For each cluster, we average the elliptical parameters at each radius from 1' out to R₅₀₀. We then follow the same procedures as in Section 3: we fit profiles extracted from each image (where we adopt the average parameters of the ellipse(s) across EPIC cameras with the center taken to be the circular centroid used in Section 3), create fractional residual images, and calculate power spectra within the same regions (rings) as before. For these elliptical models, we assume the cluster is a prolate ellipsoid with its major axis in the plane of the sky. We further assume that the deprojection can be calculated with a weighted average of N (N is defined in Equation (12)) based on the semimajor axis (rather than circular radius; see Appendix B for further details.)

We find that the choice of center and ellipticity will impact the recovered 3D amplitude spectra as evidenced in Figure 12. Accounting for ellipticity does reduce the recovered power spectra as expected, and Figure 13 shows a trend that the larger the axial ratio used, the more reduced the amplitude. The reduction in amplitude is not uniform for the two clusters; that is, we cannot write a single equation relating the reduction in amplitude (at a given wavenumber) to the axial ratio that holds for both clusters. Furthermore, we note that, for a given cluster, the shapes of the amplitude spectra are quite similar among the circular or elliptical profiles fitted. That is, this inferred injection scale does not appear to be sensitive to the geometry used. However, these conclusions are drawn from a sample of merely two clusters which are not grossly elongated. This may extend to other clusters which are not particularly elongated, though we may anticipate that for larger axial ratios the impact is more severe and we should expect it to modulate the shape of the amplitude spectra, especially at lower *k*.

We briefly digress into a discussion about the impact of elongation along the line of sight on the projection of power spectra; in particular, for the case that an axis is aligned with the line of sight, we can simply scale *W* (from Equations (9) and (10); we omit subscripts as the same scaling applies to both window functions) and *z*. In particular, let *z'* = *z*\**c*, i.e., the cluster is elongated by a factor *c* along the line of sight. In this case, *W*'(θ, *z*') = *W*(*θ*, *z*)/*c*. Their respective Fourier transforms yield the relation $\tilde{W}'(k_z/c) = \tilde{W}(k_z)$. Integrating the

<footnote n="23">²³ This, again, accounts for point-source masking and the exposure map. Rotation angles are the degrees counterclockwise from west (such that 90° is north). The fitting of the elliptical parameters uses all (none masked) pixels in the defined circle.</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 12.** Recovered 3D magnitude spectra for SPT-CLJ0232-4421 (left) and SPT-CLJ0638-5358. The solid black and blue curves are our reference (circular model centered on the centroid) model, the dashed curves correspond to circular fits about the X-ray peak, and the solid gray and cyan curves are for the various elliptical fits (about the centroid). The black and gray curves correspond to the inner ring (i.e., circle) and the blue and cyan curves correspond to the outer ring (annulus).</figcaption></figure>

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

<figure><img src="image_20.png" /><figcaption>**Figure 13.** A comparison of the recovered magnitude spectra of Ring 2 (ordinate) relative to the recovered magnitude spectra of Ring 1 (abscissa) for the two lowest wavenumbers (squares and stars; see legends) color coded by axial ratio.</figcaption></figure>

One of the main goal of this work is also to derive hydrostatic mass biases, so we proceed from our 3D spectra as in Section 6 and derive hydrostatic biases for the elliptical cases. To do this, we have adopted the same equivalences as in Section 6 between the Mach number and the peak of the amplitude spectra and the integrated power spectrum (i.e., variance). Moreover, we use the same logarithmic pressure slopes as in the circular case (see Section 6 and Table 13). We maintain the $2\sigma$ threshold when selecting a peak, while $\langle\eta_p\rangle$ takes values of 1.0 and 1.3 for in-between and unrelaxed clusters with an ellipsoidal geometry, respectively. In principle, the formulation for total pressure balance with gravitational potential should be reworked in each ellipsoidal coordinates. Though, we also lack a convention for assigning $R_{500}$ to an ellipsoidal cluster (i.e., do we still take some sphere, or do we pick one of the axes as our length metric?). Notwithstanding these caveats, our calculations suggest that even with reduced amplitude spectra (and hence reduced Mach numbers), there is less of a clear trend in hydrostatic mass biases. For SPT-

square of $|\tilde{W}'(k_z)|$ over $dk_z$ (Equation (12)), it is evident that $N'$ (the approximate scaling between $P_{2D}$ and $P_{3D}$ in the case of elongation along the line of sight) is reduced by a factor of $c$ relative to $N$. The relative bias between the full projection along the line of sight (Equation (7)) and the approximate equation (Equation (11)) as shown in Appendix B.2 remains unchanged with elongation. In short, the primary effect of elongation along the line of sight would be to boost the recovered $P_{3D}$ by a factor of $c$. If we consider the axial ratios (Table 8) in the plane of the sky serve as a naive upper bound for $c$, we find that the recover $A_{3D}$ will at most be boosted by ~10%, and thus we suspect that elongation along the line of sight is not a significant concern for these two clusters.

For the circular case, the impact of center choice (peak versus centroid) is mixed for Ring 1, but yields reduced amplitudes in Ring 2. Note that choosing the center to be the peak will also shift (translate) the center of the rings so that the spectra are taken from slightly different regions of the sky. All spectra from elliptical fits are derived from the same regions as the circular centroid case of $\phi = (1 + \sqrt{5})/2$.
<table><caption>**Table 8**
Hydrostatic Mass Biases as Derived for Different Ellipticities</caption><thead><tr><th rowspan="2">R<sub>ell</sub></th><th colspan="4">SPT-CLJ0232-4421</th><th colspan="4">SPT-CLJ0638-5358</th></tr><tr><th>a/b</th><th>θ<sub>R</sub></th><th>Ring</th><th>-b<sub>M,peak</sub></th><th>-b<sub>M,int</sub></th><th>a/b</th><th>θ<sub>R</sub></th><th>Ring</th><th>-b<sub>M,peak</sub></th><th>-b<sub>M,int</sub></th></tr></thead><tbody><tr><td rowspan="2">φ = √3</td><td rowspan="2">1</td><td rowspan="2">0</td><td>Ring 1</td><td>0.21 ± 0.39</td><td>0.25 ± 0.20</td><td rowspan="2">1</td><td rowspan="2">0</td><td>Ring 1</td><td>-0.20 ± 0.23</td><td>-0.02 ± 0.04</td></tr><tr><td>Ring 2</td><td>0.25 ± 0.17</td><td>0.06 ± 0.04</td><td>Ring 2</td><td>0.03 ± 0.38</td><td>0.05 ± 0.06</td></tr><tr><td rowspan="2">φ = (1 + √5)/2</td><td rowspan="2">1</td><td rowspan="2">0</td><td>Ring 1</td><td>0.34 ± 0.20</td><td>0.26 ± 0.26</td><td rowspan="2">1</td><td rowspan="2">0</td><td>Ring 1</td><td>-0.42 ± 0.65</td><td>-0.01 ± 0.09</td></tr><tr><td>Ring 2</td><td>0.13 ± 0.10</td><td>0.05 ± 0.04</td><td>Ring 2</td><td>0.02 ± 0.73</td><td>0.08 ± 0.10</td></tr><tr><td rowspan="2">φ = √2</td><td rowspan="2">1</td><td rowspan="2">0</td><td>Ring 1</td><td>...</td><td>0.15 ± 0.13</td><td rowspan="2">1</td><td rowspan="2">0</td><td>Ring 1</td><td>...</td><td>0.10 ± 0.17</td></tr><tr><td>Ring 2</td><td>...</td><td>0.04 ± 0.04</td><td>Ring 2</td><td>...</td><td>0.03 ± 0.09</td></tr><tr><td rowspan="2">Peak</td><td rowspan="2">1</td><td rowspan="2">0</td><td>Ring 1</td><td>...</td><td>0.42 ± 0.25</td><td rowspan="2">1</td><td rowspan="2">0</td><td>Ring 1</td><td>-0.21 ± 0.34</td><td>-0.00 ± 0.06</td></tr><tr><td>Ring 2</td><td>...</td><td>0.03 ± 0.02</td><td>Ring 2</td><td>0.10 ± 0.38</td><td>0.06 ± 0.07</td></tr><tr><td>1.0</td><td>1.14</td><td>140</td><td>Ring 1</td><td>0.28 ± 0.19</td><td>0.07 ± 0.06</td><td>1.15</td><td>164</td><td>Ring 1</td><td>-0.22 ± 0.34</td><td>-0.01 ± 0.07</td></tr><tr><td></td><td></td><td></td><td>Ring 2</td><td>0.07 ± 0.08</td><td>0.01 ± 0.02</td><td></td><td></td><td>Ring 2</td><td>0.05 ± 0.45</td><td>0.06 ± 0.08</td></tr><tr><td>1.5</td><td>1.07</td><td>141</td><td>Ring 1</td><td>0.28 ± 0.19</td><td>0.07 ± 0.06</td><td>1.17</td><td>158</td><td>Ring 1</td><td>-0.20 ± 0.30</td><td>-0.01 ± 0.07</td></tr><tr><td></td><td></td><td></td><td>Ring 2</td><td>0.09 ± 0.08</td><td>0.02 ± 0.02</td><td></td><td></td><td>Ring 2</td><td>0.06 ± 0.40</td><td>0.06 ± 0.08</td></tr><tr><td>2.0</td><td>1.06</td><td>139</td><td>Ring 1</td><td>0.30 ± 0.19</td><td>0.07 ± 0.06</td><td>1.17</td><td>150</td><td>Ring 1</td><td>-0.22 ± 0.30</td><td>-0.01 ± 0.06</td></tr><tr><td></td><td></td><td></td><td>Ring 2</td><td>0.10 ± 0.08</td><td>0.02 ± 0.02</td><td></td><td></td><td>Ring 2</td><td>0.04 ± 0.43</td><td>0.06 ± 0.07</td></tr><tr><td>2.5</td><td>1.05</td><td>123</td><td>Ring 1</td><td>0.30 ± 0.20</td><td>0.08 ± 0.07</td><td>1.18</td><td>145</td><td>Ring 1</td><td>-0.20 ± 0.29</td><td>-0.01 ± 0.06</td></tr><tr><td></td><td></td><td></td><td>Ring 2</td><td>0.11 ± 0.09</td><td>0.02 ± 0.02</td><td></td><td></td><td>Ring 2</td><td>0.06 ± 0.39</td><td>0.06 ± 0.07</td></tr><tr><td>3.0</td><td>1.07</td><td>101</td><td>Ring 1</td><td>0.27 ± 0.18</td><td>0.07 ± 0.06</td><td>1.19</td><td>142</td><td>Ring 1</td><td>-0.20 ± 0.29</td><td>-0.01 ± 0.07</td></tr><tr><td></td><td></td><td></td><td>Ring 2</td><td>0.11 ± 0.08</td><td>0.02 ± 0.02</td><td></td><td></td><td>Ring 2</td><td>0.06 ± 0.39</td><td>0.06 ± 0.08</td></tr><tr><td>3.5</td><td>1.09</td><td>95</td><td>Ring 1</td><td>0.27 ± 0.19</td><td>0.06 ± 0.06</td><td>1.20</td><td>146</td><td>Ring 1</td><td>-0.19 ± 0.28</td><td>-0.01 ± 0.06</td></tr><tr><td></td><td></td><td></td><td>Ring 2</td><td>0.10 ± 0.08</td><td>0.02 ± 0.02</td><td></td><td></td><td>Ring 2</td><td>0.06 ± 0.38</td><td>0.06 ± 0.07</td></tr><tr><td>4.0</td><td>1.11</td><td>92</td><td>Ring 1</td><td>0.27 ± 0.19</td><td>0.07 ± 0.06</td><td>1.21</td><td>145</td><td>Ring 1</td><td>-0.17 ± 0.27</td><td>-0.01 ± 0.07</td></tr><tr><td></td><td></td><td></td><td>Ring 2</td><td>0.09 ± 0.08</td><td>0.02 ± 0.02</td><td></td><td></td><td>Ring 2</td><td>0.06 ± 0.37</td><td>0.06 ± 0.08</td></tr><tr><td>4.5</td><td>1.12</td><td>91</td><td>Ring 1</td><td>0.23 ± 0.18</td><td>0.05 ± 0.05</td><td>1.21</td><td>146</td><td>Ring 1</td><td>-0.16 ± 0.27</td><td>-0.00 ± 0.07</td></tr><tr><td></td><td></td><td></td><td>Ring 2</td><td>0.09 ± 0.08</td><td>0.02 ± 0.02</td><td></td><td></td><td>Ring 2</td><td>0.08 ± 0.36</td><td>0.06 ± 0.09</td></tr><tr><td>5.0</td><td>1.11</td><td>95</td><td>Ring 1</td><td>0.26 ± 0.18</td><td>0.06 ± 0.06</td><td>1.21</td><td>144</td><td>Ring 1</td><td>-0.18 ± 0.29</td><td>-0.01 ± 0.08</td></tr><tr><td></td><td></td><td></td><td>Ring 2</td><td>0.09 ± 0.08</td><td>0.02 ± 0.02</td><td></td><td></td><td>Ring 2</td><td>0.06 ± 0.38</td><td>0.06 ± 0.09</td></tr><tr><td>5.5</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>1.20</td><td>144</td><td>Ring 1</td><td>-0.18 ± 0.29</td><td>-0.01 ± 0.07</td></tr><tr><td></td><td></td><td></td><td>...</td><td>...</td><td>...</td><td></td><td></td><td>Ring 2</td><td>0.06 ± 0.38</td><td>0.06 ± 0.08</td></tr><tr><td>6.0</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>1.19</td><td>142</td><td>Ring 1</td><td>-0.16 ± 0.30</td><td>-0.00 ± 0.08</td></tr><tr><td></td><td></td><td></td><td>...</td><td>...</td><td>...</td><td></td><td></td><td>Ring 2</td><td>0.09 ± 0.38</td><td>0.06 ± 0.09</td></tr></tbody></table>

**Note.** R<sub>ell</sub> is the radius at which the ellipticity parameters are determined (as in Figure 11). a/b is then the axis ratio and θ<sub>R</sub> is the rotation angle (in degrees). The first four rows show the results from circular cases and either different ring extents or different centers; entries for φ = (1 + √5)/2 are the reprint of the hydrostatic biases in Table 7. All numbers for elliptical fits use φ = (1 + √5)/2 to define the ring extents. Empty entries for SPT-CLJ0638-5358 arise from a lack of nodes (in Ring 2) about 2σ; empty entries for SPT-CLJ0232-4421 reflect that R<sub>500</sub> = 5.3 and thus elliptical parameters are not pursued beyond R<sub>500</sub>.

CLJ0232, the elliptical results all suggest a lower hydrostatic bias; this is not true for SPT-CLJ0638-5358.

It is hopefully evident that when considering elliptical models, the parameter space is immediately inflated: geometrically one must choose which ellipticity is appropriate and what to assume about the third ellipsoidal axis and inclination. Similar to previous discussions (e.g., Zhuravleva et al. 2015; Romero et al. 2023), it is unclear what would be correct, though some choice is taken in Dupourqué et al. (2023) and Zhang et al. (2023), where we note Zhang et al. (2023) models a circular β-model for the core and uses an elliptical β-model for the broader emission. There are further considerations in terms of deprojection, and deriving a hydrostatic mass bias for an ellipsoidal potential and how to define R<sub>500</sub> in this framework. Certainly, for our sample we conclude that adopting a spherical cluster model centered on the X-ray surface brightness centroid is the most appropriate choice.

## 7.2. Choice of Region for Extracting Spectra

There are many factors to consider in choosing the regions within which to calculate the power spectra. Some obvious limitations come about from the depth of observations (i.e., the

noise in the images) and the angular resolution of the instrument(s) used. Similarly, if one wants to constrain fluctuations on large scales, then the region ought to be large enough to well sample those scales.

In our case, we wish to measure fluctuations out to R<sub>500</sub>, both in the radial sense and in the Fourier sense. Our desire to sample scales close to R<sub>500</sub> motivates an inner circle of at least r = R<sub>500</sub>/2. Conversely, the resolution of the instruments can also set a minimum width of annuli (1.25 in our case). Even with these bounds in place, for the two clusters here, there are a range of options for two annuli (i.e., inner circle and outer annulus) that fit these criteria, though these criteria effectively rule out using three (or more) annuli.

A simple geometric argument would be to divide the area equally between the inner circle and outer annulus. While this meets the above criteria, we again note that noise should be considered. That is, the noise for spectral calculations is not the map noise but the map noise divided by the surface brightness model. Thus, regions at larger radii will have larger intrinsic noise.

From here we see the choice of regions as a vague endeavor. To limit the parameter space searched, we consider three
<figure><img src="image_21.png" /><figcaption>**Figure 14.** A comparison of the recovered magnitude spectra as a function of the inner circle radius, noted in the legend of the left panel as the divisor of $R_{500}$. The same convention relating divisors and linestyles is used in the right panel. Top row: amplitude spectra for SPT-CLJ0232-4421. Top left: amplitude spectra of pressure fluctuations. Top right: amplitude spectra of density fluctuations. Bottom row: amplitude spectra for SPT-CLJ0638-5358. Bottom left: amplitude spectra of pressure fluctuations. Bottom right: amplitude spectra of density fluctuations.</figcaption></figure>

options. We define $R_1$ as the radius of the inner circle (equally the inner radius of our outer annulus) with $R_1 = R_{500}/\phi$; as such we try three values of $\phi$: [$\sqrt{2}$, $(1+\sqrt{5})/2$, $\sqrt{3}$]. Figure 14 shows that indeed the region choice is mostly inconsequential. However, for SPT-CLJ0638-5358, we do see that there is a significant difference in the amplitude spectra of Ring 2, where $\phi = \sqrt{2}$ results in a lower spectrum than the other two values of $\phi$. This shows that there is a residual feature between $\phi = \sqrt{2}$ and $\phi = (1+\sqrt{5})/2$ that is relevant. In particular, we infer that this feature is the plausible shock feature(s) to the northeast of the cluster center (highlighted in Figure 5). The impact of $\phi$ on the resultant hydrostatic mass bias is difficult to infer with some values missing due to data quality (see Table 8). The values presented in Table 8 suggest that this choice is not particularly important, especially for these two clusters.

## 8. Conclusions

In this work we have outlined a methodology to study surface brightness fluctuations of both SZ and X-ray data that will be applied to a sample of SPT-selected galaxy clusters which have archival XMM-Newton data in future work. Our goal is to constrain the turbulent properties in the ICM. Given the quality of the data that we analyze, several specific goals arise: constraining the amplitude of fluctuations, inferring turbulent Mach numbers, and deriving a hydrostatic mass bias from the Mach numbers. To pilot this methodology we have studied two massive ($M_{500} \approx 1 \times 10^{15} M_{\odot}$) clusters at redshifts $0.2 < z < 0.3$. A shock is known to exist in one system (SPT-CLJ0638-5358; Botteon et al. 2018), but the direct impact of that shock on the fluctuations is subdominant relative to other fluctuations in the corresponding region. For our baseline
results, we adopt a circular surface brightness model (spherical ICM model) and subsequently we explore the impact of using elliptical surface brightness models.

The amplitude of fluctuations we recover in these two clusters is larger than those seen in several other studies (e.g., Churazov et al. 2012; Gaspari & Churazov 2013; Zhuravleva et al. 2015; Arévalo et al. 2016; Hofmann et al. 2016; Eckert et al. 2017; Hernández-Lang et al. 2023), but in line with other studies which probe similar physical scales (e.g., Khatri & Gaspari 2016; Dupourqué et al. 2023). We estimate the Mach number for each ring based on three methods (scaling relations): (1) relating the peak of the amplitude spectrum to the Mach number, (2) relating the integrated power spectrum to the Mach number, and (3) combining the two previous methods by taking the average of their respectively derived Mach numbers. The Mach numbers inferred from the two separate scaling relations are generally in agreement, except for Ring 2 of SPT-CLJ0638-5358 which has a high peak in its amplitude spectrum of density fluctuations and a comparatively low integrated power spectrum of density fluctuations.

We calculate hydrostatic biases for each set of Mach numbers (i.e., for each method of deriving a Mach number) and find similar bias values across the three methods (as calculated for Ring 2) in each cluster, which we interpret as the hydrostatic bias for $M_{500}$ of each cluster. While the two clusters studied here, SPT-CLJ0232-4421 and SPT-CLJ0638-5358, are ostensibly dynamically relaxed and disturbed, respectively, they both have fairly low hydrostatic biases. Moreover, the disturbed system has the lower hydrostatic bias. The disturbed cluster has greater nonthermal pressure than the relaxed cluster, as inferred from the larger Mach numbers. However, the steep logarithmic slope of the Mach numbers in the disturbed cluster yields a lower hydrostatic bias. If indeed the hydrostatic mass bias is low, especially for SPT-CLJ0638-5358, we interpret this as an endorsement of the notion that hydrostatic mass estimates can themselves be (transiently) boosted from a corresponding (transient) boost in the underlying SZ or X-ray signal.

The inferred hydrostatic mass biases were relatively insensitive to the extent of Ring 1 relative to Ring 2, with the exception of SPT-CLJ0638-5358, where toggling the ring extents shows that prominent fluctuations appear beyond $R_{500}/\sqrt{2}$, notably in the northwestern quadrant. When fitting elliptical models, our inferred fluctuations decrease by 5%–40%. However, our hydrostatic mass bias estimates did not change in a corresponding fashion. This is due to the fact that the logarithmic slope of the Mach numbers did not change dramatically and it plays a critical role in modulating the derived value of the hydrostatic mass bias. This may suggest that the hydrostatic mass bias can be robustly calculated across different assumed elliptical geometries, though this needs to be assessed across a broader sample, as will be done in our future study of the full sample of SPT-selected clusters with archival XMM-Newton data.

## Acknowledgments

The authors would like to thank the anonymous referee for comments which have improved this work. C.R. supported by NASA ADAP grant 80NSSC19K0574 and Chandra grant G08-19117X. E.B. acknowledges financial support from the European Research Council (ERC) Consolidator grant under the European Unions Horizon 2020 research and innovation program (grant agreement CoG DarkQuest No. 101002585).

M.G. acknowledges partial support by HST GO-15890.020/ 023-A, the BlackHoleWeather program, and NASA HEC Pleiades (SMD-1726). R.K. acknowledges support from the Smithsonian Institution, the Chandra High Resolution Camera Project through NASA contract NAS8-03060, and NASA Grants 80NSSC19K0116, GO1-22132X, and GO9-20109X. P. N. was supported by NASA contract NAS8-03060. The South Pole Telescope program is supported by the National Science Foundation (NSF) through award OPP-1852617. Partial support is also provided by the Kavli Institute of Cosmological Physics at the University of Chicago. Work at Argonne National Lab is supported by UChicago Argonne LLC, Operator of Argonne National Laboratory (Argonne). Argonne, a U.S. Department of Energy Office of Science Laboratory, is operated under contract No. DE-AC02-06CH11357.

*Facility:* SPT XMM-Newton.

*Software:* astropy (Astropy Collaboration et al. 2013; The Astropy Collaboration 2018; Astropy Collaboration et al. 2022), emcee (Foreman-Mackey et al. 2013), pyprofit (Eckert et al. 2017), and ESAS (Snowden et al. 2008).

## Appendix A
### The Importance of Modeling a Cool Core

In Section 3, we note that the β-model leaves residuals in the radial profile of SPT-CLJ0232-4421 owing to its cool core. Here, we investigate the impact on the resultant spectra when adopting a model with more free parameters. In particular, we adopt a double β-model as parameterized in pyprofit:

$$S(r) = S_0[(1 + (r/r_{c,1})^2)^{-3\beta+0.5} + R(1 + (r/r_{c,2})^2)^{-3\beta+0.5}] + B, \quad (A1)$$

where $S_0$, $\beta$, and $B$ are as in Equation (4), and if we make a restriction that $r_{c,2} < r_{c,1}$, then we may consider $r_{c,2}$ to be the characteristic radius of the cool core (and $r_{c,1}$ could be considered the extended core radius). The $R$ parameter allows for a ratio between the normalizations of the two β-models. Note, the model makes no restriction of which $r_c$ is smaller; the choice is arbitrary.

We perform the profile fits as in Section 3 and find that the residuals, as seen in the residual radial profiles (Figure 15), are notably reduced. While Figure 15 shows just profiles for ObsID 0042340301 in the 0.4–1.25 keV band, the results (especially the residuals) are very similar for both ObsIDs and both energy bands. We note that $\beta \sim 0.8$ across the CCDs, ObsIDs, and two energy bands with the double β-model. In comparison, the $\beta$ values were ~0.56 when fitting the single β-model. Although we do not wish to scrutinize the fitted parameter values, we highlight the larger $\beta$-value as it indicates steeper slopes, which is highly relevant. The lower panel in Figure 15 shows the resultant amplitude spectra (comparable to Figure 9) when using the double β-model for SPT-CLJ0232-4421 (the SZ-derived spectrum remains unchanged). Table 9 clearly quantifies the relevant spectral parameters; when comparing with Table 5, it is clear that both rings show an increase in fluctuations. How can this be?

The larger $\beta$-value noted earlier is the fundamental reason for this increase in fluctuations. Of course, it is not the sole culprit; the large region plays a role as well. Figure 16 shows several profiles which help to elucidate the mechanisms at play. The top panel shows the profiles of the surface brightness profiles with the background subtracted. This is the profile by
<figure><img src="image_22.png" /><figcaption>**Figure 15.** Left: surface brightness profiles for SPT-CLJ0232-4421 by CCD for ObsID 0042340301 in the 0.4–1.25 keV band. Right: resultant amplitude spectra when adopting a double β-model; see Figure 9.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 16.** Profiles from SPT-CLJ0232-4421. In all panels, the blue, orange, and green colorings correspond to the MOS1, MOS2, and PN CCDs (as in Figure 15). Top: a comparison of single β-model vs. double β-model surface brightness profiles without the contribution from the background. Middle: normalized variance within annuli, i.e., Var(δS/S) within the same annuli used to derive the observed surface brightness profiles. Var$_{2β}$ = Var(δS/S) for the case of S̄ is defined by the double β-model, and Var$_{β}$ is the corresponding normalized variance for the single β-model. The black curve denotes the area per annulus. Bottom: a comparison of the cumulative, area-weighted, variance of the normalized residuals from the double β-model relative to that of the single β-model. The magenta dashed line denotes the division between Ring 1 and Ring 2; the black dashed line denotes R$_{500}$ (and the outer edge of Ring 2).</figcaption></figure>

which the residual image is divided to obtain the normalized fluctuation image. Consequentially, for equal fluctuations (δS) at large radii, the normalized fluctuations (δS/S̄) will be larger for the profile with the steeper slope. This is seen in the middle panel, where the dashed lines correspond to the normalized fluctuations for the double β-model. The bottom panel compares the cumulative sum of the product of the colored

<table><caption>**Table 9**
Key Properties of the Amplitude Spectra of SPT-CLJ0232-4421 When the Surface Brightness Profile Model is a Double β-model</caption><thead><tr><th></th><th>A<sub>3D</sub>(k<sub>peak</sub>)</th><th>σ<sub>3D</sub></th><th>σ<sub>ln</sub></th><th>k<sub>peak</sub></th><th>λ<sub>peak</sub> (kpc)</th></tr></thead><tbody><tr><td>Ring 1</td><td>0.18 ± 0.04</td><td>0.21 ± 0.03</td><td>0.20 ± 0.03</td><td>0.49</td><td>516</td></tr><tr><td>Ring 2</td><td>0.33 ± 0.17</td><td>0.19 ± 0.09</td><td>0.19 ± 0.10</td><td>0.30</td><td>837</td></tr></tbody></table>

Note. The $k_{peak}$ column is in units of arcmin$^{-1}$.

<table><caption>**Table 10**
Derived Hydrostatic Biases When Using the Double β-model for SPT-CLJ0232-4421</caption><thead><tr><th></th><th>-b<sub>M,peak</sub></th><th>-b<sub>M,int</sub></th><th>-b<sub>M,comb</sub></th></tr></thead><tbody><tr><td>Ring 1</td><td>-0.41 ± 1.30</td><td>0.12 ± 0.19</td><td>-0.06 ± 0.70</td></tr><tr><td>Ring 2</td><td>0.30 ± 0.28</td><td>0.09 ± 0.09</td><td>0.20 ± 0.29</td></tr></tbody></table>

curves and the black curve in the middle panel. That is, the cumulative area-weighted, normalized variance serves as a proxy for the power in the fluctuations. This proxy is more reflective of the power at smaller scales. Although this proxy does not fully reveal how the power changes at difference scales, it nicely demonstrates how an excess in power can occur despite the (unnormalized) residual profile (falsely) suggesting smaller fluctuations.

If we continue our analysis to the derivation of hydrostatic biases (see Table 10), we find that adoption of the double β-model makes a more substantial difference than the elliptical considerations if looking at the recovered values and ignore the uncertainties. However, the uncertainties are still large and thus statistically we should not be consider the differences significant.

## Appendix B: Deprojection Details

We note two concerns related to the deprojection of our 2D spectra to 3D: (1) how do we deal with the dependence of the window function on projected cluster-centric radius, $\theta$ and (2) how valid is the approximation (given in Equation (11)) to Equation (7).
<figure><img src="image_24.png" /><figcaption>**Figure 17.** Left: integrated window functions, N(θ), and weighted averages within rings for SPT-CLJ0232-4421. Right: same as the left, but for SPT-CLJ0638-5358. For the X-ray data N(θ) is plotted for each CCD and each ObsID by energy band indicated as by their range in eV.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Table 11**
Projection Bias Values at $k_{500}$</figcaption></figure>

<table><thead><tr><th rowspan="3">Data Set</th><th rowspan="3">α</th><th colspan="4">B<sub>A</sub>(k<sub>500</sub>)</th></tr><tr><th colspan="2">SPT-CLJ0232-4432</th><th colspan="2">SPT-CLJ0638-5531</th></tr><tr><th>Ring 1</th><th>Ring 2</th><th>Ring 1</th><th>Ring 2</th></tr></thead><tbody><tr><td rowspan="3">SZ</td><td>2</td><td>1.02</td><td>1.01</td><td>1.01</td><td>1.00</td></tr><tr><td>3</td><td>1.03</td><td>1.01</td><td>1.02</td><td>1.01</td></tr><tr><td>4</td><td>1.04</td><td>1.02</td><td>1.03</td><td>1.01</td></tr><tr><td rowspan="3">X-ray</td><td>2</td><td>1.09</td><td>1.02</td><td>1.09</td><td>1.02</td></tr><tr><td>3</td><td>1.13</td><td>1.03</td><td>1.13</td><td>1.03</td></tr><tr><td>4</td><td>1.16</td><td>1.04</td><td>1.16</td><td>1.04</td></tr></tbody></table>

### B.1. Averaging within a Region

If we, for the moment, take the approximation given by Equation (11) to be valid, we still have the issue that N(θ) is a function of θ (see Figure 17) and we have a range of values of θ within our annuli. We posit that the appropriate way to average N(θ) is to weight the values by the area subtended by an annulus of radius θ and width dθ. In practice, we sample 16 equally spaced lines of sight in the range [0, R$_{500}$] and weight by area between these bins, given the line of sight at θ = 0 zero weight. For each ring we calculate an effective N (shown as stars in Figure 17; hereafter N$_{eff}$) as the

weighted average of N(θ) as above for those θ which lie within the ring.

It is clear that N can vary substantially from the inner edge of an annulus (or R = 0 for the inner circle) to the outer edge of the respective annulus. As noted in Zhuravleva et al. (2015), the maximal differences can be used to set the maximal uncertainty in N$_{eff}$.

### B.2. Validity of Approximation

To test the validity of Equation (11) we can perform the full integration of Equation (7) with an assumed 3D power spectrum, which we take to be a power law with a cutoff at $k_c$:

$$P_{3D} = P_0 e^{k_c/k} k^{-\alpha}, \qquad (B1)$$

where $P_0$ is the normalization (for our purposes, this is arbitrary) and α is the spectral index. We take $k_c = 1/(5 R_{500})$, which is well away from the largest scales that we sample. As discussed in Section 5, the slope α = 3 is of notable concern when deriving the peak in the magnitude spectra, but we should also be concerned with slopes shallower and steeper than it. Figure 18 shows the bias curves for SZ and X-ray when assuming α = 3. Table 11 lists the corresponding biases induced by the projection approximation at $k_c = 1/R_{500}$ when considering α = 2, 3, or 4.
<figure><img src="image_26.png" /></figure>

**Figure 18.** Left: bias in magnitude spectra for SPT-CLJ0232-4421; the top panel is with respect to the SZ (SPT) window functions and the bottom panel is with respect to X-ray (XMM-Newton) window functions. Right: same as the left, but for SPT-CLJ0638-5358. In all panels the thin and faint lines are for a specific line of sight and assuming an underlying 3D spectral index of $\alpha = 3$. The thick and darker lines are the weighted averages within their respective rings. For the top (SZ) plots, the solid thick lines are the weighted averages with an underlying 3D spectral index of $\alpha = 3$; the dotted line corresponds to $\alpha = 2$ and the dashed corresponds to $\alpha = 4$. In the bottom panels the solid and dashed lines correspond to our LE (0.4–1.25 keV) and HE (2.0–5.0 keV) bands, respectively, and are nearly identical. The legend indicates the bias values at scales of $R_{500}$, which is marked in all plots as the red dotted vertical line.

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

## Appendix C
### Masking Effects

The $\Delta$-variance method presented in Arévalo et al. (2012, hereafter A12), like other $\Delta$-variance methods, is intended to deal with gaps in an image and thus is a major motivation for using this method on the XMM-Newton data which have chip gaps and require masking point sources. For simpler masks (e.g., a circular mask for our inner ring) an alternative is to apodize the circular mask (e.g. Koch et al. 2019), though there remains the choice of apodization and whatever choice is made, there is still an effect on the recovered power spectrum, principally in normalization.<fnref n="24" /> Apodization can similarly be applied to annular masks and in all cases the effect should be controlled for the specific mask and apodization choice for an assumed underlying power spectrum or power spectra. Apodization is not a panacea and ultimately the inclusion of chip gaps and point sources would become too much for this approach.

Thus, we return again to the $\Delta$-variance method presented in Arévalo et al. (2012) and opt to use this for both our XMM-Newton and SPT data sets. Again, while it is designed to work with arbitrary masks, and was shown in Arévalo et al. (2012) to recover the power spectra shape very well for arbitrary masks for a range of spectral indices, we attempt to account for any biases in the recovered spectra due to the masking. For instance, it should not be surprising that if a region (mask) does not cover a large spatial scale, then the recovered power spectrum will be underestimated at that scale (this can be seen in Appendix B of Romero et al. 2023). Here, we check our ability to recover power spectra of varying spectral indices, maintaining $1/k_c = 5 R_{500}$, as in the previous section, given our

choice of annuli by generating image realizations described by the various power spectra. We further test the recovery of these power spectra when we include chip gaps and point-source masking. Though we might expect the scale recovery to be independent of pixel size, we decided to test with pixelizations corresponding to those in the X-ray and SZ data sets.

It should be noted that for a given realization for a given scale $l = 1/k$ which may be “fully” sampled, as there is at least one circle of diameter $l$, the sample variance will be notable. Indeed, even when a scale is sampled by several circles of corresponding diameter there will be some variance. For this reason, we iterate over 200 realizations, sampling to larger scales than reported for our data in Section 5. The bias is taken as the average recovered power spectra divided by the input power spectra and accounting for the known spectral index normalization bias (Arévalo et al. 2012). Prior to spectral measurements we do not smooth our images by the instrument PSF; thus we do not need to be concerned with beam corrections (Romero et al. 2023). For a given input spectrum the sample dispersion ($\sigma$) is taken as the standard deviation of the spectral measurements across realizations and normalized by the input power spectrum at its normalization bias. We report biases ($B_{A12}$; see Figure 19) and dispersions ($\sigma_{A12}$;<fnref n="25" /> see Figure 20 relative to the (average) recovered magnitude spectra).

Figures 19 and 20 indicate that both the bias and dispersion depend on the underlying spectral index. Interestingly, there appears to be a dependence on the pixel size. This is likely due to the coarser mask with the coarser pixelization; the radii of concern are between 12 and 25 pixels in radius for the SPT

<footnote n="24">²⁴ See, for example, https://turbustat.readthedocs.io/en/latest/tutorials/applying_apodizing_functions.html</footnote>

<footnote n="25">²⁵ Note that true sample variance is independent of the method of power spectrum estimation used. However, the reported dispersion is itself measured with the A12 method, so we give it this subscript.</footnote>
<figure><img src="image_28.png" /></figure>

**Figure 19.** Colors indicate rings as in Figure 18. Left: bias in the magnitude spectra of SPT-CLJ0232-4421 due to masking of annuli in all cases. For the X-ray lines, the legend indicates if point sources and chip gaps are also masked ("masked") or not ("unmasked"). Right: same as the left, but for SPT-CLJ0638-5358. In all panels the thin and faint lines are for spectra indices $2 \le \alpha \le 4$; the bold lines indicate $\alpha = 3$. The scale (inverse) of $R_{500}$ is indicated by the red dotted vertical line.

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

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

**Figure 20.** Colors indicate rings as in Figure 18. Left: sample variance, or rather standard deviation $\sigma$, of the magnitude spectra for SPT-CLJ0232-4421 due to the masking of annuli in all cases. For the X-ray lines, the legend indicates if point sources and chip gaps are also masked ("masked") or not ("unmasked"). Right: same as the left, but for SPT-CLJ0638-5358. In all panels the thin and faint lines are for spectra indices $2 \le \alpha \le 4$; the bold lines indicate $\alpha = 3$. The scale (inverse) of $R_{500}$ is indicated by the red dotted vertical line.

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

(SZ) pixelization (0.25). We also see that accounting for the masking of chip gaps and point sources, while not dominant relative to the annular masking, is important. Though Figure 19 presents the sample variance relative to the scales, we can recast this to account for the area in each ring and find (see Figure 21) that the sample variance scales with the number of sampling elements, as it should.

While the sampling bias is dependent on the underlying spectral index and the sample variance appears independent of the underlying spectral index, both the sampling bias and sample variance should depend on the choice of rings. Though not discussed in Section 7, these biases should also be considered when choosing regions. We present sample bias and variance for both clusters using different rings in Figure 22.

The numbers presented in the legend refer to the factor, $\phi$, which is used to define the separation between the inner ring (circle) and outer annulus, $R_1 = R_{500}/\phi$. As suggested in Section 7, within some basic constraints, the exact choice of regions (annuli) does not seem to be critical when considering the recovery of a range of scales on noiseless (synthetic) data.

In this section, we have analyzed the recovery of mock spectra which are effectively simple power laws so as to remain relatively agnostic about the spectral shape. When exploring other spectral shapes, especially with values $k_c$ consistent with our inferred injection scales, we find small differences relative to our current approach. In the case of a higher-S/N spectrum, correcting for both spectral shape and masking may be warranted.
<figure><img src="image_32.png" /><figcaption>**Figure 21.** Rather than scaling the dispersion against scale, we scale against $\sqrt{N_k}$, where $N_k$ is the number of sampling elements (if taken as a square with length $1/(2\pi k)$). This definition of $N_k$ provides the relation $\sigma_{A12} = 1/\sqrt{N_k}$.</figcaption></figure>

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

<figure><img src="image_34.png" /><figcaption>**Figure 22.** Left: bias (upper plot) and sample variance, or rather standard deviation $\sigma$, (lower plot) of the magnitude spectra for SPT-CLJ0232-4421 due to the masking of annuli in all cases and masking chip gaps and point sources for the X-ray cases. Right: same as the left, but for SPT-CLJ0638-5358. The scale (inverse) of $R_{500}$ is indicated by the magenta dotted vertical line.</figcaption></figure>
## Appendix D
SPT-CLJ0638-5358 by Quadrant

The shock to the southwest of the cluster center (Botteon et al. 2018) in SPT-CLJ0638-5358 would in itself suggest there should be elongation of the cluster along the northeast-southwest axis, and indeed we find this (Figure 11). In Section 7 we discussed the impact of fitting and subtracting elliptical surface brightness models. Here we investigate keeping a circular model but further subdividing regions (see Figure 23) to allow for more localized interrogations of fluctuations.

<figure><img src="image_35.png" /><figcaption>**Figure 23.** The same as Figure 6 for SPT-CLJ0638 but now annotated to show the quadrants.</figcaption></figure>

In particular, we opt to subdivide our rings into quadrants and also revisit the choice of radius between Ring 1 and Ring 2. As is evidenced by some empty entries in Table 12, the data are not deep enough to produce significant (here taken as just a $2\sigma$ threshold) for all cases and quadrants explored. Even so, this analysis finds larger Mach numbers in the southwest and northeast quadrants. While this is true for Ring 1 it is more strikingly so in Ring 2, but the uncertainties are also larger in Ring 2. We find it notable that (1) the northeast quadrant, i.e., quadrant 4, recovers larger Mach numbers than quadrant 2 (wherein lies the known shock; Botteon et al. 2018) and (2) that the Mach number drops sharply when adopting $\phi = \sqrt{2}$ relative to the other ring separations. This indicates that there are prominent fluctuations between $3'1 < r < 3'3$ (potentially spanning a larger range than that). This could indicate some structure more aligned with the dashed chartreuse curve in Figure 6 than the solid chartreuse curve.

We note that the residuals responsible for the dominant fluctuations in quadrant 4 are thus at larger radii than the southeast shock found in Botteon et al. (2018); thus their search via Gaussian gradient magnitude (GGM) may not have yielded a detection in part due to the decrease in surface brightness at larger radii. Admittedly the fluctuations do not appear as a canonical brightness jump, which would also reduce the prominence in a GGM image. Indeed, we might take the ambiguity from the SPT residuals as being suggestive that the gas dynamics are more complicated than a classical merger shock, which even in favorable geometries such as in Abell 2146 (e.g., Russell et al. 2012, 2022), the heating of shocks can be quite complicated (Chadayammuri et al. 2022).

<table><caption>**Table 12**
Mach Numbers and Corresponding Hydrostatic Mass Biases; Empty Fields Arise Due to a Lack of Points in Our Amplitude Spectra above $2\sigma$</caption><thead><tr><th rowspan="2">Quadrant</th><th rowspan="2">φ</th><th colspan="2">Ring 1</th><th colspan="2">Ring 2</th></tr><tr><th>M<sub>3D,peak</sub></th><th>-b<sub>M</sub></th><th>M</th><th>-b<sub>M</sub></th></tr></thead><tbody><tr><td rowspan="3">Q1</td><td>√3</td><td>0.57 ± 0.08</td><td>-0.02 ± 0.22</td><td>0.88 ± 0.40</td><td>0.20 ± 0.19</td></tr><tr><td>(1 + √5)/2</td><td>0.68 ± 0.10</td><td>...</td><td>...</td><td>...</td></tr><tr><td>√2</td><td>0.74 ± 0.17</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">Q2</td><td>√3</td><td>0.77 ± 0.20</td><td>0.06 ± 0.36</td><td>1.07 ± 0.45</td><td>0.30 ± 0.23</td></tr><tr><td>(1 + √5)/2</td><td>0.82 ± 0.25</td><td>...</td><td>...</td><td>...</td></tr><tr><td>√2</td><td>0.97 ± 0.12</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">Q3</td><td>√3</td><td>0.41 ± 0.12</td><td>-0.10 ± 0.18</td><td>0.92 ± 0.46</td><td>0.10 ± 0.23</td></tr><tr><td>(1 + √5)/2</td><td>0.50 ± 0.12</td><td>...</td><td>...</td><td>...</td></tr><tr><td>√2</td><td>0.60 ± 0.21</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">Q4</td><td>√3</td><td>0.95 ± 0.19</td><td>-1.38 ± 3.25</td><td>2.26 ± 0.83</td><td>0.33 ± 0.53</td></tr><tr><td>(1 + √5)/2</td><td>0.93 ± 0.27</td><td>-1.36 ± 4.17</td><td>2.19 ± 0.99</td><td>0.25 ± 0.85</td></tr><tr><td>√2</td><td>1.07 ± 0.49</td><td>0.21 ± 0.72</td><td>1.30 ± 0.53</td><td>0.42 ± 0.31</td></tr></tbody></table>
## Appendix E
### Intermediate Data Products

Although we do make inferences from the logarithmic pressure and Mach slopes, we present them in Table 13 for

transparency in the values use to derive the hydrostatic mass biases that are reported in Table 7. In particular, these values are used in Equations (16) and (17).

<table><caption>**Table 13**
The Logarithmic Slopes as Inferred for the Circular Cluster Model and Presented in Section 6</caption><thead><tr><th></th><th></th><th>d ln P/d ln r</th><th>d ln M<sub>peak</sub>/d ln r</th><th>d ln M<sub>int</sub>/d ln r</th><th>d ln M<sub>comb</sub>/d ln r</th></tr></thead><tbody><tr><td rowspan="2">SPT-CLJ0232-4421</td><td>Ring 1</td><td>-0.75 ± 0.01</td><td>-0.511 ± 0.70</td><td>-0.91 ± 1.07</td><td>-0.66 ± 1.07</td></tr><tr><td>Ring 2</td><td>-3.42 ± 0.09</td><td></td><td></td><td></td></tr><tr><td rowspan="2">SPT-CLJ0638-5358</td><td>Ring 1</td><td>-1.21 ± 0.01</td><td>1.50 ± 0.83</td><td>0.73 ± 0.72</td><td>1.26 ± 1.23</td></tr><tr><td>Ring 2</td><td>-3.04 ± 0.11</td><td></td><td></td><td></td></tr></tbody></table>

**Note.** In particular these values are used in Equations (16) and (17).
## Appendix F
### Mass Estimates for Our Pilot Clusters

The clusters chosen for our pilot survey are quite massive ($M_{500} \sim 10^{15} M_{\odot}$) and consequently have been observed by many facilities and have many mass estimations. We do not intend to review the details of each mass estimation method used here, but broadly characterize the mass estimation as either a total mass estimate or a hydrostatic mass (whether explicitly calculating mass via HE or building of some relation which was established using the assumption of HE). Our goal is then to contextualize our hydrostatic mass bias values obtained through our fluctuation analysis with hydrostatic mass bias values obtained by other methods. In addition to masses available in the literature, we perform another version of our pressure profile

fitting (Section 4.2) where we fit a UPP (Arnaud et al. 2010) such that it is a function of $M_{500}$ and redshift, $z$, where the former is allowed to vary and the latter is fixed. We report our obtained values of $M_{500}$ in Table 14 and note that by assuming a fixed profile shape, the associated uncertainties are artificially reduced.

We have tried to identify groups of mass estimates within which we consider the comparison of hydrostatic masses (denoted with HE) to total masses to be appropriate and have color coded these green, teal, or red; in the red sample we also include the gas masses which can be used to estimate the hydrostatic bias with an accompanying HE mass. Furthermore, in the red sample, our desire is to have a total mass estimate from lensing data. Unfortunately SPT-CLJ0638-5358 appears to be without a lensing mass estimate in the literature. In lieu of this, we have used the mass as determined by a weak-lensing

<table><caption>**Table 14**
Obtained Values of $M_{500}$</caption><thead><tr><th>Reference</th><th>Facility(ies)</th><th>Cosmology</th><th>Quantity</th><th>Type</th><th>SPT-CLJ0232-4421</th><th>SPT-CLJ0638-5358</th></tr></thead><tbody><tr><td>B15</td><td>SPT (A)</td><td>Fiducial (a), (b)</td><td><i>M</i><sub>500</sub></td><td>Total</td><td>12.01<sup>+1.80</sup><sub>-1.80</sub></td><td>12.01<sup>+1.81</sup><sub>-1.81</sub></td></tr><tr><td>H21</td><td>ACT</td><td>Fiducial (a)</td><td><i>M</i><sub>500,uncorr</sub></td><td>HE</td><td>10.12<sup>+2.12</sup><sub>-1.75</sub></td><td>10.47<sup>+2.33</sup><sub>-1.91</sub></td></tr><tr><td>H21</td><td>ACT</td><td>Fiducial (a)</td><td><i>M</i><sub>500,c</sub></td><td><b>HE</b></td><td>8.73<sup>+1.71</sup><sub>-1.43</sub></td><td>8.91<sup>+1.85</sup><sub>-1.53</sub></td></tr><tr><td>H21</td><td>ACT<sup>(d)</sup></td><td>Fiducial (a)</td><td><i>M</i><sub>500,cal</sub></td><td><b>Total</b></td><td>12.30<sup>+2.70</sup><sub>-1.43</sub></td><td>12.54<sup>+2.88</sup><sub>-1.53</sub></td></tr><tr><td>B19</td><td>SPT (B)</td><td>Fiducial (a), (b)</td><td><i>M</i><sub>500</sub></td><td>Total</td><td>11.30<sup>+1.11</sup><sub>-1.36</sub></td><td>11.29<sup>+1.10</sup><sub>-1.36</sub></td></tr><tr><td>Bu19</td><td>SPT (C)</td><td>Fiducial (c)</td><td><i>M</i><sub>500</sub></td><td><b>Total</b></td><td>9.45<sup>+1.16</sup><sub>-1.10</sub></td><td>9.42<sup>+1.18</sup><sub>-1.09</sub></td></tr><tr><td>Salvati22</td><td>Planck+SPT</td><td>Fiducial (d)</td><td><i>M</i><sub>500,fixed</sub></td><td>Total</td><td>11.16<sup>+1.52</sup><sub>-1.14</sub></td><td>9.51<sup>+1.81</sup><sub>-1.11</sub></td></tr><tr><td>B15</td><td>SPT (A)</td><td>Planck2014 (e)</td><td><i>M</i><sub>500</sub></td><td>Total</td><td>14.67 ± 2.29</td><td>14.67 ± 2.30</td></tr><tr><td>Planck16</td><td>Planck</td><td>Planck2016 (f)</td><td><i>M</i><sub>500</sub></td><td>HE</td><td>7.54</td><td>6.83</td></tr><tr><td>Melin21</td><td>Planck+SPT</td><td>Planck2020 (g)</td><td><i>M</i><sub>500</sub></td><td>HE</td><td>8.29<sup>+0.23</sup><sub>-0.23</sub></td><td>7.72<sup>+0.23</sup><sub>-0.23</sub></td></tr><tr><td>Salvati22</td><td>Planck+SPT</td><td>νΛCDM (h)</td><td><i>M</i><sub>500,free</sub></td><td>Total</td><td>10.35<sup>+2.72</sup><sub>-2.63</sub></td><td>9.24<sup>+2.44</sup><sub>-2.40</sub></td></tr><tr><td>Tarrío19</td><td>Planck+ROSAT</td><td>Fiducial (a)</td><td><i>M</i><sub>500</sub></td><td>HE</td><td>7.82<sup>+0.90</sup><sub>-0.96</sub></td><td>8.61<sup>+0.64</sup><sub>-0.67</sub></td></tr><tr><td>Piffaretti11</td><td>ROSAT</td><td>Fiducial (a)</td><td><i>M</i><sub>500</sub></td><td>HE</td><td>6.13</td><td>6.88</td></tr><tr><td>M10</td><td>ROSAT+Chandra</td><td>Fiducial (a)</td><td><i>M</i><sub>500</sub></td><td><b>HE</b></td><td>12.7 ± 2.5</td><td>10.3 ± 1.4</td></tr><tr><td>Z08</td><td>XMM</td><td>Fiducial (a)</td><td><i>M</i><sub>500</sub><sup>a</sup></td><td>HE</td><td>8.43 ± 2.48</td><td>...</td></tr><tr><td>Z08</td><td>XMM</td><td>Fiducial (a)</td><td><i>M</i><sub>500</sub><sup>b</sup></td><td>HE</td><td>7.66 ± 2.20</td><td>...</td></tr><tr><td>Fox22</td><td>HST</td><td>Fiducial (a)</td><td><i>M</i><sub>500,sl</sub><sup>c</sup></td><td><b>Total</b></td><td>7.54<sup>+0.33</sup><sub>-0.32</sub></td><td>...</td></tr><tr><td>CoMaLitV</td><td>ROSAT</td><td>Fiducial (a)</td><td><i>M</i><sub>500,X,wl</sub><sup>d</sup></td><td><b>Total</b></td><td>7.41 ± 0.59</td><td>7.96 ± 0.67</td></tr><tr><td>Klein19</td><td>WFI (E)</td><td>Fiducial (a)</td><td><i>M</i><sub>500</sub><sup>e</sup></td><td>Total</td><td>5.13<sup>+1.94</sup><sub>-1.69</sub></td><td>...</td></tr><tr><td>Klein19</td><td>WFI (E)</td><td>Fiducial (a)</td><td><i>M</i><sub>500</sub><sup>f</sup></td><td>Total</td><td>4.08<sup>+1.67</sup><sub>-1.31</sub></td><td>...</td></tr><tr><td>This work</td><td>SPT</td><td>Fiducial (a)</td><td><i>M</i><sub>500,UPP</sub><sup>g</sup></td><td><b>HE</b></td><td>8.36 ± 0.05</td><td>7.70 ± 0.05</td></tr><tr><td>Bu19</td><td>XMM</td><td>Fiducial (a)</td><td><i>M</i><sub>gas,500</sub></td><td><b>Gas</b></td><td>1.67<sup>+0.08</sup><sub>-0.08</sub></td><td>0.97<sup>+0.21</sup><sub>-0.21</sub></td></tr><tr><td>M10</td><td>ROSAT+Chandra</td><td>Fiducial (a)</td><td><i>M</i><sub>gas,500</sub></td><td><b>Gas</b></td><td>1.45 ± 0.25</td><td>1.18 ± 0.13</td></tr><tr><td>Z08</td><td>XMM</td><td>Fiducial (a)</td><td><i>M</i><sub>gas,500</sub></td><td><b>Gas</b></td><td>0.89 ± 0.09</td><td>...</td></tr></tbody></table>

**Note.** Cosmologies (third column): (a) Lambda cold dark matter (ΛCDM) with $\Omega_m = 0.3$, $\Omega_\Lambda = 0.7$, and $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$. (b) Added stipulation that $\sigma_8 = 0.8$. (c) ΛCDM with $\Omega_m = 0.3$, $\Omega_\Lambda = 0.7$, and $H_0 = 67.74$ km s$^{-1}$ Mpc$^{-1}$ and $\sigma_8 = 0.82$. (d) Flat $\nu\Lambda$CDM with $\Omega_m = 0.3$, $\sigma_8 = 0.8$, and $(1-b)_{SZ} = 0.58$. (e) Planck2014 cosmology (Planck Collaboration et al. 2014): $H_0 = 67.3$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_m = 0.315$, $\Omega_\Lambda = 1 - \Omega_m$, and $\sigma_8 = 0.84$; (f) Planck2016 cosmology (Planck Collaboration et al. 2016b): $H_0 = 67.5$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_m = 0.312$, $\Omega_\Lambda = 0.688$, and $\sigma_8 = 0.815$; (g) Planck2020 cosmology (Planck Collaboration et al. 2020): $H_0 = 67.4$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_m = 0.315$, $\Omega_\Lambda = 0.685$, and $\sigma_8 = 0.811$. (h) $\Omega_m = 0.29$, $\Omega_\Lambda = 0.71$, $H_0 = 61.3$ km s$^{-1}$ Mpc$^{-1}$, and $\sigma_8 = 0.76$. Facilities (second column): (A) additional facilities are used for cluster confirmation in the catalog; X-ray data from ROSAT has a small weight in cluster masses via the fitted $Y_X-M$ relation. (B) As with (A), but Chandra is used rather than ROSAT data; data from Magellan and the Hubble Space Telescope (HST) are also used for weak-lensing measurements. (C) External facilities are as in H16, except that XMM data are used instead of Chandra. (D) Richness-based weak-lensing calibration from the Dark Energy Survey (McClintock et al. 2019). (E) The Wide-Field Imager at the 2.2 m MPG/ESO telescope. Citation keys: B15: Bleem et al. (2015); B19: Bocquet et al. (2019); Bu19: Bulbul et al. (2019); H21: Hilton et al. (2021); Melin21: Melin et al. (2021); Salvati22: Salvati et al. (2022); Planck16: Planck Collaboration et al. (2016c); Tarrío19: Tarrío et al. (2019); Piffaretti11: Piffaretti et al. (2011); M10: Mantz et al. (2010); Z08: Zhang et al. (2008); Fox22: Fox et al. (2022); CoMaLitV: Sereno & Ettori (2017); and Klein19: Klein et al. (2019). Boldface indicates type used in the calculation of a hydrostatic mass bias (Figure 10).

a HE mass profile from $n_e.T_e$ profiles.

b Assumes an $L_X-M$ relation derived from the work's sample.

c $sl$ = strong lensing.

d Masses come from MCXC estimates (Piffaretti et al. 2011) and have been scaled using a weak-lensing calibration.

e Background galaxies selected by detection optimization with distance and purity cuts.

f Taking the conservative background selection without a prior on the concentration.

g Assumes the UPP shape is fixed to that in Arnaud et al. (2010).
calibrated $L_X-M$ relation (Sereno & Ettori 2017). The distinction between the green and teal groups is the assumed cosmology, where we deem the different Planck cosmologies to be sufficiently close to each other for comparison. For the fiducial cosmologies we exclude the uncorrected (“uncorr”) mass estimates from Hilton et al. (2021), which do not account for up-scatter in mass estimates with respect to the underlying cluster mass function, though it should be evident that a hydrostatic mass bias derived with that value will be lower than the bias value obtained with the corrected mass estimate. Additionally, we apply the small correction based on $h$, $(0.6774/0.7)^{-1}$, to the masses from Bulbul et al. (2019) for comparison with similar mass estimates.

## ORCID iDs

Charles E. Romero † https://orcid.org/0000-0001-5725-0359

Massimo Gaspari † https://orcid.org/0000-0003-2754-9258

Gerrit Schellenberger † https://orcid.org/0000-0002-4962-0740

Bradford A. Benson † https://orcid.org/0000-0002-5108-6823

Lindsey E. Bleem † https://orcid.org/0000-0001-7665-5079

Esra Bulbul † https://orcid.org/0000-0002-7619-5399

Matthias Klein † https://orcid.org/0000-0002-8248-4488

Ralph Kraft † https://orcid.org/0000-0002-0765-0511

Paul Nulsen † https://orcid.org/0000-0003-0297-4493

Christian L. Reichardt † https://orcid.org/0000-0003-2226-9169

Taweewat Somboonpanyakul † https://orcid.org/0000-0003-3521-3631

Yuanyuan Su † https://orcid.org/0000-0002-3886-1258

## References

Allen, S. W., Rapetti, D. A., Schmidt, R. W., et al. 2008, MNRAS, 383, 879

Applegate, D. E., von der Linden, A., Kelly, P. L., et al. 2014, MNRAS, 439, 48

Arévalo, P., Churazov, E., Zhuravleva, I., Hernández-Monteagudo, C., & Revnivtsev, M. 2012, MNRAS, 426, 1793

Arévalo, P., Churazov, E., Zhuravleva, I., Forman, W. R., & Jones, C. 2016, ApJ, 818, 14

Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92

Astropy Collaboration, Price-Whelan, A. M., & Lim, P. L. 2022, ApJ, 935, 167

Astropy Collaboration, Robitaille, T. P., & Tollerud, E. J. 2013, A&A, 558, A33

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

Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112

Bleem, L. E., Crawford, T. M., Ansarinejad, B., et al. 2022, ApJS, 258, 36

Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27

Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55

Botteon, A., Gastaldello, F., & Brunetti, G. 2018, MNRAS, 476, 5591

Bulbul, E., Chiu, I. N., Mohr, J. J., et al. 2019, ApJ, 871, 50

Chadayammuri, U., ZuHone, J., Nulsen, P., et al. 2022, MNRAS, 509, 1201

Churazov, E., Vikhlinin, A., Zhuravleva, I., et al. 2012, MNRAS, 421, 1123

de Haan, T., Benson, B. A., Bleem, L. E., et al. 2016, ApJ, 832, 95

Dupourqué, S., Clerc, N., Pointecouteau, E., et al. 2023, A&A, 673, A91

Eckert, D., Gaspari, M., Vazza, F., et al. 2017, ApJL, 843, L29

Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40

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

Fox, C., Mahler, G., Sharon, K., & Remolina González, J. D. 2022, ApJ, 928, 87

Gaspari, M., Brighenti, F., & Ruszkowski, M. 2013, AN, 334, 394

Gaspari, M., & Churazov, E. 2013, A&A, 559, A78

Gaspari, M., Churazov, E., Nagai, D., Lau, E. T., & Zhuravleva, I. 2014, A&A, 569, A67

Gaspari, M., Tombesi, F., & Cappi, M. 2020, NatAs, 4, 10

Gianfagna, G., Rasia, E., Cui, W., et al. 2023, MNRAS, 518, 4238

Hernández-Lang, D., Klein, M., Mohr, J. J., et al. 2023, MNRAS, 525, 24

Hilton, M., Sifón, C., Naess, S., et al. 2021, ApJS, 253, 3

Hincks, A. D., Acquaviva, V., Ade, P. A. R., et al. 2010, ApJS, 191, 423

Hitomi Collaboration, Aharonian, F., Akamatsu, H., et al. 2016, Natur, 535, 117

Hoekstra, H., Herbonnet, R., Muzzin, A., et al. 2015, MNRAS, 449, 685

Hofmann, F., Sanders, J. S., Nandra, K., Clerc, N., & Gaspari, M. 2016, A&A, 585, A130

Hudson, D. S., Mittal, R., Reiprich, T. H., et al. 2010, A&A, 513, A37

Hurier, G., & Angulo, R. E. 2018, A&A, 610, L4

Jennings, F., & Davé, R. 2023, MNRAS, 526, 1367

Kale, R., Parekh, V., Rahaman, M., et al. 2022, MNRAS, 514, 5969

Kale, R., Shende, K. M., & Parekh, V. 2019, MNRAS, 486, L80

Khatri, R., & Gaspari, M. 2016, MNRAS, 463, 655

Klein, M., Israel, H., Nagarajan, A., et al. 2019, MNRAS, 488, 1704

Koch, E. W., Rosolowsky, E. W., Boyden, R. D., et al. 2019, AJ, 158, 1

Krause, E., Pierpaoli, E., Dolag, K., & Borgani, S. 2012, MNRAS, 419, 1766

Lau, E. T., Gaspari, M., Nagai, D., & Coppi, P. 2017, ApJ, 849, 54

Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129

Lovisari, L., Forman, W. R., Jones, C., et al. 2017, ApJ, 846, 51

Mahdavi, A., Hoekstra, H., Babul, A., et al. 2013, ApJ, 767, 116

Mantz, A., Allen, S. W., Ebeling, H., Rapetti, D., & Drlica-Wagner, A. 2010, MNRAS, 406, 1773

Mantz, A. B., Allen, S. W., Morris, R. G., & Schmidt, R. W. 2016, MNRAS, 456, 4020

Martino, R., Mazzotta, P., Bourdin, H., et al. 2014, MNRAS, 443, 2342

Maughan, B. J., Giles, P. A., Rines, K. J., et al. 2016, MNRAS, 461, 4182

McClintock, T., Varga, T. N., Gruen, D., et al. 2019, MNRAS, 482, 1352

Melin, J. B., Bartlett, J. G., Tarrío, P., & Pratt, G. W. 2021, A&A, 647, A106

Menanteau, F., González, J., Juin, J., et al. 2010, ApJ, 723, 1523

Mohapatra, R., Federrath, C., & Sharma, P. 2021, MNRAS, 500, 5072

Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1

Nelson, K., Lau, E. T., & Nagai, D. 2014, ApJ, 792, 25

Parekh, V., Laganá, T. F., & Kale, R. 2021, MNRAS, 504, 610

Peacock, J. A. 1999, Cosmological physics (Cambridge: Cambridge Univ. Press), 1999

Piffaretti, R., Arnaud, M., Pratt, G. W., Pointecouteau, E., & Melin, J. B. 2011, A&A, 534, A109

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A16

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016a, A&A, 594, A24

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016b, A&A, 594, A13

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016c, A&A, 594, A27

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6

Planelles, S., Borgani, S., Dolag, K., et al. 2013, MNRAS, 431, 1487

Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, SSRv, 215, 25

Rasia, E., Meneghetti, M., Martino, R., et al. 2012, NJPh, 14, 055018

Romero, C. E., Gaspari, M., Schellenberger, G., et al. 2023, ApJ, 951, 41

Ruppin, F., Sembolini, F., De Petris, M., et al. 2019, A&A, 631, A21

Russell, H. R., McNamara, B. R., Sanders, J. S., et al. 2012, MNRAS, 423, 236

Russell, H. R., Nulsen, P. E. J., Caprioli, D., et al. 2022, MNRAS, 514, 1477

Salvati, L., Saro, A., Bocquet, S., et al. 2022, ApJ, 934, 129

Sarazin, C. L. 1988, X-ray emission from clusters of galaxies (Cambridge Astrophysics Series (Cambridge: Cambridge Univ. Press))

Sereno, M., & Ettori, S. 2017, MNRAS, 468, 3322

Siegel, S. R., Sayers, J., Mahdavi, A., et al. 2018, ApJ, 861, 71

Simonte, M., Vazza, F., Brighenti, F., et al. 2022, A&A, 658, A149

Snowden, S. L., Mushotzky, R. F., Kuntz, K. D., & Davis, D. S. 2008, A&A, 478, 615

Sunyaev, R. A., & Zel'dovich, Y. B. 1972, CoASP, 4, 173

Tarrio, P., Melin, J. B., & Arnaud, M. 2019, A&A, 626, A7

The Astropy Collaboration 2018, astropy v3.0.5: a core python package for astronomy, 3.0.5, Zenodo, doi:10.5281/zenodo.1461536

Wicker, R., Douspis, M., Salvati, L., & Aghanim, N. 2023, A&A, 674, A48

Wik, D. R., Sarazin, C. L., Ricker, P. M., & Randall, S. W. 2008, ApJ, 680, 17

Wittor, D., & Gaspari, M. 2020, MNRAS, 498, 4983

Wittor, D., & Gaspari, M. 2023, MNRAS, 521, L79

Zhang, X., Simionescu, A., Gastaldello, F., et al. 2023, A&A, 672, A42

Zhang, Y. Y., Finoguenov, A., Böhringer, H., et al. 2008, A&A, 482, 451

Zhuravleva, I., Chen, M. C., Churazov, E., et al. 2023, MNRAS, 520, 5157

Zhuravleva, I., Churazov, E., Arévalo, P., et al. 2015, MNRAS, 450, 4184

Zhuravleva, I., Churazov, E., Kravtsov, A., & Sunyaev, R. 2012, MNRAS, 422, 2712

Zhuravleva, I., Churazov, E. M., Schekochihin, A. A., et al. 2014, ApJL, 788, L13