© The Authors 2024

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⁷ − 10⁸ 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 galactive 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 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₅₀₀. 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₅₀₀. 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₅₀₀, 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 ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7. Errors are given at a 68% confidence level, unless otherwise stated. Whenever we use the notation M₅₀₀ 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₅₀₀ (ranging from 2.9 to 13.5 arcmin) and Galactic absorption (N_H spanning one order of magnitude from 1.03 × 10²⁰ cm⁻² to 1.02 × 10²¹ cm⁻²) 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.

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.

We confirm our sample definition by checking that the DR1 sample is representative of the original CHEX-MATE 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₅₀₀ 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₅₀₀) the temperatures estimated with different analysis are consistent within 1−2%. Differences raise to about 10% at R₅₀₀, 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, 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¹, showing that no major change has been applied in the meantime, with the only exception of the optional cross-calibration with NUSTAR² 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³ 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 σ = 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 different (R₅₀₀ ranging from 2.9′ to 13.5′). For all the clusters with R₅₀₀ ≤ 9′, we define an annulus between R₂₀₀ = R₅₀₀/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₅₀₀ > 12′ for which we requested an offset observation for this purpose (CHEX-MATE Collaboration 2021). Conversely, for the six clusters with 9′< R₅₀₀ < 12′ we define a close background (closebkg) region as an annulus between 1.1 R₅₀₀ 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).

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.

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 χ² 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 Δ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σ 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. 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

$$\begin{array}{c}\hfill {\mathrm{outFOV}}_{\mathrm{PN}}=A_{\mathrm{CRPB}}{\mathrm{outFOV}}_{\mathrm{MOS}2}+A_{\mathrm{SP}}(\mathrm{inFOV}-\mathrm{outFOV}).\end{array}$$(1)

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

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

$$\begin{array}{c}\hfill {\mathrm{CRPB}}_{\mathrm{PN}}=A_{\mathrm{CRPB}}{\mathrm{outFOV}}_{\mathrm{MOS}2}.\end{array}$$(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_QSPE^−Γ with a fixed slope (Γ = 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⁴. 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 emission (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⁵. 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₅₀₀, 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).

4.3.1. 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.

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 R500, 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.

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₅₀₀ < 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 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₅₀₀ < 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 normalization by 1 − 4%, which is a much smaller error than the 17% due to the use of XMM-only data for extended clusters. Moreover 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 background 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 properties of our physical background model, how we built and calibrated 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 latter 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 parameters 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 analysis).

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 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_PN and outFOV_MOS2 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 indicator 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₅₀₀, hampering the goal of reliably measuring the temperature profile up to R₅₀₀ for some clusters. As highlighted by the color-coding, most of them are located in extended clusters (R₅₀₀ > 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.

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 (R500 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 R500. Most of the regions below the threshold SOU/BKG = 0.2 at R < R500 belong to very extended clusters and are thus located in poorly exposed regions of the FOV.

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.

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_ML − T_MCMC)/T_MCMC = ( − 0.7 ± 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_ML − T_MCMC)/T_MCMC = ( − 0.4 ± 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 ± 2.5)% in regions with 0.2 < SOU/BKG < 0.5 and to ( − 7.5 ± 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σ for R < 0.4 R₅₀₀ and 0.4 − 0.8 R₅₀₀ and measure a difference of −5.6 ± 2.7% in the last radial bin (R > 0.8 R₅₀₀). 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 (≤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.

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.

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 uncertainties 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₅₀₀ 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⁵. After the rerun, the error bars increased significantly. We decided not to rerun all our fits with longer MCMC chains and that 5 × 10⁴ 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_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_HI) and computing the relative difference of the best fit values obtained with and without the molecular component: Δ_TNH = (T_NHI − T_NHtot)/T_NHtot. We then divided our sample into three subsamples: a low absorption regime (N_HI < 2 × 10²⁰ cm⁻²), an intermediate regime (2 < N_HI/10²⁰ cm⁻² < 5), and a (relatively) high absorption sample (N_HI > 5 × 10²⁰ cm⁻²). We then computed the median value of Δ_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 CHEX-MATE 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 Δ_TNH = (1.5 ± 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_ab) with respect to the values obtained with the reference table in our analysis (T_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 with N_H = 0 show negligible variations (< 1%) in the temperature estimate. Conversely, when we use an intermediate absorption, N_H = 5 × 10²⁰ cm⁻², 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₂₀₀ 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.

Fig. 7. Ratio of the difference between the temperature values obtained using XMM only in the background region (TNORASS) and that obtained with our standard pipeline (TRASS, joint fit with RASS in the background region) over TRASS, 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 R500 > 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 R500 > 9′, significantly different than the mean with the same criteria.

First, we focus on the six clusters of our sample with R₅₀₀ > 9′, for which we do not have a clean region in the FOV and we used the “closebkg” annulus (1.1 R₅₀₀ − 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% 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 ± 4.5)% (see Fig. 7). We 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₅₀₀ > 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_RASS with T_NORASS also for the remaining 22 clusters with R₅₀₀ < 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_NORASS − T_RASS)/T_RASS = (0.2 ± 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)_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)_M2 < 0.04 cts s⁻¹) and six with a significant residual contamination (highSP with 0.04 < (inFOV − outFOV)_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_NOSP − T_SP)/T_SP = (14 ± 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 ± 5%.

Fig. 8. Ratio of the difference between the temperature values obtained neglecting the soft proton component (TNOSP) and the one of our pipeline (TSP, including a SP component) over TSP 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−1), red are for observations with a low contamination. The magenta triangles shows the median value for clusters with low contamination, when the median is significantly different from the mean.

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_s, highSP − T_s, lowSP)/T_s, lowSP, where the notation T_s marks the scaled temperature T/T_{[0.15 − 0.75] R500}) as a function of the scaled radius R/R₅₀₀ 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.

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.

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 σ, in different SOU/BKG or radial ranges. This means that the true systematic errors in our method are smaller than the 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.

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₅₀₀, 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 σ_T/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.

Fig. 10. Cumulative distribution of clusters where the temperature measurement in the [0.8 − 1.2] R500 bin reaches a nominal error as a function of the relative statistical error σT/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%).

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.

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.

We rescaled the profiles according to self-similar predictions in the right panel of Fig. 11. We divide the radial coordinate by R₅₀₀, 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₅₀₀, 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.

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.

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 ± dr_i which falls in the radial bin R ± 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₅₀₀ and start to show some differences only at larger radii. While up to 0.8 R₅₀₀ differences are not significant, in the radial bin between 0.8 and 1 R₅₀₀, 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σ) 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₅₀₀, which is one of the goals of the CHEX-MATE project.

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.

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₅₀₀.

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 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₅₀₀ 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₅₀₀, 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₅₀₀: most T_i/T_{[0.15 − 0.75] R500} 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%).

Fig. 14. Distributions of the rescaled temperatures as a function of the SOU/BKG in the extraction region for the measurements encompassing 0.9 R500. The vertical dashed lines mark the thresholds used in Fig. 13, with the same color scheme.

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₅₀₀ the difference with the median is maximal, with the mean returning an intermediate value. Conversely, out of the 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.

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.

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₅₀₀ and radii by R₅₀₀. 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 profile of the intrinsic scatter by modeling the individual scaled $T_{\mathrm{i}}^{\mathrm{s}}=T_i/T_{[0.15-0.75]R500}$ 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₅₀₀, as shown in Sect. 6.1, while for the black one we used an independent estimate of T₅₀₀, derived from the SZ M₅₀₀, using Eq. (10) in Ghirardini et al. (2019). We immediately notice that the scatter profile with the T_{[0.15 − 0.75] R500} scaling is significantly lower than the one with the external mass scaling in the radial range [0.1 − 0.6] R₅₀₀. 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.

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] R500 annulus (red points) and the temperature T500 estimated from the mass (black points).

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₅₀₀ 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 CHEX-MATE 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 universal 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 CHEX-MATE 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] R500} = N(R/R₅₀₀)^−α. We decided to fit the profiles for R > 0.3 R₅₀₀ 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 α = 0.37 ± 0.04 and the intrinsic scatter σ_int = 0.11 ± 0.01. If we exclude the regions with SOU/BKG < 0.2, the slope becomes slightly flatter α = 0.29 ± 0.03 and the scatter reduces to σ_int = 0.07 ± 0.01. We also split the radial range in two bins (0.3 < R/R₅₀₀ < 0.55 and 0.55 < R/R₅₀₀ < 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₅₀₀ radial bin we do not have points with SOU/BKG < 0.2, therefore we performed a single fit and measured α = 0.25 ± 0.05. In the [0.55 − 1] R₅₀₀, the profile steepens and we find α = 0.70 ± 0.15 for the full sample and 0.44 ± 0.15 for SOU/BKG > 0.2 (see also Table 5).

Fig. 17. Fit of the temperature profile with a power law for R > 0.3 R500. 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).

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 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₅₀₀, as we do for our data, with the only difference that R₅₀₀ in the comparison samples is obtained from M_YX 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₅₀₀ in [1 − 8]×10¹⁴ M_⊙ (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₅₀₀, with respect to the other samples. The only significant difference is with REXCESS in the [0.8 − 1] R₅₀₀ 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.

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] R500. 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 T500 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.

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₅₀₀. We note that in all cases we used T₅₀₀ 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.

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.

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σ 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₅₀₀, 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 (α = −0.34 ± 0.18 for 0.56 < R/R₅₀₀ < 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] R500}, 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). 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₅₀₀, 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 profileshown in Fig. 20. The mean scaled profile, obtained by rescaling temperatures by the value measured in the annulus [0.15 − 0.75] R₅₀₀, shows a nice agreement with the observed profile at R > 0.3 R₅₀₀. 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₅₀₀ radial range, consistent with the value for the DR1 sample (see Table 5). It is interesting to note that for R > 0.55 R₅₀₀ 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.

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.

To compute the profile of intrinsic scatter, also shown in Fig. 20, both for observations and simulations we rescale the temperature by T₅₀₀ 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 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₅₀₀ (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] R500} 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₅₀₀ 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₅₀₀, 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.

Acknowledgments

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 EPICSSC 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

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)∝E^−Γ, 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 Γ = 0.6 that is nearly independent of inFOV − outFOV.

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.

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_QSPE^−Γ with a fixed slope Γ = 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 take 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,

$$\begin{array}{c}\hfill \frac{N_{\mathit{QSP}}}{\mathit{PROTSCAL}}={\mathrm{A}}_{\mathrm{SP}}{\left(\frac{inFOV-outFOV}{0.05\mathrm{cts}/\mathrm{s}}\right)}^{{\mathrm{B}}_{\mathrm{SP}}}\end{array}$$(A.1)

Fig. A.2. Normalization NQSP 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σ error envelope around the model, whereas the cyan shaded area shows the intrinsic scatter of the relation.

and a free log-normal intrinsic scatter. For MOS1 we measure A_SP, MOS1 = (4.73 ± 0.29)×10⁻⁵ and B_SP, MOS1 = 1.15 ± 0.10 with an intrinsic scatter ${\sigma }_{logN}=0.{16}_{-0.05}^{+0.07}$ dex. We obtain very similar values for MOS2, with A_SP, MOS2 = (4.53 ± 0.31)×10⁻⁵, B_SP, MOS2 = 1.06 ± 0.09 and an intrinsic scatter ${\sigma }_{logN}=0.{18}_{-0.05}^{+0.06}$ 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:

$$\begin{array}{c}\hfill {(inFOV-outFOV)}_{\mathit{PN}}=C{R}_{\mathit{ann}}-{\mathrm{A}}_{\mathrm{CRPB}}outFO{V}_{MOS2}\end{array}$$(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 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 ± 0.08)×10⁻⁴, B_SP, PN = 0.98 ± 0.05, and ${\sigma }_{logN}=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.

Fig. A.3. Normalization NQSP 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σ 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)

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⁶.

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⁶ 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 × 10²² cm⁻² 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⁷, plus Lodders et al. 2009 (hereafter lodd09) that is not available in XSPEC and is loaded with the command abund file lodders09.tab, where “lodder09.tab” is a one-column file with the number density relative to hydrogen of the elements with atomic number Z_elem from 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).

Fig. B.1. Ratio of the abundances of the different elements (indicated with their atomic number Zelem) 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 (Zelem = 2) and iron (Zelem = 26).

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 × 10²²cm⁻² 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.

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 NH = 0.05 × 1022 cm−2; (right) with NH = 0.

We emphasize that the estimates performed here for N_H ≠ 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.

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₅₀₀ in arcminutes and the mean temperature T_{[0.15 − 0.75]R500}, estimated in the annulus between 0.15 and 0.75R₅₀₀. 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.

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 R500.

Fig. C.2. Same as Fig. C.1 for cluster PSZ2G041.45+29.10.

Fig. C.3. Same as Fig. C.1 for cluster PSZ2G042.81+56.61.

Fig. C.4. Same as Fig. C.1 for cluster PSZ2G046.88+56.48.

Fig. C.5. Same as Fig. C.1 for cluster PSZ2G050.40+31.17.

Fig. C.6. Same as Fig. C.1 for cluster PSZ2G056.77+36.32.

Fig. C.7. Same as Fig. C.1 for cluster PSZ2G056.93-55.08.

Fig. C.8. Same as Fig. C.1 for cluster PSZ2G057.78+52.32.

Fig. C.9. Same as Fig. C.1 for cluster PSZ2G057.92+27.64.

Fig. C.10. Same as Fig. C.1 for cluster PSZ2G066.41+27.03.

Fig. C.11. Same as Fig. C.1 for cluster PSZ2G072.62+41.46.

Fig. C.12. Same as Fig. C.1 for cluster PSZ2G077.90-26.63.

Fig. C.13. Same as Fig. C.1 for cluster PSZ2G083.86+85.09.

Fig. C.14. Same as Fig. C.1 for cluster PSZ2G113.29-29.69.

Fig. C.15. Same as Fig. C.1 for cluster PSZ2G113.91-37.01.

Fig. C.16. Same as Fig. C.1 for cluster PSZ2G114.79-33.71.

Fig. C.17. Same as Fig. C.1 for cluster PSZ2G149.39-36.84.

Fig. C.18. Same as Fig. C.1 for cluster PSZ2G195.75-24.32.

Fig. C.19. Same as Fig. C.1 for cluster PSZ2G207.88+81.31.

Fig. C.20. Same as Fig. C.1 for cluster PSZ2G224.00+69.33.

Fig. C.21. Same as Fig. C.1 for cluster PSZ2G238.69+63.26.

Fig. C.22. Same as Fig. C.1 for cluster PSZ2G243.15-73.84.

Fig. C.23. Same as Fig. C.1 for cluster PSZ2G243.64+67.74.

Fig. C.24. Same as Fig. C.1 for cluster PSZ2G277.76-51.74.

Fig. C.25. Same as Fig. C.1 for cluster PSZ2G287.46+81.12.

Fig. C.26. Same as Fig. C.1 for cluster PSZ2G313.33+61.13.

Fig. C.27. Same as Fig. C.1 for cluster PSZ2G313.88-17.11.

Fig. C.28. Same as Fig. C.1 for cluster PSZ2G324.04+48.79.

Fig. C.29. Same as Fig. C.1 for cluster PSZ2G340.94+35.07.

Fig. C.30. Same as Fig. C.1 for cluster PSZ2G349.46-59.95.

All Tables

All Figures

	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.
In the text

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.

In the text

	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)).
In the text

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 R500, 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.

In the text

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 (R500 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 R500. Most of the regions below the threshold SOU/BKG = 0.2 at R < R500 belong to very extended clusters and are thus located in poorly exposed regions of the FOV.

In the text

	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.
In the text

Fig. 7. Ratio of the difference between the temperature values obtained using XMM only in the background region (TNORASS) and that obtained with our standard pipeline (TRASS, joint fit with RASS in the background region) over TRASS, 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 R500 > 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 R500 > 9′, significantly different than the mean with the same criteria.

In the text

Fig. 8. Ratio of the difference between the temperature values obtained neglecting the soft proton component (TNOSP) and the one of our pipeline (TSP, including a SP component) over TSP 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−1), red are for observations with a low contamination. The magenta triangles shows the median value for clusters with low contamination, when the median is significantly different from the mean.

In the text

	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.
In the text

	Fig. 10. Cumulative distribution of clusters where the temperature measurement in the [0.8 − 1.2] R500 bin reaches a nominal error as a function of the relative statistical error σT/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%).
In the text

	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.
In the text

	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.
In the text

	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.
In the text

	Fig. 14. Distributions of the rescaled temperatures as a function of the SOU/BKG in the extraction region for the measurements encompassing 0.9 R500. The vertical dashed lines mark the thresholds used in Fig. 13, with the same color scheme.
In the text

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.

In the text

	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] R500 annulus (red points) and the temperature T500 estimated from the mass (black points).
In the text

Fig. 17. Fit of the temperature profile with a power law for R > 0.3 R500. 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).

In the text

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] R500. 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 T500 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.

In the text

	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.
In the text

	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.
In the text

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.

In the text

Fig. A.2. Normalization NQSP 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σ error envelope around the model, whereas the cyan shaded area shows the intrinsic scatter of the relation.

In the text

	Fig. A.3. Normalization NQSP 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σ envelope, while the light blue shows the intrinsic scatter.
In the text

	Fig. B.1. Ratio of the abundances of the different elements (indicated with their atomic number Zelem) 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 (Zelem = 2) and iron (Zelem = 26).
In the text

	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 NH = 0.05 × 1022 cm−2; (right) with NH = 0.
In the text

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 R500.

In the text

	Fig. C.2. Same as Fig. C.1 for cluster PSZ2G041.45+29.10.
In the text

	Fig. C.3. Same as Fig. C.1 for cluster PSZ2G042.81+56.61.
In the text

	Fig. C.4. Same as Fig. C.1 for cluster PSZ2G046.88+56.48.
In the text

	Fig. C.5. Same as Fig. C.1 for cluster PSZ2G050.40+31.17.
In the text

	Fig. C.6. Same as Fig. C.1 for cluster PSZ2G056.77+36.32.
In the text

	Fig. C.7. Same as Fig. C.1 for cluster PSZ2G056.93-55.08.
In the text

	Fig. C.8. Same as Fig. C.1 for cluster PSZ2G057.78+52.32.
In the text

	Fig. C.9. Same as Fig. C.1 for cluster PSZ2G057.92+27.64.
In the text

	Fig. C.10. Same as Fig. C.1 for cluster PSZ2G066.41+27.03.
In the text

	Fig. C.11. Same as Fig. C.1 for cluster PSZ2G072.62+41.46.
In the text

	Fig. C.12. Same as Fig. C.1 for cluster PSZ2G077.90-26.63.
In the text

	Fig. C.13. Same as Fig. C.1 for cluster PSZ2G083.86+85.09.
In the text

	Fig. C.14. Same as Fig. C.1 for cluster PSZ2G113.29-29.69.
In the text

	Fig. C.15. Same as Fig. C.1 for cluster PSZ2G113.91-37.01.
In the text

	Fig. C.16. Same as Fig. C.1 for cluster PSZ2G114.79-33.71.
In the text

	Fig. C.17. Same as Fig. C.1 for cluster PSZ2G149.39-36.84.
In the text

	Fig. C.18. Same as Fig. C.1 for cluster PSZ2G195.75-24.32.
In the text

	Fig. C.19. Same as Fig. C.1 for cluster PSZ2G207.88+81.31.
In the text

	Fig. C.20. Same as Fig. C.1 for cluster PSZ2G224.00+69.33.
In the text

	Fig. C.21. Same as Fig. C.1 for cluster PSZ2G238.69+63.26.
In the text

	Fig. C.22. Same as Fig. C.1 for cluster PSZ2G243.15-73.84.
In the text

	Fig. C.23. Same as Fig. C.1 for cluster PSZ2G243.64+67.74.
In the text

	Fig. C.24. Same as Fig. C.1 for cluster PSZ2G277.76-51.74.
In the text

	Fig. C.25. Same as Fig. C.1 for cluster PSZ2G287.46+81.12.
In the text

	Fig. C.26. Same as Fig. C.1 for cluster PSZ2G313.33+61.13.
In the text

	Fig. C.27. Same as Fig. C.1 for cluster PSZ2G313.88-17.11.
In the text

	Fig. C.28. Same as Fig. C.1 for cluster PSZ2G324.04+48.79.
In the text

	Fig. C.29. Same as Fig. C.1 for cluster PSZ2G340.94+35.07.
In the text

	Fig. C.30. Same as Fig. C.1 for cluster PSZ2G349.46-59.95.
In the text