# Investigating the turbulent hot gas in X-COP galaxy clusters

Simon Dupourqué, Nicolas Clerc, Etienne Pointecouteau, Dominique Eckert, Stefano Ettori, Franco Vazza

► To cite this version:

Simon Dupourqué, Nicolas Clerc, Etienne Pointecouteau, Dominique Eckert, Stefano Ettori, et al. Investigating the turbulent hot gas in X-COP galaxy clusters. Astronomy and Astrophysics - A&A, 2023, 673, pp.A91. 10.1051/0004-6361/202245779. hal-04076308

HAL Id: hal-04076308

https://hal.science/hal-04076308

Submitted on 24 Aug 2023

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
# Investigating the turbulent hot gas in X-COP galaxy clusters

S. Dupourqué<fnref n="1" />, N. Clerc<fnref n="1" />, E. Pointecouteau<fnref n="1" />, D. Eckert<fnref n="2" />, S. Ettori<fnref n="3" />, and F. Vazza<fnref n="4" />,<fnref n="5" />,<fnref n="6" />

Received 23 December 2022 / Accepted 17 March 2023

## ABSTRACT

**Context.** Turbulent processes at work in the intracluster medium perturb this environments, impacting its properties, displacing gas, and creating local density fluctuations that can be quantified via X-ray surface brightness fluctuation analyses. Improved knowledge of these phenomena would allow for a more accurate determination of the mass of galaxy clusters, as well as a better understanding of their dynamic assembly.

**Aims.** In this work, we aim to set constraints on the structure of turbulence using X-ray surface brightness fluctuations. We seek to consider the stochastic nature of this observable and to constrain the structure of the underlying power spectrum.

**Methods.** We propose a new Bayesian approach, relying on simulation-based inference to account for the whole error budget. We used the X-COP cluster sample to individually constrain the power spectrum in four regions and within $R_{500}$. We spread the analysis on the entire set of 12 systems to alleviate the sample variance. We then interpreted the density fluctuations as the result of either gas clumping or turbulence.

**Results.** For each cluster considered individually, the normalisation of density fluctuations correlate positively with the Zernike moment and centroid shift, but negatively with the concentration and the Gini coefficient. The spectral index within $R_{500}$ and evaluated over all clusters is consistent with a Kolmogorov cascade. The normalisation of density fluctuations, when interpreted in terms of clumping, is consistent within $0.5R_{500}$ with the literature results and numerical simulations; however, it is higher between 0.5 and $1R_{500}$. Conversely, when interpreted on the basis of turbulence, we deduce a non-thermal pressure profile that is lower than the predictions of the simulations within $0.5 R_{500}$, but still in agreement in the outer regions. We explain these results by the presence of central structural residues that are remnants of the dynamical assembly of the clusters.

**Key words.** X-rays: galaxies: clusters – galaxies: clusters: intracluster medium – turbulence

## 1. Introduction

The intracluster medium (ICM) is the primary baryonic component of galaxy clusters. The hot gas ($T \sim 10^7-10^8$ K) is governed by many physical processes that introduce perturbations on various scales. In the inner parts of galaxy clusters, the strong cooling and accreting gas trigger feedback from the active galactic nucleus (AGN) hosted by the central galaxy (McNamara & Nulsen 2012; Voit et al. 2017). In the outer parts, merger events and accretion from the cosmic web generate shocks, adiabatic compression, and turbulent cascades (Nelson et al. 2012). Throughout the history of the dynamic assembly of galaxy clusters, these perturbations will play a predominant role in the non-thermal heating of the ICM (Bennett & Sijacki 2022), generating non-thermal pressure, which drives the hydrostatic mass bias when it is not accounted for in the mass estimation of massive halos when assuming their hydrostatic equilibrium (Piffaretti & Valdarnini 2008; Lau et al. 2009; Meneghetti et al. 2010; Nelson et al. 2014; Biffi et al. 2016; Pratt et al. 2019).

Turbulence is a stochastic process that occurs in high Reynolds number flows. It can be interpreted qualitatively as the transport of kinetic energy from a large injection scale to a

small viscous scale, which then dissipates as heat in the medium. Putting constraints on the turbulence that occurs in the intracluster medium is of great interest when it comes to characterising non-thermal heating, since turbulent motions can be responsible for more than 80% of these effects – and potentially all of it (Vazza et al. 2018; Angelinelli et al. 2020). The level of baseline turbulence implied by most of the modern cosmological simulations also appears to be in line with the implications of the latest discoveries of very extended radio emission, namely, out to at least $R_{500}$ in clusters of galaxies classified as ‘mega-halos’ (Cuciti et al. 2022). The detection of such unprecedentedly large volumes of cluster environments filled with relativistic particles and magnetic fields is currently best explained by the ubiquitous re-acceleration of electrons by turbulence, which is expected to be present at a similar level in all clusters on these radii, at a level compatible with simulations (see also Botteon et al. 2022).

Direct observations of this phenomenon is possible via studies of the centroid shifts and broadening of the ICM emission lines to derive the bulk and turbulent motions, respectively. Upper limits on the fraction of energy in the form of turbulent motions have been derived using the XMM-Newton spectrometer for samples of clusters (Sanders et al. 2011; Pinto et al. 2015). A novel approach using the EPIC-pn detector has

<footnote n="1">¹ IRAP, Université de Toulouse, CNRS, CNES, UT3-PS, Av. du Colonel Roche 9, 31400 Toulouse, France
e-mail: sdupourque@irap.omp.eu</footnote>

<footnote n="2">² Department of Astronomy, University of Geneva, Ch. d'Ecogia 16, 1290 Versoix, Switzerland</footnote>

<footnote n="3">³ INAF, Osservatorio di Astrofisica e Scienza dello Spazio, via Piero Gobetti 93/3, 40129 Bologna, Italy</footnote>

<footnote n="4">⁴ Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, 40122 Bologna, Italy</footnote>

<footnote n="5">⁵ Hamburger Sternwarte, University of Hamburg, Gojenbergsweg 112, 21029 Hamburg, Germany</footnote>

<footnote n="6">⁶ Istituto di Radioastronomia, INAF, Via Gobetti 101, 40122 Bologna, Italy</footnote>
allowed for the direct measurement of bulk flows in the Coma and Perseus clusters (Sanders et al. 2020), showing consistent results with Hitomi and was also applied to the Virgo and Centaurus clusters (Gatuzz et al. 2022a,b). How these measurements translate into a correction for the hydrostatic bias is illustrated in, for instance, Ota et al. (2018), Ettori & Eckert (2022). The first credible spatially resolved measurements of this type in the X-rays were achieved with the Hitomi observatory (The Hitomi Collaboration 2016) at the centre of the Perseus cluster. Future direct measurements will be obtained with the coming of the new generation of X-ray integral fields unit, that is XRISM/Resolve in the coming years (XRISM Science Team 2020) and Athena/X-IFU in the long term (Barret et al. 2020).

It is also possible to characterise these phenomena using indirect observables. The X-ray surface brightness fluctuations are mostly due to density fluctuations in the ICM, and have been studied for the first time in the Coma cluster by Schuecker et al. (2004), and later by Churazov et al. (2012), deriving constraints on the density fluctuation power spectrum on scales ranging between 30 kpc and 500 kpc. Various theoretical works have suggested a strong link between turbulent velocities and density fluctuations (e.g., Zhuravleva et al. 2014; Gaspari et al. 2014; Mohapatra et al. 2020, 2021; Simonte et al. 2022), indicating that the study of brightness fluctuations could constrain the turbulent processes that occur in the ICM. This methodology was applied to Perseus (Zhuravleva et al. 2015), and in the cool cores of ten galaxy clusters (Zhuravleva et al. 2018). These studies constrained turbulent velocities to $\le 150$ km s$^{-1}$ at scales smaller than 50 kpc, which is consistent with direct measurements from Hitomi.

Turbulent processes, and thereby the resulting density fluctuations, originate from chaotic processes that can be assimilated to random fields observed in spatially finite regions. The stochastic nature of this observable coupled with the finite size of the observations means it intrinsically carries an additional variance, which is later referred to as the sample variance. This effect has been studied for the structure function of turbulent velocities for XRISM (ZuHone et al. 2016) and Athena X-IFU mock observations (Clerc et al. 2019; Cucchetti et al. 2019), and it has been dominant in the error budget at spatial scales $\ge 50$ kpc. It is expected to play a significant role in fluctuation analyses, leading to an underestimation of the error budget when it is not accounted for. In this work, we propose a novel approach based on the forward modelling of observables related to surface brightness fluctuations, which allows, for the first time, the full error budget associated with their stochastic nature to be considered. In this work, we apply this methodology to the XMM-Newton Cluster Outskirts Project (X-COP, Eckert et al. 2017) sample and derive individual constraints on the statistical properties and spatial distribution of the density fluctuations via the characterisation of the 3D power spectrum of the associated random field. We combine these constraints under the assumption of comparable dynamics within the clusters to obtain stronger constraints on the whole sample, which can be related to the gas clumping and the turbulent processes occurring in the ICM.

In Sect. 2, we describe our methodology as well as the validation process. In Sect. 3, we present the individual results and the joint constraints on density fluctuations. In Sect. 4, we interpret these constraints as being gas clumping or coming from turbulent processes and contributing to non-thermal heating. We also inspect the correlations between the obtained parameters and the dynamic state of the clusters and discuss the various limitations of our approach. Throughout this paper, we assume a flat $\Lambda$CDM cosmology with $H_0 = 70$ km s$^{-1}$ and $\Omega_m = 1 - \Omega_\Lambda = 0.3$.

Scale radii are defined according to the critical density of the Universe at the corresponding redshift. The Fourier transform conventions are highlighted in Appendix B.

## 2. Data and method

### 2.1. X-COP Sample

The XMM-Newton Cluster Outskirts Project (X-COP, Eckert et al. 2017) is an XMM Very Large Program dedicated to the study of the X-ray emission of cluster outskirts. This sample comprises 12 massive (~3–10 × $10^{14}$ M$_\odot$) and local (0.04 < z < 0.1) galaxy clusters, for a total exposure time of ~2 Ms. The targets were selected from the first Planck catalogue (Planck Collaboration XXIX 2014) as (i) resolved sources regarding Planck's spatial resolution (i.e. all clusters have $R_{500} > 10'$), (ii) with high signal-to-noise ratio ($S/N > 12$), and (iii) in directions that ensure low hydrogen column densities to prevent soft-X ray absorption ($N_H < 10^{21}$ cm$^{-2}$). The thermodynamical properties of the X-COP sample have been characterised (Ghirardini et al. 2019). The sample was further used to derive the first observational constraints on the amount of non-thermal pressure in the outskirts of clusters (Eckert et al. 2019).

### 2.2. Data preparation

This project is based on the publicly available X-COP data products (Ghirardini et al. 2019; Ettori et al. 2019)$^1$. Here, we recall the main steps of the data reduction, image extraction, and point source subtraction. We refer to Ghirardini et al. (2019), Eckert et al. (2019) for further details on the data preparation procedure.

First, the data were reduced using the XMMSAS v13.5 software package and the Extended Source Analysis Software (ESAS) procedure; namely, the raw data were screened using the `emchain` and `epchain` executables to extract raw event files. We then extracted the light curves of the entire field of view to exclude time periods affected by soft proton flares.

From the cleaned event files, we extracted count maps from each observation in the [0.7–1.2] keV from the three detectors of the EPIC instrument (MOS1, MOS2, and pn). We then computed the corresponding exposure maps including the telescope vignetting using the `eexpmap` task. Finally, we created non X-ray background maps from the filter-wheel-closed event files, appropriately rescaling the filter-wheel-closed data to match the count rates observed in the unexposed corners of the individual CCDs. The contribution of residual soft protons was estimated by measuring the ratio between the high-energy count rates inside and outside field of view (IN-OUT) and its relation to the residual soft proton component, calibrated over ~500 blank-sky pointings (Ghirardini et al. 2018).

For each cluster, the individual images, exposure maps, and background maps were stacked and concatenated into joint EPIC mosaic images. Point sources were extracted by joining the [0.5–2] keV and [2–7] keV energy bands using the XMMSAS tool `ewavelet` with a $S/N = 5$ threshold and the detected point sources were masked when their count rate exceeded a threshold set by the peak of the logN-logS distribution, such that the source exclusion threshold is homogeneous over the entire field.

<footnote>$^1$ https://dominiqueeckert.wixsite.com/xcop</footnote>
## 2.3. Modelling surface brightness

To characterise the surface brightness fluctuations of X-COP clusters, we choose first to subtract the bulk of the main cluster emission. We determined an average emissivity profile, assuming radial symmetry, and biaxiality in directions orthogonal to the line of sight. This was achieved by fitting an average surface brightness profile on the X-ray images in the [0.7–1.2] keV energy band. In this section, we note our use of the $\mathbf{r} = (x, y, \ell)$ in the 3D parametrisation and $\rho = (x, y)$ in the 2D parametrisation, with $\ell$ as the coordinate along the line of sight and $(x, y)$ as the coordinates mapping the dimensions on the image. The true surface brightness $\Sigma_X$ is given by:

$$ \Sigma_X(\rho) = \frac{1}{4\pi(1+z)^4} \int_{-\infty}^{+\infty} \Lambda(T, Z) n_c^2(\mathbf{r}) d\ell, \quad (1) $$

where $\Lambda$ is the intra-cluster gas cooling function in the [0.7–1.2] keV band, $T$ is the 3D-temperature, $Z$ is the metallicity, and $n_c$ is the 3D electron density. The surface brightness we observe is convolved with the XMM-Newton response functions, and partially absorbed by the Galactic hydrogen column density. We then define the observed surface brightness $S_X$ as:

$$ S_X(\rho) = \int_{-\infty}^{+\infty} \Psi(\mathbf{r}) n_c^2(\mathbf{r}) d\ell + B, \quad (2) $$

where $\Psi(\mathbf{r})$ encompasses the cooling function, the cosmological dimming, the Galactic absorption, and the convolution with XMM-Newton response functions, while $B$ is a constant surface brightness background left as a free parameter. The 3D density profile is modelled using a modified version of the Vikhlinin functional form (Vikhlinin et al. 2006), defined as follows:

$$ n_c^2(\mathbf{r}) = n_{c,0}^2 \left(\frac{r}{r_c}\right)^{-\alpha} \frac{(1+r^2/r_c^2)^{-3\beta+\alpha/2}}{(1+r^\gamma/r_s^\gamma)^{\epsilon/\gamma}}, \quad (3) $$

where $n_{c,0}$ is the central density, $r_c$ and $r_s$ are two scale radii, and $\alpha, \beta, \gamma, \epsilon$ are parameters that define the slopes and smoothness of transition between the power laws in this equation. We set $\alpha = 0$, as done by Shi et al. (2016), to remove the central singularity, which causes problems due to low values of $r$ on the image. We fix $\gamma = 3$ to prevent degeneracies with the $\epsilon$ parameter, as the data do not allow for a simultaneous determination of $\epsilon$ and $\gamma$.

We computed the cooling rate $\Lambda$ and the Galactic absorption using a functional approximation fitted to the count rate derived from a “PhAbs*APEC" model under XSPEC 12.11.1, as detailed in Appendix C. Galactic absorption admits non-negligible fluctuations (~10%) at the scale of the large X-COP mosaics. To take them into account, we used the $N_H$ data from the HI4PI survey (HI4PI Collaboration 2016), which provide a mapping of the Galactic column density at the spatial scale of 12–16 arcmin. We performed the projection along the line of sight with a double exponential quadrature (Takahasi & Mori 1973; Mori & Sugihara 2001). The centre position on the model image is left free and fitted as a Gaussian perturbation of the cluster centre, as provided in the Planck catalogue (Planck Collaboration XXIX 2014). The effective centre is then parametrised as:

$$ \begin{pmatrix} x_c \\ y_c \end{pmatrix} = \begin{pmatrix} N(1, 0.5) & 0 \\ 0 & N(1, 0.5) \end{pmatrix} \begin{pmatrix} x_{\text{Planck}} \\ y_{\text{Planck}} \end{pmatrix}. \quad (4) $$

To take non-sphericity into account, we fit for an elliptical brightness distribution rather than a spherical one. The positions

($x, y$) are distorted in a surface-conserving way into a new configuration $(\tilde{x}, \tilde{y})$, introducing two additional free parameters to describe the ellipse angle $\theta$ and eccentricity $e$:

$$ \begin{pmatrix} \tilde{x} \\ \tilde{y} \end{pmatrix} = \begin{pmatrix} \frac{1}{\sqrt{1-e^2}} & 0 \\ 0 & \sqrt{1-e^2} \end{pmatrix} \begin{pmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{pmatrix} \begin{pmatrix} x - x_c \\ y - y_c \end{pmatrix}. \quad (5) $$

This is formally equivalent to assuming a biaxiality of the 3D distribution of the cluster density over the image directions. The distributions of the free parameter posteriors are obtained using Bayesian inference. We define a Poisson likelihood on images rebinned within $R_{500}$ with a Voronoi tessellation for a target count of 100 counts per bin (Cappellari & Copin 2003). As an example, for A3266, the median characteristic bin size and its 14th–86th percentile is $10^{-4}$ arcsec. The model count for each bin is approximated by the product of the surface brightness $S_X$ estimated in the centre times the area and exposure of the bin. This is a reasonable approximation, since the size of the bins is small over the regions where the surface brightness varies a lot (see Appendix D for further discussions). The total number of photon counts is obtained by adding the estimated background to the model counts, as defined in Eq. (D.1). We mask the regions associated with point sources and also remove identified substructures and groups such as in A1644, A2142, A644, A85, and RXC1825. The posterior distributions are sampled with the No U-Turn Sampler (Hoffman & Gelman 2014), as implemented in the numpyro library (Bingham et al. 2019; Phan et al. 2019). The prior distributions are displayed in Appendix E.

## 2.4. Fluctuations and power spectrum

We want to quantify the surface brightness fluctuations under the assumption that they come exclusively from intrinsic density fluctuations. In this framework, the density profile can be decomposed into a rest component $n_c^2$ and relative fluctuations $\delta$ whose amplitude is smaller than unity:

$$ S_X(\rho) = \int_{-\infty}^{+\infty} \Psi(\mathbf{r}) n_c^2(\mathbf{r}) (1 + \delta(\mathbf{r}))^2 d\ell. \quad (6) $$

By linearising the Eq. (6) in $\delta$, it follows that, at first order, the raw image corresponds to the sum of an unperturbed image $S_{X,0}$ and a surface brightness fluctuation map:

$$ S_X(\rho) \approx S_{X,0}(\rho) + 2 \int_{-\infty}^{+\infty} \epsilon_0(\mathbf{r}) \delta(\mathbf{r}) d\ell, \quad (7) $$

where $\epsilon_0(\mathbf{r}) = \Psi(\mathbf{r})n_{c,0}^2(\mathbf{r})$ is the unperturbed emissivity of the cluster. The fluctuations can be defined in a theoretically equivalent way by taking the difference or ratio of the perturbed to the unperturbed surface brightness. However, from an observational perspective, both approaches have their advantages and disadvantages, which we quantitatively discuss in Sect. 2.6. We define the absolute and relative fluctuations, $\Delta_{\text{abs/rel}}$, as follows:

$$ \Delta_{\text{abs}}(\rho) \stackrel{\text{def}}{=} \frac{S_X(\rho) - S_{X,0}(\rho)}{2} \approx \int_{-\infty}^{+\infty} \epsilon_0(\mathbf{r}) \delta(\mathbf{r}) d\ell, \quad (8) $$

$$ \Delta_{\text{rel}}(\rho) \stackrel{\text{def}}{=} \frac{S_X(\rho)}{S_{X,0}(\rho)} \approx 1 + 2 \frac{\int_{-\infty}^{+\infty} \epsilon_0(\mathbf{r}) \delta(\mathbf{r}) d\ell}{\int_{-\infty}^{+\infty} \epsilon_0(\mathbf{r}) d\ell}. \quad (9) $$

In practice, we use the definition of Zhang et al. (2023) as a variant of Eq. (9) for relative fluctuation maps, which performs better in low statistic regions while having a similar qualitative behaviour. The difference (as presented in e.g.,
<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparison between the absolute fluctuations $\Delta_{\text{absolute}}$ from Eq. (8) (left) and the relative fluctuations $\Delta_{\text{relative}}$ from Eq. (9) (right) using the surface brightness data and best fit model image of A2319. For display purposes, the images are filtered by a Gaussian kernel of 7.5". The successive contours represent the distance to the centre in units of $R_{500}$.</figcaption></figure>

Khatri & Gaspari 2016) generates maps of fluctuations whose amplitude decreases with distance from the centre. This method defines fluctuations everywhere on the image, but biases their amplitude by over-representing fluctuations in the brightest areas of the image. The ratio (e.g., Churazov et al. 2012) instead defines well-scaled fluctuations, at the expense of intensifying the noise in the low-brightness regions. In Fig. 1, we compare the fluctuation maps obtained using the absolute and relative definition for A2319. It should be noted that both methods suffer from the lack of signal over regions of the order of $R_{500}$, here the photon statistics is insufficient to perform an efficient fluctuation analysis. To study the contribution of each spatial scale to surface brightness fluctuations, it is natural to perform the analysis in Fourier space. We define the power spectrum of the fluctuations, $\mathcal{P}_{2D,\Delta}(k_\rho)$, as follows:

$$ \mathcal{P}_{2D,\Delta}(k_\rho) = \frac{1}{2\pi} \int |\hat{\Delta}(k_\rho)|^2 d\varphi, \qquad (10) $$

with $\hat{\Delta}$ the Fourier transform of the two-dimensional (2D) map $\Delta$ and $\varphi$ the azimuthal angle on the image. The numerical evaluation of $P_{2D,\Delta}$ is performed using the method described by Arévalo et al. (2012), which computes the variance of images filtered by Mexican hats on a characteristic scale to estimate the azimuthal average of the power spectrum. A detailed explanation of the method is provided in Appendix F. Other approaches found in the literature compute the structure function, $\mathcal{SF}_{2D}(\rho)$, of the image instead of the $\mathcal{P}_{2D}(k)$ (e.g., Roncarelli et al. 2018; Clerc et al. 2019; Cucchetti et al. 2019). We emphasise that these two approaches are fully equivalent, since both $\mathcal{SF}_{2D}$ and $\mathcal{P}_{2D}$ are measures of the second order correlation of the fluctuation field and are related by the following bijection, which can be derived using the Wiener–Khinchin theorem:

$$ \mathcal{SF}_{2D}(\rho) = 2 \int (1 - e^{2\pi k_\rho \cdot \rho}) \mathcal{P}_{2D}(k_\rho) dk_\rho. \qquad (11) $$

We stuck to the formulation in power spectrum as it is more practical in our chosen framework, in particular, to generate the density fluctuation random fields (see Sect. 2.5).

## 2.5. Density fluctuations as a Gaussian random field

The next step in the analysis is to properly relate $\mathcal{P}_{2D,\Delta}$ to the 3D power spectrum of density fluctuations. Previous work based on this method uses a proportional relationship between the 2D and 3D spectra, involving the Fourier transform of the assumed emissivity (see Eq. (11) in Churazov et al. 2012). This approximation is valid as long as the surface brightness can be considered constant, but is no longer valid when considering large areas with strong gradients and important variations in the surface brightness. Moreover, this approach does not consider the stochastic nature of this observable and the associated sample variance. To address this, we chose to consider density fluctuations as a random field. Since many works suggest deep connections between density fluctuations and turbulent processes (Zhuravleva et al. 2014; Gaspari et al. 2014; Simonte et al. 2022), we modelled the density fluctuation as a Gaussian random field (GRF) with a 3D power spectrum, $\tilde{\mathcal{P}}_{3D,\delta}$, noted with a bar to distinguish it from the 3D power spectrum measured from a single realisation of the random field. By noting $\langle \cdot \rangle$ the averaging operator over realisations, it is defined as:

$$ \tilde{\mathcal{P}}_{3D,\delta}(k_r)\delta^D(k_r + k'_r) = \langle \tilde{\delta}(k_r)\tilde{\delta}(-k'_r) \rangle = \langle |\tilde{\delta}(k_r)|^2 \rangle. \qquad (12) $$

The link between turbulence and density fluctuations suggests both their power spectra take a similar form. As such, we adopted the simplest possible turbulent model as a Kolmogorov cascade. We chose the following functional form as proposed by ZuHone et al. (2016), which has no direction dependence under the isotropic hypothesis:

$$ \tilde{\mathcal{P}}_{3D,\delta}(k) = \sigma^2 \frac{e^{-(k/k_{\text{dis}})^2} e^{-(k_{\text{inj}}/k)^2} k^{-\alpha}}{\int 4\pi k^2 dk e^{-(k/k_{\text{dis}})^2} e^{-(k_{\text{inj}}/k)^2} k^{-\alpha}}, \qquad (13) $$

where $k_{\text{inj}}$ and $k_{\text{disp}}$ are, respectively, the injection and dissipation scales, $\alpha$ is the inertial range spectral index, and $\sigma_\delta^2$ is the variance of fluctuations. For ease of understanding, we used spatial scales of injection rather than frequency scales, defined as $\ell_{\text{inj,dis}} = 1/k_{\text{inj,dis}}$ and expressed in units of $R_{500}$. Since the dissipation scale (of order $\ge 10^{-3} R_{500}$, Brüggen et al. 2015) is expected to be much lower than the spatial resolution of our images ($\sim 10^{-2} R_{500}$ at the average redshift of the X-COP clusters), we set it to $10^{-3} R_{500}$ for all X-COP clusters. This choice
<table><caption>**Table 1.** Regions used in the analysis, minimum and maximum scale achievable in Fourier space according to Nyquist–Shannon theorem.</caption><thead><tr><th>Region</th><th>Radius</th><th>Min scale</th><th>Max scale</th></tr></thead><tbody><tr><td>(I)</td><td>0 &lt; r &lt; R<sub>500</sub>/10</td><td>0.02R<sub>500</sub></td><td>0.10R<sub>500</sub></td></tr><tr><td>(II)</td><td>R<sub>500</sub>/10 &lt; r &lt; R<sub>500</sub>/4</td><td>0.02R<sub>500</sub></td><td>0.23R<sub>500</sub></td></tr><tr><td>(III)</td><td>R<sub>500</sub>/4 &lt; r &lt; R<sub>500</sub>/2</td><td>0.02R<sub>500</sub></td><td>0.43R<sub>500</sub></td></tr><tr><td>(IV)</td><td>R<sub>500</sub>/2 &lt; r &lt; R<sub>500</sub></td><td>0.02R<sub>500</sub></td><td>0.86R<sub>500</sub></td></tr></tbody></table>

has little effect on the other parameters of the spectrum, since the exact value of the dissipation length will very marginally affect the normalisation of the spectrum and, therefore, $\sigma_\delta$. The GRF hypothesis is valid as long as the fluctuations within the clusters can be represented by an isotropic homogeneous field. However, as clusters are dynamical objects, they are subject to, for instance, merging, accretion, or sloshing events at the centre, which will induce residuals that cannot be described by a GRF. This assumption is discussed further in Sect. 4.6.

## 2.6. Optimal definition and observable

As seen in Sect. 2.3, the definition of fluctuation and aperture size (and, more generally, the shape of the mask we choose) has various implications with regard to the signal we measure. For example, the fluctuation map, $\Delta_{rel}$, obtained through a ratio is expected to over-represent the Poisson noise for radii $\sim R_{500}$, and, conversely, for the maps, $\Delta_{abs}$, obtained with a subtraction, the very low significance of the signal at this distance to the centre does not provide additional information. To quantify the comparison between absolute and relative approaches, we used the numerical simulations (further explained in Sect. 2.8) to study the signal-to-noise ratio (S/N) that can be obtained in each case.

We simulated 100 mock surface brightness maps for density fluctuations with $\sigma_\delta = 0.32$, $\ell_{inj} = 0.05R_{500}$, $\alpha = 11/3$ (as defined in Eq. (13)) and computed the associated $\mathcal{P}_{2D}$ in the regions defined in Table 1. The maximum scale was defined using Nyquist-Shannon criterion, which is half of the highest scale accessible in the mask. For a ring of inner radius, $R_{min}$, and outer radius, $R_{max}$, it is given by $(R_{max}^2 - R_{min}^2)^{1/2}$. To quantify the valuable information in the power spectrum, we compared it to the power spectrum that is obtained by switching off density fluctuations, $N_{2D}$. This quantity represents the fluctuations which are expected when observing a perfect cluster at rest, without any density fluctuations. It can be obtained by forward modelling the power spectrum, but with a zero-normalisation density fluctuation field, accounting for the dispersion due to the mean model fit and the Poisson noise by computing it for models drawn from the posterior distribution and making various realisation of the count image. We define the S/N of the power spectrum as follows for each realisation of the mock density fluctuations:

$$S/N_{\Delta} \stackrel{\text{def}}{=} \frac{\mathcal{P}_{2D,\Delta}}{N_{2D,\Delta}}; \quad \Delta \in \{\text{abs, rel}\}. \qquad (14)$$

This quantity can trace the excess of surface brightness fluctuation when compared to what is expected with the cluster emissivity and Poisson noise. In Fig. 2, we show the comparison between the S/N obtained with the absolute method and the relative method, for the four regions of analysis and for several spatial scales. For the most central regions, the S/N metric values are high and both methods perform equivalently well. For the outermost regions, the absolute method performs better than the

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Comparison of the S/N values obtained using the absolute and relative definition for surface brightness fluctuations, for different spatial scales and in the four regions defined in Table (1). The comparison is made for 100 mock spectra using the best-fit model and the exposure map of A3266, with density fluctuation parameters fixed at $\sigma_\delta = 0.32$, $\ell_{inj} = 0.05R_{500}$ and $\alpha = 11/3$.</figcaption></figure>

relative method for large spatial scales, and both perform equally well on smaller scales. With this in mind, we favour the absolute method in the following.

## 2.7. Error budget

Constraining density fluctuations requires controlling the error budget associated with surface brightness fluctuations. The Bayesian characterisation of the cluster emission (see Sect. 2.3) allows us to estimate the uncertainties related to the mean surface brightness best-fit model, in the form of posterior distribution of the model parameters. The other major source of error is the Poissonian nature of the image counts, which introduces constant power at all spatial scales.

Observables related to the velocity (or density fluctuation) field are subject to an additional variance intrinsically linked to their stochastic nature. The measured statistics of the density fluctuation field $\delta$ will vary from one realisation to another, due to finite size effects introduced by the cluster itself and the limited field of view. Any quantity that depends on $\delta$ will therefore be affected by this ‘sample variance’, including the surface brightness fluctuation power spectrum, $\mathcal{P}_{2D,\Delta}$. This variance increases at higher spatial scales due to the lack of large-scale modes and becomes predominant in the error budget at spatial scales of $\ge 50$ kpc in the structure functions of Cucchetti et al. (2019). In the following, we implement a forward modelling approach which, for the first time, accounts for the sample variance induced by the finite sizes of galaxy clusters and limited fields of view.

## 2.8. Simulation-based inference

The Bayesian framework is well suited to the integration of the error budget, which should be reflected in the likelihood of our problem. In a classical Bayesian inference, we seek to obtain the
posterior distribution $p(\theta|x)$ of the parameters $\theta$ of our model given an observation $x$, by inverting the likelihood distribution $p(x|\theta)$, namely, the probability distribution of an observation given the parameters. In practice, the posterior distribution is sampled by estimating the likelihood for many parameters with Markov chain Monte Carlo (MCMC) methods. This is the process we describe in Sect. 2.3, assuming a Poisson distribution in each pixel, which allows us to define a corresponding likelihood. However, in the case of surface brightness fluctuations, there is no simple way to define an analytical and closed form for the likelihood function of a spectrum $P_{2D,\delta}$ given the density fluctuation parameters, as we cannot grasp the underlying distribution. Therefore, the Bayesian inference of the $\hat{P}_{3D,\delta}$ parameter distributions falls under the scope of approximate Bayesian computation methods, which only require the ability to model observables $x$ (including all the sources of variance one wishes to consider) for any set of parameters $\theta$ instead of an analytical likelihood function.

In this work, we chose to use a neural network (NN) that learns the likelihood of our problem, based on simulated observations. First, we mocked many couple of observables $x_i$ using parameters $\theta_i$ drawn from the prior distribution $p(\theta)$ selected for the inference. We then trained the neural network to learn an estimator $q(x|\theta)$ of the likelihood $p(x|\theta)$ using the couples $(x_i, \theta_i)$. This approximation is accomplished by adjusting a normalising flow, which is formally an invertible bijection between two probability distributions. The normalising flow itself is built using several layers of density estimators, implemented by masked autoencoders (Germain et al. 2015), so we ended up training a masked autoregressive flow (Papamakarios et al. 2017, 2019). Once the training had converged, the posterior distribution was sampled by evaluating the likelihood with the previously trained flow, sampling for $p(\theta|x) \propto q(x|\theta)p(\theta)$ with a classical MCMC approach.

We used the sbi library (Tejero-Cantero et al. 2020) implementation to perform this simulation-based inference. The prior distribution used for the parameters of Eq. (13) are shown in Table A.1. We generated mock X-ray images with surface brightness fluctuations by projecting an emissivity field with an additional density fluctuation field, as highlighted in Eq. (6). To achieve this, we defined emissivity cubes dimensioned as the X-COP data in the $(x, y)$ directions and with the same spatial resolution, but expanded to $\pm 5R_{500}$ along the line of sight. We used the same 3D models as in Sect. 2.3, that is: density, temperature, cooling, and ellipticity, along with their best-fit parameters, as well as the individual properties of each observation, that is, the exposure maps, background maps, and $N_H$ maps, to mock the expected rest surface brightness of each cluster in the X-COP sample. The density fluctuations $\delta$ are generated by drawing a single realisation of a Gaussian random field, assuming a Kolmogorov-like power spectrum, as defined in Eq. (13). These fluctuations are co-added to the modelled emissivity and projected along the line of sight in the same standard as Eq. (7). We then choose to use the $S/N_\delta$, as defined in Eq. (14) as an observable, since this quantity reduces the proper contribution of each cluster by dividing by the spectrum of the expected emission without density fluctuations and best represents the excess power due to their presence. It is evaluated on 20 logarithmically spaced scales between the minimum and maximum scale accessible in each region, as defined in Table 1. Generating 300 000 realisations of the fluctuation field allows us to properly learn the likelihood function for each region of analysis. The posterior parameters are then sampled using the NUTS sampler, as implemented in the pyro library (Bingham et al. 2019). We

checked that the number of simulations is sufficient by testing the neural network asserting proper convergence of the recovered parameters for a NN trained with an arbitrary number of simulations. We do not observe any significant improvement in training for more than 100 000 simulations.

## 2.9. Validation

To validate our methodology, we applied it to mock clusters with known power spectrum and statistical properties similar to the actual X-COP clusters with known density fluctuation parameters and compared them with the parameters inferred via the method outlined in Sect. 2.8, in the multiple regions defined in Table 1. We performed a validation by producing, for each mock cluster, a realisation of the density fluctuations with $\sigma_\delta = 0.32$, $\ell_{\text{inj}} = 0.05 R_{500}$ and $\alpha = 11/3$ as the parameters to recover, along with 300000 simulations for parameters drawn from the priors in Table A.1. We tackled the effect of sample variance by defining a joint likelihood for all our clusters, which is equivalent to stating that each cluster is an individual observation of the same fluctuation process. Since all spatial coordinates are scaled to $R_{500}$, we can define a joint likelihood as follows:

$$ \log \mathcal{L}_{\text{joint}} = \sum_{\text{cluster}} \log \mathcal{L}_{\text{cluster}}, \qquad (15) $$

where $\mathcal{L}_{\text{cluster}}$ is the likelihood estimated with the previously trained NN for each cluster in the X-COP sample. The effectiveness of this approach is illustrated in Fig. 3. In this case, we randomly drew N clusters ten times to alleviate cluster selection effects on the reconstructed parameters. We computed the joint parameters for each draw and average them. The error on each parameter decreases with the number of clusters used to define the likelihood, dependence, as expected when adding independent observations.

# 3. Results

The data accompanying this analysis is available online on the repository<fnref n="2" /> associated with this article.

## 3.1. Mean profile and fluctuation maps

The marginalised parameter mean and $1\sigma$ dispersion are displayed in Tables E.1 and E.2. We ensure proper convergence of the Markov chains by computing the $\hat{R}$ statistic, as proposed in (Vehtari et al. 2021) and checking that $\hat{R} < 1.01$. The parameters so found are compatible with those obtained by Ghirardini et al. (2019), which are representative of X-COP clusters. The corresponding fluctuation maps computed with Eq. (8) are shown in Fig. 4. Most of these fluctuation maps show non-Gaussian features in their residual. Sub-mergers or post-mergers can be identified in A2319 and A3266. Spiral-shaped structures in the inner regions of A85, A2029, or A2142 indicate the presence of gas sloshing.

## 3.2. Density fluctuation power spectrum parameters

We calculated, for each fluctuation map, the $S/N$ between the surface brightness fluctuations and the fluctuations expected when only Poisson noise is present (see Sect. 2.6) for the four different aforementioned regions, which is shown as a plain line

<footnote n="2">² https://github.com/renecotyfanboy/turbulence_xcop</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Estimation of the posterior parameters and qualitative behaviour of the associated errors for mock observations with known parameters with an increasing number of clusters in the joint likelihood. For each point, we randomly draw N clusters ten times, compute the posterior parameter distributions and average the mean and standard deviation to alleviate the selection effect. The parameters are in each of the four regions as defined in Table (1). Left panel: mean and standard deviation of the parameters estimated for an increasing number cluster in the joint likelihood. The black line represents the true parameters. Right panel: error on each parameter compared to the number of cluster used to find the parameters, rescaled so that it is equal to 1 for N = 12. The black line shows the expected N⁻¹/² behaviour.</figcaption></figure>

in Fig. 5. Separation of the fluctuation maps into several regions provides us a way to discriminate between the different processes at work in the ICM. The central region is expected to be dominated by the presence or absence of a cool core, as well as by AGN-feedback. The second ring should be rather sensitive to the presence or not of sloshing in the core of the cluster. It is worth noting that a sloshing extending at least out to R₅₀₀ has been observed in various clusters, in particular A2142 (Rossetti et al. 2013). The two outer rings should be related to the fluctuations induced by the larger dynamical assembly of the cluster, that is, by the accretion from the cosmic web and infalling halos. We also computed the S/N for the whole region inside R₅₀₀ to get an average statistic on the clusters. Using the methodology presented in Sect. 2.8, we constrained the parameters of the density fluctuation power spectrum.

The prior distribution and posterior median and 16th–84th percentiles are highlighted in Table A.1. The posterior probabilities for the free parameters of the 3D power spectrum derived from the joint fit over the whole X-COP sample are shown in the ‘corner plot’ of Fig. 6. We also checked the proper convergence by assessing that $\dot{R} < 1.01$ for each parameter. By folding the posterior parameter distributions into our model, we obtain the resulting 2D observables accounting for the complete error budget (including the sample variance), which are shown in Fig. 5. The inner regions are marked by a high slope and low injection scale, indicating the preponderance of finite-sized structures in these regions, which corroborates with the presence of sloshing and assemblage artefacts in our fluctuation maps in Fig. 4. The estimation for all fluctuations inside R₅₀₀ should be more resilient to the central region artefacts, since these structures can be described as marginal realisations of the random field when compared to the fluctuations in the external regions. The parameters in R₅₀₀ jointly estimated on the whole sample converge to a normalisation of $\sigma_\delta \sim 0.18$, an injection of the

order of $\ell_{\text{inj}} \sim 0.4R_{500}$, with a slope of $\alpha \sim 11/3$, compatible with a pure hydrodynamical Kolmogorov cascade. This is comparable to what was determined for the Coma cluster by Schuecker et al. (2004), studying the pressure fluctuation at spatial scales between 40 and 90 kpc, as well as in Zhuravleva et al. (2015) for the Perseus cluster.

# 4. Discussions

## 4.1. Correlation with morphological indicators

The processes that lead to density fluctuations are expected to be related to the dynamical state of the clusters. Indeed, merger events or sloshing in the centre are intrinsically linked to the dynamic assembly of clusters and can be characterised by morphological and dynamical indicators. We consider various indicators that are well correlated with the dynamic state of the cluster (Lovisari et al. 2017; Campitiello et al. 2022) and we seek correlations with the density fluctuations parameters. We chose to use the following indicators: (1) the concentration parameter $c_{SB}$ compares the central emission of the cluster to its total emission. We used the definition from Campitiello et al. (2022), namely, the ratio between the surface brightness inside $0.15R_{500}$ and $R_{500}$; (2) the centroid shift $w$ relates the variation of the distance of the peak of luminosity and the centroid of the emission for a varying aperture. We used the definition and determination of this parameter from Eckert et al. (2022); (3) the Gini coefficient $G$ measures the disparity of surface brightness in the pixels of the image. This coefficient ranges from 1 to 0 for a perfectly homogeneous and perfectly inhomogeneous distribution, respectively, and is computed using the following formula:

$$G = \frac{\sum_i^n \sum_j^n |s_i - s_j|}{2n^2 \bar{s}},$$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Surface brightness fluctuation maps, $\Delta_{abs}$, from the absolute method (see Sect. 2.6) for the X-COP cluster sample. The successive contours represent the annular regions I, II, III, and IV with outer radii of 0.1, 0.25, 0.5 and $1R_{500}$, respectively. For display purposes, the images are filtered by a Gaussian kernel of 7.5". Each image is ~2.5 Mpc on a side. The images are ranked in ascending order of disturbance for the cluster as measured by the $C_Z$ coefficient (see Sect. 4.1).</figcaption></figure>

where $s_{i,j}$ is the surface brightness in the $i$-th pixel, $\bar{s}$ is the mean surface brightness and $n$ the total number of pixels used to compute $G$; (4) finally, the Zernike moments, $Z_i$, relate the decomposition of the image on a basis of orthogonal polynomials that allow us to efficiently describe the morphology of the cluster. Here, we specifically use the Zernike coefficient, $C_Z$, that reflects the asymmetry of the clusters:

$$C_Z = \sum_{n,m \neq 0} \sqrt{|Z_n^m|},$$

where $Z_n^m$ are coefficients obtained from scalar products of the image and a given Zernike polynomial, as described in Capalbo et al. (2021). The gallery of clusters in Fig. 4 is sorted by increasing $C_Z$.

We computed these morphological indicators using all the pixels in $R_{500}$, excluding the point sources and assuming spher-

ical symmetry. We accounted for the shot noise by drawing several Poisson realisations of the measured image, with centre position drawn from our best-fit posterior distribution. The resultant morphological indicators are displayed in Table 2. We quantified the correlations by drawing 2000 values of each morphological indicators and computing the Spearman correlation coefficient with the same number of density fluctuation parameters drawn from the posterior distributions. In Fig. 7, we show the correlation matrix between the density fluctuation parameters measured inside $R_{500}$ and the morphological indicators, as measured by the Spearman correlation coefficient. We observe that the best correlations appear with $\sigma_{\delta}$ and our 4 morphological indicators. The $\sigma_{\delta}$ drives the standard deviation and thus the broadness of the density fluctuation distribution. The concentration and Gini are correlated negatively with $\sigma_{\delta}$, which can be understood as the fact that less concentrated clusters are
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Power spectrum S/N as defined in Eq. (14) and measured on the surface brightness fluctuation maps with the absolute method (see Sect. 2.6) in the four regions as defined in Table (1) and inside $R_{500}$ for each cluster in the X-COP sample. Dotted line and shaded envelopes show the posterior distribution median and 16th–84th percentiles, accounting for the complete error budget for each cluster. The plain line represents the measured observable for each cluster in X-COP. The plots are ranked in ascending order of disturbance for the cluster, as measured by the $C_Z$ coefficient (see Sect. 4.1).</figcaption></figure>

generally more disturbed and therefore admit a less even distribution of its surface brightness, which goes hand in hand with the emergence of surface brightness fluctuations and thus density fluctuations. Conversely, the positive correlation with the centroid shift and the Zernike moment (which both trace the deviation of the surface brightness from a spherically symmetrical distribution) can be understood as the emergence of density and surface brightness fluctuations will disturb the symmetry of the system and therefore increase the value of these two indicators. We plot the correlation between $\sigma_\delta$ and these parameters along with the best fit as a power-law scaling in Appendix G. The range of morphological indicators we use here cannot represent the structure of the fluctuations, so it was to be expected that there would be no significant correlation with the injection scale and the spectral index. Nevertheless, they do reflect the disturbance of the surface brightness, and, in

this sense, the correlation with the normalisation of the spectrum, related to the standard deviation of the fluctuations, is self-consistent.

## 4.2. Interpretation as gas clumping

Gas clumping refers to the local deviation of the gas density from the average density expected at that location. Throughout this paper, we prefer this definition to the more specific definition of ‘accreted substructures’, so clumping defines any form of density deviation from the smooth profile. One way of interpreting our density fluctuation measure is in the form of clumping. From the 3D power spectrum of density fluctuations, we can estimate the clumping factor in our analysis regions. Similarly to Zhuravleva et al. (2015), we define the clumping factor C as
<figure><img src="image_7.png" /><figcaption>Fig. 6. Posterior distributions (left) and radial evolution (right) of the standard deviation, $\sigma_{\delta}$, the injection scale, $\ell_{\text{inj}}$, and the spectral index $\alpha$ of the density fluctuation power spectrum parameters, jointly evaluated on the whole X-COP cluster sample in the four regions as defined in Table (1) and within $R_{500}$. The colour scale matches the four regions of interest, and the turquoise distribution represents the entire $R_{500}$ region. The black dashed line represents the expected 11/3 index from Kolmogorov-Oboukhov theory.</figcaption></figure>

follows:

$$C = \frac{\langle n_e^2 \rangle}{\langle n_e \rangle^2} = 1 + \sigma_{\delta}^2, \qquad (16)$$

where $\sigma_{\delta}$ is the normalisation of the density fluctuation power spectrum as defined in Eq. (13). In Fig. 8, we show the clumping factor we obtained for the regions (II), (III), and (IV) along with a comparison with the approximate clumping factor derived for 35 clusters in Eckert et al. (2015), the clumping factor estimated in the Perseus cluster with a model-free approach by de Vries et al. (2023), and the clumping factor estimated by Angelinelli et al. (2021) for the Itasca simulated cluster sample. We see a good agreement between the three results in the inner regions. We attribute our lower values to the circularly symmetric profiles used in Eckert et al. (2015). Furthermore, we see a good agreement with the factor from de Vries et al. (2023), which was computed from surface brightness fluctuations determined without a spatial surface brightness model (in contrast to what is done in this paper). This may point to the fact that accounting for the ellipticity of the surface brightness minimises the biases introduced by the arbitrary choice of model (discussed further in Sect. 4.5). However, we observe a tension in the outer region (radii ranging from 0.5 to 1 $R_{500}$), which is presumably due to a change in the nature of the measured fluctuations. This is further discussed in Sects. 4.6 and 4.8.

### 4.3. Interpretation as turbulent motions

Assuming that turbulence in the ICM is the main cause of measured density fluctuations, Zhuravleva et al. (2014) and Gaspari et al. (2014) have shown that the characteristic amplitude of the one-component velocity, which can be related to the 1D Mach number, $M_{1D}$, is proportional to the characteristic amplitude of density fluctuations. This relationship was

also studied in the case of turbulence in a box and for stratified atmospheres (Mohapatra & Sharma 2019; Mohapatra et al. 2020, 2021) and is discussed further in the case of stratified turbulence in Sect. 4.7. Simonte et al. (2022) derived a similar relationship based on cosmological simulations. This relation is meant to link density fluctuations entirely to the random gas motions, which is more or less arguable depending on the relaxed or perturbed nature of the clusters. This is further discussed in Sect. 4.8. In this case, the density fluctuation dispersion $\sigma_{\delta}$ can be linearly related to the velocity dispersion, $\sigma_{v}$, which we refer to in terms of the 3D Mach number $M_{3D} = \sigma_{v}/c_{s}$, where $c_{s}$ is the speed of sound in the ICM:

$$M_{3D} \approx (0.63 \pm 0.04) \times \sigma_{\delta}. \qquad (17)$$

We use the relation derived for the whole sample of Simonte et al. (2022; including relaxed and disturbed clusters). This correlation was initially derived in terms of $\sigma_{v}$ instead of the Mach number, which is equivalent in the limit of an isothermal ICM, but not entirely true in our study involving large scales. The turbulent Mach numbers determined for the X-COP sample are provided in Table A.1 for the four regions of analysis and in $R_{500}$. The level of turbulence we measure is clearly subsonic, which agrees with the non-thermal pressure support estimated by Eckert+19 using the assumption of a universal gas fraction (Eckert et al. 2019), and also with direct measurements of spectral lines broadening (Sanders et al. 2011; The Hitomi Collaboration 2016). It is also in agreement with previous results from other works based on the statistics of fluctuation in density (Zhuravleva et al. 2015, 2018) and in thermodynamical quantities (Hofmann et al. 2016) in various galaxy clusters. We compute the one-component velocity, $V_{1,k}$, as defined by Zhuravleva et al. (2018) in Eq. (4) and use the same scaling relation as proposed in Zhuravleva et al. (2014) to convert the density fluctuation to velocity:

$$V_{1,k} = \frac{c_s}{(1 \pm 0.3)} \times \sqrt{4\pi k^3 \bar{P}_{3D,\delta}(k)}. \qquad (18)$$
<table><caption>**Table 2.** Concentration, $c_{SB}$, centroid shift, $w$, Gini coefficient, $G$, and Zernkie moment, $C_Z$, computed for the X-COP cluster sample, with variance from the best-fit elliptical radii and Poisson noise.</caption><thead><tr><th>Name</th><th>c<sub>SB</sub></th><th>w(×10<sup>3</sup>)</th><th>G</th><th>C<sub>Z</sub></th></tr></thead><tbody><tr><td>A1644</td><td>0.114<sup>+0.002</sup><sub>-0.001</sub></td><td>13.761<sup>+0.085</sup><sub>-0.079</sub></td><td>0.699<sup>+0.001</sup><sub>-0.001</sub></td><td>1.0<sup>+0.034</sup><sub>-0.026</sub></td></tr><tr><td>A1795</td><td>0.52<sup>+0.001</sup><sub>-0.001</sub></td><td>2.36<sup>+0.019</sup><sub>-0.019</sub></td><td>0.823<sup>+0.001</sup><sub>-0.001</sub></td><td>0.477<sup>+0.007</sup><sub>-0.008</sub></td></tr><tr><td>A2029</td><td>0.496<sup>+0.001</sup><sub>-0.001</sub></td><td>0.85<sup>+0.009</sup><sub>-0.01</sub></td><td>0.792<sup>+0.001</sup><sub>-0.001</sub></td><td>0.437<sup>+0.011</sup><sub>-0.01</sub></td></tr><tr><td>A2142</td><td>0.408<sup>+0.001</sup><sub>-0.001</sub></td><td>4.51<sup>+0.03</sup><sub>-0.028</sub></td><td>0.771<sup>+0.001</sup><sub>-0.001</sub></td><td>0.391<sup>+0.004</sup><sub>-0.004</sub></td></tr><tr><td>A2255</td><td>0.078<sup>+0.002</sup><sub>-0.002</sub></td><td>31.391<sup>+0.269</sup><sub>-0.261</sub></td><td>0.694<sup>+0.002</sup><sub>-0.002</sub></td><td>0.389<sup>+0.011</sup><sub>-0.011</sub></td></tr><tr><td>A2319</td><td>0.215<sup>+0.001</sup><sub>-0.001</sub></td><td>33.069<sup>+0.179</sup><sub>-0.176</sub></td><td>0.68<sup>+0.001</sup><sub>-0.001</sub></td><td>0.946<sup>+0.012</sup><sub>-0.013</sub></td></tr><tr><td>A3158</td><td>0.212<sup>+0.001</sup><sub>-0.001</sub></td><td>5.87<sup>+0.039</sup><sub>-0.04</sub></td><td>0.708<sup>+0.001</sup><sub>-0.001</sub></td><td>0.523<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>A3266</td><td>0.17<sup>+0.001</sup><sub>-0.001</sub></td><td>30.837<sup>+0.166</sup><sub>-0.165</sub></td><td>0.665<sup>+0.001</sup><sub>-0.001</sub></td><td>0.79<sup>+0.016</sup><sub>-0.015</sub></td></tr><tr><td>A644</td><td>0.348<sup>+0.002</sup><sub>-0.002</sub></td><td>20.971<sup>+0.151</sup><sub>-0.147</sub></td><td>0.772<sup>+0.001</sup><sub>-0.001</sub></td><td>0.457<sup>+0.007</sup><sub>-0.008</sub></td></tr><tr><td>A85</td><td>0.417<sup>+0.001</sup><sub>-0.001</sub></td><td>3.85<sup>+0.02</sup><sub>-0.019</sub></td><td>0.767<sup>+0.001</sup><sub>-0.001</sub></td><td>0.711<sup>+0.008</sup><sub>-0.009</sub></td></tr><tr><td>RXC1825</td><td>0.153<sup>+0.001</sup><sub>-0.001</sub></td><td>8.031<sup>+0.061</sup><sub>-0.062</sub></td><td>0.625<sup>+0.001</sup><sub>-0.001</sub></td><td>0.333<sup>+0.007</sup><sub>-0.007</sub></td></tr><tr><td>ZW1215</td><td>0.205<sup>+0.001</sup><sub>-0.002</sub></td><td>3.731<sup>+0.03</sup><sub>-0.03</sub></td><td>0.703<sup>+0.001</sup><sub>-0.001</sub></td><td>0.328<sup>+0.009</sup><sub>-0.009</sub></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Comparison of the clumping factor obtained with Eq. (16) for the joint fit across the X-COP sample in each annular region compared to the clumping factor obtained by Eckert et al. (2015) for a sample of 35 clusters, by de Vries et al. (2023) for the Perseus cluster and by Angelinelli et al. (2021) for 9 clusters in the ITASCA simulated clusters. The plain line and envelope represent the median and 16th–84th percentiles of the clumping factor profile.</figcaption></figure>

Chandra energy bands that are combined, whereas we use a single narrow XMM-Newton band that minimises absorption.

### 4.4. Non-thermal pressure support

Numerical simulations predict that turbulent motions should be the dominant non-thermal pressure component in galaxy clusters (Vazza et al. 2018; Angelinelli et al. 2020). From this perspective, the aim is to characterise the ratio between the non-thermal pressure of turbulent origin and the total pressure $P_{\text{NTH}}/P_{\text{Total}}$. This ratio can be expressed as a function of the 3D turbulent Mach number, $M_{3D}$ (Eckert et al. 2019), and is expressed as:

$$ \frac{P_{\text{NTH}}}{P_{\text{TOT}}} = \frac{M_{3D}^2 \gamma}{M_{3D}^2 \gamma + 3}, \qquad (19) $$

where $\gamma = 5/3$ is the polytropic index. The $P_{\text{NTH}}/P_{\text{Total}}$ ratio determined for the X-COP sample are provided in Table A.1 for the four regions of analysis and inside $R_{500}$. The value of $P_{\text{NTH}}/P_{\text{Total}}$ derived for the joint $M_{3D}$ is shown in Fig. 10 as a function of radius. It is compared to the previous estimations from Eckert et al. (2019), and theoretical predictions from numerical simulation by Angelinelli et al. (2020) and Gianfagna et al. (2021). The level of non-thermal support we obtain is consistent with previous measurements on X-COP clusters at $R_{500}$ (Eckert et al. 2019). Similarly to what these authors found, we determined a value that is a fraction of between 0.5 and $1R_{500}$ lower than what is predicted by numerical simulations. This discrepancy is enhanced towards the centre of clusters up to a factor of 10. In our study, the centre is filled by residual structures such as sloshing spirals and cool-cores, which dictate our measurements in these regions. This is further discussed in Sects. 4.6 and 4.8. This also could point to the fact that real clusters may be more thermalised than those predicted in numerical simulations (e.g., due to an incomplete implementation of the

<caption>**Fig. 7.** Correlation matrix as measured by the Spearman coefficient between the morphological indicators and the density fluctuation power spectrum parameters. We display the median and difference with the 16th–84th percentiles. The morphological parameters are defined in Sect. 4.1 and the density fluctuation parameters are defined in Eq. (13) and evaluated inside $R_{500}$ for each cluster in the X-COP sample.</caption>

We compare it to the results from Zhuravleva et al. (2018) for the clusters in common with X-COP in Fig 9, where we see a reasonable agreement. The source of differences between our measurements may originate from the different approach used to calculate the fluctuations, as Zhuravleva et al. (2018) used wider
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Comparison of the one component velocity $V_{1,k}$ from Zhuravleva et al. (2018) in the inner and outer part of the cool-cores of each cluster, and our determination in the central region (I) as defined in Table 1. The envelope represents the 16th–84th percentiles of $V_{1,k}$.</figcaption></figure>

gas physics). In addition, other physical processes, such as cosmic rays, or magnetic fields could contribute to the non-thermal pressure support (Ruszkowski et al. 2017). However, the combination of the available radio data for clusters with radio emission and/or enough sources to study Faraday rotation allows for constraints to be set on the magnetic field pressure to $\le 1-2\%$, while upper limits drawn from the non-detection of hadronic $\gamma$-rays allows putting a strong upper limit to the level of $\le 1\%$ within the virial radius of clusters. Combined, these two non-thermal pressure sources should account for $\le 2-3\%$ of the non-thermal pressure support on the ICM within the virial radius. Numerical estimations on the radial profile of these two components can be found in, for instance, Vazza et al. (2016). A new paper by Botteon et al. (2022) on A2255 may link turbulence at $R_{500}$, or beyond, to the non-thermal pressure. With a first-order estimate, and limited to this perturbed cluster, it is argued that ~10% of turbulent energy (compared to thermal energy) at $R_{500}$ should be enough to explain the diffuse emission detected by LOFAR (assuming a Fermi II acceleration process).

### 4.5. Density model dependency

The density model chosen arbitrarily necessarily induces fluctuations of its own, thus introducing a bias in the statistics of brightness fluctuations. We assessed the impact of the choice of model on the statistics of surface brightness fluctuations across the X-COP clusters. We compared three models corresponding to increasing level of fidelity: a circular $\beta$-model (Cavaliere & Fusco-Femiano 1976), an elliptical $\beta$-model, and our mean-model (defined in Sect. 2.3). Assuming that the surface brightness fluctuations are distributed according to a log-normal distribution (Kawahara et al. 2007), we fit a dispersion $\delta_{S_x}$ inside $R_{500}$:

$$ \log S_x \sim N(\log S_{x,0}, \delta_{S_x}). $$

We assume that the mean surface brightness, $S_{x,0}$, in each pixel is given, respectively, by the best-fit circular and elliptic $\beta$-model and our best-fit model (see Sect. 3 or Table A.1). The results are shown in Fig. 11. There is a systematic reduction in dispersion when a more accurate model is used, showing that the simplest way to reduce arbitrary fluctuations is to account for the elliptical shape of the surface brightness. The best improvement is

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Comparison of the fraction of turbulent to the hydrostatic pressure in the ICM as measured for the X-COP sample (black points), compared to the previous determination at $R_{500}$ by Eckert et al. (2019, green point). The blue and green line and associated shade envelope shows the predictions for non-thermal pressure support from the numerical simulations by Gianfagna et al. (2021) with the MUSIC clusters and Angelinelli et al. (2020) with the ITASCA clusters.</figcaption></figure>

the transition between the circular and elliptical $\beta$-model. This is consistent with the results by Zhuravleva et al. (2023) on numerical simulations, showing that accounting for the ellipticity of clusters can reduce the measured density fluctuations by a factor of up to 2. Moreover, we expect the underlying potential to be intrinsically triaxial (e.g., Lau et al. 2021) which motivates the smooth, elliptical surface brightness model. Including triaxiality could further improve the density modelling, but it cannot be constrained by simply using X-ray images in the [0.7–1.2] keV band. There are very few studies constraining the triaxiality using a combination of lensing, X-ray and/or Sunyaev-Zel’dovich distortion (e.g., Filippis et al. 2005; Sereno et al. 2006, 2017; Sayers et al. 2021). Applying this methodology to
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Surface brightness deviation for a circular β-model, an elliptic β-model, and the complete model from Sect. 2.3 for each cluster in the X-COP sample</figcaption></figure>

X-COP clusters would require a multi-wavelength modelling efforts that are beyond the scope of this paper. Generally speaking, structural residuals induced by modelling flaws affect the parameters of the power spectrum from Eq. (13). They tend to increase the normalisation as more fluctuations means more variance, change the slope due to the presence of sharp structures and affect the scale of injection, depending on the size of the residuals.

## 4.6. Limitations of the Gaussian random field hypothesis

Our methodology is sensitive to density fluctuations, assumed to be homogeneous and isotropic. As such, they can be described simply by a GRF. This assumption is valid as long as the density fluctuations arise from isotropic and homogeneous turbulence, and the density fluctuations are linearly related to the turbulent velocity. In the case of the X-COP clusters, the presence of cool-cores, submergers, and sloshing causes the emergence of structures in the surface brightness fluctuations that are poorly described by a GRF. Even if these structures (e.g., sloshing) cause turbulence in the ICM (ZuHone et al. 2012), we are likely to be more sensitive to fluctuations in surface brightness generated by the spiral structure itself rather than to turbulence emerging at the interface of the two phases. The same applies to the presence of hot plasma bubbles from the central AGN feedback, which causes local gas under densities in the central regions with a characteristic size $\le 10$ kpc (Zhang et al. 2022). Such spatial scales are out of the reach of XMM-Newton. These rising bubbles could induce density fluctuations, exciting g-modes and sound waves. Although simulations predict that they could contribute up to 20% of the heat dissipated in the ICM, they would remain elusive in surface brightness fluctuations (Choudhury & Reynolds 2022). In contrast, in the outer regions, turbulent processes can be expected to emerge from the dynamic assembly of structures. These processes could be sustained over long periods due to the large relaxation time frames that characterise the cluster outskirts

and the propensity for the emergence of magnetohydrodynamic instabilities (Perrone & Latter 2022).

We therefore expect the centre of the images to represent the clumping of the gas, that is, in our case, the overdensities present due to the dynamic assembly. Conversely, we can expect the outer regions to represent fluctuations due to turbulent phenomena, since they are not close enough to the centre of the cluster for substructure accretion to bias our results and the clearly identified substructures are masked as part of our analysis. To get rid of the presence of the structures in the centre would require either a fine modelling of the sloshing spiral, as well as the modelling of cool-cores and mergers in the surface brightness residues, or by considering a departure from Gaussianity and informed phase distribution for the random field model we fit, which is beyond the scope of this paper.

## 4.7. Implications of non-ideal turbulence

In galaxy clusters, the gravity influence is expected to stratify the turbulent flows. This effect will tend to suppress parallel motion and stretch the eddies in the direction orthogonal to the gravity field. Depending on the strength of the stratification, the assumption of an isotropic Gaussian field for the velocity can be rejected, since significant differences will appear between the parallel and orthogonal component (Mohapatra et al. 2020).

Furthermore, stratification will introduce a new mode of energy cascade. In addition to cascading from high to low spatial scales, kinetic energy can be converted in either direction into gravitational potential energy via buoyancy (e.g., Bolgiano 1962). By moving hot gas to regions of lower density, an additional density contrast appears and increases the total density and surface brightness fluctuations. If stratification is not negligible, not taking this effect into account may imply that the Mach number we calculate may be overestimated. The $M_{3D}/\sigma_\delta$ relationship (see Eq. (17)) proposed by Simonte et al. (2022) considers these effects in cosmological simulations and adds another degree of anisotropy by including the effect of radial accretion.

## 4.8. Limitations of $\sigma_\delta - M_{3D}$ equivalency

As in the previous section, it is expected that the density fluctuations observed in the centre of clusters are not necessarily representative of the turbulent phenomena taking place there. For instance, in Fig. 6, we observe both a low injection scale and high spectral index, which can be interpreted in a straightforward way as the presence of small-scale, sharp features in the fluctuation maps. The scale of injection then increases and the slope tends to reduce with the radius, reflecting the fact that these structures are less prevalent away from the centre. Previous works on turbulence based on numerical simulations shows that this radial evolution should be much lower for the velocity power spectrum (Vazza et al. 2012). Simonte et al. (2022) showed that the spectrum of density fluctuations was steeper in the centre of clusters than in the outer regions, and that its normalisation increases with the radius, while the velocity spectrum remained fairly invariant. At the same time, they showed a better correlation between density fluctuations and turbulent velocity fluctuations when using only relaxed clusters, similarly to Zhuravleva et al. (2023), while the correlation worsens when disturbed clusters, associated with a larger number of clumps in the central regions, are concerned. All together, this suggests that the density fluctuations below $R_{500}/2$, are probably more sensitive to the presence of dense substructures and clumps and relatively less to the compression by turbulent motions. On the other
hand, our measurement suggests that the region between $R_{500}/2$ and $R_{500}$ gives an overall better correlation between density and velocity fluctuations.

## 5. Conclusions

We performed an analysis of the spatial statistics of surface brightness fluctuations across the twelve clusters of the X-COP sample. We derived the underlying power spectrum of density fluctuations. For the first time, we have accounted for the stochastic nature of this observable by estimating the associated sample variance, thus providing a complete error budget.

- Using a simulation-based inference approach and modelling the density fluctuations as a Kolmogorov cascade, we can constrain the normalisation, slope, and scale of injection from the surface brightness fluctuations. The increasing trend in normalisation and injection scale with distance from the cluster core is consistent, respectively, with the general idea that cluster virialisation decreases and that fluctuation-generating processes have dynamic scales that increase further away from the core.

- We correlate the parameters of the density fluctuations with morphological indicators for X-COP clusters to find connections between dynamical state and apparent fluctuations. By properly propagating the uncertainties on each of our parameters, we found a clear correlation between the normalisation of the density fluctuation field and the four parameters of interest, suggesting a link between the measured normalisation of density fluctuations and the dynamical state of clusters. The normalisation of density fluctuations within $R_{500}$ is positively correlated with the Zernike coefficient and the centroid shift, and negatively correlated with the Gini coefficient and concentration parameter. All these correlations can be interpreted as being due to the fact that the amount of fluctuation increases when the surface brightness of the cluster is disturbed, which is self-consistent in our analysis.

- We further interpreted these density fluctuations as gas clumping and turbulent motions. We constrained the clumping factor, Mach number, and non-thermal pressure support in the different regions of interest of our clusters. We observe good agreement with clumping data from other works in the central parts, but not in the external parts, where it is ~10% higher. On the other hand, the non-thermal pressure profile agrees with numerical simulations in the external parts and a previous determination at $R_{500}$ for our sample, but is underestimated in the central parts up to a factor of 10. Due to the presence of assembly artefacts in the central regions, we discuss the notion that the density fluctuations in the central regions are dominated by pure clumping and density fluctuation originating from residual substructures, while turbulent motions dominate the outer fluctuations.

Considering the sample variance obviously increases the total error in the modelling and parametrisation of the power spectrum of fluctuations. However, this component of the error budget contributes significantly at most scales. Neglecting it could lead to over or false physical interpretations. For physical processes that are universal, such as turbulence, this sample variance can be tackled by increasing the number of sources. It is then assumed that each cluster of a sample would present an individual realisation of the same stochastic physical process, possibly rescaled with its typical size. Thus, we ought to extend our work to a larger sample to allow for a further refinement of the constraints on the velocity dispersion, slope, and injection scales,

or even test a more evolved turbulence model including, such as stratification or multiple injection scales (from e.g., AGN, sloshing, dynamical assembly), while retaining sufficient statistics to study relaxed and disturbed subsets to further characterise different behaviours as a function of the source dynamical state. For instance, in a future work, we will present the application of our method to the CHEX-MATE sample (Chex-Mate Collaboration 2021).

Furthermore, upcoming and future direct measurements of turbulence will be a key step in this work and will become available in the coming years, first with the Resolve instrument on board of the XRISM missions (Terada et al. 2021) and later with the X-IFU instrument on board the Athena mission (Barret et al. 2020).

**Acknowledgements.** We would like to thank Giulia Gianfagna for sharing the MUSIC cluster non-thermal pressure profile from Gianfagna et al. (2021), Irina Zhuravleva for sharing the turbulent velocity power spectra from Zhuravleva et al. (2018), Matteo Angelinelli and Thomas Jones for sharing the clumping data and non-thermal pressure profile of the ITASCA cluster sample from Angelinelli et al. (2020, 2021) and Rajsekhar Mohapatra for the insightful discussions about the impact of stratification. This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2022-22052. F.V. acknowledges financial support from the Hizon 2020 program under the ERC Starting Grant magcow, no. 714196. This work used various open-source packages such as matplotlib (Hunter 2007), astropy (Astropy Collaboration 2013, 2018), ChainConsumer (Hinton 2016), cmasher (Velden 2020), sbi (Tejero-Cantero et al. 2020), pyro (Bingham et al. 2019), jax (Bradbury et al. 2018), haiku (Hennigan et al. 2020), numpyro (Bingham et al. 2019; Phan et al. 2019).

## References

- Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197

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

- Angelinelli, M., Ettori, S., Vazza, F., & Jones, T. W. 2021, A&A, 653, A171

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

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

- Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

- Barret, D., Decourchelle, A., Fabian, A., et al. 2020, Astron. Nachr., 341, 224

- Bennett, J. S., & Sijacki, D. 2022, MNRAS, 514, 313

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

- Bingham, E., Chen, J. P., Jankowiak, M., et al. 2019, J. Mach. Learn. Res., 20, 1

- Bolgiano, R. Jr. 1962, J. Geophys. Res. (1896-1977), 67, 3015

- Botteon, A., van Weeren, R. J., Brunetti, G., et al. 2022, Sci. Adv., 8, eabq7623

- Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: Composable Transformations of Python+NumPy Programs

- Brüggen, M., & Vazza, F. 2015, Astrophys. Space Sci. Lib., 407, 599

- Campitiello, M. G., Ettori, S., Lovisari, L., et al. 2022, A&A, 665, A117

- Capalbo, V., De Petris, M., De Luca, F., et al. 2021, MNRAS, 503, 6155

- Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345

- Cavaliere, A., & Fusco-Femiano, R. 1976, A&A, 500, 95

- Chex-Mate Collaboration (Arnaud, M., et al.) 2021, A&A, 650, A104

- Choudhury, P. P., & Reynolds, C. S. 2022, MNRAS, 514, 3765

- 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

- Cucchetti, E., Clerc, N., Pointecouteau, E., Peille, P., & Pajot, F. 2019, A&A, 629, A144

- Cuciti, V., de Gasperin, F., Brüggen, M., et al. 2022, Nature, 609, 911

- de Vries, M., Mantz, A. B., Allen, S. W., et al. 2023, MNRAS, 518, 2954

- Eckert, D., Roncarelli, M., Ettori, S., et al. 2015, MNRAS, 447, 2198

- Eckert, D., Ettori, S., Pointecouteau, E., et al. 2017, Astron. Nachr., 338, 293

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

- Eckert, D., Ettori, S., Robertson, A., et al. 2022, A&A, 666, A41

- Ettori, S., & Eckert, D. 2022, A&A, 657, L1

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

- Filippis, E. D., Sereno, M., Bautz, M. W., & Longo, G. 2005, ApJ, 625, 108

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

- Gatuzz, E., Sanders, J. S., Canning, R., et al. 2022a, MNRAS, 513, 1932

- Gatuzz, E., Sanders, J. S., Dennerl, K., et al. 2022b, MNRAS, 511, 4511

- Germain, M., Gregor, K., Murray, I., & Larochelle, H. 2015, in Proceedings of the 32nd International Conference on Machine Learning (PMLR), 881
Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7

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

Hennigan, T., Cai, T., Norman, T., & Babuschkin, I. 2020, Haiku: Sonnet for JAX

HI4PI Collaboration 2016, A&A, 594, A116

Hinton, S. 2016, J. Open Source Softw., 1, 45

Hoffman, M. D., & Gelman, A. 2014, J. Mach. Learn. Res., 15, 1593

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

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Kawahara, H., Suto, Y., Kitayama, T., et al. 2007, ApJ, 659, 257

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

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

Lau, E. T., Hearin, A. P., Nagai, D., & Cappelluti, N. 2021, MNRAS, 500, 1029

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

McNamara, B. R., & Nulsen, P. E. J. 2012, New J. Phys., 14, 055023

Meneghetti, M., Rasia, E., Merten, J., et al. 2010, A&A, 514, A93

Mohapatra, R., & Sharma, P. 2019, MNRAS, 484, 4881

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

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

Mori, M., & Sugihara, M. 2001, J. Comput. Appl. Math., 127, 287

Nelson, K., Rudd, D. H., Shaw, L., & Nagai, D. 2012, ApJ, 751, 121

Nelson, K., Lau, E. T., Nagai, D., Rudd, D. H., & Yu, L. 2014, ApJ, 782, 107

Ota, N., Nagai, D., & Lau, E. T. 2018, PASJ, 70, 51

Papamakarios, G., Pavlakou, T., & Murray, I. 2017, Advances in Neural Information Processing Systems (Curran Associates, Inc.), 30

Papamakarios, G., Sterratt, D., & Murray, I. 2019, in Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics (PMLR), 837

Perrone, L. M., & Latter, H. 2022, MNRAS, 513, 4605

Phan, D., Pradhan, N., & Jankowiak, M. 2019, ArXiv e-prints [arXiv:1912.11554]

Piffaretti, R., & Valdarnini, R. 2008, A&A, 491, 71

Pinto, C., Sanders, J. S., Werner, N., et al. 2015, A&A, 575, A38

Planck Collaboration XXIX. 2014, A&A, 571, A29

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

Roncarelli, M., Gaspari, M., Ettori, S., et al. 2018, A&A, 618, A39

Rossetti, M., Eckert, D., Grandi, S. D., et al. 2013, A&A, 556, A44

Ruszkowski, M., Yang, H.-Y. K., & Reynolds, C. S. 2017, ApJ, 844, 13

Sanders, J. S., Fabian, A. C., & Smith, R. K. 2011, MNRAS, 410, 1797

Sanders, J. S., Dennerl, K., Russell, H. R., et al. 2020, A&A, 633, A42

Sayers, J., Sereno, M., Ettori, S., et al. 2021, MNRAS, 505, 4338

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

Sereno, M., Filippis, E. D., Longo, G., & Bautz, M. W. 2006, ApJ, 645, 170

Sereno, M., Ettori, S., Meneghetti, M., et al. 2017, MNRAS, 467, 3801

Shi, X., Komatsu, E., Nagai, D., & Lau, E. T. 2016, MNRAS, 455, 2936

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

Takahasi, H., & Mori, M. 1973, Publ. Res. Inst. Math. Sci., 9, 721

Tejero-Cantero, A., Boelts, J., Deistler, M., et al. 2020, J. Open Source Softw., 5, 2505

Terada, Y., Holland, M., Loewenstein, M., et al. 2021, J. Astron. Telesc. Instrum. Syst., 7, 037001

The Hitomi Collaboration 2016, Nature, 535, 117

Vazza, F., Roediger, E., & Brüggen, M. 2012, A&A, 544, A103

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

Vazza, F., Angelinelli, M., Jones, T. W., et al. 2018, MNRAS, 481, L120

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. 2021, Bayesian Anal., 16, 667

Velden, E., & v. d., 2020, J. Open Source Softw., 5, 2004

Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691

Voit, G. M., Meece, G., Li, Y., et al. 2017, ApJ, 845, 80

Weisstein, E. W. 1995, Fourier Transform (Wolfram Research, Inc.)

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

Zhang, C., Zhuravleva, I., Gendron-Marsolais, M.-L., et al. 2022, MNRAS, 517, 616

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

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

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

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

ZuHone, J. A., Markevitch, M., Brunetti, G., & Giacintucci, S. 2012, ApJ, 762, 78

ZuHone, J. A., Markevitch, M., & Zhuravleva, I. 2016, ApJ, 817, 110
Appendix A: Additional table

<table><caption>**Table A.1.** Prior distribution, inferred median and difference with 16th-84th percentiles for the density fluctuation parameters, and deduced parameters, for each cluster in the X-COP sample and for a joint fit over the whole sample. The first section represents the density fluctuation parameters defined in Sec. 2.5. The second section corresponds to the parameters deduced from the density fluctuations, namely, the clumping factor $\sqrt{C}$ (Sect. 4.2), Mach number $M$ (Sect. 4.3) and ratio of non-thermal pressure to total pressure $P_{\text{NTH}}/P_{\text{TOT}}$ (Sect. 4.4). The parameters of the density fluctuations as well as the inferred parameters are calculated in the four regions defined in Table (1) and within $R_{500}$. We ensured the proper convergence of the MCMC sampling by assessing that $\hat{R} < 1.01$ for each parameter.</caption><thead><tr><th>Name</th><th>Prior</th><th>A1644</th><th>A1795</th><th>A2029</th><th>A2142</th><th>A2255</th><th>A2319</th><th>A3158</th><th>A3266</th><th>A644</th><th>A85</th><th>RXC1825</th><th>ZW1215</th><th>Joint</th></tr></thead><tbody><tr><td>σ<sub>δ</sub> (I)</td><td>10<sup>2</sup>u(-2.1)</td><td>0.54<sup>+0.27</sup><sub>-0.18</sub></td><td>0.2<sup>+0.08</sup><sub>-0.05</sub></td><td>0.17<sup>+0.09</sup><sub>-0.06</sub></td><td>0.12<sup>+0.05</sup><sub>-0.03</sub></td><td>0.06<sup>+0.16</sup><sub>-0.04</sub></td><td>0.19<sup>+0.1</sup><sub>-0.05</sub></td><td>0.28<sup>+0.28</sup><sub>-0.14</sub></td><td>0.42<sup>+0.28</sup><sub>-0.13</sub></td><td>0.25<sup>+0.2</sup><sub>-0.11</sub></td><td>0.09<sup>+0.04</sup><sub>-0.02</sub></td><td>0.28<sup>+0.24</sup><sub>-0.1</sub></td><td>0.05<sup>+0.24</sup><sub>-0.04</sub></td><td>0.16<sup>+0.02</sup><sub>-0.01</sub></td></tr><tr><td>σ<sub>δ</sub> (II)</td><td>10<sup>2</sup>u(-2.1)</td><td>0.31<sup>+0.23</sup><sub>-0.09</sub></td><td>0.18<sup>+0.1</sup><sub>-0.04</sub></td><td>0.09<sup>+0.04</sup><sub>-0.03</sub></td><td>0.24<sup>+0.12</sup><sub>-0.06</sub></td><td>0.05<sup>+0.13</sup><sub>-0.04</sub></td><td>0.27<sup>+0.08</sup><sub>-0.05</sub></td><td>0.11<sup>+0.09</sup><sub>-0.04</sub></td><td>0.36<sup>+0.09</sup><sub>-0.05</sub></td><td>0.06<sup>+0.03</sup><sub>-0.02</sub></td><td>0.15<sup>+0.05</sup><sub>-0.03</sub></td><td>0.31<sup>+0.17</sup><sub>-0.1</sub></td><td>0.1<sup>+0.07</sup><sub>-0.04</sub></td><td>0.23<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>σ<sub>δ</sub> (III)</td><td>10<sup>2</sup>u(-2.1)</td><td>0.18<sup>+0.09</sup><sub>-0.05</sub></td><td>0.22<sup>+0.21</sup><sub>-0.09</sub></td><td>0.13<sup>+0.13</sup><sub>-0.05</sub></td><td>0.18<sup>+0.07</sup><sub>-0.04</sub></td><td>0.09<sup>+0.49</sup><sub>-0.07</sub></td><td>0.27<sup>+0.09</sup><sub>-0.06</sub></td><td>0.21<sup>+0.1</sup><sub>-0.05</sub></td><td>0.23<sup>+0.06</sup><sub>-0.04</sub></td><td>0.32<sup>+0.28</sup><sub>-0.11</sub></td><td>0.28<sup>+0.12</sup><sub>-0.08</sub></td><td>0.28<sup>+0.43</sup><sub>-0.15</sub></td><td>0.18<sup>+0.07</sup><sub>-0.04</sub></td><td>0.23<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>σ<sub>δ</sub> (IV)</td><td>10<sup>2</sup>u(-2.1)</td><td>0.74<sup>+0.78</sup><sub>-0.35</sub></td><td>0.42<sup>+0.03</sup><sub>-0.03</sub></td><td>0.33<sup>+0.1</sup><sub>-0.07</sub></td><td>0.57<sup>+0.25</sup><sub>-0.21</sub></td><td>0.23<sup>+0.84</sup><sub>-0.19</sub></td><td>0.45<sup>+0.21</sup><sub>-0.16</sub></td><td>0.45<sup>+0.21</sup><sub>-0.14</sub></td><td>0.77<sup>+0.21</sup><sub>-0.15</sub></td><td>0.48<sup>+0.1</sup><sub>-0.08</sub></td><td>0.41<sup>+0.38</sup><sub>-0.13</sub></td><td>0.39<sup>+0.26</sup><sub>-0.11</sub></td><td>0.25<sup>+0.17</sup><sub>-0.09</sub></td><td>0.59<sup>+0.08</sup><sub>-0.07</sub></td></tr><tr><td>σ<sub>δ</sub> in R<sub>500</sub></td><td>10<sup>2</sup>u(-2.1)</td><td>0.37<sup>+0.11</sup><sub>-0.09</sub></td><td>0.19<sup>+0.05</sup><sub>-0.04</sub></td><td>0.14<sup>+0.07</sup><sub>-0.04</sub></td><td>0.19<sup>+0.05</sup><sub>-0.03</sub></td><td>0.17<sup>+0.18</sup><sub>-0.06</sub></td><td>0.44<sup>+0.13</sup><sub>-0.08</sub></td><td>0.28<sup>+0.19</sup><sub>-0.12</sub></td><td>0.32<sup>+0.04</sup><sub>-0.04</sub></td><td>0.14<sup>+0.08</sup><sub>-0.03</sub></td><td>0.14<sup>+0.05</sup><sub>-0.04</sub></td><td>0.41<sup>+0.16</sup><sub>-0.14</sub></td><td>0.06<sup>+0.03</sup><sub>-0.02</sub></td><td>0.18<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>ℓ<sub>inj</sub>[R<sub>500</sub>] (I)</td><td>10<sup>2</sup>u(-2.1)</td><td>0.46<sup>+1.92</sup><sub>-0.33</sub></td><td>0.44<sup>+1.2</sup><sub>-0.28</sub></td><td>1.41<sup>+3.49</sup><sub>-1.05</sub></td><td>0.15<sup>+0.43</sup><sub>-0.05</sub></td><td>0.06<sup>+0.38</sup><sub>-0.04</sub></td><td>0.19<sup>+0.55</sup><sub>-0.08</sub></td><td>1.77<sup>+4.76</sup><sub>-1.49</sub></td><td>0.34<sup>+0.5</sup><sub>-0.15</sub></td><td>1.37<sup>+2.68</sup><sub>-0.8</sub></td><td>0.81<sup>+3.07</sup><sub>-0.65</sub></td><td>0.49<sup>+2.16</sup><sub>-0.31</sub></td><td>0.12<sup>+1.69</sup><sub>-0.1</sub></td><td>0.26<sup>+0.12</sup><sub>-0.07</sub></td></tr><tr><td>ℓ<sub>inj</sub>[R<sub>500</sub>] (II)</td><td>10<sup>2</sup>u(-2.1)</td><td>0.6<sup>+2.02</sup><sub>-0.35</sub></td><td>0.49<sup>+0.95</sup><sub>-0.22</sub></td><td>3.34<sup>+3.72</sup><sub>-2.34</sub></td><td>0.47<sup>+0.64</sup><sub>-0.2</sub></td><td>0.53<sup>+3.01</sup><sub>-0.47</sub></td><td>0.34<sup>+0.51</sup><sub>-0.12</sub></td><td>0.78<sup>+2.07</sup><sub>-0.43</sub></td><td>0.29<sup>+0.27</sup><sub>-0.09</sub></td><td>1.99<sup>+4.74</sup><sub>-1.55</sub></td><td>0.77<sup>+1.67</sup><sub>-0.42</sub></td><td>1.29<sup>+2.84</sup><sub>-0.87</sub></td><td>1.38<sup>+3.03</sup><sub>-0.9</sub></td><td>2.55<sup>+1.09</sup><sub>-0.76</sub></td></tr><tr><td>ℓ<sub>inj</sub>[R<sub>500</sub>] (III)</td><td>10<sup>2</sup>u(-2.1)</td><td>0.47<sup>+2.1</sup><sub>-0.22</sub></td><td>0.61<sup>+1.05</sup><sub>-0.31</sub></td><td>1.16<sup>+2.52</sup><sub>-0.62</sub></td><td>0.91<sup>+1.65</sup><sub>-0.44</sub></td><td>0.1<sup>+1.32</sup><sub>-0.08</sub></td><td>2.75<sup>+3.06</sup><sub>-1.55</sub></td><td>0.54<sup>+0.96</sup><sub>-0.24</sub></td><td>0.55<sup>+0.69</sup><sub>-0.21</sub></td><td>1.82<sup>+4.6</sup><sub>-1.14</sub></td><td>2.87<sup>+3.4</sup><sub>-1.63</sub></td><td>1.14<sup>+4.14</sup><sub>-0.71</sub></td><td>1.08<sup>+1.48</sup><sub>-0.52</sub></td><td>1.65<sup>+0.55</sup><sub>-0.4</sub></td></tr><tr><td>ℓ<sub>inj</sub>[R<sub>500</sub>] (IV)</td><td>10<sup>2</sup>u(-2.1)</td><td>2.48<sup>+3.91</sup><sub>-2.36</sub></td><td>0.12<sup>+0.03</sup><sub>-0.02</sub></td><td>1.22<sup>+1.37</sup><sub>-0.45</sub></td><td>5.88<sup>+2.63</sup><sub>-2.88</sub></td><td>0.35<sup>+4.45</sup><sub>-0.33</sub></td><td>4.73<sup>+3.34</sup><sub>-2.51</sub></td><td>1.71<sup>+3.41</sup><sub>-1.18</sub></td><td>1.19<sup>+2.15</sup><sub>-0.63</sub></td><td>0.32<sup>+0.43</sup><sub>-0.14</sub></td><td>2.27<sup>+4.33</sup><sub>-1.64</sub></td><td>2.59<sup>+3.23</sup><sub>-1.6</sub></td><td>1.15<sup>+3.12</sup><sub>-0.71</sub></td><td>6.53<sup>+1.45</sup><sub>-1.21</sub></td></tr><tr><td>ℓ<sub>inj</sub>[R<sub>500</sub>] in R<sub>500</sub></td><td>10<sup>2</sup>u(-2.1)</td><td>0.32<sup>+0.19</sup><sub>-0.12</sub></td><td>1.31<sup>+1.58</sup><sub>-0.61</sub></td><td>0.32<sup>+0.39</sup><sub>-0.13</sub></td><td>0.52<sup>+0.38</sup><sub>-0.16</sub></td><td>0.57<sup>+0.74</sup><sub>-0.26</sub></td><td>5.6<sup>+2.35</sup><sub>-1.82</sub></td><td>2.17<sup>+3.65</sup><sub>-1.52</sub></td><td>0.35<sup>+0.12</sup><sub>-0.08</sub></td><td>0.5<sup>+1.15</sup><sub>-0.32</sub></td><td>0.73<sup>+2.25</sup><sub>-0.44</sub></td><td>2.87<sup>+4.98</sup><sub>-2.4</sub></td><td>0.73<sup>+0.99</sup><sub>-0.28</sub></td><td>0.48<sup>+0.06</sup><sub>-0.05</sub></td></tr><tr><td>α (I)</td><td>U(0, 7)</td><td>4.0<sup>+1.24</sup><sub>-0.85</sub></td><td>2.83<sup>+0.61</sup><sub>-0.47</sub></td><td>3.18<sup>+0.61</sup><sub>-0.41</sub></td><td>5.17<sup>+1.19</sup><sub>-1.31</sub></td><td>4.02<sup>+2.12</sup><sub>-2.73</sub></td><td>4.89<sup>+1.22</sup><sub>-1.07</sub></td><td>2.9<sup>+1.37</sup><sub>-0.73</sub></td><td>5.21<sup>+1.09</sup><sub>-0.96</sub></td><td>5.86<sup>+0.76</sup><sub>-0.82</sub></td><td>3.06<sup>+0.98</sup><sub>-0.63</sub></td><td>4.97<sup>+1.41</sup><sub>-1.33</sub></td><td>1.81<sup>+2.2</sup><sub>-1.22</sub></td><td>4.08<sup>+0.35</sup><sub>-0.28</sub></td></tr><tr><td>α (II)</td><td>U(0, 7)</td><td>3.81<sup>+1.39</sup><sub>-0.73</sub></td><td>5.94<sup>+0.73</sup><sub>-0.91</sub></td><td>3.79<sup>+1.31</sup><sub>-0.7</sub></td><td>6.01<sup>+0.7</sup><sub>-1.09</sub></td><td>4.73<sup>+1.47</sup><sub>-2.56</sub></td><td>4.81<sup>+1.35</sup><sub>-0.88</sub></td><td>3.69<sup>+1.5</sup><sub>-0.95</sub></td><td>4.94<sup>+0.88</sup><sub>-0.65</sub></td><td>4.97<sup>+1.35</sup><sub>-1.44</sub></td><td>4.1<sup>+1.11</sup><sub>-0.61</sub></td><td>3.37<sup>+0.92</sup><sub>-0.51</sub></td><td>4.53<sup>+1.65</sup><sub>-1.34</sub></td><td>3.68<sup>+0.18</sup><sub>-0.14</sub></td></tr><tr><td>α (III)</td><td>U(0, 7)</td><td>4.85<sup>+1.52</sup><sub>-1.54</sub></td><td>2.96<sup>+1.03</sup><sub>-0.52</sub></td><td>3.4<sup>+1.13</sup><sub>-0.64</sub></td><td>5.13<sup>+1.15</sup><sub>-1.01</sub></td><td>2.27<sup>+2.45</sup><sub>-1.59</sub></td><td>3.54<sup>+0.59</sup><sub>-0.37</sub></td><td>4.3<sup>+1.63</sup><sub>-1.18</sub></td><td>3.97<sup>+1.15</sup><sub>-0.64</sub></td><td>4.59<sup>+1.59</sup><sub>-1.16</sub></td><td>2.96<sup>+0.37</sup><sub>-0.31</sub></td><td>2.85<sup>+1.47</sup><sub>-0.58</sub></td><td>4.77<sup>+1.22</sup><sub>-1.1</sub></td><td>3.15<sup>+0.14</sup><sub>-0.12</sub></td></tr><tr><td>α (IV)</td><td>U(0, 7)</td><td>3.72<sup>+2.12</sup><sub>-1.35</sub></td><td>5.69<sup>+0.74</sup><sub>-0.81</sub></td><td>4.85<sup>+1.49</sup><sub>-1.1</sub></td><td>3.31<sup>+0.89</sup><sub>-0.37</sub></td><td>3.54<sup>+1.99</sup><sub>-2.01</sub></td><td>4.01<sup>+1.44</sup><sub>-0.6</sub></td><td>3.05<sup>+0.99</sup><sub>-0.37</sub></td><td>2.72<sup>+0.26</sup><sub>-0.21</sub></td><td>4.48<sup>+1.51</sup><sub>-1.25</sub></td><td>4.1<sup>+2.32</sup><sub>-1.15</sub></td><td>3.93<sup>+0.81</sup><sub>-0.58</sub></td><td>4.68<sup>+1.68</sup><sub>-1.56</sub></td><td>2.9<sup>+0.1</sup><sub>-0.09</sub></td></tr><tr><td>α in R<sub>500</sub></td><td>U(0, 7)</td><td>3.4<sup>+1.01</sup><sub>-0.5</sub></td><td>3.07<sup>+0.38</sup><sub>-0.33</sub></td><td>2.9<sup>+0.66</sup><sub>-0.6</sub></td><td>4.05<sup>+0.52</sup><sub>-0.42</sub></td><td>3.5<sup>+1.22</sup><sub>-0.87</sub></td><td>3.2<sup>+0.22</sup><sub>-0.22</sub></td><td>2.83<sup>+0.83</sup><sub>-0.44</sub></td><td>4.26<sup>+0.45</sup><sub>-0.38</sub></td><td>3.32<sup>+2.45</sup><sub>-0.72</sub></td><td>3.88<sup>+1.66</sup><sub>-0.86</sub></td><td>2.81<sup>+0.37</sup><sub>-0.26</sub></td><td>4.89<sup>+1.29</sup><sub>-1.24</sub></td><td>3.67<sup>+0.12</sup><sub>-0.12</sub></td></tr><tr><td>√C (I)</td><td></td><td>1.14<sup>+0.15</sup><sub>-0.07</sub></td><td>1.02<sup>+0.02</sup><sub>-0.01</sub></td><td>1.02<sup>+0.02</sup><sub>-0.01</sub></td><td>1.01<sup>+0.01</sup><sub>-0.01</sub></td><td>1.0<sup>+0.02</sup><sub>-0.0</sub></td><td>1.02<sup>+0.02</sup><sub>-0.01</sub></td><td>1.04<sup>+0.11</sup><sub>-0.03</sub></td><td>1.09<sup>+0.14</sup><sub>-0.02</sub></td><td>1.03<sup>+0.07</sup><sub>-0.02</sub></td><td>1.0<sup>+0.0</sup><sub>-0.0</sub></td><td>1.04<sup>+0.09</sup><sub>-0.02</sub></td><td>1.0<sup>+0.04</sup><sub>-0.04</sub></td><td>1.01<sup>+0.0</sup><sub>-0.0</sub></td></tr><tr><td>√C (II)</td><td></td><td>1.05<sup>+0.09</sup><sub>-0.02</sub></td><td>1.02<sup>+0.02</sup><sub>-0.01</sub></td><td>1.0<sup>+0.0</sup><sub>-0.0</sub></td><td>1.03<sup>+0.03</sup><sub>-0.0</sub></td><td>1.0<sup>+0.02</sup><sub>-0.0</sub></td><td>1.04<sup>+0.02</sup><sub>-0.01</sub></td><td>1.01<sup>+0.01</sup><sub>-0.0</sub></td><td>1.06<sup>+0.03</sup><sub>-0.02</sub></td><td>1.0<sup>+0.0</sup><sub>-0.0</sub></td><td>1.01<sup>+0.01</sup><sub>-0.0</sub></td><td>1.05<sup>+0.06</sup><sub>-0.03</sub></td><td>1.01<sup>+0.01</sup><sub>-0.0</sub></td><td>1.03<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>√C (III)</td><td></td><td>1.02<sup>+0.02</sup><sub>-0.01</sub></td><td>1.02<sup>+0.06</sup><sub>-0.02</sub></td><td>1.01<sup>+0.03</sup><sub>-0.0</sub></td><td>1.02<sup>+0.02</sup><sub>-0.01</sub></td><td>1.0<sup>+0.15</sup><sub>-0.0</sub></td><td>1.04<sup>+0.03</sup><sub>-0.01</sub></td><td>1.02<sup>+0.03</sup><sub>-0.01</sub></td><td>1.03<sup>+0.01</sup><sub>-0.01</sub></td><td>1.05<sup>+0.12</sup><sub>-0.03</sub></td><td>1.04<sup>+0.04</sup><sub>-0.02</sub></td><td>1.04<sup>+0.19</sup><sub>-0.03</sub></td><td>1.02<sup>+0.01</sup><sub>-0.01</sub></td><td>1.03<sup>+0.01</sup><sub>-0.0</sub></td></tr><tr><td>√C (IV)</td><td></td><td>1.24<sup>+0.58</sup><sub>-0.17</sub></td><td>1.09<sup>+0.01</sup><sub>-0.01</sub></td><td>1.05<sup>+0.03</sup><sub>-0.02</sub></td><td>1.15<sup>+0.14</sup><sub>-0.09</sub></td><td>1.03<sup>+0.44</sup><sub>-0.02</sub></td><td>1.1<sup>+0.1</sup><sub>-0.05</sub></td><td>1.09<sup>+0.05</sup><sub>-0.05</sub></td><td>1.26<sup>+0.14</sup><sub>-0.09</sub></td><td>1.11<sup>+0.05</sup><sub>-0.03</sub></td><td>1.08<sup>+0.19</sup><sub>-0.04</sub></td><td>1.07<sup>+0.12</sup><sub>-0.04</sub></td><td>1.03<sup>+0.06</sup><sub>-0.02</sub></td><td>1.16<sup>+0.05</sup><sub>-0.04</sub></td></tr><tr><td>√C in R<sub>500</sub></td><td></td><td>1.07<sup>+0.04</sup><sub>-0.03</sub></td><td>1.02<sup>+0.01</sup><sub>-0.01</sub></td><td>1.01<sup>+0.01</sup><sub>-0.0</sub></td><td>1.02<sup>+0.01</sup><sub>-0.01</sub></td><td>1.01<sup>+0.05</sup><sub>-0.01</sub></td><td>1.09<sup>+0.06</sup><sub>-0.03</sub></td><td>1.04<sup>+0.07</sup><sub>-0.03</sub></td><td>1.05<sup>+0.01</sup><sub>-0.01</sub></td><td>1.01<sup>+0.02</sup><sub>-0.0</sub></td><td>1.01<sup>+0.01</sup><sub>-0.0</sub></td><td>1.08<sup>+0.07</sup><sub>-0.05</sub></td><td>1.0<sup>+0.0</sup><sub>-0.0</sub></td><td>1.02<sup>+0.0</sup><sub>-0.0</sub></td></tr><tr><td>M (I)</td><td></td><td>0.34<sup>+0.18</sup><sub>-0.11</sub></td><td>0.13<sup>+0.05</sup><sub>-0.03</sub></td><td>0.11<sup>+0.06</sup><sub>-0.04</sub></td><td>0.08<sup>+0.03</sup><sub>-0.02</sub></td><td>0.04<sup>+0.1</sup><sub>-0.02</sub></td><td>0.12<sup>+0.06</sup><sub>-0.03</sub></td><td>0.17<sup>+0.09</sup><sub>-0.09</sub></td><td>0.27<sup>+0.18</sup><sub>-0.08</sub></td><td>0.16<sup>+0.13</sup><sub>-0.07</sub></td><td>0.06<sup>+0.02</sup><sub>-0.01</sub></td><td>0.18<sup>+0.15</sup><sub>-0.07</sub></td><td>0.03<sup>+0.15</sup><sub>-0.02</sub></td><td>0.1<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>M (II)</td><td></td><td>0.19<sup>+0.15</sup><sub>-0.06</sub></td><td>0.11<sup>+0.06</sup><sub>-0.03</sub></td><td>0.06<sup>+0.03</sup><sub>-0.02</sub></td><td>0.15<sup>+0.08</sup><sub>-0.04</sub></td><td>0.03<sup>+0.09</sup><sub>-0.02</sub></td><td>0.17<sup>+0.05</sup><sub>-0.03</sub></td><td>0.07<sup>+0.06</sup><sub>-0.02</sub></td><td>0.23<sup>+0.06</sup><sub>-0.04</sub></td><td>0.04<sup>+0.02</sup><sub>-0.01</sub></td><td>0.1<sup>+0.03</sup><sub>-0.02</sub></td><td>0.19<sup>+0.11</sup><sub>-0.06</sub></td><td>0.06<sup>+0.05</sup><sub>-0.02</sub></td><td>0.15<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>M (III)</td><td></td><td>0.11<sup>+0.06</sup><sub>-0.03</sub></td><td>0.14<sup>+0.13</sup><sub>-0.06</sub></td><td>0.08<sup>+0.08</sup><sub>-0.03</sub></td><td>0.11<sup>+0.05</sup><sub>-0.03</sub></td><td>0.05<sup>+0.3</sup><sub>-0.04</sub></td><td>0.17<sup>+0.06</sup><sub>-0.04</sub></td><td>0.13<sup>+0.06</sup><sub>-0.03</sub></td><td>0.15<sup>+0.04</sup><sub>-0.02</sub></td><td>0.2<sup>+0.18</sup><sub>-0.07</sub></td><td>0.18<sup>+0.07</sup><sub>-0.05</sub></td><td>0.18<sup>+0.27</sup><sub>-0.09</sub></td><td>0.11<sup>+0.05</sup><sub>-0.02</sub></td><td>0.15<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>M (IV)</td><td></td><td>0.47<sup>+0.5</sup><sub>-0.23</sub></td><td>0.27<sup>+0.03</sup><sub>-0.03</sub></td><td>0.21<sup>+0.06</sup><sub>-0.04</sub></td><td>0.36<sup>+0.16</sup><sub>-0.13</sub></td><td>0.14<sup>+0.53</sup><sub>-0.12</sub></td><td>0.28<sup>+0.13</sup><sub>-0.1</sub></td><td>0.28<sup>+0.13</sup><sub>-0.09</sub></td><td>0.48<sup>+0.14</sup><sub>-0.1</sub></td><td>0.3<sup>+0.07</sup><sub>-0.05</sub></td><td>0.26<sup>+0.24</sup><sub>-0.08</sub></td><td>0.25<sup>+0.16</sup><sub>-0.07</sub></td><td>0.16<sup>+0.11</sup><sub>-0.06</sub></td><td>0.37<sup>+0.06</sup><sub>-0.05</sub></td></tr><tr><td>M in R<sub>500</sub></td><td></td><td>0.23<sup>+0.07</sup><sub>-0.06</sub></td><td>0.12<sup>+0.03</sup><sub>-0.02</sub></td><td>0.09<sup>+0.05</sup><sub>-0.02</sub></td><td>0.12<sup>+0.03</sup><sub>-0.02</sub></td><td>0.11<sup>+0.11</sup><sub>-0.04</sub></td><td>0.27<sup>+0.08</sup><sub>-0.05</sub></td><td>0.18<sup>+0.12</sup><sub>-0.07</sub></td><td>0.2<sup>+0.03</sup><sub>-0.03</sub></td><td>0.09<sup>+0.05</sup><sub>-0.02</sub></td><td>0.09<sup>+0.03</sup><sub>-0.02</sub></td><td>0.26<sup>+0.1</sup><sub>-0.09</sub></td><td>0.04<sup>+0.02</sup><sub>-0.01</sub></td><td>0.11<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>P<sub>NTH</sub>/P<sub>TOT</sub> (I)</td><td></td><td>5.99<sup>+6.83</sup><sub>-3.18</sub></td><td>0.87<sup>+0.83</sup><sub>-0.38</sub></td><td>0.67<sup>+0.87</sup><sub>-0.41</sub></td><td>0.32<sup>+0.34</sup><sub>-0.13</sub></td><td>0.07<sup>+0.93</sup><sub>-0.06</sub></td><td>0.83<sup>+1.04</sup><sub>-0.36</sub></td><td>1.65<sup>+4.8</sup><sub>-1.26</sub></td><td>3.79<sup>+6.12</sup><sub>-1.93</sub></td><td>1.39<sup>+2.92</sup><sub>-0.92</sub></td><td>0.19<sup>+0.2</sup><sub>-0.09</sub></td><td>1.73<sup>+3.98</sup><sub>-1.04</sub></td><td>0.06<sup>+1.78</sup><sub>-0.06</sub></td><td>0.54<sup>+0.14</sup><sub>-0.11</sub></td></tr><tr><td>P<sub>NTH</sub>/P<sub>TOT</sub> (II)</td><td></td><td>2.06<sup>+3.98</sup><sub>-1.08</sub></td><td>0.68<sup>+1.0</sup><sub>-0.3</sub></td><td>0.19<sup>+0.22</sup><sub>-0.1</sub></td><td>1.23<sup>+1.59</sup><sub>-0.54</sub></td><td>0.06<sup>+0.7</sup><sub>-0.05</sub></td><td>1.59<sup>+1.1</sup><sub>-0.52</sub></td><td>0.24<sup>+0.58</sup><sub>-0.14</sub></td><td>2.81<sup>+1.53</sup><sub>-0.77</sub></td><td>0.08<sup>+0.1</sup><sub>-0.04</sub></td><td>0.51<sup>+0.4</sup><sub>-0.19</sub></td><td>2.05<sup>+2.71</sup><sub>-1.12</sub></td><td>0.23<sup>+0.44</sup><sub>-0.14</sub></td><td>1.2<sup>+0.34</sup><sub>-0.26</sub></td></tr><tr><td>P<sub>NTH</sub>/P<sub>TOT</sub> (III)</td><td></td><td>0.71<sup>+0.87</sup><sub>-0.32</sub></td><td>1.03<sup>+2.84</sup><sub>-0.7</sub></td><td>0.38<sup>+1.12</sup><sub>-0.22</sub></td><td>0.72<sup>+0.7</sup><sub>-0.28</sub></td><td>0.17<sup>+6.57</sup><sub>-0.11</sub></td><td>1.59<sup>+1.27</sup><sub>-0.61</sub></td><td>0.97<sup>+1.09</sup><sub>-0.4</sub></td><td>1.16<sup>+0.64</sup><sub>-0.35</sub></td><td>2.22<sup>+5.17</sup><sub>-1.26</sub></td><td>1.71<sup>+1.67</sup><sub>-0.82</sub></td><td>1.75<sup>+8.31</sup><sub>-1.35</sub></td><td>0.68<sup>+0.65</sup><sub>-0.26</sub></td><td>1.19<sup>+0.33</sup><sub>-0.25</sub></td></tr><tr><td>P<sub>NTH</sub>/P<sub>TOT</sub> (IV)</td><td></td><td>10.75<sup>+23.14</sup><sub>-7.57</sub></td><td>3.79<sup>+0.76</sup><sub>-0.66</sub></td><td>2.36<sup>+1.58</sup><sub>-0.85</sub></td><td>6.53<sup>+6.4</sup><sub>-3.8</sub></td><td>1.11<sup>+19.13</sup><sub>-1.08</sub></td><td>4.19<sup>+4.51</sup><sub>-2.38</sub></td><td>4.18<sup>+4.4</sup><sub>-2.14</sub></td><td>11.58<sup>+6.14</sup><sub>-3.89</sub></td><td>4.81<sup>+2.21</sup><sub>-1.52</sub></td><td>3.53<sup>+8.66</sup><sub>-1.82</sub></td><td>3.28<sup>+5.22</sup><sub>-1.61</sub></td><td>1.39<sup>+2.42</sup><sub>-0.83</sub></td><td>7.2<sup>+2.17</sup><sub>-1.74</sub></td></tr><tr><td>P<sub>NTH</sub>/P<sub>TOT</sub> in R<sub>500</sub></td><td></td><td>2.91<sup>+1.96</sup><sub>-1.26</sub></td><td>0.81<sup>+0.52</sup><sub>-0.28</sub></td><td>0.45<sup>+0.56</sup><sub>-0.2</sub></td><td>0.79<sup>+0.5</sup><sub>-0.26</sub></td><td>0.66<sup>+2.04</sup><sub>-0.36</sub></td><td>4.02<sup>+2.63</sup><sub>-1.32</sub></td><td>1.74<sup>+2.96</sup><sub>-1.12</sub></td><td>2.14<sup>+0.68</sup><sub>-0.51</sub></td><td>0.45<sup>+0.66</sup><sub>-0.19</sub></td><td>0.44<sup>+0.37</sup><sub>-0.2</sub></td><td>3.57<sup>+3.14</sup><sub>-2.02</sub></td><td>0.09<sup>+0.11</sup><sub>-0.04</sub></td><td>0.7<sup>+0.12</sup><sub>-0.1</sub></td></tr></tbody></table>
## Appendix B: Fourier transform convention

In this paper, we define the Fourier transform with the classical signal processing convention, namely (0, -2π), see Weisstein (1995). This pairs results in the forward transform highlighted in Eq. B.1 and B.2. We use $\hat{f}$ and $\tilde{f}$ to refer, respectively, to the 2D and 3D Fourier transform of a function, $f$.

$$ \mathcal{F}\mathcal{T}_{2D}\{f\} = \int d^2\rho f(\rho)e^{-2i\pi k_{\rho}\cdot\rho} = \hat{f}(k_{\rho}), \quad (B.1) $$

$$ \mathcal{F}\mathcal{T}_{3D}\{f\} = \int d^3\mathbf{r} f(\mathbf{r})e^{-2i\pi k_{\mathbf{r}}\cdot\mathbf{r}} = \tilde{f}(k_{\mathbf{r}}). \quad (B.2) $$

## Appendix C: Functional approximation for $\Psi$

The cooling function as seen by XMM-Newton in the [0.7 - 1.2] keV energy band can be phenomenologically modelled with a smooth-broken power law for the temperature dependency and an exponential absorption for the column density:

$$ \Psi(N_H, T) \approx \Lambda_0 e^{-N_H \sigma} \left(\frac{T}{T_{\text{break}}}\right)^{-\alpha_1} \left(\frac{1}{2} + \frac{1}{2}\left(\frac{T}{T_{\text{break}}}\right)^{1/\Delta}\right)^{(\alpha_1 - \alpha_2)\Delta} . \quad (C.1) $$

To determine the best set of parameters for this analytical approximation, we minimised the least-square of the residuals between the functional given in the previous equation and the cooling function estimated with XSPEC for each cluster of the sample. The true cooling rate is determined by computing the count rate for XMM-Newton with an PhABS*APEC model. The temperature at a given radius is set using the universal profile derived from the X-COP sample by Ghirardini et al. (2019). The exposure time is arbitrarily set to 1 Ms and split among PN, MOS-1 and MOS-2 with a 60%, 20%, 20% ratio. The abundance is fixed to 0.3 Z⊙. The count-rate in an element of volume is computed in the [0.7-1.2] keV band, assuming that $n_e = 1.17 n_H$ for a fully ionised ICM (Anders & Grevesse 1989), and with the following normalisation:

$$ N = \frac{10^{-14}}{4\pi[d_A(1+z)]^2} \int d^3V n_e n_H, $$

where $d_A$ is the angular diameter distance at the redshift $z$ of the cluster. The XSPEC cooling function $\bar{\Lambda}$ is computed on a $10 \times 10$ grid with $k_B T \in [1, 10]$ keV and $N_H \in [1 \times 10^{19}, 2 \times 10^{21}]$ cm$^{-2}$ to cover the ranges expected in the X-COP clusters.

## Appendix D: Binning approximation

In Sect. 2.3, we make the approximation that the model counts in each bin is given by the product of the surface brightness, $S_X$, estimated in the centre times the effective exposure $\int_{\text{Bin}} \tau d\rho$, instead of the integral over the bin geometry:

$$ \text{Model Counts} = \int_{\text{Bin}} S_X(\rho) \tau d\rho \approx S_X(\rho_0) \times \int_{\text{Bin}} \tau d\rho, \quad (\text{D.1}) $$

where $\tau$ is the exposure time in each pixel and $d\rho$ is a surface element of the bin area. To validate this approximation, we used a A3266 image, exposure map, and binning, and we define the relative error $\epsilon$ between the true integral and the approximation as follows:

$$ \epsilon = \frac{\int_{\text{Bin}} S_X(\rho) \tau d\rho - S_X(\rho_0) \times \int_{\text{Bin}} \tau d\rho}{\int_{\text{Bin}} S_X(\rho) \tau d\rho} $$

We plot in Figure (D.1), $\epsilon$ compared to the standard deviation of $\tau d\rho$ within the bin. We observe that $\epsilon$ is contained around 0.5% for most of the bins. It rises to 4% for some bins located in areas of the image away from the centre, where the mosaic exposure varies greatly. These regions correspond to the places where the XMM-Newton mosaics overlap, and produce strong variations of exposure within the same bin. In any case, the bins with the coarsest approximation are those that contribute the least to the likelihood of our average model.

<figure><img src="image_12.png" /><figcaption>Fig. D.1. Evaluating the uncertainties induced by our approximation of a flat surface brightness in 2D spatial bins. The relative error, $\epsilon$, is compared to the standard deviation of effective exposure in each bin and the color of each point indicate the distance of the bin to the centre of the cluster image.</figcaption></figure>
# Appendix E: Mean model parameters

Here, we provide the parameters obtained for the average surface brightness profile, as presented in Sect. 2.3. The notation

$U$ represents a uniform prior distribution, while $N$ stands for a normal distribution. The density parameters are displayed in Table E.1 and the spatial parameters are displayed in Table E.2

<table><caption>**Table E.1.** Mean and standard deviation of the density model parameters from Eq. 3. The prior distributions are displayed in the first line.</caption><thead><tr><th></th><th>log<sub>10</sub>n<sup>2</sup><sub>e,0</sub></th><th>log<sub>10</sub>r<sub>c</sub></th><th>log<sub>10</sub>r<sub>s</sub></th><th>β</th><th>ε</th><th>log<sub>10</sub>B</th></tr><tr><th>Prior</th><th>U(-8, -3)</th><th>U(-2, 0)</th><th>U(-1, 1)</th><th>U(0, 5)</th><th>U(0, 5)</th><th>U(-10, -4)</th></tr></thead><tbody><tr><td>A1644</td><td>-4.057 ± 0.084</td><td>-1.934 ± 0.049</td><td>-0.091 ± 0.103</td><td>0.341 ± 0.004</td><td>3.254 ± 1.221</td><td>-7.162 ± 0.438</td></tr><tr><td>A1795</td><td>-3.574 ± 0.005</td><td>-1.475 ± 0.004</td><td>-0.254 ± 0.034</td><td>0.502 ± 0.002</td><td>2.075 ± 0.321</td><td>-7.34 ± 0.039</td></tr><tr><td>A2029</td><td>-3.241 ± 0.014</td><td>-1.912 ± 0.016</td><td>-0.988 ± 0.007</td><td>0.31 ± 0.004</td><td>1.647 ± 0.022</td><td>-7.334 ± 0.012</td></tr><tr><td>A2142</td><td>-4.267 ± 0.006</td><td>-1.274 ± 0.004</td><td>-0.08 ± 0.007</td><td>0.476 ± 0.002</td><td>4.891 ± 0.104</td><td>-7.44 ± 0.01</td></tr><tr><td>A2255</td><td>-6.09 ± 0.034</td><td>-0.6 ± 0.098</td><td>0.021 ± 0.431</td><td>0.526 ± 0.122</td><td>2.643 ± 1.49</td><td>-7.925 ± 0.65</td></tr><tr><td>A2319</td><td>-4.978 ± 0.007</td><td>-1.033 ± 0.007</td><td>-0.001 ± 0.018</td><td>0.459 ± 0.004</td><td>4.56 ± 0.365</td><td>-6.907 ± 0.019</td></tr><tr><td>A3158</td><td>-5.146 ± 0.026</td><td>-1.234 ± 0.051</td><td>-0.527 ± 0.025</td><td>0.27 ± 0.021</td><td>2.522 ± 0.115</td><td>-7.458 ± 0.055</td></tr><tr><td>A3266</td><td>-5.606 ± 0.008</td><td>-0.992 ± 0.012</td><td>-0.202 ± 0.006</td><td>0.327 ± 0.006</td><td>4.976 ± 0.025</td><td>-7.229 ± 0.011</td></tr><tr><td>A644</td><td>-4.678 ± 0.009</td><td>-1.049 ± 0.01</td><td>-0.112 ± 0.048</td><td>0.565 ± 0.009</td><td>3.914 ± 0.77</td><td>-7.486 ± 0.029</td></tr><tr><td>A85</td><td>-3.224 ± 0.008</td><td>-1.987 ± 0.006</td><td>-0.639 ± 0.008</td><td>0.357 ± 0.001</td><td>1.381 ± 0.019</td><td>-8.785 ± 0.548</td></tr><tr><td>RXC1825</td><td>-5.527 ± 0.072</td><td>-1.528 ± 0.145</td><td>-0.56 ± 0.022</td><td>0.149 ± 0.016</td><td>2.995 ± 0.098</td><td>-7.178 ± 0.023</td></tr><tr><td>ZW1215</td><td>-5.447 ± 0.019</td><td>-0.997 ± 0.037</td><td>-0.298 ± 0.059</td><td>0.412 ± 0.03</td><td>3.554 ± 0.54</td><td>-7.545 ± 0.035</td></tr></tbody></table>

<table><caption>**Table E.2.** Mean and standard deviation of the spatial parametrisation from Eq. 5. The prior distributions are displayed in the first line.</caption><thead><tr><th></th><th>θ</th><th>e</th><th>Right ascension</th><th>Declination</th></tr><tr><th>Prior</th><th>U(-π/2, +π/2)</th><th>U(0, 0.99)</th><th>RA<sub>Planck</sub> × N(1, 0.5)</th><th>DEC<sub>Planck</sub> × N(1, 0.5)</th></tr></thead><tbody><tr><td>A1644</td><td>1.148 ± 0.048</td><td>0.41 ± 0.02</td><td>12h57m10.923s ± 0.07s</td><td>-17°24′41.001” ± 0.75”</td></tr><tr><td>A1795</td><td>-0.197 ± 0.005</td><td>0.578 ± 0.002</td><td>13h48m52.778s ± 0.007s</td><td>26°35′33.404” ± 0.126”</td></tr><tr><td>A2029</td><td>-0.347 ± 0.004</td><td>0.612 ± 0.002</td><td>15h10m56.236s ± 0.005s</td><td>5°44′42.049” ± 0.087”</td></tr><tr><td>A2142</td><td>0.854 ± 0.003</td><td>0.753 ± 0.001</td><td>15h58m20.173s ± 0.013s</td><td>27°13′55.763” ± 0.17”</td></tr><tr><td>A2255</td><td>1.481 ± 0.036</td><td>0.553 ± 0.016</td><td>17h12m51.890s ± 0.397s</td><td>64°03′48.037” ± 2.149”</td></tr><tr><td>A2319</td><td>0.445 ± 0.005</td><td>0.651 ± 0.003</td><td>19h21m10.511s ± 0.038s</td><td>43°57′21.618” ± 0.527”</td></tr><tr><td>A3158</td><td>1.418 ± 0.01</td><td>0.639 ± 0.005</td><td>3h42m52.761s ± 0.09s</td><td>-53°37′44.784” ± 0.624”</td></tr><tr><td>A3266</td><td>-1.136 ± 0.006</td><td>0.58 ± 0.003</td><td>4h31m25.116s ± 0.079s</td><td>-61°25′33.411” ± 0.548”</td></tr><tr><td>A644</td><td>-0.205 ± 0.01</td><td>0.621 ± 0.005</td><td>8h17m25.447s ± 0.025s</td><td>-7°31′02.495” ± 0.459”</td></tr><tr><td>A85</td><td>0.355 ± 0.005</td><td>0.509 ± 0.002</td><td>0h41m50.414s ± 0.005s</td><td>-9°18′11.713” ± 0.074”</td></tr><tr><td>RXC1825</td><td>1.57 ± 0.0</td><td>0.624 ± 0.006</td><td>18h25m21.609s ± 0.069s</td><td>30°26′22.100” ± 0.679”</td></tr><tr><td>ZW1215</td><td>-1.094 ± 0.011</td><td>0.637 ± 0.005</td><td>12h17m41.309s ± 0.051s</td><td>3°39′30.396” ± 0.617”</td></tr></tbody></table>
Appendix F: Mexican hat filtering

We redirect the readers to the original article by Arévalo et al. (2012) for additional details on the Mexican hat filtering method. Here, we extend this formalism to arbitrary conventions for Fourier transforms (a, b), see Weisstein (1995). To do so, we first define the n-dimensional Fourier transform for a radially symmetric function f as the following Hankel transform:

$$
\mathcal{F}\mathcal{T}_{nD}\{f\} = \left(\frac{|b|}{(2\pi)^{1-a}}\right)^{n/2} \frac{(2\pi)^{n/2}}{(bk)^{n/2-1}} \int_0^\infty r^{n/2} f(r) J_{n/2-1}(bkr) dr. \quad (\text{F.1})
$$

In the following, we designate $\tilde{f}$ as the N-dimension Fourier
transform of $f$. We define the n-dimensional Gaussian kernel of
$\sigma$ standard deviation and its Fourier transform as follows:

$$
G_{\sigma,r,nD}(r) = \frac{1}{(2\pi\sigma^2)^{n/2}} \exp\left(-\frac{r^2}{2\sigma^2}\right), \quad (\text{F.2})
$$

$$
\tilde{G}_{\sigma, nD} = (2\pi)^{\frac{n(a-1)}{2}} |b|^{\frac{n}{2}} \exp\left(-\frac{b^2 k^2 \sigma^2}{2}\right). \quad (\text{F.3})
$$

The kernel of the Mexican hat is formally defined as the
difference of two Gaussian kernels, whose standard deviations
are, respectively, $\sigma_1 = \sigma / \sqrt{1 + \epsilon}$ and $\sigma_2 = \sigma \times \sqrt{1 + \epsilon}$ with
$\epsilon = 10^{-3}$. Computing the Taylor expansion of this kernel gives:

$$
\tilde{F}(k) = (\tilde{G}_{\sigma_1,nD} - \tilde{G}_{\sigma_2,nD})(k) \approx |b|^{\frac{n}{2}+2} \epsilon k^2 \sigma^2 (2\pi)^{\frac{n(a-1)}{2}} e^{-\frac{b^2 k^2 \sigma^2}{2}}. \quad (\text{F.4})
$$

The relation between the filtered frequency $k_{\max}$, which cor-
responds to the maximum value of the previous filter, and the
standard deviation $\sigma$ is given by $k_{\max} = \sqrt{2}/(|b|\sigma)$. We then
assume that we set $\sigma$ in Eq. F.4 so that $k_r$ is the frequency of
maximum value. We convolve the true image $I$ with this kernel:

$$
I_c(\mathbf{r}) = (G_{\sigma_1} - G_{\sigma_2}) * I(\mathbf{r}) \quad (\text{F.5})
$$

The variance of the convolved image $I_c$ is directly related
to the power spectrum of $I$ evaluated at $k_r$, with a simple pro-
portionality relation: $\text{Var } I_c \propto |\tilde{I}(k_r)|^2$. It can be shown using

Plancherel theorem:

$$
\int d^n \mathbf{r} |I_c(\mathbf{r})|^2 = \int d^n \mathbf{k} |\tilde{I}_c(\mathbf{k})|^2 = \int d^n \mathbf{k} |\tilde{I}(\mathbf{k})|^2 |\tilde{F}(\mathbf{k})|^2 \quad (\text{F.6})
$$

The kernel $|\tilde{F}(\mathbf{k})|^2$ is supposed to be thin enough to act like a Dirac function, and be close to zero when away from $k_r$. We can then rewrite Eq. F.6 as follows:

$$
\begin{align}
\int d^n \mathbf{r} |I_c(\mathbf{r})|^2 &\approx |\tilde{I}(k_r)|^2 \int d^n \mathbf{k} |\tilde{F}(\mathbf{k})|^2, \tag{F.7} \\
&= |\tilde{I}(k_r)|^2 \epsilon^2 \Upsilon(n) k_r^n \nonumber
\end{align}
$$

where $\Upsilon(n) = 2^{n(a-\frac{1}{2})-1}\pi^{n(a+\frac{1}{2})}n^{\frac{n}{2}+1}$ is a constant which
depends on the dimension $n$. The idea now is to incorporate an
arbitrary mask, $M$, to remove unwanted pixels on the image. We
defined the new convolved image as follows:

$$
I_c(\mathbf{r}) = \left( \frac{G_{\sigma_1} * I}{G_{\sigma_1} * M} - \frac{G_{\sigma_2} * I}{G_{\sigma_2} * M} \right) \times M. \quad (\text{F.8})
$$

Then, including these changes in Eq. F.7, we can estimate
the power spectrum of the image by computing:

$$
|\hat{I}(k_r)|^2 \approx \frac{1}{\epsilon^2 \Upsilon(n) k_r^n} \frac{\int d^n \mathbf{r} I_c(\mathbf{r})^2}{\int d^n \mathbf{r} M(\mathbf{r})}. \quad (\text{F.9})
$$

The integrals in Eq. F.9 are simply expressed as the sum over all cells (or pixels in 2D). This methodology intrinsically handles border effects as well as point source exclusions in the data. We test the effectiveness of this approach for an image generated via a GRF Φ, with the exposure masks of A3266. In Fig. F.1, we show a single realisation of this GRF along with the regions used to compute the power spectra, along with the distribution of their estimations for 1000 realisations of the GRF. We see a good agreement between the power spectrum used to generate the GRF and the reconstructed ones, with a slight bias, which is investigated in detail by Arévalo et al. (2012).
<figure><img src="image_13.png" /><figcaption>**Fig. F.1.** Assessing the efficiency of Mexican hats to compute the power spectrum of a field with arbitrary masks by computing the power spectra of 1000 realisations of a random field with known power spectrum in various ring regions. **Left panel:** regions of interest overplotted on a single realisation of the known random field with the A2319 exposure map. **Right panel:** recovered power spectra in the various regions of interest for 1000 realisations of the known random field, and comparison with the input power spectrum</figcaption></figure>
## Appendix G: Correlation between $\sigma_{\delta}$ and the morphological indicators

Here, we show the individual correlations between the normalisation of density fluctuation estimated inside $R_{500}$ and the various morphological indicators introduced in Sect. 4.1 for each cluster in the X-COP sample in Fig. G.1. We plot the best-fit distribution of a power law scaling for each correlation. This scaling is defined as $\log_{10} y = a_1 \log_{10} x + a_2$. The best-fit values of $a_1$ and $a_2$ for each correlation are displayed in Table G.1.

<table><caption>**Table G.1.** Median values of $a_1$ and $a_2$ and difference with the 16th-84th percentiles of the best fit for each scaling</caption><thead><tr><th></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th></tr></thead><tbody><tr><td>σ<sub>δ</sub> = f(c<sub>SB</sub>)</td><td>-0.198<sup>+0.097</sup><sub>-0.134</sub></td><td>0.129<sup>+0.06</sup><sub>-0.069</sub></td></tr><tr><td>σ<sub>δ</sub> = f(w)</td><td>0.123<sup>+0.057</sup><sub>-0.045</sub></td><td>0.507<sup>+0.141</sup><sub>-0.1</sub></td></tr><tr><td>σ<sub>δ</sub> = f(G)</td><td>-2.15<sup>+0.771</sup><sub>-0.988</sub></td><td>-0.05<sup>+0.101</sup><sub>-0.118</sub></td></tr><tr><td>σ<sub>δ</sub> = f(C<sub>Z</sub>)</td><td>0.37<sup>+0.174</sup><sub>-0.179</sub></td><td>0.355<sup>+0.061</sup><sub>-0.048</sub></td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. G.1.** Correlation between $\sigma_{\delta}$ evaluated in $R_{500}$ for each cluster and the morphological indicators defined in Sect. 4.1. The plain line and envelop represent the median and 16th-84th percentiles of the best fit of a power law scaling.</figcaption></figure>