# PITSZI: Probing intra-cluster medium turbulence with Sunyaev–Zel’dovich imaging

R. Adam, T. Eynard-Machet, I. Bartalucci, D. Cherouvrier, N. Clerc, L. Di Mascolo, S. Dupourqué, C. Ferrari, J.-F. Macías-Pérez, E. Pointecouteau, et al.

► To cite this version:

R. Adam, T. Eynard-Machet, I. Bartalucci, D. Cherouvrier, N. Clerc, et al., PITSZI: Probing intra-cluster medium turbulence with Sunyaev–Zel’dovich imaging. Astronomy & Astrophysics - A&A, 2025, 694, pp.A182. 10.1051/0004-6361/202452342. hal-04942288

HAL Id: hal-04942288

https://hal.science/hal-04942288v1

Submitted on 12 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# PITSZI: Probing intra-cluster medium turbulence with Sunyaev–Zel’dovich imaging

## Application to the triple merging cluster MACS J0717.5+3745

R. Adam<fnref n="1" />,*, T. Eynard-Machet<fnref n="1" />, I. Bartalucci<fnref n="2" />, D. Cherouvrier<fnref n="3" />, N. Clerc<fnref n="4" />, L. Di Mascolo<fnref n="1" />, S. Dupourqué<fnref n="4" />, C. Ferrari<fnref n="1" />, J.-F. Macías-Pérez<fnref n="3" />, E. Pointecouteau<fnref n="4" />, and G. W. Pratt<fnref n="5" />

Received 23 September 2024 / Accepted 2 January 2025

## ABSTRACT

Turbulent gas motions are expected to dominate the non-thermal energy budget of the intracluster medium (ICM). The measurement of pressure fluctuations from high angular resolution Sunyaev–Zel’dovich imaging opens a new avenue to study ICM turbulence, complementary to X-ray density fluctuation measures. We developed a methodological framework designed to optimally extract information on the ICM pressure fluctuation power spectrum statistics, and publicly released the associated software named PITSZI (Probing ICM Turbulence from Sunyaev–Zel’dovich Imaging). We applied this tool to the New IRAM KIDs Array (NIKA) data of the merging cluster MACS J0717.5+3745 to measure its pressure fluctuation power spectrum at high significance, and to investigate the implications for its non-thermal content. Depending on the choice of the radial pressure model and the details of the applied methodology, we measured an energy injection scale $L_{inj} \sim 800$ kpc. The power spectrum normalization corresponds to a characteristic amplitude reaching $A_{\delta p/p}(k_{peak}) \sim 0.4$. These results were obtained assuming that the ICM of MACS J0717.5+3745 can be described as pressure fluctuations on top of a single (smooth) halo, and were dominated by systematics due to the choice of the radial pressure model. Using simulations, we determined that fitting a radial model to the data can suppress the observed fluctuations by up to ~50%, while a poorly representative radial model can induce spurious fluctuations, which we also quantified. Assuming standard scaling relations between the pressure fluctuations and turbulence, we find that MACS J0717.5+3745 presents a turbulent velocity dispersion $\sigma_v \sim 1200$ km/s, a kinetic to kinetic plus thermal pressure fraction $P_{kin}/P_{kin+th} \sim 20\%$, and we estimate the hydrostatic mass bias to $b_{HSE} \sim 0.3-0.4$. Our results are in excellent agreement with alternative measurements obtained from X-ray surface brightness fluctuations, and in agreement with the fluctuations being adiabatic in nature.

**Key words.** turbulence – techniques: high angular resolution – galaxies: clusters: intracluster medium

## 1. Introduction

In the 1970s, X-ray observations revealed the existence of a diffuse hot gas component that dominates the baryon content of galaxy clusters, known as the intracluster medium (ICM; see Sarazin 1986, for a review). Since then, the ICM has been extensively studied using X-ray observations, but also in the millimetre via the Sunyaev–Zel’dovich (SZ) effect (Sunyaev & Zeldovich 1970, 1972). Today, we know that the physical properties of the ICM are driven by the gravitational collapse of the surrounding cosmic web, in which the kinetic energy of the infalling gas is converted into heat primarily through shocks and turbulent cascades (Kravtsov & Borgani 2012). Consequently, the ICM is essentially thermal. However, a non-negligible fraction of the energy is also channeled into a non-thermal component in the form of turbulence, magnetic fields, and cosmic rays. While the thermal state of galaxy clusters is now well understood, the non-thermal content is not. In particular, turbulent gas motions are expected to

dominate the non-thermal ICM energy budget (e.g. Vazza et al. 2016) and are presumably key to explain cosmic-ray (re-) acceleration and magnetic field amplification (Brunetti & Jones 2014; Donnert et al. 2018). Therefore, the characterization of ICM turbulence is becoming crucial to acquire a comprehensive view of the physical processes involved in the assembly of massive halos and utilizing them for high-precision cluster cosmology.

Robust and precise total mass estimates are indeed essential for the use of the galaxy cluster population as a cosmological probe (Pratt et al. 2019). As the dominant contribution to the non-thermal pressure support, turbulence is expected to play a leading contribution to the hydrostatic mass bias that affects cluster masses derived from the hydrostatic equilibrium assumption (Angelinelli et al. 2020). The hydrostatic mass bias has been widely examined in the literature through both observational data, especially via the calibration of mass scaling relations, and numerical simulations (e.g. Hoekstra et al. 2015; Smith et al. 2016; Eckert et al. 2019; Gianfagna et al. 2021; Muñoz-Echeverría et al. 2024), in particular in light of

<footnote n="1">¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France</footnote>

<footnote n="2">² INAF, IASF-Milano, Via A. Corti 12, 20133 Milano, Italy</footnote>

<footnote n="3">³ Laboratoire de Physique Subatomique et de Cosmologie, Université Grenoble Alpes, CNRS/IN2P3, 53, Avenue des Martyrs, Grenoble, France</footnote>

<footnote n="4">⁴ IRAP, CNRS, Université de Toulouse, CNES, UT3-UPS, Toulouse, France</footnote>

<footnote n="5">⁵ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette, France</footnote>

<footnote>* Corresponding author; remi.adam@oca.eu</footnote>
the tension between cosmic microwave background-based and
cluster-based cosmology (Planck Collaboration XX 2014). Nev-
ertheless, the understanding of its physical origin and its
complete characterization is still an open issue. In fact, the
hydrostatic mass bias is currently one of the major limitations
for the use of ICM-based observations as a cosmological probe
(e.g. Planck Collaboration XXIV 2016; Pratt et al. 2019).

Turbulence is also a key ingredient to explain the diffuse
radio emission that is widely observed in galaxy clusters today,
but not yet fully understood. We can differentiate between radio
relics, which are elongated structures located at the periphery
of clusters, and radio halos, which coincide spatially with the
thermal gas (van Weeren et al. 2019). While relics are believed
to trace merger shock acceleration (van Weeren et al. 2010),
turbulence is often invoked to explain the failure of the sim-
ple diffuse shock acceleration model (Fujita et al. 2015). Radio
halos are generally further divided into mini-halos (in relaxed
cool-core clusters), megaparsec-sized giant radio halos (in dis-
turbed objects), and mega-halos, which extend out to the extreme
periphery. In all cases, observations suggest that turbulence is a
key player to re-accelerate synchrotron-emitting relativistic elec-
trons (Brunetti & Jones 2014). In this scenario, mini-halos may
be powered by active galactic nuclei (AGN) feedback or slosh-
ing motions (Mazzotta & Giacintucci 2008), giant-halos by the
energy dissipated in mergers (Cassano et al. 2010), and mega-
halos by large-scale structure formation (Beduzzi et al. 2023).
However, the details of the mechanisms involved are still poorly
understood and the origin of the particles re-accelerated by tur-
bulence is debated (e.g. Pinzke et al. 2017; Adam et al. 2021).
Characterizing the connection between these non-thermal pro-
cesses and turbulence, in different environments will therefore
be key to understand the plasma physics at play.

Ideally, turbulence should be observed by directly prob-
ing the ICM velocity field (Simionescu et al. 2019). This will
become routine with high-resolution X-ray spectroscopy experi-
ments (XRISM Science Team 2020; Barret et al. 2020), follow-
ing on from the early results from the *Hitomi* satellite towards
the core of the Perseus cluster (Hitomi Collaboration 2016).
Alternatively, turbulence can be probed indirectly through fluc-
tuations of ICM thermodynamic quantities, assuming these are
being produced by the turbulent velocity field. This approach
relies on the calibration of the power spectrum of the target
ICM fluctuation with the turbulent Mach number, as studied
in numerical simulations (Zhuravleva et al. 2014; Gaspari et
2014; Mohapatra et al. 2020, 2021, 2022; Simonte et al. 2022;
Zhuravleva et al. 2023). To date, such studies have been con-
ducted through the statistics of X-ray surface brightness and
pressure fluctuations obtained relative to a smooth model (e.g.
Schuecker et al. 2004; Churazov et al. 2012; Dupourqué et al.
2023, 2024; Heinrich et al. 2024). However, the physical inter-
pretation of these fluctuations as being generated by turbulence
is not straightforward from X-rays only. They may also arise
from structures that are independent of turbulence (e.g. clump-
ing, sloshing features, and shocks; see Churazov et al. 2012;
Dupourqué et al. 2023, for discussions). Moreover, the thermo-
dynamic nature of the fluctuations may depend on the ICM
environment<fnref n="1" />. For instance, isobaric fluctuations may arise from
slow sloshing motions, gas bubbles induced by AGN feedback
will appear as isothermal fluctuations, and adiabatic fluctuations
are associated with more vigorous motions such as with weak
shocks or sonic turbulence (see, e.g. Zhuravleva et al. 2018).

These regimes are very difficult to access directly through X-ray
surface brightness fluctuations alone, since they probe essen-
tially density perturbations<fnref n="2" />.

The SZ effect offers an independent and complementary
probe of the ICM (see Birkinshaw 1999; Mroczkowski et al.
2019, for reviews). Unlike the X-ray emission, which traces gas
density and temperature (with imaging and spectroscopy, respec-
tively), the SZ effect directly measures the thermal pressure<fnref n="3" />.
Furthermore, the SZ surface brightness is independent of the
redshift, which makes it very attractive for the characterization
of distant objects, provided that sufficiently high angular resolu-
tion and sensitivity are available. As such, the SZ surface bright-
ness fluctuations are obvious tracers to investigate the impact
of turbulence on the ICM physics. To date, only three such stud-
ies have been conducted on individual objects (Khatri & Gaspari
2016; Romero et al. 2023, 2024), albeit with a limited signal-
to-noise-ratio. Romero (2024) also undertook a forecast study
on constraints from SZ surface brightness fluctuations in galaxy
clusters.

Here, we aim to develop a methodological framework that
efficiently extracts information on the ICM pressure fluctuation
power spectrum from SZ observations. We also release the cor-
responding software, named PITSZI (Probing ICM Turbulence
from Sunyaev–Zel’dovich Imaging), to the public. We apply
this tool to NIKA observations of the triple merging cluster
MACS J0717.5+3745, which probe the SZ signal at an angular
resolution of 18 arcsec, to measure its pressure fluctuation power
spectrum and explore the implications for its non-thermal energy
content.

The present article is structured as follows. Section 2 gives
an overview of the PITSZI software. In Section 3, we present the
NIKA data of the cluster MACS J0717.5+3745, which we use
as a test case throughout the paper. The modelling of the pres-
sure distribution and its SZ observable, accounting for instru-
mental effects and noise, is discussed in Section 4. Section 5
presents the methodology used to infer constraints on the pres-
sure profile and the pressure fluctuation power spectrum. We dis-
cuss the results highlighting the implications for the non-thermal
ICM physics of MACS J0717.5+3745 in Section 6. Section 7
gives a summary of the paper and details our conclusions. Sev-
eral appendices complement the paper by quantifying the relevant
sources of uncertainty and systematic effects, and discuss how
background and foreground contamination may be accounted
for. Throughout the paper, we assume a flat ΛCDM cosmology
with H₀ = 70 km s⁻<fnref n="1" /> Mpc⁻<fnref n="1" /> and Ωm = 0.3. At the redshift
of MACS J0717.5+3745 (z = 0.546), 1 arcmin corresponds to
395 kpc.

## 2. General overview and structure of the PITSZI software

PITSZI is a python-based code designed to model the ICM
pressure distribution in galaxy clusters. It allows us to generate
SZ image mocks and derive constraints on the power spectrum of
the pressure fluctuations and its implications to the non-thermal

<footnote n="2">² The soft X-ray band primarily reflects density fluctuations, while the hard X-ray band captures fluctuations influenced by both density and temperature. Consequently, X-ray surface brightness fluctuations can, in principle, provide insights into the thermodynamic nature of these variations. However, achieving the necessary photon statistics in the hard band is challenging, which restricts the range of observable scales.</footnote>

<footnote n="3">³ It may also measure the ICM velocity from the kinetic SZ effect (kSZ) in individual systems, but this remains extremely challenging (Sayers et al. 2013; Adam et al. 2017b).</footnote>

<footnote n="1">¹ Pressure and density perturbations are related via $\delta P/P = \Gamma \delta \rho / \rho$, with $\Gamma = [0, 1, 5/3]$ for isobaric, isothermal, and adiabatic fluctuations.</footnote>
ICM physics. Before entering into the details of the modelling and analysis framework, this section gives a brief overview of the code.

## 2.1. Overview and key ingredients

Modelling of the pressure distribution requires a description of the ICM pressure profile, to which pressure fluctuations are added via the modelling of their power spectrum. Given a 3D grid sampling along the sky coordinates and the line of sight, the pressure is integrated to produce a mock SZ image. The SZ image may then be convolved with an instrument response function (beam smoothing and data reduction transfer function) and realistic noise may be added to mimic the data obtained from a given telescope. PITSZI was developed in the context of the analysis of NIKA and NIKA2 data (Adam et al. 2014; Catalano et al. 2014; Adam et al. 2018a), but it can in principle be applied to data obtained from any SZ imaging instrument (e.g. Planck, ACT, SPT, Bolocam, MUSTANG, MUSTANG2).

To extract constraints on the pressure fluctuation power spectrum in three dimensions, one must first compute the power spectrum of the SZ image in two dimensions after accounting for the radial pressure model and given a region of interest (e.g. within $R_{500}$). This power spectrum is then compared to a model, which can be obtained via different methods, accounting for projection effects, the instrument response function, the masking of bad pixels, the noise, and possible astrophysical contaminants. To avoid noise bias and mitigate systematic effects in the data, one may use the cross power spectra of two independent data (sub)sets, if available.

With the constraints on the pressure fluctuation power spectrum in hand, PITSZI implements several approaches from the literature to derive the non-thermal properties of the ICM. This includes the turbulent velocity or the non-thermal to thermal energy ratio.

## 2.2. Code structure

The structure of the code is illustrated in Figure 1. PITSZI consists of four main classes.

1. The class *Model* is designed to describe the cluster object and the pressure sampling. In addition to basic properties such as sky coordinates, mass and redshift, it includes the description of the pressure profile and the pressure fluctuations as the key physical properties.

2. The class *Data* can be used to construct a data object, defined according to an SZ image, its astrometry information and the instrument response function and noise properties. This class also implements functions to handle the noise properties and generate noise Monte Carlo realizations.

3. The classes *InferenceRadial* and *InferenceFluctuation* are dedicated to extraction of constraints on the radial and fluctuation model parameters. Several methods to compute the SZ fluctuation power spectrum from a given pressure model are implemented, in addition to two different fitting methods based on non-linear least squares or Markov Chain Monte Carlo.

4. The class *Physics* derives physical constraints on the non-thermal ICM physics given the pressure fluctuation power spectrum. It implements various results from the literature that connect turbulence to pressure fluctuations.

In addition, libraries of utility functions, including those related to power spectrum analysis, are used by the different classes.

## 2.3. Sampling

The ICM pressure is sampled in real space by defining a box centered at the cluster redshift and on a reference map coordinate. This setup assumes that the cluster extent is sufficiently small on the sky to approximate the box as a simple pixelated rectangular parallelepiped. The line of sight is gridded given the pixel resolution and the box size along the line of sight (typically $3 \times R_{500}$). The projected grid may be provided via a header (e.g. corresponding to real data) or by defining the pixel resolution and the size of the field of view along longitude and latitude directions.

Fast Fourier transform definitions of the spatial frequency space, $k$, follow that of the ` numpy `. ` fft` package<fnref n="4" />, which is used in PITSZI. We note that the $k$-space does not include modes for which $|k| > \min(\max(k_1), \max(k_2), \max(k_3))$, where $k_1, k_2, k_3$ are the spatial frequency along the plane of the sky and line of sight, respectively, so that the fluctuations would not be isotropic at these scales (see Ponthieu et al. 2011, and in particular their Figure 1), but they are usually well below the limit given by the telescope angular resolution so that this is not relevant in practice.

# 3. Data: MACS J0717.5+3745 as observed with NIKA

In this paper, we apply the PITSZI code to the cluster MACS J0717.5+3745, as imaged by the NIKA camera. This section discusses the target choice, presents the NIKA data, and reviews the possible contaminants that are accounted for in the analysis.

## 3.1. MACS J0717.5+3745

MACS J0717.5+3745 is one of the most complex and dynamically active galaxy clusters known. It is a triple merging system at redshift $z = 0.546$, involving at least four subclusters (Ma et al. 2009), and one of the most massive clusters in the Universe ($M(<1\,\text{Mpc}) \sim 2 \times 10^{15}\,\text{M}_\odot$, Jauzac et al. 2018). MACS J0717.5+3745 is characterized by the presence of very large amounts of both thermal and non-thermal energy. Its ICM was found to be extremely hot, reaching up to $25\,	ext{keV}$ (Adam et al. 2017a), and it hosts one of the most powerful known radio halos (Rajpurohit et al. 2021a) and a complex radio relic (Rajpurohit et al. 2021b). Recently, the analysis of ICM density fluctuations with *Chandra* data in a sample of 80 clusters found that MACS J0717.5+3745 displayed the highest characteristic velocities in the sample (Heinrich et al. 2024). In addition, MACS J0717.5+3745 is currently the only individual cluster in which a kinetic SZ signal has been detected (Mroczkowski et al. 2012; Sayers et al. 2013; Adam et al. 2017b). As a well-studied system, the SZ images of MACS J0717.5+3745 are among the best available, both in terms of signal-to-noise ratio and angular resolution. Given these properties, this cluster is an ideal test and benchmark target to search for pressure fluctuations via SZ mapping. However, as will be further discussed in Section 6, we note that the interpretation of the results in terms of turbulence may not be straightforward given the complexity of the system. Nevertheless, it remains an excellent target for addressing the methodology that we develop here, which is the main goal of the paper.

<footnote n="4">⁴ https://numpy.org/doc/stable/reference/routines.fft.html</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Overview of the code structure, user interfaces, and links between the different components. The yellow boxes display the main classes and the inner orange boxes indicate the dependence on the main sub-classes, with class names indicated at the top together with the name of the main file where this is defined. The grey box gives the name of the utility function files. The blue arrows indicate the interfaces with the user. The orange arrows indicate how the different modules are inputs to one another. The grey arrows show how utility functions are used by the classes.</figcaption></figure>

## 3.2. Data

We use the 150 GHz NIKA data obtained from the IRAM 30m telescope under projects 237-13 and 222-14 (displayed in Figure 2). These correspond to about 13 h of unflagged data taken under good atmospheric conditions. Calibration uncertainties were estimated to be 7% at 150 GHz and the Gaussian beam FWHM (full width at half maximum) was measured to be 18.2 arcsec<fnref n="5" />. The data processing induced transfer function, that filters scales larger than the field-of-view (~2 arcmin) was estimated following Adam et al. (2015). Noise Monte Carlo realizations were generated following the method described in Adam et al. (2016). These data have been made publicly available in Adam et al. (2018b); see Section 7. We refer to Adam et al. (2017b,a) for further details regarding the data reduction.

## 3.3. Contaminants

In addition to the thermal SZ signal and the target pressure fluctuations, a few astrophysical contaminants may affect the

data. We account for these contaminants as follows. The atmospheric and instrumental noise will be discussed separately in Section 4.

First of all, MACS J0717.5+3745 is known to be contaminated by kinetic SZ signal. We use the best fit template obtained in Adam et al. (2017b) as a correction. However, following their recommendation, since the model is degenerate and presents large uncertainties, we only used it to test the impact of the kinetic SZ signal as a systematic effect. We refer to Appendix D for details regarding the quantification of this effect.

The NIKA 150 GHz images are affected by resolved individual radio and infrared sources. We used the point source model built in Adam et al. (2017b) to correct for these. Additionally, we built a conservative point source mask by thresholding the point source template at 0.1 mJy/beam. We will apply this mask to the data to test the systematic effects due to possible remaining point source residuals.

In addition to the sources directly detected in the data, we expect the SZ signal to be affected by a contribution from unresolved galaxies below the detection threshold, which constitute the cosmic infrared background (CIB). We followed the methodology developed by Béthermin et al. (2017) to generate 1000 mock CIB realizations which were used to account for this contaminant in the following. We refer to Appendix A.1 for further details regarding the construction of these mocks.

<footnote n="5">⁵ This corresponds to 0.3 arcmin FWHM. For comparison, the Planck, ACT, SPT, Bolocam, and MUSTANG2 beam FWHM are about 10 arcmin, ~2 arcmin, ~1.5 arcmin, 0.97 arcmin and 0.16 arcmin, respectively, in the context of SZ observations.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. NIKA Compton parameter (y) data. From left to right, we show the NIKA MACS J0717.5+3705 input data, its best-fit radial model map, a jackknife map (half difference between two equal data subsets) and the residual between the data and the radial model, δy. The data were smoothed with a 15 arcsec gaussian kernel (FWHM) for visual purposes. The contours are showing the S/N in units of 3σ. The comparison with mock data can be observed in Figure 3.</figcaption></figure>

The contribution from undetected lower mass halos in the field of view may also potentially affect the data. However, we verify in Appendix A.2, using simulations, that this signal is negligible in the present case.

Finally, the diffuse galactic emission and the primary CMB anisotropies may also affect the data and should be accounted for if relevant. In the case of the NIKA data, however, both of these contaminants were shown to be negligible given the image sensitivity and the angular scales involved (Adam et al. 2016, 2017b).

# 4. Modelling the pressure fluctuations and its Sunyaev-Zel'dovich observable

This section discusses how the pressure distribution is modelled in PITSZI, how it is used to generate mock SZ images, and how the observables are related to the underlying pressure fluctuation power spectrum.

## 4.1. Pressure distribution and projection

Let us first consider the modelling of the pressure and its projection to obtain the requested SZ observable and the power spectrum of its fluctuations.

### 4.1.1. The Sunyaev-Zel'dovich effect

The thermal SZ effect is due to the inverse Compton scattering of CMB photons by the hot ICM (Sunyaev & Zeldovich 1970, 1972). The SZ surface brightness $\Delta I_v$, at frequency $v$, relative to the CMB specific intensity $I_0$, can be expressed as (see Mroczkowski et al. 2019, for a review)

$$ \frac{\Delta I_v}{I_0} = f(v, T) \times y, \qquad (1) $$

where $f(v, T)$ is the characteristic SZ spectrum and $y$ the Compton parameter, which gives the amplitude of the SZ effect. The SZ spectrum slightly depends on the temperature $T$ for very hot gas in the relativistic regime (Pointecouteau et al. 1998; Itoh et al. 1998). We used the coefficients listed in Adam et al. (2017b), together with their XMM-Newton X-ray spectroscopic temperature map, to convert from surface brightness to Compton parameter as a function of temperature.

The Compton parameter is related to the electron pressure $P_e$, via integration along the line of sight $\ell$, such that

$$ y = \frac{\sigma_T}{m_e c^2} \int P_e d\ell. \qquad (2) $$

The parameter $\sigma_T$ is the Thomson cross section, $m_e$ is the electron rest mass, and $c$ is the speed of light.

In PITSZI, we decompose the pressure into a mean radial component and pressure fluctuations such that

$$ \begin{aligned} P_e(x_1, x_2, x_3) &= \bar{P}_e(r) + \delta P_e(x_1, x_2, x_3) \\ &= \bar{P}_e(r) \left( 1 + \frac{\delta P_e(x_1, x_2, x_3)}{\bar{P}_e(r)} \right). \end{aligned} \qquad (3) $$

Here, $x_1, x_2$, and $x_3 \equiv \ell$ are the coordinates in the plane of the sky and along the line of sight, and $r = \sqrt{x_1^2 + x_2^2 + x_3^2}$ is the 3D radius.

### 4.1.2. The mean pressure profile

The choice of the mean radial distribution is key for fluctuation analyses (see discussions in, e.g. Romero et al. 2023; Dupourqué et al. 2023). Here, we discuss how the pressure profile can be modelled in PITSZI and we refer to Appendix C for a study of the impact of the choice of the model on the results.

**Mean pressure profile.** PITSZI was built taking advantage of the developments made for the MINOT software$^6$. As such, it implements the same profile libraries, including: the $\beta$, double-$\beta$, Navarro-Frenk-White (NFW), generalized Navarro-Frenk-White (gNFW), simplified Vikhlinin model (SVM), or any user-defined profile to be interpolated in log-space (see Adam et al. 2020, for details). It is also possible to set the pressure profile given the mass and the redshift of the cluster, according to various universal models available from the literature (e.g. Arnaud et al. 2010; Planck Collaboration Int. V 2013; Ghirardini et al. 2019; Melin & Pratt 2023).

**Triaxiality.** The expected shape of dark matter halos is triaxial rather than spherical. Even though the gas should be rounder, its structure is still driven by the dark matter collapse (Limousin et al. 2013). It is therefore mandatory to allow for triaxiality in the radial pressure modelling in order not to confuse pressure fluctuations with deviations from sphericity.

In PITSZI, the shape of the ellipsoid is parametrized via the minor-to-major and intermediate-to-major axis ratio, $q_{min}$ and $q_{int}$, respectively. We define the ellipsoid radius as

$$ \xi = \frac{\sqrt{\left(\frac{x_1}{q_{\min}}\right)^2 + \left(\frac{x_2}{q_{\int}}\right)^2 + (x_3)^2}}{\sqrt{\left(\frac{x_1}{q_{\min}}\right)^2 + \left(\frac{x_2}{q_{\int}}\right)^2 + (x_3)^2}} \qquad (4) $$
The pressure profile is then given by $\bar{P}_e(r) \rightarrow \bar{P}_e(\xi)$.

The axes are initially defined so that the major axis is aligned with the line of sight, the intermediate axis aligned with the latitude coordinates, and the minor axis aligned with the longitude. The cluster is then rotated according to the three Euler angles, $\phi_1, \phi_2, \phi_3$, to allow for any orientation, using the scipy.spatial transform.Rotation package. The rotation follows the so-called ZXZ sequence, meaning that $\phi_1$ corresponds to a first rotation around the line of sight axis ($x_3$), $\phi_2$ corresponds to a second rotation around the longitude axis ($x_1$), and $\phi_3$ corresponds to a third rotation around the line of sight axis ($x_3$).

### 4.1.3. Description of the pressure fluctuations

The pressure fluctuations relative to the radial pressure profile $(\frac{\delta P_e(x_1, x_2, x_3)}{P_e(r)})$ are modelled as an isotropic random field. This field is described by its power spectrum and a probability density function. While radial evolution of the fluctuations are usually expected (Shaw et al. 2010; Battaglia et al. 2012; Nelson et al. 2014), this is not currently implemented in PITZSI, which assumes a single pressure fluctuation power spectrum over the whole ICM.

For the time being, only one power spectrum model is available, a cutoff powerlaw function, but any other function can easily be implemented. This 3D model power spectrum, which we use in the following, is parameterized as

$$ \mathcal{P}_{\delta P/\bar{P}}(k_{3D}) = \sigma_{\bar{P}}^2 \frac{k_{3D}^{\alpha} \exp\left(-\frac{1}{k_{3D}^2 L_{\text{inj}}^2}\right) \exp\left(-k_{3D}^2 L_{\text{dis}}^2\right)}{\int 4\pi k_{3D}^{\alpha+2} \exp\left(-\frac{1}{k_{3D}^2 L_{\text{inj}}^2}\right) \exp\left(-k_{3D}^2 L_{\text{dis}}^2\right) dk_{3D}}. \quad (5) $$

Here, the notation $k_{3D}$ specifies that we consider the power spectrum in three dimensions (contrary to the 2D power spectrum of SZ images discussed later) and $k_{3D}$ is the norm of the wavenumber $k = (k_1, k_2, k_3)$. The parameter $\sigma_{\bar{P}}$ represents the standard deviation of the pressure fluctuations and gives the normalization of the power spectrum$^7$. The parameter $\alpha$ is the slope of the spectrum in the inertial range. Its canonical value is $\alpha = -11/3$ for a Kolmogorov cascade (Kolmogorov 1941). The parameters $L_{\text{inj}}$ and $L_{\text{dis}}$ are the injection scale and the dissipation scale, respectively. While the model given in Equation (5) has been largely used in the literature, the true power spectrum could be more complicated, and this choice represents one of the strong assumption of the present work.

With the power spectrum in hand, it is also useful to consider the unitless characteristic spectrum, defined as

$$ A_{\delta P/\bar{P}}(k_{3D}) = \sqrt{4\pi k_{3D}^3 \mathcal{P}_{\delta P/\bar{P}}(k_{3D})}. \quad (6) $$

The pressure fluctuations are also described given a probability distribution function, for which either a gaussian or a lognormal distribution is available. We use the formalism described in Greiner & Enßlin (2015) to generate lognormal fluctuations that follow the required power spectrum.

### 4.1.4. Relation between pressure fluctuation (3D) and SZ surface brightness fluctuations (2D) power spectra

The pressure fluctuation power spectrum is related to the SZ surface brightness fluctuation power spectrum by (e.g. Churazov et al. 2012; Khatri & Gaspari 2016, for the formalism in X-ray and SZ, respectively)

$$ \mathcal{P}_{\delta y/\bar{y}}(k_{2D}) = \int \mathcal{P}_{\delta P/\bar{P}}(k_{3D}) |\tilde{W}(k_1, k_2, k_3)|^2 dk_3. \quad (7) $$

The quantity $\mathcal{P}_{\delta y/\bar{y}}(k_{2D})$ is the power spectrum of the residual SZ image between the radial model $\bar{y}(R)$ and the total SZ image $y(x_1, x_2)$, such that

$$ \frac{\delta y}{\bar{y}} = \frac{y(x_1, x_2) - \bar{y}(R)}{\bar{y}(R)}, \quad (8) $$

with $R = \sqrt{x_1^2 + x_2^2}$ the projected radius and $k_{2D}$ the norm of the wavenumber in the plane of the sky. The quantity $\tilde{W}(k_1, k_2, k_3)$ is the Fourier transform of the window function given by

$$ W(x_1, x_2, x_3) = \frac{\sigma_{\mathrm{T}}}{m_e c^2} \frac{\bar{P}_e(r)}{\bar{y}(R)}. \quad (9) $$

In practice, the window function quickly drops for $k_3$ above a given cutoff, such that Equation (7) is well approximated by

$$ \mathcal{P}_{\delta y/\bar{y}}(k_{2D}) \approx \mathcal{P}_{\delta P/\bar{P}}(k_{3D}) \int |\tilde{W}(k_1, k_2, k_3)|^2 dk_3. \quad (10) $$

Using this simplification implies that the map

$$ C_{3D \rightarrow 2D} = \int |\tilde{W}(k_1, k_2, k_3)|^2 dk_3 \quad (11) $$

provides a simple normalization to convert the 3D to 2D power spectrum. In practice, one needs to average the conversion factor over the effective area over which the power spectrum is computed to obtain

$$ C_{3D \rightarrow 2D}^{(\text{eff})} = \frac{\sum_{\text{map pixels}} \omega(x_1, x_2)^2 C_{3D \rightarrow 2D}(x_1, x_2)}{\sum_{\text{map pixels}} \omega(x_1, x_2)^2}, \quad (12) $$

where $\omega(x_1, x_2)$ is a weight map (e.g. data mask, region of interest, or any weighting function that is applied to the data prior computing the power spectrum; see Section 4.2.2). In the end, Equation (7) becomes

$$ \mathcal{P}_{\delta y/\bar{y}}(k_{2D}) = C_{3D \rightarrow 2D}^{(\text{eff})} \times \mathcal{P}_{\delta P/\bar{P}}(k_{3D}). \quad (13) $$

This approximation was introduced by Churazov et al. (2012), who studied the X-ray surface brightness fluctuations in the core of the Coma cluster with XMM-Newton and Chandra. However, it may significantly break at large scales, especially for 3D power spectra with low power on large scales. Additionally, the approximation holds better for flatter pressure profiles since the window function in this case would be sharper (see also Clerc et al. 2019 for more insights). The limits of this approximation are quantified in Appendix B, where we see that it is in fact fairly good in the vast majority of practical cases. In the present analysis, it implies a bias on the power spectrum model that remains within about 10% over the range of scales that we probe.

<footnote>$^7$ In practice, the output standard deviation of a pressure fluctuation box obtained by generating a random field that follows this spectrum may be lower than $\sigma_{\bar{P}}$ if the scales below and above those sampled by the model are not available.</footnote>
## 4.2. Power spectrum measurement and modelling accounting for the instrumental response, data weighting and noise

After calibration, the observed Compton parameter map obtained from SZ imaging instruments can usually be expressed as

$$y_{\text{observed}} = H \otimes B \otimes y_{\text{true}} + n, \quad (14)$$

where $y_{\text{true}}$ is the true incoming signal. The quantity B represents the beam smoothing induced by the limited telescope angular resolution and H is the transfer function induced by the data reduction that usually filters the scales larger than the field of view (Adam et al. 2015; Romero et al. 2017; Ruppin et al. 2018). The term n is the noise. The map used to extract the 2D SZ power spectrum is defined as

$$S = \omega \times \frac{y_{\text{observed}} - y_{\text{radial model}}^{(\text{num})}}{y_{\text{radial model}}^{(\text{den})}}, \quad (15)$$

with $\omega$ the weight map. The radial component model in the numerator and denominator, $y_{\text{model}}^{(\text{num, den})}$, may slightly differ (see below).

Let us now consider the measurement of the power spectrum itself and its modelling in the presence of these instrumental effects. This includes beam smoothing and transfer function filtering due to the mapmaking procedure, data masking or weighting, and noise.

### 4.2.1. Instrumental response function

Leaving the weight and the noise aside for the moment, in Fourier space, application of the instrument response function translates into

$$\mathcal{P}_{\text{observed}}(k_{2\text{D}}) \approx B(k_{2\text{D}})^2 H(k_{2\text{D}})^2 \mathcal{P}_{\delta y/\bar{y}}(k_{2\text{D}}). \quad (16)$$

Equation (16) is not a strict equality due the fact that the convolution of the instrument response with the observable S is not the same as the ratio of the convolved numerator and denominator. While $y^{(\text{num})}$ accounts for the instrumental response to match the observed data, we therefore define $y^{(\text{den})}$ such that it does not account for the transfer function, but is smoothed by the beam. This implies that deviations from equality only appears on scales below the beam cutoff, where the information on fluctuations is lost anyway.

In the case where the cross spectra between two maps a and b are used, the beam and transfer function become $B^2(k_{2D}) \rightarrow B^{(a)}(k_{2D}) \times B^{(b)}(k_{2D})$ and $H^2(k_{2D}) \rightarrow H^{(a)}(k_{2D}) \times H^{(b)}(k_{2D})$, respectively.

For a Gaussian beam with FWHM = $2\sqrt{2\ln(2)}\sigma_{\text{beam}}$, the transmission reads

$$B(k_{2\text{D}}) = \exp(-2\pi^2 k_{2\text{D}}^2 \sigma_{\text{beam}}^2). \quad (17)$$

The exact shape of the transfer function will depend on the details of the data reduction and shall be provided as part of the data products for any quantitative analysis.

### 4.2.2. Data weighting

In practice, the $\delta y/\bar{y}$ image is not used to extract the power spectrum directly. Indeed, one needs to account for possible masking, or more generally data weighting. This may arise due to

the presence of contaminant point sources that cannot be subtracted perfectly, because we might want to extract the power spectrum over a given region of interest (e.g. within $R_{500}$), or in the case where data weighting is required (e.g. to account for noise inhomogeneity). In PITSZI, we implemented the methodology developed by Ponthieu et al. (2011) for the POKER software in order to account for data masking and weighting in the power spectrum modelling.

Assuming the SZ fluctuation image to be weighted as $S \equiv \omega \times \delta y/\bar{y}$, its pseudo-Fourier coefficients, for each mode $m_1, m_2$, are given by

$$\tilde{S}_{m_1, m_2} = \sum_{\mu_1, \mu_2} \omega_{\mu_1, \mu_2} \left(\frac{\delta y}{\bar{y}}\right)_{\mu_1, \mu_2} e^{-2i\pi(\mu_1 m_1 / N_1 + \mu_2 m_2 / N_2)}, \quad (18)$$

with $N_1, N_2$ the number of pixels along latitude and longitude, indexed $\mu_1$ and $\mu_2$, respectively. We can show that (Ponthieu et al. 2011, their Appendix A, albeit with different Fourier transform normalization)

$$\tilde{S}_{m_1, m_2} = \sum_{m'_1, m'_2} \tilde{S}_{m_1, m_2}^{\text{unweighted}} \frac{K_{m_1, m'_1}^{m_2, m'_2}}{N_1 N_2}, \quad (19)$$

with $\tilde{S}_{m_1, m_2}^{\text{unweighted}}$ the pseudo-Fourier coefficients of the unweighted image $\delta y/\bar{y}$ and $K_{m_1, m'_1}^{m_2, m'_2}$ a matrix that depends on the pseudo-Fourier coefficients of the weight map. The power spectrum of the weighted image is given by $\mathcal{P}_{m_1, m_2} = |\tilde{S}_{m_1, m_2}|^2$ for each mode $m_1, m_2$ and $\mathcal{P}(k_{2D})$ is obtained by binning $\mathcal{P}_{m_1, m_2}$ in $k_{2D}$ assuming isotropy. We refer to Ponthieu et al. (2011) for more insight into the mathematical details.

The comparison between the power spectrum of weighted data and a model of SZ fluctuations should therefore either convolve the model to account for the effect of the weights (using Equation (19) together with binning), or deconvolve the data from the weights. Moreover, it is advisable to embed the observed sky patch into a larger patch and pad it with zeros to prevent aliasing from scales larger than the observed sky. The observation patch should be apodized to minimize large-scale aliasing. Other methods have been used in the literature to estimate the power spectrum in this context. In particular, most results have used the method from Arévalo et al. (2012), which differs from ours. Appendix B quantifies the underlying uncertainties.

### 4.2.3. Noise

The contribution from the noise is obtained by applying the same process to noise Monte Carlo realizations. In this case, we use the image given by

$$S_{\text{noise},i} = \omega \times \frac{n_i}{y_{\text{radial model}}^{(\text{den})}}, \quad (20)$$

for each Monte Carlo i, to obtain the power spectra $\mathcal{P}_{\text{noise},i}(k_{2\text{D}})$. The mean across all the Monte Carlo realizations provides the expected noise bias $\mathcal{P}_{\text{noise}}(k_{2\text{D}}) = \langle \mathcal{P}_{\text{noise},i}(k_{2\text{D}}) \rangle$ and the noise covariance matrix is given by $C_{\text{noise}} = \langle (\mathcal{P}_{\text{noise},i}(k_{2\text{D}}) - \mathcal{P}_{\text{noise}}(k_{2\text{D}}))(\mathcal{P}_{\text{noise},i}(k_{2\text{D}}) - \mathcal{P}_{\text{noise}}(k_{2\text{D}})) \rangle$. It will be an essential ingredient for the likelihood estimation and power spectrum estimation (Section 5). We note that model uncertainties are not accounted in Equation (20) because they will be treated as a systematic effect.
When two independent data (sub)sets are used to compute the power spectrum, the noise is estimated accordingly. We thus use the cross spectra between two independent noise realizations $S_{\text{noise},i}^{(a)}$ and $S_{\text{noise},i}^{(b)}$. In this case, the mean noise power spectrum is expected to be zero.

The exact same procedure allows us to extract the contribution from the CIB, to obtain $\mathcal{P}_{\text{CIB}}(k_{2\text{D}})$ and $C_{\text{CIB}}$, or any astrophysical contaminant in general.

### 4.3. Summary, discussion, and example

In PITSZI, the power spectrum of the data is computed as the modulus squared of the fast Fourier transform of the quantity $S$ defined in Equation (15), which is then binned in $k_{2\text{D}}$ to obtain $\mathcal{P}_{\text{observed}}(k_{2\text{D}})$. It is also possible to use the cross spectra between two independent data (sub)sets, $S^{(a)}$ and $S^{(b)}$, if available. The weight map accounts for the region of interest and discarded pixels. In practice, the weights are somewhat arbitrary. For instance, one might set $\omega$ to $y_{\text{model}}^{(\text{den})}$ to avoid a dramatic increase in the noise with radius. However, this will be at the cost of obtaining a power spectrum estimate that will be weighted towards the cluster center.

To obtain a mock or an actual model power spectrum that can be compared to the data, one must account for the instrument response function and the weighting scheme. Starting from a 3D pressure fluctuation power spectrum model, PITSZI implements two different approaches. The projection approach computes the model as follows: 1) Equation (13) is used to compute the theoretical 2D power spectrum; 2) Equation (16) is applied to account for the instrument response function; 3) the power spectrum is projected in $k$-space given the sampled mode $m_1, m_2$ and Equation (19) is applied to account for data weighting; 4) the power spectrum is binned in $k_{2\text{D}}$ to obtain $\mathcal{P}_{\text{ICM}}(k_{2\text{D}})$. Alternatively, a brute force approach may be used in which full SZ image Monte Carlo simulations are generated and their power spectra calculated as for the data. This is done by generating mock pressure cubes including the radial profile and the fluctuations; obtaining SZ images via Equation (2); convolving the image with the instrument response function; applying the weights; and finally extracting the 2D power spectrum. As in the case of the noise (Section 4.2.3), the statistics of the power spectra obtained in this manner allow us to determine the mean power spectrum $\mathcal{P}_{\text{ICM}}(k_{2\text{D}}) = \langle \mathcal{P}_{\text{Monte Carlo},i}(k_{2\text{D}}) \rangle$ and its covariance matrix $C_{\text{ICM}}$. The brute force approach is significantly more computationally expensive than the projection approach, and is not used to fit the data.

To illustrate the methodology discussed in this Section, we now compare the MACS J0717.5+3745 NIKA data with three mock SZ maps constructed from different pressure fluctuation models with PITSZI (model 1, model 2 and model 3, see Table 1). These are designed to qualitatively mimic the NIKA data, but we refer to Section 5 for a full discussion of model inference. The radial model is the same for the three cases, and is obtained by fitting an elliptical gNFW model to the NIKA data

$$ \bar{P}(r) = \frac{P_0}{\left(\frac{r}{r_p}\right)^c \left(1 + \left(\frac{r}{r_p}\right)^a\right)^{\frac{b-c}{a}}}. \qquad (21) $$

The slope parameters $a$, $b$ and $c$ are fixed to the morphologically disturbed model of Arnaud et al. (2010), the centre is free to vary, and the mass is fit via a single normalization parame-

<table><caption>**Table 1.** Summary of the test model properties used for illustration.</caption><thead><tr><th></th><th>σ<sub>ρ</sub></th><th>L<sub>inj</sub></th><th>α</th><th>Statistics</th></tr></thead><tbody><tr><td>Test model 1</td><td>0.5</td><td>1 Mpc</td><td>-11/3</td><td>Gaussian</td></tr><tr><td>Test model 2</td><td>0.5</td><td>1 Mpc</td><td>-11/3</td><td>Lognormal</td></tr><tr><td>Test model 3</td><td>0.2</td><td>250 kpc</td><td>-11/3</td><td>Gaussian</td></tr></tbody></table>

ter that also defines the scale radius<fnref n="8" />. The first fluctuation model is built with $\sigma_{\rho} = 0.5$, $L_{\text{inj}} = 1$ Mpc, and $\alpha = -11/3$, with Gaussian pressure fluctuations. The second model has the same power spectrum parameters but with lognormal pressure fluctuations (and we use the same seed to generate random fluctuations). The third model is built with $\sigma_{\rho} = 0.2$, $L_{\text{inj}} = 250$ kpc, and $\alpha = -11/3$, with Gaussian pressure fluctuations. In all cases, $L_{\text{dis}} = 1$ kpc, but the exact value of this parameter is irrelevant here because it lies well below the instrumental resolution.

In Figure 2, we can see that, given the elliptical radial model, SZ fluctuations are observed by eye in the NIKA residual map. The half-difference (Jackknife) map gives a visual comparison of the noise contribution to the data. In Figure 3, the SZ mocks derived from model 1 and 2 visually resemble the data, at least qualitatively. When including the noise Monte Carlo, it would be nearly impossible to distinguish mock images from true data. This indicates that our description of the ICM as a radial model plus isotropic pressure perturbations mimics the real data well, at least given the sensitivity and the scales that are probed with NIKA. We defer discussion of the limitations of the modelling approach to Section 6. Contrary to models 1 and 2, we can already observe that model 3 presents too few fluctuations and its residual is nearly consistent with the half difference map of the true data. Visually, models 1 and 2 are nearly identical, and look statistically similar.

In Figure 4, we show that indeed, the fluctuation statistics (lognormal or Gaussian) cannot be distinguished for the pressure fluctuation models 1 and 2. While the pressure fluctuation distributions present significant differences, when seen as SZ fluctuation models 1 and 2 are nearly compatible, even without including noise in the data, owing to projection effects. The Gaussian versus lognormal models become significantly different only when the fluctuation amplitude is significantly increased. When including noise, it is not possible to distinguish the two within uncertainties. The fluctuations obtained for MACS J0717.5+3745 are compatible with those expected in both models 1 and 2.

In Figure 5, we show the input power spectrum of the pressure fluctuations for the different models. The scales accessible with NIKA roughly correspond to the shaded area, since the signal is progressively filtered for scales larger than the field of view, and on scales smaller than the beam. We thus expect sensitivity to the normalization, the injection scale, and possibly the slope if a sufficient signal-to-noise ratio is available. The dissipation scales, on the other hand, is expected to be far beyond the reach of NIKA in the case of MACS J0717.5+3745. We also compute the power spectrum of the corresponding SZ fluctuations with our different methods and include the different contributions to the data. Here, the weight map, shown in the bottom right of the right panel, was defined so that pixels beyond $\theta_{500}/2 \approx 2$ arcmin were masked and the ICM was weighted

<footnote n="8">⁸ In Equation (21), the characteristic radius $r_p = R_{500}/c_{500}$ and the normalization $P_0$ can be related to the mass assuming that clusters follow the same scaled pressure profile. See Arnaud et al. (2010) for details.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Simulated SZ mocks to be compared with the NIKA data (Figure 2). The smoothing kernel and the contours are the same as Figure 2. *From left to right*, we show a noiseless SZ realization, the same image with noise, the noiseless residual between the mock and the radial model, and the residual including noise. The first, second and third rows correspond to model 1, model 2 and model 3, respectively. See Table 1 for the definition of the models.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison between the pressure and 1D SZ fluctuations probability distributions for model 1 (gaussian) and model 2 (lognormal), and MACS J0717.5+3745 data. *Left*: probability density function of the pressure fluctuations $1 + \delta P / \bar{P}$ for model 1 (gaussian fluctuations) and model 2 (lognormal fluctuations). *Middle*: probability density function of the SZ fluctuations $1 + \frac{\delta y}{y}$ for model 1 and model 2, without noise. *Right*: probability density function of the SZ fluctuations $1 + \frac{\delta y}{y}$ for model 1 and model 2, including noise, and comparison with MACS J0717.5+3745 NIKA data. The data and model were smoothed with a 10 arcsec FWHM gaussian kernel to reduce the noise. The mean distributions and error bars were obtained from the mean and standard deviations, respectively, of 100 Monte Carlo realizations.</figcaption></figure>

(i.e. the $\omega$ term) according to the smooth radial model (with inner slope set to zero) and smoothed with a 20 arcsec FWHM Gaussian kernel to minimize aliasing. We can see that the projection approach matches very well the mean brute force approach, indicating that the approximations discussed in this Section are robust over the full range of scales. The contribution of the CIB is expected to be lower than the noise, but it is still non-negligible. While model 1 is expected to be well detected, and in fact match the data relatively well, model 3 is below the noise and the CIB contributions. Model 2 is not represented because

it is very similar to model 1. The shaded regions provide the standard deviation of the different power spectra, including that of the pressure fluctuation model, which arise from the fact that they intrinsically correspond to stochastic processes. The data show both the auto spectrum of the full data and the cross spectrum of two independent equivalent data subsets. In both cases, a clear excess associated with the injected signal is observed in the range $k \sim 10^{-3} - 3 \times 10^{-3}$ kpc$^{-1}$. The auto-spectrum converges to the expected noise mean at high $k$, where it dominates, while the cross-spectrum is scattered around zero. The auto- and
<figure><img src="image_6.png" /><figcaption>Fig. 5. Input pressure fluctuation power spectrum and its 2D projection. Left: input 3D pressure fluctuations power spectra for model 1, 2, and 3. The grey area roughly corresponds to the observable modes with the NIKA data. Right: 2D SZ fluctuations power spectra for MACS J0717.5+3745 NIKA data, the noise contribution, the CIB contribution, and models 1 and 3 including the projection and brute force approaches. The shaded area corresponds to the measured standard deviation of the respective spectra. For the data, error bars indicate the contribution from the noise and CIB only in black, and account for the intrinsic scatter of model 1 in grey. We show both the auto spectrum of the full data, and the cross spectra of two independent data subsets. For the later, the power spectrum may be negative, as indicate in the legend. The uncertainties due to the noise are $\sqrt{2}$ larger for the cross spectrum than for the auto spectrum.</figcaption></figure>

cross-spectra agree very well on scales where the signal dominates. The noise contribution dominates the uncertainties at high $k$, while the intrinsic model scatter dominates at low $k$. In the case of MACS J0717.5+3745, the CIB contribution to the uncertainties is subdominant at all scales, although its mean contribution is not negligible.

## 5. Pressure profile and pressure fluctuation inference

In Section 4, we discussed how to compute the SZ fluctuations 2D power spectrum, and how to derive the corresponding model given a 3D pressure fluctuation power spectrum and radial model. Here we discuss the methodology implemented in PITSZI to infer constraints on the radial model and the pressure fluctuation power spectrum, and then apply it to MACS J0717.5+3745 data.

### 5.1. Pressure profile

The fitting approach used to constrain the radial pressure profile consists of comparing the model, including projection effects and instrumental effects as described in Section 4.1, to the data. While the ICM triaxiality can be modelled in three dimensions, SZ data alone cannot constrain the ICM shape along the line of sight (see, e.g. Kim et al. 2024, for such an analysis). Therefore, in the fitting process, only a single axis ratio (assuming $q_{int} = q_{min}$) and the rotation angle $\phi_3$ (assuming $\phi_1 = 0$ and $\phi_2 = 90$ deg) are fit. This corresponds to a fit of the ellipticity as projected onto the plane of the sky. At this stage, PITSZI also accounts for the map zero level as a nuisance parameter, since it is generally unconstrained from SZ observations. In the end, PITSZI allows for the fit of any pressure profile model parameters, the cluster centre, the cluster ellipticity, and the map zero level.

Two different algorithms are implemented to constrain the model: sampling of the parameter space via Markov Chain Monte Carlo (MCMC) using the `emcee` package (Foreman-Mackey et al. 2013), and the use of a non-linear least

square fit via the `scipy CURVE_fit` package<fnref n="9" />. The MCMC fit uses a Gaussian log likelihood computed as

$$ \ln \mathcal{L} \propto \sum_{m,n} (D - M)_m (C^{-1})_{m,n} (D - M)_n \quad (22) $$

where $D$ represents the unmasked data map, $M$ the radial model, and $C$ the noise covariance matrix. In the present paper, we neglect pixel-to-pixel correlations and only account for the noise standard deviation in the fit. This is because obtaining the full noise covariance matrix is computationally expensive, and because its use essentially affects the posterior uncertainties on the parameters which are dominated by systematic uncertainties associated to the choice of radial model. In the case where the posterior is well approximated by a multivariate Gaussian distribution, the MCMC and non-linear least square methods give fully consistent results. In the present analysis, we will neglect the statistical uncertainties associated with the fit parameters of the radial model, because they are much smaller than the systematic uncertainties induced by the choice of the radial model.

Indeed, the choice of radial model is known to be key for surface brightness fluctuation analysis (e.g. Dupourqué et al. 2023; Romero et al. 2023) and is responsible for most of the systematic uncertainties. Therefore, we explore the stability of the results by using the following radial models (RM):

1. RM₁: centre fixed to the X-ray peak (roughly consistent with X-ray barycenter) and a spherical profile. The shape of the pressure profile is fixed to the morphologically disturbed model from Arnaud et al. (2010), and the mass $M_{500}$ is the only free parameter of the pressure profile. The model has two free parameters (including the map zero level);

2. RM₂: free centre, spherical profile, $M_{500}$ fitted as the normalization of the morphologically disturbed Arnaud et al. (2010) pressure profile ($2 + 1 + 1 = 4$ parameters);

3. RM₃: free centre, free elliptical shape, $M_{500}$ fitted as the normalization of the Arnaud et al. (2010) morphologically disturbed pressure profile ($2 + 2 + 1 + 1 = 6$ parameters);

<footnote n="9">⁹ All parameters may be fixed, sharp limits can be imposed, and in the case of MCMC, Gaussian priors can be used.</footnote>
<table><caption>**Table 2.** Best-fit parameters obtained for the different radial models.</caption><thead><tr><th></th><th>RA<br>(deg)</th><th>Dec<br>(deg)</th><th>q<sub>int, maj</sub><br>(–)</th><th>φ<sub>3</sub><br>(deg)</th><th>M<sub>500</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>P<sub>0</sub><br>(10<sup>-2</sup> keV cm<sup>-3</sup>)</th><th>r<sub>p</sub><br>(kpc)</th><th>a<br>(–)</th><th>b<br>(–)</th><th>c<br>(–)</th></tr></thead><tbody><tr><td>RM<sub>1</sub></td><td>109.3806</td><td>37.7583</td><td>1</td><td>–</td><td>18.1</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td>RM<sub>2</sub></td><td>109.3872</td><td>37.7530</td><td>1</td><td>–</td><td>21.7</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td>RM<sub>3</sub></td><td>109.3878</td><td>37.7528</td><td>0.77</td><td>23</td><td>26.5</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td>RM<sub>4</sub></td><td>109.3874</td><td>37.7525</td><td>0.74</td><td>24</td><td>–</td><td>11.7</td><td>852</td><td>–</td><td>–</td><td>–</td></tr><tr><td>RM<sub>5</sub></td><td>109.3870</td><td>37.7543</td><td>0.71</td><td>21</td><td>–</td><td>10.5</td><td>463</td><td>6.6</td><td>7.4</td><td>0.06</td></tr></tbody></table>

**Notes.** Uncertainties on the radial model parameters are negligible compared to the systematic effects associated with the choice of the model, and are thus omitted. The results shown here were obtained using the non-linear least square method.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Residual SZ maps between the data and the radial models under consideration. Contours are multiples of $3\sigma$, and the maps have been smoothed with a 15 arcsec FWHM gaussian kernel for visual purposes.</figcaption></figure>

4. RM$_4$: free centre, free elliptical shape, free gNFW pressure profile normalization $P_0$ and scale radius $r_p$ ($2+2+2+1=7$ parameters);

5. RM$_5$: free centre, free elliptical shape, free gNFW pressure profile ($2+2+5+1=10$ parameters).

In the case of elliptical profiles, we stress that the mass determination is driven by the shape and normalization of the pressure profile along the major axis, so that its meaning should be taken with caution (see Section 4.1.2). In the case of RM$_5$, the profile parameters being highly degenerate, we impose a flat prior on the slopes of the gNFW model such that $a \in [0, 5]$, $b \in [2, 8]$, $c \in [0, 2]$ given the results from Planck Collaboration Int. V (2013). On one hand, if the model is too simplistic, the cluster potential well is not well represented, and we would expect to observe excess residual features unrelated to pressure fluctuations. On the other hand, models that are too complex would be expected to overfit the pressure fluctuations. See Appendix C for more discussion. Accordingly, we expect that our baseline model RM$_3$ should provide a fairly good compromise.

In Table 2, we give the best fit results obtained for the five models. In Figure 6, we show the SZ residual maps. As expected, the amount of SZ surface brightness residuals diminishes as a function of increasing model complexity. While fitting for the centre largely improves the residual, fitting for the ellipticity has a smaller impact. The most complex model is able to account for a substantial fraction of the fluctuations, and no large features are observed visually in the residual. In Appendix D, we also show that attempting to correct for the kSZ signal implies larger residuals, in particular for model RM$_5$.

## 5.2. Pressure fluctuations

With the radial model in hand, PITSZI can constrain the pressure fluctuation model parameters given the data. The power spectrum of the data is extracted from the map $S$ given by Equation (15), or from the cross spectra between two independent maps. Data and model are then compared by processing the model following the formalism discussed in Section 4.

Again, two fitting methods have been implemented: the MCMC parameter space sampling, and a non-linear least-squares fit. For the MCMC, the log-likelihood is defined, as for the profile, under the Gaussian approximation$^{10}$

$$ \ln \mathcal{L} \propto \sum_{m,n} (D - M)_m (C^{-1})_{m,n} (D - M)_n, \quad (23) $$

where $D$ corresponds to the power spectrum of the data, $M$ is the model power spectrum, and $C$ is the covariance matrix. The model is given for both methods by the sum of the different contributions to the power spectrum

$$ M = M(k_{2D}) = \mathcal{P}_{\text{ICM}}(k_{2D}) + A_{\text{noise}}\mathcal{P}_{\text{noise}}(k_{2D}) + A_{\text{CIB}}\mathcal{P}_{\text{CIB}}(k_{2D}). \quad (24) $$

Here, the quantity $\mathcal{P}_{\text{ICM}}(k_{2D})$ refers to the ICM pressure model, which depends on the (fixed) pressure profile and on the pressure fluctuation parameters. The nuisance parameters $A_{\text{noise}}$ and $A_{\text{CIB}}$ allow PITSZI to marginalize over the amplitude of the noise and the CIB, to account for possible differences between the predictions and these components in the real data, and to check for degeneracies with the pressure fluctuation model. The power spectrum covariance matrix is given by the sum of the contributions from the noise and the CIB. The contribution from the pressure fluctuation can be added to account for sample variance using a reference spectrum that matches the data. In practice, we first fit the model without accounting for sample variance to obtain a reference model, and then refit the data accounting for the sample variance using the pressure fluctuation model covariance given the best-fit model from the first step. We check that the results are stable after only 1 iteration, given the relatively high signal-to-noise ratio for MACS J0717.5+3745. The sample variance is thus accounted for in the total error budget.

We apply this method to the NIKA data on MACS J0717.5+3745. We fit for the pressure fluctuation

<footnote>$^{10}$ In practice, the probability distribution followed by the data is not strictly Gaussian, especially for the first bins if the number of modes that are averaged to extract the spectrum is low.</footnote>
<table><caption>**Table 3.** Constraints on the pressure fluctuation parameters.</caption><thead><tr><th>Analysis</th><th>Model</th><th>σ<sub>ρ</sub><br>(-)</th><th>L<sub>inj</sub><br>(kpc)</th><th>A<sub>noise</sub><br>(-)</th><th>A<sub>CIB</sub><br>(-)</th></tr></thead><tbody><tr><td>Reference<sup>†</sup></td><td>RM<sub>1</sub></td><td>1.01<sup>+0.08</sup><sub>-0.11</sub></td><td>3128<sup>+1322</sup><sub>-1792</sub></td><td>1.06<sup>+0.04</sup><sub>-0.04</sub></td><td>1.8<sup>+1.8</sup><sub>-1.3</sub></td></tr><tr><td>Reference<sup>†</sup></td><td>RM<sub>2</sub></td><td>0.63<sup>+0.03</sup><sub>-0.03</sub></td><td>859<sup>+63</sup><sub>-62</sub></td><td>1.01<sup>+0.04</sup><sub>-0.04</sub></td><td>2.4<sup>+1.7</sup><sub>-1.6</sub></td></tr><tr><td>Reference<sup>†</sup></td><td>RM<sub>3</sub></td><td>0.58<sup>+0.04</sup><sub>-0.04</sub></td><td>873<sup>+74</sup><sub>-75</sub></td><td>1.04<sup>+0.04</sup><sub>-0.04</sub></td><td>2.4<sup>+1.7</sup><sub>-1.7</sub></td></tr><tr><td>Reference<sup>†</sup></td><td>RM<sub>4</sub></td><td>0.68<sup>+0.07</sup><sub>-0.07</sub></td><td>653<sup>+78</sup><sub>-78</sub></td><td>1.03<sup>+0.04</sup><sub>-0.04</sub></td><td>2.4<sup>+1.7</sup><sub>-1.7</sub></td></tr><tr><td>Reference<sup>†</sup></td><td>RM<sub>5</sub></td><td>0.27<sup>+0.11</sup><sub>-0.11</sub></td><td>1010<sup>+410</sup><sub>-418</sub></td><td>1.05<sup>+0.04</sup><sub>-0.04</sub></td><td>2.5<sup>+1.7</sup><sub>-1.7</sub></td></tr><tr><td>Cross-spectrum</td><td>RM<sub>3</sub></td><td>0.64<sup>+0.04</sup><sub>-0.04</sub></td><td>812<sup>+71</sup><sub>-72</sub></td><td>0 (fixed)</td><td>2.5<sup>+1.7</sup><sub>-1.7</sub></td></tr><tr><td>kSZ template subtracted</td><td>RM<sub>3</sub></td><td>0.70<sup>+0.03</sup><sub>-0.03</sub></td><td>852<sup>+49</sup><sub>-51</sub></td><td>1.03<sup>+0.04</sup><sub>-0.04</sub></td><td>2.4<sup>+1.7</sup><sub>-1.6</sub></td></tr><tr><td>Point sources masked</td><td>RM<sub>3</sub></td><td>0.61<sup>+0.04</sup><sub>-0.04</sub></td><td>851<sup>+72</sup><sub>-69</sub></td><td>1.03<sup>+0.04</sup><sub>-0.04</sub></td><td>2.4<sup>+1.7</sup><sub>-1.6</sub></td></tr><tr><td>No radial model based weight</td><td>RM<sub>3</sub></td><td>0.70<sup>+0.04</sup><sub>-0.04</sub></td><td>822<sup>+62</sup><sub>-62</sub></td><td>0.98<sup>+0.04</sup><sub>-0.04</sub></td><td>2.4<sup>+1.7</sup><sub>-1.7</sub></td></tr><tr><td>MCMC fit</td><td>RM<sub>3</sub></td><td>0.58<sup>+0.03</sup><sub>-0.03</sub></td><td>869<sup>+67</sup><sub>-60</sub></td><td>1.04<sup>+0.04</sup><sub>-0.04</sub></td><td>1.1<sup>+1.7</sup><sub>-0.8</sub></td></tr></tbody></table>

**Notes.** The central value gives the median of the distribution and the uncertainties correspond to the 68% confidence intervals. The reference values († symbol) were obtained with the following analysis choices: auto-spectrum, no kSZ correction, point sources subtracted, radial model based weighting plus 2 arcmin radius region of interest, non-linear least square fit, and projection approach.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Constraints on the fluctuation model obtained for the reference analysis. Left: comparison between the data and model for the 2D SZ fluctuation power spectrum, in the case of radial model RM<sub>3</sub>. The residual shows the offset to the best-fit model relative to the uncertainties, accounting or not for intrinsic model variance (sample variance). Shaded area provide 68% and 95% confidence intervals. Right: constraints on the 3D pressure fluctuation power spectrum for the different radial model, with 68% confidence interval shown. The sample variance is included in the error budget.</figcaption></figure>

normalization σ<sub>ρ</sub>, the injection scale L<sub>inj</sub>, and marginalize over the noise (A<sub>noise</sub>) and CIB (A<sub>CIB</sub>) normalization. Given the scales sampled by the data, it is not feasible to fit for the slope without large degeneracies with other parameters, and we therefore fix it to α = −11/3 (see also Dupourqué et al. 2023, 2024). The impact of this choice on the results is discussed in Appendix E. We use flat priors on parameters, as σ<sub>ρ</sub> ∈ [0, 3], L<sub>inj</sub> ∈ [50, 5000], A<sub>noise</sub> ∈ [0, 5] and A<sub>CIB</sub> ∈ [0, 5]. However, this only affects the results for the CIB normalization since it is poorly constrained, so that we marginalize over its uncertain amplitude within a large range. We use RM<sub>3</sub> as our baseline radial model, but since this choice is the main source of uncertainty, we also provide the results for the other models. The pressure fluctuations are assumed to follow Gaussian statistics, and trace the velocity field. In terms of methodology, we make the following choices as a reference: the auto-spectrum is used (we note that when using the cross-spectrum between two data subsets, the noise amplitude is fixed to zero since the two sets

are independent); no kSZ correction is performed; the point sources are subtracted but not masked; the weight map (ω) is defined as being proportional to the radial model Compton parameter map multiplied by the binary mask (2 arcmin radius region of interest, plus point sources if relevant) and smoothed with a 20 arcsec FWHM Gaussian kernel, giving more weight to the inner part of the cluster; the non-linear least square fit is used; the projection approach is used. Nevertheless, we provide the results when varying all of these choices. For the non linear least square fit, the best fit parameters and their covariance matrix is used to produce multivariate Gaussian distributions as the constraint in the parameter space, accounting for sharp limits on the parameters.

Numerical results of the fits are listed in Table 3. In Figure 7, we show the constraints in terms of the 2D power spectrum and the corresponding constraints in terms of the 3D power spectrum, for our reference analysis choice RM<sub>3</sub>. Figure 8 provides the posterior constraints on the parameter space, focusing only
<figure><img src="image_9.png" /><figcaption>Fig. 8. Top: constraints in the parameter space for models RM$_{2}$, RM$_{3}$ and RM$_{4}$. Other extreme models are not shown for clarity. Contours correspond to the 68 and 95% confidence intervals. Bottom: comparison between the constraints obtained on RM$_{3}$ for the full MCMC sampling and the non linear least square approximation methods, in the case of RM$_{3}$.</figcaption></figure>

on the non-extreme models RM$_{2}$, RM$_{3}$ and RM$_{4}$ for clarity. As expected, the value of $\sigma_{\rho}$ and $L_{\text{inj}}$ tends to decrease with model complexity, from RM$_{1}$ to RM$_{5}$, although the value of $L_{\text{inj}}$ is nearly unconstrained for RM$_{5}$, which explains its high central value. Excluding the two extremes (RM$_{1}$ and RM$_{5}$), the results are relatively stable with respect to the choice of the radial model. This can also be observed in the 3D power spectrum constraint and the parameter space constraint.

The statistical uncertainty on the normalization $\sigma_{\rho}$ is of the order of 10%, while the injection scale is also well con-

strained, with similar precision. If the injection scale is poorly constrained, for instance because the power spectrum peak falls beyond the scales accessible from the data (RM$_{1}$) or because the fluctuations have been absorbed in the radial model (RM$_{5}$), the uncertainties increase significantly. The systematic uncertainty due to the radial model is much larger than the statistical uncertainty when including all radial models (about 50%), but these uncertainties are comparable when excluding the two extreme radial models (see also Appendix C for an attempt to quantify the systematics induced by these models). We also show that the analysis choices have a minor impact on the results as they are all nearly compatible within the statistical uncertainties. The main changes are observed in the case of the kSZ template subtraction and the change of the weighting scheme, which effectively slightly changes the power spectrum extraction region under consideration. They both imply a slightly higher pressure fluctuation amplitude (by about $3\sigma$). Given the relatively high signal-to-noise detection of the pressure fluctuations, the parameter posterior distribution is nearly Gaussian, so that the MCMC fit and the non-linear least squares method are fully consistent for $\sigma_{\rho}$, $L_{\text{inj}}$ and $A_{\text{noise}}$. The CIB normalization is compatible with unity, but is poorly constrained. As its posterior distribution significantly deviates from that of a Gaussian, the CIB constraints given in Table 3 for the non-linear least squares and MCMC fits differ. The noise normalization is well constrained and is compatible with unity, indicating that it describes the data well.

## 6. Implications for the non-thermal ICM physics

The constraints on the pressure fluctuation power spectrum give information on the non-thermal physics of the ICM, assuming the fluctuations are directly related to the turbulence. In this Section, we explore the physical consequences of the measurements obtained above, compare the results with other work, and discuss limitations and caveats of the method.

In PITSZI, the physical interpretation is implemented in the class Physics (see Figure 1), essentially as a list of functions that encode the equations described in detail below, and the results from the literature used here. As will be discussed further later, we stress that the results derived here are subject to several assumptions and are limited by systematic uncertainties, in particular due to the choice of the radial model.

### 6.1. From the power spectrum to non-thermal physics

The power spectrum of the pressure fluctuations is related to the turbulent 3D Mach number, as characterized by numerical simulations (e.g. Gaspari et al. 2014; Mohapatra et al. 2020, 2021, 2022; Zhuravleva et al. 2023). For instance, assuming adiabatic perturbations, the relation obtained by Gaspari et al. (2014) for density perturbations becomes, when expressed as pressure perturbations

$$ \mathcal{M}_{3D} = \frac{4}{\gamma} A_{\delta P / \bar{P}}(k_{\text{peak}}) \left( \frac{L_{\text{inj}}}{500 \text{ kpc}} \right)^{-s}, \qquad (25) $$

where $A_{\delta P / \bar{P}}(k_{\text{peak}})$ is the peak of the characteristic amplitude of the pressure fluctuation power spectrum (Equation (6)). The slope $s \approx 0.25$ gives a small correction as a function of injection scale. In contrast, Zhuravleva et al. (2023) provide a scaling of the Mach number with the integrated power spectrum ($\propto \sigma_{\rho}$), for pressure perturbations, accounting for a radial dependence and for different dynamical states. These relations read

$$ \frac{\mathcal{M}_{3D}}{\sqrt{3}} = a + b \times \frac{r}{R_{500}} + c \times \frac{\delta\xi}{\xi}, \qquad (26) $$
<table><caption>**Table 4.** Non-thermal ICM physics derived from our reference analysis.</caption><thead><tr><th>Model</th><th><i>M</i><sub>3D</sub><br>(-)</th><th>σ<sub>v</sub><br>(km/s)</th><th><i>P</i><sub>kin</sub>/<i>P</i><sub>kin</sub>+<i>P</i><sub>th</sub><br>(-)</th><th><i>b</i><sub>HSE,Δ=2500</sub><sup>(†)</sup><br>(-)</th><th><i>b</i><sub>HSE,Δ=500</sub><sup>(†)</sup><br>(-)</th></tr></thead><tbody><tr><td>RM<sub>1</sub></td><td>0.82<sup>+0.16</sup><sub>-0.16</sub> (+0.04<br>-0.06 stat.)</td><td>1618<sup>+334</sup><sub>-323</sub> (+74<br>-115 stat.)</td><td>0.27<sup>+0.08</sup><sub>-0.08</sub> (+0.02<br>-0.03 stat.)</td><td>0.33<sup>+0.07</sup><sub>-0.06</sub> (+0.007<br>-0.010 stat.)</td><td>0.41<sup>+0.05</sup><sub>-0.05</sub> (+0.005<br>-0.008 stat.)</td></tr><tr><td>RM<sub>2</sub></td><td>0.62<sup>+0.11</sup><sub>-0.11</sub> (+0.02<br>-0.02 stat.)</td><td>1221<sup>+220</sup><sub>-230</sub> (+41<br>-45 stat.)</td><td>0.17<sup>+0.05</sup><sub>-0.05</sub> (+0.01<br>-0.01 stat.)</td><td>0.30<sup>+0.06</sup><sub>-0.06</sub> (+0.004<br>-0.004 stat.)</td><td>0.38<sup>+0.05</sup><sub>-0.05</sub> (+0.003<br>-0.003 stat.)</td></tr><tr><td>RM<sub>3</sub></td><td>0.59<sup>+0.20</sup><sub>-0.19</sub> (+0.02<br>-0.03 stat.)</td><td>1176<sup>+404</sup><sub>-385</sub> (+43<br>-52 stat.)</td><td>0.16<sup>+0.10</sup><sub>-0.08</sub> (+0.01<br>-0.01 stat.)</td><td>0.29<sup>+0.07</sup><sub>-0.06</sub> (+0.004<br>-0.005 stat.)</td><td>0.38<sup>+0.05</sup><sub>-0.05</sub> (+0.003<br>-0.004 stat.)</td></tr><tr><td>RM<sub>4</sub></td><td>0.64<sup>+0.21</sup><sub>-0.21</sub> (+0.04<br>-0.05 stat.)</td><td>1261<sup>+421</sup><sub>-413</sub> (+87<br>-90 stat.)</td><td>0.18<sup>+0.10</sup><sub>-0.09</sub> (+0.02<br>-0.02 stat.)</td><td>0.30<sup>+0.06</sup><sub>-0.06</sub> (+0.008<br>-0.009 stat.)</td><td>0.39<sup>+0.05</sup><sub>-0.05</sub> (+0.006<br>-0.007 stat.)</td></tr><tr><td>RM<sub>5</sub></td><td>0.28<sup>+0.17</sup><sub>-0.14</sub> (+0.10<br>-0.11 stat.)</td><td>545<sup>+332</sup><sub>-273</sub> (+205<br>-226 stat.)</td><td>0.04<sup>+0.06</sup><sub>-0.03</sub> (+0.04<br>-0.03 stat.)</td><td>0.24<sup>+0.07</sup><sub>-0.06</sub> (+0.021<br>-0.023 stat.)</td><td>0.34<sup>+0.06</sup><sub>-0.05</sub> (+0.017<br>-0.021 stat.)</td></tr></tbody></table>

**Notes.** The reported values correspond to the median of the distribution, and the uncertainties account for both the statistical uncertainties in the power spectrum measurement and the scatter in the relation derived by Zhuravleva et al. (2023). The uncertainties obtained when ignoring the scatter are given in parenthesis. (†) These results were obtained with the scaling of Equation (26) (sum) but no significant change is observed when using Equation (27) (product) instead.

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Constraints on the non-thermal ICM parameters. *From left to right:* constraints on the 3D Mach number, the velocity dispersion, the kinetic to total pressure, and the hydrostatic mass bias at $R_{500}$. The constraints for the different radial models are shown with different colours, as indicated in the legend. The top panels provide the constraints obtained without accounting for the scatter in the relation by Zhuravleva et al. (2023), while in the bottom panels this scatter is accounted for.</figcaption></figure>

or

$$ \frac{M_{3D}}{\sqrt{3}} = a \times \left(\frac{r}{R_{500}}\right)^b \times \left(\frac{\delta\xi}{\xi}\right)^c, \qquad (27) $$

where *a*, *b* and *c* are coefficients that depend on the dynamical state and the choice of scaling relation (see their Table 1 and their Fig. 6, for the pressure). Here, we have converted our normalization to the natural logarithm based normalization used by

Zhuravleva et al. (2023), $\sigma_{\text{ln}} = \sqrt{\log((1 + \sqrt{1 + 4\sigma_{\text{P}}^2})/2)}$, and

we have used their notation, $\delta\xi/\xi = 2\sqrt{2\log(2)/\log(10)}\sigma_{\text{ln}}$. Applying the results from Zhuravleva et al. (2023), ignoring the radial dependence and accounting for the ellipticity of our models (their Fig. 6), we obtain the 3D Mach number values listed in Table 4; Figure 9 shows the corresponding probability distribution of the Mach number. For the intermediate model RM3, $M_{3D} \sim 0.6$. Error bars were computed by propagating the uncertainty on the power spectrum parameters and including the scatter from the relation by Zhuravleva et al. (2023) via Monte Carlo realizations. As we can see, the scatter in the relation dominates

over the statistical error, but it remains smaller than the systematic effect due to the choice of the radial model.

The 3D Mach number is related to the velocity dispersion $\sigma_v$ via

$$ M_{3D} = \frac{\sigma_v}{c_s}, \qquad (28) $$

where $c_s$ is the sound speed, which depends on the gas temperature $T_{\text{gas}}$ as

$$ c_s = \sqrt{\frac{\gamma k_B T_{\text{gas}}}{\mu_{\text{gas}} m_p}}. \qquad (29) $$

Here $\gamma = 5/3$ is the gas adiabatic index, $k_B$ the Boltzmann constant, $\mu_{\text{gas}} \approx 0.61$ the gas mean molecular weight and $m_p$ the mass of the proton. Here we assume a mean temperature $T_{\text{gas}} = 15$ keV (Adam et al. 2017a). This implies that $c_s = 1981$ km/s. We list the corresponding values in Table 4 and show the probability distribution of this parameter in Figure 9. For the intermediate model RM3, we obtain $\sigma_v \sim 1200$ km/s.
The ratio between kinetic and thermal ICM pressure is related to the Mach number as

$$X = \frac{P_{\text{kin}}}{P_{\text{th}}} = \frac{1}{2}\gamma(\gamma - 1)M_{\text{3D}}^2. \quad (30)$$

We use this equation to derive the kinetic to kinetic plus thermal pressure ratio. Our findings are listed in Table 4, and are displayed in Fig. 9. We obtain $P_{\text{kin}}/P_{\text{th+kin}} \sim 0.16$ for RM$_3$.

The hydrostatic equilibrium assumption is often used to obtain estimates of the mass of galaxy cluster via (Pratt et al. 2019)

$$M_{\text{HSE}}(<r) = -\frac{r^2}{G\rho_{\text{gas}}(r)} \frac{dP_{\text{th}}(r)}{dr}, \quad (31)$$

using the gas density, $\rho_{\text{gas}}$, and thermal pressure estimates, $P_{\text{th}}$, derived from SZ or X-ray observations, with G being Newton's constant. However, turbulent gas motions will contribute to balance the gravitational field and will induce a bias in the hydrostatic mass if they are non-negligible. Turbulent motions expected to dominate the hydrostatic mass bias (Pratt et al. 2019). Thanks to Equation (30), we can compute a corrected mass profile, $M_{\text{corr}}$, by substituting the thermal pressure with the thermal plus kinetic pressure: $P_{\text{th}}(r) \rightarrow P_{\text{th}}(r) + P_{\text{kin}}(r) = (1 + X(r)) P_{\text{th}}(r)$. The mass bias profile is then obtained from

$$b_{\text{HSE}} = 1 - \frac{M_{\text{HSE}}(<r)}{M_{\text{corr}}(<r)}. \quad (32)$$

The mass bias at an overdensity $\Delta$ is given by $b_{\text{HSE},\Delta} = 1 - \frac{M_{\text{HSE}}(R_{\Delta,\text{HSE}})}{M_{\text{corr}}(R_{\Delta,\text{corr}})}$, where $R_{\Delta}$ is estimated by integrating the mass profile such that the enclosed density within $R_{\Delta}$ reaches $\Delta$ times the critical density of the Universe. In the following, we perform our calculations by using as a reference the gas density profile obtained from XMM-Newton observation used in Adam et al. (2017a) and the pressure profile of model RM$_1$, as it is centered on the X-ray peak. The kinetic pressure correction, X, is obtained via Equation (30) using the Mach number derived from the scaling by Zhuravleva et al. (2023) given in Equation (26), as a function of radius. The statistical uncertainties arising from the pressure fluctuation power spectrum are also accounted for, but these are subdominant. Our estimates of the kinetic pressure correction, except for model RM$_1$, use radial pressure models for which the centre is not fully consistent with that for the thermal model (density and pressure). However, the results obtained for the hydrostatic mass bias are barely affected by the exact choice of thermal model. This is because to first order the bias corresponds to a relative difference with respect to the perfect hydrostatic equilibrium model, such that changes in the exact reference model are not important, at least compared to other uncertainties. Our findings are listed in Table 4 and shown in Figure 9. We obtain $b_{\text{HSE},\Delta=2500} \sim 0.3$ and $b_{\text{HSE},\Delta=500} \sim 0.4$.

## 6.2. Discussions: Comparison with previous work

Recently, Heinrich et al. (2024) used Chandra data to infer the ICM turbulence from X-ray surface brightness density fluctuations. In their sample of 80 clusters, they found the highest characteristic velocities for MACS J0717.5+3745. For this cluster, they obtain measurable density perturbations in the range of scale of about $k \in [2, 10]$ Mpc$^{-1}$, which presents a good overlap with our measurement, as we have a significant signal at scales of about $k \in [0.8, 4]$ Mpc$^{-1}$. They obtained a Mach number $M_{3D} \simeq 0.47$, a velocity dispersion in the range $\sim 1000-1300$ km/s, and

a kinetic to total pressure fraction $P_{\text{kin}}/(P_{\text{kin}} + P_{\text{th}})$ in the range $\sim 7.8-18.9\%$. These results are in excellent qualitative agreement with our findings (Table 4), indicating that we are indeed probing the same physical process with these SZ observations, and that systematic effects are limited. In fact, combining the SZ pressure and X-ray density fluctuations, it is possible to constrain the thermodynamic nature of the turbulence itself. The pressure ($P$) and density ($\rho$) perturbations are related via

$$\frac{\delta P}{\bar{P}} = \Gamma \frac{\delta \rho}{\bar{\rho}}, \quad (33)$$

with $\Gamma = [0, 1, \frac{5}{3}]$ for isobaric, isothermal, and adiabatic fluctuations, respectively. Heinrich et al. (2024) measured $(\delta\rho/\bar{\rho})_k \equiv A_{\delta\rho/\bar{\rho}}(k_{3D}) = \sqrt{4\pi k_{3D}^3 P_{\rho,3D}(k_{3D})}$ at $R_{2500} = 507$ kpc for MACS J0717.5+3745, obtaining $A_{\delta\rho/\bar{\rho}}(1/R_{2500}) = 0.329^{+0.164}_{-0.110}$. Combined with our measurement at the same scale, $A_{\delta\rho/\bar{\rho}}(1/R_{2500}) = [0.47^{+0.04}_{-0.02}, 0.39^{+0.02}_{-0.02}, 0.37^{+0.02}_{-0.02}, 0.40^{+0.03}_{-0.04}, 0.15^{+0.06}_{-0.06}]$ we estimate $\Gamma = [1.43^{+0.72}_{-0.50}, 1.18^{+0.59}_{-0.40}, 1.12^{+0.56}_{-0.38}, 1.22^{+0.62}_{-0.42}, 0.46^{+0.29}_{-0.25}]$ for models RM$_1$ to RM$_5$, respectively. These values agree with those expected for adiabatic or isothermal perturbations, but exclude isobaric perturbations (except for RM$_5$, which sits between isobaric and isothermal). Nevertheless, we note that the two analyses are not necessarily self-consistent since the radial model may be incompatible (e.g. the cluster centre differs). In the future, a joint analysis of SZ and X-ray perturbations at the surface brightness level, with joint radial modelling, may be effective in better constraining the nature of the ICM fluctuations.

We also compared the velocity dispersion obtained via the pressure fluctuations measurement to that obtained directly using the kinetic SZ effect (Sayers et al. 2013; Adam et al. 2017b). We used the ICM line of sight velocity maps derived by Adam et al. (2017b) under the assumptions of their two ICM density models F1 and F2. The image standard deviation, extracted within their mask, which roughly matches our region of interest, gives $\sigma_{kSZ} \simeq 1500$ km/s and $\sigma_{kSZ} \simeq 1300$ for their models F1 and F2, respectively. Although the comparison is limited because of projection effects, noise contamination in the kSZ velocity image, and because the kSZ is only sensitive to the line of sight velocity component, the two measurements compare well. This could indicate that the velocity measured via kSZ, which is essentially dominated by the bulk velocity of the sub-cluster components, corresponds well to the injection of turbulence on the largest scales, that will eventually cascade and dissipate at smaller scales.

The kinetic to thermal pressure support in galaxy clusters has been predicted using numerical simulations (Shaw et al. 2010; Battaglia et al. 2012; Nelson et al. 2014). Using the relation between Mach number and the pressure fluctuation power spectrum obtained by Zhuravleva et al. (2023), we can estimate this quantity via Equation (30) and compare our results with expectations. The comparison is shown in Figure 10, where we account for both statistical uncertainties in the power spectrum constraint and the scatter in the relation by Zhuravleva et al. (2023) that we use. As we can see, our estimates compare well with expectations, although we are limited by the uncertainty in the shape of the radial dependence both in the predictions from numerical simulations and in the scaling relation that we use. Apart from model RM$_5$, our results are generally slightly above the predictions, which might be anticipated since MACS J0717.5+3745 is one of the most dynamically active clusters in the Universe.
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Kinetic to thermal plus kinetic pressure ratio and hydrostatic mass bias as a function of radius. Constraints are obtained using the scaling relation of Zhuravleva et al. (2023) given in Equation (26) (sum, left panel) and Equation (27) (product, right panel). Error contours account for the statistical uncertainties in the power spectrum measures and the scatter in the relation by Zhuravleva et al. (2023). Top: kinetic to thermal plus kinetic pressure ratio constraint and comparison to expectations from numerical simulations (Shaw et al. 2010; Battaglia et al. 2012; Nelson et al. 2014). Here we assumed $M_{500} = 26.5 \times 10^{14} M_{\odot}$ according to the result from model RM$_3$, but this choice does not significantly affect our results. Bottom: hydrostatic mass bias profile computed as $b \equiv 1 - \frac{M_{HSE}}{M_{corr}}$. The results from RM$_2$ and RM$_4$ were omitted for clarity, but they lie between models RM$_1$ and RM$_3$, and RM$_3$ and RM$_5$, respectively.</figcaption></figure>

In Figure 10, we also derived a constraint on the hydrostatic mass bias profile, for which we observed a strong radial shape dependence on the choice of the scaling relation used to obtain the constraint. The two converge at $r \gtrsim 300$ kpc and we obtain a bias of ~0.15–0.3. When comparing directly the recovered masses at overdensity $\Delta$, this increases to $b_{HSE,\Delta=2500} \sim 0.3$ and $b_{HSE,\Delta=500} \sim 0.4$. Such a value corresponds to that generally assumed to reconcile CMB cosmological constraints with cluster counts from Planck (e.g. Planck Collaboration XX 2014). However, it is obtained for a single cluster exhibiting extreme dynamical activity, that is not representative of the full population. Application of the same analysis to other clusters is likely to lead to much lower values for the bias. Moreover, these results are affected by strong assumptions and limitations, which we discuss next, and which prevent us from drawing any strong conclusions.

## 6.3. Discussion: Limitations, approximations, and systematic effects

The results presented in this paper are intrinsically limited by choices in the analysis and systematic effects, in addition to assumptions regarding the nature of the signal that

we are actually probing. It should again be stressed that MACS J0717.5+3745 is an extremely complex system. While it is an excellent target to develop the methodology presented in this paper, the signal that we are extracting is likely to be affected not only by turbulent motions, but also by bulk motions, clumping, or a residual kSZ signal.

Considering MACS J0717.5+3745 as a single cluster is also a strong approximation: it could be thought as a single object with strong internal activity and strong pressure fluctuations, as we did, or could be considered as a set of still nearly independent sub-clusters. Consequently, we could wonder what part of the signal should be included in the hydrostatic component and what part should go to the fluctuations. This question reflects in the issue of the choice of the radial model, which has been intensively discussed in the literature. It is an obvious issue for MACS J0717.5+3745, but it should affect all clusters to some extent, not only strong mergers (e.g. Romero et al. 2023, for discussion). Here we have attempted to estimate the associated uncertainty by applying a number of radial models and bracketing our results with two extreme models. Despite these limitations, our results provide a step forward in the attempt of characterizing the non-thermal physics of the ICM from resolved SZ data, and we can still compare our finding with alternative observations.
The model which we use to describe the power spectrum of the pressure fluctuations consists in a single power law that is truncated with an exponential cutoff at the injection scale and at the dissipation scale. However, the true power spectrum is likely to be more complicated, especially in a merging cluster with multiple components such as MACS J0717.5+3745. Indeed, different physical processes could inject turbulence at different scales, with different power spectra, and in different regions of the cluster. The present model will therefore encapsulate all these contributions into a single global component. The large-scale cutoff may also differ from an exponential. Although it would be very difficult to constrain its shape from the data, because uncertainties are strongly increasing at these scales, this choice could affect the results. We expect that a sharper (softer) cutoff would lead to a smaller (larger) Mach number since the latter scales with the integral of the power spectrum in our analysis. The choice of the power spectrum model is likely to be one of the strongest assumptions of the present work.

In the present paper, we have implemented an analysis framework based on the projection of the 3D pressure fluctuation power spectrum to the 2D SZ fluctuation power spectrum. This involves a linear scaling between the 2D and the 3D spectra, which is completely accurate only for a flat surface brightness and at scales that are sufficiently small, particularly compared to the injection scale. While we have shown that this method leads to only ≤10% bias on the model predictions in our context, future work will likely benefit from avoiding such approximations to obtain a more reliable model in any configuration and at all scales. This could be done via the use of simulation based inference techniques, as applied in (e.g. Dupourqué et al. 2023, 2024).

The pressure fluctuation power spectrum measure is directly obtained from SZ data. When deriving the non-thermal ICM physical properties in Section 6, we relied on scaling relations, calibrated from numerical simulations, to connect the pressure fluctuation power spectrum to the turbulence via the 3D Mach number. Here we have focused on the relations calibrated by Zhuravleva et al. (2023), but testing the implications of different relations based on different simulations could be valuable. We also noted that the scatter in the relation between the pressure fluctuation power spectrum and the 3D Mach number was one of the main sources of uncertainty when extracting non-thermal information on the ICM.

In the present work, we did not attempt to measure the evolution of the power spectrum in different radial bins because we were already limited in terms of the spatial scales accessible with the NIKA data. This implied that our measurement was only used to fix the normalization of the 3D Mach number, and we relied on the radial trend given by the scaling relations from Zhuravleva et al. (2023) when deriving the non-thermal radial profiles. In the future, the NIKA2 data should allow us to probe larger scales, owing to the larger field of view of its camera (Adam et al. 2018a), thus allowing us measure pressure fluctuations further out of the cluster core and split the data into different radial bins to directly constrain the radial evolution from the data.

We used a fixed mean gas temperature across the cluster when extracting the turbulent velocity dispersion. However, MACS J0717.5+3745 contains complex temperature structure which significantly varies over the core region (Adam et al. 2017a). This implies that even when using a fixed Mach number, the velocity dispersion should present significant spatial variations, which we did not consider here.

Finally, and as already mentioned in Section 6.2, the combination of X-ray and SZ fluctuation measures is in principle very

valuable to constrain the nature of the fluctuations. However, we only used the measured constraints on the power spectra and combined pressure and density fluctuation estimates a posteriori. In addition to the coherence of the radial model that is not guaranteed here, such a combination does not allow one to truly test that the measured fluctuations arise from the same physical origin. In the future, applying cross correlation techniques between the different observables should be efficient in testing the true nature of the fluctuations, and in limiting systematic effects that are not correlated between SZ and X-ray observables.

## 7. Summary and conclusions

In this paper, we presented new developments made to measure the pressure fluctuation power spectrum of the diffuse gas in galaxy clusters using resolved SZ data. This led to the construction of the PITSZI software suite, which is publicly available. PITSZI allows one to model the ICM pressure distribution, accounting for the radial (possibly triaxial) component and the fluctuations (gaussian or lognormal) described by a power spectrum. The software can be used to produce mock data and to constrain SZ pressure fluctuations via different fitting methods. PITSZI has been developed and extensively tested with NIKA and NIKA2 data, but it can in principle work with any SZ imaging data (e.g. Planck, SPT, ACT, Bolocam, MUSTANG, MUSTANG2, NIKA, NIKA2).

The code was applied to the triple merging cluster MACS J0717.5+3745. We summarize our main findings below.

- The methodological framework developed in PITSZI improves over previous analyses. In particular, the power spectrum model prediction was shown to be accurate to within less than 10% for realistic ICM models. Additionally, PITSZI allows us directly to constrain the 3D pressure fluctuation power spectrum parameters using a Bayesian approach.

- We reported a high significance detection of pressure fluctuations within ~$\theta_{500}/2$. By modelling and fitting the power spectrum with an exponential cutoff powerlaw with a single injection scale, we obtained a large amplitude, $\sigma_{\rho} \approx 0.6$. The injection scale, $L_{\text{inj}} \approx 800$ kpc corresponds to expectations from merger driven energy injection. Given the scales sampled by the data and the sensitivity that was achieved, it was not possible to constrain the slope, which was fixed to the canonical value of $\alpha = -11/3$.

- Because of projection effects, the instrument response and the noise in the data, it was not possible to distinguish Gaussian from lognormal perturbations despite the large amplitude of the fluctuations. As projection effects dominate, we anticipate that it will be extremely challenging to differentiate the two, even with improved observations. This result contrasts with the measurement reported in Khatri & Gaspari (2016) for the Coma cluster with Planck.

- The origin and nature of the detected fluctuations may be questionable in such a complex merging system, particularly given the presence of a known kSZ signal, and the fact that the choice of the radial (smooth) model plays a crucial role regarding this issue. We have provided a first attempt to quantify the induced systematic effects. We found that fitting a radial model to the data filters the fluctuations by a factor that increases with increasing model complexity, reaching 80% for the most complex test models. On the other hand, a model that is too simplistic was found to induce spurious fluctuations that could overshoot the true signal, even for large fluctuations for the most simple model. These
systematic effect were found to act mainly on scales $k \le 1/L_{\text{inj}}$. By testing radial models with different complexity, we have attempted to bracket the systematics linked to the choice of radial model, but we conclude that this issue is still the most critical when performing quantitative fluctuation analyses.

- Assuming that the measured pressure fluctuations are related to turbulence, the estimated power spectrum of MACS J0717.5+3745 implies a 3D Mach number of about $M_{3D} \sim 0.6$. This corresponds to a kinetic to thermal pressure fraction $P_{\text{kin}}/P_{\text{th+kin}} \sim 0.2$. Given the cluster mean temperature, this implies a very large velocity dispersion $\sigma_v \sim 1200$ km/s. Our measurement aligns very well with previous constraints from X-ray surface brightness fluctuations based on the ICM density. This may indicate that the ICM perturbations are likely adiabatic in nature. Interestingly, we note that the turbulent velocity dispersion constraint also qualitatively matches the direct measurement obtained from kinetic SZ measures (Sayers et al. 2013; Adam et al. 2017b).

- Assuming that turbulence is the main driver of the hydrostatic mass bias, we derive a value of $b_{\text{HSE}} \sim 0.3-0.4$, depending on radius and on the considered model.

Future high-resolution X-ray spectroscopic data will allow us to calibrate the relation between the pressure fluctuation spectrum and the true turbulent velocity spectrum using real data. As high-resolution X-ray spectroscopy will remain much more expensive than deep SZ imaging, especially at high redshift, pressure fluctuation analyses are likely to become an important tool for investigating ICM turbulence. In fact, the methodology presented in this paper should be applied to any other cluster sample. PITSZI is in continuous development. In particular, the implementation of simulation based inference methods (Dupourqué et al. 2023, 2024) are ongoing, with the aim of accounting for sample variance in a more efficient way, which is crucial when dealing with cluster samples. In addition, joint constraints from SZ plus X-ray images are also being considered, and could provide further interesting insights into the nature of the perturbations at play.

## Data availability

The data used in this work are publicly available at https://lpsc.in2p3.fr/NIKA2LPSZ/nika2sz.release.php. See Adam et al. (2018b) for the corresponding article. The PITSZI code, together with several example notebooks, is available at https://github.com/remi-adam/pitszi.

**Acknowledgements.** We are thankful to the anonymous referee for useful comments, which helped improve the quality of the paper. This work is based on observations carried out under project number 237-13 and 222-14 with the NIKA camera at the IRAM 30 m Telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This work was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence of Université Côte d'Azur under reference number ANR-15-IDEX-01. This work has been supported by the French government, through the UCA¹,E.D.¹. Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01. This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. EP and NC acknowledge support from the French Agence Nationale de la Recherche (ANR), under grant ANR-22-CE31-0010. GWP acknowledges support from CNES, the French space agency. We would like to thank the IRAM staff for their support during the campaigns. The NIKA dilution cryostat has been designed and built at the Institut Néel. In particular, we acknowledge the crucial contribution of the Cryogenics Group, and in particular Gregory Garde, Henri Rodenas, Jean Paul Leggeri, Philippe Camus. This work has been partially funded by the Foundation Nanoscience Grenoble, the LabEx FOCUS ANR-11-LABX-0013 and the ANR under the

contracts 'MKIDS' and 'NIKA'. This work has benefited from the support of the European Research Council Advanced Grants ORISTARS and M2C under the European Union's Seventh Framework Programme (Grant Agreement nos. 291294 and 340519). We thank Instituto de Astrofísica de Andalucía (IAA-CSIC), Centro de Supercomputación de Galicia (CESGA) and the Spanish academic and research network (RedIRIS) in Spain for hosting Uchuu DR1, DR2 and DR3 in the Skies & Universes site for cosmological simulations. The Uchuu simulations were carried out on Aterui II supercomputer at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan, and the K computer at the RIKEN Advanced Institute for Computational Science. The Uchuu Data Releases efforts have made use of the skunIAA_RedIRIS and skun6IAA computer facilities managed by the IAA-CSIC in Spain (MICINN EU-Feder grant EQC2018-004366-P). This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013), in addition to NumPy (van der Walt et al. 2011). SciPy (Jones et al. 2001), and Ipython (Pérez & Granger 2007). Figures were generated using Matplotlib (Hunter 2007).

## References

Adam, R., Comis, B., Macías-Pérez, J. F., et al. 2014, A&A, 569, A66

Adam, R., Comis, B., Macías-Pérez, J. F., et al. 2015, A&A, 576, A12

Adam, R., Comis, B., Bartalucci, I., et al. 2016, A&A, 586, A122

Adam, R., Arnaud, M., Bartalucci, I., et al. 2017a, A&A, 606, A64

Adam, R., Bartalucci, I., Pratt, G. W., et al. 2017b, A&A, 598, A115

Adam, R., Adane, A., Ade, P. A. R., et al. 2018a, A&A, 609, A115

Adam, R., Hahn, O., Ruppin, F., et al. 2018b, A&A, 614, A118

Adam, R., Goksu, H., Leingártner-Goth, A., et al. 2020, A&A, 644, A70

Adam, R., Goksu, H., Brown, S., Rudnick, L., & Ferrari, C. 2021, A&A, 648, A60

Angelinielli, M., Vazza, F., Giocoli, C., et al. 2020, MNRAS, 495, 864

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

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

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Barret, D., Decourchelle, A., Fabian, A., et al. 2020, Astronomische Nachrichten, 341, 224

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

Beduzzi, L., Vazza, F., Brunetti, G., et al. 2023, A&A, 678, L8

Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89

Béthermin, M., Gkogkou, A., Van Cuyck, M., et al. 2022, A&A, 667, A156

Birkinshaw, M. 1999, Phys. Rep., 310, 97

Brunetti, G., & Jones, T. W. 2014, International Journal of Modern Physics D, 23, 1430007P

Cassano, R., Ettori, S., Giacintucci, S., et al. 2010, ApJ, 721, L82

Catalano, A., Calvo, M., Ponthieu, N., et al. 2014, A&A, 569, A9

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

Clerc, N., Cucchetti, E., Pointecouteau, E., & Peille, P. 2019, A&A, 629, A143

Donnert, J., Vazza, F., Brüggen, M., & ZuHone, J. 2018, Space Sci. Rev., 214, 122

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

Dupourqué, S., Clerc, N., Pointecouteau, E., et al. 2024, A&A, 687, A58

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

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

Fujita, Y., Takizawa, M., Yamazaki, R., Akamatsu, H., & Ohno, H. 2015, ApJ, 815, 116

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

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

Gianfagna, G., De Petris, M., Yepes, G., et al. 2021, MNRAS, 502, 5115

Greiner, M., & Enßlin, T. A. 2015, A&A, 574, A86

Heinrich, A., Zhuravleva, I., Zhang, C., et al. 2024, MNRAS, 528, 7274

Hitomi Collaboration (Aharonian, F., et al.) 2016, Nature, 535, 117

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

Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90

Ishiyama, T., Prada, F., Klypin, A. A., et al. 2021, MNRAS, 506, 4210

Itoh, N., Kohyama, Y., & Nozawa, S. 1998, ApJ, 502, 7

Jauzac, M., Eckert, D., Schaller, M., et al. 2018, MNRAS, 481, 2901

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python http://www.scipy.org/

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

Kim, J., Sayers, J., Sereno, M., et al. 2024, A&A, 686, A97

Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301

Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353

Limousin, M., Morandi, A., Sereno, M., et al. 2013, Space Sci. Rev., 177, 155

Ma, C.-J., Ebeling, H., & Barrett, E. 2009, ApJ, 693, L56
Mazzotta, P., & Giacintucci, S. 2008, ApJ, 675, L9

Melin, J. B., & Pratt, G. W. 2023, A&A, 678, A197

Mohapatra, R., Federrath, C., & Sharma, P. 2020, MNRAS, 493, 5838

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

Mohapatra, R., Jetti, M., Sharma, P., & Federrath, C. 2022, MNRAS, 510, 3778

Mroczkowski, T., Dicker, S., Sayers, J., et al. 2012, ApJ, 761, 47

Mroczkowski, T., Nagai, D., Basu, K., et al. 2019, Space Sci. Rev., 215, 17

Muñoz-Echeverría, M., Macías-Pérez, J. F., Pratt, G. W., et al. 2024, A&A, 682, A147

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

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

Pérez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21

Pinzke, A., Oh, S. P., & Pfrommer, C. 2017, MNRAS, 465, 4800

Planck Collaboration XX. 2014, A&A, 571, A20

Planck Collaboration XXIV. 2016, A&A, 594, A24

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

Pointecouteau, E., Giard, M., & Barret, D. 1998, A&A, 336, 44

Ponthieu, N., Grain, J., & Lagache, G. 2011, A&A, 535, A90

Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25

Rajpurohit, K., Brunetti, G., Bonafede, A., et al. 2021a, A&A, 646, A135

Rajpurohit, K., Wittor, D., van Weeren, R. J., et al. 2021b, A&A, 646, A56

Romero, C. E. 2024, ApJ, 975, 197

Romero, C. E., Mason, B. S., Sayers, J., et al. 2017, ApJ, 838, 86

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

Romero, C. E., Gaspari, M., Schellenberger, G., et al. 2024, ApJ, 970, 73

Ruppin, F., Mayet, F., Pratt, G. W., et al. 2018, A&A, 615, A112

Sarazin, C. L. 1986, Reviews of Modern Physics, 58, 1

Sayers, J., Mroczkowski, T., Zemcov, M., et al. 2013, ApJ, 778, 52

Schuecker, P., Finoguenov, A., Miniati, F., Böhringer, H., & Briel, U. G. 2004, A&A, 426, 387

Shaw, L. D., Nagai, D., Bhattacharya, S., & Lau, E. T. 2010, ApJ, 725, 1452

Simionescu, A., Nakashima, S., Yamaguchi, H., et al. 2019, Space Sci. Rev., 215, 24

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

Smith, G. P., Mazzotta, P., Okabe, N., et al. 2016, MNRAS, 456, L74

Sunyaev, R. A., & Zeldovich, Y. B. 1970, Ap&SS, 7, 3

Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments on Astrophysics and Space Physics, 4, 173

Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22

van Weeren, R. J., Röttgering, H. J. A., Brüggen, M., & Hoeft, M. 2010, Science, 330, 347

van Weeren, R. J., de Gasperin, F., Akamatsu, H., et al. 2019, Space Sci. Rev., 215, 16

Vazza, F., Wittor, D., Brüggen, M., & Gheller, C. 2016, Galaxies, 4, 60

XRISM Science Team 2020, ArXiv e-prints [arXiv:2003.04962]

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

Zhuravleva, I., Allen, S. W., Mantz, A., & Werner, N. 2018, ApJ, 865, 53

Zhuravleva, I., Chen, M. C., Churazov, E., et al. 2023, MNRAS, 520, 5157
<figure><img src="image_12.png" /><figcaption>**Fig. A.1.** Example of a CIB realization at 150 GHz obtained using the PYSIDES code (Béthermin et al. 2022).</figcaption></figure>

# Appendix A: Contribution from diffuse backgrounds and foregrounds

## A.1. Cosmic infrared background (CIB)

The CIB is one of the major contaminants in SZ fluctuation analyses, at least in the case of the scales probed with NIKA data. We use the method developed in Béthermin et al. (2017) to construct the mock realizations used in this paper. We apply the PYSIDES code (Béthermin et al. 2022) to the UCHUU simulations (Ishiyama et al. 2021). We obtain a source catalogue with sky positions, redshift and the NIKA 150 GHz fluxes over a 2 deg² region. We then split the catalogue into 1000 independent regions of 10 × 10 arcmin². From each of them, a CIB map was generated by summing the contribution from the sources in the catalogue, which we modelled given their flux, relative position with respect to the centre, and the NIKA 150 GHz Gaussian beam. Catalogue entries for which the flux is above the NIKA detection threshold (taken as 0.4 mJy, corresponding to 2σ, Adam et al. 2017b) were removed to avoid including sources that would be subtracted or masked in real data. In Figure A.1, we show an example of one of the CIB realizations. These Monte Carlo realizations were used to model and account for the CIB along the analysis.

## A.2. Undetected low mass halos

We also considered the contribution from undetected low mass galaxy clusters in the field. We estimated the expected number of clusters in a 30 × 30 arcmin² field by integrating the Tinker et al. (2008) mass function across bins of redshift ($0 < z < 3$) and mass ($2 \times 10^{13} M_{\odot} < M_{500} < 3 \times 10^{15} M_{\odot}$). We then drew the number of clusters from a Poisson distribution based on this estimate, and assigned each cluster a mass and redshift drawn from the mass function, along with a random position in the field. We repeated this process to obtain a 1000 mock cluster catalogues. Assuming spherical symmetry, we modelled the cluster's electronic pressure profile with a gNFW profile (Nagai et al. 2007), which was then integrated along the line of sight to obtain the Compton parameter profile of the cluster. By summing each cluster's contribution we obtained simulated maps of the tSZ

signal. The modelling was undertaken using the minot Python package (Adam et al. 2020). The maps were convolved with the NIKA instrument response function (beam and transfer function) and used to compute its associated power spectrum as described in Section 4.2.3. We found that given the present analysis framework and the scales that we probe, the contribution from undetected low mass clusters is about two orders of magnitude lower than that of the CIB, and therefore fully negligible.

# Appendix B: Uncertainties and systematics in the recovered power spectrum due to the methodology

In this Appendix, we quantify the systematic uncertainties in the power spectrum modelling discussed in Section 4. We first compared the SZ fluctuation power spectrum measured after performing Monte Carlo simulations to that predicted using the framework discussed in Section 4.1.4, in order to test the projection assumptions (Equation 13). We then compared the power spectrum obtained from Monte Carlo realizations, including the instrument response function and data weighting, to that predicted from the modelling methodology discussed in Section 4 in order to test the overall uncertainties associated with the modelling of the pressure fluctuation power spectrum.

Figure B.1 presents the comparison between the Monte Carlo and the model for projection only (top) and accounting for instrumental and weighting effects (bottom). The comparison is performed assuming a morphologically disturbed Arnaud et al. (2010) pressure profile, but we also show the impact of changing this model to a cool-core profile. The fluctuations were set to $\sigma_P = 0.5$, $\alpha = -11/3$, $L_{\text{dis}} = 1$ kpc, and $L_{\text{inj}}$ was varied to change the power spectrum shape since this is key for the comparison. As we can see, the model predictions are accurate within ~10% at scales above or near the peak, $k \ge 1/L_{\text{inj}}$. At large scales, we observe larger deviations for $k \ll 1/L_{\text{inj}}$, as expected (see Section 4.1.4). These deviations are larger for more peaked profiles, such as the cool-core profile in the test shown here. However, they appear in a regime where the power spectrum drops significantly with respect to the peak value and where it cannot be measured in practice.

We conclude that systematic uncertainties in the power spectrum modelling are negligible compared to other uncertainties in the present work. While we only show the case of a reference spectrum and profile, we have checked that these results are robust against the choice of the model and analysis parameters. We also note that the methodology developed here significantly improves over currently used derivations of the power spectrum, which essentially rely on the delta variance framework from Arévalo et al. (2012). For more details, we point the reader to the extensive discussions in Romero (2024), who show that power spectrum modelling uncertainties can reach a factor of a few, depending on the input power spectrum.

# Appendix C: Systematic uncertainties associated with the radial model

The systematic effects associated with the radial model are twofold. Firstly, the chosen radial model may not adequately represent the cluster, such that differences will induce spurious SZ fluctuations. Secondly, fitting the radial model to the data is likely to absorb part of the fluctuations. In this Appendix, we attempt to estimate the magnitude of these systematic effects.
<figure><img src="image_13.png" /><figcaption>**Fig. B.1.** Comparison between the expected and the recovered (mean of Monte Carlo realizations) pressure fluctuation power spectra. The bias, shown at the bottom of each panel, is defined as the ratio between the two and we show the ±10% interval as a dashed line. Error bars correspond to the uncertainties on the mean associated with the variability of Monte Carlo realizations. **Top:** comparison after projection only, i.e. we are testing the validity of Equation 13. The bias is larger than a factor of 1.3 only for the lowest bin for the blue curve, at scales well above the peak. **Bottom:** comparison when including projection effects, the instrument response function and masking, i.e. we are testing the framework discussed in Section 4. The labels A10MD and A10CC refer to the morphologically disturbed and cool-core pressure profile models of Arnaud et al. (2010), respectively.</figcaption></figure>

Fluctuation removal from the radial model fit The first test consists of simulating SZ images, fitting the radial component according to the different models used in this paper, and comparing the resulting fluctuation power spectra to the one expected. Since any change in the radial model will also affect the conversion from the SZ 2D to the pressure 3D power spectrum, the SZ fluctuation power spectrum has to be normalized by $C_{3D→2D}^{(eff)}$. Instrumental and masking effects are neglected because they affect all the simulated data in the same way. The induced bias is

<figure><img src="image_14.png" /><figcaption>**Fig. C.1.** Estimate of the bias induced by the radial model as a function of wavenumber. **Top:** bias induced by the removal of fluctuations implied by the fit of the radial model (Equation C.1). **Bottom:** excess fluctuations induced by the mismodelling, relative to a reference pressure fluctuation power spectrum ($\sigma_p = 0.5$, $L_{inj} = 750$ kpc, $\alpha = -11/3$ and $L_{dis} = 1$ kpc). The different colours represent the different models, according to Section 5.1. In the right panel, the different line styles refer to the different input radial model. In both panels, the vertical dashed line gives $k = 1/L_{inj}$.</figcaption></figure>

thus computed as

$$ \text{bias} = \frac{\left( \frac{P_{\omega\delta y/\bar{y}}}{C_{3D\to2D}^{(\text{eff})}} \right)_{\text{fitted RM}}}{\left( \frac{P_{\omega\delta y/\bar{y}}}{C_{3D\to2D}^{(\text{eff})}} \right)_{\text{true RM}}} \quad (C.1) $$

In the top panel of Figure C.1, we show the mean bias induced by radial model fitting to the data, as a function of $k$, obtained by averaging 100 mock realizations. The model used to generate the noiseless mocks was set to a spherical, morphologically disturbed (Arnaud et al. 2010) pressure profile ($M_{500} = 2.5 \times 10^{15}$ M$_\odot$), and the pressure fluctuations were defined as lognormal with $\sigma_p = 0.5$, $L_{inj} = 750$ kpc, $\alpha = -11/3$ and $L_{dis} = 1$ kpc. The fitted radial models are those discussed in Section 5.1. We note that in the context of this test, the models should perfectly fit the mocks in the absence of fluctuations. We can see that increasing the model complexity implies an increasing removal of fluctuations. Most of the filtering appears on large scales since the models are relatively smooth. Up to about 50% filtering is observed near the power spectrum peak and up to
80% on larger scales, for the most complex models. The bias converges to ~ 1 on small scales.

**Incorrect radial model-induced bias** Quantification of the impact of an incorrect radial model choice on the results is more difficult, since the true model of MACS J0717.5+3745 is unknown. Moreover, the excess fluctuation bias should be quantified relative to a given fluctuation model. We thus proceeded as follows. We simulated SZ data given a radial model, without fluctuations, fit the image with another (different) radial model, and computed the SZ fluctuation power spectrum of the residuals. We then compare the induced spurious surface brightness fluctuations to a reference pressure fluctuation model given by $\sigma_{\rho} = 0.5$, $L_{\text{inj}} = 750$ kpc, $\alpha = -11/3$ and $L_{\text{dis}} = 1$ kpc. The input radial models were defined in terms of increasing complexity: 1) a centered, spherical, morphologically disturbed Arnaud et al. (2010) profile; 2) as input model 1, but miscenteried by 20 arcsec along R.A. and 20 arcsec along Dec. (this compares well with the difference between the X-ray and SZ peak for MACS J0717.5+3745); 3) as input model 2, but with an elliptical profile with $q_{\text{int,maj}} = 0.7$; and 4) as input model 3, but changing the slopes of the pressure profile with $a \to a + 0.1$, $b \to b - 0.7$, $c \to c - 0.2$, which typically corresponds to the difference one may have between different clusters (e.g. Arnaud et al. 2010; Planck Collaboration Int. V 2013).

In Figure C.1, bottom panel, we can observe that the excess power generally increases with input model complexity, and decreases with fit model complexity. This is expected since more complex models are better able to account for a mismatch with the data. The bias is also more prominent on large scales, at least given the tested models. As expected, no bias is observed when the model is able to fit the data perfectly. For instance, this is the case for the first input model, regardless of the fit model.

According to this test, RM₁ may lead to an excess in the power spectrum larger than 10 times the reference power spectrum on large scales, because of miscentering. RM₂ may lead to a bias which is of the same order of magnitude as the reference spectrum, essentially because it does not account for ellipticity. RM₃ may lead to up to 40% bias due to the mismatch of the shape of the pressure profile (given the different shapes for the input and the fit used here). RM₄ is able to mitigate the mismatch in the shape by fitting $P_0$ and $r_\rho$ so that the bias reaches 1% at most. Finally, RM₅ is able to perfectly match all the input pressure profiles so that it does not lead to any bias. Although these tests give an insight into the level of bias induced by the choice of the radial model, we stress that they strongly depend on the choice of the input models, and should only be taken as a qualitative estimate.

Nevertheless, these two tests indicate that systematic effects in the choice of the radial model are likely to alter the measured power spectrum on scales comparable to or larger to the peak, and may thus alter the recovery of the peak power spectrum (e.g. $L_{\text{inj}}$, $A(k_{\text{peak}})$). Given the results obtained here, RM₃ appears as a reasonable baseline choice since the induced excess power expected compares well to the filtering, so that in the end the associated systematic effects may be minimal.

## Appendix D: Systematic uncertainties associated with the kinetic Sunyaev-Zel'dovich signal

MACS J0717.5+3745 is known to host significant kSZ signal. Following the recommendation from the public release of NIKA SZ data (see Section 7), we used the kSZ corrected map to

test the stability of the results. Figure D.1 displays the residuals between the data and the considered models, as in Figure 6, but accounting for the kSZ correction prior the fit. Similarly, Figure D2 shows the 3D pressure fluctuation power spectra accounting for the kSZ correction. The residual and recovered spectra are comparable to those of the baseline results, although model RM₅ is not able to fit the map as well as in the earlier version, and its resulting power spectra aligns better with RM₂, RM₃ and RM₄. The results of the power spectrum constraints are listed in Table 3.

## Appendix E: Impact of the choice of the pressure fluctuation power spectrum slope

The results presented in this paper have been obtained using a fixed value of the slope of the pressure fluctuation power spectrum, $\alpha = -11/3$. This choice corresponds to the canonical value for a Kolmogorov cascade. However, the true power spectrum slope may slightly differ, and the power spectra associated with merging clusters are often rather steep. While the slope of the power spectrum cannot be constrained from the data, due to the combination of beam smoothing and noise that dominates on small scales, this appendix quantifies the dependence on the results of our baseline choice.

We use the reference analysis framework together with model RM₃ and reproduce the results assuming $A_{\delta P/P}(k_{3D}) \propto k^{-1/5}$ (shallow) and $A_{\delta P/P}(k_{3D}) \propto k^{-4/5}$ (steep), which bracket reasonable values of the slope (Gaspari et al. 2014). This corresponds to $\alpha = -3.4$ (shallow) and $\alpha = -4.6$ (steep), respectively, to be compared with $\alpha = -11/3 \approx -3.7$ in the canonical case. As we can see in Figure E1, the data are not sufficiently sensitive to clearly distinguish between these different models, although we note that a slightly better residual is obtained for the steep spectrum.

We report the changes of the main analysis parameters in Table E.1 when varying the slope of the power spectrum. Because the model parameters are degenerate, using the shallow (steep) spectrum implies an increase (decrease) of the normalization $\sigma_{\rho}$ of about 25%. The injection scale is more stable and remains stable within 10%. The derived Mach number, velocity dispersion, kinetic to thermal pressure and hydrostatic mass bias all decrease (increase) when increasing (decreasing) the steepness of the power spectrum. The changes are smaller than for the amplitude of the power spectrum due to the anticorrelation between $\alpha$ and $\sigma_{\rho}$ and because we are using the integrated spectrum as a proxy for the Mach number. While the results of the analysis are affected by the choice of the slope of the pressure fluctuation power spectrum, fixing this parameter within a reasonable range does not alter our conclusions.

<table><caption>**Table E.1.** Change in the values of the key parameters derived from the reference analysis when varying the slope of the pressure fluctuations power spectrum.</caption><thead><tr><th>Slope (α)<br>(-)</th><th>σ<sub>ρ</sub><br>(-)</th><th>L<sub>inj</sub><br>(kpc)</th><th>M<sub>3D</sub><br>(km/s)</th><th>σ<sub>v</sub><br>(-)</th><th>P<sub>kin</sub>/P<sub>kin+P<sub>th</sub></sub><br>(-)</th><th>b<sub>HSE,Δ=2500</sub></th></tr></thead><tbody><tr><td>Shallow (-3.4)</td><td>0.74</td><td>880</td><td>0.68</td><td>1350</td><td>0.20</td><td>0.31</td></tr><tr><td>Reference (≈ -3.7)</td><td>0.58</td><td>873</td><td>0.59</td><td>1176</td><td>0.16</td><td>0.29</td></tr><tr><td>Steep (-4.6)</td><td>0.45</td><td>790</td><td>0.47</td><td>925</td><td>0.11</td><td>0.27</td></tr></tbody></table>

Notes. These results are obtained using the radial model RM₃.
<figure><img src="image_15.png" /><figcaption>Fig. D.1. Same as Figure 6 but with kSZ correction included.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. D2. Same as Figure 7 right panel but with kSZ correction included.</figcaption></figure>

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

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

<figure><img src="image_19.png" /><figcaption>Fig. E1. Same as Figure 7 left panel but when fixing the slope of the power spectrum of pressure fluctuations to $\alpha = -3.4$ (top), $\alpha \approx -3.7$ (middle), and $\alpha = -4.6$ (bottom).</figcaption></figure>