© The Authors 2024

1. Introduction

Lying at the nodes of the cosmic web, massive clusters grow mainly through mergers of smaller mass units (groups and poor clusters) and through the continuous accretion of matter and field galaxies along filaments. Their matter content reflects that of the Universe (∼85% dark matter, ∼15% baryons), making them unique laboratories for testing models of gravitational structure formation and studying the thermodynamics of the intra-cluster medium (ICM). The hot (T ∼ 10⁷ − 10⁸ K) tenuous plasma, which accounts for the majority (∼85%) of the baryonic content, is responsible for the X-ray light through thermal bremsstrahlung and line emission and has been the focus of many investigations since the launch of the first X-ray satellites.

It is often assumed that, after the collapse of the main halo progenitor, the cluster gas settles in hydrostatic equilibrium into a spherically symmetric potential well. Under this assumption, all thermodynamic properties (e.g., temperature, density, and pressure) depend only on the distance from the center and are therefore homogeneous within a narrow radial shell. However, both X-ray observations and hydrodynamical simulations show that the gas is continuously perturbed by mergers and cosmological accretion (e.g., Jones & Forman 1999; Mathiesen et al. 1999; Markevitch & Vikhlinin 2007; Churazov et al. 2012; Vazza et al. 2012, 2017; Simonte et al. 2022). In addition, turbulent motions due to a large variety of other astrophysical ICM processes (e.g., sloshing, active galactic nuclei (AGN) feedback, thermal instabilities; see, e.g., Simionescu et al. 2019; Gaspari et al. 2020) can also significantly alter the hydrostatic equilibrium at different scales, from the inner core to the cluster outskirts, potentially in an anisotropic manner. Bulk motions may also evolve into turbulence and there have been several observational evidence over the last few years (e.g., Liu et al. 2015, 2016; Hitomi Collaboration 2016; Sanders et al. 2020) Thus, the complex physics of galaxy clusters may result in a high level of temperature and density substructures (i.e., fluctuations) which are tied to the dynamical history of the clusters (e.g., Simionescu et al. 2019; ZuHone & Su 2022). The level of inhomogeneities in the ICM may also depend on the thermal conductivity and the viscosity of the gas (e.g., Dolag et al. 2004; Sijacki & Springel 2006; Gaspari et al. 2014; ZuHone et al. 2015). Hence, a quantitative characterization of the thermodynamic structures in the ICM is needed before understanding a comprehensive picture of the physical processes at work in galaxy clusters.

Perturbations can appear as isobaric (e.g., slow-motion regime, i.e., low Mach numbers, or gas cooling), adiabatic (sound waves and/or weak shocks), and isothermal (characterized by strong conduction or global perturbations of gravitational potential). The nature of the perturbations is often evaluated with the “effective” equation of state δT/T = (γ − 1)δρ/ρ where γ = 0 for isobaric perturbations (i.e., temperature and density fluctuations are anticorrelated), γ = 5/3 for adiabatic perturbations (i.e., density and temperature fluctuations are positively correlated); and γ = 1 for isothermal perturbations (an unstable narrow regime). Slow motions (i.e., ℳ < 0.5) tend to be in pressure equilibrium with the medium, while stronger turbulence (i.e., 0.5 < ℳ < 1) overcomes the cluster stratification and is associated with an increased level of density fluctuations compared to the level of temperature fluctuations (i.e., higher effective γ). This behavior is likely related to the relative importance of gravity waves and sound waves with the latter showing an increasing role with the radius (e.g., Zhuravleva et al. 2013; Gaspari et al. 2014). High-resolution simulations show that plasma perturbations are tightly related to the dynamical state of the cluster and that there is an interplay between different fluctuations (e.g., Gaspari & Churazov 2013; Gaspari et al. 2014).

Regardless of their nature, all the turbulent processes are expected to generate fluctuations in ICM thermodynamic properties, that should be detectable in the related observables (see, e.g., Simionescu et al. 2019 and references therein). Various theoretical works have suggested a strong link between turbulence and thermodynamic fluctuations (e.g., Zhuravleva et al. 2014, 2023; Gaspari et al. 2014; Mohapatra et al. 2020, 2021; Simonte et al. 2022). From the observational point of view, Schuecker et al. (2004) were the first to investigate the relations between the thermodynamic fluctuations in pressure (as seen in projection from X-ray observations) and the turbulence in the ICM. More recent observational studies constrained the turbulence level by measuring the fluctuations in surface brightness or gas density in the inner regions of several galaxy clusters (e.g., Churazov et al. 2012; Sanders & Fabian 2012; Arévalo et al. 2016; Zhuravleva et al. 2016, 2018; de Vries et al. 2023) pointing to the isobaric nature of most of the total variance of perturbations. However, we should note that most of these AGN-related studies are mainly focused on the inner regions of galaxy clusters. By comparing perturbations in the central regions and in the outer regions Hofmann et al. (2016) found hints for a change in the thermodynamic state from the isobaric to the adiabatic regime.

The temperature structure can provide further information on the physics of shock-heated gas in merging events and on the role of turbulence and gas sloshing. In fact, simulations by Gaspari et al. (2014) indicate that, in the low ℳ regime, the fluctuations in temperature are as clear and robust tracers as density fluctuations, and even better than pressure fluctuations. With higher ℳ numbers (e.g., ℳ > 0.5) the amplitudes of the temperature fluctuations are a bit lower than the amplitudes of the other thermodynamic properties but are still a reasonably competitive tracer.

Thermodynamic 2D distributions can be used to identify cluster substructures and in turn to trace the merging process and the mass assembly history. By visual inspection of the thermodynamic maps of a large sample of clusters, Laganá et al. (2019) classified the clusters as relaxed or disturbed and correlated this classification with several standard cool-core diagnostic parameters. Their finding shows that the standard cool-core diagnostics are often too simplistic to account for the overall cluster dynamics because it is not rare that the cluster core can appear dynamically relaxed, while the outskirts can be dynamically unrelaxed. This is not surprising given that the relevant timescales in low density regions are much longer than in the core and that any denser infalling clump would be more prominent than in the inner regions. The thermodynamic maps can reveal in detail the complex structure of a cluster. However, something more quantitative is now needed to link the morphology and the amplitude of the temperature substructures to the cluster dynamical state.

The presence of substructures may induce biases in the determination of cluster properties, such as global gas temperature or total mass (and therefore impact the scaling relations; see Zhang et al. 2009). For instance, a disturbance (e.g., an infalling structure or a cooler or hotter spot) in a cluster may appear as a peak and/or a discontinuity in the radial profile of the scatter of the fluctuations. Unrelaxed clusters may show strong fluctuations and significant correlations between, for example, temperature and gas density fluctuations. Thus, any diagnostic of substructures, asymmetries, and turbulence at any scale is expected to be directly linked to the scatter of the scaling relations, and to the bias in X-ray hydrostatic equilibrium (HE) masses. We note that hydrodynamical simulations showed that structures in the ICM temperature distribution are, indeed, the main sources of systematic bias in HE mass estimates from X-ray data (see, e.g., Rasia et al. 2006, 2014; Ansarifard et al. 2020). Thus, measuring the level of the ICM inhomogeneity and turbulent motions is crucial to estimate the hydrostatic mass bias (e.g., Roncarelli et al. 2013; Angelinelli et al. 2020).

In this work, we study the temperature map of 28 clusters and investigate which kind of information can be extracted from their analysis. We compare them with the electron density maps with the same resolution (although in principle a significantly better resolution can be achieved) to map the inhomogeneities at the same scale. The paper is organized as follows. In Sect. 2, we present the selected sample, the data available, and the spectral analysis performed to reconstruct the 2D distribution of the gas temperature. The results on the general properties of the temperature maps and the characteristics of the observed fluctuations are presented in Sect. 3. The discussion of our main findings and the conclusions are described in Sects. 4 and 5, respectively.

Throughout the paper, we assume a flat ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹. All uncertainties are 1σ confidence intervals unless stated otherwise. We notice that several symbols are considered throughout the paper and to help the reader we provide a summary of the definitions in Table 1.

2. Data analysis

2.1. Sample

The Cluster HEritage project with XMM-Newton – Mass Assembly and Thermodynamics at the Endpoint of structure formation (CHEX-MATE¹; CHEX-MATE Collaboration 2021) is a ground-breaking multi-wavelength investigation of 61 local clusters (Tier-1: 0.05 < z < 0.2; 2 × 10¹⁴ M_⊙ < M₅₀₀ < 9 × 10¹⁴ M_⊙) and 61 massive clusters (Tier-2: z < 0.6 with M₅₀₀ > 7.25 × 10¹⁴ M_⊙); four clusters belong to both Tiers. This project aims at covering, with homogenous XMM-Newton exposures, the ICM emission up to at least R₅₀₀, for this minimally biased, signal-to-noise-limited sample of objects detected by Planck through the Sunyaev–Zeldovich effect. One of the goals of the project is the robust assessment of the level of any systematic error affecting the X-ray analysis of the cluster gas density, temperature, and mass profiles. Indeed, gas inhomogeneities prevent a robust determination of the global cluster properties and need to be properly understood to use clusters in astrophysical and cosmological studies.

We selected a pilot sample of 28 clusters (corresponding to ∼1/4 of the CHEX-MATE sample) where R₅₀₀ is completely covered by one single XMM-Newton pointing (i.e., simplifying the analysis and speeding up the fitting process) to (i) investigate how the thermodynamical maps can be used to investigate the fluctuations in the ICM and (ii) complement the information about the dynamical state of the clusters coming from the standard morphological analysis of the X-ray images presented in Campitiello et al. (2022). For this reason, the sample has been selected to cover a large range of morphological properties as estimated in CHEX-MATE (see the distribution of the concentration, c, and centroid-shift, w, parameters in Fig. 1). Although we did not attempt to be representative in terms of mass and redshift, the selected sample roughly covers the CHEX-MATE ranges. The list of clusters is provided in Table 2 and the gallery is presented in Appendix A.

Fig. 1. Distribution of the clusters in our sample in the concentration-centroid-shift (i.e., c − w) space. The values are taken from Campitiello et al. (2022). The dashed lines indicate the median values. The most relaxed clusters are in the bottom-right corner while the most disturbed in the upper-left corner.

2.2. Data reduction

Observation data files (ODFs) were retrieved from the XMM-Newton archive and reprocessed with the XMMSAS v19.0.0 with the latest calibration information available in October 2021. We used tasks emchain and epchain to generate calibrated event files from raw data. We only considered single to quadruple pixel MOS events (i.e., PATTERN < 13), and single to double pixel pn events (i.e., PATTERN < 5). Moreover, we considered only the high-quality MOS (i.e., #XMMEA_EM) and pn (i.e., FLAG==0) events. The data were cleaned for periods of high background due to the soft protons using the XMM-ESAS tools mos-filter and pn-filter, respectively. Additionally, we also removed the CCDs in the so-called ‘anomalous state’ (see Kuntz & Snowden 2008 for more details). Finally, we also obtain a list of out-of-time events, which are then subtracted from images and spectra after rescaling them based on the pn observation mode.

Point sources, identified by running the SAS wavelet detection tool ewavelet, were excluded from the analysis. To ensure a constant CXB flux across the cluster volume we excluded only the sources higher than a threshold value in the Log N-pLog S distribution. More details can be found in Bartalucci et al. (2023). We note that only point sources were masked, while clumps or large scale substructures were not removed. However, since the detection algorithm accounts for the increasing PSF at large radii we do not expect that point sources can be confused with clumps or large scale substructures in the outer regions of the clusters.

2.3. Spectral fitting

The modeling of the background was done following Lovisari & Reiprich (2019) with a few changes that we highlight here. As in Lovisari & Reiprich (2019) the cosmic X-ray background (CXB) was modeled by fitting simultaneously the XMM-Newton spectra with the ROSAT All-Sky Survey (RASS) spectra extracted from a region beyond the virial radius using the available tool² (Sabol & Snowden 2019) at the HEASARC webpage. However, we now use the counts-based spectrum, implemented with the v3.0.0, which allows us to use properly the cstat statistic during the fit. The non-vignetted quiescent particle background (QPB) was estimated using the filter wheel closed (FWC) observations after renormalizing them to match the observation of interest following the procedure presented in Zhang et al. (2009). In Lovisari & Reiprich (2019), the renormalizing values were obtained for each detector individually. However, Marelli et al. (2021) showed that the corners (i.e., the out-of-field of view regions) of the pn detector are also exposed to photons and particles concentrated by the X-ray telescope, leading to wrong renormalization values. Luckily, the authors showed that there is a very tight and linear correlation between the particle backgrounds level detected in pn and MOS2 cameras. Therefore, we applied to pn the renormalization values determined for MOS2. Marelli et al. (2021) showed that also some MOS out-of-field-of-view regions are exposed to celestial photons and provide the mask to exclude them from MOS2 data. Since the same mask cannot be applied to MOS1, we used MOS2 values as a proxy also to renormalize the MOS1 FWC datasets. Finally, we added an extra broken power-law (folded only with the RMF), to account for a residual soft proton contamination which affects many observations even after filtering the flare events (e.g., Leccardi & Molendi 2008). The slopes are fixed to 0.4 (below 5 keV) and 0.8 (above 5 keV) as described in Leccardi & Molendi (2008). Since this component may be different for MOS and pn detectors, the normalization is left free to vary in the three detectors and in all the regions of interest (this accounts, in first approximation, for the proton vignetting).

To perform the spectral analysis of our sample, we use the XSPEC fitting package (Arnaud 1996) v12.11.0k. We fit all our spectra with an APEC thermal plasma model (Smith et al. 2001) with an absorption fixed at the values provided by Bourdin et al. (2023)³. The best fits were obtained by minimizing the C-statistics (i.e., a modified Cash statistics; Cash 1979) assuming the metal abundances provided by Asplund et al. (2009). Our analysis slightly differs from the standard CHEX-MATE pipeline (described in a forthcoming work) which is partially using some ESAS tools (see Snowden et al. 2008) not optimized to work with a large number of regions, but returns consistent results.

2.4. Temperature maps

We determined the maps using the Weighted Voronoi Tesselation (WVT) method (a comparison with the results from the curvelet analysis, see Bourdin et al. 2015, is shown in Appendix B). The cells of the Voronoi map were generated using the WVT algorithm provided by Diehl & Statler (2006), which is a generalization of the Cappellari & Copin (2003) Voronoi binning algorithm, and requiring a signal-to-noise ratio S/N ∼ 30 in the 0.3–7 keV band. As shown in Lovisari & Reiprich (2019), this S/N allows us to estimate the temperatures with a statistical uncertainty of 10–20% (68% c.l., see the relative error maps in Appendix A). This choice may impact some of our results, therefore, we also produced the maps with S/N ∼ 50 (i.e., larger cells and smaller statistical uncertainties) that we use for comparison purposes. For each region of the map, we obtain a measure of the projected temperature and electron density. The projected density value for each cell of the map is given by

$$\begin{array}{c}\hfill n_{\mathrm{e}}=\sqrt{\frac{N4\pi D_{\mathrm{a}}^2{(1+z)}^2}{f_{n_{\mathrm{e}}/n_{\mathrm{H}}}{10}^{-14}V}},\end{array}$$(1)

where N is the APEC normalization, D_a the angular diameter distance, f_ne/nH is the fraction between the electron and proton density and depends on the metallicity value, and V the volume of the emitting region. The normalization N is computed as a linear combination of the individual detectors’ normalization, each weighted for the exposed detector area and the count rate contribution in the energy band of the spectral fit. The volume was determined as $V=2A\sqrt{R_{500}^2-X^2-Y^2}$, where A is the area of the region, and X and Y are the projected distances in the east-west and north-south directions, and we assumed that the properties of the material in each region are homogeneous and that there is no other material projected onto them. In Appendix A, we show the recovered maps for all the clusters.

2.5. Reconstruction of the 1D ICM profiles of temperature and density

From the 2D distributions, we calculate the projected temperature and density profiles by properly weighting each spatial cell.

The projected temperature in each annulus j is computed as

$$\begin{array}{c}\hfill T_{1\mathrm{D},j}=\frac{\sum T_{2\mathrm{D},i}{w}_i}{\sum w_i},\end{array}$$(2)

where i identifies the cell number in the 2D distributions, and w_i are the weights that take into account both the area A_ij covered by the cell i in a given annulus j and the emissivity (more details about the weighting are provided in Appendix C)

$$\begin{array}{c}\hfill w_i=A_{\mathit{ij}}{S}_{X,i}{T}_{2\mathrm{D},i}^{\alpha },\end{array}$$(3)

with S_X, i being the surface brighness in the cell i, and α = −0.75 (see, e.g., Mazzotta et al. 2004). The latter component of w_i is needed because when the overall spectrum is given by the superposition of several single-temperature spectra (e.g., the different cells), the cooler gas components are relatively more important in determining the temperature resulting from the spectral fit because the shape of the XMM-Newton responses (e.g., Mazzotta et al. 2004; Vikhlinin 2006). We note that such a value of α is strictly valid only in the cluster regime, and, if not included in the calculation, the recovered 1D temperatures tend to be overestimated (see Appendix C). In Fig. 2 we show the comparison between the profiles recovered from the maps and the profiles obtained with the direct fitting of the spectra extracted for each annulus.

Fig. 2. Comparison between the profile derived with the direct spectral fitting of each annulus j (i.e., Tspec, j). and the temperatures recovered from the maps using Eq. (2) (i.e., T1D, j) in the same spatial regions. Each blue point corresponds to an annulus of the 28 clusters, while the black points represents the average values with the shadow area including the 16th and 84th percentiles.

We also associate to each measurement of T_1D, j an uncertainty ϵ_T1D, j from the propagation of the error ϵ_T2D, i (we assume symmetric errors) on the spectral measurement T_2D, i

$$\begin{array}{c}\hfill {ϵ}_{T_{1\mathrm{D},j}}=\frac{\sqrt{\sum w_i^2{ϵ}_{T_{2\mathrm{D},i}}^2{[(1+\alpha )\sum w_i-(\alpha /T_{2\mathrm{D},i})\sum T_{2\mathrm{D},i}{w}_i]}^2}}{{(\sum w_i)}^2},\end{array}$$(4)

and a standard deviation σ_Tj (i.e., the total scatter or dispersion around the average value) defined as

$$\begin{array}{c}\hfill {\sigma }_{T_j}=\sqrt{\frac{\sum w_i{(T_{2\mathrm{D},i}-T_{1\mathrm{D},j})}^2}{\sum w_i}}.\end{array}$$(5)

The measured temperature dispersion in a given annulus is a combination of the intrinsic variation (σ_Tj, int) of the ICM temperature distribution (due to turbulent motions and mergers), of the statistical uncertainty ϵ_{T1D, j, stat} on T_1D, j (see Eq. (4)), and of the statistical uncertainty ϵ_Ti, stat of all the cells in the region of interest (as value we used the weighted mean of the statistical errors in a given region). Thus, the intrinsic temperature fluctuations are:

$$\begin{array}{c}\hfill {\sigma }_{T_{j,\mathrm{int}}}=\sqrt{{\sigma }_{T_j}^2-{ϵ}_{T_{1\mathrm{D},j,\mathrm{stat}}}^2-{ϵ}_{{\mathrm{T}}_{i,\mathrm{stat}}}^2}.\end{array}$$(6)

The gas density radial profile is recovered from the geometrical deprojection of the normalization of the spectral model (see Eq. (1)). Since the length along the line of sight is the same at a given radius (under the assumption of spherical symmetry), we can relate the line of sight averaged electron density at the annulus j to the ones measured in each cell i as

$$\begin{array}{c}\hfill n_j^2{A}_j={\displaystyle \sum }{n}_i^2{A}_{\mathit{ij}}.\end{array}$$(7)

From Eq. (7), we can then write

$$\begin{array}{cc}& n_{\mathrm{e},j}=\sqrt{\frac{\sum n_{\mathrm{e},i}^2{A}_{\mathit{ij}}}{\sum A_{\mathit{ij}}}},\hfill \\ \hfill & {ϵ}_{n_{\mathrm{e},j}}=\frac{\sqrt{\sum n_{\mathrm{e},i}^2{A}_{\mathit{ij}}^2{ϵ}_{n_{\mathrm{e},i}}^2}}{n_{\mathrm{e},j}\sum A_{\mathit{ij}}},\hfill \\ \hfill & {\sigma }_{n_{\mathrm{e},j}}^2=\frac{\sum A_{\mathit{ij}}{(n_{\mathrm{e},i}-n_{\mathrm{e},j})}^2}{\sum A_{\mathit{ij}}}.\hfill \end{array}$$(8)

We repeat the process after excluding the cells where the temperature deviates by more than 1σ from the corresponding azimuthally averaged value T_1D, j. This choice allow to remove also some of the noise in the maps which is related to the S/N criteria used to produce the maps and can be changed or improved with deeper observations. This σ-clipping is applied to the quantity

$$\begin{array}{c}\hfill s_i=\frac{T_{1\mathrm{D},j}-T_{2\mathrm{D},i}}{{({ϵ}_{T_{1\mathrm{D},j}}^2+{ϵ}_{T_{2\mathrm{D},i}}^2)}^{1/2}},\end{array}$$(9)

with T_1D, j (computed using Eq. (2)) being the temperature in the annulus encompassing the cell of interest. An example of such clipping is shown in Fig. 3. If these regions are associated to cold clumps, turbulence, and/or bulk motions as suggested by simulations, then the recovered profiles will be more representative of a component in a nearly hydrostatic equilibrium. We discuss in Appendix D the impact of the Voronoi binning in detecting the temperature inhomogeneities.

Fig. 3. Example of clipping for the cluster G041.45+29.10. In the top panel we show the distribution of si values with the regions in red and blue being the cells deviating more than 1σ from the azimuthal value. In the bottom panel, we show their distribution in the temperature map derived within R500.

3. Results on the temperature maps

In this section, we present the results of the analysis done by estimating both the global values (i.e., within R₅₀₀) of the mean temperature and of the associated scatter, and of the same quantities resolved radially. We study how these quantities distribute, how they relate among them, and if they can be used as a proxy of the dynamical state. The relations between scatter and central values of the gas temperature (and density) are then used to investigate the physical properties of the local variations.

3.1. Properties of the 2D distribution

The number of regions in the maps determined with S/N = 30 varies from ∼50 to ∼300 per cluster (median ∼150; see Table 3) depending on the quality of the data and the brightness of the cluster. The size of the cells as a function of the distance from the X-ray center (in units of R₅₀₀) is shown in Fig. 4. The distribution of the temperature and electron density values can provide some basic information about the physics of the ICM, like the constraining of the turbulent velocity in galaxy clusters. Several numerical studies suggest that the temperature and electron density distributions are approximately log-normal (e.g., Kawahara et al. 2007; Zhuravleva et al. 2013; Frank et al. 2013; Rasia et al. 2014; Gaspari et al. 2014; Towler et al. 2023). We tested whether the observed temperature distributions for our clusters differ from a normal distribution. For simplicity, we neglect the variations between the different rings and we fit the distribution for the whole cluster within R₅₀₀. In Fig. 5 (left panel), we show the measured skewness of the distribution of T_2D, i/T_1D, j (where T_1D, j is the temperature in the annulus encompassing the cell i) that should be close to zero for normal distributions. For all the clusters, we measure a positive skewness indicating a tail of higher temperature values that could be associated with heating events (e.g., shocks) that did not have time to thermalize. However, the chosen sampling of the map (i.e., the number of cells) can impact the measured values (i.e., the skewness measured in the presence of a low number of cells can be quite uncertain). To test the uncertainty associated with the number of cells present in the map for each cluster, we simulate 1000 distributions of N values (where N is equal to the number of cells in the map) randomly selected from a normal distribution. We then measure the skewness for each distribution and evaluate the dispersion around zero (i.e., around the expected value). The result is shown with green points in the left panel of Fig. 5, where typically the measured values are significantly larger than the uncertainty that can be associated with the limited numbers of cells available. That suggests that indeed the distributions have a tail of high-temperature values. The skewness values of the observed distribution do not change significantly even if we remove the core regions (i.e., < 0.15R₅₀₀), typically cooler, that potentially drive the distribution. We note that, since the temperatures T_1D, j are recovered from the measured T_2D, i values (and not from a direct spectral fitting of the region), potential biased measurements in the latter quantity will affect as well the former (i.e., it is not expected to impact significantly the skewness).

Fig. 4. Radius of the cells for the maps obtained with S/N = 30, in fraction of R500, as a function of the distance from the center. The feature observed at r/R500 ⪆ 0.65 is the result of the large cell sizes in the outer regions.

A way to graphically show the deviation from a normal distribution is through the quantile–quantile (QQ) plots, where a quantile is the fraction of points below a given value. If two sets come from a population with the same distribution, then the points should fall approximately along the 45-deg reference line. The greater the departure from this reference line, the greater the evidence that the two data sets have come from populations with different distributions. In the middle panel of Fig. 5, we compare the temperature distribution of the T_2D, i/T_1D, j values of each cluster with the normal distribution, and we see that the QQ plot for most of the clusters appears curved above the line. This shape indicates a distribution that is heavily right-skewed (i.e., hot clumps) with a light left tail (only partially associated with the cooler cells in the core region because some are associated with substructures that have merged). When comparing the logarithm of the temperatures with the lognormal distribution, we see that the QQ plots appear as roughly a straight line, suggesting that indeed the projected temperatures are roughly log-normally distributed. The Anderson-Darling test (more sensitive to the tails of the distribution with respect to the Kolmogorov–Smirnov test which is instead more sensitive to the center of distribution) indicates that the hypothesis of the normal distribution can be rejected at a significance level of 5% (1%) for 20 (14) clusters while the log-normal distribution only for 5 (4) clusters. We have also verified that similar results are obtained by excluding the innermost regions (i.e., the effect is not driven by the core where we have most of the map cells) or the outermost regions (where the cluster temperatures decline). By clipping out the regions with |s_i|> 1 the distribution becomes closer to the normal distribution, especially, if the core regions are not considered. Similar results, although less significant, are obtained with coarser but higher S/N maps (see the inset plots of Fig. 5).

Fig. 5. Observed properties of the 2D temperature distribution. Left: Measured values (blue circles) of the skewness in the temperature distribution of T2D, i/T2D, j (being T2D, j the temperature of the annulus encompassing the cell i) compared with a typical uncertainty (green squares) that can be associated with the number of cells on the temperature map. Middle: QQ plot, assuming a normal distribution, of the T2D, i/T1D, j values for the maps with S/N = 30. Right: QQ plot assuming a log-normal distribution. In the inset plots, we show the results when using the maps with S/N = 50. The units of the QQ plots are given in z-score (i.e., subtracting the mean from each data point and dividing it by the standard deviation).

3.2. Temperature dispersion and dynamical state

Radiative cooling, AGN feedback, and mergers leave clear imprints on the thermodynamic properties (e.g., temperature) of the ICM at different scales. However, the two-dimensional structures observed in the temperature maps have been mainly used so far for qualitative analysis (e.g., the determination of a shock or a cold front position) leaving their full potential unexploited.

To investigate possible connections between the temperature dispersion and the cluster dynamical state, we use the morphological parameters determined from the imaging analysis (see Campitiello et al. 2022). In particular, we consider the morphological parameters c (i.e., the ratio of the emission within two different apertures; more sensitive to the core properties), w (i.e., the variance of the projected separation between the X-ray peak and the centroid of the emission obtained within 10 apertures; more sensitive to substructures), P₂₀ (i.e., a measure of the ellipticity based on the second moment of the power ratios which consist of a 2D multipole decomposition of the surface brightness distribution), and their combination M_all. We refer to Campitiello et al. (2022) for the definition of these parameters.

For each temperature map cell, we computed s_i given in Eq. (9). We then computed for each cluster an average $\overline{|s_i|}$ value and the standard deviation of s_i. The average values are given in Table 3. In Table 4, we show the correlation of these values, in different cluster regions, with morphological parameters coming from the analysis of the X-ray images tabulated in Table 2 (i.e., computed within R₅₀₀). Apart from a moderate correlation of the spectroscopic parameters with w (more related to the gas inhomogeneities than c, which is more related to the presence of strong cores) and M_all in the intermediate regions of the clusters (i.e., 0.15–0.5R₅₀₀), there is no (or very weak) correlation between the parameters from the spectral and imaging analyses. Although there is a mild anti-correlation (i.e., r ∼ 0.4) between the spectroscopic derived parameters and the cluster temperature that could potentially hide a moderate correlation between s_i and the morphological parameters, it seems that the impact of the dynamical state on the temperature dispersion is subdominant with respect to other effects. The lack of correlation (or the difficulty to detect it) is probably related to the small range of radial temperature values (∼90% of the T_2D, i cell values lie in the range 0.5–1.5T₅₀₀), or to the size of the cells that may be sensitive only to certain fluctuations (e.g., large scale). This is different from the gas density which varies by ∼3 orders of magnitude in the same radial range and is, therefore, more sensitive to local fluctuations induced from, for example, displacement of the gas, as occurs during mergers that are expected to be the primary source of the injection of extra-energy and turbulence. This reinforces the evidence that gas density-based proxies provide more powerful estimators of the dynamical state of massive halos (thanks also to the higher resolution that can be achieved).

For each cluster we also obtained the mean relative scatter $\overline{{\sigma }_{{\mathrm{T}}_{i,\mathrm{int}}}/T_j}$ of the temperature (see Table 3) using Eq. (5) where the azimuthal value, T_j at the radius of the cell is computed using Eq. (2). Again, there is no clear correlation with the morphological parameters.

We repeated the calculation after removing all the cells with |s_i|> 1 and we report the values in Table 3. We note that, when considering the full maps, the parameters span a broad range of values, while after clipping there is a convergence toward constant values (see also Fig. 6). In particular, $\overline{{\sigma }_{T_{i,\mathrm{int}}}/T_j}$ goes to zero as expected. The average values of $\overline{|s_i|}$ also decreases by a factor of ∼2 with a much smaller standard deviation and the removal of the moderate correlation with the cluster temperature. A similar effect is seen in the values of std(s_i). It is important to note that the values decrease significantly in all clusters independently of their morphological appearance and that after clipping both $\overline{|s_i|}$ and std(s_i) tend to converge to similar values independently of the resolution of the maps (i.e., S/N = 30 vs. S/N = 50).

Fig. 6. Spectral parameters derived before (left panels) and after (right panels) clipping for the regions deviating more than 1σ from the azimuthal values. In blue and red we show the results obtained with maps of S/N = 30 and S/N = 50, respectively. In each panel we provide the correlation (i.e., Spearman rank test value r) between the spectral parameter and T, the median value of the distribution μ, and the standard deviation σ.

3.3. Impact of the temperature fluctuations in the recovered azimuthal profile

Removing the temperature inhomogeneities (associated with cold and hot clumps and substructures) should provide profiles that are closer to the expectation for gas in hydrostatic equilibrium with the gravitational potential.

The quantity s (Eq. (9)) can be used to exclude regions that deviate from the azimuthal value. To highlight the impact, we excluded all the cells with |s_i|> 1 before reconstructing the temperature profile. The fraction of cells that deviate more than 1σ from the azimuthal value show no significant correlation with the dynamical state; the Spearman rank test between the fraction of scattered temperature cells and the M_all parameter does not show evidence of a strong correlation (r ∼ 0.32; p-value = 0.10). Even morphologically relaxed systems (i.e., high c and low w or M_all) do not have T_1D, j, c/T_1D, j ∼ 1.

In Fig. 7, we show the impact on the temperature measurements after removing the cells with |s_i|> 1. In many central annuli, the new temperatures are slightly larger than the original value (a sign that we are slightly preferentially removing the, easier to detect, cold clumps). In general, within ∼0.3–0.4R₅₀₀, the effect is somehow small (i.e., < 5%) but at larger radii, the effect can be of 10–20% or more in some particular regions.

Fig. 7. Ratio between T1D, j, c (i.e., the temperature estimated in each annulus j after clipping the most deviating regions) and T1D, j (i.e., the temperature estimated using all the temperature cells in each annulus) as a function of radius. The black points represents the average values, while the shadow area includes the 16th and 84th percentiles.

This could have an obvious impact, for instance, on the mass reconstruction under the HE assumption and on the use of the overall cluster temperature in scaling relations (e.g., M_tot − T or M_tot − Y_X). In fact, the recovered total mass depends linearly on the temperature at R₅₀₀ but also on its gradient. We are not in the position to quantify the real impact with our data because at large radii the binning of the maps is quite large and prevents a proper recovering of the inhomogeneities. However, in Fig. 8 (top panel), we show how the global temperature changes after removing cells with |s_i|> 1. There are several clusters showing a difference of ∼0.5 keV, corresponding to a relative change of ∼5%. Assuming a self-similar M_tot − T scaling, this relative change in temperature would correspond to a change in the total mass of ∼7–8%. In Fig. 8 (bottom panel), we show the change in slope of the temperature profile when fitting a power-law in the region 0.15–0.75R₅₀₀ (beyond 0.75R₅₀₀ the map is too coarse). There is a very mild correlation (i.e., r = −0.24) between the change in slope and in temperature.

Fig. 8. Impact of the temperature fluctuations to the cluster overall temperature and azimuthal profile. Top: Comparison between the global temperature estimated using all the cells in the 2D maps and the temperature estimated after excluding the cells deviating more than 1σ with respect to the azimuthal value. In the inset plots we show the absolute difference between T500 and T500, c (bottom-right inset) and its fractional variation (top-left inset). Bottom: Change in the gradient of the temperature profile fitted in the region 0.15–0.75R500 before and after the clipping. In the inset plot, we show its fractional variation.

3.4. Properties of the fluctuations

Since the pioneering work by Schuecker et al. (2004) that detected a scale-invariant pressure fluctuation spectrum in the range between 40 and 90 kpc in the Coma cluster well described by a projected Kolmogorov/Oboukhov-type turbulence spectrum, several works have tried to model the expected signal associated with the propagation of large-scale eddies due to the diffusion of energy injected from dynamical phenomena like mergers and mass accretion, and have put interesting constraints from the observed fluctuations in both the density distribution of the X-ray emitting gas (e.g., Churazov et al. 2012; Arévalo et al. 2012; Gaspari & Churazov 2013; Gaspari et al. 2014; Zhuravleva et al. 2014, 2016; Hofmann et al. 2016) and in the SZ-based pressure maps (e.g., Khatri & Gaspari 2016).

In the present work, we complement previous studies with the analysis of the temperature fluctuations measured in the (projected) temperature maps, under the assumption that the projected fluctuations provide a proxy for the turbulence field (as done in, e.g., Schuecker et al. 2004; Hofmann et al. 2016; Khatri & Gaspari 2016). To support this assumption, we computed also the second-order structure function (SF) of the temperature. The SF (see a definition in, e.g., ZuHone et al. 2016) is directly related to the power spectrum of the projected quantity of interest, and has the advantage that it can be computed from the 2D distributions independently from the spatial shape of the regions under consideration. We focused on small scales (r < 150 kpc) to avoid the flattening due to the large-scale plateau (see, e.g., Roncarelli et al. 2018), and measured an average slope of the SF of ∼0.8–1. Although this is flatter than the slope of 5/3 expected for a Kolmogorov-like power spectrum, it is roughly in agreement with the finding of Roncarelli et al. (2018), where the overall shape of the SFs derived in constrained hydrodynamical simulations are less steep than the expectations from the adopted power spectra. It is known that measuring temperature fluctuations is more complicated than measuring density fluctuations (from surface brightness) and a larger S/N is required. Clearly, this translates to a coarser temperature mapping, so that different physical scales are probed with respect to the SB maps. In this case, the spectroscopically estimated T_2D, i is the result of the emission-weighted⁴ average along the line of sight of the three-dimensional temperature T (see, e.g., Mazzotta et al. 2004). For a polytropic gas, fluctuations in the latter quantity, $\delta T_i=T_i-\overline{T}$, should satisfy the relation

$$\begin{array}{c}\hfill \frac{\delta T}{T}=(\gamma -1)\frac{\delta n}{n},\end{array}$$(10)

with the fluctuations in the gas density n, where γ is the polytropic index and is equal 0 and 5/3 for the isobaric and adiabatic case, respectively. Following the arguments in Schuecker et al. (2004), we can convert these three-dimensional fluctuations into the corresponding projected quantities once the temperature is considered as a weighted quantity along the line of sight

$$\begin{array}{c}\hfill \frac{\delta T_{2\mathrm{D}}}{T_{2\mathrm{D}}}\approx \frac{\int \mathrm{d}lw\delta T}{\int \mathrm{d}lwT}\approx \frac{\gamma -1}{2}\frac{\delta n_{2\mathrm{D}}^2}{n_{2\mathrm{D}}^2},\end{array}$$(11)

where w ≈ n²T^a, $n_{\text{2D}}^2=\int \text{d}l{n}^2$, and $\delta n_{\text{2D}}^2=2n_{\text{2D}}\delta n_{\text{2D}}$. Here, we are assuming that within each cell, a single temperature and density, with their associated fluctuations, are present. It is worth noting that projection effects might be more complex than the assumption adopted here, implying the calculation of a proper window function in the band of interest for each object (see, e.g., Khatri & Gaspari 2016). However, we show in Appendix E that the relation between Eqs. (10) and (11) works pretty well for a sample of simulated objects.

In our analysis, we identify as fluctuations in temperature the dispersion (see Eq. (5)) of the measured spectral temperature around the average cluster temperature in the region of interest. Therefore, we fitted the relation:

$$\begin{array}{c}\hfill \frac{{\sigma }_{T_{j,\mathrm{int}}}}{T_j}=\frac{\gamma -1}{2}(\frac{\sum n_{\mathrm{e},i}^2{A}_{i,j}}{n_{\mathrm{e},j}^2}-1).\end{array}$$(12)

In the left panel of Fig. 9, we show the radial trend of the temperature scatter measured by fitting a power-law to the radial cells of each cluster with a size of 0.1R₅₀₀ up to a radius where the S/N = 30 maps have enough resolution (typically 0.6R₅₀₀; see Fig. 4). The temperature dispersion profiles are similar to one another and are in general quite flat (i.e., slope close to zero for most of the clusters). This result stands both for relaxed and disturbed clusters (i.e., independent of the dynamical state) and does not depend on the map resolution.

Fig. 9. Slope values of the intrinsic scatter of the temperature (left panel) and electron density (middle panel) measured within ∼0.6R500. In the right panel, we show the polytropic index γ (see Eq. (10)) derived for each cluster. Blue points are from maps with S/N = 30 and red points are from maps with S/N = 50. The solid lines represent the median value μ (also reported in red on the top-left corners) while the shaded area includes the 16th and 84th percentiles.

The fluctuations in density have been estimated similarly to the ones in temperature (see Eqs. (5) and (8)) using the line of sight averaged electron density values estimated with Eq. (1). The shape of the density distribution (e.g., the QQ plots) is similar to the one in temperature and to what is found in simulations (e.g., Zhuravleva et al. 2013) with a heavy tail in the high-density regime (i.e., well-described by a log-normal distribution) usually associated with cold regions (e.g., substructures). The absolute values of the scatter in density are quite large and are mainly driven by the poor resolution of the electron density maps (derived with the temperature map binning) paired with steep density gradients. In Appendix D, using surface brightness maps, we show how the fluctuations change as a function of the resolution. However, independent of the resolution, the scatter on the density for many clusters slightly increases with radius (see the slope values in the middle panel of Fig. 9) in agreement with the finding by simulations (see Rasia et al. 2014; Towler et al. 2023).

By fitting Eq. (12), we can make a measure of γ for each cluster. The recovered values are shown in the right panel of Fig. 9. They do not depend on the map resolution (although the S/N = 50 maps allow a better constraint thanks to their lower statistical uncertainties) or on the dynamical state (i.e., r ∼ 0.1 and r ∼ 0.2 with concentration and centroid-shift, respectively). They lie, on average, between the isobaric and the adiabatic expectations. The 1σ interval of derived γ is large and also overlaps with γ = 1 (although in our sample there are only 4 clusters with γ in the range [0.8–1.2]). However, the pure isothermal state is a very narrow, hard to achieve regime, since turbulence is continuously regenerated (see details in Gaspari et al. 2014). Obviously, the perturbations in real clusters are a mix of the different regimes and the γ values reflect such mixing. The values of γ do not show any correlation with morphological parameters, total mass, or redshift.

3.5. Global intrinsic scatter

In the previous section, we discussed the radial scatter (i.e., the dispersion in different annuli with size 0.1R₅₀₀). Here, we discuss the total dispersion within R₅₀₀ (i.e., the dispersion of all the cells around the temperature of the cluster, T₅₀₀). On top of the intrinsic variation of the ICM temperature distribution (due to turbulent motions and mergers) and the statistical uncertainties, we have also to account for the contribution due to the intrinsic radial variation over the same region. In fact, since the ICM is not isothermal, an underlying temperature profile will result in a distribution of spectroscopic measurements over the map with a given dispersion. Therefore, we compute the true intrinsic dispersion in the kT distribution as

$$\begin{array}{c}\hfill {\sigma }_{T_{\mathrm{int}}}=\sqrt{{\sigma }_{T_{\mathrm{tot}}}^2-{\sigma }_{T_{\mathrm{prof}}}^2-{ϵ}_{T_{1\mathrm{D},j,\mathrm{stat}}}^2-{ϵ}_{T_{2\mathrm{D},i,\mathrm{stat}}}^2},\end{array}$$(13)

where σ_Ttot and ϵ_{T1D, j, stat} are computed using Eqs. (4) and (5) replacing T_2D, j with the overall cluster temperature T₅₀₀, and ϵ_{T2D, i, stat} is the average statistical uncertainty of the T_2D, i measurements. For each object with a given M₅₀₀ and redshift (see Table 2), we recovered the expected projected temperature profile, σ_T, prof, by assuming the universal profile in Ghirardini et al. (2019).

In Fig. 10, we show that most of the scatter in the maps is associated with real inhomogeneities in the gas and not with the dispersion caused by the underlying profile or to the statistical uncertainties. Our results are in quite good agreement with the finding by Frank et al. (2013), although with moderately lower intrinsic scatter values. The fact that there is a significant scatter in the σ_Tint–T relation may suggest that the specific details of the cluster history (e.g., mergers, gas cooling) play a significant role in the temperature distribution. The measured scatter is weakly dependent on the number of cells in the maps (r = 0.36, p = 0.06). The positive correlation observed in Fig. 10 (left panel) may suggest that the stronger shocks and turbulence of the more massive systems are driving the large range in values of the scatter. However, the most massive clusters are also those where we expect the largest impact of substructures because of the large differences between the temperature of the infalling component with respect to the temperature of the main cluster. Nevertheless, the relative intrinsic scatter is much less dependent of the overall temperature (i.e., r ∼ −0.16, p = 0.44; see also right panel of Fig. 10) and, has values ranging between 0.11 and 0.41 (excluding the two clusters for which the statistical errors prevents a proper determination of the intrinsic scatter), suggesting that the relative intrinsic scatter is almost independent of the mass (confirmed by the low value of the Spearman coefficient, r = −0.12, p = 0.58). The relative intrinsic scatter values are also independent of the dynamical state (i.e., correlation with M_all gives r = −0.00, p = 0.98).

Fig. 10. Total temperature dispersion within R500 as a function of the overall cluster temperature. Left: Relation between the mean cluster temperature T500 and the intrinsic scatter σTint within R500. The dotted orange line represents the best-fit relation to the measured points and is compared to the observations by Frank et al. (2013; black) and simulations with (green) and without (cyan) thermal conduction by Rasia et al. (2014). Right: Relative scatter as function of the cluster temperature. In blue, we show the observed values (i.e., total scatter σTtot/T500), in green we plot the relative intrinsic scatter (i.e., σTint/T500) computed as given in Eq. (13) where the statistical errors (i.e., ϵT1D, j, stat) are shown in pink and are computed using Eq. (4), in magenta we show the scatter associated with the underlying temperature profile (i.e., σTprof/T500), and in red the average statistical uncertainty of the T2D, i measurements (i.e., ϵT2D, i, stat).

4. Discussion

Gas inhomogeneities are ubiquitous in the ICM and affect various scales (see, e.g., Nagai & Lau 2011; Roncarelli et al. 2013; Vazza et al. 2013; Towler et al. 2023). Hydrodynamical simulations show that, even for relaxed clusters, the electron density in a given radial shell can be well described by a log-normal distribution plus a possible tail (Zhuravleva et al. 2013). The former component, which accounts for a large fraction of the volume, can be considered as nearly hydrostatic and can be used, for instance, to infer the total hydrostatic mass. A similar effect is present in the temperature distribution. It is thus crucial to understand and model the temperature inhomogeneities for a better characterization of the temperature profile. In fact, if multi-phase components coexist with the virial temperature, then a spectroscopic measurement, obtained assuming a single-temperature thermal model, can be biased low with respect to the “physically-motivated” gas mass-weighted temperature within the same radial shell.

Our results show that a log-normal distribution is a good description of the 2D spectroscopic measurements of the temperatures in the ICM for most of the clusters. There are hints for a tail of high-temperature values as observed for the electron density. This is in agreement with numerical simulations (see, e.g., Kawahara et al. 2007; Zhuravleva et al. 2013) which however often use 3D information. It is worth noticing that the reconstruction of the 3D distribution from observations is limited by the geometry of the gas halo and from projection effects (significantly unknown for disturbed clusters).

Since a disturbance in the ICM is expected to contribute to the observed fluctuations and to the estimated scatter we explored if spectroscopic parameters (e.g., $\overline{|s_i|}$, std(s_i), $\overline{{\sigma }_{T_{i,\mathrm{int}}}/T_j}$) derived from the temperature maps can be used as a quantitative measure of the dynamical state of clusters, similar to what is done by morphological parameters as the surface brightness concentration and centroid-shift (among the most robust morphological parameters; see Lovisari et al. 2017; Campitiello et al. 2022). Our results provide no evidence for such correlation, possibly because the amplitude of the clumps in the ICM temperature is somewhat smaller than the ones in gas density, but also because the current resolution does not allow us to identify small features. Despite the lack of correlation between the spectrally derived parameters and the standard morphological estimators, in some cases, the spectral information can still provide complementary information about the dynamical state of clusters. A good example is G266.04−21.25 (aka “Bullet” cluster): the standard morphological parameters are unable to identify it as a disturbed system (e.g., M_all = 0.02); on the contrary, the values of the spectroscopically derived parameters are among the highest observed in our sample (see Table 3).

The measurements of the radial behavior of the temperature and density fluctuations in our sample suggest that there is a mix of isobaric and adiabatic fluctuations for most of the clusters (see Fig. 9). For comparison, the temperature and density gradients measured by Schuecker et al. (2004) in Coma suggest that the substructures are closer to the adiabatic case (see also Arévalo et al. 2016), while the finding by, for instance, Zhuravleva et al. (2018) suggests that the perturbations in the core are mainly isobaric. By comparing perturbations within and outside 100 kpc, Hofmann et al. (2016) found an indication for a change from isobaric to adiabatic perturbations. Since our maps cover the full volume within R₅₀₀ (although most of the cells are within 0.5R₅₀₀), our study of the temperature fluctuations is in agreement with previous studies and may confirm that isobaric processes dominate the core while adiabatic processes are often localized in the outskirts.

The observed temperature inhomogeneities affect the global and radial temperatures which in turn can affect the mass bias level. Rasia et al. (2012) showed that in simulations ∼10–15% of the total mass bias can be attributed to such inhomogeneities, therefore by removing them we may indeed reduce the mass bias. Because of the current quality of the data, we could only show the impact of removing 1σ level inhomogeneities (Fig. 7). With deeper exposures, providing smaller statistical uncertainties, we can relax such strong threshold. It is useful to note that since the observed inhomogeneities affect the derived integrated properties, these results not only may provide complementary information about the dynamical state of a system (e.g., the Bullet cluster), but, more importantly, can provide a statistical tool to evaluate how fluctuations impact the general behavior of the relations among these integrated quantities. The estimated level of fluctuations, like the scatter of their distribution, should be included in the statistical analysis for a better understanding of the scaling laws.

4.1. Comparison with hydrodynamical simulations

By comparing the level of temperature inhomogeneities in our cluster sample with those reconstructed in objects extracted from hydrodynamical simulations, we can give insights into the physical mechanisms at work in the ICM. In particular, we compare our results with those obtained from different simulations (FLASH, see Gaspari et al. 2014; The Three Hundred carried out with the GADGET code, The300-GX – see Cui et al. 2018; and The Three Hundred carried out with the GIZMO code, The300-GIZMO – see Cui et al. 2022) and different physics. The cosmological simulations include the structure assembly and inflows, but have low resolution. Hydrodynamical simulations test pure turbulence without substructures but have high resolution (i.e., a few kpc).

Following Rasia et al. (2014), we calculate the logarithmic gas density and temperature distributions in equal radial shells with a size of 0.1R₅₀₀. We refer to ς_ne and ς_T as the standard deviations of Gaussian distributions of the logarithmic values. Apart from some differences in smoothed-particle-hydrodynamics (SPH) and adaptive-mesh-refinement (AMR) codes (see discussion in Rasia et al. 2014 for all the ICM physics prescription, including non-radiative/NR, cooling-star-formation without/CSF, and with BH growing and AGN feedback/AGN) the agreement between observation and simulations is quite good (see Fig. 11). For the values of ς_ne and ς_T covered by our observations, it seems that the agreement is better with SPH simulations which are characterized by a higher level of temperature fluctuations.

Fig. 11. Correlation between Gaussian standard deviations ςT and ςne measured in different radial shells. Each data point corresponds to the values from a single shell of a cluster. The orange line represents the best-fit relation of the data points and is compared with the best-fits obtained by Rasia et al. (2014).

In Fig. 12, we compare the radial behavior of the ratio between the intrinsic dispersion and the median value of the ICM temperature as obtained for this CHEX-MATE subsample and a set of different hydrodynamical simulations. We also verified the impact of the Voronoi tesselation (specifically, we produce temperature and emissivity maps from the model and use the Voronoi binning adopted in the spectral analysis to recover the expected contribution of the model on σ_Tj, int/T) and found that is negligible. There is a general agreement on the trend with the radius between the constraints obtained by observations and simulations, with a very weak increase ratio moving outwards. There is also a good agreement in the level of the temperature dispersion measured in observations and simulations. We also note that, in constrained hydrodynamical simulations, a relatively high turbulence level (but still subsonic) is required to reach values closer to the observed data (and to the results of the cosmological simulations).

Fig. 12. Comparison of σTj, int/T estimated in cells of 0.1R500 as a function of radius between the CHEX-MATE objects (blue; the 16th and 84th percentiles were computed excluding the annuli for which the statistical errors prevents a proper determination of the intrinsic scatter) and results extracted from various hydrodynamical simulations (see text for details in Sect. 4.1). For cosmological simulations (i.e., The 300 GX and GIZMO), we represent the distribution for 45 objects with the median and a dispersion equal to the inter-quartile range divided by 1.35. For the FLASH simulations, we present the results obtained for 3 different projections (with different line styles) of the same object simulated with two different levels of turbulence (i.e., different Mach).

High-resolution simulations have shown a connection between such fluctuations and the Mach number of gas motion in the ICM (e.g., Gaspari & Churazov 2013; Gaspari et al. 2014). In the low Mach number regime (i.e., ℳ < 0.5) perturbations are mainly isobaric and therefore one expects σ_T/T ∼ σ_ne/n_e, while for high Mach numbers (i.e., 0.5 < ℳ < 1) the perturbations shift to the adiabatic regime implying (e.g., for γ = 5/3) σ_T/T ∼ 0.67σ_ne/n_e. The fluctuations in the ICM move from isobaric to adiabatic from the ℳ_3D = 0.25 to ℳ_3D = 0.75. This is true for subsonic turbulence but a compressional components may also lead to weak shocks and sounds waves (i.e., adiabatic fluctuations) even in low-Mach case (e.g., Mohapatra et al. 2022). However, in a seminal work, Ryu et al. (2008) showed that turbulence in clusters is largely solenoidal with subsonic velocities (see also Miniati 2014; Vazza et al. 2017). We also note that the low/high turbulence simulations by Gaspari & Churazov (2013) and Gaspari et al. (2014) mimic a relaxed/unrelaxed cluster (or merger-like vs. internal turbulence). The temperature features become more washed out in the simulated clusters with lower turbulence (i.e., relaxed system). Even though in the relaxed clusters shocks are weak (essentially sound waves, as shown by the absence of thin sheets), some extended rolls/eddies/rarefactions are still visible in the spectroscopic-like temperature map.

4.2. Role of the thermal conduction in the ICM

The level of ICM inhomogeneity may also depend on diffusive processes, such as thermal conduction. In fact, thermal conduction makes the gas more isothermal by smoothing out the temperature substructure in the ICM and therefore reducing the values of σ_T. This should be particularly true at high temperatures where conductivity, being strongly temperature dependent (k ∝ T^5/2), becomes more efficient (e.g., Dolag et al. 2004). Thus, the slope of the σ_T–T relation can be used to qualitatively constrain the degree of thermal conduction in the intracluster plasma. The observed relation stands between the values obtained from cosmological simulations with and without thermal conduction using a conductivity of 1/3 of the Spitzer value (see Fig. 10).

The link between thermal conduction and turbulence is particularly important, as turbulence can re-orient the magnetic field and thus reduce or restore heat flow to a region. Thermal conduction has been mainly investigated via theoretical studies (e.g., Dolag et al. 2004; ZuHone et al. 2013; Gaspari et al. 2014; Biffi & Valdarnini 2015) because, from an observational point of view, it is very challenging to resolve local features. However, the observations of sharp temperature gradients linked to cold fronts (Ghizzardi et al. 2010; ZuHone et al. 2013) and the presence of cool cores favor a scenario where the conduction is highly suppressed (Molendi et al. 2023, and refs therein). Our results suggest that the level of thermal conductivity, although non-zero, is smaller than 1/3 of the Spitzer value (i.e., the value assumed in the cosmological simulations used for the comparison). The low value of thermal conduction is also suggested by the hydrodynamical simulations. In fact, the red lines in Fig. 12 which are obtained with a relatively high Mach number are already unable to fit the observed data. Any extra conduction will further suppress the scatter down to lower values requiring even higher Mach numbers to match the observed trend.

4.3. Connections with turbulence and mass bias b

Our work focuses on the variations we resolve in maps of the temperature measurements that we obtain from the fitting of counts-based spectra, integrated along the line of sight. We discuss here the physical implications of such fluctuations in the spectral temperature distribution. By interpreting most of these fluctuations as generated from turbulence, induced by the mass accretion in the process of cluster formation, we can infer the ratio between the energy in turbulence and the thermal energy, and translate this ratio in terms of a predicted value of the hydrostatic mass bias b = 1 − M_HE/M_tot (see, e.g., Khatri & Gaspari 2016; Ettori & Eckert 2022).

In the adiabatic regime, we have σ_T/T ∼ 0.67σ_ne/n_e, while in the isobaric regime σ_T/T ∼ σ_ne/n_e (see Gaspari et al. 2014). Since our results suggest a mix of the two regimes (see Fig. 9), in first approximation we can take the mean of them, i.e., σ_T/T ∼ 0.83σ_ne/n_e. We can now estimate the Mach number by applying the relation ${\mathcal{M}}_{1\mathrm{D}}\sim {\sigma }_{n_{\mathrm{e}}}/n_{\mathrm{e}}\sim 1.2\ast {\sigma }_T/T$ which gives ${\mathcal{M}}_{3\mathrm{D}}=\sqrt{(}3)\ast {\mathcal{M}}_{1\mathrm{D}}\sim 2.1\ast {\sigma }_{\mathrm{T}}/T$. We measured ${\sigma }_{T,\mathrm{int}}/T=0.{17}_{-0.05}^{+0.08}$ (see Sect. 3.5 and Fig. 10). Therefore, ${\mathcal{M}}_{3\mathrm{D}}=0.{36}_{-0.09}^{+0.16}$. Hydrodynamical simulations of galaxy clusters usually find E_turb/E_th in the range 5–50% (e.g., Vazza et al. 2009; Lau et al. 2009; Gaspari et al. 2012). For $E_{\mathrm{turb}}/E_{\mathrm{th}}=0.5\gamma (\gamma -1){\mathcal{M}}_{3\mathrm{D}}^2$ that translates into a ℳ_3D from simulations in the range 0.30–0.95, consistent with our observational constraints. Our result is also in agreement with the finding of Hofmann et al. (2016) who found ℳ_3D = 0.16 − 0.40. It is worth noticing that since we do not remove the substructures it is possible that the measured temperature fluctuations are not only associated with turbulence. Therefore, the derived constraints should be considered as upper limits. However, given that the measured temperature scatter does not correlate with the dynamical state of the systems, the presence of significant substructures (which is true only for a very few systems in our sample) may not be playing a significant role. Moreover, if cosmological substructures were to dominate, ℳ numbers would appear supersonic (mimicking shocks) in the majority of systems because of the substantial boost in the density maps.

The measured ℳ_3D corresponds to $E_{\mathrm{turb}}/E_{\mathrm{th}}\sim 0.{07}_{-0.03}^{+0.09}$ which is in agreement with the observational finding of Eckert et al. (2019) and Dupourqué et al. (2023), and of the theoretical work by Angelinelli et al. (2020). As discussed in Ettori & Eckert (2022; see also Khatri & Gaspari 2016), we can then relate the quantity E_turb/E_th to the hydrostatic mass bias b as

$$\begin{array}{c}\hfill b={(\frac{E_{\mathrm{th}}}{E_{\mathrm{turb}}}+1)}^{-1}.\end{array}$$(14)

In Fig. 13, we show the estimated ratios between the turbulent and thermal energy within R₅₀₀ which points, for the studies sample, to a mass bias of $b=0.{06}_{-0.03}^{+0.07}$, comparable to that derived by the X-COP collaboration for a small sample of massive, nearby, mostly relaxed clusters (see, e.g., Ettori et al. 2019; Eckert et al. 2019).

Fig. 13. Mass bias and ratio between turbulent and thermal energy in the clusters in our sample. Top: Ratio between the turbulent and the thermal energy computed based on the σT, int/T values and under the assumption that the isobaric fluctuations are dominant. The dashed horizontal lines are computed for different Mach numbers using the following equation: $E_{\mathrm{turb}}/E_{\mathrm{th}}=0.5\gamma (\gamma -1){\mathcal{M}}_{3D}^2$. In the right y-axis, we provide the mass bias b estimated within R500 through Eq. (14). The empty squares indicate the values of Eturb/Eth (and b) after the correction accounting for the underlying power spectrum (see the text for more details). Bottom: Distribution of the values of the hydrostatic bias. In gray (blue), we show the b values before (after) accounting for the integration of the power spectrum between some characteristic scales (see text in Sect. 4.3). Dotted lines indicate the means of the distributions.

The variance in the maps is related to the turbulent energy E_turb, which is proportional to the integral of the power spectrum P(k) = k^α:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{T}}^2\propto E_{\mathrm{turb}}\propto {\displaystyle \int }P(k)k^2\mathrm{d}k={\displaystyle \int }E(k)\mathrm{d}k,\end{array}$$(15)

with the energy spectrum E(k) = k^α + 2. Equation (15) is fully generic and can be applied to any mode of fluctuation. For Kolmogorov-like turbulence (i.e., α = −11/3 leading to E(k)∝k^−5/3), we expect σ_T to increase at larger scales (i.e., for decreasing k). In our calculations, being the dispersion measured on some characteristic scales regulated by the resolution of the maps and the overall volume sampled, we have to estimate a correction to E_turb by integrating the power spectrum between a dissipation scale and an injection scale. We take these to be 10 kpc and 0.5R₅₀₀ (see Gaspari et al. 2014), respectively. The shape of the power-spectrum is assumed to be Kolmogorov-like. With the correction for the power spectrum we obtain $E_{\mathrm{turb}}/E_{\mathrm{th}}\sim 0.{12}_{-0.08}^{+0.16}$ and $b\sim 0.{11}_{-0.07}^{+0.11}$ (see Fig. 13) which tend to agree better to what derived from X-ray vs lensing studies (e.g., Sereno & Ettori 2015; Herbonnet et al. 2020; Lovisari et al. 2020a), and to recent hydrodynamical simulations (e.g., Barnes et al. 2021; Gianfagna et al. 2023). It is worth noticing that the injection and dissipation scales from simulations are not well known and depends on sub-grid models.

5. Conclusions

We investigate the level of inhomogeneities in the gas temperature (and gas density) 2D distribution for a sample of 28 CHEX-MATE galaxy clusters. The temperature maps show clear structures that we associate with fluctuations in the ICM of heterogeneous origin (mostly due to the ongoing accretion processes of, for instance, cold clumps and subhalos, and energy diffusion related to the feedback). The distributions of the temperature bins are well described by a log-normal function (see Fig. 5). Once these inhomogeneities are removed, the reconstructed temperature profiles vary, with the intensity of the variation increasing with the radius (see Fig. 7). That has indeed the effect of changing the absolute temperature and gradient of the profiles with an obvious impact on the estimate of the fundamental integrated properties (e.g., the total mass).

The overall level of turbulence, resolved in our spectral analysis, suggests that, on average, there is a mix of isobaric and adiabatic fluctuations (see Fig. 9). We constrained the average 3D Mach number in our sample to ${\mathcal{M}}_{3\mathrm{D}}=0.{36}_{-0.09}^{+0.16}$. The global effect translates into a level of energy in turbulence that is about 7% of the thermal energy (see Fig. 13), implying an estimated hydrostatic bias of approximately 6% (ranging between 0 and 29%) in this representative subsample of CHEX-MATE clusters. However, when we apply the correction to E_turb by integrating the power spectrum between a dissipation scale and an injection scale, we obtain a median hydrostatic bias of b ∼ 11%, covering a range between 0 and 37%.

Once resolved as a function of the radius, the estimates of σ_Tj, int/T match the ones from state-of-art hydrosimulations (see Fig. 12). Also, the observed σ_Tint − T relation is slightly steeper than the one obtained with hydrodynamical simulations including thermal conduction (which drastically smooths temperature variations and homogenizes the ICM; see Fig. 10) at 1/3 of the Spitzer value, but it is significantly shallower than the one derived without thermal conduction.

The knowledge acquired by this work on the level of fluctuations present in the gas temperature maps, together with the evaluation of the dynamical state obtained by the estimate of the morphological parameters determined from the analysis of the X-ray images (Campitiello et al. 2022), will allow us both to guide the interpretation on, and to quantify, the amount of scatter that will be resolved both in the thermodynamic radial profiles and in the integrated quantities of the CHEX-MATE clusters. Given the encouraging results presented in this work, the next step will be to extend this kind of analysis to the entire CHEX-MATE sample and to investigate the link between the observed gas T, gas density, and surface brightness fluctuations to the scatter characterizing the distribution of the integrated quantities in the scaling relations (e.g., Pratt et al. 2009; Lovisari et al. 2020b). An investigation of the statistics of the fluctuations in the surface brightness maps of the CHEX-MATE objects, and relating the density fluctuations and the turbulent motions, will be presented in a forthcoming paper.

Acknowledgments

We thank the anonymous referee for her/his review of the manuscript. We also thank M. Roncarelli for useful discussion. L.L. and S.E. acknowledge financial contribution from the contracts ASI-INAF Athena 2019-27-HH.0, “Attività di Studio per la comunità scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASI-INAF n. 2017-14-H.0), and from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158). L.L. also acknowledges financial contribution from the INAF grant 1.05.12.04.01. W.F. and C.J. acknowledge support from the Smithsonian Institution, the Chandra High Resolution Camera Project through NASA contract NAS8-03060, and NASA Grants 80NSSC19K0116, GO1-22132X, and GO9-20109X. B.J.M. acknowledges support from STFC grant ST/V000454/1. G.W.P. acknowledges support from CNES, the French space agency. J.S. and J.K. were supported by NASA Astrophysics Data Analysis Program (ADAP) Grant 80NSSC21K1571. This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project #565 (Multi-Wavelength Studies of the Culmination of Structure Formation in the Universe). In this work we made use of the WVT binning algorithm by Diehl & Statler (2006), which is a generalization of Cappellari & Copin (2003) Voronoi binning algorithm. We also made use of the following packages: FTOOLS (Blackburn et al. 1995), DS9 (Joye et al. 2003), Astropy (Astropy Collaboration 2013, 2018).

References

Appendix A: Gallery

In Fig. A.1 we show the recovered surface brightness, electron density, and temperature maps within R₅₀₀ for all the clusters analyzed in this work.

Fig. A.1. From left to right: (a) bkg-subctracted and exposure-corrected images in the 0.3-7 keV band; The size of the boxes corresponds to R500; (b) binned surface brightness maps; (c) projected electron density maps (d) projected temperature maps; (e) relative temperature error maps. The white circles in the voronoi maps corresponds to 0.15R500, 0.5R500, and R500.

Fig. A.2. From left to right: (a) bkg-subctracted and exposure-corrected images in the 0.3-7 keV band; The size of the boxes corresponds to R500; (b) binned surface brightness maps; (c) projected electron density maps (d) projected temperature maps; (e) relative temperature error maps. The white circles in the voronoi maps corresponds to 0.15R500, 0.5R500, and R500.

Fig. A.3. From left to right: (a) bkg-subctracted and exposure-corrected images in the 0.3-7 keV band; The size of the boxes corresponds to R500; (b) binned surface brightness maps; (c) projected electron density maps (d) projected temperature maps; (e) relative temperature error maps. The white circles in the voronoi maps corresponds to 0.15R500, 0.5R500, and R500.

Fig. A.4. From left to right: (a) bkg-subctracted and exposure-corrected images in the 0.3-7 keV band; The size of the boxes corresponds to R500; (b) binned surface brightness maps; (c) projected electron density maps (d) projected temperature maps; (e) relative temperature error maps. The white circles in the voronoi maps corresponds to 0.15R500, 0.5R500, and R500.

Appendix B: Comparison with the curvelet maps

Thermodynamical 2D maps have been extensively used in the study of galaxy clusters, thanks to their great potential to characterize the dynamical state of a system. For this reason, a number of spectral-mapping methods have been developed. Based on the approach, these methods can be divided into three main categories: (i) hardness ratio method which combines X-ray imaging in several energy bands (e.g., Ferrari et al. 2006; ii) wavelet (see Bourdin et al. 2004; Bourdin & Mazzotta 2008) or curvelet (see Bourdin et al. 2015) analysis using the wavelet or curvelet coefficients to couple a multi-scale spectroscopic analysis with a structure detection scheme; (iii) spatially resolved spectroscopy of independent cells sampled by a required S/N or a minimum number of counts for which at least three different binning strategies have been proposed in this context: adaptive binning (e.g., Sanders & Fabian 2001, Lovisari et al. 2009, 2011, O’Sullivan et al. 2011); contour binning (Sanders 2006, Hofmann et al. 2016); Weighted Voronoi Tesselation (WVT) method (e.g., Cappellari & Copin 2003). In general, the Voronoi algorithms define the smallest possible regions for the spectral analysis, given the statistics of the observation. It is based on the distribution of the signal and noise in a given band and it puts together pixels to reach a required S/N, regardless of the “temperature” that would be measured in that cell. The wavelet or curvelet analysis has the advantage of being not driven by the distribution of the brightness of the emission. Conversely, the advantage of Voronoi algorithms is that (i) regions can be treated as independent from each other (when neglecting the effect of the PSF) and (ii) it does not try to follow the gradients of the ICM emission.

For ten objects in our sample, we compared the temperature map obtained through the spectral analyses of the Voronoi tessellation regions with the one obtained through the curvelet analysis. For each Voronoi region, we computed the weighted temperature from the curvelet map and we then compared the two T values. In Fig. B.1, we show a summary plot to demonstrate that the temperatures recovered between the two methods are in good agreement, with an average value of the ratio (T_this work − T_cur)/ϵ of 0.005 (median: 0.06; 1st-3rd quartile: -0.42, 0.53).

Fig. B.1. Comparison of the temperatures measured in the maps reconstructed in this work and the ones obtained from the curvelet technique. The histogram represents the distribution of the difference in spectroscopic measurements normalized to the statistical error $ϵ=\sqrt{{ϵ}_{\mathrm{this}\mathrm{work}}^2+{ϵ}_{\mathrm{cur}}^2}$ for 1,916 spatial cells obtained in 10 objects. The vertical dashed lines indicate 0 (red) and -3, 3 (green).

We also compared the radial temperature fluctuations estimated from the Voronoi and curvelet temperature maps. We found a remarkably good agreement between the σ_T/T recovered with the two methods in the regions between 0.2-0.8R₅₀₀. At large radii (i.e., > 0.8R₅₀₀) the comparison is not straightforward because of the coarse mapping with the Voronoi method and of the masking of low signal-to-noise regions in the curvelet analysis. Within 0.2R₅₀₀, the σ_T/T values in the Voronoi maps are slightly larger than what was found by the curvelet analysis. A detailed investigation is required to understand the cause of this slight difference but it is beyond the scope of the paper. However, this does not affect any of the main conclusions.

Appendix C: Temperature cell weights

In Fig. C.1, we show how well the overall cluster temperatures recovered from the maps (i.e., T_1D, 500) match the fitted values using one single extraction region (i.e., T_spec). We used six different weightings, each of them including a proxy of the X-ray emissivity (i.e., either the value of the surface brightness of the cells or the electron density obtained using Eq. 1) multiplied by: (i) the area of the cell, (ii) the area of the cell multiplied by T^α (with α=-0.75 from Mazzotta et al. 2004), and (iii) the area multiplied by the relative errors of the cells. Since the relative errors are smaller for lower temperatures, both the latter approaches weight more the cells with lower temperatures.

Fig. C.1. Comparison between the temperature derived from the maps (i.e., T1D, 500) and the global temperature (i.e., Tspec) obtained by fitting a single spectrum extracted within R500. In each subplot the text in magenta indicates the weighting used to recover T1D, 500 from the maps, while the values of μ and σ indicate the median of the distribution and its dispersion. In blue we show the results using the maps obtained with S/N=30 while in red we show the ones obtained with S/N=50.

In general, using SB instead of n_e as a proxy for the X-ray emissivity returns T_1D, 500 in better agreement with T_spec. This is possible because there is no assumption on the scale of the line of sight (i.e., the volume to be used to recover the electron density). If so, in principle this result could be used to provide a constraint on the elongation along the line of sight (helpful for a measure of the triaxiality), but this is beyond the scope of this paper.

Using only the area of the cells as weights is clearly not sufficient (see top panels of Fig. C.1) because it does not account for the energy dependence of the effective area. Using the relative errors of each cell in the weighting (bottom panels) returns T_1D, 500 in much better agreement with T_spec (e.g., compare the μ and σ values of the different subplots). Weighting each cell by T^α works also relatively well (middle panels), although there is an increasing deviation at lower temperatures. However, since for comparison with the simulations using the factor T^α is the best option (e.g., Mazzotta et al. 2004), and after verifying that the impact on the results of the paper is not significant, we used this latter weighting as our default method.

Appendix D: Voronoi tessellation and map resolution

Given a certain image, the binning mask obtained with WVT is deterministic (i.e., by running the algorithm multiple times, one obtains always the same result). Depending of the resolution of the mask, it is possible that some inhomogeneities are missed because they are spanning multiple cells, with the net result that the feature is washed out. In order to investigate how the choice of the map center impacts the obtained map and our results, we derived maps centered offset from the X-ray peak. The center of the new maps was computed from a grid with the X values obtained as X_center=X_peak+C×R₅₀₀ and the Y values as Y_center=Y_peak+C×R₅₀₀, where C can assume 3 values: -0.1, 0, +0.1. As an example, in the top panel of Fig. D.1, we show the results for G041.45+29.10. Despite the different binning, the same features are observed in all the maps. In the bottom panel we show that the recovered profiles from the individual maps differ only by a few percent (and are always consistent within the uncertainties) as far as the map resolution is good enough (typically within ∼0.6R₅₀₀). At large radii, when the binning becomes quite coarse there is an increase of the scatter, although the data points are still consistent within the 1σ errors. Similar results are obtained with the other clusters in the sample, so that the results of the paper are not affected by the deterministic nature of the Voronoi binning.

Fig. D.1. Impact of the Voronoi tessellation on the 2D temperature distribution. Top: temperature maps obtained for G041.45+29.10 by offsetting the center of the map wrt to the X-ray peak. The green circles identify R500 centered on the peak. Bottom: ratio between the temperature recovered from the nine maps (using Eq. 2) and the temperature estimated using the spectra extracted from azimuthal annuli. Each color (from blue to red) represents one of the maps in the top panel (from top-left to bottom-right).

The results presented in the current analysis are based on the binning of temperature maps obtained with S/N=30. Such a choice is set by the requirement to have typical temperature uncertainties of the order of 10-20%. However, the resolution of the maps can hide some relevant features and reduce or increase the observed scatter, in particular in the electron density maps used to infer the type of fluctuations. This is probably due to the increasing size of the binning in the outer regions, amplified by the presence of steep surface brightness gradients. To verify that, we produced SB maps with S/N=15, 30, and 50 and we derived the scatter as a function of radius. In Fig. D.2 we show the results with the electron densities calculated under the assumption of $n_e=\sqrt{\mathit{SB}}$. Indeed, the ratio σ_ne/n_e increases from low to high S/N. However, it is worth noticing that the trend remains the same with an increase as a function of the radius (as discussed in Sect. 3.4 and shown in Fig. 9), although there is a flattening beyond 0.3-0.4R₅₀₀.

Fig. D.2. Impact of the map resolution on the density profile in the case of G041.45+29.10. We truncate the plot at 0.6R500 because beyond that radius we have a poor resolution (in particular with S/N=50). Each point represents the value of σne/ne computed at the radius of the cells in the map.

Appendix E: Projection effects to the density and temperature perturbations

By analyzing a set of 3D cosmological simulations (Vazza et al. 2017; Simonte et al. 2022) we tested how well Eq. 11 (derived following Schuecker et al. 2004) provides a good representation of the correlation between the projected density and temperature fluctuations (i.e., the relation of the 2D perturbations being a factor of ∼2 flatter than the relation between 3D density and temperature perturbations). The simulations are performed with the code ENZO (Bryan et al. 2014) and accurately capture the dynamical evolution of the ICM as galaxy clusters undergo matter accretion throughout their evolutionary process. Although the employed simulations are classified as non-radiative, previous studies have indicated that the influence of non-gravitational effects on the larger scales of interest (> few tenths of kpc) is comparatively limited, in contrast to the significant impact of mergers and accretion phenomena in shaping the properties of ICM (e.g., Vazza et al. 2012; Valdarnini 2019). We refer the reader to Simonte et al. (2022, and references therein) for further details. The 3D ICM perturbations are estimated in spherical shells with a size of 0.1R₅₀₀. In Fig. E.1 we show as red points the fitted slope (i.e., γ-1) between the 3D density and temperature perturbations and as blue diamonds the 2D slopes from three different projections. It is clear that projection effects tend to flatten the relation as expected. The flattening is roughly a factor of 2 as derived following Schuecker et al. (2004; see the green squares).

Fig. E.1. For each of the 8 simulated clusters, we show the mean ratio between the radial values of δT3D/T3D and δn/n (red points), δT2D/T2D and $\delta n_{2D}^2/n_{2D}^2$ in three different projections (blue diamonds), and the latter values corrected by the factor of 2 from Schuecker et al. (2004; the green squares represent the average values while the error bars represent the minimum and maximum values from the different projections).

All Tables

All Figures

	Fig. 1. Distribution of the clusters in our sample in the concentration-centroid-shift (i.e., c − w) space. The values are taken from Campitiello et al. (2022). The dashed lines indicate the median values. The most relaxed clusters are in the bottom-right corner while the most disturbed in the upper-left corner.
In the text

	Fig. 2. Comparison between the profile derived with the direct spectral fitting of each annulus j (i.e., Tspec, j). and the temperatures recovered from the maps using Eq. (2) (i.e., T1D, j) in the same spatial regions. Each blue point corresponds to an annulus of the 28 clusters, while the black points represents the average values with the shadow area including the 16th and 84th percentiles.
In the text

	Fig. 3. Example of clipping for the cluster G041.45+29.10. In the top panel we show the distribution of si values with the regions in red and blue being the cells deviating more than 1σ from the azimuthal value. In the bottom panel, we show their distribution in the temperature map derived within R500.
In the text

	Fig. 4. Radius of the cells for the maps obtained with S/N = 30, in fraction of R500, as a function of the distance from the center. The feature observed at r/R500 ⪆ 0.65 is the result of the large cell sizes in the outer regions.
In the text

Fig. 5. Observed properties of the 2D temperature distribution. Left: Measured values (blue circles) of the skewness in the temperature distribution of T2D, i/T2D, j (being T2D, j the temperature of the annulus encompassing the cell i) compared with a typical uncertainty (green squares) that can be associated with the number of cells on the temperature map. Middle: QQ plot, assuming a normal distribution, of the T2D, i/T1D, j values for the maps with S/N = 30. Right: QQ plot assuming a log-normal distribution. In the inset plots, we show the results when using the maps with S/N = 50. The units of the QQ plots are given in z-score (i.e., subtracting the mean from each data point and dividing it by the standard deviation).

In the text

Fig. 6. Spectral parameters derived before (left panels) and after (right panels) clipping for the regions deviating more than 1σ from the azimuthal values. In blue and red we show the results obtained with maps of S/N = 30 and S/N = 50, respectively. In each panel we provide the correlation (i.e., Spearman rank test value r) between the spectral parameter and T, the median value of the distribution μ, and the standard deviation σ.

In the text

	Fig. 7. Ratio between T1D, j, c (i.e., the temperature estimated in each annulus j after clipping the most deviating regions) and T1D, j (i.e., the temperature estimated using all the temperature cells in each annulus) as a function of radius. The black points represents the average values, while the shadow area includes the 16th and 84th percentiles.
In the text

Fig. 8. Impact of the temperature fluctuations to the cluster overall temperature and azimuthal profile. Top: Comparison between the global temperature estimated using all the cells in the 2D maps and the temperature estimated after excluding the cells deviating more than 1σ with respect to the azimuthal value. In the inset plots we show the absolute difference between T500 and T500, c (bottom-right inset) and its fractional variation (top-left inset). Bottom: Change in the gradient of the temperature profile fitted in the region 0.15–0.75R500 before and after the clipping. In the inset plot, we show its fractional variation.

In the text

Fig. 9. Slope values of the intrinsic scatter of the temperature (left panel) and electron density (middle panel) measured within ∼0.6R500. In the right panel, we show the polytropic index γ (see Eq. (10)) derived for each cluster. Blue points are from maps with S/N = 30 and red points are from maps with S/N = 50. The solid lines represent the median value μ (also reported in red on the top-left corners) while the shaded area includes the 16th and 84th percentiles.

In the text

Fig. 10. Total temperature dispersion within R500 as a function of the overall cluster temperature. Left: Relation between the mean cluster temperature T500 and the intrinsic scatter σTint within R500. The dotted orange line represents the best-fit relation to the measured points and is compared to the observations by Frank et al. (2013; black) and simulations with (green) and without (cyan) thermal conduction by Rasia et al. (2014). Right: Relative scatter as function of the cluster temperature. In blue, we show the observed values (i.e., total scatter σTtot/T500), in green we plot the relative intrinsic scatter (i.e., σTint/T500) computed as given in Eq. (13) where the statistical errors (i.e., ϵT1D, j, stat) are shown in pink and are computed using Eq. (4), in magenta we show the scatter associated with the underlying temperature profile (i.e., σTprof/T500), and in red the average statistical uncertainty of the T2D, i measurements (i.e., ϵT2D, i, stat).

In the text

	Fig. 11. Correlation between Gaussian standard deviations ςT and ςne measured in different radial shells. Each data point corresponds to the values from a single shell of a cluster. The orange line represents the best-fit relation of the data points and is compared with the best-fits obtained by Rasia et al. (2014).
In the text

Fig. 12. Comparison of σTj, int/T estimated in cells of 0.1R500 as a function of radius between the CHEX-MATE objects (blue; the 16th and 84th percentiles were computed excluding the annuli for which the statistical errors prevents a proper determination of the intrinsic scatter) and results extracted from various hydrodynamical simulations (see text for details in Sect. 4.1). For cosmological simulations (i.e., The 300 GX and GIZMO), we represent the distribution for 45 objects with the median and a dispersion equal to the inter-quartile range divided by 1.35. For the FLASH simulations, we present the results obtained for 3 different projections (with different line styles) of the same object simulated with two different levels of turbulence (i.e., different Mach).

In the text

Fig. 13. Mass bias and ratio between turbulent and thermal energy in the clusters in our sample. Top: Ratio between the turbulent and the thermal energy computed based on the σT, int/T values and under the assumption that the isobaric fluctuations are dominant. The dashed horizontal lines are computed for different Mach numbers using the following equation: $E_{\mathrm{turb}}/E_{\mathrm{th}}=0.5\gamma (\gamma -1){\mathcal{M}}_{3D}^2$. In the right y-axis, we provide the mass bias b estimated within R500 through Eq. (14). The empty squares indicate the values of Eturb/Eth (and b) after the correction accounting for the underlying power spectrum (see the text for more details). Bottom: Distribution of the values of the hydrostatic bias. In gray (blue), we show the b values before (after) accounting for the integration of the power spectrum between some characteristic scales (see text in Sect. 4.3). Dotted lines indicate the means of the distributions.

In the text

	Fig. A.1. From left to right: (a) bkg-subctracted and exposure-corrected images in the 0.3-7 keV band; The size of the boxes corresponds to R500; (b) binned surface brightness maps; (c) projected electron density maps (d) projected temperature maps; (e) relative temperature error maps. The white circles in the voronoi maps corresponds to 0.15R500, 0.5R500, and R500.
In the text

	Fig. A.2. From left to right: (a) bkg-subctracted and exposure-corrected images in the 0.3-7 keV band; The size of the boxes corresponds to R500; (b) binned surface brightness maps; (c) projected electron density maps (d) projected temperature maps; (e) relative temperature error maps. The white circles in the voronoi maps corresponds to 0.15R500, 0.5R500, and R500.
In the text

	Fig. A.3. From left to right: (a) bkg-subctracted and exposure-corrected images in the 0.3-7 keV band; The size of the boxes corresponds to R500; (b) binned surface brightness maps; (c) projected electron density maps (d) projected temperature maps; (e) relative temperature error maps. The white circles in the voronoi maps corresponds to 0.15R500, 0.5R500, and R500.
In the text

	Fig. A.4. From left to right: (a) bkg-subctracted and exposure-corrected images in the 0.3-7 keV band; The size of the boxes corresponds to R500; (b) binned surface brightness maps; (c) projected electron density maps (d) projected temperature maps; (e) relative temperature error maps. The white circles in the voronoi maps corresponds to 0.15R500, 0.5R500, and R500.
In the text

Fig. B.1. Comparison of the temperatures measured in the maps reconstructed in this work and the ones obtained from the curvelet technique. The histogram represents the distribution of the difference in spectroscopic measurements normalized to the statistical error $ϵ=\sqrt{{ϵ}_{\mathrm{this}\mathrm{work}}^2+{ϵ}_{\mathrm{cur}}^2}$ for 1,916 spatial cells obtained in 10 objects. The vertical dashed lines indicate 0 (red) and -3, 3 (green).

In the text

Fig. C.1. Comparison between the temperature derived from the maps (i.e., T1D, 500) and the global temperature (i.e., Tspec) obtained by fitting a single spectrum extracted within R500. In each subplot the text in magenta indicates the weighting used to recover T1D, 500 from the maps, while the values of μ and σ indicate the median of the distribution and its dispersion. In blue we show the results using the maps obtained with S/N=30 while in red we show the ones obtained with S/N=50.

In the text

Fig. D.1. Impact of the Voronoi tessellation on the 2D temperature distribution. Top: temperature maps obtained for G041.45+29.10 by offsetting the center of the map wrt to the X-ray peak. The green circles identify R500 centered on the peak. Bottom: ratio between the temperature recovered from the nine maps (using Eq. 2) and the temperature estimated using the spectra extracted from azimuthal annuli. Each color (from blue to red) represents one of the maps in the top panel (from top-left to bottom-right).

In the text

	Fig. D.2. Impact of the map resolution on the density profile in the case of G041.45+29.10. We truncate the plot at 0.6R500 because beyond that radius we have a poor resolution (in particular with S/N=50). Each point represents the value of σne/ne computed at the radius of the cells in the map.
In the text

Fig. E.1. For each of the 8 simulated clusters, we show the mean ratio between the radial values of δT3D/T3D and δn/n (red points), δT2D/T2D and $\delta n_{2D}^2/n_{2D}^2$ in three different projections (blue diamonds), and the latter values corrected by the factor of 2 from Schuecker et al. (2004; the green squares represent the average values while the error bars represent the minimum and maximum values from the different projections).

In the text