CHEX-MATE: Characterization of the intra-cluster medium temperature distribution

We study the perturbations in the temperature (and density) distribution for 28 clusters selected from the CHEX-MATE sample to evaluate and characterize the level of inhomogeneities and the related dynamical state of the ICM. We use these spatially resolved 2D distributions to measure the global and radial scatter and identify the regions that deviate the most from the average distribution. During this process, we introduce three dynamical state estimators and produce clean temperature profiles after removing the most deviant regions. We find that the temperature distribution of most of the clusters is skewed towards high temperatures and is well described by a log-normal function. There is no indication that the number of regions deviating more than 1$\sigma$ from the azimuthal value is correlated with the dynamical state inferred from morphological estimators. The removal of these regions leads to local temperature variations up to 10-20% and an average increase of $\sim$5% in the overall cluster temperatures. The measured relative intrinsic scatter within $R_{500}$, $\sigma_{T,int}/T$, has values of 0.17$^{+0.08}_{-0.05}$, and is almost independent of the cluster mass and dynamical state. Comparing the scatter of temperature and density profiles to hydrodynamic simulations, we constrain the average Mach number regime of the sample to $M_{3D}$=0.36$^{+0.16}_{-0.09}$. We infer the ratio between the energy in turbulence and the thermal energy, and translate this ratio in terms of a predicted hydrostatic mass bias $b$, estimating an average value of $b\sim$0.11 (covering a range between 0 and 0.37) within $R_{500}$. This study provides detailed temperature fluctuation measurements for 28 CHEX-MATE clusters which can be used to study turbulence, derive the mass bias, and make predictions on the scaling relation properties.


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 7 -10 8 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. 2012Vazza et al. , 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. 2015Liu et al. , 2016;;Hitomi Collaboration et al. 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., M<0.5) tend to be in pressure equilibrium with the medium, while stronger turbulence (i.e., 0.5<M<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;Gaspari et al. 2014;Mohapatra et al. 2020Mohapatra et al. , 2021;;Simonte et al. 2022;Zhuravleva et al. 2023).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. 2016Zhuravleva et al. , 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 M regime, the fluctuations in temperature are as clear and robust tracers as density fluctuations, and even better than pressure fluctuations.With higher M numbers (e.g., M>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 Notes.When this value is associated with a single cell on the maps we add the subscript i to these symbols, when it is associated to an annulus we add the subscript j, while if the subscript is omitted, the value refers to the global property (i.e., within R 500 ).
temperature distribution are, indeed, the main sources of systematic bias in HE mass estimates from X-ray data (see, e.g., Rasia et al. 2006Rasia et al. , 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 Sect. 4 and 5, respectively.
Throughout the paper, we assume a flat ΛCDM cosmology with Ω m = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 .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.

Sample
The Cluster HEritage project with XMM-Newton -Mass Assembly and Thermodynamics at the Endpoint of structure formation (CHEX-MATE 1 ; CHEX-MATE Collaboration et al. 2021) is a ground-breaking multi-wavelength investigation of 61 local clusters (Tier-1: 0.05 < z < 0.2; 2×10 14 M ⊙ < M 500 < 9×10 14 M ⊙ ) and 61 massive clusters (Tier-2: z<0.6 with M 500 > 7.25 × 10 14 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 500 , for this minimally biased, signal-tonoise-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 500 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.

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 LogN-LogS distribution.More details can be found in Bartalucci et al. (2023).
Article number, page 3 of 24 Notes.Cluster indices and names, X-ray peak coordinates, and redshifts are shown in columns 1-4.In columns 5-6 we provide the list of XMM-Newton observations that we investigated, and the corresponding clean exposure times (in ks) for MOS1, MOS2, and pn.The value of R 500 (in kpc) used to extract the maps provided in column 7 are taken from Planck Collaboration et al. (2016).In columns 8-11, we provide the concentration, centroid-shift, power-ratio, and M all parameters estimated by Campitiello et al. (2022).
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.

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 tool2 (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 powerlaw (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 Bour-Article number, page 4 of 24 din et al. (2023) 3 .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.

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 where N is the APEC normalization, D a the angular diameter distance, f n e /n H 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 de- , 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.

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 where i identifies the cell number in the 2D distributions, and w i are the weights that take into account both the area A i j covered by the cell i in a given annulus j and the emissivity (more details about the weighting are provided in Appendix C) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 r/R 500 0.8 0.9 Comparison between the profile derived with the direct spectral fitting of each annulus j (i.e., T spec, j ). and the temperatures recovered from the maps using Eq. 2 (i.e., T 1D, 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.
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.
We also associate to each measurement of T 1D, j an uncertainty ϵ T 1D, j from the propagation of the error ϵ T 2D,i (we assume symmetric errors) on the spectral measurement T 2D,i and a standard deviation σ T j (i.e., the total scatter or dispersion around the average value) defined as The measured temperature dispersion in a given annulus is a combination of the intrinsic variation (σ T j,int ) of the ICM temperature distribution (due to turbulent motions and mergers), of the statistical uncertainty ϵ T 1D, j,stat on T 1D, j (see Eq. 4), and of the statistical uncertainty ϵ T i,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: The gas density radial profile is recovered from the geometrical deprojection of the normalization of the spectral model (see Article number, page 5 of 24 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 From Eq. 7, we can then write 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 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.

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 500 ) 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.

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 500 ) 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 500 .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 Article number, page 6 of 24   4 3 2 1 0 1 2 3 4 Theoretical Quantiles Fig. 5. Observed properties of the 2D temperature distribution.Left: measured values (blue circles) of the skewness in the temperature distribution of T 2D,i /T 1D, j (being T 1D, 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 T 2D,i /T 1D, 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).
Table 3. Parameters determined from the temperature maps with S/N=30.

Name
N cells T 1D,500 Notes.Cluster names are shown in column 1.In column 2, we provide the number of temperature cells for each cluster.In column 3, the overall temperature in keV recovered using Eq. 2 in the annulus 0-R 500 .In columns 4-7, we provide the estimated spectroscopic parameters.Column 8 indicates the number of cells left after removing the ones deviating more than 1σ.The corresponding temperature and spectroscopic parameters are quoted in columns 9-13.
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 Article number, page 7 of 24 A&A proofs: manuscript no.2Dmaps Notes.The coefficient r is the non-parametric statistical measure used to study how well the relationship between two variables can be described using a monotonic function.The p-value for a hypothesis test whose null hypothesis is that two sets of data are linearly uncorrelated.
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 hightemperature values.The skewness values of the observed distribution do not change significantly even if we remove the core regions (i.e., < 0.15R 500 ), 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).
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-degree 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).

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 20 (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 |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 500 ).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.5R500 ), 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.5T500 ), 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 σ T i,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 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 σ.
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, σ T i,int /T j goes to zero as expected.The average values of |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 |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).

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 Fig. 7. Ratio between T 1D, j,c (i.e., the temperature estimated in each annulus j after clipping the most deviating regions) and T 1D, 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.
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.4R500 , 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.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 500 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 powerlaw in the region 0.15-0.75R500 (beyond 0.75R 500 the map is too coarse).There is a very mild correlation (i.e., r=-0.24) between the change in slope and in temperature.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 T 500 and T 500,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.

Properties of the fluctuations
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. 2014Zhuravleva et al. , 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-weighted4 average along the line of sight of the threedimensional temperature T (see, e.g., Mazzotta et al. 2004).For a polytropic gas, fluctuations in the latter quantity, δT i = T i − T , should satisfy the relation 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 where w ≈ n 2 T a , n 2 2D = dl n 2 , and δn 2 2D = 2n 2D δn 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 Eq. 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: 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 500 up to a radius where the S/N=30 maps have enough resolution (typically 0.6R 500 ; 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.
The fluctuations in density have been estimated similarly to the ones in temperature (see Eq. 5 and Eq.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.

Global intrinsic scatter
In the previous section, we discussed the radial scatter (i.e., the dispersion in different annuli with size 0.1R 500 ).Here, we discuss the total dispersion within R 500 (i.e., the dispersion of all the cells around the temperature of the cluster, T 500 ).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 where σ T tot and ϵ T 1D, j,stat are computed using Eq. 5 and Eq. 4 replacing T 2D, j with the overall cluster temperature T 500 , and ϵ T 2D,i,stat is the average statistical uncertainty of the T 2D,i measurements.For each object with a given M 500 and redshift (see Table 2), we recovered the expected projected temperature profile, σ T,pro f , 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 σ T int -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).

Discussion
Gas inhomogeneities are ubiquitous in the ICM and affect various scales (see, e.g., Nagai & Lau 2011;Roncarelli et al. 2013; Article number, page 11 of 24 A&A proofs: manuscript no.2Dmaps  T1D,j,stat T2D,i,stat Fig. 10.Total temperature dispersion within R 500 as a function of the overall cluster temperature.Left: relation between the mean cluster temperature T 500 and the intrinsic scatter σ T int within R 500 .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 /T 500 ), in green we plot the relative intrinsic scatter (i.e., σ T int /T 500 ) computed as given in Eq. 13 where the statistical errors (i.e, ϵ T 1D, 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., σ T pro f /T 500 ), and in red the average statistical uncertainty of the T 2D,i measurements (i.e, ϵ T 2D,i,stat ).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 "physicallymotivated" 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., |s i |, std(s i ), σ T i,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 mor-phological 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 500 (although most of the cells are within 0.5R 500 ), 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.

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 500 .We refer to ς ne and ς T as the standard deviations of Gaussian distributions of the logarithmic values.Apart from some differences in smoothed-particlehydrodynamics (SPH) and adaptive-mesh-refinement (AMR) codes (see discussion in Rasia et al. 2014  physics prescription, including non-radiative / NR, cooling-starformation 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.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 σ T j,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).
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., M<0.5) perturbations are mainly isobaric and therefore one expects σ T /T ∼ σ n e /n e , while for high Mach numbers (i.e., 0.5<M<1) the perturbations shift to the adiabatic regime implying (e.g. for γ=5/3) σ T /T ∼ 0.67σ n e /n e .The fluctuations in the ICM move from isobaric to adiabatic from the M 3D =0.25 to M 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 mergerlike 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.

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.

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).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 Eq.: E turb /E th = 0.5 γ (γ − 1) M 2 3D .In the right y-axis, we provide the mass bias b estimated within R 500 through Eq. 14.The empty squares indicate the values of E turb /E th (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 adiabatic regime, we have σ T /T ∼ 0.67σ n e /n e , while in the isobaric regime σ T /T ∼ σ n e /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σ n e /n e .We can now estimate the Mach number by applying the relation M 1D ∼ σ n e /n e ∼ 1.2 * σ T /T which gives M 3D = √ (3) * M 1D ∼ 2.1 * σ T /T .We measured σ T,int /T =0.17 +0.08  −0.05 (see Sect. 3.5 and Fig. 10).Therefore, M 3D = 0.36 +0.16  −0.09 .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 turb /E th = 0.5 γ (γ − 1) M 2 3D that translates into a M 3D from simulations in the range 0.30-0.95,consistent with our observational con-Article number, page 14 of 24 straints.Our result is also in agreement with the finding of Hofmann et al. (2016) who found M 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, M numbers would appear supersonic (mimicking shocks) in the majority of systems because of the substantial boost in the density maps.
The measured M 3D corresponds to E turb /E th ∼ 0.07 In Fig. 13, we show the estimated ratios between the turbulent and thermal energy within R 500 which points, for the studies sample, to a mass bias of b=0.06 +0.07 −0.03 , 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).
The variance in the maps is related to the turbulent energy E turb , which is proportional to the integral of the power spectrum with the energy spectrum E(k) = k α+2 .Eq. 15 is fully generic and can be applied to any mode of fluctuation.For Kolmogorovlike 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 500 (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 turb /E th ∼ 0.12 +0.16 −0.08 and b ∼0.11 +0.11   −0.07 (see Fig. 13) which tend to agree better to what derived from Xray 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.

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 M 3D =0.36 +0.16  −0.09 .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 σ T j,int /T match the ones from state-of-art hydrosimulations (see Fig. 12).Also, the observed σ T int -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.Article number, page 24 of 24

1Fig. 1 .
Fig. 1.Distribution of the clusters in our sample in the concentrationcentroid-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.

Fig. 3 .
Fig. 3. Example of clipping for the cluster G041.45+29.10.In the top panel we show the distribution of s i 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 R 500 .

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

TFig. 8 .
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 T 500 and T 500,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.

Fig. 9 .
Fig. 9. Slope values of the intrinsic scatter of the temperature (left panel) and electron density (middle panel) measured within ∼0.6R 500 .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.

Fig. 11 .
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 byRasia et al. (2014).

Fig
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 ϵ = ϵ 2 this work + ϵ 2 cur for 1,916 spatial cells obtained in 10 objects.The vertical dashed lines indicate 0 (red) and -3, 3 (green).
Fig. C.1.Comparison between the temperature derived from the maps (i.e., T 1D,500 ) and the global temperature (i.e., T spec ) obtained by fitting a single spectrum extracted within R 500 .In each subplot the text in magenta indicates the weighting used to recover T 1D,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.

Table 4 .
Spearman correlation coefficient r and p − value between the spectral and morphological parameters.
for all the ICM Comparison of σ T j,int /T estimated in cells of 0.1R 500 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., The300 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).