# CHEX-MATE: Robust reconstruction of temperature profiles in galaxy clusters with XMM-Newton

M Rossetti, D Eckert, F Gastaldello, E Rasia, G.W Pratt, S Ettori, S Molendi, M Arnaud, M Balboni, I Bartalucci, et al.

► To cite this version:

M Rossetti, D Eckert, F Gastaldello, E Rasia, G.W Pratt, et al.. CHEX-MATE: Robust reconstruction of temperature profiles in galaxy clusters with XMM-Newton. Astronomy & Astrophysics - A&A, 2024, 686, pp.A68. 10.1051/0004-6361/202348853. hal-04523136

HAL Id: hal-04523136

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

Submitted on 26 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
CHEX-MATE: Robust reconstruction of temperature profiles in
galaxy clusters with XMM-Newton*

M. Rossetti¹, D. Eckert², F. Gastaldello¹, E. Rasia³,⁴, G. W. Pratt⁵, S. Ettori⁶,⁷, M. Arnaud⁵, M. Balboni¹,⁸, I. Bartalucci¹, R. M. Batalha⁵, S. Borgani³,⁴,⁹, H. Bourdin¹⁰,¹¹, S. De Grandi¹², F. De Luca¹⁰,¹¹, M. De Petris¹³, W. Forman¹⁴, M. Gaspari¹⁵, S. Ghizzardi¹, A. Iqbal⁵, S. Kay¹⁶, L. Lovisari¹,¹⁴, B. J. Maughan¹⁷, P. Mazzotta¹⁰,¹¹, E. Pointecouteau¹⁸, G. Riva¹,¹⁹, J. Sayers²⁰, and M. Sereno⁶,⁷

(Affiliations can be found after the references)

Received 5 December 2023 / Accepted 26 February 2024

ABSTRACT

The “Cluster HEritage project with *XMM-Newton*: Mass Assembly and Thermodynamics at the End point of structure formation” (CHEX-MATE) is a multi-year heritage program to obtain homogeneous *XMM-Newton* observations of a representative sample of 118 galaxy clusters. The observations are tuned to reconstruct the distribution of the main thermodynamic quantities of the intra-cluster medium up to $R_{500}$ and to obtain individual mass measurements, via the hydrostatic-equilibrium equation, with a precision of 15–20%. Temperature profiles are a necessary ingredient for the scientific goals of the project and it is thus crucial to derive the best possible temperature measurements from our data. This is why we have built a new pipeline for spectral extraction and analysis of *XMM-Newton* data, based on a new physically motivated background model and on a Bayesian approach with Markov chain Monte Carlo methods, which we present in this paper for the first time. We applied this new method to a subset of 30 galaxy clusters representative of the CHEX-MATE sample and show that we can obtain reliable temperature measurements up to regions where the source intensity is as low as 20% of the background, keeping systematic errors below 10%. We compare the median profile of our sample and the best-fit slope at large radii with literature results and we find a good agreement with other measurements based on *XMM-Newton* data. Conversely, when we exclude the most contaminated regions, where the source intensity is below 20% of the background, we find significantly flatter profiles, in agreement with predictions from numerical simulations and independent measurements with a combination of Sunyaev–Zeldovich and X-ray imaging data.

**Key words.** galaxies: clusters: general – galaxies: clusters: intracluster medium

# 1. Introduction

Galaxy clusters represent the end point of structure formation in the Universe: they are the most massive virialized structures to have formed thus far and are located at the nodes of the cosmic web, continuously accreting matter from the filaments. The intra-cluster medium (ICM) constitutes their main baryonic component and is heated to X-ray-emitting temperatures ($10^7-10^8$ K) during the formation and accretion processes driven by dark matter. As such, the ICM carries important information on the physical processes of structure formation and on the global cluster properties: ICM thermodynamic properties correlate well with the total gravitational mass and properly scaled radial profiles of thermodynamic quantities are nearly universal, reflecting the properties of the underlying dark matter structure (Kaiser 1986). Deviations from these relations do exist and are a clear indication of important effects such as cooling, nongravitational feedback from active galactic nuclei (AGN), or supernovae, bulk motions, and turbulence induced by accretion processes (see Lovisari & Maughan 2022 for a recent review, Gaspari et al. 2014; Bulbul et al. 2019; Sereno et al. 2020; Poon et al. 2023).

The ICM temperature is one of the key observable quantities that can be derived from X-ray observations of galaxy clusters, along with density and metal abundance. It is derived from the analysis of X-ray spectra, mainly through the position of the Bremsstrahlung exponential cutoff. In the last 30 years, the previous and present generation of X-ray satellites have allowed for spatially resolved spectral measurements, mapping the distribution of ICM temperatures from the cores to the external regions. Temperature profiles in radial annuli are a necessary ingredient to derive the radial distributions of thermodynamic quantities, such as pressure and entropy, and to reconstruct the total mass through the hydrostatic equilibrium equation. Early temperature profiles with ASCA and BeppoSAX (e.g., Fabian et al. 1994; De Grandi & Molendi 2002) unambiguously demonstrated that some clusters do show a temperature decline in their inner regions, consistent with a short cooling time of the high-density central regions and the prediction of the cooling flow model (Fabian et al. 1984; Fabian 1994). When high resolution spectra with XMM-Newton RGS showed that the gas does not cool below a critical value (e.g., Peterson et al. 2001, 2003), causing the fall of the cooling flow model, these systems were dubbed cool cores (CCs, Molendi & Pizzolato 2001) and were shown to be typically associated with clusters in a relaxed dynamical state. Conversely, in the external regions, ASCA and BeppoSAX observations, lead to controversial results on the presence of a temperature gradient (e.g., Markevitch et al. 1998; Irwin et al. 1999; White 2000; De Grandi & Molendi 2002). This is due to

<footnote>* Full Table C.1 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/686/A68</footnote>
technical challenges in measuring temperature profiles in the external regions of galaxy clusters, where the intensity of the source is low with respect to contamination from other brighter cluster regions because of the large point spread function (PSF) and several background components (see Ettori & Molendi 2011 for a review). The former issue has been resolved with the advent of the current generation of X-ray telescopes *XMM-Newton* and *Chandra*, which unambiguously showed a declining trend in the temperature profiles with radius in the external regions (e.g., Vikhlinin et al. 2006; Pratt et al. 2007; Snowden et al. 2008; Leccardi & Molendi 2008). However, the low intensity of the cluster emission with respect to the instrumental and celestial background and foreground is still a major challenge, hampering reliable temperature measurements up to the virial radius of galaxy clusters and thus leaving a large fraction of the cluster volume poorly explored. In the last decade, the *Suzaku* satellite, benefiting from a low-earth orbit and thus a lower particle background with respect to *XMM-Newton* and *Chandra*, allowed us to extend temperature measurements up to the virial radius for a limited sample of a dozen clusters. The steep decrease in the temperature profile and flattening in the entropy at the virial radius observed in some clusters is still under debate (see the review by Walker et al. 2019 and references therein).

An independent method that avoids the issues of X-ray spectroscopy for deriving temperature profiles is the combination of the spatially resolved Sunyaev–Zeldovich (SZ, Sunyaev & Zeldovich 1972) signal, which is proportional to the integral of the ICM pressure along the line of sight, with the gas density, which can be derived from X-ray images and is less dependent on the background treatment than the temperature. This approach has been attempted for a few clusters (e.g., Basu et al. 2010; Eckert et al. 2013; Tchernin et al. 2016; Pratt et al. 2016; Adam et al. 2017; Sereno et al. 2018; Ghirardini et al. 2021). In the X-COP project, Ghirardini et al. (2019) compared the temperature profiles obtained with the combination of *Planck* and *XMM-Newton* data with those obtained with the purely spectroscopic XMM analysis and found them to be consistent in the radial range where the two methods overlap, that is within $R_{500}$. While this approach proved very successful for measuring temperatures in the outskirts of galaxy clusters, thus far it has been limited to clusters with a high enough SZ signal and extension in the sky matching the capabilities of the SZ instruments (nearby extended objects for low resolution instruments such as *Planck* or compact systems for high resolution instruments such as NIKA2 and ALMA). Moreover, the spatial resolution of spectroscopic X-ray temperature profiles is still superior to that obtained with the joint SZ-X-ray analysis.

The use of radial profiles for studying the temperature structure of the ICM is of course an approximation: several authors have shown two-dimensional maps of the ICM temperatures, with significant non-radial variations (e.g., Bourdin & Mazzotta 2008; Frank et al. 2013; Laganá et al. 2019). Nonetheless, the spherical symmetry approximation is necessary for reconstructing the total mass with the hydrostatic equilibrium equation and in conditions of limited statistics, such as in faint and distant clusters or in the external regions. Recently, Lovisari et al. (2024) quantified the impact of temperature inhomogeneities on the reconstructed radial profiles to be below 5% on average, reaching larger values (10–20%) in external and peculiar regions.

The Cluster HEritage project with *XMM-Newton*: Mass Assembly and Thermodynamics at the End point of structure formation (CHEX-MATE, CHEX-MATE Collaboration 2021) is

a multi-year heritage program with *XMM-Newton* awarded in 2017. It consists of a homogeneous set of observations of a representative sample of 118 galaxy clusters, extracted from the *Planck* SZ catalog (PSZ2, Planck Collaboration XXVII 2016). It aims to derive the total masses, through the hydrostatic equilibrium equation, and characterize the thermodynamic properties of the full sample (see CHEX-MATE Collaboration 2021, for more details). The observations are tailored to reach a precision of 15–20% on the hydrostatic masses within $R_{500}$. In this context, the temperature profile is a fundamental ingredient of the CHEX-MATE analysis and the desired precision on the mass translates to a 15% precision on the temperature measurement at $R_{500}$, that is, the radius within which the mean density is 500 times the critical density of the Universe at the cluster’s redshift. Any issue affecting the temperature profiles will inevitably be propagated into the final CHEX-MATE results, and it is thus crucial for our project to obtain the most accurate temperature profiles from our data and assess the impact of our choices and of possible systematics on the temperature measurements. This is why we developed a new pipeline for spectral extraction and analysis (Sect. 3), with a special focus on the background model (Sect. 4), which we present in this paper for the first time. We applied this method to a sample of 30 objects (Sect. 2), which we used to test it and evaluate the impact of our analysis on the temperature measurements (Sect. 5). We present the first results on the temperature profiles for this sample in Sect. 6 and we discuss our findings in Sect. 7.

In our paper, we assume a flat $\Lambda$CDM cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.3$, and $\Omega_\Lambda = 0.7$. Errors are given at a 68% confidence level, unless otherwise stated. Whenever we use the notation $M_{500}$ in this paper, we are referring to the masses derived from the *Planck* SZ signal with the method described in Planck Collaboration XXVII (2016), using the MMF3 algorithm (Melin et al. 2006), as discussed in CHEX-MATE Collaboration (2021).

## 2. The sample

We built our sample to fulfill two different requirements: (i) being a “technical” sample, allowing us to test our pipeline and new methods with a limited number of objects under different analysis conditions (source angular extension, background levels, etc.); (ii) being representative of the full CHEX-MATE sample, in terms of its selection quantities (mass, redshift, and *Planck* S/N). We dubbed it “Data Release 1” (DR1 hereafter) because it is the first sample for which results have been internally released to the CHEX-MATE collaboration. To build the sample, we selected clusters in the $M–z$ plane (Fig. 1), looking also at the distribution of other physical parameters, such as the extension in the sky as measured by $R_{500}$ (ranging from 2.9 to 13.5 arcmin) and Galactic absorption ($N_H$ spanning one order of magnitude from $1.03 \times 10^{20} \text{ cm}^{-2}$ to $1.02 \times 10^{21} \text{ cm}^{-2}$) with different levels of molecular contribution (Bourdin et al. 2023). We also looked at the distribution of the quality assessment indicators that we computed for each observation (see Sect. 3.1), making sure that the selected observations encompass different levels of cosmic-ray induced particle background (Sect. 4.1), following the solar cycle (Marelli et al. 2021; Gastaldello et al. 2022), and of the residual focused component (Sect. 4.2, Salvetti et al. 2017). We also checked that the clean exposure time of the selected observations allows us to meet the feasibility requirements of the CHEX-MATE project.

We confirm our sample definition by checking that the DR1 sample is representative of the original CHEX-MATE
<table><caption>**Table 1.** Clusters of the DR1 sample.</caption><thead><tr><th>PSZ2 name</th><th>Redshift</th><th><i>M</i><sub>500</sub><br>(10<sup>14</sup> <i>M</i><sub>☉</sub>)</th><th><i>R</i><sub>500</sub><br>(arcmin)</th><th>Tier</th><th>XMM obsid<br>center/offset</th></tr></thead><tbody><tr><td>PSZ2 G008.31-64.74</td><td>0.312</td><td>7.42</td><td>4.51</td><td>2</td><td>0827010901</td></tr><tr><td>PSZ2 G041.45+29.10</td><td>0.178</td><td>5.41</td><td>6.48</td><td>1</td><td>0601080101</td></tr><tr><td>PSZ2 G042.81+56.61</td><td>0.072</td><td>4.22</td><td>13.50</td><td>1</td><td>0202080201/0827361101</td></tr><tr><td>PSZ2 G046.88+56.48</td><td>0.115</td><td>5.10</td><td>9.40</td><td>1</td><td>0827010601</td></tr><tr><td>PSZ2 G050.40+31.17</td><td>0.164</td><td>4.22</td><td>6.40</td><td>1</td><td>0827040101</td></tr><tr><td>PSZ2 G056.77+36.32</td><td>0.095</td><td>4.34</td><td>10.53</td><td>1</td><td>0740900101</td></tr><tr><td>PSZ2 G056.93-55.08</td><td>0.447</td><td>9.49</td><td>3.70</td><td>2</td><td>0503490201</td></tr><tr><td>PSZ2 G057.78+52.32</td><td>0.065</td><td>2.32</td><td>12.15</td><td>1</td><td>0827040301/0827041801</td></tr><tr><td>PSZ2 G057.92+27.64</td><td>0.076</td><td>2.66</td><td>11.08</td><td>1</td><td>0827030301</td></tr><tr><td>PSZ2 G066.41+27.03</td><td>0.575</td><td>7.69</td><td>2.88</td><td>2</td><td>0827320601</td></tr><tr><td>PSZ2 G072.62+41.46</td><td>0.228</td><td>11.43</td><td>6.73</td><td>2</td><td>0605000501</td></tr><tr><td>PSZ2 G077.90-26.63</td><td>0.147</td><td>4.99</td><td>7.46</td><td>1</td><td>0827020101</td></tr><tr><td>PSZ2 G083.86+85.09</td><td>0.183</td><td>4.74</td><td>6.04</td><td>1</td><td>0827030701</td></tr><tr><td>PSZ2 G113.29-29.69</td><td>0.107</td><td>3.57</td><td>8.86</td><td>1</td><td>0827021201</td></tr><tr><td>PSZ2 G113.91-37.01</td><td>0.371</td><td>7.58</td><td>3.96</td><td>2</td><td>0827021001</td></tr><tr><td>PSZ2 G114.79-33.71</td><td>0.094</td><td>3.79</td><td>10.20</td><td>1</td><td>0827320401</td></tr><tr><td>PSZ2 G149.39-36.84</td><td>0.170</td><td>5.35</td><td>6.71</td><td>1</td><td>0827030601</td></tr><tr><td>PSZ2 G195.75-24.32</td><td>0.203</td><td>7.80</td><td>6.53</td><td>2</td><td>0201510101</td></tr><tr><td>PSZ2 G207.88+81.31</td><td>0.353</td><td>7.44</td><td>4.09</td><td>2</td><td>0827020301</td></tr><tr><td>PSZ2 G224.00+69.33</td><td>0.190</td><td>5.11</td><td>5.99</td><td>1</td><td>0827020901</td></tr><tr><td>PSZ2 G238.69+63.26</td><td>0.169</td><td>4.17</td><td>6.21</td><td>1</td><td>0500760101</td></tr><tr><td>PSZ2 G243.15-73.84</td><td>0.410</td><td>8.09</td><td>3.75</td><td>2</td><td>0827011301</td></tr><tr><td>PSZ2 G243.64+67.74</td><td>0.083</td><td>3.62</td><td>11.23</td><td>1</td><td>0827010801</td></tr><tr><td>PSZ2 G277.76-51.74</td><td>0.438</td><td>8.65</td><td>3.65</td><td>2</td><td>0674380301</td></tr><tr><td>PSZ2 G287.46+81.12</td><td>0.073</td><td>2.56</td><td>11.33</td><td>1</td><td>0149900301</td></tr><tr><td>PSZ2 G313.33+61.13</td><td>0.183</td><td>8.77</td><td>7.42</td><td>2</td><td>0093030101</td></tr><tr><td>PSZ2 G313.88-17.11</td><td>0.153</td><td>7.86</td><td>8.37</td><td>2</td><td>0692932001</td></tr><tr><td>PSZ2 G324.04+48.79</td><td>0.452</td><td>10.58</td><td>3.81</td><td>2</td><td>0112960101</td></tr><tr><td>PSZ2 G340.94+35.07</td><td>0.236</td><td>7.80</td><td>5.76</td><td>2</td><td>0827311201</td></tr><tr><td>PSZ2 G349.46-59.95</td><td>0.347</td><td>11.36</td><td>4.77</td><td>2</td><td>0504630101</td></tr></tbody></table>

**Notes.** We quote: the PSZ2 name, the redshift, the nominal integrated mass from Planck MMF3 data, the cluster size $R_{500}$ in arcminutes, as derived from the mass with Eq. (9) in Planck Collaboration XX (2014), the Tier to which each cluster belongs (see CHEX-MATE Collaboration 2021 for details) and the XMM-Newton observations used in the analysis.

selection (CHEX-MATE Collaboration 2021), by performing Kolmogorov-Smirnov (KS) tests on the distribution of mass, redshift, or Planck S/N. Since we excluded from our selection double clusters that would require an ad hoc strategy in the definition of regions for the radial profile (Sect. 3.2), we further checked that we did not preferentially select relaxed clusters, by performing a KS test on the distribution of the *M* morphological parameter, which combines the information from the light concentration and centroid shift parameters and from the second and third order power ratios (see Campitiello et al. 2022 for details). The test returns a very high probability (>99%) that DR1 is representative of the morphological distribution of the original CHEX-MATE sample.

The final DR1 sample is composed of 30 clusters, listed in Table 1, along with their main physical properties and the XMM-Newton observations that we used in our analysis. We show their distribution in the mass-redshift plane along with the parent CHEX-MATE sample in Fig. 1. We recall that masses here are estimated from the Planck SZ data as described in Planck Collaboration XXVII (2016) using the $Y_{SZ}-M_{500}$ scaling relation calibrated with XMM-Newton data and are not corrected for any hydrostatic bias.

We show the images of all clusters in our sample in Appendix C.

## 3. The CHEX-MATE pipeline for spectral results

In this section and in Sect. 4, we present in details the pipeline that we developed to extract and fit spectra in radial profiles for the CHEX-MATE project. We tested it with other pipelines available in the collaboration and used for other CHEX-MATE projects (Bourdin et al. 2023; Lovisari et al. 2024; De Luca et al., in prep.). In the internal regions of galaxy clusters ($R < 0.5 R_{500}$) the temperatures estimated with different analysis are consistent within 1–2%. Differences raise to about 10% at $R_{500}$, where the details of the background modeling have a large impact on the temperature measurements.

### 3.1. Data reduction

We performed data reduction and carried out our analysis making use of the XMM-Newton Science Analysis System (SAS), version 16.1 and the Extended Source Analysis Software (ESAS, Snowden et al. 2008; Kuntz & Snowden 2008) embedded within SAS. We keep the calibration database (CalDB) updated to ensure that we apply the latest calibration files to our data, when we reprocess them starting from the raw ODF files retrieved from the archive. We tested more recent SAS versions (19.1 and 20.0), but the tasks are either much slower than with SAS 16.1,
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Distribution of the CHEX-MATE (black circles) and DR1 (orange open squares) objects in the mass-redshift plane, showing that our subsample is representative of the original selection.</figcaption></figure>

thus requiring a too long computing time for a large data-set as CHEX-MATE, either show unresolved bugs. We are supported in our choice by the EPIC calibration status document<fnref n="1" />, showing that no major change has been applied in the meantime, with the only exception of the optional cross-calibration with NUSTAR<fnref n="2" /> that we are not using in our pipeline.

We refer to Bartalucci et al. (2023) for details on the first steps of the data reduction, including calibration, standard pattern cleaning, removal of noisy CCDs in the MOS cameras and light-curve filtering. We emphasize that for our analysis we used the same selection of contaminating sources as in Bartalucci et al. (2023), and thus we mask both point sources and extended substructures. The latter include known merging substructures, such as in PSZ 2G046.88+56.48, a.k.a. A2069, where the northern subcluster A2069B (e.g., Owers et al. 2009) has been masked in our analysis, and PSZ 2G195.75–24.32, a.k.a. A520, where we masked the merging subcluster, dubbed as “the foot” in Wang et al. (2016) and leading to the formation of a shock front (Markevitch et al. 2005), and other extended features whose origin is still unclear (in PSZ 2G056.77+36.32, PSZ 2G057.78+52.31, PSZ 2G057.92+27.64). Masked sources and regions for all clusters are shown in Appendix C.

While reducing each observation, we also compute two indicators that allow us to estimate the level of instrumental background affecting them. The first one is the count-rate in the region of the MOS2 detector not exposed to the sky (i.e., outside the Field Of View, outFOV), which has been shown to be a good proxy of the background due to the interaction of high-energy (>10 MeV) cosmic ray particles, modulated by the solar cycle (Sect. 4.1, see Gastaldello et al. 2022). We calculate this parameter following the prescriptions (detector, energy band and region definition) provided in Marelli et al. (2017). The second indicator (inFOV–outFOV) is the difference between the count-rate in the region exposed to X-ray photons of the MOS2 detector

(inFOV) and the outFOV value, both computed in a hard band<fnref n="3" /> to minimize the celestial contribution. This allows us to estimate the excess background in the FOV with respect to the background level predicted by the outFOV value (e.g., Salvetti et al. 2017). We refer to this component as the residual focused component, as it is at least partly due to contamination from quiescent soft protons not filtered in the light curve cleaning but may also have a different origin (see Salvetti et al. 2017). To compute the inFOV–outFOV indicator we could not use the count rate in the full MOS2 FOV as in Marelli et al. (2017) and Salvetti et al. (2017) because of the presence of the cluster emission in our observation. We thus used an external annulus between 12 and 14 arcmin from the aimpoint which, in combination with the hard energy band, minimizes the contribution of the cluster emission. These indicators have been very useful to assess the quality of our observations. We subsequently discuss in Sect. 4 how we used them to build our background model.

## 3.2. Spectral extraction in annuli

We define regions for the spectral extraction as concentric annuli centered on the peak of the X-ray image. We automatically define the position of the peak as the brightest pixel in the soft band image (0.7–1.2 keV), after correcting for the exposure map, masking the point sources, and smoothing with a Gaussian function with $\sigma = 5$ pixels (see Bartalucci et al. 2023 for more details). We compute the width of each annulus with an adaptive binning method to reach a constant signal to noise ratio ($S/N = 50$ in the 0.3–2 keV energy range) in each spectral bin (e.g., Pratt et al. 2010; Chen et al. 2023).

For each region, we extract spectra, redistribution matrix files (RMF), and ancillary response files (ARF) with the ESAS tools *mos-spectra* and *pn-spectra*. Following the ESAS procedure (Kuntz & Snowden 2008; Snowden et al. 2008), these tools also compute the count rate in the unexposed corners (outFOV) and select the most similar (in terms of magnitude and hardness ratio, see Kuntz & Snowden 2008 for details) filter-wheel-closed observations in the *XMM-Newton* calibration database. A spectrum for the cosmic-ray induced particle background (CRPB) is then produced for each region with the *mos-back* and *pn-back* tools. We use these ESAS products to build our background model as shown in Sect. 4.1.

The regions for spectral extraction can be very small (width as small as 3.3′′) in the central regions of galaxy clusters with a peaked surface brightness profile. In these conditions, we noticed an error in the *arfgen* and *rmfgen* tools called by *mos/pn-spectra* producing null ARFs and RMFs. To fix this issue, we extract RMFs, ARFs in a circular region with radius 30′′ around the peak and we assign them to all regions in the inner 30′′. Since the statistical quality of the background files in small regions is very poor, we also use the appropriately rescaled CRPB spectra in the circular region with $R = 30′′$ for all central regions. Both the effective area and the particle background intensity do not show significant variations on scales smaller than 1′.

Finally we extracted a spectrum from an external “background” region, ideally free of cluster emission, which is needed to calibrate our model of the sky background components (Sect. 4.3). The definition of this region is not trivial in a sample where the extension of the cluster emission in the FOV is very

<footnote n="1">¹ https://xmmweb.esac.esa.int/docs/documents/CAL-TN-0018.pdf</footnote>

<footnote n="2">² https://xmmweb.esac.esa.int/docs/documents/CAL-SRN-0388-1-4.pdf</footnote>

<footnote n="3">³ As described in Marelli et al. (2017), we use events with $E > 7$ keV, excluding the energy ranges [9.4–10] keV and [11–12] keV, affected by instrumental lines.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Examples of our spectral fits for the cluster PSZ2 G349.46−59.95 in different regions: the external region where we fit the sky background (left), a central region dominated by the cluster emission (middle), and an external background-dominated region with SOU/BKG = 0.4 (right). In all panels, the black, red, and green crosses mark the MOS1, MOS2, and pn spectra, while the lines represent the corresponding models. With the same color scheme, the continuous lines mark the total resulting model for each detector, the dotted lines represent the model for particle background, while the dashed lines the sky background. The dotted cyan line (when visible) mark the contribution of the pn OoT events. In the middle and right panel, we highlight in blue the cluster component. In the left panel, magenta crosses represent the data from the RASS spectrum, which is jointly fitted to the EPIC data in our baseline pipeline. In all plots, the bottom panels mark the ratio of the data with respect to the total model.</figcaption></figure>

different ($R_{500}$ ranging from 2.9' to 13.5'). For all the clusters with $R_{500} \le 9'$, we define an annulus between $R_{200} = R_{500}/0.7$ and the edge of the FOV. This is also possible for the two objects (PSZ2 G057.78+53.31 and PSZ2 G042.81+56.61) with $R_{500} > 12'$ for which we requested an offset observation for this purpose (CHEX-MATE Collaboration 2021). Conversely, for the six clusters with $9' < R_{500} < 12'$ we define a close background (closebkg) region as an annulus between $1.1 R_{500}$ and the FOV edge. We separate the residual cluster emission from the celestial background components, making use of ancillary data from the ROSAT All Sky Survey (RASS) diffuse background map (Snowden et al. 1995, 1997), as described in detail in Sect. 4.3. We subsequently discuss the impact of this strategy in Sect. 5.4.1.

### 3.3. Spectral fitting

For each region, we joint-fit the spectra of the three detectors and eventually of multiple observations covering the same region within XSPEC (version 12.11). We show some examples of our fits in Fig. 2. We link together the parameters of the models which are applied to the three instruments and to all observations, taking into account the different areas by using a constant factor. Indeed, all the normalization terms in our fits are per unit area (in arcmin²). Unless otherwise stated, we perform our fits in the band 0.5–12 keV for the MOS detectors and 0.5–14 keV for the pn. From these bands, we further exclude the energy ranges corresponding to the more prominent fluorescent lines, namely Al and Si in the range 1.2–1.9 keV for the MOS detectors and Al in 1.2–1.7 keV and Cu in 7–9.2 keV for the pn. As in Leccardi & Molendi (2008) and Ghirardini et al. (2019), we do not subtract any background component in our fit but we model them with a physical model described in Sect. 4 that we apply to the spectra in each region. We estimate the parameters of the celestial components by fitting the spectra in the external background region (see Sect. 3.2), with the model described in Sect. 4.3. Also the pn Out-of-Time (OoT) events are treated as a further background component: we extract a spectrum from the OoT event files in each region as described in Sect. 3.2 and we fit it with a *phabs(bapec+powerlaw)* model, to mimic the cluster contribution (but with a degraded energy resolution) and the CXB, and another power law for the particle background. We

then properly rescale the best fit model of the OoT contribution, taking into account the operating mode (Full Frame or Extended Full Frame).

The baseline method for our fits makes use of Bayesian statistics and of the Markov chain Monte Carlo (MCMC) algorithm (Goodman & Weare 2010) within XSPEC. To assess and test it, we also perform more standard fits with a standard Maximum Likelihood approach (ML). In both cases, we use the Cash statistics (Cash 1979) within XSPEC, as it is more appropriate than $\chi^2$ for Poissonian data, especially in conditions of high background and poor statistics (Leccardi & Molendi 2007; Humphrey et al. 2009; Kaastra 2017). In the ML approach, we build our “best” model for the background components (Sect. 4) and we freeze it, leaving as the only free parameters in our fits those related to the cluster emission (temperature, metal abundance and normalization). In this framework, the best fit parameters are the values minimizing the Cash statistics, and the uncertainties are computed imposing $\Delta C_{stat.} = 1$. Conversely, in the Bayesian MCMC we estimate priors for the parameters associated to the background components (Sect. 4) that are then joint-fitted with the parameters associated to the cluster emission. This method allows us to propagate the uncertainties on the background parameters to the scientific results. In this approach, we compute the best fit values as the mean of the marginal posterior distribution of each parameter and the $1\sigma$ errors by sorting the chain values and then taking the upper and lower limit of the central 68% of the values. We run the MCMC fit within XSPEC, with the *chain* command, using the Goodman–Weare algorithm (Goodman & Weare 2010) with 26 walkers (twice the number of free parameters). We tried different lengths both for the burn-in and for the fit phase on a few typical cases and we noticed to have a good convergence and stable results discarding the first 50 000 steps (burn-in phase) and running the chain for other 50 000 iterations.

Concerning the source component, we model the cluster emission with *apec* (Astrophysical Plasma Emission Code Smith et al. 2001) and the Galactic absorption with *phabs*, assuming *aspl* abundance table (Asplund et al. 2009) and *vern* photo-electric cross sections (Verner et al. 1996). We discuss in Appendix B the impact of the abundance table on the parameters of cluster emission. We fix the redshift at the values provided in CHEX-MATE Collaboration (2021) and the column density
to the value $N_{Htot}$ which includes the contribution of both the atomic hydrogen ($N_{HI}$) and the molecular one (see Bourdin et al. 2023 for more details). We explain how we tested the impact of using the total $N_H$ with respect to $N_{HI}$ in Sect. 5.2.

## 4. A physical model for the background components

The typical accuracy of temperature profiles in cluster outskirts is largely determined by a telescope's background level and our ability at predicting the normalization and spectral shape of background components. Here we introduce a physically motivated model for the XMM-Newton background which allows us to accurately predict its properties. The XMM-Newton background is made of several components of different physical origin, in the form of the cosmic-ray induced particle background (CRPB), the residual focused component (RFC), and photons (cosmic X-ray background and foreground emission). We modeled individually these components and calibrated them on a large set of over 500 blank-sky observations from the XMM XXL survey (Pierre et al. 2016), totaling more than 6 Ms of data. The adopted procedure, its calibration and an estimation of the corresponding background subtraction accuracy are described in the following.

### 4.1. Cosmic-ray particle-induced background

The primary background component in XMM-Newton data originates from the interaction between Galactic cosmic rays and the structure of the spacecraft, generating secondary electrons. The secondary electrons are recorded by the CCDs and cannot be individually distinguished from photon-induced events, which produces a background signal that will be referred to as the cosmic-ray particle induced background. The normalization of this component is known to be modulated by the Solar cycle (Gastaldello et al. 2022) and its spectral shape resembles a hard power law with a number of fluorescence lines (Kuntz & Snowden 2008). In the MOS detectors, the CRPB intensity can be monitored on-the-fly using the unexposed corners of the detectors, which are shielded from the telescope's FOV and thus record only CRPB-induced events. On the other hand, a large collection of data with the filter wheel set in closed position (hereafter filter-wheel-closed data, FWC) were accumulated by XMM-Newton throughout the mission. Again, FWC observations only record CRPB-induced events, allowing us to determine the spectral shape and the variation of the CRPB across the detector with high accuracy. For each observation, we determine the current CRPB level using the MOS corners and renormalize the FWC data such that they match the level observed in the corner data using the ESAS tasks mos-spectra and mos-back. We then use the rescaled FWC data to extract the CRPB spectrum from the region of interest, and we fit the resulting spectrum with a phenomenological model made of two power laws and a combination of Gaussian fluorescence lines to create an analytical model for the CRPB.

In the case of the pn detector, the situation is rendered more complicated by the fact that the device does not include an area that is fully shielded from the telescope's FOV (Marelli et al. 2021). Therefore, it is not possible to monitor the CRPB rate in the same way as for the MOS, especially in time frames that are significantly affected by soft proton (SP) flares. Following Marelli et al. (2021), we define an outer region ($R > 905''$ from the aimpoint) which is reasonably shielded from the telescope.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Relation between the pn and MOS2 high-energy (10–14 keV) corner count rates in a large number of blank-sky observations with a low level of residual focused contamination. The color code indicates the level of contamination as traced by the inFOV-outFOV indicator and the dotted line is the best-fit relation (Eq. (1)).</figcaption></figure>

We then take advantage of the simultaneous MOS data to predict the level of CRPB in the pn. To this end, we make use of the inFOV and outFOV count rates, as introduced in Sect. 3.1. We recall that inFOV-outFOV can be used as an indicator of the remaining SP contamination (Salvetti et al. 2017). We then select a subset of our blank-sky observations with low contamination (inFOV-outFOV < 0.04 cts s⁻¹), for which we expect the pn corner data to be unaffected by soft protons and representative of the CRPB only. We then correlated the [10–14] keV pn count rate in the corners with the MOS2 corner count rate and found an excellent correlation between the two quantities, as already reported by Marelli et al. (2021).

In Fig. 3 we show the relation between the pn and MOS2 corner count rates for observations with inFOV-outFOV < 0.04 cts s⁻¹. The tight correlation demonstrates that for uncontaminated observations the pn corner count rates provide an accurate estimate of the CRPB level. However, the color code implemented in Fig. 3 shows that the scatter of the points around the mean relation correlates with the level of contamination, with the observations falling above the relation exhibiting systematically higher values of inFOV-outFOV than the relations falling below it. To alleviate this issue, we fit a simple model to the data in which the pn high-energy count rate is described as

$$ \text{outFOV}_{\text{PN}} = A_{\text{CRPB}} \text{outFOV}_{\text{MOS2}} + A_{\text{SP}}(\text{inFOV} - \text{outFOV}). \quad (1) $$

In other words the pn outFOV rate is the sum of a CRPB component that is related to the MOS2 corner count rate with a proportionality constant $A_{CRPB}$ and of a SP contamination term that is proportional to inFOV-outFOV with a proportionality constant $A_{SP}$. We then fitted the relation shown in Fig. 3 with this model and determined the best-fitting values of $A_{CRPB}$ and $A_{SP}$. We are then able to predict the level of the CRPB-induced pn rate per unit area as

$$ \text{CRPB}_{\text{PN}} = A_{\text{CRPB}} \text{outFOV}_{\text{MOS2}}. \quad (2) $$

For each observation and each spectrum, we extract the local CRPB spectrum from the FWC data using pn-spectra and renormalize it such that the high-energy normalization matches the prediction of Eq. (2).
While this method has proven to be very effective in reproducing the shape and normalization of the continuum emission, it does not reproduce well the properties of the emission lines, whose central energy and width are sensitive to gain variations and charge-transfer-inefficiency (see Kuntz & Snowden 2008; Sanders et al. 2020 for more details). This is why we decided to exclude the energy ranges of the most prominent spectral lines in our fit, as described in Sect. 3.3.

## 4.2. Residual focused component

The XMM-Newton orbit is known to host a large number of clouds of low-energy protons (~100 MeV) that are trapped within the Earth's magnetosphere. When the trajectory of the spacecraft crosses such a cloud, soft protons (SP) can be funneled through the telescope and focused in a similar way as actual photons, albeit with a slightly different vignetting curve (Kuntz & Snowden 2008). Soft protons can generate very bright flares of the background intensity, with the background count rates increasing by more than two orders of magnitude. While the light curve filtering process introduced in Sect. 3.1 allows us to exclude time periods that are strongly affected by soft proton flares, in some cases there remains a low level of SP contamination even after the most affected time intervals have been clipped, which constitutes a significant part of the residual focused component. We attempted to create a predictive model of the residual SP contamination from the set of ~500 blank-sky pointings discussed above, as described in Appendix A. The SP component is modeled as a power law $I(E) = N_{QSP}E^{-\Gamma}$ with a fixed slope ($\Gamma = 0.6$, Appendix A.1), while its normalization $N_{QSP}$ is predicted from the inFOV-outFOV indicator, using the calibration expressions described in Appendix A.2, one for each instrument. We take into account the variation of the SP component across the field-of-view with its vignetting function, as provided by the ESAS task protonscale. In this way, we can predict the value of $N_{QSP}$ for each region, detector, and observation.

## 4.3. Sky background and foreground components

The sky contribution to our spectra is mainly due to three components (e.g., Kuntz & Snowden 2000): the residual Cosmic X-ray Background (CXB), and the foreground emission of the Galactic Halo (GH), and the emission of the Local Hot Bubble (LHB). For the LHB, we use an unabsorbed apec thermal model with temperature fixed to 0.11 keV, solar metal abundance, and $z = 0$, while for the GH we used an absorbed apec with temperature allowed to vary in the range 0.1–0.6 keV (flat prior in the MCMC fit) as in McCammon et al. (2002) and solar metal abundance<fnref n="4" />. In both cases, the normalizations are free parameters of the fit. We model the CXB as an absorbed power-law, with a fixed slope of 1.46 (Moretti et al. 2009) and free normalization. For the absorption, we always use the phabs model and the same $N_H$ value used to fit the cluster emission (Sect. 3.3). We derive the free parameters of this model by fitting the spectra extracted in an external annulus, ideally free of cluster emis-

sion (Sect. 3.2). While running fits with the Bayesian MCMC approach, we also used this approach to estimate the sky background: we derived posteriors for the four free parameters that we modeled with Gaussian distributions to be given as priors in the subsequent fits.

As discussed in Sect. 3.2, for six clusters in our sample we could not define a region not contaminated by the cluster emission in the FOV and offset observations are not available. In these conditions, we need to take into account a residual cluster component in the “closebkg” region, which is very degenerate with the sky background parameters that we need to estimate. To break this degeneracy and provide better constraints to the free parameters in our model, we follow the method suggested by Snowden et al. (2008) and we make use of an ancillary data set: a spectrum extracted from the ROSAT All Sky Survey (RASS) diffuse background map (Snowden et al. 1995, 1997) in an annulus between one and two degrees from the cluster center, using the HEASARC X-ray background tool<fnref n="5" />. We jointly fit the RASS spectrum with the XMM one, taking into account a cross-correlation factor between the two instruments of 15% (Eckert et al. 2011) and binding together the parameters of the sky background components. The normalization of the cluster component is a free parameter in the XMM spectrum while it is fixed to zero in the external RASS data, which essentially do not contain cluster emission. During the fit of the sky background spectra, we fix the temperature of the cluster component beyond $R_{500}$, since leaving it free often lead us to unphysical values for this parameter. The fit is not very sensitive to the exact value of this parameter and we chose to fix it to one third of the mean expected temperature of the cluster, as estimated by the $M-T$ relation in Arnaud et al. (2005), taking into account the typical decline of temperature profiles in the literature. We perform the fits in the background region both with and without the RASS ancillary dataset, to study its impact both on the sky parameters estimates and on the temperature profiles (Sect. 5).

## Impact of RASS joint fit on sky background parameters

We note that the joint fit with the RASS affects significantly the best fit values of some sky components: the normalization of the LHB, which is only poorly constrained in the XMM-only fit, the normalization and temperature of the Galactic Halo (see Fig. 4). This variation is systematic: with the RASS+XMM joint fit we find lower LHB normalization (16%, computed as the weighted mean of the ratio XMM+RASS/XMMonly), higher normalization (55%) and lower temperature (16%) of the GH, reflecting the strong degeneracies between these parameters.

For the CXB normalization, we should consider two different effects. For the more extended clusters, where we do not have a clean region in the field of view to estimate the sky background, the fit without RASS should overestimate the CXB normalization, since the residual cluster emission will be attributed to the background components. Indeed, if we consider only the 6 extended clusters ($9' < R_{500} < 12'$) where we used a “closebkg” region (purple points in Fig. 4), the normalization of the fit without RASS is larger than the one with the joint fit, with a mean difference of 17%. Besides that, there is another fundamental effect: since the source detection limit in the XMM-Newton observations is lower than in the RASS data, the fraction of the resolved CXB is different in the two spectra, leading to different normalizations for this component, while in our joint

<footnote n="4">⁴ Recently, Ponti et al. (2023) reported a much lower value for the abundance of the hot circum-galactic medium with eROSITA data. The value of the metal abundance is degenerate with the other free parameters, especially the normalization of the GH and LHB. We performed simulations of a typical sky spectrum with $Z = 0.06 Z_⊙$ and found that we can always retrieve a good fit with our procedure, even assuming $Z = Z_⊙$. This suggests that we can still use our baseline model as a good phenomenological representation of the sky background components.</footnote>

<footnote n="5">⁵ https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/xraybg/xraybg.pl</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of best-fit values for sky background parameters obtained with the RASS joint fit (x-axis) and using XMM-Newton only (y-axis). We show the normalization of the LHB (upper left), and of the CXB (upper right), while in the bottom panels we report the GH normalization (left) and temperature (right). Points are color-coded by the extension of the cluster, estimated by R<sub>500</sub>, in arcminutes. The dashed line marks the equality line. The LHB and GH components show large systematic variations, reflecting the degeneracies between the components. The CXB normalization shows only a moderate variation on the full sample, with larger difference for extended clusters.</figcaption></figure>

fit we impose them to be the same. We can quantify this effect,
by comparing the best fit normalization with the joint fit and
the XMM-only, in the objects where we have a large enough
cluster-free region in the *XMM-Newton* observations. Indeed, for
the clusters with R<sub>500</sub> < 6', the normalization in the joint fit is
larger by 1.4% (mean, while 3.7% is the median difference). We
can thus conclude that with the joint fit we overestimate the nor-
malization by 1–4%, which is a much smaller error than the 17%
due to the use of XMM-only data for extended clusters. More-
over this error is smaller than the typical statistical errors on this
parameter (≈4%), over which we marginalize during the MCMC
fit.

We address the impact of the variations of the sky back-
ground parameters on the temperature estimates in Sect. 5.4.1.

4.4. Application to the spectral fit

In the previous subsections, we have described in detail the prop-
erties of our physical background model, how we built and cal-

ibrated it. Here, we focus on how we apply it to the spectra
during the fitting process, both in the Maximum-Likelihood and
in the MCMC approach. As presented in Sect. 3.3, we joint fit
the spectra of the three instruments and to all of them we fit
the cluster component and the sky background model. The lat-
ter is expressed as *apec+phabs*(*apec+powerlaw*), where the
temperature of the Galactic Halo and the normalizations of all
three components are estimated from the fit in the external region
(Sect. 4.3) and appropriately rescaled for the area. All parame-
ters are kept fixed in the standard ML fit, while in the MCMC the
normalizations and the GH temperature are treated as Gaussian
priors, using the statistical errors on the best-fit parameters for
the width (see Table 2, where we summarize the properties of
the priors for the background components in the MCMC analy-
sis).

For the CRPB and the RFC we build a separate model for
each detector. The former is expressed as a combination of two
power laws for the continuum emission and several Gaussians to
represent the fluorescent lines (Sect. 4.1), overall multiplied by a
<table><caption>**Table 2.** Properties of the priors set in the MCMC analysis for the background parameters.</caption><thead><tr><th>Parameter</th><th>Shape</th><th>Central value</th><th>Width</th></tr></thead><tbody><tr><td>CRPB norm. (MOS)</td><td>Gaussian</td><td>Best Fit of <i>mos</i>-back spectra</td><td>2% intrinsic scatter</td></tr><tr><td>CRPB norm. (pn)</td><td>Gaussian</td><td>Best fit of <i>pn</i>-back spectra renormalized by Eq. (2)</td><td>6% intrinsic scatter</td></tr><tr><td>RFC norm. (MOS)</td><td>Uniform</td><td>inFOV-outFOV, with Eq. (A.1)</td><td>Intrinsic scatter</td></tr><tr><td>RFC norm. (pn)</td><td>Uniform</td><td>inFOV<sub>PN</sub> and outFOV<sub>MOS2</sub> with Eqs. (A.1) and (A.2)</td><td>Intrinsic scatter</td></tr><tr><td>LHB norm.</td><td>Gaussian</td><td>Best fit in background region</td><td>1&sigma; errors</td></tr><tr><td>GH temp.</td><td>Gaussian</td><td>Best fit in background region</td><td>1&sigma; errors</td></tr><tr><td>GH norm.</td><td>Gaussian</td><td>Best fit in background region</td><td>1&sigma; errors</td></tr><tr><td>CXB norm.</td><td>Gaussian</td><td>Best fit in background region</td><td>1&sigma; errors</td></tr></tbody></table>

**Notes.** For each of them, we report the shape of the prior, how we computed its central value, and its width.

constant which sets the intensity of this component. This model is first fitted to the spectra of FWC data produced by the ESAS tools *mos/pn-back* and then applied, with its best fit parameters, to the spectra extracted from the observation. When performing the ML fit all parameters are kept fixed, while in the MCMC fit the multiplicative constant is allowed to vary within its uncertainty and treated as a Gaussian prior (Table 2). For the pn detector, the procedure is the same but we renormalize the overall constant according to Eq. (2). Concerning the Residual Focused component, we model it with a power law with fixed slope -0.6 and normalization predicted by the inFOV-outFOV value, measured for each detector, with Eq. (A.1). For the pn, we estimate it by combining the inFOV<sub>PN</sub> and outFOV<sub>MOS2</sub> as described in Eq. (A.2). While doing the MCMC fit, we treat the normalization as a flat prior, within the intrinsic scatter of the relation between the Normalization of the SP component and the inFOV-outFOV indicator (Figs. A.2 and A.3).

## 5. Impact of spectral analysis on temperature measurements

In this section, we evaluate the robustness of our pipeline and the impact of our decisions about the spectral fitting methods on the temperature measurements. While some of these decisions (such as the spectral model, the assumed $N_H$, and the techniques for the spectral fitting) affect all our measurements, other assumptions on the background model will impact mainly the outer cluster regions. Indeed, the effects of the latter strongly depend on the relative intensity of the source with respect to the background. Briefly, even a large error on the background model will not affect significantly the best fit parameters in a region where the source outshines the background, such as in the central regions of galaxy clusters. Conversely, in the external regions where the source intensity is comparable to, or often even smaller than, the background level, even a few percent error in the estimates of the background may affect significantly our temperature measurements. To quantify this, following Leccardi & Molendi (2008), we define the indicator SOU/BKG for all DR1 clusters and for all regions where we extracted spectra. We compute it as (OBS-BKG)/BKG, where OBS is the observed count rate and BKG is the predicted count-rate by the best fit model of all the background components. We measure all count rates in the energy band 0.7–10 keV ( comparable to the band we used in our spectral fit, see Sect. 3.3), to be consistent with Leccardi & Molendi (2008). We are aware that the cluster count rate in the 0.7–10 keV band depends on temperature, inducing a slight dependence also on the SOU/BKG indicator. We estimated its impact by computing our indica-

tor also in a soft band (0.7–1.2 keV), where the count rate is mostly independent of temperature for $kT > 2$ keV. We checked that all the results presented later in this paper (e.g., Figs. 13 and 17) are robust with respect to the choice of the band in which we compute our indicator and we decided to keep using the full spectral band, to allow quantitative comparisons with Leccardi & Molendi (2008). We compute this indicator for each camera but we refer to the values obtained with MOS2.

In Fig. 5 we show the distribution of the SOU/BKG indicator as a function of radius for all the regions, spanning four orders of magnitudes. While most regions are clearly source dominated, we have a significant fraction (~20%) of regions where the background intensity is larger than the source, justifying our need for the most accurate background modeling (Sect. 4). As we subsequently discuss in Sect. 6.2, we trust our temperature measurements for regions with SOU/BKG > 0.2 and we thus excluded about 30 regions with SOU/BKG below this threshold in some of the analysis presented later. In Fig. 5, we note that about half of the regions below the threshold are located at $R < R_{500}$, hampering the goal of reliably measuring the temperature profile up to $R_{500}$ for some clusters. As highlighted by the color-coding, most of them are located in extended clusters ($R_{500} > 9'$), which by construction of our sample correspond to low-redshift and low-mass objects. These regions are located at relatively large off-axis angles from the aimpoint, where the exposure is lower than in the center, suppressing the source count-rate more than the background one, whose dominant CRPB component is not affected by vignetting.

### 5.1. MCMC vs. maximum-likelihood fitting

As discussed in Sect. 3.3, to assess the MCMC framework we have utilized for our analysis, which marginalizes over uncertainties in the background components, we also perform a second set of fits using the standard ML approach, where the background model is fixed. It is possible in principle that the best fit parameters obtained with the two methods could differ from one another, since their definition is different. While the best fit temperature in the standard ML fit is the value which minimizes the Cash statistic, in the MCMC it is the mean of the posterior distribution. These values could be different, especially in the case of a skewed distribution. Moreover, during the marginalization process of the MCMC the background parameters can move with respect to the fixed values in the ML fit, and this can lead to a change in the best fit temperature, especially in background-dominated regions. We expect also an effect on the error bars of the best fit parameters, which should be larger in the MCMC fit, encompassing also the uncertainties on the background modeling.
<figure><img src="image_6.png" /></figure>

**Fig. 5.** Distribution of the SOU/BKG indicator for MOS2 as a function of scaled radius, for all the 354 regions from which we extracted spectra. The color coding refers to the extension in the sky ($R_{500}$ in arcmin) of the corresponding cluster. The horizontal continuous line shows the value where the source intensity gets lower than the background line, while the dashed line marks the threshold SOU/BKG = 0.2 used in the analysis of mean and median profiles (Sect. 6.2). The vertical line marks $R_{500}$. Most of the regions below the threshold SOU/BKG = 0.2 at $R < R_{500}$ belong to very extended clusters and are thus located in poorly exposed regions of the FOV.

In the left panel of Fig. 6, we compare the best fit temperature values obtained with these two methods: the agreement is very good, with a weighted mean of the difference $(T_{\text{ML}} - T_{\text{MCMC}})/T_{\text{MCMC}} = (-0.7 \pm 0.3)\%$. We expect the MCMC approach to have a small impact in source dominated regions and indeed the difference is consistent with zero, $(T_{\text{ML}} - T_{\text{MCMC}})/T_{\text{MCMC}} = (-0.4 \pm 0.4)\%$ when SOU/BKG > 10. Conversely, in background dominated regions, the MCMC approach has a larger impact on the best-fit parameters, since it allows the background parameters to move from their value. Hence, temperatures in this regime are very sensitive to the background models. Indeed, the differences increase to ($-4.6 \pm 2.5)\%$ in regions with $0.2 < \text{SOU/BKG} < 0.5$ and to ($-7.5 \pm 3.5)\%$ for SOU/BKG < 0.2. If we split the regions in three radial bins, we find that the difference is consistent with zero at less than $1-2\sigma$ for $R < 0.4 R_{500}$ and $0.4-0.8 R_{500}$ and measure a difference of $-5.6 \pm 2.7\%$ in the last radial bin ($R > 0.8 R_{500}$). While the results are only moderately significant, we can conclude that the standard ML approach leads to an underestimate of the temperature values of about 3–7% in background-dominated regions. We note that this comparison depends on the estimator used in the MCMC fit to derive the best fit parameters. As described in Sect. 3.3, we use the mean of the posterior distribution, but we could use also the median or the mode. These values do not differ significantly from the mean on the full sample ($\le 1\%$) but can vary by a few percent in background-dominated regions, where the posterior distributions become more skewed, moving closer to the ML estimate. We should thus consider the value provided above as an upper limit for this effect.

In the right panel of Fig. 6, we compare the 68% relative symmetric error on the temperature measurement with the ML and MCMC approach. Again, in source-dominated regions the error on the measurement is dominated by the statistical quality of the data and is very similar in both approaches. Conversely, in background-dominated regions the error on the ML reflects only the statistical quality on the data, while the one on the MCMC should be larger because it propagates also the uncer-

tainties on the background models. This is indeed observed in Fig. 6, where most of the points with low SOU/BKG lie above the equality line. On average, the MCMC error is 25% (31%) larger in regions with SOU/BKG < 1 (SOU/BKG < 0.5) with respect to the standard ML approach. Conversely, the few points with very low SOU/BKG that lie significantly below the line may indicate that the MCMC did not converge there. We individually investigated the three most deviant outliers: they correspond to radial bins whose center is beyond $R_{500}$ and with a very low SOU/BKG < 0.1, corresponding to regions that we later excluded from our analysis (see Sect. 6.2). The visual inspection of the chains confirms that the MCMC had not converged in these cases and we thus reran the MCMC fitting, increasing the number of steps to $10^5$. After the rerun, the error bars increased significantly. We decided not to rerun all our fits with longer MCMC chains and that $5 \times 10^4$ steps is a good compromise between convergence in most cases and computation speed.

## 5.2. Impact of total $N_H$

As presented in Sect. 3.3, we use the total $N_{\text{Htot}}$ value, comprising the contribution of both the atomic and molecular hydrogen, computed in Bourdin et al. (2023). We tested the impact of this choice on the temperature measurements by refitting all spectra with the standard atomic absorption ($N_{\text{HI}}$) and computing the relative difference of the best fit values obtained with and without the molecular component: $\Delta_{\text{TNH}} = (T_{\text{NHI}} - T_{\text{NHtot}})/T_{\text{NHtot}}$. We then divided our sample into three subsamples: a low absorption regime ($N_{\text{HI}} < 2 \times 10^{20} \text{ cm}^{-2}$), an intermediate regime ($2 < N_{\text{HI}}/10^{20} \text{ cm}^{-2} < 5$), and a (relatively) high absorption sample ($N_{\text{HI}} > 5 \times 10^{20} \text{ cm}^{-2}$). We then computed the median value of $\Delta_{\text{TNH}}$ on the three subsamples. In the low absorption regime differences in the temperature values are negligible and the median difference is only 0.2%. In the intermediate regime, we start to see some differences, especially at high temperature, but the median difference is still 1%. Finally, in the high absorption sample we systematically find higher temperatures when fitting with the atomic component only, with a median difference of 6.8%. These results are consistent with those for the full CHEXMATE sample described in Bourdin et al. (2023).

We use the mean difference in the derived temperatures as an estimate of the systematic associated to the assumption of a fixed $N_H$ in our fit. On the full DR1 sample, we estimated $\Delta_{\text{TNH}} = (1.5 \pm 0.3)\%$ and found a consistent value at low SOU/BKG and in different radial ranges. This is indeed expected since it should not depend on the background level. We emphasize that $N_H$ is a fixed parameter in all our fits and thus the uncertainty on its measurement is not propagated in the temperature estimate even in the MCMC fit.

## 5.3. Impact of abundance table

Abundance tables are used in XSPEC to compute plasma emission and photoelectric emission models and thus can impact all the free parameters in our fits, not only the metal abundance. The relative changes can be estimated making use of XSPEC simulations, as described in Appendix B. Here we briefly recall the results concerning the temperature, comparing the ones measured with a given table ($T_{\text{ab}}$) with respect to the values obtained with the reference table in our analysis ($T_{\text{aspl}}$, Asplund et al. 2009, see Sect. 3.3). The abundance table is used both in the apec model and in the absorption phabs model, and it is the latter which causes the most significant differences in the temperature measurements. Indeed, all the simulations performed
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Comparison of best-fit temperature values (left panel) and relative symmetric errors (right panel) for all clusters and all regions obtained with the standard ML technique (x-axis) and with the MCMC approach (y-axis). The color coding refers to the SOU/BKG ratio in each region, the dashed line represents the equality line.</figcaption></figure>

with $N_H = 0$ show negligible variations (<1%) in the temperature estimate. Conversely, when we use an intermediate absorption, $N_H = 5 \times 10^{20} \text{ cm}^{-2}$, we find significant variations ($T_{ab} - T_{aspl})/T_{ab}$ of -3 ± 1% using Lodders et al. (2009) and -6 ± 3% using Anders & Grevesse (1989). We interpret this result as being due to the different abundance of helium in the solar tables.

## 5.4. Impact of background assumption in low SOU/BKG regions

In this section, we assess the impact of different background treatments on the temperature measurements, especially in outer regions associated with a low SOU/BKG ratio. Since we want to address the impact of background parameters that are free to vary within their priors in the MCMC fit, in this section we show results obtained with the ML fit, where these parameters are kept fixed.

### 5.4.1. Impact of joint fit with RASS

In Sect. 4.3 we introduced our strategy to derive the parameters of the sky background components by performing a joint fit of the *XMM-Newton* spectrum beyond $R_{200}$ with the RASS spectrum around the cluster and discussed the impact on the background parameters themselves. Here we want to assess the impact on temperature measurements, by comparing the best fit values obtained with our pipeline ($T_{RASS}$) with those obtained turning off the joint fit with RASS data ($T_{NORASS}$), using XMM data only to estimate the sky background parameters. We compare the temperature ratios as a function of SOU/BKG in Fig. 7. We emphasize that we use the nomenclature $T_{RASS}$, but we never use RASS data for measuring the cluster temperature, only for determining the sky background model in a joint fit with *XMM-Newton*. The differences in the temperature that we discuss here are the consequence of the different sky background parameters, presented in Sect. 4.3 and Fig. 4.

First, we focus on the six clusters of our sample with $R_{500} > 9'$, for which we do not have a clean region in the FOV and we used the “closebkg” annulus (1.1 $R_{500} - 14'$), likely contaminated by cluster emission, as a background region (Sect. 3.2). In Sect. 4.3 we noticed that the introduction of the joint fit with RASS reduces the CXB normalization by a factor of 17%

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Ratio of the difference between the temperature values obtained using XMM only in the background region ($T_{NORASS}$) and that obtained with our standard pipeline ($T_{RASS}$, joint fit with RASS in the background region) over $T_{RASS}$, in percent units, as a function of the SOU/BKG ratio in the spectral region. Filled points show the individual measurements (error bars not shown for clarity), while the open squares represent the weighted mean ratio and its error in different intervals of SOU/BKG. Blue points refer to measurements in clusters with $R_{500} > 9'$, red are for the remaining objects. The cyan triangle shows the median value for the regions with SOU/BKG < 0.2 and clusters with $R_{500} > 9'$, significantly different than the mean with the same criteria.</figcaption></figure>

for these clusters. Concerning the ICM temperature, the effect of a reduced CXB normalization is negligible in regions with SOU/BKG > 1: temperatures obtained with the two methods are always consistent within the error bars and the mean difference is less than 1%. In background-dominated regions (SOU/BKG < 1), individual measurements may be significantly different, with variations up to 20–50%, larger than the statistical errors on the measurements. Nonetheless, there is not a clear trend (see blue points in Fig. 7): in most cases we find $T_{RASS} > T_{NORASS}$ but in others we find the opposite. We computed the weighted mean variations in different SOU/BKG regimes and we find a significant difference only in the bin with SOU/BKG < 0.2, where we find $(T_{NORASS} - T_{RASS})/T_{RASS} = (-22 \pm 4.5)\%$ (see Fig. 7). We
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Ratio of the difference between the temperature values obtained neglecting the soft proton component ($T_{\text{NOSP}}$) and the one of our pipeline ($T_{\text{SP}}$, including a SP component) over $T_{\text{SP}}$ in percent units, as a function of the SOU/BKG ratio in the spectral region. Filled points show the individual measurements (error bars not shown for clarity), while the open squares represent the weighted mean ratio and its error in different intervals of SOU/BKG. Blue points refer to measurements in observations with a significant residual contamination (inFOV-outFOV > 0.04 cts s⁻¹), red are for observations with a low contamination. The magenta triangles show the median value for clusters with low contamination, when the median is significantly different from the mean.</figcaption></figure>

note however that this mean value is driven by one single measurement with small uncertainty but with a variation of ≈50% and may return an overestimate for this effect. The median value (also shown in Fig. 7) on the same sample is significantly lower (−4.6±4.7%) and may better reflect the properties of the distribution. We note that, by construction of the CHEX-MATE sample, the clusters with $R_{500} > 9'$ are nearby low-mass objects, where the temperature in the external regions lies in the range 2–3 keV. It is possible that the variations of the sky-background parameters, especially the CXB normalization, may have a larger impact in different temperature regimes.

We also checked that the joint fit with RASS did not introduce a bias in our temperature measurements, possibly due to the cross-calibration of the two instruments and to the different parameters of the foreground components, with the systematic underestimate of 1–4% in the CXB normalization discussed in Sect. 4.3. To assess this, we compare $T_{\text{RASS}}$ with $T_{\text{NORASS}}$ also for the remaining 22 clusters with $R_{500} < 9'$ (red points in Fig. 7): in all cases temperatures are consistent within the two methods and start to differ only at SOU/BKG < 1, but remain consistent within the error bars in most cases. We computed the mean weighted ratio to be $(T_{\text{NORASS}} - T_{\text{RASS}})/T_{\text{RASS}} = (0.2 \pm 0.4)%$ on the full sample. The difference increases in regions with low SOU/BKG (3.5±3.3% in the range 0.2–0.5 and 2.3±4.7% when SOU/BKG < 0.2), remaining nonetheless consistent with zero. We can conclude that if we introduce a systematic bias with our method, this is less than 5%.

### 5.4.2. Impact of the residual focused component

As discussed in Sect. 4.2, we introduced a novel approach to model the residual focused component basing on the IN-

OUT indicator. This component, which we attribute to residual soft protons, is typically subdominant with respect to other background components in all DR1 observations ((IN-OUT)$_{\text{M2}}$ < 0.1 cts s⁻¹), and becomes comparable to the intensity of the sky background and to the CRPB only in the cluster with the most contaminated observation (PSZ2 G195.75–24.32). Still, it can have an important impact in background dominated regions. To estimate this effect, we divide this sample in two subsamples: 24 objects with low residual contamination (lowSP with (inFOV-outFOV)$_{\text{M2}}$ < 0.04 cts s⁻¹) and six with a significant residual contamination (highSP with 0.04 < (inFOV-outFOV)$_{\text{M2}}$ < 0.1 cts s⁻¹) and we perform the fit turning off the SP correction. We then compare the temperatures obtained with this method with those obtained with our baseline pipeline and we show their ratio in Fig. 8 as a function of the SOU/BKG in the fitting region. Individual measurements are typically consistent in source-dominated regions and start to show large variations (up to more than 100%) when SOU/BKG < 1. This effect is particularly large in the subsample with large residual contamination, where we measure a weighted mean $(T_{\text{NOSP}} - T_{\text{SP}})/T_{\text{SP}} = (14 \pm 6)%$ in regions with 0.2 < SOU/BKG < 1 and (23 ± 11)% when SOU/BKG < 0.2. Nonetheless, the effect is visible also on the individual measurements of clusters with a low contamination. While the weighted mean value for the highSP sample does not seem to report significant differences (Fig. 8), probably driven low by a few measurements with a small statistical error, the median difference (magenta square in Fig. 8) in the SOU/BKG < 0.2 range is $16 \pm 5%$.

The large impact of including or neglecting the SP component on the temperature measurements in external regions of galaxy clusters highlights the importance of modeling this component, even in “clean” observations. We emphasize here however that with this comparison we are not measuring the real systematics associated to our SP modeling but rather measuring the potential effect if we neglect this component in our analysis, which we do not.

An alternative method of measuring the systematic effect associated with the SP modeling is to compare the median scaled temperature profiles of the subsample with high contamination and those of the subsample with low contamination. In principle the profiles should not differ and any systematic difference could be interpreted as due to the SP modeling. We show their difference $((T_{\text{s,highSP}} - T_{\text{s,lowSP}})/T_{\text{s,lowSP}}$, where the notation $T_s$ marks the scaled temperature $T/T_{[0.15-0.75]R_{500}}$) as a function of the scaled radius $R/R_{500}$ in Fig. 9. In some radial bins, we measure a significant difference, which is, however, consistent with the observed dispersion of the profiles, preventing its interpretation as a systematic bias. We expect to derive better estimates when we analyze the full CHEX-MATE sample. For the moment, we can use the dispersion level (10% outside the core) as an upper limit to the systematics related to SP modeling.

### 5.5. Error budget

In the previous sections, we have quantified the impact of our choices for spectral analysis on the temperature estimate. In some cases, such as the choice of $N_H$ and the abundance table, we are able to measure a significant effect which allows us to associate a possible systematic error in our measurement. In other cases, the mean differences between the temperature under different assumptions are consistent with zero at one or two $\sigma$, in different SOU/BKG or radial ranges. This means that the true systematic errors in our method are smaller than the sta-
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Ratio of the difference between the median scaled temperature profiles for the subsample with high contamination with respect to the low contamination sample. The points show the ratio between the median value in each bin, while the dashed lines mark the observed dispersion, which is typically larger than the measured difference, preventing its interpretation as a systematic bias.</figcaption></figure>

<table><caption>**Table 3.** Potential systematic errors in temperature measurements in percent units.</caption><thead><tr><th rowspan="2">Effect</th><th colspan="2">Percent impact in SOU/BKG range</th></tr><tr><th>0.2-0.5</th><th>&lt;0.2</th></tr></thead><tbody><tr><td>N<sub>H</sub></td><td>1.5%</td><td>1.5%</td></tr><tr><td>Ab. table</td><td>3-6%</td><td>3-6%</td></tr><tr><td>ML-MCMC</td><td>&lt;5%</td><td>&lt;7%</td></tr><tr><td>RASS</td><td>&lt;6.6%</td><td>&lt;9.4%</td></tr><tr><td>SP<sup>(*)</sup></td><td>&lt;10%</td><td>&lt;10%</td></tr></tbody></table>

**Notes.** Upper limits are at 95% c.l. (*)Using the difference between the scaled profile in the high and low contamination subsamples.

Statistical uncertainty for our sample of 30 clusters and we can only provide an upper limit for each of them. We estimated the upper limit at 95% confidence level as twice the uncertainty on the weighted mean temperature difference and we report all of them in Table 3. As discussed before, in some cases, related to the background modeling, our estimate of systematics depends on the SOU/BKG of the regions we are considering, with larger effects in more background dominated regions. Conversely, in other cases the effect is independent of the intensity of the source. In Table 3, we report the values in two ranges, for SOU/BKG < 0.2 (where, as we subsequently show in Sect. 6.2, temperature measurements are probably biased) and 0.2 < SOU/BKG < 0.5, that is, the most background dominated regions where we consider we have reliable temperatures. In both regimes, possible systematics related to the choices of our spectral analysis are on the order of a few percent, with upper limit on the order of 10% in the most background dominated regions. We note that this is the same order of magnitude of the cross-calibration between the three EPIC cameras, leading to significantly different temperatures when estimated with single detectors (e.g., Schellenberger et al. 2015; Nevalainen & Molendi 2023). We checked that for the DR1 sample, the mean difference between the value measured with a single detector and our baseline joint fit is 5%, consistent with literature results.

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Cumulative distribution of clusters where the temperature measurement in the [0.8–1.2] R<sub>500</sub> bin reaches a nominal error as a function of the relative statistical error σ<sub>T</sub>/T. The black continuous line shows the goal in the CHEX-MATE feasibility (15%), while the dashed line shows the median value on the sample (12.7%).</figcaption></figure>

Besides the estimate of systematic errors induced by our method, it is also important to assess the statistical quality of our measurements. In Fig. 6 we show the distribution of relative errors in the MCMC and ML fit: the bulk of the values is below 15%, but in a few cases errors can reach values of 30–40%. The size of the errors depends on several factors, such as the signal in the radial bin, the SOU/BKG ratio and even the temperature value itself (with low temperature values featuring smaller errors, even in relative terms, than higher temperatures). In the CHEX-MATE XMM-Newton proposal preparation, we tuned the exposure time to reach an error of 15% in the annulus [0.8–1.2] R<sub>500</sub>, based on an empirical relation between the counts in the soft band and the error on temperature (see CHEX-MATE Collaboration 2021, for more details). We can use our sample to verify this method a posteriori: we extracted spectra in this annulus for each DR1 clusters and fitted them with our pipeline. We used the ML fit, since we did not foresee using the MCMC approach at the time of the proposal. In Fig. 10 we show the cumulative distribution of clusters with a symmetric relative temperature error smaller than a given σ<sub>T</sub>/T. We reach the requirement of 15% error for 22 clusters and the median error on the full sample is 12.7%. The maximum fractional error in this radial range is 23%. We do not find any clear correlation between the size of the errors and the physical properties of the clusters.

# 6. Results

## 6.1. The temperature profiles

We show the temperature profiles of the 30 DR1 clusters in Fig. 11. The measured temperatures range from 2 to more than 14 keV, reflecting the large spread in mass of the CHEX-MATE and DR1 sample. The radial range of the profiles also reflects the presence of some nearby and very extended objects and of high-redshift and compact clusters in the sample. The color coding by total mass in the left panel of Fig. 11 highlights these dependencies.

We rescaled the profiles according to self-similar predictions in the right panel of Fig. 11. We divide the radial coordinate by
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Individual temperature profiles of the 30 DR1 clusters, color-coded by the total mass. The left panel shows the temperature as a function of the projected position of the annuli in arcminutes, highlighting the large range of measured temperatures and extension of the DR1 sample. The right panel shows the rescaled profiles.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Rescaled temperature profiles of the DR1 sample: clusters belonging to the most relaxed, the most disturbed, and the mixed class are shown in blue, red, and gray, respectively. We use the logarithmic scale on the x-axis to focus on the more central regions.</figcaption></figure>

$R_{500}$, estimated from the Planck SZ mass, reported in Table 1. We underline that these masses can be biased with respect to total masses and may scatter around the hydrostatic masses that would be derived from the X-ray data. As discussed in Bartalucci et al. (2023), this may have an impact in the shape of the profiles and should be taken into account when comparing with simulations and with other samples. Concerning the y-axis, we rescaled the temperature by a mean value obtained by fitting a spectrum extracted in an annulus in the range [0.15–0.75] $R_{500}$, that is, excluding the core and the outer background-dominated regions. This scaling is not completely independent from the data in the profiles as it is extracted from the same data set. An alternative

approach (which we use for the scatter profile in Sect. 6.3 and for some literature comparison in Sect. 7.1) is to rescale by an external value derived from the mass, through scaling relations. The rescaling of the profiles reduces significantly the scatter in the profiles and the dependency on the cluster mass, as highlighted in the comparison of the two plots in Fig. 11, both color-coded by cluster mass. Nonetheless, the scatter among the profiles is non-negligible, reflecting both the uncertainties on the measurements and the cluster to cluster variations (Sect. 6.3). The temperature profiles in the central regions of galaxy clusters have been shown to correlate with the dynamical state of the cluster (e.g., Hudson et al. 2010; Leccardi et al. 2010; Pratt et al. 2010). In Campitiello et al. (2022), we identified the most relaxed and the most disturbed clusters in the CHEX-MATE sample, according to the *M* parameter, which is a combination of several dynamical indicators (surface brightness concentration, centroid shift, and power ratios). Bartalucci et al. (2023) showed that this classification segregates clusters with peaked emission measure profiles in the center from those showing a flatter slope. In Fig. 12, we show the temperature profiles of the most relaxed, most disturbed and mixed clusters in the DR1 sample. We do not find a clear segregation of clusters, contrary to the emission measure profile, suggesting that our classification based on the *M* parameter is not effective in separating clusters with different temperature profiles. This may be due to the definition of the parameter itself, which combines indicators on different scales and is not tailored to measure the properties of the core. However, we should also note that we have only three clusters in DR1 belonging to the most relaxed class according to Campitiello et al. (2022) and none of them show very low values of the *M* parameter. It is also, thus, possible that we are lacking the most extreme clusters which would drive the segregation, and this result should be verified in the full CHEX-MATE sample. Indeed, the KS test, used to check the consistency with the parent sample (Sect. 2), is not very sensitive to the tails of the distribution.
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Comparison of weighted mean (left panel) and median (right panel) temperature radial profiles applying different thresholds in the SOU/BKG indicator. We plot only values with more than 15 (half sample) individual measurements.</figcaption></figure>

## 6.2. The mean and median temperature profile

Several methods can be applied to stack the rescaled temperature profiles of individual clusters (Sect. 6.1) into a single average profile of the sample under study. For instance, Leccardi & Molendi (2008) report a weighted mean profile, obtained by weighting the $T_i$ individual measurements by their absolute errors $dT_i$. They apply a further weight by considering the fraction of the original radial bin $r_i \pm dr_i$ which falls in the radial bin $R \pm dR$ used in the output average profile. This method has the advantage to account for the position and the extent of the original radial binning of the individual measurements with respect to the output mean profiles. Moreover, it considers the different errors of individual measurements, giving more weight to more precise estimates. However, it is important to note that the error bars on temperature measurements do not reflect only the statistical quality of the spectra from which they were derived, but also depend strongly on the temperature itself. Indeed, it is easier to measure low temperatures than high temperatures with XMM-Newton, mainly because of the shape of the effective area and the higher intensity of the source with respect to the particle background at low energy. This is reflected into smaller uncertainties (also in relative terms) for low temperature values than for higher temperature measurements. It is easy to infer that the error weighting can bias low the average values in situations where we have a large spread in the temperature measurements, such as in the innermost bins, due to the presence of cool cores, and also in the external regions where measurements in background-dominated regions can lead to (biased) low-temperature values (Leccardi & Molendi 2008). An alternative approach used by other authors (e.g., Arnaud et al. 2010; Lovisari & Reiprich 2019) is to compute the median of the individual measurements at a given radius rather than the weighted mean (see the comparison in Lovisari & Reiprich 2019 for the metal abundance profile). By construction, the median is less sensitive to the tails of the distribution and therefore should be less biased in the situations described above with low temperature values. Conversely, as it is not weighted, we do not account for the position of the initial radial binning with respect to the one of the final profiles. In our analysis, we estimate the error bars on the median profile by performing a randomization of the individual $T_i$ measurements

in each bin within their error bars and recompute the median 1000 times.

In Fig. 13, we show the average profile, both with the weighted mean (taking into account the radial bin fraction) and with the median. In both cases, following Leccardi & Molendi (2008), we show the profiles obtained considering only temperature measurements in regions with the SOU/BKG indicator above different thresholds (namely 2, 1, 0.5, 0.3, 0.2, 0.1 and no threshold), plotting only the points where the average is obtained with more than 15 independent measurements (half sample). Focusing on the weighted mean profile (left panel), we note that profiles with different thresholds are almost identical up to 0.5 $R_{500}$ and start to show some differences only at larger radii. While up to 0.8 $R_{500}$ differences are not significant, in the radial bin between 0.8 and $1R_{500}$, the mean values of the profile without applying any threshold and the one with SOU/BKG > 0.1 (magenta and purple in the left panel of Fig. 13) are significantly lower (at more than $3\sigma$) than the other values. This means that the mean profile is likely biased low, by a factor of about 15%. Conversely, the results obtained with higher thresholds do not deviate from one another, suggesting that a threshold as low as SOU/BKG = 0.2 can be safely used to measure temperature profiles up to $R_{500}$, which is one of the goals of the CHEX-MATE project.

We underline that the value we consider here as a safe threshold for our measurements, SOU/BKG > 0.2, is three times lower than the safe value suggested by Leccardi & Molendi (2008), who introduced this diagnostic plot and also used the weighted mean for their profiles. This implies that the CHEX-MATE pipeline, allows to obtain reliable temperature profiles at much lower SOU/BKG limit, where the cluster intensity is only 20% of the background. We investigated if this improvement could be due to the MCMC fitting approach used here, but results obtained with the standard ML fitting do not differ significantly from what is shown in Fig. 13. The improvement is most likely due to the more advanced physical background model introduced in the CHEX-MATE pipeline and to the high quality of our observations, tailored to reach a good accuracy in temperature measurements at $R_{500}$.

The median profile, shown in the right panel of Fig. 13, does not show the steepening of the profiles obtained with low SOU/BKG thresholds that was apparent in the weighted mean
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Distributions of the rescaled temperatures as a function of the SOU/BKG in the extraction region for the measurements encompassing $0.9 R_{500}$. The vertical dashed lines mark the thresholds used in Fig. 13, with the same color scheme.</figcaption></figure>

profile: except for the last bin, all measurements are very consistent within each other at a given radius. Even in the bin $[0.8-1]R_{500}$ the difference is less than 10% and only appears if we include regions with SOU/BKG < 0.1. This is qualitatively expected since the median is less sensitive to the low temperature values in the tails of the distribution, as it is apparent also in the inner two bins of the profile, where the median profile shows only a moderate decrease, while the weighted mean profile shows a much larger decrease, due to the presence of a few low temperature values, with small error bars, in cool-core clusters.

To understand the effect of including low SOU/BKG measurements in our mean and median profiles we looked at the distribution of the individual measurements $T_i$ and (SOU/BKG)$_i$ values in different radial bins. Up to $0.6 R_{500}$, the distribution of $T_i$ as a function of SOU/BKG is rather flat, therefore the exclusion of a few points at low SOU/BKG does not alter significantly the distribution and its mean and median are not affected by the cut. However, at larger radii, we start to see a correlation between temperature and SOU/BKG, with lower temperatures measured in more background dominated regions. An example is shown in Fig. 14 for the bin at $0.9 R_{500}$: most $T_i/T_{[0.15-0.75]R_{500}}$ values below 0.7 are measured in regions with SOU/BKG < 0.2. Therefore, if we exclude these regions the median changes from 0.72 to 0.79 (~9%) while the weighted mean moves from 0.63 to 0.78 (~18%).

In Fig. 15 we compare the weighted mean, the median and the mean fit profile. The latter is computed as the best-fit mean value in each radial bin, approximating the distribution with a Gaussian function, and is a byproduct of the scatter computation described in Sect. 6.3, to which we refer for details. We show both the profiles obtained with all measurements and those excluding the regions with SOU/BKG < 0.2 (shaded areas). In the central regions, the weighted mean is biased low because of the presence of cool cores, while the mean and median return consistent results. If we consider the values obtained with all measurements (filled circles), the weighted mean tend to be always lower than the other methods and in the last bin at $0.9 R_{500}$ the difference with the median is maximal, with the mean returning an intermediate value. Conversely, out of the

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Comparison of the average profiles obtained with different methods: median (black points and gray shaded area), weighted mean (red), and fit mean (blue). Filled points refer to values derived by all measurements, while shaded regions show the values obtained when applying the SOU/BKG > 0.2 selection. Out of the core, the different methods return consistent results after excluding regions with low SOU/BKG, while if we use all measurements the weighted mean always returns lower values than the other methods.</figcaption></figure>

<table><caption>**Table 4.** Weighted mean and median temperature profiles (rescaled by $T_{[0.15-0.75]R_{500}}$) and their ratio, after excluding regions with SOU/BKG < 0.2.</caption><thead><tr><th>R/R<sub>500</sub></th><th>ΔR/R<sub>500</sub></th><th>Weighted mean</th><th>Median</th><th>Ratio</th></tr></thead><tbody><tr><td>0.025</td><td>0.025</td><td>0.822 ± 0.005</td><td>1.015 ± 0.017</td><td>0.810</td></tr><tr><td>0.075</td><td>0.025</td><td>0.983 ± 0.005</td><td>1.085 ± 0.015</td><td>0.906</td></tr><tr><td>0.125</td><td>0.025</td><td>1.067 ± 0.006</td><td>1.109 ± 0.013</td><td>0.962</td></tr><tr><td>0.175</td><td>0.025</td><td>1.079 ± 0.007</td><td>1.102 ± 0.014</td><td>0.978</td></tr><tr><td>0.225</td><td>0.025</td><td>1.062 ± 0.008</td><td>1.079 ± 0.015</td><td>0.985</td></tr><tr><td>0.275</td><td>0.025</td><td>1.057 ± 0.009</td><td>1.067 ± 0.017</td><td>0.991</td></tr><tr><td>0.350</td><td>0.050</td><td>1.027 ± 0.007</td><td>1.046 ± 0.013</td><td>0.982</td></tr><tr><td>0.450</td><td>0.050</td><td>0.969 ± 0.009</td><td>0.994 ± 0.017</td><td>0.975</td></tr><tr><td>0.550</td><td>0.050</td><td>0.919 ± 0.012</td><td>0.941 ± 0.021</td><td>0.977</td></tr><tr><td>0.700</td><td>0.100</td><td>0.838 ± 0.012</td><td>0.847 ± 0.026</td><td>0.990</td></tr><tr><td>0.900</td><td>0.100</td><td>0.783 ± 0.023</td><td>0.813 ± 0.036</td><td>0.964</td></tr></tbody></table>

core all profiles are consistent when we exclude the regions with SOU/BKG < 0.2. We report the weighted mean and the median profiles under this condition in Table 4). Since the median is the quantity less affected by the presence of low SOU/BKG, from now on we shall use it for the comparison with other samples in the literature (Sect. 7.1), for which it is not possible to apply the SOU/BKG selection.

## 6.3. The intrinsic scatter

In Fig. 11, we qualitatively show that the scatter in the profiles reduces significantly when we rescale temperatures by their mean value in the range $[0.15-0.75]R_{500}$ and radii by $R_{500}$. The residual scatter after scaling is due both to the statistical uncertainties on the individual measurements and to the intrinsic scatter between the profiles, which reflects the variations from cluster to cluster and can shed light on the non self-similar processes affecting the ICM. We computed the radial
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Profile of the intrinsic scatter in the scaled temperature measurements of DR1 clusters. The two profiles refer to the two different methods for scaling temperature: a global temperature measured in the [0.15–0.75] $R_{500}$ annulus (red points) and the temperature $T_{500}$ estimated from the mass (black points).</figcaption></figure>

profile of the intrinsic scatter by modeling the individual scaled $T_i^s = T_i/T_{[0.15–0.75]R_{500}}$ measurements in a given radial bin with a Gaussian, whose width is the sum of the squares of the statistical and intrinsic scatter. We then fit them in a Bayesian framework to derive the parameters of the Gaussian: the mean (used in Sect. 6.2 and Fig. 15) and the intrinsic scatter (see Ghirardini et al. 2019; Bartalucci et al. 2023 for more details). The resulting intrinsic scatter profile for the DR1 temperatures is shown in Fig. 16. We show two profiles, obtained with two different scalings of the temperature values. For the red one, we used the mean temperature estimated in the annulus [0.15–0.75] $R_{500}$, as shown in Sect. 6.1, while for the black one we used an independent estimate of $T_{500}$, derived from the SZ $M_{500}$, using Eq. (10) in Ghirardini et al. (2019). We immediately notice that the scatter profile with the $T_{[0.15–0.75]R_{500}}$ scaling is significantly lower than the one with the external mass scaling in the radial range [0.1–0.6] $R_{500}$. Indeed, the covariance between the measured points and the scaling quantities has been shown to induce suppression of the scatter (Pratt et al. 2022). Conversely, the scaling with the temperature derived from the mass relation is independent of the measured values, since we use the mass derived from the SZ signal, and the scatter profile is indeed larger (see Sect. 7.1). For this reason, from now on, we shall use the mass scaling when measuring the scatter.

The scatter profile shows the typical behavior of thermodynamic quantities, highlighted also in Ghirardini et al. (2019): it shows a high value in the center, reflecting the difference between cool-core and non cool-core clusters, then reaches a minimum between 0.2 and $0.5 R_{500}$ and starts increasing again at large radii. We note that the scatter in the temperature profiles is low, around 10% at its minimum value, smaller than the values observed in the Emission Measure profiles for the full CHEXMATE sample (Bartalucci et al. 2023). Indeed, Ghirardini et al. (2019) showed that temperature is the thermodynamic quantity with the smallest intrinsic scatter.

## 6.4. Fit of the temperature profile

The shape of the temperature profiles of galaxy clusters has been widely discussed in the literature, looking for a uni-

<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Fit of the temperature profile with a power law for $R > 0.3 R_{500}$. Orange (blue) points represent regions where SOU/BKG is greater (smaller) than 0.2. The blue line and dark envelope represent the best fit function and its errors on the full dataset, while the light shaded area marks the intrinsic scatter. The red line and envelope show the best-fit for regions with SOU/BKG > 0.2 (intrinsic scatter not shown for clarity).</figcaption></figure>

versal function (e.g., Baldi et al. 2012; Ghirardini et al. 2019). The most used function has been introduced by Vikhlinin et al. (2006), originally applied to the three-dimensional deprojected profiles, but often used also to model the projected profiles with different choices of the parameters (e.g., Baldi et al. 2012; Ghirardini et al. 2019). This model predicts a temperature gradient in the central regions, which however is typically present only in cool-core clusters, and it is thus customary to fit separately the profiles of the two classes of objects. As discussed in Sect. 6.1 we do not have a clear separation between the two classes in our sample and the number of objects is relatively small. We thus prefer to postpone this analysis to the full CHEXMATE sample. Nonetheless, beyond the core regions, the temperature profiles of cool core and non-cool core clusters have been shown to be quite similar, thus allowing us to fit the entire sample with a single function. Moreover, in this radial regime the profile can be well approximated with a power law with $T/T_{[0.15–0.75]R_{500}} = N(R/R_{500})^{-\alpha}$. We decided to fit the profiles for $R > 0.3 R_{500}$ and derive the best fit parameters for the slope, the normalization and the intrinsic scatter, first using the full dataset and then considering only measurements with SOU/BKG > 0.2. In both cases, we perform our fit in a Bayesian framework, as presented in Sect. 6.3. The points and best fit functions are shown in Fig. 17. For the full dataset, we measure the slope $\alpha = 0.37 \pm 0.04$ and the intrinsic scatter $\sigma_{int} = 0.11 \pm 0.01$. If we exclude the regions with SOU/BKG < 0.2, the slope becomes slightly flatter $\alpha = 0.29 \pm 0.03$ and the scatter reduces to $\sigma_{int} = 0.07 \pm 0.01$. We also split the radial range in two bins ($0.3 < R/R_{500} < 0.55$ and $0.55 < R/R_{500} < 1$) and fit each of them with a separate power law (see e.g., Ghirardini et al. 2019). In the [0.3–0.55] $R_{500}$ radial bin we do not have points with SOU/BKG < 0.2, therefore we performed a single fit and measured $\alpha = 0.25 \pm 0.05$. In the [0.55–1] $R_{500}$, the profile steepens and we find $\alpha = 0.70 \pm 0.15$ for the full sample and $0.44 \pm 0.15$ for SOU/BKG > 0.2 (see also Table 5).

The steepening of the profile while moving to external regions has been already observed in the literature (e.g., Ghirardini et al. 2019) and our results are in qualitative agreement (see also Sect. 7.1). We note however that the steepening
<table><caption>**Table 5.** Best fit slope in different radial range for our sample and comparison samples.</caption><thead><tr><th>Sample</th><th>[0.3–1]</th><th>[0.3–0.55]</th><th>[0.55–1]</th></tr></thead><tbody><tr><td>DR1 all</td><td>0.37 ± 0.04</td><td>0.25 ± 0.05</td><td>0.70 ± 0.15</td></tr><tr><td>DR1 high S/B</td><td>0.29 ± 0.03</td><td>0.25 ± 0.05</td><td>0.44 ± 0.15</td></tr><tr><td>REXCESS</td><td>0.50 ± 0.05</td><td>0.20 ± 0.05</td><td>0.99 ± 0.19</td></tr><tr><td>ESZ</td><td>0.38 ± 0.03</td><td>0.16 ± 0.04</td><td>0.73 ± 0.14</td></tr><tr><td>X-COP spectro</td><td>0.45 ± 0.05</td><td>0.16 ± 0.14</td><td>0.61 ± 0.16</td></tr><tr><td>Simulations</td><td>0.34 ± 0.03</td><td>0.24 ± 0.07</td><td>0.41 ± 0.08</td></tr></tbody></table>

is less pronounced if we consider only measurements with SOU/BKG > 0.2 and that even in the full radial range the best fit profile is flatter than when using all measurements (Fig. 17). This is in agreement with the results already shown for the mean and median profile in Sect. 6.2 and it is likely due to temperature measurements biased low in heavily background-dominated regions (Fig. 14). We checked the results presented in this section also using a SOU/BKG estimator computed in the soft band and we find consistent results.

# 7. Discussion

## 7.1. Comparison with previous results

In this section, we compare our average profiles with results from other samples with different selections. We tried to be consistent in our comparisons, computing the median from the individual profiles for the literature samples with the same code used for DR1 analysis and applying a similar scaling on temperature whenever possible. Since we cannot apply SOU/BKG thresholds in the comparison samples, we do not apply it to the DR1 sample either. We focus only on cluster samples observed by *XMM-Newton* and for which we had access to the individual profiles.

In the left panel of Fig. 18, we compare our median temperature profile with the one of REXCESS (Pratt et al. 2010) and ESZ (a sample of 62 SZ-selected clusters described in Planck Collaboration V 2013). For these samples, we could rescale the temperature by the value measured in an annulus [0.15–0.75] $R_{500}$, as we do for our data, with the only difference that $R_{500}$ in the comparison samples is obtained from $M_{Y_x}$ while in DR1 from $M_{SZ}$. Moreover the binning scheme used to derive the temperature profiles in both samples is very similar to the one we used here. REXCESS is a purely X-ray selected sample (Böhringer et al. 2007), covering a redshift range 0.05–0.18 and a mass range $M_{500}$ in [1–8]×10$^{14}$ $M_{\odot}$ (Pratt et al. 2010). Conversely, ESZ is the intersection of the 189 clusters detected by Planck in the Early SZ sample (Planck Collaboration VIII 2011) and the *XMM-Newton* archive as of mid-2011 (Planck Collaboration XI 2011). The DR1 median profile is consistent, at less than one sigma, with the other measurements in almost all radial bins. We note that the error bars of the DR1 profiles are always smaller than in the other samples although with a similar sample size. This may be related to the better quality of the CHEX-MATE data, tailored to reach a 15% error on the temperature measurements at $R_{500}$, with respect to the other samples. The only significant difference is with REXCESS in the [0.8–1] $R_{500}$ bin, where temperature measurements are more sensitive to the details of the background treatment. Moreover, the different selection as well as the lower median mass in the sample could also have a role in the difference.

In the right panel of Fig. 18 we compare the DR1 median profile with that of the X-COP sample, both from the standard X-ray analysis (blue squares) and with the joint fit with SZ data (red triangles, see Ghirardini et al. 2019 for more details). Here, the background treatment of the *XMM-Newton* data is more similar to what we used, but the binning and observation strategy is different. Moreover, we could not scale the temperature by an external value in the same region, so we relied on rescaling by a predicted temperature from the mass, as introduced in Sect. 6.3, using in both cases Eq. (10) in Ghirardini et al. (2019). This could induce some difference, since for X-COP we used the hydrostatic masses rather than $M_{SZ}$. In the external regions, the agreement between the DR1 profile and the spectroscopic X-COP is very good. Only in the last radial bin, the joint X-SZ measurement is slightly larger than the DR1 profile, but it is interesting to note that it is consistent with our measurement once we exclude low SOU/BKG regions (gray empty circles in Fig. 18). This is further indirect indication of a bias in measuring X-ray temperatures in strongly background dominated spectra. There are some differences in the more internal regions which are at least partly due to the different fraction of cool cores in the two samples and to the different binning.

In Fig. 19, we compare the scatter profile of CHEX-MATE DR1 with the literature samples discussed above. The DR1 scatter profile is consistent with that of REXCESS and ESZ at all radii and only marginally smaller than the one in X-COP between [0.15–0.5] $R_{500}$. We note that in all cases we used $T_{500}$ derived from the mass to rescale the temperature measurements, but only for the DR1 sample the scaling value is independent, as it is derived from the SZ masses. Indeed, for X-COP we used the hydrostatic masses, while for ESZ and REXCESS we used the mass derived from $Y_x$. In both cases they are not independent from the temperature measurements and this may in principle induce a covariance which suppresses the scatter.

We fit the profile with a power law also for the comparison samples in the same radial ranges used for DR1 (Sect. 6.4) and report the best fit slopes in Table 5. For the comparison samples, we always use all available measurements in our fits. The slope values are typically consistent within $2\sigma$ with the values measured for our sample using all our measurements. It is interesting to note that, especially in the radial bin [0.55–1] $R_{500}$, all profiles show a steep slope and only our measurement with SOU/BKG > 0.2 is significantly flatter, which is however consistent with the slope measured from the combination of X-ray and SZ data in X-COP by Ghirardini et al. (2019) in a similar radial range ($\alpha$ = -0.34 ± 0.18 for 0.56 < $R/R_{500}$ < 0.86). This suggests that once we exclude heavily background-contaminated regions, our analysis allows us to reproduce the shape of the temperature profiles as derived by other independent methods, in principle less affected by background systematic than X-ray measurements. During the fit procedure, we scale our profiles, as well as REXCESS and ESZ, with $T_{[0.15–0.75]R_{500}}$, while the X-COP profiles are rescaled based on the mass.

## 7.2. Comparison with numerical simulations

In this section, we compare our results with predictions from numerical cosmological simulations. More specifically, we used the GADGET-X version of the simulations from THE THREE HUNDRED collaboration (Cui et al. 2018), which contains a relatively large set of massive clusters and is thus well-suited to compare with the CHEX-MATE sample. The baryon physics adopted in these simulations includes the star formation models and stellar and AGN feedback, described in Rasia et al. (2015).
<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Comparison of median temperature profiles of the CHEX-MATE DR1 sample (black circles) with other samples studied in the literature. Left panel: comparison with ESZ and REXCESS median profiles, where we have scaled temperature values by the best fit in a region [0.15–0.75] $R_{500}$. Right panel: comparison with X-COP, both the standard results from the X-ray analysis (blue squares) and the joint fit with SZ data (red triangles). The gray empty circles show the DR1 median profile after applying the SOU/BKG > 0.2 selection. Here the temperature scaling $T_{500}$ is computed from the mass with the expression in Ghirardini et al. (2019), both for X-COP and for CHEX-MATE data. Dashed lines show the total scatter in the DR1 profile.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Comparison of the intrinsic scatter profiles in CHEX-MATE DR1 (black points) with other samples in the literature (colored shaded areas), showing consistent results at all radii.</figcaption></figure>

We built a simulated counterpart for the DR1 sample using a procedure similar to the one adopted in Bartalucci et al. (2023), to which we refer for details. Briefly, we match each observed cluster with a simulated one, in the closest redshift snapshot and with similar mass, taking into account a 20% possible underestimate of the mass $M_{SZ}$ of DR1 clusters due to hydrostatic bias. For each simulated cluster, we extract spectroscopic-like temperature (Mazzotta et al. 2004) profiles in fixed bins of $R_{500}$, projected in three orthogonal directions. We compute a temperature and an error in each radial bin for each object, using the mean and dispersion of the three projections.

We then analyze the profiles of simulated clusters as we do for observations (Sect. 6.3), to obtain the mean profile shown in Fig. 20. The mean scaled profile, obtained by rescaling temperatures by the value measured in the annulus [0.15–0.75] $R_{500}$, shows a nice agreement with the observed profile at $R > 0.3 R_{500}$.

This is confirmed by the fit with a power law of the simulated profiles, which returns a slope 0.34 ± 0.03 in the [0.3–1] $R_{500}$ radial range, consistent with the value for the DR1 sample (see Table 5). It is interesting to note that for $R > 0.55 R_{500}$ we obtain a slope of 0.41 ± 0.07, which is similar to the value obtained using only the regions with SOU/BKG > 0.2, thus supporting our choice of excluding the more background-dominated regions. We see some differences in the central regions, which may be due to the different morphological mix between the observed and simulated sample, but also to the difficulties of simulations to fully reproduce the complex physics of cool core clusters.

To compute the profile of intrinsic scatter, also shown in Fig. 20, both for observations and simulations we rescale the temperature by $T_{500}$ derived from the mass to minimize the covariance. For simulations, we reduce the true masses by a factor of 20% to mimic “hydrostatic-like” masses. The scatter in simulated clusters shows a very nice agreement with the observations, suggesting that our simulations and the DR1 sample are consistently sampling the variety of properties of the clusters population. Moreover, it suggests that these simulations capture the main physical processes (gravity-driven accretion and feedback processes) driving the shape of the profiles and their diversity.

# 8. Summary and conclusions

In this paper, we present the pipeline that will be used in the CHEX-MATE collaboration for spectral extraction in radial annuli, for their fit, and for the production of projected radial profiles. We also applied it for the study of the temperature profile of a subsample of 30 objects that are representative of the CHEX-MATE sample.

The pipeline combines the best practices developed during previous projects, namely REXCESS (Pratt et al. 2007, 2010), Planck (Planck Collaboration XI 2011; Planck Collaboration V 2013), X-COP (Tchernin et al. 2016; Ghirardini et al. 2019), and M2C (Bartalucci et al. 2018, 2019), taking the peculiar characteristics of the CHEX-MATE sample (including both nearby extended objects and more distant compact clusters
<figure><img src="image_22.png" /><figcaption>**Fig. 20.** Comparison of the mean temperature (left panel) and intrinsic scatter (right panel) of the DR1 sample with those obtained from numerical simulations, using both spectroscopic-like temperatures.</figcaption></figure>

CHEX-MATE Collaboration 2021) into account. The main novelties of our pipeline, which have been presented in this paper for the first time, are the construction of a physical model for the particle background (both induced by high-energy Galactic cosmic rays and soft protons) for all EPIC detectors and the application of a Bayesian MCMC framework, allowing for the propagation of uncertainties on the background parameters up to the final results of the spectral analysis.

We applied the pipeline to a subset of 30 CHEX-MATE clusters (DR1), which was built to cover the different observation properties of the full sample (background level and soft proton contamination) and, at the same time, to be representative of the physical properties of the parent sample, including the different angular sizes of the targets. We measured temperatures ranging from ~2 to ~14 keV in different clusters and regions and we could map the profiles from the cores to the external regions around $R_{500}$ (Fig. 11). We have shown that, with the CHEX-MATE pipeline and data quality, we are able to recover temperatures in regions where the source intensity is as low as 20% of the background in the 0.7–10 keV band, which is a significant improvement with respect to past analysis (Leccardi & Molendi 2008 considered reliable-only measurements in regions where SOU/BKG > 0.6). We estimated that the choices that we inevitably had to make during the spectral fitting and background modeling impact the temperature measurements for less than 10%, as long as we focus on regions with SOU/BKG > 0.2.

We used our results to measure the weighted mean and median temperature profile and found that the latter is less affected by regions with low SOU/BKG. We thus used the median profile to compare our results with literature samples, where we could not apply the SOU/BKG cut, and found consistent results. We also measured the profile of the intrinsic scatter and showed that scaling by a mean temperature $T_{[0.15-0.75]R_{500}}$ suppresses the scatter with respect to the scaling for an external quantity, such as the gravitational temperature derived from the mass (Fig. 16). With the external scaling, the agreement is very good, both with other samples (Fig. 19) and with simulations (Fig. 20). We modeled the temperature profile at $R > 0.3 R_{500}$ with a piece-wise power law and noticed that the profile steepens going outward, in agreement with other samples in the literature. It is interesting to note that in the radial range $[0.55-1] R_{500}$, the slope of the profile is steep and consistent with all other X-ray measurements in the literature. However, once we excluded

heavily background-contaminated regions with SOU/BKG < 0.2, our analysis returned a flatter profile, reproducing the shape of the temperature profiles derived by the joint fit of SZ and X-ray measurements (Ghirardini et al. 2019), which in principle are less affected by background systematic than X-ray spectroscopic measurements. Moreover, our profile matches the predictions of numerical simulations.

The temperature profiles presented in this paper will be used in combination with the surface brightness measurements shown in Bartalucci et al. (2023) to reconstruct the three-dimensional radial profiles of most thermodynamic quantities of the ICM, the total mass profile through the hydrostatic equilibrium equation, and all of the related integrated quantities, with the purpose of addressing the main scientific question of the CHEX-MATE projects (CHEX-MATE Collaboration 2021). We will apply the CHEX-MATE pipeline presented in this paper for spectral extraction and fitting to the full CHEX-MATE sample of 118 clusters to extract radial temperature profiles for all of them. This will allow us to extend the analysis presented in this paper to the full sample and thus to perform the statistical comparison that was not possible with only 30 objects. The full sample will allow us to compare the temperature profiles in different mass and redshift ranges, as well as in different dynamical states.

**Acknowledgements.** We thank the referee for their timely report and for providing constructive comments. Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. M.R., F.G., S.E., S.M., I.B., H.B., S.D.G., F.D.L., S.G., L.L., P.M., and G.R. acknowledge the financial contribution from the contracts Prin-MUR 2022, supported by Next Generation EU (M4.C2.1.1, n.20227RNLY3, the concordance cosmological model: stress-tests with galaxy clusters), 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), H.B., P.M., and F.D.L. acknowledge also the financial contribution from INFN through the InDark initiative, from Tor Vergata Grant “SUPERMASSIVE-Progetti Ricerca Scientifica di Ateneo 2021”, and from Fondazione ICSC, Spoke 3 Astrophysics and Cosmos Observations, National Recovery and Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR) Project ID CN_00000013 “Italian Research Center on High-Performance Computing, Big Data and Quantum Computing” funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di “campioni nazionali di R&S (M4C2-19)” – Next Generation EU (NGEU). L.L. also acknowledges support from INAF mini grant 1.05.12.04.01. M.S. acknowledges financial contributions from contract ASI-INAF n.2017-14-H.0
and INAF Theory Grant 2023 “Gravitational lensing detection of matter distribution at galaxy cluster boundaries and beyond” (1.05.23.06.17). R.M.B. acknowledges support by Centre National d'Études Spatiales (CNES) under the EPIC-SSC project for *XMM-Newton*. G.W.P. acknowledges support from CNES, the French space agency. J.S. was 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”. We also made use of the following packages: FTOOLS (Blackburn 1995), DS9 (Joye & Mandel 2003), and several python libraries, such as numpy (Harris et al. 2020), matplotlib (Hunter 2007), Astropy (Astropy Collaboration 2013, 2018), PyMC (Abril-Pla et al. 2023).

## References

Abril-Pla, O., Andreani, V., Carroll, C., et al. 2023, PeerJ Comput. Sci., 9, e1516

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

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

Arnaud, M., Pointecouteau, E., & Pratt, G. W. 2005, A&A, 441, 893

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

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

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

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

Baldi, A., Ettori, S., Molendi, S., & Gastaldello, F. 2012, A&A, 545, A41

Bartalucci, I., Arnaud, M., Pratt, G. W., & Le Brun, A. M. C. 2018, A&A, 617, A64

Bartalucci, I., Arnaud, M., Pratt, G. W., Démoclès, J., & Lovisari, L. 2019, A&A, 628, A86

Bartalucci, I., Molendi, S., Rasia, E., et al. 2023, A&A, 674, A179

Basu, K., Zhang, Y. Y., Sommer, M. A., et al. 2010, A&A, 519, A29

Blackburn, J. K. 1995, in Astronomical Data Analysis Software and Systems IV, eds. R. A. Shaw, H. E. Payne, & J. J. E. Hayes, ASP Conf. Ser., 77, 367

Böhringer, H., Schuecker, P., Pratt, G. W., et al. 2007, A&A, 469, 363

Bourdin, H., & Mazzotta, P. 2008, A&A, 479, 307

Bourdin, H., De Luca, F., Mazzotta, P., et al. 2023, A&A, 678, A181

Bulbul, E., Chiu, I. N., Mohr, J. J., et al. 2019, ApJ, 871, 50

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

Cash, W. 1979, ApJ, 228, 939

Chen, C. M. H., Arnaud, M., Pointecouteau, E., Pratt, G. W., & Iqbal, A. 2023, A&A, submitted [arXiv:2311.10397]

CHEX-MATE Collaboration (Arnaud, M., et al.) 2021, A&A, 650, A104

Cui, W., Knebe, A., Yepes, G., et al. 2018, MNRAS, 480, 2898

De Grandi, S., & Molendi, S. 2002, ApJ, 567, 163

Eckert, D., Molendi, S., & Paltani, S. 2011, A&A, 526, A79

Eckert, D., Molendi, S., Vazza, F., Ettori, S., & Paltani, S. 2013, A&A, 551, A22

Ettori, S., & Molendi, S. 2011, Mem. Soc. Astron. It. Suppl., 17, 47

Ettori, S., Ghirardini, V., & Eckert, D. 2020, Astron. Nachr., 341, 210

Fabian, A. C. 1994, in Cosmological Aspects of X-Ray Clusters of Galaxies, ed. W. C. Seitter, NATO Adv. Study Inst. (ASI) Ser. C, 441, 163

Fabian, A. C., Nulsen, P. E. J., & Canizares, C. R. 1984, Nature, 310, 733

Fabian, A. C., Arnaud, K. A., Bautz, M. W., & Tawara, Y. 1994, ApJ, 436, L63

Frank, K. A., Peterson, J. R., Andersson, K., Fabian, A. C., & Sanders, J. S. 2013, ApJ, 764, 46

Gaspari, M., Brighenti, F., Temi, P., & Ettori, S. 2014, ApJ, 783, L10

Gastaldello, F., Marelli, M., Molendi, S., et al. 2022, ApJ, 928, 168

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

Ghirardini, V., Bulbul, E., Kraft, R., et al. 2021, ApJ, 910, 14

Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hudson, D. S., Mittal, R., Reiprich, T. H., et al. 2010, A&A, 513, A37

Humphrey, P. J., Liu, W., & Buote, D. A. 2009, ApJ, 693, 822

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

Irwin, J. A., Bregman, J. N., & Evrard, A. E. 1999, ApJ, 519, 518

Joye, W. A., & Mandel, E. 2003, in Astronomical Data Analysis Software and Systems XII, eds. H. E. Payne, R. I. Jedrzejewski, & R. N. Hook, ASP Conf. Ser., 295, 489

Kaastra, J. S. 2017, A&A, 605, A51

Kaiser, N. 1986, MNRAS, 222, 323

Kuntz, K. D., & Snowden, S. L. 2000, ApJ, 543, 195

Kuntz, K. D., & Snowden, S. L. 2008, A&A, 478, 575

Laganá, T. F., Durret, F., & Lopes, P. A. A. 2019, MNRAS, 484, 2807

Leccardi, A., & Molendi, S. 2007, A&A, 472, 21

Leccardi, A., & Molendi, S. 2008, A&A, 486, 359

Leccardi, A., Rossetti, M., & Molendi, S. 2010, A&A, 510, A82

Lodders, K., Palme, H., & Gail, H. P. 2009, Landolt-Börnstein, 4B, 712

Lovisari, L., & Maughan, B. J. 2022, in Handbook of X-ray and Gamma-ray Astrophysics, eds. C. Bambi, & A. Santangelo, 65

Lovisari, L., & Reiprich, T. H. 2019, MNRAS, 483, 540

Lovisari, L., Ettori, S., Rasia, E., et al. 2024, A&A, 682, A45

Marelli, M., Salvetti, D., Gastaldello, F., et al. 2017, Exp. Astron., 44, 297

Marelli, M., Molendi, S., Rossetti, M., et al. 2021, ApJ, 908, 37

Markevitch, M., Forman, W. R., Sarazin, C. L., & Vikhlinin, A. 1998, ApJ, 503, 77

Markevitch, M., Govoni, F., Brunetti, G., & Jerius, D. 2005, ApJ, 627, 733

Mazzotta, P., Rasia, E., Moscardini, L., & Tormen, G. 2004, MNRAS, 354, 10

McCammon, D., Almy, R., Apodaca, E., et al. 2002, ApJ, 576, 188

Melin, J. B., Bartlett, J. G., & Delabrouille, J. 2006, A&A, 459, 341

Molendi, S., & Pizzolato, F. 2001, ApJ, 560, 194

Moretti, A., Pagani, C., Cusumano, G., et al. 2009, A&A, 493, 501

Nevalainen, J., & Molendi, S. 2023, A&A, 676, A142

Owers, M. S., Nulsen, P. E. J., Couch, W. J., & Markevitch, M. 2009, ApJ, 704, 1349

Peterson, J. R., Paerels, F. B. S., Kaastra, J. S., et al. 2001, A&A, 365, L104

Peterson, J. R., Kahn, S. M., Paerels, F. B. S., et al. 2003, ApJ, 590, 207

Pierre, M., Pacaud, F., Adami, C., et al. 2016, A&A, 592, A1

Planck Collaboration VIII. 2011, A&A, 536, A8

Planck Collaboration XI. 2011, A&A, 536, A11

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

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

Planck Collaboration XXVII. 2016, A&A, 594, A27

Ponti, G., Zheng, X., Locatelli, N., et al. 2023, A&A, 674, A195

Poon, H., Okabe, N., Fukazawa, Y., Akino, D., & Yang, C. 2023, MNRAS, 520, 6001

Pratt, G. W., Böhringer, H., Croston, J. H., et al. 2007, A&A, 461, 71

Pratt, G. W., Arnaud, M., Piffaretti, R., et al. 2010, A&A, 511, A85

Pratt, G. W., Pointecouteau, E., Arnaud, M., & van der Burg, R. F. J. 2016, A&A, 590, L1

Pratt, G. W., Arnaud, M., Maughan, B. J., & Melin, J. B. 2022, A&A, 665, A24

Rasia, E., Borgani, S., Murante, G., et al. 2015, ApJ, 813, L17

Salvetti, D., Marelli, M., Gastaldello, F., et al. 2017, Exp. Astron., 44, 309

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

Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30

Sereno, M., Umetsu, K., Ettori, S., et al. 2018, ApJ, 860, L4

Sereno, M., Umetsu, K., Ettori, S., et al. 2020, MNRAS, 492, 4528

Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91

Snowden, S. L., Freyberg, M. J., Plucinsky, P. P., et al. 1995, ApJ, 454, 643

Snowden, S. L., Egger, R., Freyberg, M. J., et al. 1997, ApJ, 485, 125

Snowden, S. L., Mushotzky, R. F., Kuntz, K. D., & Davis, D. S. 2008, A&A, 478, 615

Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments Astrophys. Space Phys., 4, 173

Tchernin, C., Eckert, D., Ettori, S., et al. 2016, A&A, 595, A42

Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487

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

Walker, S., Simionescu, A., Nagai, D., et al. 2019, Space Sci. Rev., 215, 7

Wang, Q. H. S., Markevitch, M., & Giacintucci, S. 2016, ApJ, 833, 99

White, D. A. 2000, MNRAS, 312, 663

<footnote n="1">¹ INAF, IASF-Milano, Via A. Corti 12, 20133 Milano, Italy e-mail: mariachiara.rossetti@inaf.it</footnote>

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

<footnote n="3">³ INAF, Osservatorio di Trieste, Via Tiepolo 11, 34131 Trieste, Italy</footnote>

<footnote n="4">⁴ IFPU, Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy</footnote>

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

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

<footnote n="7">⁷ INFN, Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy</footnote>

<footnote n="8">⁸ DiSAT, Università degli Studi dell'Insubria, Via Valleggio 11, 22100 Como, Italy</footnote>

<footnote n="9">⁹ Astronomy Unit, Department of Physics, University of Trieste, Via Tiepolo 11, 34131 Trieste, Italy</footnote>

<footnote n="10">¹⁰ Università degli Studi di Roma “Tor Vergata”, Via della Ricerca Scientifica, 1, 00133 Roma, Italy</footnote>
11 INFN, Sezione di Roma 2, Università degli studi di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy

12 INAF, Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807 Merate, Italy

13 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 5, 00185 Roma, Italy

14 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

15 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

16 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The University of Manchester, Alan Turing Building, Manchester M13 9PL, UK

17 HH Wills Physics Laboratory, University of Bristol, Tyndall Ave, Bristol BS8 1TL, UK

18 IRAP, Université de Toulouse, CNRS, CNES, UPS, 9 Av. du Colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France

19 Dipartimento di Fisica, Università degli Studi di Milano, Via G. Celoria 16, 20133 Milano, Italy

20 California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA
# Appendix A: Calibration of the background model

This section presents the calibration of the physically motivated background model introduced in Sect. 4 on a large set of more than 500 individual blank-sky pointings. The selected blank fields were mainly collected in the framework of the XMM-XXL survey (Pierre et al. 2016) and following surveys. The typical duration of the observations is in the range 10–40 ks, which matches well the typical duration of CHEX-MATE observations. The blank-sky pointings were reduced using the exact same procedure as devised in Sect. 3.1, including the extraction of CRPB spectra (Sect. 4.1) and of the SP contamination indicator inFOV-outFOV.

## A.1. Stacked residual spectra

We extracted spectra from four annuli centered on the telescope’s aim point covering the radial ranges [0’–5’], [5’–9’], [9’–12’], and [12’–15’]. For each spectrum, we modeled the CRPB contribution using the method devised in Sect. 4.1 and used the ROSAT all-sky survey spectra to model the sky background (see Sect. 4.3). These two components are well characterized and their origin is well understood, such that the difference between the measured spectra and the predicted model gives us a handle on the residual focused component, at least partly due to SP. We then sorted the observations in term of their residual SP contamination as traced by the inFOV-outFOV indicator and stacked the residual spectra in bins of inFOV-outFOV.

In Fig. A.1 we show the stacked residuals in the [12’–15’] annulus for low (inFOV-outFOV < 0.03 cts/s) and high SP contamination (0.08 < inFOV-outFOV < 0.18 cts/s). We see a significant excess over the model in all cases, although the amplitude of the excess strongly depends on inFOV-outFOV. At low energy (E < 2 keV), where the dominant background component is the sky emission, the residuals are essentially independent of inFOV-outFOV. This points toward a systematic calibration offset between XMM-Newton and ROSAT, as already shown by Eckert et al. (2011). The ROSAT spectral normalization appears on average 14% lower than the XMM-Newton one, which we correct for the remainder of the analysis by applying a renormalization factor of 1/0.88 to the ROSAT spectra (see Sect. 4.3). Conversely, beyond 2 keV the residual spectra strongly depend on inFOV-outFOV, indicating that the model residuals scale well with the SP contamination. Therefore, we fit the residual spectra beyond 2 keV with a power law $I(E) \propto E^{-\Gamma}$, which is indicated with the solid lines in Fig. A.1. The spectral shape of the residual SP component appears universal, with a slope of $\Gamma = 0.6$ that is nearly independent of inFOV-outFOV.

## A.2. Relation between SP normalization and inFOV-outFOV for MOS

Given the finding that the excess component seems to correlate with SP contamination and that its spectral shape is nearly universal, we fit the residual spectra of each observation with a power law $I(E) = N_{QSP}E^{-\Gamma}$ with a fixed slope $\Gamma = 0.6$. We then study how the normalization of the SP component, $N_{QSP}$, can be predicted on the basis of inFOV-outFOV. Assuming that the additional component considered here can be entirely attributed to a residual SP contamination, we expect this component to vary across the detector with a vignetting curve that is consistent with the known SP vignetting, which is flatter than that of photons. To take the SP vignetting into account, for each spectrum we used the SAS task protonscale to determine an effective region area

(hereafter PROTSCAL) and normalized the fitted SP normalization by the corresponding PROTSCAL value, such that in case the excess component is induced by a residual SP component its pattern across the detector is correctly taken into account.

In Fig. A.2 we show the fitted SP normalization scaled by the SP vignetting factor PROTSCAL as a function of inFOV-outFOV for the two MOS detectors. We observe a clear correlation between the two quantities, with a Pearson correlation coefficient of 0.82 (MOS1) and 0.88 (MOS2). For comparison, we also show the normalization of the excess component obtained when fitting stacked spectra in bins of inFOV-outFOV, which allows us to better highlight the relation. We then fit the relation with a power law,

$$ \frac{N_{QSP}}{\text{PROTSCAL}} = A_{\text{SP}} \left( \frac{\text{inFOV} - \text{outFOV}}{0.05 \text{cts/s}} \right)^{B_{\text{SP}}} \quad (\text{A.1}) $$

and a free log-normal intrinsic scatter. For MOS1 we measure $A_{SP,MOS1} = (4.73 \pm 0.29) \times 10^{-5}$ and $B_{SP,MOS1} = 1.15 \pm 0.10$ with an intrinsic scatter $\sigma_{logN} = 0.16^{+0.07}_{-0.05}$ dex. We obtain very similar values for MOS2, with $A_{SP,MOS2} = (4.53 \pm 0.31) \times 10^{-5}$, $B_{SP,MOS2} = 1.06 \pm 0.09$ and an intrinsic scatter $\sigma_{logN} = 0.18^{+0.06}_{-0.05}$ dex. The consistent results obtained with the two MOS detectors are probably unsurprising, given that the detector and telescope properties are very similar.

## A.3. Relation between SP normalization and inFOV-outFOV for pn

As in the previous section, we attempted to correlate the normalization of the excess component with the inFOV-outFOV indicator for the pn detector as well. While we also found a relation between the pn excess component and inFOV-outFOV (computed for the MOS2 detector as in Marelli et al. 2017), the relation appears to have a much larger intrinsic scatter of ~0.4 dex, such that it is much more difficult to predict the value of $N_{QSP}$ for the pn data on the basis of inFOV-outFOV. The main reason for the large intrinsic scatter is the difference in the selection of good time intervals between the MOS and pn detectors. Indeed, given the higher sensitivity of the pn to SP, the light curve filtering procedure (see Sect. 3.1) typically excludes a higher fraction of the total observation time for pn than for MOS, such that the selection of good time intervals need not coincide. Given the highly variable nature of the SP component, the relation between the excess component and the inFOV-outFOV indicator based on MOS2 data exhibits a substantial level of scatter. To alleviate this issue, we considered the possibility of defining its own inFOV-outFOV indicator for pn as well. This is rendered difficult by the lack of shielded area in the pn. Following the procedure devised in Sect. 4.1, we predict the expected intensity of the pn CRPB on the basis of the MOS2 outFOV value, which can be done in spite of the different selection of good time intervals as the CRPB is known to vary over timescales that are typically much longer than our observation. We can then determine the high-energy count rate inside the pn FoV, scale this value per unit area, and compute an equivalent of inFOV-outFOV in the following way:

$$ (\text{inFOV-outFOV})_{\text{PN}} = CR_{\text{ann}} - A_{\text{CRPB}} \text{outFOV}_{\text{MOS2}} \quad (\text{A.2}) $$

with $A_{CRPB}$ the CRPB proportionality constant defined in Eq. 2 and $CR_{ann}$ the [10-14] keV count rate in the [12–15] arcmin annulus centered on the aim point. The use of the external annulus allows us to exclude any potential photon contribution in case a bright, hot cluster is observed on axis. We then attempted
<figure><img src="image_23.png" /></figure>

**Fig. A.1.** Stacked residual spectra in the [12' – 15'] annulus for two sets of blank-sky observations sorted by their values of *inFOV – outFOV*. The left-hand panel shows the low-contamination dataset (*inFOV – outFOV* < 0.03 cts/s) whereas in the right-hand panel we are selecting only highly contaminated observations (0.08 < *inFOV – outFOV* < 0.18 cts/s). The residuals were fit beyond 2 keV with a power law indicated by the solid lines. The bottom panels show the residuals compared to the model.

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

**Fig. A.2.** Normalization $N_{QSP}$ of the excess component as a function of *inFOV – outFOV* for MOS1 (left) and MOS2 (right). The gray points show individual blank-sky observations whereas the red points show the stacked results in bins of *inFOV – outFOV*. The blue solid line and shaded area show the fit to the data and the $1\sigma$ error envelope around the model, whereas the cyan shaded area shows the intrinsic scatter of the relation.

to correlate the new indicator defined in Eq. A.2 with the normalization of the pn excess component. The resulting relation is shown in Fig. A.3. We found a much tighter correlation than with the MOS-based quantity, with a correlation coefficient of 0.85. Fitting the data with the same model as in Eq. A.1 albeit with (*inFOV – outFOV*)PN instead of *inFOV – outFOV*, we found $A_{SP,PN} = (1.85 \pm 0.08) \times 10^{-4}$, $B_{SP,PN} = 0.98 \pm 0.05$, and $\sigma_{\log N} = 0.10_{-0.04}^{+0.05}$ dex. We can see that the normalization of the excess component for pn is about 4 times higher than that of the MOS detectors. In all cases, the slope of the relation is close to one.

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

**Fig. A.3.** Normalization $N_{QSP}$ of the excess component as a function of the pn *inFOV – outFOV* indicator, computed as in Eq. A.2. As in the previous plot, gray points represent individual measurements and red points are stacked results, the blue solid line and area show the best fit relation and its $1\sigma$ envelope, while the light blue shows the intrinsic scatter.
## Appendix B: Effects of the different abundance tables on the best-fit spectral measurements of the parameter of the phabs(apec)

<figure><img src="image_26.png" /><figcaption>**Fig. B.1.** Ratio of the abundances of the different elements (indicated with their atomic number $Z_{elem}$) between the values quoted in *angr* and *lodd09* and the values of reference in *aspl*. In brown (*lodd09*) and filled diamonds (*angr*), we show the relative abundance of helium ($Z_{elem} = 2$) and iron ($Z_{elem} = 26$).</figcaption></figure>

We have estimated the effects of changing the adopted abundance table in XSPEC on the best-fit spectral parameters of a thermal model, either absorbed or not. We refer to the spectral model phabs(apec) set as a reference and in which both components depend on the assumed table for the metal abundance<fnref n="6" />.

We simulate spectra in XSPEC convolved with the *XMM* pn response using the *fake* command with no background, and defining an exposure time ($t_{exp} = 10^6$ sec) and a normalization of the thermal model ($K_{apec} = 1$) high enough that the count statistic is not a limitation. We evaluate the impact of the assumed table for the metal abundance in two cases, one with $N_H = 0.05 \times 10^{22} \text{ cm}^{-2}$ and fixing $N_H = 0$. The latter condition allows us to quantify the impact on the thermal model only. We simulate the spectrum with a metallicity of $Z = 0.3$ times the solar values in *aspl* (Asplund et al. 2009), and investigate how the best-fit parameters change by changing the abundance table with the XSPEC command *abund*. We consider all the abundance tables available in XSPEC<fnref n="7" />, plus Lodders et al. 2009 (hereafter *lodd09*) that is not available in XSPEC and is loaded with the command *abund file ladders09.tab*, where “ladder09.tab” is a one-column file with the number density relative to hydrogen of the elements with atomic number $Z_{elem}$ from

<table><caption>**Table B.1.** Variation in percentage with respect to *aspl*. The quoted values indicate the minimum, mean and maximum over the range investigated in temperature and redshift. The best-fit parameters are the apec normalization *K*, the temperature *T* and the metal abundance *Z*, and are estimated with different $n_H$.</caption><thead><tr><th></th><th colspan="2">$N_H = 0.05 \times 10^{22} \text{ cm}^{-2}$</th></tr><tr><th></th><th>angr</th><th>lodd09</th></tr></thead><tbody><tr><th>K</th><td>-0.5, 1.4, 3.8</td><td>-2.7, -2.3, -1.7</td></tr><tr><th>T</th><td>-8.7, -5.7, -3.0</td><td>-3.8, -2.7, -1.8</td></tr><tr><th>Z</th><td>-38.5, -32.6, -25.1</td><td>-8.5, -3.1, 2.0</td></tr></tbody><tbody><tr><th></th><th colspan="2">$N_H = 0$</th></tr><tr><th></th><th>angr</th><th>lodd09</th></tr><tr><th>K</th><td>-3.9, -3.5, -2.5</td><td>-4.7, -4.2, -3.9</td></tr><tr><th>T</th><td>0.0, 0.3, 1.0</td><td>-0.4, -0.3, -0.2</td></tr><tr><th>Z</th><td>-30.3, -29.3, -26.5</td><td>-2.4, -0.7, 0.2</td></tr></tbody></table>

1 to 30. We run the procedure using a range of temperatures (2, 5, 8 keV) and redshifts (0.05, 0.2, 0.6) defined to enclose the properties of the systems that will be studied in the CHEX-MATE analysis. Each spectrum is rebinned to have at least 1 count in each bin and fitted with Cash statistics over the range 0.5–8 keV.

We focus our analysis on the changes due to the use of the abundance tables *angr* and *lodd09* in fitting simulated spectra produced with the abundance table *aspl*. The values in the tabulated metal abundances are plotted in Fig. B.1 as a function of the atomic number and normalized to the ones in *aspl*. The main differences with respect to *aspl* are the higher (by ~40%) Fe abundance in *angr* and the lower (by 12%) He abundance in *lodd09* (see also Ettori et al. 2020).

The outputs of the spectral analysis are saved and the results are plotted in Fig. B.2 and summarized in Table B.1. We conclude that the differences are more severe (with an average bias of -30%) on the estimates of the metal abundance in *angr* because of its higher reference abundance of Fe, whereas they are in the order of few percent on the measurements of *K* and *T*. When an $N_H = 0.05 \times 10^{22} \text{ cm}^{-2}$ is considered, the mean bias in normalization *K* is -2(+2)% for *lodd09* (*angr*). On the temperature *T*, it is -3% for *lodd09* and -6% for *angr*. The metallicity *Z* is biased by -3% (-30%) in *lodd09* (*angr*). When no Galactic absorption is considered and just the effect on the thermal model is quantified, the biases using *lodd09* with respect to the adopted *aspl* are in the order of < 1% on *T*, -1% on *Z*, and about -4% on *K*; using *angr*, they are similar on *T* and *K*, and are still -30% on *Z*.

We emphasize that the estimates performed here for $N_H \neq 0$ are valid only for the phabs model, which uses the abundance table set for the XSPEC fit. Other models, such as wabs and tbabs use their own hard-coded abundance table.

<footnote n="6">⁶ about the dependence of the phabs model see https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/node249.html</footnote>

<footnote n="7">⁷ see https://heasarc.nasa.gov/xanadu/xspec/manual/ XSabund.html also for references</footnote>
<figure><img src="image_27.png" /><figcaption>**Fig. B.2.** Ratios between the best-fit parameters and the input values. Points are color-coded according to the adopted abundance table (aspl: red, lodd09: yellow, aspl: blue, grsa: green, wilm: purple). Different symbols refer to different redshifts ($z = 0.05$: open diamonds; $z = 0.2$: open circles; $z = 0.6$: filled circles). (Left) Assuming $N_H = 0.05 \times 10^{22} \text{ cm}^{-2}$; (right) with $N_H = 0$.</figcaption></figure>
## Appendix C: Data and figures of individual clusters

In this Appendix, we report the data, the plots of the temperature profiles, and the images with masked point sources and annuli used for spectral extraction for all the clusters in our sample. The full table with the temperature profiles is available in electronic format at CDS. Here we show an excerpt in Table C.1. The first column shows the name of the cluster, while the second and third ones show the sky extension $R_{500}$ in arcminutes and

the mean temperature $T_{[0.15-0.75]R_{500}}$, estimated in the annulus between 0.15 and $0.75R_{500}$. Columns 4 and 5 report the coordinates of the center used for spectral extraction (the peak in Bartalucci et al. 2023), while columns 6 and 7 show the center of the radial bin $R$ and the width $dR$ (annulus included between $R-dR$ and $R+dR$). Columns 8, 9, and 10 report the best fit temperature and its errors (downward and upward) obtained with our baseline MCMC procedure, while column 11 shows the $SOU/BKG$ ratio in the radial bin.

<table><caption>**Table C.1.** Temperature profiles for the clusters of our sample. The full table is available at CDS.</caption><thead><tr><th>Name</th><th>R<sub>500</sub><br>arcmin</th><th>T<sub>mean</sub><br>keV</th><th>R.A.<br>deg</th><th>Dec.<br>deg</th><th>radius<br>arcmin</th><th>width<br>arcmin</th><th>T<br>keV</th><th>dT<sub>down</sub><br>keV</th><th>dT<sub>up</sub><br>keV</th><th>SOU/BKG</th></tr></thead><tbody><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>0.0825</td><td>0.0825</td><td>7.74</td><td>0.91</td><td>0.82</td><td>29.69</td></tr><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>0.2475</td><td>0.0825</td><td>6.95</td><td>0.58</td><td>0.57</td><td>29.72</td></tr><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>0.4400</td><td>0.1100</td><td>6.76</td><td>0.45</td><td>0.47</td><td>18.40</td></tr><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>0.6875</td><td>0.1375</td><td>7.51</td><td>0.47</td><td>0.46</td><td>11.29</td></tr><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>0.9900</td><td>0.1650</td><td>6.55</td><td>0.38</td><td>0.36</td><td>7.28</td></tr><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>1.3750</td><td>0.2200</td><td>6.60</td><td>0.33</td><td>0.30</td><td>5.01</td></tr><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>1.8700</td><td>0.2750</td><td>8.44</td><td>0.54</td><td>0.53</td><td>3.38</td></tr><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>2.5025</td><td>0.3575</td><td>6.66</td><td>0.36</td><td>0.38</td><td>1.98</td></tr><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>3.3275</td><td>0.4675</td><td>6.44</td><td>0.50</td><td>0.51</td><td>1.01</td></tr><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>4.4000</td><td>0.6050</td><td>6.81</td><td>0.86</td><td>0.86</td><td>0.35</td></tr><tr><td>PSZ2 G008.31-64.74</td><td>4.51</td><td>6.83</td><td>344.7010</td><td>-34.8005</td><td>5.6925</td><td>0.6875</td><td>6.40</td><td>1.86</td><td>1.17</td><td>0.12</td></tr><tr><td>PSZ2 G041.45+29.10</td><td>6.48</td><td>6.35</td><td>259.4366</td><td>19.6766</td><td>0.1100</td><td>0.1100</td><td>6.14</td><td>0.73</td><td>0.71</td><td>12.85</td></tr><tr><td>PSZ2 G041.45+29.10</td><td>6.48</td><td>6.35</td><td>259.4366</td><td>19.6766</td><td>0.3025</td><td>0.08250</td><td>7.37</td><td>0.66</td><td>0.67</td><td>13.4</td></tr></tbody></table>

<figure><img src="image_28.png" /><figcaption>**Fig. C.1.** Image and temperature profile of PSZ2G008.31-64.74. The left panel shows the cluster image in the [0.7 – 1.2] keV band, with overlaid masked point sources and the annulus used for estimating the sky background. In the central panel we show the annular regions used for the extraction of the temperature profile, which is represented in the right panel, where the dashed line marks the position of $R_{500}$.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Fig. C.2. Same as Fig. C.1 for cluster PSZ2G041.45+29.10.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. C.3. Same as Fig. C.1 for cluster PSZ2G042.81+56.61.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. C.4. Same as Fig. C.1 for cluster PSZ2G046.88+56.48.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Fig. C.5. Same as Fig. C.1 for cluster PSZ2G050.40+31.17.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Fig. C.6. Same as Fig. C.1 for cluster PSZ2G056.77+36.32.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. C.7. Same as Fig. C.1 for cluster PSZ2G056.93-55.08.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Fig. C.8. Same as Fig. C.1 for cluster PSZ2G057.78+52.32.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Fig. C.9. Same as Fig. C.1 for cluster PSZ2G057.92+27.64.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Fig. C.10. Same as Fig. C.1 for cluster PSZ2G066.41+27.03.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>Fig. C.11. Same as Fig. C.1 for cluster PSZ2G072.62+41.46.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Fig. C.12. Same as Fig. C.1 for cluster PSZ2G077.90-26.63.</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Fig. C.13. Same as Fig. C.1 for cluster PSZ2G083.86+85.09.</figcaption></figure>
<figure><img src="image_41.png" /><figcaption>Fig. C.14. Same as Fig. C.1 for cluster PSZ2G113.29-29.69.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Fig. C.15. Same as Fig. C.1 for cluster PSZ2G113.91-37.01.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Fig. C.16. Same as Fig. C.1 for cluster PSZ2G114.79-33.71.</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>Fig. C.17. Same as Fig. C.1 for cluster PSZ2G149.39-36.84.</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>Fig. C.18. Same as Fig. C.1 for cluster PSZ2G195.75-24.32.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>Fig. C.19. Same as Fig. C.1 for cluster PSZ2G207.88+81.31.</figcaption></figure>
<figure><img src="image_47.png" /><figcaption>Fig. C.20. Same as Fig. C.1 for cluster PSZ2G224.00+69.33.</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>Fig. C.21. Same as Fig. C.1 for cluster PSZ2G238.69+63.26.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>Fig. C.22. Same as Fig. C.1 for cluster PSZ2G243.15-73.84.</figcaption></figure>
<figure><img src="image_50.png" /><figcaption>Fig. C.23. Same as Fig. C.1 for cluster PSZ2G243.64+67.74.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Fig. C.24. Same as Fig. C.1 for cluster PSZ2G277.76-51.74.</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>Fig. C.25. Same as Fig. C.1 for cluster PSZ2G287.46+81.12.</figcaption></figure>
<figure><img src="image_53.png" /><figcaption>Fig. C.26. Same as Fig. C.1 for cluster PSZ2G313.33+61.13.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>Fig. C.27. Same as Fig. C.1 for cluster PSZ2G313.88-17.11.</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>Fig. C.28. Same as Fig. C.1 for cluster PSZ2G324.04+48.79.</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>Fig. C.29. Same as Fig. C.1 for cluster PSZ2G340.94+35.07.</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>Fig. C.30. Same as Fig. C.1 for cluster PSZ2G349.46-59.95.</figcaption></figure>