# Forecasting constraints from surface brightness fluctuations in galaxy
clusters

Charles E. Romero E-mail<EMAIL_ADDRESS>Center for
Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138,
USA

###### Abstract

Studies of surface brightness (SB) fluctuations in the intracluster medium
(ICM) present an indirect estimate of turbulent pressure support and
associated Mach numbers. While high resolution X-ray spectroscopy offer means
to directly constrain line of sight gas motions, including those due to
turbulence, such observations are relatively expensive and will be limited to
nearby, bright clusters. In this respect, SB fluctuations are the most
economical means to constrain turbulent motions at large cluster radii across
a range of redshifts and masses.

To forecast what current and future X-ray and SZ facilities may achieve in SB
fluctuation studies, I review and synthesize matters of accuracy and precision
with respect to calculating power spectra of SB fluctuations, from which
turbulent properties are derived. Balance concerns of power spectrum accuracy
and precision across a range of spatial scales, I propose the use of three
annuli with: [0,0.4] $R_{500}$, [0.4,1] $R_{500}$, and [1,1.5] $R_{500}$.
Adopting these three regions, I calculate the uncertainty in the hydrostatic
mass bias, $\sigma_{b_{\mathcal{M}}}$, can be achieved for various instruments
in several scenarios. I find that Lynx and AtLAST are competitive in their
constraints at $R_{500}$, while AtLAST should perform better when constraining
$\sigma_{b_{\mathcal{M}}}$ at $R_{200}$.

††software:  astropy (Astropy Collaboration et al., 2013; The Astropy
Collaboration, 2018; Astropy Collaboration et al., 2022), SciPy (Virtanen et
al., 2020), numpy (Harris et al., 2020), emcee (Foreman-Mackey et al., 2013),
pyproffit (Eckert et al., 2017), SOXS (ZuHone et al., 2023)

##  1 Introduction

Galaxy clusters are filled with hot, ionized gas known as the intracluster
medium (ICM). It can be probed by X-ray (via brehmsstrahlung and line
emission) and via the Sunyaev-Zel’dovich (Sunyaev & Zel’dovich, 1972, ; SZ)
effect. These two probes are complementary as X-ray emission in strongly
related to gas density and the SZ (specifically the thermal SZ, hereafter tSZ)
effect is proportional to the (thermal) electron pressure. From density and
pressure distributions we can infer gas temperature and entropy distributions.
Despite this wealth of information, key pieces are missing.

One key piece missing is velocity information (see Simionescu et al., 2019,
for a review). Gas velocities may be inferred for some substructures (e.g.
Markevitch & Vikhlinin, 2007; Sayers et al., 2013; Adam et al., 2017; van
Weeren et al., 2019), but the gas velocity away from these substructures
remains unconstrained. In particular, turbulent motions are expected to fill
the ICM volume (Simionescu et al., 2019), and turbulence will play an
important role in heating, cooling, and gas and energy transport (Rebusco et
al., 2007; Chen et al., 2019; Gaspari et al., 2020; Wang et al., 2023)

In addition to the aforementioned motivations of studying turbulence, a major
driver of studying turbulence is assessing its contribution to the non-thermal
pressure in galaxy clusters. That is, a common approach to estimating cluster
masses has been to assume hydrostatic equilibrium with solely thermal pressure
support. It is unclear how much non-thermal pressure support derives from
turbulent pressure specifically (e.g. Battaglia et al., 2012; Lau et al.,
2013; Nelson et al., 2014; Biffi et al., 2016; Angelinelli et al., 2020).
Observationally, many approaches to estimate the hydrostatic mass bias have
been taken(e.g. Maughan et al., 2016; Andreon et al., 2017; Ettori et al.,
2020; Logan et al., 2022; Wicker et al., 2023), though most credence appears
to be afforded to avenues which contrast a hydrostatic mass estimate to a
total mass estimate from either weak lensing(e.g. Applegate et al., 2014;
Hoekstra et al., 2015; Sereno, 2015; Mantz et al., 2016; Okabe & Smith, 2016;
Sereno & Ettori, 2017; Bulbul et al., 2024) or cosmic microwave background
(CMB) lensing (e.g. Hurier & Angulo, 2018; Madhavacheril et al., 2020;
Ansarinejad et al., 2024).

X-ray spectroscopy (e.g. Sanders et al., 2010; Pinto et al., 2015; Ogorzalek
et al., 2017; Sanders et al., 2020), and especially high resolution X-ray
spectroscopy (Hitomi Collaboration et al., 2016, 2018) will directly constrain
gas motions and thus place constraints on turbulent velocities. While their
direct access to gas motions is critical, these observations, for XRISM
(Kitayama et al., 2014; ZuHone et al., 2016; Ota et al., 2018) or NewAthena
(e.g. Roncarelli et al., 2018; Zhang et al., 2024), will still be expensive
(potentially prohibitively so), especially if probing out to $R_{500}$ or
beyond.

Yet another avenue to accessing gas motions (indirectly) comes from studies of
surface brightness (SB) (e.g. Churazov et al., 2012) or thermodynamic
fluctuations (Schuecker et al., 2004; Hofmann et al., 2016; Ueda et al.,
2018). Initially these studies focused on the central (brightest) regions of
galaxy clusters (e.g. Zhuravleva et al., 2015; Arévalo et al., 2016; Romero et
al., 2023), but deeper observations have allowed studies out to larger radii
(Eckert et al., 2017; Zhuravleva et al., 2019; Sanders et al., 2020; Heinrich
et al., 2024), now with a few studies produce constraints even out to
$R_{500}$ (Khatri & Gaspari, 2016; Dupourqué et al., 2023, 2024; Romero et
al., 2024). While studies of SB fluctuations are still expensive, they offer
more affordable access to constraints on turbulent motions, especially in
lower density regions of clusters, than high resolution X-ray spectroscopy.

The utility of surface brightness fluctuations in the intracluster medium
(ICM) arises due to the optically thin nature of the ICM via either the
Sunyaev-Zel’dovich (SZ) effect or bremsstrahlung radiation in the X-rays. The
surface brightness in SZ (specifically the thermal SZ, hereafter tSZ) is
characterized by the Compton y parameter: $y\propto\int P_{e}dz$, where
$P_{e}$ is the thermal electron pressure and $dz$ indicates integration along
the line of sight. The analog for X-ray is $I\propto n_{e}^{2}dz$, where
$n_{e}$ is the electron density and the (omitted) temperature dependence is
often quite negligible, especially for hot ($kT_{e}>3$ keV) ICM gas when
observing in energy bands between roughly 1 and 2 keV.

The foundation of Fourier analyses of fluctuations in the ICM is the relation
between the power spectra of 2-dimensional surface brightness fluctuations,
$P_{\text{2D}}$ and the power spectra of the corresponding 3-dimensional
thermodynamic quantities, $P_{\text{3D}}$ (Peebles, 1993; Churazov et al.,
2012):

$P_{\text{2D}}(k_{\theta})=\int
P_{\text{3D}}(\mathbf{k})|\tilde{W}|^{2}dk_{z},$ (1)

where $k_{\theta}$ is a circular wavenumber (in 2D), $\mathbf{k}$ is a
spherical wavenumber in 3D, $k_{z}$ is the wavenumber with respect to the
scales along the line of sight and $\tilde{W}$ is the Fourier transform of the
window function, $W$. A (surface brightness) fluctuation map is taken as the
residual map divided the model of ICM surface brightness to get either a
$\delta y/y$ or $\delta S/S$ map. For clarity, I write $\delta y/\bar{y}$ and
$\delta S/\bar{S}$ so as to denote that $\bar{y}$ and $\bar{S}$ are models
(e.g. fitted radial profiles) of the surface brightness.

In practice, Equation 1 is often approximated as:

$P_{\text{2D}}(k_{\theta})\approx
P_{\text{3D}}(\mathbf{k})\int|\tilde{W}|^{2}dk_{z},$ (2)

and the integral is assigned its own variable:

$N(\theta)=\int|\tilde{W}|^{2}dk_{z}.$ (3)

###  1.1 Overview

Figure 1: A flowchart which shows the interconnected nature of choices
specific to power spectra calculation. Teal boxes indicate data products; blue
boxes indicate analysis choices, and red boxes are questions to be asked for
each choice to be made. The prime (′) notation denotes that the respective map
is convolved by the instrument beam.

With respect to SB fluctuation analyses, Figure 1 highlights the
interconnectedness of data products and analysis choices. In Section 2 I
consider various biases and address their relative importance. In Section 3, I
focus on uncertainties in measured power spectra of SB fluctuations,
especially for a circular region and show where a minimum in such
uncertainties can be obtained. Section 4 synthesizes the prior two sections
and sets the stage for investigating what constraints on the turbulent
hydrostatic mass bias may be achieved with current and future X-ray and SZ
instruments (Section 5). Idealized forecasting results are presented in
Section 6 and in Section 7 I discuss limitations and caveats of these
forecasts and I consider how biases will reduce the idealized constraints.
Conclusions are presented in Section 8.

Throughout this paper, I adopt a 3D power spectrum with the following
parameterization:

$P_{\text{3D,model}}(k)=P_{0}e^{-(k_{\text{c}}/k)^{\eta_{\rm
c}}}k^{-\alpha}e^{-(k/k_{\text{dis}})^{\eta_{\rm d}}},$ (4)

where $P_{0}$ is the normalization, $k_{\text{c}}$ is a cutoff wavenumber,
$\alpha$ is the spectral index, $k_{\text{dis}}$ is a dissipation wavenumber,
$\eta_{\rm c}$ and $\eta_{\rm d}$ regulate how quickly the cutoff and
dissipation occur. In Gaspari & Churazov (2013), approximate values for
$\eta_{\rm c}$ and $\eta_{\rm d}$ are found to be 4 and 3/2, respectively.
Except for an exploration in biases, I largely neglect the dissipation term in
this study. Note that if one defines the injection scale as where the
amplitude spectrum, $A_{\text{3D}}=\sqrt{4\pi k^{3}P_{\text{3D}}}$ peaks, then
this occurs where $d\ln P_{\text{3D}}/d\ln k=-3$. In the case of the model
spectrum (and $k_{\text{dis}}\gg k_{\text{c}}$) the injection wavenumber is
$k_{\text{inj}}=\eta_{\rm c}k_{c}/(\alpha-3)^{1/\eta_{\rm c}}$. For a
Kolmogorov spectrum ($\alpha=11/3$) this would mean
$k_{\text{inj}}=(1.5/\eta_{\rm c})^{(1/\eta_{\rm c})}k_{c}$.

I adopt a concordance cosmology: $H_{0}=70$ km s-1 Mpc-1, $\Omega_{M}=0.3$,
$\Omega_{\Lambda}=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 Addressing the accuracy of power spectra of fluctuations

In this section I address the matter of accuracy by investigating several
biases. Of principal concern are biases due to the method of power spectrum
estimation and physically induced biases for a perfectly spherical cluster.
Other biases, such as those due to gas clumping, cluster geometry, and
substructure (e.g. Dupourqué et al., 2023; Heinrich et al., 2024) are better
investigated suited to investigation with hydrodynamical or cosmological
simulations (e.g. Zhuravleva et al., 2023), which are beyond the scope of this
work.

###  2.1 Choice of Estimator

In the introduction I already indicated a focus on the power spectra of
normalized fluctuation maps. There are direct approaches to estimating the
Fourier transform or power spectrum of a region within an image (e.g. Koch et
al., 2019), several of which are designed to handle arbitrary masking, such a
delta-variance methods (Ossenkopf et al., 2008; Arévalo et al., 2012). When a
map is mostly empty one can employ a structure function (e.g ZuHone et al.,
2016; Seta et al., 2023), or even if a map is filled but the distances between
cells is irregular as with Voronoi tessellation (Roncarelli et al., 2018;
Zhang et al., 2024). Although much of the work of studies of fluctuations is
done for power spectra of a flat-sky image, if one is working with a spherical
sky (e.g. a HEALpix map: Górski et al. (2005)), the spherical harmonic
transforms (or Pseudo-$C_{\ell}$s) can be calculated via packages like
PolSpice (Chon et al., 2004), NaMaster (Alonso et al., 2019), or anafast
(Zonca et al., 2019). The spherical harmonic transforms can be converted to a
flat-sky power spectrum via its approximation (e.g. Khatri & Gaspari, 2016).
It is also possible to project a HEALpix map into flat-sky, though this
projection may leave artifacts in the power spectrum (e.g. Romero et al.,
2024, where the artifacts were on scales smaller than the instrument beam and
ultimately inconsequential.)

I do not endeavour to compare the accuracy of the various methodologies in
this work. I will adopt the delta variance method as proposed in Arévalo et
al. (2012), hereafter A12, as they validated its use in arbitrarily masked
maps and showed improved accuracy over the delta variance method proposed in
Ossenkopf et al. (2008). As advised in Koch et al. (2019), it is best to
perform one’s own validation tests to see how well a method performs in a
particular scenario.

####  2.1.1 Biases within A12 Delta Variance Method

The foundation of delta-variance approaches lie in the relation between the
variance in the real domain and the Fourier domain. In the case of the Mexican
hat filter proposed in A12, they note two equations for the variance for a
filtering at a spatial frequency of $k_{r}$, denoted $V_{k_{r}}$ in their
appendix. Specifically, they write for a general $n$-dimensional case:

$\displaystyle V_{k_{r}}$ $\displaystyle=\int(F*I)^{2}d^{n}x$
$\displaystyle=\int P(k)|\hat{F}_{k_{r}}(k)|^{2}d^{n}k,$

where $I$ is the image, $F$ is the Mexican hat filter, with $F*I$ denoting the
convolution of the two terms and $x$ denoting spatial dimensions; $P(k)$ is
the power spectrum of the image, and $\hat{F}$ is the Fourier transform of $F$
and $k$ denoting spatial frequencies. As $P(k)$ is the quantity of interest,
the method in A12 makes the approximation that this term can come out of the
integral. While this makes this method relatively easy to use, the drawback is
that this approximation can induce biases. For the case of a power-law
spectrum with normalization $P_{0}$ and spectral index $\alpha$:

$P_{\text{inj}}=P_{0}k^{-\alpha},$ (5)

Arévalo et al. (2012) find that there will be a normalization bias in the
recovered spectrum $P_{\text{meas}}$:

$\frac{P_{\text{meas}}}{P_{\text{inj}}}=2^{\alpha/2}\frac{\Gamma(n/2+2-\alpha/2)}{\Gamma(n/2+2)}.$
(6)

Of course, the underlying power spectrum may not be a single power-law, and
there have been two additional expressions derived for additional Gaussian
filtering. In particular, Zhuravleva et al. (2015) derive an expression for
the influence when one subtracts Gaussian patches from an image, and Romero et
al. (2023) derive an expression for the bias imparted by a point spread
function that can be characterized as a sum of Gaussians.

Figure 2: Several biases to be aware of for the A12 method. Top: Biases due to
the spectral shape, notably the spectral turnover. Here we have adopted
$\eta_{c}=1$ and $k_{\text{dis}}\gg 1$. Middle: Biases due to the spectral
shape, but regarding an exponential cutoff at high $k$. In particular, we have
set $k_{c}=1e-5$ and vary $k_{\text{dis}}$; $\eta_{d}=1$. Bottom: an example
of a masking bias. We adopt $k_{c}$ and $k_{\rm dis}$ such that our input
spectrum is effectively a pure power law with known normalization bias denoted
by the horizontal black dashed line. The legend denotes the radius of the
circular mask; the vertical dashed lines indicate the wavenumber corresponding
to the diameter of the mask.

Speaking more generally, we see that the magnitude and shape of the bias being
induced depends on the shape of the underlying spectrum. The top and middle
panels in Figure 2 show how the biases behave for some example spectra,
notably different $k_{c}$ and $k_{\text{dis}}$ for $\eta_{c}=1$ and
$\eta_{d}=1$. Note that the bias at the dissipation end is more significant as
the assumed slope is already steep ($\alpha=11/3$); that is, from the
normalization bias, we could expect any further steeping to induce larger
biases, as happens. Conversely, at the cutoff end, the power spectrum must
first pass through its peak (slope equal to zero), where the normalization
bias is small, before the slope becomes steeply positive and incurs strong
biases.

$\frac{\tilde{P}}{P}(k_{r})=\frac{\Gamma(n/2+2-\alpha/2)}{\Gamma(n/2+2)}2^{\alpha/2}\times\\\
\left[\sum_{i=1}^{N}\sum_{j=1}^{N}c_{i}c_{j}\left(\frac{2x_{i}^{2}x_{j}^{2}+x_{i}^{2}+x_{j}^{2}}{2x_{i}^{2}x_{j}^{2}}\right)^{-(n/2+2-\alpha/2)}\right].$
(7)

The adopted mask will induce its own bias; thus termed a mask bias. This mask
bias has a spectral dependency too; therefore, to emphasize the bias is due to
the mask, I choose parameters for our power spectrum model (Equation 4) such
that I effectively model a pure power law. I then measure the power spectrum
in various circles (about the center of the image). The bottom panel of Figure
2 shows the mask biases incurred. In particular, there is decent power
recovery up until scales corresponding to the diameter of the circular masks
(vertical dashed lines), at which point the recovery falls off dramatically,
which coincides with scales where one has few samples and sample variance (see
Section 3) will dominate. I note that there is also a gradual overestimation
in power from small scales up until the diameter of the mask.

###  2.2 A note on (de)projection and its associated bias

A so-called projection bias is induced due to the approximation (Equation 2)
of Equation 1. Given that the projection dependence on $W(\theta,z)$ (in
Fourier space, $|\tilde{W}(\theta,k_{z})|$) reveals a dependence on the
projected radius, this needs to be taken into account. In particular, the
projection bias for spectra measured over a given area can be calculated as
the area-averaged projection bias.

By extension, the variable $N(\theta)$ (Equation 3) will have a corresponding
area-averaged (effective) value, deemed $N_{\text{eff}}$. Alternatively, one
can produce a map of $\sqrt{N(\theta)}$ and divide a $\delta y/\bar{y}$ or
$\delta S/\bar{S}$ by this map. In doing so, the measured spectrum, while
being measured in 2D space will correspond to the deprojected (3D) spectrum.
For clarity, I explicitly label these power (or amplitude) spectra with
$\delta y/\bar{y}/\sqrt{N}$ or $\delta S/\bar{S}/\sqrt{N}$ in the subscript.

###  2.3 Recovery of single underlying power spectrum

To perform a numerical check on the projection bias over a large region of a
cluster, I inject pressure fluctuations, $\delta P/\bar{P}$, described by a
single power spectrum ($P_{\text{3D}}$). For a given pressure profile model,
$\bar{P}$, I integrate the quantity $\bar{P}*(1+\delta P/\bar{P})$ along the
line of sight to produce a Compton $y$ map. Here, I take $\bar{P}$ to be a
spherical pressure distribution described by the A10 (Arnaud et al.,
2010)111This may also be refered to as a UPP model. The parameters are those
derived from the full sample in Arnaud et al. (2010). pressure profile for a
cluster at $z=0.3$ with $M_{500}=6\times 10^{14}M_{\odot}$. An integration of
the UPP model without fluctuations yields a $\bar{y}$ map. From the two maps I
produce a $\delta y/\bar{y}$ map, on which power spectra are measured, via FFT
(no masking; over the entire image) and via A12 out to $R_{500}$.

Figure 3: Dividing by a map of $\sqrt{N}$ uncovers a well-normalized image
with the underlying $P_{3D}$ power spectrum. The recovered spectra suffer from
distortions due to projection, but otherwise this is good.

Here, I am primarily concerned with the fidelity at large scales. As such, I
adopt the model parameterization (Equation 4) with five different values of
$k_{c}$: [5e-4, 1e-3, 2e-3, 5e-3, 1e-2] kpc-1, all with $\alpha=11/3$,
$\eta_{\rm H}=4$, $P_{0}=1e{-4}{~{}\rm kpc}^{3}$, and I ignore the dissipation
term (in practice $k_{\text{dis}}=1e3$ kpc-1). For each set of parameters, 100
realizations are created of which power spectra are calculated. The recovered
power spectra at taken as the average across these realizations for each set
of parameters.

Figure 3 shows the recovered amplitude spectra on images of $\delta
y/\bar{y}/\sqrt{N}$ using both the FFT and A12 methods. The top panel shows
the uncorrected amplitude spectra, while the bottom panel shows the corrected
spectra. The corrected spectra are in excellent agreement with the input
spectra, and only at extremely steep slopes does numerical resolution reveal a
limitation of this check.

Although the shape of spectra explored here is rather uniform, it may be of
interest to consider some recovered quantities before and after correcting for
biases. Excluding the injected spectrum with $k_{c}=5e-4$, the uncorrected
spectra yield $k_{\text{inj}}$ ($A(k_{\text{inj}}$) 4% (3%) and 70% (1%)
higher for the FFT and A12 methods, respectively. For the corrected spectra,
the FFT and A12 methods recover both $k_{\text{inj}}$ and $A(k_{\text{inj}}$
within 1%. I find that these results are nearly identical if one calculates
the spectra on $\delta y/\bar{y}$ images and subsequently divides by
$N_{\text{eff}}$.

###  2.4 Which residual map to use

The formalism presented in the introduction and checked in just the previous
section is calculated in the case that there is no beam convolution. The
canonical correction to beam convolution is to simply divide the inferred
power spectrum (after noise debiasing) by the power spectrum of the beam.
However, this correction is only approximate as it is not the map of $\delta
I/\bar{I}$ that is convolved by the beam. Rather, $I$ is smoothed by the beam,
call this $I^{\prime}$, and the residual can be labeled $\delta I^{\prime}$. I
can consider the analogous model variables $\bar{I}$ and $\bar{I}^{\prime}$.

Noting that I do not have access to $\delta I/\bar{I}$, one may ask: should I
prefer $\delta I^{\prime}/\bar{I}$ or $\delta I^{\prime}/\bar{I}^{\prime}$? An
implied question is how accurately do these fractional residuals recover the
true spectrum?

To address this question I build on the infrastructure previously used to
generate pressure fluctuations and generate $y$ and $\bar{y}$ maps. I adopt
Gaussian beams with full-width half-maxima (FWHM) of $5^{\prime\prime}$,
$10^{\prime\prime}$, and $20^{\prime\prime}$ to generate corresponding
$y^{\prime}$ and $\bar{y}^{\prime}$ maps. I adopt a single power spectrum of
fluctuations with $k_{c}=1e-3$ kpc-1 and the remaining parameters as set
previously.

As the importance of the smoothing should depend on the angular size of the
cluster relative to the beam, I explore a range of cluster sizes.
Specifically, I choose to sample a grid of clusters with masses
$M_{500}=[2,4,6]\times 10^{14}{\rm M}_{\odot}$ at redshifts
$z=[0.1,0.3,0.5,0.7,0.9]$. For the assumed cosmology, this translates to
$1.6^{\prime}<R_{500}<12.4^{\prime}$.

Now, I simply need to compare the power spectra of these images. However, as
indicated in Section 2.1, measuring a power spectrum is not necessarily
trivial. There is also the matter of the region within which to measure the
power spectrum, discussed further in Section 4. As will be motivated in
Section 4, I calculate the power spectra in three regions:
$\theta<0.4R_{500}$, $0.4R_{500}<\theta<1.0R_{500}$, and
$1.0R_{500}<\theta<1.5R_{500}$. I employ the delta-variance method detailed in
Arévalo et al. (2012) and account for the induced PSF bias explained in
Section 2.1.1.

Figure 4 shows the recovered power spectra of either the $\delta
y^{\prime}/\bar{y}^{\prime}$ image relative to the power spectra of $\delta
y/\bar{y}$ (left column) or the power spectra of $\delta y^{\prime}/\bar{y}$
relative to $\delta y^{\prime}/\bar{y}^{\prime}$ (right column). To reiterate,
the (standard) PSF corrections (discussed in Section 2.1.1) are already
accounted for in $P_{\delta y^{\prime}/\bar{y}}$ and $P_{\delta
y^{\prime}/\bar{y}^{\prime}}$ and thus the motivation for this investigation.
The column on the right indicates that there is not a significant difference
between taking power spectra of $\delta y^{\prime}/\bar{y}^{\prime}$ versus
$\delta y^{\prime}/\bar{y}$.

Figure 4: The left column compares the recovered $P_{\delta
y^{\prime}/\bar{y}^{\prime}}$ relative to $P_{\delta y/\bar{y}}$, while the
right column compares $P_{\delta y^{\prime}/\bar{y}}$ relative to $P_{\delta
y/\bar{y}}$. Both $P_{\delta y^{\prime}/\bar{y}^{\prime}}$ and $P_{\delta
y^{\prime}/\bar{y}}$ have been corrected for canonical PSF bias via Equation
7. The shapes of of the relative curves in the two columns are very similar.
For each region, the darkest shade corresponds to smoothing by a beam with
FWHM = $5^{\prime\prime}$, and the lightest shade as a smoothing with FWHM =
$20^{\prime\prime}$.

I define biases, $b_{y^{\prime}/\bar{y}}$ as $1-P_{\delta
y^{\prime}/\bar{y}}/P_{\delta y/\bar{y}}$ if $P_{\delta
y^{\prime}/\bar{y}}/P_{\delta y/\bar{y}}>1$ otherwise as $1-P_{\delta
y/\bar{y}}/P_{\delta y^{\prime}/\bar{y}}$. I define
$b_{y^{\prime}/\bar{y}^{\prime}}$ analogously. Secondarily, I calculate the
ratio between the appropriate FWHM and $\theta_{500}$, the angular size of
$R_{500}$. I find that the maxima of both $b_{y^{\prime}/\bar{y}}$ and
$b_{y^{\prime}/\bar{y}^{\prime}}$ are well approximated as
$\text{FWHM}/\theta_{500}$.

####  2.4.1 Remarks on smoothing

I find no clear preference between the use of $\delta y^{\prime}/\bar{y}$ and
$\delta y^{\prime}/\bar{y}^{\prime}$ images. It is possible that different
region selections and assumed underlying spectra change the outcome; though in
the latter case I did not discern a strong dependence on the underlying
spectrum for the range I investigated. This bias may also change for X-ray
analyses given the steeper drop in surface brightness profile. However, given
the comparatively small PSF of current and future X-ray instruments, I do not
perform an analogous investigation.

While this bias can ostensibly be corrected if the underlying spectral shape
is known, the the scaling of the bias indicates that situation (FWHM relative
to $\theta_{500}$) where the bias is greatest is where the constraints on the
the spectral shape will be worst as the scales accessible will be most
limited. This suggests that in practice this bias will be quite difficult to
combat. If one adopts target uncertainties near $\lesssim 10$%, then this
investigation suggest that fluctuations studies are best suited to cases where
$\theta_{500}\gtrsim 10$FWHM.

###  2.5 Multiple underlying spectra

The evidence that there is not a single underlying $P_{3D}$ within the ICM for
any thermodynamic quantity is substantial. Of those studies of surface
brightness fluctuations that have investigated multiple regions (or apertures
about a center) have found that the inferred $P_{3D}$ vary with region or
aperture size (e.g. Zhuravleva et al., 2015; Khatri & Gaspari, 2016; Romero et
al., 2023; Dupourqué et al., 2023). This is not only expected, but it is in
fact a requirement if turbulent non-thermal pressure increases with cluster-
centric radius (e.g Battaglia et al., 2012; Nelson et al., 2014; Biffi et al.,
2016; Angelinelli et al., 2020).

####  2.5.1 A plausible scenario

Figure 5: The resultant $A_{3D}$ spectra as measured by the A12 method
relative to the volume-weighted $A_{3D}$ per shell.

In this section I explore the noiseless effects on projected spectra when
there exists a distribution of power spectra throughout a cluster. In addition
to the above notes on biases and spectral weighting, Appendix A revisits the
formalism of a projected power spectrum in the case of multiple underlying
$P_{\text{3D}}$.

Figure 6: Left: The amplitude spectra as recovered via the A12 method with no
bias corrections applied. $A_{\text{3D,wtd}}$ refers to the pressure-weighted
average of the input spectra (within the corresponding 3D shell); shown as
dotted curves. The solid curves are the recovered spectra when assuming
observations with an instrument that has a Gaussian FWHM of
$10^{\prime\prime}$, while the dashed curves are the recovered spectra with no
PSF distortion. Right: Best fit amplitude spectra to the raw spectra as
observed with an instrument with FWHM of $10^{\prime\prime}$, with bias
correction included. Fits are performed only with measured points at scales
larger than the FWHM of $10^{\prime\prime}$ (indicated by the dashed red
line). Shaded regions indicate $1\sigma$ intervals about the best fit. In all
panels, blue, orange, and green curves correspond to Ring 1, Ring 2, and Ring
3, respectively.

I adopt a schema wherein I have ten different $P_{\text{3D}}$ with
$\alpha=11/3$ and no dissipation term. The ten power spectra occupy spherical
shells (with tapered transitions) and a peak chosen such that the non-thermal
pressure follows the radial profile found in Nelson et al. (2014); see Section
5 for a more details on the non-thermal pressure profile. As the injection
scale is expected to increase with cluster-centric radius ($r$; e.g. Dupourqué
et al., 2023; Romero et al., 2024), I assume that $l_{\rm inj}\propto r$. Once
again, I assume an A10 profile of a $z=0.3$ and $M_{500}=6\times
10^{14}M_{\odot}$ cluster. The input power spectra are shown as amplitude
spectra in Figure 5.

As was done in Section 2.4, I adopt regions (rings) with edges of
[[0.4,1.0,1.5] $R_{500}$ (see also Section 4). Spectra from maps of $\delta
y/\bar{y}/\sqrt{N}$ and $\delta y^{\prime}/\bar{y}/\sqrt{N}$ are calculated.
In the left panel of Figure 6 are raw (measured) amplitude spectra insofar as
they are the values recovered by the A12 method without any bias corrections
applied and averaged over realizations. If considered from the perspective of
not knowing the input spectra, we may hope to fit a parameterized power
spectrum model to our recovered spectra, calculate biases, and infer bias-
corrected spectra. This should be done iteratively until convergence.

It is computationally easy to calculate bias corrections for spectral shape,
but more expensive to correct for a mask bias. That is, for each fit
performed, a bias should be calculated. In practice, I use
scipy.optimize.curve_fit to fit for parameters in Equation 4, where I put
strong priors on $\alpha$, $\eta_{c}$, and $\eta_{d}$ and include shape bias
corrections with each spectral model. After each iteration of
scipy.optimize.curve_fit, I calculate a mask correction for the obtained power
spectrum. I repeat these two steps until convergence. While I am working with
noiseless synthetic data, the standard deviation of the power spectra across
realizations imparts sample variance, which is used in the fitting procedures.

Table 1: Peaks of Amplitude Spectra

Region |  | $\ell_{\text{inj}}$ (kpc) |  $A_{\text{3D}}(k_{\text{peak}})$ |  $\sigma_{\text{3D}}$
---|---|---|---|---
Ring1 | Expected | 515 | 0.24 | 0.33
Raw (10) | 359 | 0.25 | 0.35
Raw | 359 | 0.25 | 0.38
Corrected | 559 | 0.25 | 0.37
Ring2 | Expected | 746 | 0.30 | 0.40
Raw (10) | 479 | 0.34 | 0.47
Raw | 479 | 0.34 | 0.51
Corrected | 673 | 0.33 | 0.47
Ring3 | Expected | 1003 | 0.46 | 0.58
Raw (10) | 479 | 0.44 | 0.60
Raw | 479 | 0.44 | 0.65
Corrected | 890 | 0.43 | 0.58

Note. — Peaks of A${}_{\text{3D}}$.

In this example, where the cutoff scale is within the measured range, I find
that the spectral shape bias is dominant, though the mask bias is still
important as it is not unity. Given the convergence from
scipy.optimize.curve_fit, I fix the mask bias to its converged value. I
subsequently run spectral fits again, but through emcee over 10000 steps
without priors on spectral parameters. From the MCMC chains I recover
amplitude spectra shown in the right panel of Figure 6.

These (bias-corrected) fits show amplitude spectra that are much better
matched to the expected amplitude spectra than the raw spectra are matched to
the expected spectra. The enlarged uncertainty envelopes (shaded regions) at
high $k$ is due to the degeneracy between the the PSF bias correction and the
slope $\alpha$. Conversely, at low $k$, there is also a degeneracy between the
(shape) bias correction and the rate of exponential cutoff ($k_{c}$ and
$\eta_{c}$). When including measurement errors, it is likely that use of
priors will be warranted (especially for terms like $\eta_{c}$.

I investigated incorporating projection biases and found that attempting to
correct for projection biases assuming the recovered (corrected) spectrum
within a ring generally produced worse results. The problem is that a proper
calculation of the projection bias requires knowing the detailed distribution
of power spectra (cf Appendix A).

Additionally, I investigated onion peeling the power spectra and found that
this does improve the recovered spectrum. For brevity, I present the findings
of onion peeling in an investigation of power spectrum recovery with added
measurement noise in Appendix C. In the scenario here, $A_{\text{3D},\delta
y^{\prime}/\bar{y}/\sqrt{N}}$ matches $A_{\text{3D},\delta
y/\bar{y}/\sqrt{N}}$ at low to moderate $k$, reinforcing the inference from
Section 2.4 that for this cluster ($z=0.3$, $M_{500}=6\times
10^{14}\text{M}_{\odot}$), the (secondary) bias smoothing due to a
$10^{\prime\prime}$ beam will be small.

In the context of relating spectra of fluctuations to turbulent velocities,
Table 1 reports on key spectral quantities for the input spectra and recovered
spectra. As seen in Section 2.3, the raw amplitude spectra decently capture
the amplitude of the peak, but it will underestimate the injection scale (in
this scenario by $\sim 50$%). Table 1 also includes the dispersion in
fluctuations (taken as the square root of the integrated spectra).

###  2.6 Remarks

In my investigation, I have found that the most important biases for accurate
spectral recover are the spectral shape and masking biases, which are induced
by the A12 method. While exploring these biases for other methods is beyond
the scope of this paper, such biases should clearly be checked for whichever
method is adopted.

In the plausible scenario, the radially evolving $P_{\text{3D}}(r)$ also
produces a ”contamination” bias, where $P_{\text{3D}}(r>\theta)$ contaminate
spectra within an annulus out to $\theta$. This appears to be minimal, but for
sufficiently sensitive observations, applying an onion peeling approach may be
helpful (see Appendix C).

Additional biases noted in this section include the projection bias and
smoothing bias. However, I found that, for my plausible scenario, attempts to
correct for the projection bias were detrimental. Given the relatively small
impact of the smoothing bias expected given $\theta_{500}=4^{\prime}.3$ and a
Gaussian beam of $10^{\prime\prime}$, no treatment is applied here.

##  3 Uncertainties in power spectra

Uncertainties in power spectra have long been calculated assuming either an
underlying Gaussian or Poisson random field (e.g Tukey, 1956; Blackman &
Tukey, 1958; Hoyng, 1976). In Section 2.1, I adopted the A12 method of
calculating power spectra. In Arévalo et al. (2012), they provide a formula
assuming that the data correspond to a Gaussian random field. They find that
$\sigma_{P}(k)\propto P(k)/\sqrt{N}$ for $N$ total data points.

Reframing the what is presented in Arévalo et al. (2012), I find:

$\sigma_{P}(k)\simeq P(k)\sqrt{2/N_{k}},$ (8)

where $N_{k}=\Omega\pi k^{2}$ sampling elements for an sampling area $\Omega$
at a given frequency $k$. Note that this $N_{k}$ corresponds to the number of
Gaussian smoothing kernels (those used in the Mexican hat filter) fit within a
region of area $\Omega$. This scaling matches that found in Arévalo et al.
(2012) and is found to match results from previous works (see Appendix D).
Note that this scaling implies that the uncertainty in the power spectrum will
be steeper than the power spectrum by 1 for a two dimensional dataset, i.e. in
images.

It is useful to write the total measured power spectrum, $P_{\rm tot,m}$ as
the sum of the measured power spectrum due to the signal, $P_{\rm s,m}$ and
the measured noise component $P_{\rm n,m}$:

$P_{\rm tot,m}(k)=P_{\rm s,m}(k)+P_{\rm n,m}(k).$ (9)

We are interested in the calculation of the signal power spectra and its
uncertainty. We do not have direct access to $P_{\rm s,m}(k)$, but access it
through $P_{\rm s,m}(k)=P_{\rm tot,m}(k)-P_{\rm n,m}(k)$. Let us consider
$P_{\rm s,e}(k)$ to be the expected signal power spectrum. From this, we can
calculate the expected sample variance, $\sigma_{P_{\rm s,e}(k)}$. The total
uncertainty on the inferred signal spectrum is then given by:

$\sigma_{P_{\rm s,m}(k)}^{2}=\sigma_{P_{\rm n,m}(k)}^{2}+\sigma_{P_{\rm
s,e}(k)}^{2}.$ (10)

###  3.1 White noise approximation

Note that sample variance, $\sigma_{P_{\rm s,e}(k)}$, does not improve with
observational sensitivity. However, $\sigma_{P_{\rm n,m}(k)}$ will improve (by
definition) with the sensitivity of observations. I will focus on the
uncertainty due to noise, $\sigma_{P_{\rm n,m}(k)}$.

With X-ray observations, the noise is Poisson, which equates to white noise.
In SZ observations, the noise can be pink, but often at the scales of interest
one tends to be in a regime where the noise is sufficiently white (e.g. Romero
et al., 2023, 2024). For the case of white noise, $P_{n}(k)=\text{constant}$,
in a 2D image Parseval’s theorem gives us

$P_{n}(k)\pi k_{\rm max}^{2}=V,$ (11)

where $V$ is the variance in the image and $k_{\rm max}$ is the highest
frequency being sampled. In practice, the highest frequency is set by the
(binned) pixel size for X-ray observations and by the FWHM of the beam for SZ
observations.

###  3.2 Uncertainties within a circular region

Let us assume that an observational facility points at the center of a
(spherical) cluster. The variance in the image for both X-ray and SZ
facilities can be reasonably well approximates as a function of $\theta$, the
projected cluster radius. For SZ observations, let us denote the variance and
standard deviation in the observed images as $V_{y}$ and $\sigma_{y}$,
respectively. The variance in the normalized residual ($\delta y/\bar{y}$) map
is then give by:

$V_{\delta
y/\bar{y}}(\theta)=\left(\frac{\sigma_{y}(\theta)}{\bar{y}(\theta)}\right)^{2},$
(12)

In order to measure a power spectrum, one needs to do so over an area. The
desired quantity is then the area-averaged variance. Let us consider the case
of a circular region about the cluster center. The variance within a circle of
radius $\theta$ can be expressed as:

$\langle V_{\delta y/\bar{y}}\rangle(\theta)=\frac{\int_{0}^{\theta}V_{\delta
y/\bar{y}}2\pi\theta^{\prime}d\theta^{\prime}}{\pi\theta^{2}}.$ (13)

Consequently, from equations 8 and 11 the uncertainty due to noise in a
residual spectrum will be given by:

$\sigma_{P_{\text{n,m},\delta y/\bar{y}}}(k,\theta)=\frac{\langle V_{\delta
y/\bar{y}}\rangle(\theta)}{\pi k_{\rm
max}^{2}}\sqrt{\frac{2}{N_{k}(\theta)}}.$ (14)

If we note that the dependence on $k$ simply governs a normalization factor
for $\sigma_{P_{\text{n,m},\delta y/\bar{y}}}(k,\theta)$, at a given $\theta$,
we can arbitrarily choose a $k$ to investigate the behavior over $\theta$.

In the case of X-ray observations, the variance is given simply by the
expected counts, such that for a surface brightness rate model ($\bar{S}_{\rm
tot}$ with units counts/s/area), one can write:

$V_{\text{counts}}=\bar{S}_{\rm tot}E\Omega_{\text{pix}},$ (15)

where $E$ is the exposure map (in seconds), and $\Omega_{\text{pix}}$ is the
pixel area. The variance in the rate map is given by
$V_{\text{counts}}/(E\Omega_{\text{pix}})^{2}$ and the variance in the $\delta
S/\bar{S}$ image is given by:

$V_{\delta
S/\bar{S}}=\frac{\bar{S}_{tot}}{E\Omega_{\text{pix}}\bar{S}_{ICM}^{2}}.$ (16)

Having arrived at an expression for $V_{\delta S/\bar{S}}$, the calculation of
$\langle V_{\delta S/\bar{S}}\rangle$ and $\sigma_{P_{\text{n,m},\delta
S/\bar{S}}}(\theta)$ follow as in Equations 13 and 14 respectively.

At this point it is instructive to consider an example cluster with
implications that can largely be generalized under the assumption of self-
similarity. I take again a cluster at $z=0.3$ with $M_{500}=6\times
10^{14}M_{\odot}$ and the A10 pressure profile. I also explore various other
pressure profiles (sets of gNFW parameters) presented in Arnaud et al. (2010)
and McDonald et al. (2014), which have corresponding profiles in grey in
Figure 7.222Because I chose $z=0.3$, I adopt the parameters in the low-$z$
sample of McDonald et al. (2014).

Figure 7: Top: surface brightness profiles for an example cluster; grey lines
are gNFW pressure profiles with parameter sets other than that of the A10.
Bottom: the variance within a circle of radius $r$ (dashed lines) as expressed
in Equations 13 and 16; for the SZ profiles ($\langle V_{\delta
y/\bar{y}}\rangle$ and $\sigma_{\langle V_{\delta y/\bar{y}}\rangle}$) have
coloring matching that in the top panel. Each star marks the minimum of the
profile upon which it is marked.

For the X-ray $\beta$-model:

$I=\frac{I_{0}}{(1+(r/r_{c})^{2})^{x}},$ (17)

I adopt $\beta=0.66$, and $r_{c}=0.1R_{500}$ and normalized to be consistent
with the luminosity and temperature scaling relation presented in Bulbul et
al. (2019). These parameters give a qualitatively good agreement with surface
brightness profiles presented in, for example, Bartalucci et al. (2023);
Lyskova et al. (2023); Romero et al. (2024). A background rate is taken that
is similar to that seen in Romero et al. (2023, 2024).

Adopting an XMM-like exposure (in the EPIC cameras) and a uniform sensitivity
map for SZ ($V_{y}(\theta)=\text{constant}$), the bottom panel of Figure 7
show the variance profiles and associated $\sigma_{P_{\text{n,m}}}(k,\theta)$
for SZ and X-ray for arbitrary $k$ and exposure times, as neither $k$ nor
exposure times affect the shape of the profiles.

###  3.3 Incorporating biases into uncertainties

From Section 2, we’ll recall that the measured power spectra can suffer from
multiple biases, which we can aggregate into a bias for each component. That
is, $B_{\rm s}$ and $B_{\rm n}$ can account for the total biases on the
measured signal and noise power spectra such that $P_{\rm s,m}(k)=P_{\rm
s,t}(k)B_{\rm s}$ and $P_{\rm n,m}(k)=P_{\rm n,t}(k)B_{\rm s}$, for true
(unbiased) spectra $P_{\rm s,t}(k)$ and $P_{\rm n,t}(k)$. Rewriting Equation
9, we have:

$P_{\rm s,t}(k)=\frac{P_{\rm tot,m}(k)-P_{\rm n,t}(k)B_{\rm n}}{B_{\rm s}}.$
(18)

For $B_{\rm n}\approx 1$, one can infer a bias-corrected uncertainty on the
signal power spectrum as $\sigma_{P_{\rm s,t}(k)}=\sigma_{P_{\rm
s,m}(k)}/B_{\rm s}$.

##  4 Synthesizing accuracy and uncertainty concerns

The previous two sections have addressed matters of accuracy and uncertainty
in recovered power spectra of surface brightness fluctuations. Here, I
synthesize the previous two sections with respect to choice of regions. In
particular, I focus on the case circular (or elliptical) annuli about the
cluster center, where the goal(s) align with the determination of properties
of turbulence as a function of radius at to (at least) $R_{500}$. Moreover, I
want the choice of regions to accommodate SB fluctuation analysis via multiple
instruments, so that analyses, at least with respect to region choices, are
more homogeneous.

When determining properties of turbulence as a function of radius, there is a
trade-off between radial resolution (how many annular regions are used) and
the dynamic range that is probed with those regions, due to masking bias and
sample variance. Moreover, recovering the injection scale suggests avoiding
thin annuli as injection scales are likely to occur at scales where sample
variance is important (e.g. Dupourqué et al., 2023; Heinrich et al., 2024). If
we want to infer some radial trend, out to $R_{500}$, we want a minimum of two
regions within $R_{500}$. Adopting two regions will be most accommodating to
homogeneous analyses with instruments of various spatial resolutions.

Section 3 strongly motivates choosing an inner region to be a circle of radius
$R_{1}$, such that $0.1R_{500}<R_{1}<0.3R_{500}$ as that range yields a
minimum in uncertainties in a power spectrum of SZ fluctuations, while
$R_{1}\approx 0.1R_{500}$ yields a minimum in the uncertainties of a
corresponding X-ray fluctuation power spectrum. Conversely, fluctuations
within $0.1R_{500}$ may be susceptible to SB model choices and thus omitted
(e.g. Dupourqué et al., 2024), though given the nature of area-averaging, such
concerns can be mitigated with a larger central region (e.g. Romero et al.,
2024).

While $R_{1}=0.3R_{500}$ may be indicated (and $R_{2}=R_{500}$) for the above
motivations, I propose adopting $R_{1}=0.4R_{500}$. The increase in
uncertainty (relative to $R_{1}=0.3R_{500}$), especially for SZ measurements,
is not particularly consequential. In comparison, $R_{1}=0.4R_{500}$
accommodates homogeneous analysis with a wider range of instrument resolutions
and attenuates the impact of fluctuations in the central most regions. As an
added bonus, $0.4R_{500}$ is a decent approximation to $R_{2500}$ (e.g. Arnaud
et al. (2010) takes $R_{2500}\sim 0.44R_{500}$).

Finally, one may consider measurements beyond $R_{500}$, for which I propose
$R_{200}\sim 1.5R_{500}$ (e.g. Romero et al., 2020) as it satisfies the above
concerns of biases and dynamic range recoverable, and of course $R_{200}$ is a
useful point of comparison.

##  5 Forecasting constraints on the hydrostatic mass bias: approach

In the remainder of the paper, I will turn my attention to the what
constraints on the (turbulent) hydrostatic mass bias we may expect with some
current and proposed future facilities. In particular, I wish to determine the
redshifts and masses of clusters for which constraints of a given level can be
achieved. To this end, I explore clusters between $0.05\leq z\leq 1.65$ and
$0.08<(M_{500}/10^{14}M_{\odot}<12.6$. I adopt the annular regions suggested
in Section 4. Because the dominant biases when using the A12 method depend on
the spectral shape, but I do not wish to assume a spectral shape nor explore a
large parameter space, I calculate power spectral values assuming perfect
accounting of biases. Of more concern are the associated uncertainties
(detailed in Appendix A.1), where I do incorporate a correction to the PSF
bias in the spectral uncertainties.

###  5.1 Formulation of Hydrostatic Bias

To relate the amplitude spectrum to the turbulent Mach number, $\mathcal{M}$,
I use the (adapted) relation from Gaspari & Churazov (2013):

$\mathcal{M}=4A_{\rho}\left(\frac{l_{\rm
inj}}{0.4R_{500}}\right)^{-\alpha_{H}}=2.4A_{P}\left(\frac{l_{\rm
inj}}{0.4R_{500}}\right)^{-\alpha_{H}},$ (19)

where $A_{\rho}$ and $A_{P}$ are the peaks of the density and pressure
amplitude spectra, respectively. The injection scale is denoted as $l_{\rm
inj}$ and $\alpha_{H}\simeq 0.25$ is an empirically derived scaling for the
hydrodynamical cascade of fluctuations. This scaling relation has been found
to be in agreement with with the relation from Zhuravleva et al. (2012); see,
e.g. Romero et al. (2023); Dupourqué et al. (2024).

To determine the input turbulent Mach numbers I employ its relation to the
non-thermal pressure fraction due to turbulent motions:

$\frac{P_{\rm NT}}{P_{\rm th}}=(\gamma/3)\mathcal{M}^{2},$ (20)

where I take $\gamma=5/3$ to be the adiabatic index. With this relation, I use
the radial profile of non-thermal pressure as given in Nelson et al. (2014):

$\frac{P_{\rm NT}}{P_{\rm th}+P_{\rm
NT}}=1-A\left(1+exp\left[-\left(\frac{r/r_{200m}}{B}\right)^{\gamma_{P}}\right]\right),$
(21)

where they found values of $A=0.452\pm 0.001$, $B=0.841\pm 0.008$, and
$\gamma_{P}=1.628\pm 0.019$. While Nelson et al. (2014) and Battaglia et al.
(2012) consider the non-thermal pressure to be completely from random (i.e.
turbulent) motions, other sources of non-thermal pressure may not be
negligible (e.g. Angelinelli et al., 2020).

Figure 8: Mach profiles, as inferred from Equation 20 and non-thermal pressure
profiles in Battaglia et al. (2012) and Nelson et al. (2014).

Khatri & Gaspari (2016) provide a relation between the hydrostatic bias and
$\mathcal{M}_{\text{3D}}$ (and attach a corresponding subscript to denote the
method of calculation, $b_{\mathcal{M}}$):

$b_{\mathcal{M}}=\frac{-\gamma\mathcal{M}_{\text{3D}}^{2}}{3}\frac{d\ln
P_{\text{NT}}}{d\ln
P_{\text{th}}}\left(1+\frac{\gamma\mathcal{M}_{\text{3D}}^{2}}{3}\frac{d\ln
P_{\text{NT}}}{d\ln P_{\text{th}}}\right)^{-1},$ (22)

where $\gamma$ is the adiabatic index, taken to be 5/3 for the ICM, and
$b_{\mathcal{M}}$ is implicitly a function of $r$. NB that as defined in
Khatri & Gaspari (2016) $b_{\mathcal{M}}\equiv M_{x}/M_{\text{tot}}-1$.
Following the recasting performed in Khatri & Gaspari (2016), I find:

$\frac{d\ln P_{\text{NT}}}{d\ln P_{\text{th}}}=\frac{d\ln P_{\text{NT}}/d\ln
r}{d\ln P_{\text{th}}/d\ln r}=1+2\frac{d\ln\mathcal{M}_{\text{3D}}/d\ln
r}{d\ln P_{\text{th}}/d\ln r}.$ (23)

I adopt the rings suggested in Section 4; those with edges at [0.4, 1.0, 1.5]
$R_{500}$. The effective radii for these rings is roughly [0.3, 0.8, 1.3]
$R_{500}$. To calculate $b_{\mathcal{M}}(R_{2500})$,
$b_{\mathcal{M}}(R_{500})$, and $b_{\mathcal{M}}(R_{200})$, I need to estimate
$\mathcal{M}(R_{2500})$, $\mathcal{M}(R_{500})$, and $\mathcal{M}(R_{200})$.
From the inner two rings, or all three rings, I fit a power law slope to
$\mathcal{M}(r)$ at the effective radii of each ring and use the fitted power
law to calculate $\mathcal{M}(R_{2500})$, $\mathcal{M}(R_{500})$, and
$\mathcal{M}(R_{200})$ (and consequently $b_{\mathcal{M}}(R_{2500})$,
$b_{\mathcal{M}}(R_{500})$, and $b_{\mathcal{M}}(R_{200})$), propagating all
uncertainties, including those from sample variance (see Appendix A.1).

Equation 23 also requires calculation of $d\ln P_{\text{th}}/d\ln r$, for
which I assume that the logarithmic slope from the A10 profile is recovered.
Relative to the uncertainties on other quantities in Equation 23, this should
not be the limiting factor, even for X-ray observation. Given that the
uncertainty for this quantity in the case SZ observations is straightforward
(see below), I include uncertainties in the SZ case. In the X-ray case, this
uncertainty is ignored.

###  5.2 Toy setup

In order to quantify the expected uncertainties on turbulent hydrostatic mass
bias, I require a SB model, a model of the noise (spatially and spectrally),
and (3) a defined region and scale being sampled for a given observation
(Section 3).

####  5.2.1 SZ modelling

To forecast all SZ observations, I again adopt the A10 pressure profile and I
assume white noise, often with a radial dependence. Mapping speeds can be
calculated from the various time sensitivity calculators provided for
MUSTANG-2333https://greenbankobservatory.org/science/gbt-observers/mustang-2/,
NIKA2444https://publicwiki.iram.es/Continuum/TimeEstimatorScriptGuideW2023,
TolTEC555https://toltec.lmtgtm.org/toltec_sensitivity_calculator, and
AtLAST666https://github.com/ukatc/AtLAST_sensitivity_calculator. WIKID is a
proposed successor to MUSTANG-2; I have produced mapping speed estimates using
MUSTANG-2 data with the projected increased detectors and improved detector
and readout noise. I assume that the instruments perform scans which are
roughly equivalent to Lissajous daisy scans of $5{\arcmin}$ (NIKA2),
$3{\arcmin}$ (TolTEC), $3{\arcmin}.5$ (MUSTANG-2), and $3{\arcmin}.5$ (WIKID),
on a single pointing.

AtLAST777https://www.atlast.uio.no/ is intended to have a camera with over a
degree field of view; currently the aim is 2 degrees. Assuming an optimal
internal linear combination (ILC) of sensitivites in ALMA-like bands, I find a
mapping speed of $1.6\times 10^{-6}$-root hour in Compton y, similar to that
found in Di Mascolo et al. (2024). For comparison, SPT-SZ (Bleem et al., 2015)
would achieve an optimal ILC sensitivity of $\sim 3\times 10^{-6}$ in Compton
$y$-arcmin, but their minimum variances actually achieve $1.4\times 10^{-6}$
in Compton $y$-arcmin. This corresponds to a degradation factor of $4.8$.
While the sensitivities of each band and angular resolution will affect this
degradation, I adopt this factor and arrive at a Compton $y$ mapping speed of
$8\times 10^{-6}$-root hour at beam scales.

SPT-3G is expected to be $\sim$ten times more sensitive than SPT-SZ (Benson et
al., 2014). For both SPT-3G and AtLAST will employ a survey strategy such that
the noise is uniform and the transfer function is effectively unity at scales
being probed. I utilize a parametric model of the transfer function fit to the
MUSTANG-2 transfer function (Romero et al., 2020), and I scale this model to
the FOV of NIKA2, TolTEC, and WIKID.

####  5.2.2 X-ray modelling

For the case of the X-ray observations, I assume an isothermal $\beta$ model
(Equation 17) with $\beta=2/3$, $r_{c}=0.1R_{500}$ and $I_{0}$ taken to
satisfy the $L_{X,cin}{–}M500{–}z$ and $T_{X,cin}{–}M500{–}z$ relations in
Bulbul et al. (2019), which is seen to be qualitatively similar to to profiles
in the literature (as mentioned in Section 3). To properly normalize the SB
profiles, I employ SOXS888https://hea-
www.cfa.harvard.edu/soxs/index.html(ZuHone et al., 2023) which calculates the
rest-frame luminosity based on temperature, and from this I derive the
appropriate SB normalization.

I then use SOXS to simulate observations, with a single pointing, in the 0.5
to 2.0 keV range of our X-ray (imaging) instruments on Lynx, NewAthena, AXIS,
and LEM. I assume a modest column density (4e18 cm${}^{-}2$) on all
observations. I apply instrumental background and foreground to one
observation per instrument to calculate a background count rate. I use
pyproffit to extract SB profiles, and fix $\beta=2/3$ and $r_{c}=0.1R_{500}$.

For XMM, I use results obtained in Romero et al. (2024) to establish a scaling
relation with self-similar scaling: $L_{X}\propto M_{500}E(z)^{2}$ that
corresponds to the soft energy band (0.5 to 1.25 keV) used in Romero et al.
(2024). As the background, I take an average of the total backgrounds in the
0.5 to 1.25 keV and assume this background is constant for our forecasting
constraints. This creates a slight mismatch, but I do not intend to focus on
the forecast results for XMM; I include them to show consistency, in a very
broad sense, between my forecasts and results which have been published.

####  5.2.3 Remaining parameters

In order to employ Equation 19, I must adopt an injection scale. However, as
the literature has not borne out a clear expectation for what the injection
scales will be across a large range of redshifts and masses. Observational and
numerical simulations tends to find injection scales of a few hundreds of kpc
(in the masses, redshifts probed; e.g. Gaspari et al., 2013; Khatri & Gaspari,
2016; Sanders et al., 2020; Romero et al., 2023; Dupourqué et al., 2023,
2024). Additionally, the injection scale is expected increase with cluster
radius. Therefore, I consider a range of injections scales:
$[0.005,0.008,0.013,0.025,0.05,0.1,0.2,0.4,0.6,0.8]\times[R_{500},1$ Mpc],
where injection scales $<0.1\times[R_{500},1$ Mpc] are included for the
purposes of cataloging the uncertainties in the recovered spectra at those
scales (see Section 7). I consider plausible injection scales to be $\geq
0.1\times R_{500}$.

Finally, to account for the noise in maps with targeted instruments, I must
assume an exposure time. Although AtLAST is likely to perform a survey (Di
Mascolo et al., 2024), the depth of that survey is yet to be determined. As
such, I explore various exposure times for it as well. In particular, I
consider three exposure times: 50, 250, and 1000 ks. While 50 ks is still not
a short exposure, in the regime of SB fluctuations, it should be considered a
short exposure. 250 and 1000 ks thus correspond to a modest and deep exposure.

##  6 Forecasting constraints on the hydrostatic mass bias: results

Figure 9: Mass thresholds for achieving $\sigma_{b_{\mathcal{M}}}=0.1$ for
MUSTANG-2 (solid lines) and XMM (dashed lines). Above these lines the
uncertainties on $b_{\mathcal{M}}$ will be less than $0.1$. The various panels
in this figure show how the constraints change based on assumptions and
exposure time. Top left: The injection scale is assumed to scale with
$R_{500}$ and all results take a fixed exposure time (on-source) of 50 ks. Top
right: same as top left, but injection scales are at specified physical
scales. Bottom left: Injection scale is held at a fixed fraction of $R_{500}$
and exposure times are varied. Bottom right: Same as bottom left except the
injection scale is held at a fixed physical scale.

As discussed in the previous section, SB profiles have been assumed for
clusters between $0.05\leq z\leq 1.65$ and
$0.08<(M_{500}/10^{14}M_{\odot}<12.6$. Accounting for expected image noise,
imaging fields of view (for a single pointing), transfer functions, and a
Gaussian PSF assumption, I can calculate the mass threshold (as a function of
redshift) above which instruments can achieve a given uncertainty in the
hydrostatic mass bias, $\sigma_{b_{\mathcal{M}}}$ or better. Given that values
for the hydrostatic mass bias, $b$, are expected to be $0.1\lesssim b\lesssim
0.3$(e.g. Hurier & Angulo, 2018), a comparable (or better) can be taken as a
target uncertainty. I take $\sigma_{b_{\mathcal{M}}}=0.1$ as my target
uncertainty (on a single cluster) as this appears achievable in the presented
approach, while uncertainties much less than this may not be possible due to
sample variance.

Figure 10: Mass thresholds (as in Figure 9 for achieving
$\sigma_{b_{\mathcal{M}}}=0.1$ when $l_{\rm inj}=0.4R_{500}$ for select
current and proposed X-ray facilities (top) and SZ facilities (bottom with 50
ks on-source time. Here, I assume imaging cameras are used. Limitations at the
low-redshift end come from $R_{500}$ exceeding the field of view. The inferred
SPT-3G results assume the end-of-survey expected sensitivity.

Figure 9 shows how the mass threshold to achieve
$\sigma_{b_{\mathcal{M}}}=0.1$ changes as a function of $l_{\text{inj}}$ and
exposure (integration) time, $T_{\rm exp}$. The improvement in constraints as
the injection scale increases can be inferred as a consequence of how the
noise varies with scale (Secton 3) and how quickly the signal varies with
scale (e.g. $\alpha\gtrsim 3$).

Figure 11: Mass thresholds for achieving $\sigma_{b_{\mathcal{M}}}=0.1$ when
$l_{\rm inj}=0.4R_{500}$ for some of the more sensitive X-ray or SZ
facilities. The grey dots are clusters in ACT, SPT, Planck, MCXC, or eRosita
catalogs and serve to contextualize the mass threshold curves. Moreover, the
blue boxes and chartreuse diamonds indicate masses and redshifts of clusters
(or elliptical galaxies) studied via SB fluctuations or spectroscopy to infer
gas velocities, but generally not out to $R_{500}$.

While I explored the consequence of assuming injections scales at fixed
physical scales, I anticipate that injection scales that scale with $R_{500}$
are likely to be the better descriptor of fluctuations in clusters. As such,
Figure 10 adopts the case that $l_{\text{inj}}=0.4R_{500}$ and
$T_{\text{exp}}=50$ ks, with the exception of SPT-3G, whose sensitivity
(noise) is the expected noise at the end of the survey.

Not surprisingly, one of the more striking results is that SZ instruments can
perform much better at higher redshift. Conversely, several SZ instruments are
quite hampered at low redshift due to their instantaneous field of view, which
hurts both the imageable field of view and most consequentially, the transfer
function. While SZ survey instruments have sufficiently large fields of view,
their resolution and sensitivity limit the constraints they can provide. I
note that results for SO and CMB-S4 are similar as SPT-3G as the expected
sensitivities are similar with SO and CMB-S4 having marginally poorer angular
resolution than SPT-3G.

Figure 11 retains forecasts of $\sigma_{b_{\mathcal{M}}}$ for just three
instruments. For context, I have added clusters detected in several catalogs.
For additional context, I have added clusters, groups, and/or BCGs for which
turbulent velocities have been measured or inferred. Most of studies of
turbulent velocities, either by spectroscopy or SB fluctuations, have focused
on the central regions of clusters, and several clusters have been studied
multiple times. That is to say that the constraints denoted by dotted curves
in Figure 11 signify a substantial improvement in the existing constraints of
fluctuations out to $R_{500}$.

##  7 Discussion

Section 6 presented constraints on the hydrostatic mass bias assuming perfect
bias corrections and that a single measurement at the injection peak was
sufficient to infer the turbulent velocities and subsequently the hydrostatic
mass bias. While decent bias corrections can be achieved, doing so requires
access to nearly a decade in scales (Section 2.5.1). This requirement, along
side the secondary smoothing bias (Section 2.4.1) suggest the need for the
FWHM of an instrument to be less than $0.1R_{500}$ for a targeted cluster.
Figure 12 adds curves denoting what fraction of $R_{500}$ corresponds to
$10^{\prime\prime}$. Additionally, constraints for the three integration times
are shown for AtLAST and Lynx. Finally, in Figure 12, I have trimmed the
studies in the literature to those with constraints (even if marginal) out to
$R_{500}$. Figure 13 replicates Figure 12, but shows constraints for
$\sigma_{b_{\mathcal{M}}}=0.1$ at $R_{200}$.

Figure 12: Similar to Figure 11 with thresholds for achieving
$\sigma_{b_{\mathcal{M}}}=0.1$ on $M_{500}$ when $l_{\rm inj}=0.4R_{500}$ by
instrument and exposure time. Here I have limited markers to those which
indicate some (often poor) constraints on turbulent velocities out to
$R_{500}$ from SB fluctuation analyses, some of which have corresponding
hydrostatic mass bias estimates.

Two prominent findings from Section 2.5.1 were that (1) assuming the injection
scale is within the scales probed, good constraints on the bias-corrected
spectral shape benefit greatly when $k_{\rm peak,raw}\lesssim k_{\rm max}/4$
and (2) the peaks of the raw amplitude spectra, $A_{\rm 3D}(k_{\rm peak})$,
and the overall density or pressure fluctuations, $\sigma_{\rm 3D}$, as
recovered in the raw spectra are not substantially biased. In the case that
$l_{\rm inj}=0.4R_{500}$, then the goal is to sufficiently recover
fluctuations at scales of $0.1R_{500}$. In this respect, resolution should
never be a limiting factor Lynx. However, for instruments such as WIKID or
AtLAST, with FWHM of $10^{\prime\prime}$, then resolution is a concern, where
the grey line marked $0.1R_{500}=10^{\prime\prime}$ suggests a a threshold
beyond which AtLAST or WIKID may not provide robust constraints.

Figure 13: Similar to Figure 12, but curves represent thresholds for achieving
$\sigma_{b_{\mathcal{M}}}=0.1$ on $M_{200}$.

On account of the injection scale, the requirement to obtain fluctuations down
to $0.1R_{500}$ may be alleviated if the injection scale is larger than
$0.4R_{500}$, which is quite plausible for Rings 2 and 3. However, there
remains the challenge of recovering down to $0.1R_{500}$ in Ring 1 (assuming
$l_{\rm inj}\simeq 0.4R_{500}$). If the secondary bias due to beam smoothing
(Section 2.4.1) can be satisfactorily addressed, then the principal remaining
concern is what sensitivity is required at $4k_{\rm peak,raw}$ to
satisfactorily recover the unbiased power spectrum. Appendix C suggest that
$\sigma_{A_{\text{3D}}}\simeq 0.05$ is sufficient.

Figure 14: Scales relative to $R_{500}$ for which either Lynx (left, green) or
AtLAST (right, blue) achieves an uncertainty in $A_{3D}$ of 0.05 or smaller
with a 1 Ms exposure. The top row is taken for Ring 1 (0-0.4$R_{500}$); the
middle row for Ring 2 (0.4-1.0 $R_{500}$); and the bottom row for Ring 3
(1.0-1.5$R_{500}$). Maximum scale at which the shading starts is $0.4R_{500}$.
Dotted curves correspond to thresholds if no PSF corrections are needed (i.e.
the PSF is infinitesimal); dashed curves correspond to thresholds if the FWHM
of the PSF were half of the nominal FWHM. Black curves denote constraints at
scales which could plausibly be at $l_{\rm inj}/4$ (see text). Grey dots are
detected clusters from surveys, as in Figure 11. The red dashed curve denotes
$0.1R_{500}=10^{\prime\prime}$.

Figure 14 shows the mass thresholds where $\sigma_{A_{\rm 3D}}(k=1/l)=0.05$
for various scales. The black curves indicate the thresholds for the (modestly
conservative) case that $l_{\rm inj,2}\gtrsim 0.6R_{500}$ in Ring 2 and
$l_{\rm inj,3}\gtrsim 0.8R_{500}$ in Ring 3.999Recall from Section 5, the
effective radii for the rings are roughly [0.3, 0.8, 1.3] $R_{500}$. That
threshold is fairly consistent across the three rings (at $M_{500}\approx
3\times 10^{14}\text{M}_{\odot}$) for AtLAST, while the threshold shows a
marked evolution across the three rings for Lynx.

Figure 14 also shows the constraints if the noise in the images were the same
and the FWHM of the PSF were half of its nominal value (dashed curves) or were
infinitesimal (dotted curves). For Lynx, the improvement (on the scales shown)
is minimal, while for AtLAST the improvements are not trivial. Even so, the
nature of a cascading power spectrum simply makes measurements at very small
scales difficult even with an infinitesimal PSF.

###  7.1 Caveats

In the course of making these forecasts many assumptions were made, all of
which will have some impact on the results. One of the most fundamental
choices was to adopt the A12 method, which has two dominant biases. Even so,
its flexibility seems likely to keep it as a preferred method.

The choice of Mach relation (that from Gaspari et al. (2013)) was made out of
simplicity. Adopting a different relation (e.g. those from Zhuravleva et al.,
2012, 2023) are unlikely to dramatically change the results. That is, both
relations depend on decent recovery of spectra over a similar range of scales
to ensure accuracy. Moreover, these two relations are found to be in good
agreement (e.g. Romero et al., 2023; Dupourqué et al., 2024) with respect to
inferred Mach numbers from observed spectra.

To model the turbulent non-thermal pressure profile, one looks to simulations
(e.g Nelson et al., 2014; Shi et al., 2015; Biffi et al., 2016; Angelinelli et
al., 2020), where parameterizations were provided in Battaglia et al. (2012);
Nelson et al. (2014). Battaglia et al. (2012) and Nelson et al. (2014) are in
decent agreement (see Figure 8), while Angelinelli et al. (2020) finds a
shallower slope similar to Eckert et al. (2019). While I took the profile from
Nelson et al. (2014) as my fiducial turbulent non-thermal pressure profile, I
find similar results if I use the profiles in Battaglia et al. (2012) or
Angelinelli et al. (2020).

Another assumption falls into the assumed SB profiles (or equally, pressure
profiles, for SZ observations). Figure 7 showed that the A10 pressure profile
is similar to other pressure profiles; conversely, all sample-averaged
profiles compared are quite similar. The adopted X-ray SB profile was also
found to be consistent with those in Bartalucci et al. (2023); Lyskova et al.
(2023).

####  7.1.1 Additional Noise Concerns

Another critical assumption corresponds to the noise characteristics assumed.
Both SZ and X-ray are susceptible to contamination from (1) compact sources,
(2) intracluster contaminants, (3) extended foreground, and (4) extended
background contaminants. In the X-ray, compact sources are generally masked,
while in the SZ they may be subtracted (at the cost of additional uncertainty)
or masked. In the literature, the treatment of intracluster contaminants is
currently ad hoc, especially as there is not an unambiguous definition of what
constitutes an intracluster contaminant (e.g. Zhuravleva et al., 2015; Romero
et al., 2023). Notwithstanding a lack of clear definition, the approach is
generally to mask these contaminants.

It may be possible that some extended foregrounds or backgrounds can be fit
out (e.g. with some linear combination method). Alternatively, the may
”simply” be an additional source of noise. In the case of AtLAST observations,
I have opted to not take the optimal linear combination sensitivity to
heuristically account for cosmic backgrounds (e.g. CMB and Cosmic Infrared
Background, CIB) and dust. There will also be some floor to the noise
achievable: the confusion limit. For the case of SZ observations, previous
calculations for confusion limits have focused on what is likely to be found
for CMB surveys, with instruments with resolution (FWHM of their beams) worse
than 1 arcminute. Recently, Raghunathan (2022) found that one may expect a
confusion limit $<5\times 10^{-7}$ in Compton $y$ for CMB-HD (FWHM of 0.25
arcminutes at 150 GHz Sehgal et al., 2019). In particular, much of the
confusion would arise from halos with $M_{200}<10^{13}$ M⊙, which is the
expected mass threshold for AtLAST (Di Mascolo et al., 2024).

With these additional noise concerns in mind, it’s clear that the constraints
on hydrostatic mass bias due to turbulent motions will degrade relative to
what was presented in Section 6 and results presented in Figure 14 may also
degrade. Further study quantifying these additional sources of noise is
warranted, but outside the scope of this work.

###  7.2 Implications for AtLAST

With my assumed mapping speed, AtLAST achieves a sensitivity of $\lesssim
5\times 10^{-7}$ in Compton $y$ with 1 Ms on source integration time. Assuming
a 2 deg2 FOV for AtLAST, and a $\sim 20$% telescope downtime over a year, one
could survey to this depth over $\sim 250$ deg2 over 5 years. Di Mascolo et
al. (2024) finds that this sensitivity could be achieved over 4000 deg2 (over
5 years) if the noise adheres to the optimal internal linear combination.
Given cluster catalogs from ACT (Hilton et al., 2018, 2021) and SPT (Bleem et
al., 2015, 2022), the number of clusters with $M_{500}\gtrsim 3.5\times
10^{14}\text{M}_{\odot}$ in $\sim 250$ deg2 and 4000 deg2 will be $\sim 30$
and $\sim 500$, respectively.

AtLAST has the potential to excel at high redshift, especially with
constraints out to $R_{200}$. AtLAST’s beam should not be a hindrance to
achieving significant constraints on $\sigma_{b_{\mathcal{M}}}$ at $R_{500}$
and $R_{200}$ for many clusters, but for smaller clusters (i.e. at higher
redshifts), it presents a potential limitation.

###  7.3 Impact of results

Non-thermal pressure support may include other (non-turbulent) sources of
pressure support which are not negligible. From a cosmological context, the
total mass of clusters is of fundamental interest (e.g. Allen et al., 2011).
Measurements of total mass can be derived from caustics, weak lensing, and CMB
lensing. This may be the preferred avenue for cosmological constraints based
on cluster abundance.

However, access to the total non-thermal pressure from total mass measurements
does not necessarily reveal which mechanisms provide the non-thermal pressure
support. There is thus still a interest in assessing the non thermal pressure
due to turbulent motions. In addition to the use of SB fluctuations, high
resolution X-ray spectroscopy (such as XRISM or NewAthena) can access
turbulent velocities and it can do so directly. However, comparable velocity
constraints from spectroscopic measurements are more expensive than from SB
fluctuations. Therefore, SB fluctuations will remain a critical tool in
assessing turbulent gas velocities at large cluster radii, especially at
moderate to high redshifts.

##  8 Conclusions

In this work I investigated concerns of the accuracy of power spectra
estimation on noiseless synthetic data. Separately I utilized properties of
variance (and power spectra) to derive estimates of uncertainties in power
spectra assuming radial surface brightness profiles and noise profiles. These
two investigations laid the groundwork to forecast what constraints can be
produced with various current and upcoming SZ and X-ray instruments.

On matters of accuracy. In Section 2 I revisited known biases and unveiled a
secondary manner in which beam convolution (which acts directly on $\delta y$
or $\delta S$ and not $\delta y/\bar{y}$ or $\delta S/\bar{S}$) can bias
results. I find that accounting for biases due to the method of power spectrum
calculation, i.e. shape and mask biases, especially for the (Arévalo et al.,
2012, ;A12) method are of primary concern. Conversely, the secondary bias due
to beam convolution and projection biases are of less concern and it may be
advisable to not attempt to correct for a projection bias.

On matters of uncertainty and choice of regions

For plausible, circular, surface brightness profiles (and noise or exposure
profiles), a minimum in spectral uncertainty can be obtained at $\sim
0.2R_{500}$. Physical expectations and instrumental factors lead me to propose
(when possible) a scheme of annuli with boundaries: [0-0.4] $R_{500}$,
[0.4-1.0] $R_{500}$, and if data quality allows, [1.0-1.5] $R_{500}$.

Forecasting constraints for various instruments I assume universal surface
brightness profiles for both X-ray and SZ observations and non-thermal
pressure profiles to derive constraints on amplitude spectra (at various
scales) for a handful of current and proposed X-ray and SZ facilities. These
uncertainties are propagated to an uncertainty in the hydrostatic mass bias
(due to turbulent motions). In this idealized investigation, many instruments
are or will be capable of constraining to hydrostatic mass bias to
$\sigma_{b_{\mathcal{M}}}=0.1$, where Lynx and AtLAST perform best. Accounting
for requiments to sufficiently recover a fluctuation spectrum, I find that the
target of $\sigma_{b_{\mathcal{M}}}=0.1$ on $M_{500}$ can likely be achieved
for clusters with $M_{500}\gtrsim 3\times 10^{14}M_{\odot}$ at $z=0.5$ (and
lower masses for $z<0.5$) by AtLAST and Lynx. Both AtLAST and Lynx face some
challenges to perform constraints down to $M_{500}\simeq 3\times
10^{14}M_{\odot}$ at $z>0.5$. AtLAST appears to have a clear advantage in its
ability to constrain $\sigma_{b_{\mathcal{M}}}$ on $M_{200}$.

I would like to thank my colleagues for their insights and encouragement in
previous works which acted as a springboard for this work. I would like to
further thank Max Gaspari and Paul Nulsen for their theoretical insights,
Rishi Khatri for sanity checks along the way, Gerrit Schellenberger and Ralph
Kraft for their expertise with X-ray observations, and Irina Zhuravleva for
the idea to use a map of $\sqrt{N}$. CR is supported by NASA ADAP grant
80NSSC19K0574 and Chandra grant G08-19117X.

## References

  * Adam et al. (2017) Adam, R., Bartalucci, I., Pratt, G. W., et al. 2017, A&A, 598, A115, doi: 10.1051/0004-6361/201629182
  * Allen et al. (2011) Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409, doi: 10.1146/annurev-astro-081710-102514
  * Alonso et al. (2019) Alonso, D., Sanchez, J., Slosar, A., & LSST Dark Energy Science Collaboration. 2019, MNRAS, 484, 4127, doi: 10.1093/mnras/stz093
  * Andreon et al. (2017) Andreon, S., Trinchieri, G., Moretti, A., & Wang, J. 2017, A&A, 606, A25, doi: 10.1051/0004-6361/201730920
  * Angelinelli et al. (2020) Angelinelli, M., Vazza, F., Giocoli, C., et al. 2020, MNRAS, 495, 864, doi: 10.1093/mnras/staa975
  * Ansarinejad et al. (2024) Ansarinejad, B., Raghunathan, S., Abbott, T. M. C., et al. 2024, arXiv e-prints, arXiv:2404.02153, doi: 10.48550/arXiv.2404.02153
  * Applegate et al. (2014) Applegate, D. E., von der Linden, A., Kelly, P. L., et al. 2014, MNRAS, 439, 48, doi: 10.1093/mnras/stt2129
  * Arévalo et al. (2016) Arévalo, P., Churazov, E., Zhuravleva, I., Forman, W. R., & Jones, C. 2016, ApJ, 818, 14, doi: 10.3847/0004-637X/818/1/14
  * Arévalo et al. (2012) Arévalo, P., Churazov, E., Zhuravleva, I., Hernández-Monteagudo, C., & Revnivtsev, M. 2012, MNRAS, 426, 1793, doi: 10.1111/j.1365-2966.2012.21789.x
  * Arnaud et al. (2010) Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92, doi: 10.1051/0004-6361/200913416
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Bartalucci et al. (2023) Bartalucci, I., Molendi, S., Rasia, E., et al. 2023, A&A, 674, A179, doi: 10.1051/0004-6361/202346189
  * Battaglia et al. (2012) Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012, ApJ, 758, 74, doi: 10.1088/0004-637X/758/2/74
  * Benson et al. (2014) Benson, B. A., Ade, P. A. R., Ahmed, Z., et al. 2014, in Proc. SPIE, Vol. 9153, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII, 91531P, doi: 10.1117/12.2057305
  * Biffi et al. (2016) Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112, doi: 10.3847/0004-637X/827/2/112
  * Blackman & Tukey (1958) Blackman, R. B., & Tukey, J. W. 1958, Bell System Technical Journal, 37, 185, doi: https://doi.org/10.1002/j.1538-7305.1958.tb03874.x
  * Bleem et al. (2015) Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27, doi: 10.1088/0067-0049/216/2/27
  * Bleem et al. (2022) Bleem, L. E., Crawford, T. M., Ansarinejad, B., et al. 2022, ApJS, 258, 36, doi: 10.3847/1538-4365/ac35e9
  * Bulbul et al. (2019) Bulbul, E., Chiu, I. N., Mohr, J. J., et al. 2019, ApJ, 871, 50, doi: 10.3847/1538-4357/aaf230
  * Bulbul et al. (2024) Bulbul, E., Liu, A., Kluge, M., et al. 2024, arXiv e-prints, arXiv:2402.08452, doi: 10.48550/arXiv.2402.08452
  * Chen et al. (2019) Chen, Y.-H., Heinz, S., & Enßlin, T. A. 2019, MNRAS, 489, 1939, doi: 10.1093/mnras/stz2256
  * Chon et al. (2004) Chon, G., Challinor, A., Prunet, S., Hivon, E., & Szapudi, I. 2004, MNRAS, 350, 914, doi: 10.1111/j.1365-2966.2004.07737.x
  * Churazov et al. (2012) Churazov, E., Vikhlinin, A., Zhuravleva, I., et al. 2012, MNRAS, 421, 1123, doi: 10.1111/j.1365-2966.2011.20372.x
  * Di Mascolo et al. (2024) Di Mascolo, L., Perrott, Y., Mroczkowski, T., et al. 2024, arXiv e-prints, arXiv:2403.00909, doi: 10.48550/arXiv.2403.00909
  * Dupourqué et al. (2023) Dupourqué, S., Clerc, N., Pointecouteau, E., et al. 2023, A&A, 673, A91, doi: 10.1051/0004-6361/202245779
  * Dupourqué et al. (2024) —. 2024, arXiv e-prints, arXiv:2403.03064, doi: 10.48550/arXiv.2403.03064
  * Eckert et al. (2017) Eckert, D., Gaspari, M., Vazza, F., et al. 2017, ApJ, 843, L29, doi: 10.3847/2041-8213/aa7c1a
  * Eckert et al. (2019) Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40, doi: 10.1051/0004-6361/201833324
  * Ettori et al. (2020) Ettori, S., Lovisari, L., & Sereno, M. 2020, A&A, 644, A111, doi: 10.1051/0004-6361/202038586
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Gaspari et al. (2013) Gaspari, M., Brighenti, F., & Ruszkowski, M. 2013, Astronomische Nachrichten, 334, 394, doi: 10.1002/asna.201211865
  * Gaspari & Churazov (2013) Gaspari, M., & Churazov, E. 2013, A&A, 559, A78, doi: 10.1051/0004-6361/201322295
  * Gaspari et al. (2020) Gaspari, M., Tombesi, F., & Cappi, M. 2020, Nature Astronomy, 4, 10, doi: 10.1038/s41550-019-0970-1
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Heinrich et al. (2024) Heinrich, A., Zhuravleva, I., Zhang, C., et al. 2024, MNRAS, 528, 7274, doi: 10.1093/mnras/stae208
  * Hilton et al. (2018) Hilton, M., Hasselfield, M., Sifón, C., et al. 2018, ApJS, 235, 20, doi: 10.3847/1538-4365/aaa6cb
  * Hilton et al. (2021) Hilton, M., Sifón, C., Naess, S., et al. 2021, ApJS, 253, 3, doi: 10.3847/1538-4365/abd023
  * Hitomi Collaboration et al. (2016) Hitomi Collaboration, Aharonian, F., Akamatsu, H., et al. 2016, Nature, 535, 117, doi: 10.1038/nature18627
  * Hitomi Collaboration et al. (2018) —. 2018, PASJ, 70, 10, doi: 10.1093/pasj/psx127
  * Hoekstra et al. (2015) Hoekstra, H., Herbonnet, R., Muzzin, A., et al. 2015, MNRAS, 449, 685, doi: 10.1093/mnras/stv275
  * Hofmann et al. (2016) Hofmann, F., Sanders, J. S., Nandra, K., Clerc, N., & Gaspari, M. 2016, A&A, 585, A130, doi: 10.1051/0004-6361/201526925
  * Hoyng (1976) Hoyng, P. 1976, A&A, 47, 449
  * Hurier & Angulo (2018) Hurier, G., & Angulo, R. E. 2018, A&A, 610, L4, doi: 10.1051/0004-6361/201731999
  * Khatri & Gaspari (2016) Khatri, R., & Gaspari, M. 2016, MNRAS, 463, 655, doi: 10.1093/mnras/stw2027
  * Kitayama et al. (2014) Kitayama, T., Bautz, M., Markevitch, M., et al. 2014, arXiv e-prints, arXiv:1412.1176, doi: 10.48550/arXiv.1412.1176
  * Koch et al. (2019) Koch, E. W., Rosolowsky, E. W., Boyden, R. D., et al. 2019, AJ, 158, 1, doi: 10.3847/1538-3881/ab1cc0
  * Korngut et al. (2011) Korngut, P. M., Dicker, S. R., Reese, E. D., et al. 2011, ApJ, 734, 10, doi: 10.1088/0004-637X/734/1/10
  * Lau et al. (2013) Lau, E. T., Nagai, D., & Nelson, K. 2013, ApJ, 777, 151, doi: 10.1088/0004-637X/777/2/151
  * Logan et al. (2022) Logan, C. H. A., Maughan, B. J., Diaferio, A., et al. 2022, A&A, 665, A124, doi: 10.1051/0004-6361/202243347
  * Lyskova et al. (2023) Lyskova, N., Churazov, E., Khabibullin, I. I., et al. 2023, MNRAS, 525, 898, doi: 10.1093/mnras/stad2305
  * Madhavacheril et al. (2020) Madhavacheril, M. S., Sifón, C., Battaglia, N., et al. 2020, ApJ, 903, L13, doi: 10.3847/2041-8213/abbccb
  * Mantz et al. (2016) Mantz, A. B., Allen, S. W., Morris, R. G., & Schmidt, R. W. 2016, MNRAS, 456, 4020, doi: 10.1093/mnras/stv2899
  * Markevitch & Vikhlinin (2007) Markevitch, M., & Vikhlinin, A. 2007, Phys. Rep., 443, 1, doi: 10.1016/j.physrep.2007.01.001
  * Maughan et al. (2016) Maughan, B. J., Giles, P. A., Rines, K. J., et al. 2016, MNRAS, 461, 4182, doi: 10.1093/mnras/stw1610
  * McDonald et al. (2014) McDonald, M., Benson, B. A., Vikhlinin, A., et al. 2014, ApJ, 794, 67, doi: 10.1088/0004-637X/794/1/67
  * Nelson et al. (2014) Nelson, K., Lau, E. T., & Nagai, D. 2014, ApJ, 792, 25, doi: 10.1088/0004-637X/792/1/25
  * Ogorzalek et al. (2017) Ogorzalek, A., Zhuravleva, I., Allen, S. W., et al. 2017, MNRAS, 472, 1659, doi: 10.1093/mnras/stx2030
  * Okabe & Smith (2016) Okabe, N., & Smith, G. P. 2016, MNRAS, 461, 3794, doi: 10.1093/mnras/stw1539
  * Ossenkopf et al. (2008) Ossenkopf, V., Krips, M., & Stutzki, J. 2008, A&A, 485, 917, doi: 10.1051/0004-6361:20079106
  * Ota et al. (2018) Ota, N., Nagai, D., & Lau, E. T. 2018, PASJ, 70, 51, doi: 10.1093/pasj/psy040
  * Peebles (1993) Peebles, P. J. E. 1993, in Princeton Series in Physics, Princeton, NJ: Princeton University Press, —c1993.  http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1993ppc..book.....P&db_key=AST
  * Pinto et al. (2015) Pinto, C., Sanders, J. S., Werner, N., et al. 2015, A&A, 575, A38, doi: 10.1051/0004-6361/201425278
  * Raghunathan (2022) Raghunathan, S. 2022, ApJ, 928, 16, doi: 10.3847/1538-4357/ac510f
  * Rebusco et al. (2007) Rebusco, P., Churazov, E., Böhringer, H., & Forman, W. 2007, in Relativistic Astrophysics Legacy and Cosmology – Einstein’s, ed. B. Aschenbach, V. Burwitz, G. Hasinger, & B. Leibundgut (Berlin, Heidelberg: Springer Berlin Heidelberg), 340–343
  * Romero et al. (2018) Romero, C., McWilliam, M., Macías-Pérez, J.-F., et al. 2018, A&A, 612, A39, doi: 10.1051/0004-6361/201731599
  * Romero et al. (2020) Romero, C. E., Sievers, J., Ghirardini, V., et al. 2020, ApJ, 891, 90, doi: 10.3847/1538-4357/ab6d70
  * Romero et al. (2023) Romero, C. E., Gaspari, M., Schellenberger, G., et al. 2023, ApJ, 951, 41, doi: 10.3847/1538-4357/acd3f0
  * Romero et al. (2024) —. 2024, arXiv e-prints, arXiv:2404.04373, doi: 10.48550/arXiv.2404.04373
  * Roncarelli et al. (2018) Roncarelli, M., Gaspari, M., Ettori, S., et al. 2018, A&A, 618, A39, doi: 10.1051/0004-6361/201833371
  * Sanders et al. (2010) Sanders, J. S., Fabian, A. C., Smith, R. K., & Peterson, J. R. 2010, MNRAS, 402, L11, doi: 10.1111/j.1745-3933.2009.00789.x
  * Sanders et al. (2020) Sanders, J. S., Dennerl, K., Russell, H. R., et al. 2020, A&A, 633, A42, doi: 10.1051/0004-6361/201936468
  * Sayers et al. (2013) Sayers, J., Mroczkowski, T., Zemcov, M., et al. 2013, ApJ, 778, 52, doi: 10.1088/0004-637X/778/1/52
  * Schuecker et al. (2004) Schuecker, P., Finoguenov, A., Miniati, F., Böhringer, H., & Briel, U. G. 2004, A&A, 426, 387, doi: 10.1051/0004-6361:20041039
  * Sehgal et al. (2019) Sehgal, N., Aiola, S., Akrami, Y., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 6, doi: 10.48550/arXiv.1906.10134
  * Sereno (2015) Sereno, M. 2015, MNRAS, 450, 3665, doi: 10.1093/mnras/stu2505
  * Sereno & Ettori (2017) Sereno, M., & Ettori, S. 2017, MNRAS, 468, 3322, doi: 10.1093/mnras/stx576
  * Seta et al. (2023) Seta, A., Federrath, C., Livingston, J. D., & McClure-Griffiths, N. M. 2023, MNRAS, 518, 919, doi: 10.1093/mnras/stac2972
  * Shi et al. (2015) Shi, X., Komatsu, E., Nelson, K., & Nagai, D. 2015, MNRAS, 448, 1020, doi: 10.1093/mnras/stv036
  * Simionescu et al. (2019) Simionescu, A., ZuHone, J., Zhuravleva, I., et al. 2019, Space Sci. Rev., 215, 24, doi: 10.1007/s11214-019-0590-1
  * Sunyaev & Zel’dovich (1972) Sunyaev, R. A., & Zel’dovich, Y. B. 1972, Comments Astrophys. Space Phys., 4, 173
  * The Astropy Collaboration (2018) The Astropy Collaboration. 2018, astropy v3.0.5: a core python package for astronomy, 3.0.5, Zenodo, doi: 10.5281/zenodo.1461536
  * Tukey (1956) Tukey, J. W. 1956, The Annals of Mathematical Statistics, 27, 722.  http://www.jstor.org/stable/2237379
  * Ueda et al. (2018) Ueda, S., Kitayama, T., Oguri, M., et al. 2018, ApJ, 866, 48, doi: 10.3847/1538-4357/aadd9d
  * van Weeren et al. (2019) van Weeren, R. J., de Gasperin, F., Akamatsu, H., et al. 2019, Space Sci. Rev., 215, 16, doi: 10.1007/s11214-019-0584-z
  * Vikhlinin et al. (2001) Vikhlinin, A., Markevitch, M., & Murray, S. S. 2001, ApJ, 549, L47, doi: 10.1086/319126
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang et al. (2023) Wang, C., Oh, S. P., & Ruszkowski, M. 2023, MNRAS, 519, 4408, doi: 10.1093/mnras/stad003
  * Wicker et al. (2023) Wicker, R., Douspis, M., Salvati, L., & Aghanim, N. 2023, A&A, 674, A48, doi: 10.1051/0004-6361/202243922
  * Zhang et al. (2024) Zhang, C., Zhuravleva, I., Markevitch, M., et al. 2024, MNRAS, doi: 10.1093/mnras/stae1022
  * Zhuravleva et al. (2018) Zhuravleva, I., Allen, S. W., Mantz, A., & Werner, N. 2018, The Astrophysical Journal, 865, 53, doi: 10.3847/1538-4357/aadae3
  * Zhuravleva et al. (2023) Zhuravleva, I., Chen, M. C., Churazov, E., et al. 2023, MNRAS, doi: 10.1093/mnras/stad470
  * Zhuravleva et al. (2012) Zhuravleva, I., Churazov, E., Kravtsov, A., & Sunyaev, R. 2012, MNRAS, 422, 2712, doi: 10.1111/j.1365-2966.2012.20844.x
  * Zhuravleva et al. (2019) Zhuravleva, I., Churazov, E., Schekochihin, A. A., et al. 2019, Nature Astronomy, 3, 832, doi: 10.1038/s41550-019-0794-z
  * Zhuravleva et al. (2015) Zhuravleva, I., Churazov, E., Arévalo, P., et al. 2015, MNRAS, 450, 4184, doi: 10.1093/mnras/stv900
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298, doi: 10.21105/joss.01298
  * ZuHone et al. (2016) ZuHone, J. A., Markevitch, M., & Zhuravleva, I. 2016, ApJ, 817, 110, doi: 10.3847/0004-637X/817/2/110
  * ZuHone et al. (2023) ZuHone, J. A., Vikhlinin, A., Tremblay, G. R., et al. 2023

##  Appendix A Formalism with multiple underlying power spectra

Figure 15: The resultant $P_{2D}$ spectrum due to two distinct $P_{3D}$ along
the line of sight is roughly the weighted sum of the two corresponding
$P_{2D}$ (see text). Here, I use $P_{2D,1\text{w}}=w_{1}P_{2D,1}$ and
$P_{2D,2\text{w}}=w_{2}P_{2D,2}$.

In order to explore the effects of multiple $P_{3D}$ on a resultant $P_{2D}$,
I perform a simple numerical simulation wherein I assume two distinct power
spectra, $P_{3D,1}$ and $P_{3D,2}$ lie in separate regions along the line of
sight. In particular I create a 3D box where a region (e.g. $|z|<200$ kpc) is
occupied by ($\delta P/P$) fluctuations corresponding to $P_{3D,1}$ and the
remaining (complementary, e.g. $|z|>=200$ kpc)) region is occupied by ($\delta
P/P$) fluctuations corresponding to $P_{3D,2}$ with the same underlying
Gaussian random field. I numerically integrate $W(\delta P/P)$ for
$W(\theta,z)$ at a specific $\theta$ ($500$ kpc). To calculate $W$, I assume
an A10 profile of a $z=0.3$ and $M_{500}=6\times 10^{14}M_{\odot}$ cluster.

I adopt a weighting scheme for each $P_{\text{3D}}$, $t(r)$, where $t=1$ where
the fluctuations have the underlying power spectrum $P_{3D,1}$ and $t=0$
otherwise to transition between the two power spectra, giving:

$\displaystyle P_{\text{3D,1t}}(k,r)$
$\displaystyle=P_{\text{3D,1}}(k)t(r)\text{, and}$ (A1) $\displaystyle
P_{\text{3D,2t}}(k,r)$ $\displaystyle=P_{\text{3D,2}}(k)[1-t(r)].$ (A2)

When considering the projection, we can apply the weights to the window
functions $W_{1t}(\theta,z)=W(\theta,z)t(\theta,z)$ and
$W_{2t}(\theta,z)=W(\theta,z)[1-t(\theta,z)]$ such that:

$\displaystyle P_{\text{2D},1\text{t}}$ $\displaystyle=\int
P_{\text{3D}}|\tilde{W}_{1\text{t}})^{2}dk_{z}\text{, and}$ (A4)
$\displaystyle P_{\text{2D},2\text{t}}$ $\displaystyle=\int
P_{\text{3D}}|\tilde{W}_{2\text{t}})^{2}dk_{z}.$ (A5)

It is of no surprise that the discontinuity between the two regions induces
aliasing. Thus I alter the weight function from a step function to that of a
Tukey window (tapered cosine). Specifically, I found that a width of six
resolution elements produced a sufficiently smooth numerical transition; these
results can be seen in Figure 15

It is also possible to derive (another) approximation in the case of multiple
underlying $P_{\text{3D}}$. From Equation 1, I derive $P_{\text{2D},1}$ and
$P_{\text{2D},2}$ corresponding to the integration of $P_{\text{3D},1}$ and
$P_{\text{3D},2}$, respectively, over $dk_{z}$ with the same
$|\tilde{W}|^{2}$. I the following values:

$\displaystyle w_{1}$
$\displaystyle=\frac{\int(tW)^{2}dz}{\int(W^{2})dz}\text{, and}$ (A6)
$\displaystyle w_{2}$
$\displaystyle=\frac{\int([1-t]W)^{2}dz}{\int(W^{2})dz},$ (A7)

respectively, . With these weights, I find:

$P_{2D,\text{w}}\approx w_{1}P_{2D,1}+w_{2}P_{2D,2}.$ (A8)

Figure 15 shows that the approximate form does quite well, especially for
features that are likely of interest.

###  A.1 Uncertainty in hydrostatic bias

For each instrument, redshift, cluster mass, injection scale, region, and
exposure time, I calculate the inferred uncertainty in the 2D amplitude
spectrum as indicated in Section 3. The only bias corrections applied to power
spectra are those from the PSF correction, wherein I take $\alpha=3$.101010The
choice of $\alpha=3$ is motivated because this should be indicative of the
bias around the peak, where (necessarily) $d\ln P_{\text{3D}}/d\ln r=-3$. I
use the regions outlined in Section 4: rings with bounds at [0.4, 1.0, 1.5]
$R_{500}$; the innermost ring being a circle of radius $0.4R_{500}$. Values
for $N_{\rm eff}$ are calculated as in Section 2 for each ring and allow me to
deproject spectra to $P_{\rm 3D}$ and $\sigma_{P_{\rm 3D}}$.

To condense expressions, I will adopt the following nomenclatures:

$\displaystyle\alpha_{\text{NT}}$ $\displaystyle=d\ln P_{\text{NT}}/d\ln r$
$\displaystyle\alpha_{\text{th}}$ $\displaystyle=d\ln P_{\text{th}}/d\ln r$
$\displaystyle\alpha_{\mathcal{M}}$
$\displaystyle=d\ln\mathcal{M}_{\text{3D}}/d\ln r$

With the above condensations, I derive:

$\sigma_{\alpha_{\text{NT}}/\alpha_{\text{th}}}=2\left[\left(\frac{\sigma_{\alpha_{\mathcal{M}}}}{\alpha_{\text{th}}}\right)^{2}+\left(\frac{\alpha_{\mathcal{M}}\sigma_{\alpha_{\text{th}}}}{\alpha_{\text{th}}^{2}}\right)^{2}\right]^{1/2}$
(A9)

While the uncertainties on $\mathcal{M}$ at [0.4, 1.0, 1.5] $R_{500}$ is
dependent on $\alpha_{\mathcal{M}}$, I treat them as being independent for
simplicity of calculation. The uncertainty of the hydrostatic mass bias can
then be computed as:

$\sigma_{b_{\mathcal{M}}}=\frac{\gamma}{3}\left[\left(2\mathcal{M}_{\text{3D}}\sigma_{\mathcal{M}}\frac{\alpha_{\text{NT}}}{\alpha_{\text{th}}}\right)^{2}+(\mathcal{M}_{\text{3D}}^{2}\sigma_{\alpha_{\text{NT}}/\alpha_{\text{th}}})^{2}\right]^{1/2}\\\
\times\left(1+\frac{\gamma\mathcal{M}_{\text{3D}}^{2}}{3}\frac{\alpha_{\text{NT}}}{\alpha_{\text{th}}}\right)^{-2}.$
(A10)

In this computation, the fundamental uncertainties are taken as
$\sigma_{\alpha_{\mathcal{M}}}$, $\sigma_{\alpha_{\text{th}}}$, and
$\sigma_{\mathcal{M}}$. The uncertainty in the logarithmic slope of the Mach
number, $\sigma_{\alpha_{\mathcal{M}}}$ is calculated by fitting a single
power law with just two rings (Ring 1 and Ring 2) or all three rings. The
uncertainty in the logarithmic slope of thermal pressure,
$\sigma_{\alpha_{\text{th}}}$, is set to zero for X-ray instruments as a
robust estimation is beyond the intended scope of this paper.

For SZ instruments, $\sigma_{\alpha_{\text{th}}}$ is more readily inferred by
calculating the uncertainty on $d\ln y/d\ln r=\alpha_{y}$ using annuli of
[0.8,1.0] and [1.0,1.2] times a ring edge, where ring edges are $0.4R_{500}$,
$1.0R_{500}$, or $1.5R_{500}$, where I assume
$\sigma_{\alpha_{\text{th}}}\approx\sigma_{\alpha_{y}}$. This approximation is
sufficient at moderate radii and increasingly accurate at large radii, as the
pressure profile tends towards behaving as a single power law. See Appendix B
for a derivation of this relation between the Compton $y$ profile and electron
pressure, $P_{e}$ profile which is valid when $P_{e}$ behaves as a single
power law at $r>\theta$, for a projected radius $\theta$.

Lastly, $\sigma_{\mathcal{M}}$ is calculated by propagating errors through the
power law fit and corresponding estimation of $\mathcal{M}$ at [0.4, 1.0, 1.5]
$R_{500}$.

##  Appendix B Power Law deprojection

Other works have shown analytic expressions for interations along the line of
sight for emissivity or pressure which is characterized as a power law (e.g.
Vikhlinin et al., 2001; Korngut et al., 2011; Romero et al., 2018). Here, let
us consider the case where the Compton $y$ profile is taken to be a power law,
$y=y_{0}\theta^{-\alpha_{y}}$, for projected radii $\theta$. From the inverse
Abel transform, we have:

$P_{e}(r)=\frac{m_{e}c^{2}}{\sigma_{T}}\frac{1}{\pi}\int_{\infty}^{r}\frac{dy}{d\theta}\frac{d\theta}{(\theta^{2}-r^{2})^{1/2}},$
(B1)

where $\theta$ is the projected radius, $m_{e}$ is the electron mass, $c$ is
the speed of light, and $\sigma_{T}$ is the Thomson cross section. If drop
constants, we can write:

$P_{e}(r)\propto\int_{\infty}^{r}\frac{\theta^{-\alpha_{y}-1}d\theta}{(\theta^{2}-r^{2})^{1/2}}.$
(B2)

This equation can be converted into a form recognized as the Beta function by
defining a new variable, $\zeta=r/\theta$. We now derive:

$P_{e}(r)\propto
r^{-\alpha_{y}-1}\int_{0}^{1}\frac{\zeta^{\alpha_{y}}d\zeta}{(1-\zeta^{2})^{1/2}},$
(B3)

where the integral is the Beta function, and is equal to
$\Gamma(\alpha_{y}+1)\Gamma(0.5)/\Gamma(\alpha_{y}+1.5)$. The salient point is
that $P_{e}(r)\propto r^{-\alpha_{y}-1}$, such that the logarithmic slope of
the pressure profile is simply the logarithmic slope of the Compton $y$
profile minus 1.

##  Appendix C Recovery of power spectra with plausible measurement errors

Section 5 outlines that the uncertainties in amplitude spectra, thus in
inferred Mach numbers, are calculated assuming that biases are (perfectly)
accounted for. In Section 2.5.1, it is shown that this can be done quite well
for a spectrum with no measurement uncertainties (only sample variance).

Figure 16: As in Figure 16, but with measurement errors (variance) added to
the existing sample variance. In particular,
$\sigma_{A_{\text{3D}}}(k_{\text{inj}}*4)=0.05$, as a plausible target
sensitivity in Section 7.

In particular, I explore the option of restricting $\alpha=11/3$ (or adopting
priors such as $\alpha=11/3\pm 1/3$), which I consider to be appropriate as
this is not only theoretically expected (in most cases), it is well supported
observationally (e.g. Zhuravleva et al., 2018; Dupourqué et al., 2023;
Heinrich et al., 2024).

Table 2: Recovered Amplitude Spectra

Region |  |  $\ell_{\text{inj}}$ (kpc) |  $A_{\text{3D}(k_{\text{peak}})}$ |  $\sigma_{\text{3D}}$
---|---|---|---|---
Ring1 | Expected | 515 | 0.24 | 0.33
Corrected | 614 | 0.26 | 0.36
Peeled | 559 | 0.24 | 0.33
Ring2 | Expected | 746 | 0.30 | 0.40
Corrected | 739 | 0.35 | 0.47
Peeled | 739 | 0.33 | 0.44
Ring3 | Expected | 1003 | 0.46 | 0.58
Corrected | 890 | 0.45 | 0.58

Note. — Properties of recovered A${}_{\text{3D}}$ with noise.

##  Appendix D Expected vs. measured spectral uncertainties

From Section 3, a relation between a given power spectrum and its uncertainty
are given in Equation 8. Additionally, Equation 11 allows one to convert from
the variance (at a specified pixel or beam resolution) to the equivalent white
noise power spectrum.

Figure 17 acts as a check on these relations. The brown curve adopts the
measured noise power spectrum (Romero et al., 2024) and applies Equation8. The
remaining ”theory” curves (lines) adopt the white noise estimate. Both methods
require scaling by the area-average of the square of the inverse (ICM-only)
surface brightness model.

Figure 17: Verification of Equation 8 with data from SPT-CLJ0232-4421 (Romero
et al., 2024). Both measured and theoretical uncertainties are presented
without bias corrections.

In the case of deriving the Poisson power spectrum via Equation 11 for the
X-ray case, One finds $k_{\rm max}=k_{\rm pix}/\sqrt{\pi}$, where $k_{\rm
pix}=1/l_{\rm pix}$ for $l_{\rm pix}$ being the angular length of a pixel
(assuming square pixels). This estimation underestimates the reported
uncertainties because the reported uncertainties are calculated from
realizations which sample the MCMC chain of $\beta$ profile parameters for the
given cluster.

In the SZ case, $k_{\rm max}=1/(2\pi\sigma_{\rm PSF}^{2})\approx
1/\text{FWHM}^{2}$ for noise reported at beam scales.