© The Authors 2024

1. Introduction

Galaxy clusters are dark-matter-dominated astrophysical objects with total masses¹ in the range M₅₀₀ ≃ 10¹⁴ − 10¹⁵ M_⊙. Following the hierarchical scenario of structure formation, they reside at the nodes of the cosmic web and grow through accretion of matter along the filaments of the large-scale structure and through episodic mergers of small mass systems. Falling into the potential well of the dark matter, the intracluster medium (ICM) is heated to X-ray-emitting temperatures (∼10⁷ − 10⁸ K) by shocks and compression. Given the scale-free nature of gravity, the resulting X-ray cluster population is expected to be self-similar and for tight scaling relations to exist between the ICM observables and cluster mass and redshift (Kaiser 1986; Bryan & Norman 1998). However, second-order effects, mainly linked to feedback from active galactic nuclei (AGNs) and radiative cooling of the gas, act to modify the properties of the ICM and induce some degrees of departure from self-similar predictions (see e.g. Voit 2005; Pratt et al. 2009, 2010; Borgani & Kravtsov 2011; Giodini et al. 2013; Gaspari et al. 2020; Lovisari & Maughan 2022).

In recent years, spatially resolved observations have allowed for a more detailed examination of the impact of non-gravitational processes on the ICM. In particular, the study of the thermodynamic properties of the ICM has been of great interest, since it allows us to obtain useful insights into the history and level of the energy deposited in the ICM through feedback processes. Gas entropy, defined as K = T/n_e^2/3 (Ponman et al. 1999; Lloyd-Davies et al. 2000), where T is the gas temperature and n_e the electron density², plays a key role in this context, since it both determines the structure of the ICM and provides a record of the processes (both gravitational and non gravitational) that influence the properties of the ICM. Entropy is generated during the hierarchical assembly process and then is modified by any other process that can change the physical characteristics of the gas (e.g. AGN heating, cooling, and star formation). For these reasons, it is the ideal tool for investigating the thermodynamic history of the cluster population (see Voit 2005, for a review).

In a stable ICM in hydrostatic equilibrium, stratification naturally results in a radially increasing entropy profile. Tozzi & Norman (2001) used analytical modelling to study the characteristic radial distribution of the entropy, where, outside the core region, it was shown to steadily increase with radius, following a power law with a slope of K ∝ R^1.1. This conclusion was reinforced by the non-radiative (gravity-only) simulations of Voit et al. (2005), where, again outside the core, the entropy slope was found to be remarkably stable out to ∼2 R₂₀₀, independent of system mass. These predictions are regarded as robust and depend little on the adopted numerical approach outside of the central region (Mitchell et al. 2009; Gaspari et al. 2012).

Early observations, however, showed that the entropy profiles of lower-mass systems exceeded that expected from gravity alone (Lloyd-Davies et al. 2000). Subsequent observations have refined this picture. The mass dependence of the entropy excess has been confirmed, and has been shown to extend to a larger radius in lower-mass systems (Pratt et al. 2010). As a result, lower-mass systems generally exhibit a shallower entropy slope, with significantly increased core entropy relative to the expectation from simple gravitational collapse. Moving outwards, the observed profiles converge towards the non-radiative prediction, but the convergence radius is larger for less massive systems. This contributes to the observed large scatter in the central regions (e.g. Cavagnolo et al. 2009; Pratt et al. 2010; Ghirardini et al. 2019). The now-standard explanation for this behaviour is that it is due to the combined effect of non-gravitational processes. These include, but are not limited to, the following: radiative cooling of the central gas, which removes material from the hot phase, leaving only gas on a higher adiabat; stochastic feedback from the central AGN over time, which injects energy into the ICM and raises its entropy; feedback from supernovae (SNe) in the cluster galaxies, which also inject energy into the ICM; and conduction, which smooths out temperature (and therefore entropy) gradients. While numerical simulations that include these processes can go some way to explaining the observed entropy distributions, they still struggle to reproduce the observed population behaviour (see e.g. Nagai et al. 2007; Barnes et al. 2017; Altamura et al. 2023; Oppenheimer et al. 2021; Kay & Pratt 2022).

Interpretation of the observed properties of the ICM entropy is complicated by the clear dependence of the entropy distribution on the morphological properties of the cluster gas (Pratt et al. 2010). Relaxed systems generally possess cool cores, and such objects exhibit an entropy distribution that closely follows the K ∝ R^1.1 distribution down to very small radii. In contrast, disturbed systems always exhibit high core entropy (e.g. Pratt et al. 2010; Babazaki et al. 2018). However, the central surface brightness of morphologically disturbed systems is always flatter than that of cool-core systems of a similar mass. This means that at a given mass, to obtain a spectrum with a given number of counts or signal-to-noise ratio, a larger central region must be used for morphologically disturbed systems than for cool core systems. The resulting entropy profiles therefore have poorer angular resolution in the centre, further complicating the interpretation.

In this context, it is interesting to revisit the entropy distribution of the highest-mass systems. In such objects, gravity is the dominant entropy generation mechanism, and non-gravitational effects should be minimised. As such, they should provide a robust baseline for theoretical works. In addition, such systems are highly luminous in X-rays, allowing measurement of the entropy both deep into the core and out to a large radius. For this work, we used a subset of 32 very high-mass (M₅₀₀ > 7.75 × 10¹⁴ M_⊙) clusters, which we have dubbed HIGHMz, selected from the ‘Cluster HEritage project with XMM-Newton: Mass Assembly and Thermodynamics at the Endpoint of structure formation’ (CHEX-MATE; CHEX-MATE Collaboration 2021), a multiyear Heritage programme to observe 118 Planck SZ-selected clusters with XMM-Newton. CHEX-MATE observations are tailored for deriving the total masses through the hydrostatic equilibrium equation, and for characterising the thermodynamic properties of the entire sample with good precision within R₅₀₀. In this context, the very high mass of the systems in HIGHMz allowed us to further investigate the ICM entropy and its radial distribution in the gravity-dominated regime, and to compare our results with recent large-volume numerical simulations including non-gravitational sub-grid physics. The results from this work will thus be a reference point for future studies at lower-mass scales and higher redshifts.

The paper is organised as follows. In Sect. 2 we present the HIGHMz sample, together with additional samples (both observational and simulated) used for comparison; in Sect. 3 we describe our methods for the data processing and analysis; and Sect. 4 is dedicated to the results of this work where we present the measured entropy profiles and compare them to other samples. Finally, we discuss our findings and provide a conclusion of the work in Sects. 5 and 6. Throughout the paper, we assume a flat Λ cold dark matter (CDM) cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ (i.e. h₇₀ = 1), Ω_m = 0.3, and Ω_Λ = 0.7. The solar abundance table is set to Asplund et al. (2009), while all the quoted errors hereafter are at the 1σ confidence level. Whenever we use the notation M₅₀₀^Y_SZ in this paper, we refer to the masses derived from the Planck SZ (Sunyaev & Zeldovich 1972) 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. Samples

2.1. HIGHMz: High-mass clusters

For the present work, we have selected the most massive systems of the CHEX-MATE (CHEX-MATE Collaboration 2021) project. This sample, which we have dubbed HIGHMz, comprises 35 clusters at redshift z > 0.2 and with M₅₀₀^Y_SZ > 7.75 × 10¹⁴ M_⊙.

We excluded from the present analysis three galaxy clusters, namely PSZ2G107.10+65.32 (a.k.a. A1758N), PSZ2G225.93−19.99 (a.k.a. MACS J0600-20), and PSZ2G339.63−69.34 (a.k.a. Phoenix). The former two are double clusters and do not allow a simple radial analysis. The latter has a strong central AGN, whose emission, given the size of XMM-Newton’s PSF, has an impact also on regions beyond the core. Since including these three complex clusters may introduce systematics in our results, we preferred to keep them aside and to include them subsequently in the final study on the full CHEX-MATE sample, adopting an ad-hoc analysis. The final list of the 32 analysed HIGHMz clusters is presented in Table 1, together with information about cluster masses, redshifts and indicators for their morphological state (light concentration, c; centroid shift, w), as measured by Campitiello et al. (2022). High values of c and low values of w are commonly measured for clusters in a dynamically relaxed state, while low c and high w are commonly associated with disturbed systems (e.g. Santos et al. 2008; Hudson et al. 2010; Lovisari et al. 2017; Campitiello et al. 2022). The cluster masses and redshifts are also shown in Fig. 1, as black dots.

Fig. 1. Distribution of HIGHMz (black dots) clusters in the mass-redshift plane. X-COP (red squares), REXCESS (orange triangles) and ESZ (green stars) clusters are also plotted for comparison. Vertical and horizontal coloured lines mark the median redshift and mass, respectively, of the four samples.

Given its selection in mass and redshift, HIGHMz allows us to address multiple and interesting aims. In particular, it offers the opportunity to study the cluster entropy in a regime where gravitational processes are considered dominant, thus setting a valuable baseline for future studies at lower-mass scales, and to compare the results with local samples available in the literature, thus investigating potential hints of variation of the gas entropy profiles with mass and/or redshift.

2.2. Comparison samples: REXCESS, ESZ, and X-COP

We will compare our results with three cluster samples available in the literature, namely REXCESS (Böhringer et al. 2007), ESZ (Planck Collaboration XI 2011), and X-COP (Eckert et al. 2017). All of these samples have been observed with XMM-Newton and thus allow direct comparison with CHEX-MATE. The four samples together allow investigation of the gas entropy in clusters spanning a wide range both in mass and redshift, thus giving us the possibility to test predictions from self-similar scenario. REXCESS, ESZ, and X-COP are extensively described in the referenced papers; however, we provide here below a summary of their main features:

• REXCESS consists of 31 nearby (z < 0.2) galaxy clusters with temperatures in the range 2 − 9 keV, drawn from the REFLEX catalogue (Böhringer et al. 2004). Therefore, it provides a census of the local X-ray-selected cluster population. REXCESS clusters were selected in X-ray luminosity only, with no bias towards any particular morphological type;

• the original ESZ sample comprises 189 clusters, which were identified via their SZ effect in the first all-sky coverage by the Planck satellite (Planck Collaboration I 2011; Planck Collaboration VIII 2011). After cross-correlating with the MCXC catalogue (Piffaretti et al. 2011) and checking the quality of XMM-Newton archive data, the suitable sample for the X-ray analysis was reduced to 62 clusters. This final ESZ sample consists of clusters at z < 0.5, spanning one decade in mass;

• X-COP is a set of 12 massive and local (z < 0.1) clusters selected from the Planck all-sky survey of SZ sources (Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016). The aim of the X-COP project was to advance the knowledge of the physical conditions in the cluster outskirts, by combining high signal-to-noise SZ data with good X-ray data coverage at and beyond ∼R₅₀₀. The thermodynamic profiles of X-COP clusters from joint XMM-Newton and Planck data were published in Ghirardini et al. (2019). The entropy profiles used here slightly differ from the ones presented in Ghirardini et al. (2019) as they are the result of a joint non-parametric reconstruction of the deprojected thermodynamic profiles, whereby the 3D temperature profile is described as a linear combination of a large number of log-normal functions sampling the radial range of interest. For more details on the deprojection technique, we refer the reader to Eckert et al. (2022).

Masses and redshifts of REXCESS, ESZ, and X-COP samples are shown in Fig. 1 and compared to HIGHMz. In particular, we plot (i) M₅₀₀^Y_SZ masses for HIGHMz; (ii) masses measured using the M₅₀₀ − Y_X scaling (Kravtsov et al. 2006), as calibrated in Arnaud et al. (2010), for REXCESS and ESZ; and (iii) hydrostatic-equilibrium masses for X-COP, as computed in Eckert et al. (2022). The HIGHMz sample has both a higher average mass and a higher average redshift than the other samples. We note that ESZ has 10, 8, and 7 clusters in common with HIGHMz, X-COP, and REXCESS, respectively, while the other three samples do not have clusters in common with each others. We checked that the entropy profiles of these common clusters, as measured for the different samples, are in agreement within the errors.

2.3. Simulated samples: MACSIS and The300

We also compare our results with recent hydrodynamic simulations that are able to reproduce a sufficient number of massive clusters in the mass range of interest. In particular, we consider two simulated datasets, taken from the parent MACSIS (Barnes et al. 2017) and THE THREE HUNDRED (hereafter The300; Cui et al. 2018) projects. Briefly:

• MACSIS is a sample of 390 massive clusters selected from a large-volume (3.2 Gpc box size) dark-matter simulation, which were then re-simulated with full gas physics. MACSIS extends the BAHAMAS (400 h⁻¹ Mpc box size) simulation (McCarthy et al. 2017) to the most massive clusters expected to form in a ΛCDM cosmology. These were both run with a heavily modified version of the GADGET3 TreePM-SPH code, last described in Springel (2005). The code incorporates radiative cooling on an element-by-element basis; heating and cooling from UV/X-ray and cosmic microwave backgrounds; star formation; stellar evolution and enrichment from asymptotic-giant-branch (AGB) stars and Type Ia/II SNe; kinetic stellar feedback; super-massive black hole growth and thermal AGN feedback. MACSIS was run with the same sub-grid parameter choices as BAHAMAS, calibrated to match the low redshift galaxy stellar mass function and cluster gas fractions (see McCarthy et al. 2017 for further details of the baryonic physics model and calibration process);

• The300 (GADGET-X version) consists of the 324 most massive clusters identified at z = 0 within the dark-matter-only MultiDark simulation (Klypin et al. 2016) of size 1 h⁻¹ Gpc and resimulated accounting for baryons. The hydrodynamic code, implementing the SPH description as in Beck et al. (2016), is similar to what was used in Rasia et al. (2015), and includes metal-dependent radiative gas cooling, star formation, stellar feedback by AGB stars, Type-Ia and core-collapse SNe, super-massive black hole growth, and AGN feedback.

Given their properties, both datasets are therefore ideal samples to test the comparison with the massive clusters in HIGHMz.

Following Planck Collaboration XX (2014), we assumed that the mass measurements M₅₀₀^Y_SZ show an average hydrostatic bias of 20% (i.e. 1 − b = 0.8). To ensure a fair comparison with simulations, which measure the true masses of the reproduced clusters, we rescaled HIGHMz masses accordingly, whose ‘true’ values were considered to be M₅₀₀^{Y_SZ, corr.} = M₅₀₀^Y_SZ/0.8. In order to define the most suitable simulated dataset for comparison with HIGHMz, we thus imposed MACSIS and The300 cluster masses to be greater than M₅₀₀^{Y_SZ, corr.}, as also previously done in Bartalucci et al. (2023). The resulting sample contains 25 clusters from The300 and 75 from MACSIS, taken from one snapshot corresponding to redshift z = 0.333 (The300) and z = 0.34 (MACSIS), the closest one for each sample to the median redshift of HIGHMz (z = 0.39). Similarly to observations, we also used information on morphological indicators (in particular, c and w) for The300 and MACSIS, computed as detailed in Campitiello et al. (2022) and Towler et al. (2023), respectively. Figure 2 shows the comparison of both cluster masses and morphological indicators for the three samples, while the complete list of the selected simulated datasets is reported in Tables A.1 and A.2. Both simulation datasets adopt a flat ΛCDM cosmology with h = 0.678 km s⁻¹ Mpc⁻¹, Ω_m = 0.307 and Ω_Λ = 0.693 (Planck Collaboration XIII 2016). We corrected for the difference in the adopted cosmology when comparing to observations (Eq. (1)).

Fig. 2. Properties of MACSIS (light blue) and The300 (violet) clusters, in comparison to HIGHMz (black). Top panel: distribution of the cluster masses. Bottom panel: distribution in the concentration-centroid shift plane. In both panels, coloured lines show the median value of the corresponding axis for the three samples. MACSIS and The300 clusters are all at redshifts z = 0.34 and z = 0.333, respectively.

3. Analysis procedures

The X-ray data have been reduced and analysed using the CHEX-MATE pipeline, which is extensively described in Bartalucci et al. (2023) and Rossetti et al. (2024). In the following, we report the main steps of our procedure and present the emission measure and gas temperature profiles of HIGHMz clusters. Finally, we also summarise our deprojection techniques to build the three-dimensional entropy profiles and provide information on the adopted self-similar scaling.

3.1. Data reduction

The HIGHMz clusters were observed with the European Photon Imaging Camera (EPIC, Turner et al. 2001; Strüder et al. 2001) on board XMM-Newton. The datasets were reprocessed using the Extended-Science Analysis System (ESAS, Snowden et al. 2008) embedded in SAS version 16.1 (see Rossetti et al. 2024 for details of this choice). We removed flare events by using the tools mos-filter and pn-filter, by extracting the light curves in the [2.5 − 8.5] keV energy range and excising time intervals with count rates exceeding 3σ times the mean. Point sources were filtered from the analysis following the scheme detailed in Sect. 2.2.3 of Ghirardini et al. (2019) and summarised in Sect. 3.1.1 of Bartalucci et al. (2023).

For each cluster, we used the XMM-Newton observations listed in Table 1. When multiple pointings were available, we combined them to increase the available count statistics, both for the image and the spectral analysis, as described in Bartalucci et al. (2023) and Sect. 3.3 of Rossetti et al. (2024).

3.2. Projected profiles

3.2.1. Emission measure

The procedure we followed to derive HIGHMz emission measure profiles is outlined in detail in Sect. 3 of Bartalucci et al. (2023). Briefly, we first produced the EPIC image for each cluster in the energy band 0.7 − 1.2 keV, together with exposure and background maps, merging together those of the three XMM-Newton’s cameras. From the combination of these images, we then extracted both azimuthal mean and median background-subtracted and exposure-corrected surface brightness (SB) profiles, from the coordinates of the X-ray peaks. The use of the azimuthal mean allows inspection of the central regions of a cluster with higher resolution³ (e.g. Pratt et al. 2022); conversely, the technique of computing azimuthal median SB profiles was introduced by Eckert et al. (2015) to limit the impact of sub-clumps and sub-structures too faint to be identified and masked (e.g. Roncarelli et al. 2013; Zhuravleva et al. 2013). As specified later, we will make use of both these products in the computation of the ICM entropy. Finally, the SB profiles were converted to emission measures following Eq. (1) of Arnaud et al. (2002). Emission measures are directly related to the gas density of the cluster, being the integral of the density squared along the line of sight. They are therefore necessary ingredients to recover the 3D density profiles.

The azimuthal median emission measure profiles of the HIGHMz clusters are presented in Fig. 3 (left) and colour coded according to their masses M₅₀₀^Y_SZ, with blue to red indicating from less to more massive clusters. As highlighted by Bartalucci et al. (2023) for the entire CHEX-MATE sample, in the central regions of HIGHMz clusters we find considerable dispersion, with cool-core clusters exhibiting steeper profiles, while moving towards the outskirts the profiles appear more self-similar.

Fig. 3. Azimuthal median emission measure (left) and projected temperature (right) profiles of clusters in HIGHMz. The profiles are colour coded according to cluster masses (M500YSZ; Table 1). Temperature measurements with S/B < 0.2 are marked in black.

3.2.2. Gas temperature

The steps of the CHEX-MATE pipeline built for the spectral extraction and fitting, which allowed us to measure the gas temperature profiles of HIGHMz clusters, are detailed in Sects. 3 and 4 of Rossetti et al. (2024). The pipeline combines the best practices developed during previous projects (e.g. Pratt et al. 2007; Planck Collaboration XI 2011; Bartalucci et al. 2018; Ghirardini et al. 2019), with the introduction of some novelties. These include, for example, the construction of a physical model for the particle background for all EPIC detectors and the application of a Bayesian MCMC framework, allowing the propagation of the uncertainties on the background parameters up to the final results of the spectral analysis.

An important point of strength of the CHEX-MATE pipeline is the careful analysis of the underlying systematics. Based on the study of a sample of 30 clusters, selected to be representative of the entire CHEX-MATE population, Rossetti et al. (2024) showed that we can obtain reliable temperature estimates at least up to regions where the source intensity is larger than 20% of the background. Below this value (i.e. where S/B < 0.2) temperature measurements can be biased low and the exclusion of these bins resulted in a better reconstruction of the mean temperature profile in the outer regions (Rossetti et al. 2024).

In the spectral fitting, we fixed the redshift and Galactic absorption, but the temperature, metallicity, and normalisation of the cluster spectra were left free to vary, as described in Rossetti et al. (2024). The measured temperature profiles of HIGHMz clusters are presented in Fig. 3 (right). Temperatures range from ∼5 to ∼15 keV, reflecting the high masses of the clusters in the sample. The adopted colour coding is the same as for the emission measures and highlights the dependence on cluster mass. Temperature measurements where S/B < 0.2 are also shown in the figure and marked in black. Almost all of these are located beyond R₅₀₀, while CHEX-MATE observations are tailored to measure temperatures within this radius. In about four cases, these measurements exhibit very low temperature values. Including these external bins in the deprojection can introduce a bias on the final products, as discussed later.

3.3. Measuring the ICM entropy

Following the astrophysical convention, we construct the ICM entropy as K = T/n_e^2/3, using the three-dimensional profiles of gas density and temperature. The deprojection of both the emission measure and temperature profiles presented in Sect. 3.2 is therefore needed in advance.

The emission measures were deprojected and PSF-corrected using the non-parametric method described in Croston et al. (2006). As also detailed in Croston et al. (2008), the emission measure was converted to gas density by calculating a global conversion factor for each profile in Xspec v. 12.13 (Arnaud 1996), using the spectroscopic temperature in the [0.15 − 1] R₅₀₀^Y_SZ aperture. A secondary correction factor, taking into account radial variations of temperature and abundance, was also included to give the final gas density profiles. This was obtained from analytical fits to the projected quantities (see Pratt & Arnaud 2003). Both azimuthal mean and median density profiles were produced and used to calculate the ICM entropy. We will consider the entropy profiles derived from azimuthal median densities as our reference throughout the paper. However, azimuthal mean densities will be used in Sect. 4.3, for comparison with other samples.

To deproject and PSF-correct the 2D temperature profiles of HIGHMz clusters, we adopted the non-parametric-like technique described in Démoclès et al. (2010) and Bartalucci et al. (2018). Briefly, we assumed that the 3D temperature profiles can be described by a parametric model, adapted from Vikhlinin et al. (2006), that is convolved with a response matrix which simultaneously takes into account projection and PSF redistribution. The projection procedure additionally took into account the bias introduced by fitting isothermal models to multi-temperature plasma (Mazzotta et al. 2004; Vikhlinin 2006). The final 3D temperature profile is then estimated at the weighted radii corresponding to the 2D annular binning scheme. Low and external temperature measurements with S/B < 0.2, often accompanied with small statistical errors, have an impact on the deprojection. This is particularly true for the cluster outskirts, but also for the internal regions. We will show in Sect. 4.1, and discuss in more detail in Sect. 5.1, the impact of including these low-S/B measurements on the shape of the entropy profiles. Unless otherwise stated, in the following sections we will consider the entropy profiles constructed from temperature measurements with S/B > 0.2 as our reference.

3.4. Scaling and self-similar predictions

Among the products of our deprojection pipeline are also integrated masses within R₅₀₀, which we adopted for rescaling the entropy profiles. In particular, temperature, density, and their gradients were used to compute both hydrostatic masses (hereafter $M_{500}^{\mathrm{HE}}$) and the quantity Y_X = M_gas × T (Kravtsov et al. 2006), which was used as a proxy for the total mass (hereafter M₅₀₀^Y_X), using the scaling calibrated in Arnaud et al. (2010). Hydrostatic masses were measured using azimuthal median densities and adopted in the rescaling of the ICM entropy in comparison with X-COP clusters only (Sect. 4.3); masses from Y_X were computed both from median and mean densities and used in the following sections depending on the relevant study.

According to the self-similar scenario, the physical properties of galaxy clusters should coincide, once they are correctly scaled for cluster mass and redshift. For the ICM entropy, the predicted dependence on mass and redshift from pure gravitational collapse is K ∝ M^2/3E(z)^−2/3 (Voit 2005). In this work we adopted the characteristic entropy K₅₀₀, as computed by Pratt et al. (2010):

$$\begin{array}{c}\hfill {\mathrm{K}}_{500}=106(\frac{M_{500}}{h_{70}^{-1}{10}^{14}{\mathrm{M}}_{\odot }}{)}^{2/3}E{(z)}^{-2/3}f{{}_b}^{-2/3}{h}_{70}^{-4/3}{\mathrm{keV}\mathrm{cm}}^2,\end{array}$$(1)

where f_b is the Universal baryon fraction, which we take to be 0.16 (Planck Collaboration XIII 2016). We will show, through the combined sample of HIGHMz, ESZ, X-COP, and REXCESS clusters, that the adopted self-similar scaling may not be adequate to describe the properties of galaxy clusters, due to the impact of non-gravitational processes on the ICM. In Sect. 5.3 and Appendix B, we will provide updated dependencies on mass and redshift that go beyond self-similar predictions.

Predictions from cosmological simulations that implement gravitational processes only (e.g. Voit et al. 2005) can be considered as a baseline to assess the impact of non-gravitational processes on the ICM. These simulations are commonly referred to as ‘non-radiative’ and predict that, in the gravity-dominated regime, cluster entropy outside the core regions should follow a power-law increase out to ∼2 R₂₀₀. Any other physical processes capable of injecting energy into the ICM would result in an overestimate of the ICM entropy with respect to the expected power-law behaviour. Voit et al. (2005) presented a simple formula that provides both the normalisation and slope of this expected power law, normalised for an overdensity Δ = 200. In this work, we adopted the equivalent equation computed by Pratt et al. (2010), written for an overdensity Δ = 500:

$$\begin{array}{c}\hfill K/K_{500}=1.42{(R/R_{500})}^{1.1}.\end{array}$$(2)

Given their notably high masses, one would expect the properties of HIGHMz clusters to be mainly driven by gravitational processes at large radii and so their external entropy profiles to closely follow the prediction by Voit et al. (2005). We will verify if this is the case in the following.

4. Results

4.1. Entropy profiles of HIGHMz clusters

We present in Fig. 4 the entropy profiles of the HIGHMz clusters, derived from azimuthal median densities. Those obtained from azimuthal mean densities have similar properties and are presented in Appendix C. To ensure a better examination of the core regions, for those clusters with poor central temperature resolution, we computed the entropy on the density radial grid, assuming a constant central temperature, as already done in previous works (Donahue et al. 2005; Cavagnolo et al. 2009; Pratt et al. 2010). These additional measurements, if any, are shown using dotted lines in Fig. 4.

Fig. 4. Entropy profiles of HIGHMz clusters, derived using azimuthal median densities. The profiles are shown in physical units in the upper panels, while they are scaled using K500YX in the centre and bottom. In the upper and central panels, they are colour coded according to the mass, while in the bottom we show the dependence on the concentration (left) and centroid shift (right), as measured by Campitiello et al. (2022). The two labels ‘Without S/B cut’ and ‘S/B > 0.2’ indicate whether all the temperature bins were used in the deprojection or we limited to measurements with S/B > 0.2, respectively. In the central right panel, we also plot the median profile (solid black) and the measured intrinsic dispersion (grey shaded area). At small radii, dotted lines denote profiles derived assuming a constant core temperature on the density radial grid. Black dash-dotted lines in the central and lower panels are predictions from non-radiative simulations (Voit et al. 2005).

In the upper panels of Fig. 4, the profiles are shown in physical units, with their radii normalised by R₅₀₀^Y_X. On the left, all the temperature bins have been included in the deprojection, while on the right we have limited the deprojection to temperature bins with S/B > 0.2, to reduce the impact of the uncertainty in temperature measurements at large radii discussed in Sect. 3.2.2. As already noted in previous works (e.g. Cavagnolo et al. 2009; Sun et al. 2009; Pratt et al. 2010), most profiles do not show simple power-law behaviour down to arbitrarily small radii, but some flattening is typically observed in the central regions (R ≲ 0.2 R₅₀₀). Here, considerable dispersion is observed, which then reduces moving outwards, where the measured entropy profiles increase with radius following power laws with almost the same slope. At radii R ≳ 0.6 R₅₀₀, the effect of excluding temperature measurements with S/B < 0.2 is visible observing the shape of the profiles: some of the clusters exhibit significant flattening in the outskirts when no cut for S/B is applied (top left); conversely, including only temperature bins with S/B > 0.2 in the deprojection allows the shape of the entropy profiles to be regularised at all radii, especially in the outskirts, where the entropy profiles resume a steadily increasing trend (top right). As the inclusion of low-S/B temperature measurements may introduce a bias into the ICM entropy reconstruction at large radii, as mentioned above, we therefore consider profiles with a cut for S/B > 0.2 as our reference.

In Fig. 4 (upper panels), gas entropy profiles are colour coded according to their values of M₅₀₀^Y_SZ, with blue to red for least to most massive clusters. This shows that a clear mass dependence is present beyond R ∼ 0.2 R₅₀₀^Y_X, with the most massive clusters exhibiting higher entropy, when this is not rescaled for any global quantity. By adopting the self-similar rescaling for K₅₀₀^Y_X (Eq. (1) using M₅₀₀^Y_X), and thus taking into account the dependence on both cluster mass and redshift, the correlation with cluster mass disappears at large radii (central panels). In the outskirts, we observe perhaps even more clearly the impact of excluding temperature measurements with S/B < 0.2 from the analysis on the shape of the profiles. After this correction, the entropy profiles of massive clusters resume a steadily increasing trend with radius also at ∼R₅₀₀ and more closely follow the predictions from pure gravitational collapse (Eq. (2), Voit et al. 2005), represented in Fig. 4 as the black dash-dotted lines. Conversely, in the core regions, notable scatter is still observed.

In the lower panels of Fig. 4, the entropy profiles of HIGHMz clusters are colour coded according to their morphological indicators, such as the concentration (c, left) and centroid shift (w, right), as measured by Campitiello et al. (2022) and reported in Table 1. In the central regions we find a strong (anti-)correlation with (c) w, meaning that disturbed clusters tend to have higher central entropy, likely due to residual merger activity, while cool-core clusters exhibit lower central entropy, which is thought to arise from a balance between feedback and cooling processes (Pratt et al. 2010). Not surprisingly, the dependence of the ICM entropy on the two morphological indicators is similar, although they are sensitive to different cluster scales; indeed, the correlation with c, which is directly linked to the core properties, is stronger. Finally, we notice that moving outwards the dependence on morphological parameters is reduced and appears even reversed at ∼R₅₀₀, with cool-core clusters featuring higher, and steeper, entropy profiles, reflecting the behaviour of the gas density at these radii (e.g. Maughan et al. 2012). A more quantitative estimate of the dependence of the gas entropy on both c and w will be provided in Sect. 4.4.2, in comparison with simulated datasets.

We also computed the median and the intrinsic scatter of the profiles, following the procedure detailed in Sect. 6.1 of Bartalucci et al. (2023). As a first step, we interpolated each scaled profile linearly in the log-log plane, on a common grid of twelve bins in the [0.02 − 1] R₅₀₀ radial range. The medians were then derived through a Monte Carlo technique: for each radial bin, we generated 1000 random realisations of each entropy measurement, normally distributed around its nominal value and adopting the statistical error as sigma. From the final stacked distribution at each radius, we then computed the median and associated errors as the 50th and [16th–84th] percentiles, respectively. To derive the intrinsic scatter of the entropy measurements as a function of radius, we fitted the interpolated data points at each radial bin in a Bayesian framework, using PyMC v.5.6.1 (Abril-Pla et al. 2023), with the model:

$$\begin{array}{c}\hfill K/K_{500}^{Y_{\mathrm{X}}}=A·\exp (\pm {\sigma }_{\mathrm{int}}),\end{array}$$(3)

where A is a constant and σ_int the intrinsic scatter. We note that the best-fitting values of A and σ_int are different at each considered bin. We report our results in Table 2 and plot the median profile both in the central right panel of Fig. 4 and also in Fig. C1, where we compare it with other scalings adopted throughout the paper.

4.2. Analytical fits

To provide further information on the properties of HIGHMz entropy profiles, including their normalisation, slope and intrinsic scatter at all radii, we adopted a constant-plus-power-law model, as first introduced by Donahue et al. (2005), together with a radially dependent scatter (Eq. (5) of Ghirardini et al. 2019). This parameterisation is able to capture both the power-law increase of the profiles at large radii and the flattening observed in the central regions. By using PyMC, we jointly fitted the scaled entropy profiles in the radial range [0.01 − 1] R₅₀₀^Y_X, with the model:

$$\begin{array}{c}\hfill K(x)/K_{500}^{Y_{\mathrm{X}}}=(K_0+K_1{x}^{\alpha })·\exp [\pm {\sigma }_{\mathrm{int}}(x)],\end{array}$$(4)

where x = R/R₅₀₀^Y_X and K₀ is the central entropy excess above the power law, with normalisation K₁ and slope α, that is measured at large radii; σ_int(x) is the intrinsic scatter⁴ of the profiles and varies with the radius following the quadratic functional form:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{int}}(x)={\sigma }_0+{\sigma }_1{\log }^2(x/x_0),\end{array}$$(5)

with σ₁ the width of the log-parabola, and x₀ and σ₀ the location and the intercept of the minimum of the log-parabola, respectively. The adopted priors and the resulting best-fitting parameters are reported in Table 3 (top); in Fig. 5, we show our best-fitting curve (top), with particular focus on its normalisation (centre) and slope (bottom), in comparison with predictions from Voit et al. (2005).

Fig. 5. Joint fit of the HIGHMz entropy profiles using the constant-plus-power-law parameterisation (Eqs. (4) and (5)). Top panel: best-fitting model (black), super-imposed to observational measurements (grey dots). Black dashed lines mark the intrinsic scatter, while red dash-dotted line is the prediction from non-radiative simulations (Voit et al. 2005). Central panel: ratio between our best-fitting model and predictions from non-radiative simulations. Grey shaded area is the associated statistical error. Bottom panel: local slope of the best-fitting analytical model (black), together with the associated statistical error (grey). In blue are the slopes obtained using fitting with piece-wise power laws (Eq. (6)). Red dash-dotted line is the canonical 1.1 slope.

Within ∼0.1 R₅₀₀^Y_X the best-fitting analytical profile becomes progressively flatter and the measured dispersion is large, as already discussed; moving out to the cluster outskirts, the intrinsic scatter is reduced and the best-fitting profile approaches the predictions from pure gravitational collapse (Voit et al. 2005), albeit with a gentler slope at all radii (α ∼ 0.9 at R₅₀₀^Y_X). The measured best-fitting profile is above the theoretical predictions for R ≲ 0.7 R₅₀₀^Y_X, while it becomes lower at larger radii. In particular, the measured gas entropy is ∼90% of that predicted by Voit et al. (2005) at R₅₀₀^Y_X.

Due to the rigidity of the adopted analytical model, any possible variations of the best-fitting slope with the radius would be smoothed away in favour of a mean slope. To address this issue, we also tested the fit with piece-wise power laws (as described in Eq. (4) of Ghirardini et al. 2019), which allow for greater flexibility. We divided the profiles into four radial bins, in the range [0.1 − 1] R₅₀₀^Y_X, and used PyMC again to build the model:

$$\begin{array}{c}\hfill K(x)/K_{500}^{Y_{\mathrm{X}}}=A·x^B·\exp (\pm {\sigma }_{\mathrm{int}}),\end{array}$$(6)

where A is the normalisation of the power law, with slope B, and σ_int is the intrinsic scatter in each bin. The priors and best-fitting parameters are reported in Table 3 (bottom), while the measured slopes are also shown in Fig. 5 (bottom). We do not find significant variations of the slopes with radius, which are in agreement with results from the global analytical fit. The measured piece-wise power-law slopes are flatter than the canonical 1.1 value at all radii.

4.3. Comparison with other observational samples

We compared HIGHMz entropy profiles to those of other observational samples found in the literature, namely REXCESS, ESZ, and X-COP (see Sect. 2.2). We stress that a blind comparison of the products obtained from different samples, analysed with distinct techniques, is likely to be affected by the distinct systematic uncertainties of each sample results. In the following, we take this into account and perform two separate comparisons, bringing the four samples to conditions that are as similar as possible:

• the entropy profiles of X-COP clusters were derived using azimuthal median densities and rescaled (for K₅₀₀ and R₅₀₀) using hydrostatic-equilibrium masses (Eckert et al. 2022). Consequently, we compare these data with HIGHMz entropies derived and scaled in a similar way;

• similarly, REXCESS and ESZ profiles were produced using azimuthal mean densities and rescaled using M₅₀₀^Y_X (Pratt et al. 2010). For comparison, we thus consider HIGHMz profiles derived from mean densities and adopt masses derived from the Y_X proxy in the rescaling. These profiles are also shown in Fig. C2.

Also for X-COP, REXCESS, and ESZ, we adopt a rescaling by K₅₀₀ following Eq. (1), using f_b = 0.16, as already done for HIGHMz clusters. For each of the four samples, we computed the median entropy profile and the intrinsic scatter following the method outlined in Sect. 4.1. The results are illustrated in Fig. 6, where we compare HIGHMz to X-COP (left), and to REXCESS and ESZ (right). The two HIGHMz median profiles shown in Fig. 6 are also compared in Fig. C1 to the ‘reference’ measurement, that is, entropy using median densities and M₅₀₀^Y_X in the rescaling.

Fig. 6. Comparison with other observational samples. On the left, we compare our results to X-COP (red squares), while on the right to REXCESS (yellow triangles) and ESZ (green stars) clusters. Top: median entropy profiles, with their intrinsic scatter. Centre: ratio to the HIGHMz median profile. The grey shaded area is the intrinsic scatter of the HIGHMz clusters, while the horizontal dotted lines indicate a 20% variation. Bottom: radial profiles of intrinsic scatter.

The median scaled entropy profile of HIGHMz consistently lies below those of the other samples, regardless of whether entropies are derived from azimuthal median (left) or mean (right) densities, and regardless the adopted mass in the scaling. The only exception is given by the core region, where we measure flatter entropies. This may be due to a combination of different effects, such as a different morphological distribution of HIGHMz clusters and/or a distinct resolution of the profiles due to the different redshift range. Further discussion on the distribution of the central entropies is provided in Sect. 5.2.1.

Beyond ∼0.1 R₅₀₀, the median profiles of the four samples exhibit similar shapes, as shown by the nearly constant ratios (sample/HIGHMz) with radius in the central panels of Fig. 6. The good agreement in shape with X-COP suggests that restricting the analysis to temperature measurements with S/B > 0.2 for HIGHMz has a comparable impact on the entropy profiles as a joint X/SZ analysis, as also pointed out by Rossetti et al. (2024) for the temperatures. Conversely, larger discrepancies are found when comparing normalisations. REXCESS, X-COP, and ESZ median profiles are, on average, ∼42%, ∼14% and ∼12% higher than those of HIGHMz, respectively. These differences in normalisation correlate with the median masses of the samples (Fig. 1): low-mass objects feature higher scaled entropy profiles than massive systems, which more closely approach the predictions from non-radiative simulations, consistent with the idea that non-gravitational effects play a more important role at the low-mass end. We discuss further the departures from self-similar predictions in Sect. 5.3.

In Fig. 6 (bottom), we also show the dispersion of the entropy profiles for the four samples as a function of the radius. In general, the measured intrinsic scatter profiles are in agreement at all radii. We measure considerable dispersion in the central regions, reflecting the morphological variety of the clusters, while the scatter reduces towards the external regions, where clusters appear to be more self-similar. A small increase of the intrinsic scatter is then observed at radii R ≳ 0.8 R₅₀₀, especially using the azimuthal mean, likely due to the larger impact of undetected sub-clumps. The only notable exception regards REXCESS clusters, which display a larger scatter (factor of ∼1.5 with respect to ESZ and HIGHMz) in the radial range [0.1 − 0.4] R₅₀₀, due to the wide range of masses that are present in the sample. To test this hypothesis, we divided REXCESS into two equal parts based on cluster masses and verified that the measured scatter of the more massive half is consistent with the ESZ and HIGHMz values.

4.4. Comparison with simulations

We now compare the HIGHMz entropy profiles with those of clusters taken from the MACSIS and The300 datasets (see Sect. 2.3), to test whether our theoretical understanding of the physical processes acting at the galaxy cluster scale can reproduce the properties of the observed entropy profiles.

As specified in Sect. 2.3, we assumed that Planck masses present an average 20% hydrostatic mass bias. We therefore consider HIGHMz masses corrected accordingly (i.e. M₅₀₀^Y_SZ/(1 − b), with 1 − b = 0.8), to compute both K₅₀₀ and R₅₀₀, which we then used in the rescaling. To ensure a fair comparison, we also adopted the scaling for K₅₀₀ reported in Eq. (1) for the simulations, using f_b = 0.16, as done for observations. We note that we have used mass-weighted entropy profiles for the simulated clusters, to reduce the effect of cold and dense sub-structures that are not masked in the simulations. However, in Appendix A we also present spectroscopic-like (Mazzotta et al. 2004) entropy profiles, as compared to the mass-weighted ones (Fig. A2).

4.4.1. Shape and dispersion

Figure 7 compares the scaled entropy profiles of the HIGHMz, MACSIS, and The300 samples. In general, the simulations reproduce the observed dispersion in the central regions and also the steady increase of ICM entropy with radius, with reduced scatter. Beyond the core, both MACSIS and The300 profiles increase according to a power law, with almost the same slope predicted from pure gravitational collapse (Voit et al. 2005). For a more quantitative comparison, we measured the median and the intrinsic scatter of the entropy profiles as a function of radius, following the same procedure detailed in Sect. 4.1.

Fig. 7. Comparison with simulated datasets. From the top to the bottom: (i) measured profiles for MACSIS (light blue), The300 (violet), and HIGHMz, rescaled to mimic a 20% hydrostatic bias, as described in the text (black); (ii) median entropy profiles of the three samples, together with their intrinsic scatters and the red dash-dotted line showing the prediction from Voit et al. (2005); (iii) ratio to HIGHMz, with the grey shaded area showing the intrinsic dispersion of HIGHMz entropies and horizontal dotted lines marking a 20% discrepancy; (iv) radial profile of the intrinsic dispersion for the three samples.

Beyond the core, we find a good agreement in shape with the median entropy profile of both simulated datasets, although HIGHMz profiles are flatter in the outskirts. Regarding the normalisation, HIGHMz median profile lies between those of MACSIS and The300, which slightly overestimate (∼6%) and underestimate (∼9%) the observed one, while still remaining within the measured dispersion at all radii. These differences are due to the combination of differences in temperatures and electron densities, as shown in Appendix A. Looking at the central regions (R ≲ 0.1 R₅₀₀), larger discrepancies are observed for MACSIS clusters, which predict strong cooling in the core; as a result, the simulated galaxy clusters have low central temperatures and peaked central densities (see Fig. A.1), which together give a systematically lower gas entropy than observed. This is not a new finding, as a similar behaviour has already been observed in previous works comparing simulations with observations, (e.g. Barnes et al. 2017) and also a recent study of FLAMINGO clusters (Braspenning et al. 2024), based on a similar code. Conversely, somewhat better agreement with observations is observed in the core of The300 clusters, probably favoured by the implemented artificial thermal diffusion term (Rasia et al. 2015), which allows for more efficient gas mixing in the central regions.

Regarding the dispersion of the profiles, presented in Fig. 7 (bottom), the simulations are able to reproduce the general behaviour observed in real clusters, albeit with some differences. For example, both simulated datasets overestimate the observed intrinsic dispersion in the central regions, although the difference is not statistically significant for The300, and underestimate it in the outskirts. At intermediate radii, i.e. [0.2 − 0.5] R₅₀₀, both simulated datasets reproduce well the measured dispersion for HIGHMz clusters.

4.4.2. Correlation with morphological parameters

In Fig. 4, we have highlighted that the central dispersion of HIGHMz entropy profiles is closely related to their dynamical state, as supported by the strong correlation between entropy and morphological indicators, such as c and w, and that this dependence is then reduced outside the core. In the following, we test the ability of simulations to reproduce the dependence of the entropy profiles on c and w. We approach this question in two different ways: by dividing the samples into disturbed and relaxed objects and comparing their median profiles; and by quantifying the correlation between entropy and morphological indicators as a function of radius.

We first divided HIGHMz, MACSIS, and The300 clusters into two sub-samples based on their classification as either most relaxed or most disturbed. This was achieved by sorting the clusters according to their centroid shift value, w, and selecting those with w values falling in the lower and upper tails of the distribution. Specifically, we identified the most relaxed and disturbed clusters as the 20% of the clusters with the lowest and highest w, respectively. This approach only considers clusters that are in the two tails of the distribution, thus reducing the impact of potential systematics arising from differences in the calculation of w between observations and simulations. We then computed the median entropy profiles of the two sub-samples and compared them with the median profile of all the clusters in HIGHMz, MACSIS, and The300 (Fig. 8). We find that the simulations reproduce the general behaviour seen in observations, i.e. relaxed clusters have steeper profiles at all radii than perturbed clusters. The difference between the median profiles of the two sub-samples is particularly pronounced in the central regions, where disturbed systems show higher entropy profiles than relaxed ones; the dichotomy disappears at ∼0.3 − 0.4 R₅₀₀ and is then slightly reversed in the outer regions, where we observe relaxed clusters more closely approaching predictions from Voit et al. (2005). In the outskirts, the simulations underestimate the differences between the two sub-samples, which is not surprising given the lower reproduced scatter in these regions with respect to the observations, as already noted in Sect. 4.4.1.

Fig. 8. Comparison between relaxed and disturbed systems, for HIGHMz (left), The300 (centre), and MACSIS (right). Green is used for the most relaxed clusters, while yellow for most disturbed ones. Black is used for the median profile obtained using all the clusters in each sample. In the bottom panels, we show the ratios to the medians of the full samples, with the dotted lines marking a 20% discrepancy.

We also calculated a more quantitative comparison of the dependence of the scaled entropy on the morphological parameters for the three samples. Specifically, we interpolated the entropy profiles on a dense radial grid, and for each radius we calculated (in the log-log plane) the correlation between the ICM entropy and both c and w, as quantified by the Pearson’s ρ parameter. For HIGHMz, we included information on the statistical errors in the measurements by performing a Monte Carlo simulation: for each radius, we generated 1000 random values (around the nominal measurements, with the sigma given by the statistical errors) for K/K₅₀₀, c and w, and for each realisation we calculated the correlation. Finally, from these distributions of measured Pearson’s ρ, we computed the medians and associated errors. Plotting the results as a function of radius as in Fig. 9 (left for c and right for w), it is possible to investigate how the dependence varies from the centre to the outskirts. Figure 9 reflects some of the aspects already discussed in Sect. 4.1: in the centre we observe a strong (anti-)correlation with (c) w for HIGHMz clusters, with ρ ∼ ( − 0.8) 0.6. The dependence on c and w reaches a minimum between R ∼ 0.2 − 0.3 R₅₀₀, and finally undergoes a slight inversion at larger radii. The simulations reproduce both the shape and intensity of these correlations at all radii, indicating that they can accurately replicate the spread of the entropy profiles and the dependence on dynamical state. The only difference is found at ∼R₅₀₀ for The300 where no significant correlation is found. This is consistent with the finding of Fig. 8 where the profiles of the relaxed and disturbed sub-samples of these 25 clusters from The300 are almost coincident at R₅₀₀.

Fig. 9. Pearson’s correlation coefficient between c (left) or w (right) and scaled entropy K/K500 as a function of cluster radius. Light blue, violet and black are used for MACSIS, The300, and HIGHMz, respectively. The grey shaded area is the statistical error associated with HIGHMz measurements.

5. Discussion

5.1. Flattening of the profiles in the cluster outskirts?

The shape of the entropy profiles in the cluster outskirts has been the subject of debate during the past decades, as conflicting results have been derived. According to cosmological simulations which implement gravitational processes only, entropy should exhibit a steady increase with radius (∝R^1.1) at least out to ∼2 R₂₀₀ (Tozzi & Norman 2001; Voit et al. 2005). However, observations exploiting the low particle background of Suzaku have indicated significant flattening at ∼0.6 R₂₀₀ and even at inner radii (e.g. Walker et al. 2012; Urban et al. 2014; Simionescu et al. 2017). Various explanations have been discussed to account for the observed flattening in Suzaku data, including differences between electron and ion temperatures in the outer ICM regions (Hoshino et al. 2010; Akamatsu et al. 2011) and a weakening of the accretion shock as it expands (Lapi et al. 2010; Cavaliere et al. 2011; Walker et al. 2012). Constraining the radius at which the flattening starts, and the extent of it, would represent a significant step forwards in our understanding of the physical mechanisms occurring at the accretion shock.

New insights were provided by Ghirardini et al. (2019), who reconstructed the entropy profiles of X-COP clusters, revealing a steady increase with radius out to at least ∼R₂₀₀, consistent with predictions from simulations. X-COP introduced several methodological improvements, such as the use of median instead of mean densities. In particular, the ability of median densities to exclude clumps down to scales of ∼10 − 20 kpc in the outer regions has been suggested as the main reason for discrepancies with Suzaku. The mean densities, combined with the low resolution of Suzaku (∼2 arcmin), would have prevented the exclusion of cool, over-dense structures, leading to a bias towards higher gas density and lower temperature values, and consequently resulting in an underestimation of the measured gas entropy.

With the high data quality of the CHEX-MATE clusters, we are now able to provide further conclusions about the shape of the entropy profiles in the cluster outskirts. We showed that flattening is occasionally observed at ∼R₅₀₀, both using azimuthal median (Fig. 4) and mean (Fig. C2) densities, when all the temperature bins are used in the reconstruction of the gas entropy. Conversely, correcting for the bias highlighted by Rossetti et al. (2024), i.e. excluding temperature measurements with S/B < 0.2, the entropy profiles resumed a steadily increasing trend, consistent with X-COP and with predictions from cosmological simulations. This suggests that the contribution of azimuthal median densities alone may not fully explain the differences with Suzaku, as it was suggested by Ghirardini et al. (2019), but a correction of a potential bias in the temperature measurement is also needed.

We investigate the impact of alleviating the bias in gas density (using the azimuthal median) and temperature (imposing S/B > 0.2) separately in Fig. 10. More specifically, we performed the ratio between the median entropy profiles of the sample reconstructed in two different ways:

• using azimuthal median and mean densities, with no correction for the bias in temperature (grey). This highlights the possible impact of a more efficient exclusion of cold and over-dense regions from the reconstruction of the gas density, while temperatures are left untouched;

• using azimuthal mean density profiles, combined with temperature profiles with and without low-S/B measurements (purple). This shows the impact of correcting the bias in temperature only.

Fig. 10. Relative contribution of potential density and temperature-related systematic uncertainties to the median entropy profile. We show the ratios between the medians $\stackrel{\sim }{K}$ of the scaled entropy profiles, reconstructed using: (i) median and mean densities, without cut for S/B > 0.2 (grey); (ii) mean densities, with and without cut for S/B > 0.2 (purple) and (iii) a combination of the two (red), as detailed in the text. Shaded areas are the 1σ statistical errors.

In Fig. 10, we show the ratios of the medians of the measured entropy profiles, together with their statistical errors.

The use of azimuthal median instead of mean density profiles results in a marginal increase (≲5%) of the ICM entropy at all radii, with a measured trend similar to the one shown by Eckert et al. (2015) and Bartalucci et al. (2023). Conversely, excluding temperature bins where S/B < 0.2 barely affects the reconstructed entropy profiles within ∼0.6 R₅₀₀ (although some differences are noted, due to deprojection effects), while it has larger impact in the outer regions (∼8% at R₅₀₀), especially for those clusters with significantly biased low external temperature measurements. This clearly indicates that, in the outskirts, correcting for the bias in temperature has a comparable or perhaps even larger impact than the use of azimuthal median densities.

Finally, we also show the joint effect of correcting densities (through the use of azimuthal medians) and temperatures (imposing S/B > 0.2) on the entropy profiles (red in Fig. 10). To first order, this is the net sum of the previous two measurements. We notice a radial increase of the median ratio (reaching ∼1.13 at R₅₀₀), meaning that the contribution from the systematics is larger in the outskirts than in the central regions. If these effects were not taken into account, they would have significant impact on the reconstruction of the entropy profiles, thus leading to incorrect conclusions on their shape in the cluster outskirts.

5.2. Shape of the profiles

In this section we present the results of the analytical fit to the individual HIGHMz entropy profiles using a constant-plus-power-law model of the form K = K₀ + K₁(r/100 kpc)^α (Donahue et al. 2005). This exercise allows us to discuss both the central entropy (K₀) distribution – a topic that has been debated in the past decades – and the distribution of the outer slopes α; in addition, it allows us to compare them with the results from previous works (Cavagnolo et al. 2009; Pratt et al. 2010).

We performed the fits in a Bayesian framework using PyMC, assuming flat priors in a wide range for the three parameters to be determined. The best-fitting parameters for each cluster are given in Table 4. To compare fairly with Cavagnolo et al. (2009) and Pratt et al. (2010), we used the HIGHMz entropy profiles derived from mean densities (see Appendix C). The histograms of the resulting central entropies K₀ and outer slopes α are shown in Fig. 11 (top left and top right), together with the cumulative distributions (bottom panels), which are independent of the adopted binning. In constructing the histograms, we took into account the statistical uncertainties on the best-fitting parameters, in order to correctly weight each measurement. More specifically, for each measurement of K₀ and α, we generated 1000 random realisations, normally distributed around the best-fitting value and using the statistical error as sigma, and then computed the histograms from these distributions, normalised to the total number of realisations. The histograms in Fig. 11 are colour coded according to the median centroid shift of the clusters in each bin, in order to link the shape of the profiles to their morphological state.

Fig. 11. Distribution of central entropies K0 (left) and outer slopes α (right) of the HIGHMz clusters. Plot bars are colour coded according to the median centroid shift w in each bin. Bottom panels are the cumulative distributions of the two histrograms. Blue and golden vertical lines mark the positions of the peaks identified with ACCEPT (Cavagnolo et al. 2009) and REXCESS (Pratt et al. 2010), respectively. The red dotted line on the right is the canonical slope value of 1.1.

5.2.1. Central entropy distribution

Based on an archival collection of 239 clusters observed with Chandra, known as ACCEPT, Cavagnolo et al. (2009) found a bimodal distribution of the central entropies K₀, with two peaks of similar amplitude located at K₀ ∼ 15 keV cm² and K₀ ∼ 150 keV cm². Conversely, by studying REXCESS clusters, Pratt et al. (2010) identified two tentative peaks at lower entropy (K₀ ∼ 3 and ∼75 keV cm²), with an amplitude ratio of 1 : 3. However, they could not statistically distinguish between a bimodal and a left-skewed distribution of K₀. Pratt et al. (2010) proposed various reasons for the differences observed, including a possible overestimation of the central temperature distribution by Cavagnolo et al. (2009) due to the lack of deprojection of their measured 2D temperatures. Moreover, since ACCEPT is an archive-limited sample, the clear gap identified by Cavagnolo et al. (2009) between K₀ ∼ 30 − 50 keV cm² may be the simple result of scientists’ prevailing interest in either strong-cool-cores clusters or in energetic mergers, leaving out clusters with ‘average’ properties, which will be observed in a representative sample like REXCESS. In summary, what the real distribution of central entropies in the cluster population is is still under debate. Clearly, given its close connection to the underlying physical mechanisms operating in cluster cores, the investigation of the central entropy distribution is a key issue for the cluster community.

We further investigated the topic of cluster central entropy distribution with the HIGHMz clusters. The results are shown in Fig. 11 (left). Unlike Cavagnolo et al. (2009), and more similarly to Pratt et al. (2010), we do not find evidence for a bimodal distribution of central entropies, but rather for a left-skewed distribution, which peaks at K₀ ∼ 200 keV cm², with an extended tail towards low K₀ values. This is not an artefact of the adopted binning, as supported by the cumulative distribution in the bottom panel of Fig. 11 (left), which shows a gradual increase with K₀. This finding indicates that we observe more clusters with high central entropy in the present sample than in REXCESS and ACCEPT. Similar behaviour has been noted in Sect. 4.3 in comparison to X-COP and ESZ, where the median entropy profile of HIGHMz has been found to be flatter in the central regions (Fig. 6). The small peak observed at ∼7 keV cm² is simply the result of two clusters with similar central entropy (PSZ2G073.97−27.82 and PSZ2G340.94+35.07). We note that we restricted the histogram to values of K > 0.1 keV cm²; however, four clusters have a best-fitting central entropy that is consistent with zero within 1σ.

Finally, we observe a clear correlation between K₀ and the centroid shift w, as indicated by the colour coding adopted in Fig. 11. Indeed, morphologically disturbed clusters tend to have higher central entropy, while morphologically relaxed clusters populate the left tail of the distribution, as already observed in Fig. 4. Typically, K₀ is used to classify cool-core and non-cool-core clusters from a thermodynamical point of view; although it indicates the morphological state of a cluster as measured at large scales, w shows a good correlation with K₀. As the median mass of HIGHMz is larger than both X-COP and ESZ, it may be that this very high-mass population is dominated by disturbed objects with a high central entropy. Alternatively, the pure SZ-selection of HIGHMz may contribute to the presence of a greater number of disturbed, high central entropy systems (Planck Collaboration IX 2011; Rossetti et al. 2016, 2017; Andrade-Santos et al. 2017; Lovisari et al. 2017).

Although our results seem to be in contrast with previous findings, they cannot be considered as definitive. In fact, the study of central regions is not straightforward, and the presence of additional uncertainties may bias our conclusions. For example, the physical size of the central bins is different for morphologically relaxed and disturbed clusters, a consequence of the adaptive binning method we used for the spectral extraction (Chen et al. 2024; see Sect. 3.2 of Rossetti et al. 2024), thus possibly altering to some level the measured distribution of the central entropies. In addition, differences in sample selection methods (X-rays and SZ), the adopted satellite for the study (Chandra and XMM-Newton) and redshift ranges may prevent a fair comparison between the samples. Future studies on the entire CHEX-MATE sample will provide increased statistical power to derive additional constraints on the central entropy distribution in clusters. However, while CHEX-MATE was designed to ensure homogeneous observations at R₅₀₀, the same cannot be guaranteed for the central regions. Therefore, a dedicated programme on a representative cluster sample, with the aim of reaching homogeneous spatial resolution in the central regions, may be needed to definitively shed light on this issue.

5.2.2. External slopes of the profiles

The right-hand panel of Fig. 11 shows the distribution of outer slopes, α, of the HIGHMz clusters, as from the constant plus power-law model. A wide range of values for the parameter α is observed, extending from ∼0.5 to ∼2. Similarly to the findings by Pratt et al. (2010) for REXCESS clusters, there is no indication of bimodality in the distribution of slopes. The median of the measured distribution is α ∼ 1.12, remarkably close to the self-similar nominal value of 1.1, and slightly larger than 0.98, as measured by Pratt et al. (2010), shown by the golden vertical line in Fig. 11 (right) using REXCESS. While the median slope of the REXCESS sample still falls within the large 1σ dispersion of HIGHMz (σ = 0.35), the lower median value likely reflects the presence of more lower-mass systems in the REXCESS sample. Finally, the peak value measured by individual fits of the entropy profiles (α ∼ 1.12) is larger than the best fitting value presented in Fig. 5, when all the HIGHMz entropy profiles are jointly fitted (α ∼ 0.87). Although the latter is still included within the 1σ distribution presented in this section, the observed differences may simply reflect the rigidity of the assumed model, which allows for larger variability when fitting individual profiles.

Similarly to the central entropies, the distribution of the outer slopes presented in Fig. 11 (right) is also colour coded according to the median centroid shift parameter (w) of each bin. This allows us to link the measured outer slopes to the morphological state of the clusters. From the figure, it is evident that more relaxed clusters, with lower values of w, exhibit a slope α much closer to the canonical 1.1 value, with a relatively small scatter. Conversely, more disturbed clusters, with larger values of w, are more prevalent in the tails of the distribution, displaying both lower and higher values than 1.1. This is also in line with the findings presented in Fig. 5 of Pratt et al. (2010), showing that the best fitting slopes of the sub-sample of relaxed REXCESS clusters is closer to ∼1 with a relatively small scatter, while dynamical disturbed objects are characterised by larger scatter in the α distribution.

5.3. Beyond the self-similar scenario

In Sect. 4.3, we compared entropy profiles of ESZ, X-COP, and REXCESS clusters to those from HIGHMz, noting similar shapes but significant differences for their normalisations. These discrepancies were observed to correlate with the median masses of each sample, with low-mass clusters exhibiting more pronounced deviations. These findings support the idea that non-gravitational processes leave a signature in the physical properties of the ICM, particularly so for less massive systems, and so a self-similar scaling (Eq. (1)) may not be fully representative of all the physical processes affecting the gas entropy in a cluster.

Some works already focused on alternative rescalings of the thermodynamic profiles, that go beyond the self-similar scenario and account for the effects of non-gravitational processes. For example, Pratt et al. (2010) showed a mass-dependent excess of the gas entropy in REXCESS clusters, while Pratt et al. (2022) highlighted a stronger than self-similar evolution with redshift, together with a significant residual dependence on mass, for the density profiles of ∼120 galaxy clusters. In the following, we investigate further the impact of non-gravitational processes through the combined sample of REXCESS, X-COP, ESZ, and HIGHMz clusters, to explore the possibility of mass and redshift dependencies of the gas entropy which go beyond the self-similar scenario.

We initially studied the impact of both an additional mass dependence and a modified evolution with respect to self-similar expectations. This investigation is detailed in Appendix B, where we show that we did not find strong statistical evidence for a modified redshift dependence in the scaling, likely due to the limited redshift range of the sample under consideration. We therefore assumed self-similar evolution and studied the impact of a residual mass dependence only. This is detailed in Sects. 5.3.1 and 5.3.2, where we investigate both global and radial effects, respectively. Finally, although we used the same scaling and adopted the same code to compute the medians of the different samples, some additional systematic effects, for example due to different data analysis approaches, may affect the following results to some degree. Through the future study on the entire CHEX-MATE sample, we will be able to derive definitive constraints on the departure from the self-similar scenario, minimising potential systematics through a homogeneous analysis of all the 118 clusters in the sample.

5.3.1. Global dependence on mass

We assumed that the differences between the median entropy profiles of the four samples can be entirely explained by a modified dependence on cluster mass, rather than a self-similar one. We thus introduced the parameter A_m, that quantifies the departure from self-similar predictions, and built a modified entropy rescaling:

$$\begin{array}{c}\hfill \stackrel{\sim }{K_{500}}=M_{500}^{2/3+A_m}E{(z)}^{-2/3}.\end{array}$$(7)

To identify the value of A_m that is needed to explain the differences we have observed in Sect. 4.3, we studied how the scatter of the measured profiles varies as a function of A_m. Our best-fitting residual mass dependence is then the value of A_m that minimises the dispersion of the profiles in a given radial range.

To do so, we followed a procedure similar to the one described in Sect. 4.2. Specifically, we fitted together the entropy profiles of the considered samples in the radial range [0.3 − 0.8] R₅₀₀, thus excluding the flattening that is observed in the core region and potentially biased measurements in the outskirts, using a power law and an intrinsic scatter:

$$\begin{array}{c}\hfill K(x)/\stackrel{\sim }{K_{500}}=A·x^{\alpha }·\exp (\pm {\sigma }_{\mathrm{int}}),\end{array}$$(8)

where again x = R/R₅₀₀. To investigate the dependence on the residual mass dependence, we generated 250 random values of A_m, uniformly distributed in the range [ − 0.75, +0.5]. At each iteration, the entropy profiles were rescaled using a different realisation of $\stackrel{\sim }{K_{500}}$ and the best-fitting intrinsic scatter was then plotted as a function of A_m. Similarly to what was discussed in Sect. 4.3, we specify that we considered X-COP and HIGHMz and REXCESS, ESZ, and HIGHMz separately in the analysis, to reduce the impact of systematic errors coming from a heterogeneous cluster sample. We also noticed that some clusters are in common between HIGHMz, ESZ, and REXCESS; we thus excised duplicates that may bias our results, keeping our measurements first and then, if a cluster was only in common between REXCESS and ESZ, those reported in Pratt et al. (2010). The final combined samples count 44 clusters for the combined X-COP – HIGHMz sample, and 108 for REXCESS – ESZ – HIGHMz.

The results are presented in Fig. 12 (left), where we plot the distribution of the measured dispersions as a function of the residual mass dependence A_m. We immediately notice that these describe well defined curves in the plane A_m − σ_int and that the minima of these curves are not achieved at A_m = 0, as would be the case if cluster entropies scaled self-similarly. The exact positions of the minima and the associated errors were identified through a Monte Carlo simulation; specifically, (i) we generated 1000 random realisations of each data point around its statistical error; (ii) each time we linearly interpolated the randomly generated curve and found its minimum; and (iii) we computed the median and 1σ dispersion from this distribution. The measured errors give an idea of the flattening of each curve around the position of the minimum. We measure A_m (min. scatter) = − 0.24 ± 0.03 and −0.21 ± 0.04 for the combined samples of REXCESS – ESZ – HIGHMz (brown in Fig. 12, left) and X-COP – HIGHMz (yellow), respectively. We thus find evidence that a milder mass dependence than self-similar is needed to reduce the scatter of the measured entropy profiles to its minimum in the [0.3 − 0.8] R₅₀₀ radial range, supporting the idea that non-gravitational processes break the self-similarity of the observed clusters and cause an enhancement of their gas entropy that depends on cluster mass. We note that the measured values of A_m (min. scatter) for the two combined samples are in good agreement within their error bars. This is indeed expected, as REXCESS and ESZ together cover the mass range of the X-COP sample.

Fig. 12. Constraints to the residual mass dependence Am in the entropy rescaling. Left panel: distribution of the intrinsic scatter of different samples as a function of Am, to identify the global (i.e. within [0.3 − 0.8] R500) Am correction (Sect. 5.3.1). Yellow is used for the combined X-COP – HIGHMz sample, brown for REXCESS – ESZ – HIGHMz, while blue for HIGHMz only. Vertical dashed lines and shaded areas mark the positions of the minima, together with their statistical errors. Right panel: radial dependence of the Am values which minimise the scatter of the combined REXCESS – ESZ – HIGHMz sample (Sect. 5.3.2), together with their statistical errors (black solid line and grey shaded area). Shown in brown (blue) is the global measurements using REXCESS – ESZ – HIGHMz (HIGHMz only) sample. In red are results from REXCESS only (Pratt et al. 2010). Both in the left and right panels, black dashed lines are the prediction from self-similar scenario, i.e. Am = 0.

Finally, we also repeated the same exercise using the HIGHMz clusters only (blue in Fig. 12, left), to test whether the departure from the self-similar scaling is observed also at the massive cluster scale. We find A_m (min. scatter) = − 0.04 ± 0.10, consistent with zero within the large error bar. This supports the idea that, beyond the core, the properties of massive clusters are mainly driven by the gravitational assembly, which results in a good agreement with the self-similar predictions for the rescaling.

5.3.2. Radial dependence on mass correction

In the previous section, we identified a value A_m that would explain the global differences observed between samples at different masses. However, Pratt et al. (2010) also identified a radial trend of the mass-dependent entropy excess, that is stronger at ∼R₂₅₀₀ and becomes negligible at ∼R₅₀₀ within the large observational uncertainties. The increased sample size of the present study allows a similar investigation to be conducted with an improved precision, and thus to map a radial dependence of the measured A_m values that minimise the dispersion of the profiles. To do so, we repeated the procedure presented in Sect. 5.3.1 in a slightly different fashion. Specifically, we first interpolated the profiles on a common grid of 10 equally spaced radial bins within the range [0.1 − 1] R₅₀₀ (plus an additional measurement at 0.05 R₅₀₀), similarly to Sect. 4.3. We then fitted the data points belonging to each radial bin with the model:

$$\begin{array}{c}\hfill K/\stackrel{\sim }{K_{500}}=A·\exp (\pm {\sigma }_{\mathrm{int}}),\end{array}$$(9)

and studied the distribution of the intrinsic scatter as a function of A_m at each radius. Similarly to what was done in Sect. 5.3.1, we generated 100 random values of A_m to use in the rescaling, uniformly distributed in the range [ − 0.75, 0.25]. However, for this test we used only the combined sample of REXCESS – ESZ – HIGHMz, since the larger number of objects affords better precision on the radial dependence of the mass correction. The position of the minimum value of each distribution and the associated errors were computed as in Sect. 5.3.1. These are shown in Fig. 12 (right), as a function of radius and compared to previous measurements by Pratt et al. (2010).

We find a radial trend of the A_m values minimising the ten distributions, meaning that the minimum scatter is achieved at lower A_m values at small radii, while a weaker correction is needed moving towards the cluster outskirts. More specifically, we measure A_m (min. scatter) = − 0.41 ± 0.15 at 0.05 R₅₀₀^Y_X, while this value is closer to self-similar predictions at R₅₀₀^Y_X, where we find A_m (min. scatter) = − 0.14 ± 0.06. Our measurements are in agreement with previous findings by Pratt et al. (2010) and support the idea that non-gravitational processes have larger effects on the ICM properties at small radii, while corrections to self-similar predictions are smaller moving towards the outskirts. Interestingly, at R₅₀₀^Y_X their measurement is consistent with zero within the large error bars; with increased statistics, we are able to constrain a residual mass dependence other than zero also at this external radius.

Our measured value of A_m (min. scatter) at R₅₀₀^Y_X is in agreement with predictions based on the semi-analytic model proposed by Ettori et al. (2023), i.e. A_m (min. scatter) = − 0.16 ± 0.01, although they also considered a modified dependence on redshift, rather than assuming a self-similar evolution as we have done. At the same radius, our measurement is also consistent within 1σ with results from simulations implementing full baryonic physics, as presented in Nagai et al. (2007), which predict A_m = −0.062 ± 0.036.

5.4. Rescaling using weak-lensing masses

In Sect. 4.1, we have adopted masses derived from Y_X to build the rescaling of the measured entropy profiles, due to their lower sensitivity to cluster dynamical state compared to hydrostatic equilibrium masses. However, Pratt et al. (2022) have shown that the use of M₅₀₀^Y_X leads to a systematic suppression of the true dispersion of the density profiles, for a sample of ∼120 clusters. This happens as a cluster with density higher than the average (at a fixed mass) will possess greater gas mass, resulting in a higher measured Y_X value. Consequently, this will lead to an overestimate of M₅₀₀^Y_X (and so R₅₀₀^Y_X) that, when used in the rescaling, will move the measured density profile back to the mean one, thus reducing the apparent scatter. In other words, Pratt et al. (2022) showed that the covariance between M₅₀₀^Y_X and gas density implies a suppression of the intrinsic scatter in rescaled density profiles.

In this section we investigate further the same issue using the entropy profiles of HIGHMz. In particular, to test the degree of covariance between our profiles and M₅₀₀^Y_X, we performed an alternative rescaling, using masses derived from a ground-based weak lensing analysis (hereafter $M_{500}^{\mathrm{WL}}$). These masses are not estimated from any of the ICM observables and do not rely on any assumptions about the dynamical state of the clusters; therefore they are commonly considered a good estimate of the ‘true’ mass of a cluster in terms of accuracy (Becker & Kravtsov 2011; Pratt et al. 2019; Umetsu 2020; Euclid Collaboration 2024). However, an unbiased reconstruction of the weak-lensing signal (or reduced tangential shear profile) within the cluster extension is not straightforward and is particularly challenging for individual low-mass systems (M₅₀₀ ∼ 10¹⁴ M_⊙; see Fig. 17 of Umetsu et al. 2020) due to their shallower potential well, and for distant ones, as it is more difficult to have information on a sufficient number of background galaxies. As a result, weak-lensing masses are often accompanied with large statistical errors or, in other terms, they are nearly unbiased, but not precise.

For our purpose, we took advantage of weak-lensing M₅₀₀ mass estimates available in the Amalgam2 dataset based on archival Subaru Suprime-Cam and CFHT MegaPrime/MegaCam observations (Umetsu et al., in prep. and Gavazzi et al., in prep.), produced within the CHEX-MATE collaboration. These weak-lensing masses were inferred from the reduced tangential shear profiles of individual clusters assuming the spherical NFW profile (Navarro et al. 1996) with uninformative priors on the halo mass and concentration (M₂₀₀, c₂₀₀), closely following the procedures outlined in Umetsu et al. (2020). To date, these masses have been measured for 42 CHEX-MATE clusters including a subset of HIGHMz comprising 16 clusters, that we then considered for the present study; the entropy profiles of this sub-sample, scaled using $M_{500}^{\mathrm{WL}}$, are illustrated in Fig. 13 (top).

Fig. 13. Rescaling using weak-lensing masses. From top to the bottom: (i) entropy profiles of the sub-sample of 16 clusters with weak-lensing masses available, scaled using $M_{500}^{\mathrm{WL}}$, where shaded grey lines show the impact of accounting for statistical errors on $M_{500}^{\mathrm{WL}}$ in the rescaling, as described in the text; (ii) median entropy profiles, obtained using M500YX (red dotted line) and $M_{500}^{\mathrm{WL}}$ (blue dotted line) in the rescaling; continuous lines are the measurements after including statistical errors on mass measurements; (iii) same as centre, but for the intrinsic dispersion; here, simple scalings without including errors on mass measurements are shown as shaded areas.

To test the impact of the adopted mass in the rescaling on the measured dispersion of the profiles, we employed the same fitting procedure outlined in Sect. 4.1. Specifically, we adopted the model described in Eq. (3), both considering M₅₀₀^Y_X and $M_{500}^{\mathrm{WL}}$ to build K₅₀₀ and R₅₀₀, to test the impact of the two mass rescalings on the measured intrinsic scatter.

The measured medians with this first fitting procedure are shown in Fig. 13 (centre), as dotted red (M₅₀₀^Y_X) and blue ($M_{500}^{\mathrm{WL}}$) lines, while the intrinsic scatter profiles are plotted in the bottom panel, as shaded pink and light blue envelopes, respectively. The figure shows that using weak-lensing masses leads to a systematic decrease (a flat ∼10%) of the median profile at all radii and a net increase (≳50%) of the measured dispersion beyond 0.3 R₅₀₀. However, this might not be the result of the covariance between the entropy profiles and M₅₀₀^Y_X, but it may simply reflect the impact of the larger statistical errors associated with the $M_{500}^{\mathrm{WL}}$ measurements. In fact, the median fractional statistical uncertainty on weak-lensing masses is ∼28%⁵, while that of masses derived from Y_X is ∼3%.

To explore their impact on the final measurements, we propagated the statistical errors associated with mass measurements into the scaled entropy profiles. More specifically, we generated 1000 random realisations of M₅₀₀^Y_X and $M_{500}^{\mathrm{WL}}$ for each cluster, normally distributed around the measured values using the 68% dispersion as sigma, that we used to scale the profiles. Figure 13 (top) also illustrates the spread around the measured entropy profiles due to mass randomisation, when scaling using $M_{500}^{\mathrm{WL}}$. From the random distribution around each entropy bin, we computed (as 16th and 84th percentiles) a new statistical error associated with each measurement that quantifies the impact of mass randomisation. This additional error was then used to build an updated version of the error measurements at each bin:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{incl}.\mathrm{errors}}=\sqrt{{\sigma }_{\mathrm{stat}.}^2+{\sigma }_{\mathrm{ran}.\mathrm{masses}}^2},\end{array}$$(10)

including the original measured statistical errors (σ_stat.) and the spread due to mass randomisation (σ_{ran. masses}). Adopting Eq. (3), we then recalculated the median and the intrinsic scatter of the scaled profiles, this time including the information from mass randomisation. The results of these new fitting procedures are shown in Fig. 13 (centre and bottom) as red (M₅₀₀^Y_X) and blue ($M_{500}^{\mathrm{WL}}$) solid lines.

Accounting for the statistical errors associated with weak lensing mass measurements has small impact on the measured median (lower than 5% at all radii), while it systematically reduces the dispersion (average ∼15% beyond the core), reflecting the impact of their large uncertainties. Conversely, the small errors on M₅₀₀^Y_X barely affect the median and intrinsic scatter when employing this mass rescaling. A comparison between the two rescalings suggests that employing Y_X masses leads to a ∼30% suppression of the intrinsic scatter of the profiles beyond ∼0.2 R₅₀₀, although this result is not statistically significant. In the core, no evident differences are found. This would point to the conclusion that a degree of covariance of entropies and Y_X is indeed observed, and using weak-lensing masses would allow us to recover the ‘true’ entropy scaling. As already mentioned, a similar behaviour was also observed by Pratt et al. (2022); assuming that the suppression profile of their Fig. 10 is translated into a similar one for entropies, we would find a correction of ∼35% beyond the core, similar to the one reported here.

However, the limited size of the sample adopted for the present study precludes definitive conclusions on this issue. The future study on the entropy profiles for the full CHEX-MATE sample, together with a similar investigation performed on the other thermodynamic quantities, will allow more robust conclusions on the impact of different mass measurements on the reconstruction of the intrinsic scatter of the thermodynamic profiles in galaxy clusters.

6. Summary and conclusions

The present study focused on the entropy profiles of 32 galaxy clusters, with masses M₅₀₀^Y_SZ > 7.75 × 10¹⁴ M_⊙ and located at redshifts 0.2 < z < 0.6. This sample, named HIGHMz, was selected from the parent CHEX-MATE sample (CHEX-MATE Collaboration 2021) to be representative of the most massive cluster population. The high masses of HIGHMz clusters allowed us to investigate the cluster entropy in the gravity-dominated regime, where non-gravitational processes are expected to have only minor effects, thus setting a valuable baseline for future studies at lower-mass scales and higher redshifts. The ultimate goal of this work is to provide a set of observational constraints on the entropy profiles of the most massive clusters of the Universe, which will serve as a future reference for cosmological simulations and theoretical models describing entropy production and modification mechanisms at the galaxy cluster scale.

In general, the entropy profiles of HIGHMz clusters share similar properties to those measured in previous studies (e.g. Cavagnolo et al. 2009; Pratt et al. 2010; Ghirardini et al. 2019). For example, we found that:

• notable dispersion is present in the core regions, while the profiles scale more self-similarly and approach predictions from non-radiative simulations (Voit et al. 2005) towards the cluster outskirts;

• the central dispersion reflects the morphological variety of the clusters in the sample, as revealed by the strong correlation with morphological indicators; and

• compared with the predictions from non-radiative simulations, beyond the core, the profiles are well fitted by a power law with a gentler slope at all radii (α ∼ 0.9 at R₅₀₀; Fig. 5), leading to a lower (∼10%) measured entropy in the cluster outskirts.

In contrast with Cavagnolo et al. (2009) and more similarly to Pratt et al. (2010), we found no evidence of bimodality in the central entropy distribution (Fig. 11).

We investigated the relative role of possible systematic errors in the entropy reconstruction. We showed that:

• the exclusion of temperature measurements in radial bins with a low source-to-background (S/B < 0.2, Rossetti et al. 2024) from the deprojection, which lie mostly beyond R₅₀₀, allows us to regularise the entropy profiles in the outskirts. The change in the external part of the profiles is larger than 10% only for four clusters, which are those featuring a very low temperature in the last radial bin (Fig. 3). After exclusion of these temperature measurements, these profiles more closely approach the predictions from non-radiative simulations; and

• the use of azimuthal median densities, which allows a more efficient excision of dense clumps, impacts the final shape of the entropy profiles at all radii (Fig. 10), but in the cluster outskirts the temperature measurements are the major contributor to the systematic uncertainties.

We compared our profiles with other observational samples, namely REXCESS (Pratt et al. 2010), ESZ (Planck Collaboration XI 2011), and X-COP (Eckert et al. 2022). Entropy profiles are all in good agreement in shape, but some differences in normalisation are present, with the low-mass clusters showing higher scaled entropy at all radii (Fig. 6). These discrepancies correlate with the median masses of the four samples and are a clear indication of the impact that non-gravitational processes have on the properties of the ICM. Through the combined sample of HIGHMz, ESZ, and REXCESS, we investigated the possibility of an alternative scaling of entropy with mass (Eq. (7)) that can account for the break of the self-similarity due to non-gravitational processes. In particular, we found that:

• a weaker mass dependence than self-similar (A_m ∼ −0.25) is able to minimise the dispersion of the entropy profiles in the radial range [0.3 − 0.8] R₅₀₀ (Fig. 12, left);

• the effect is radial (Fig. 12, right), with stronger deviations from self-similarity at small radii, although a non-zero mass dependence was measured also at R₅₀₀, where we found A_m = −0.14 ± 0.06;

• the same study, performed on HIGHMz clusters only, showed that the A_m value minimising the scatter is consistent with being zero, that is, self-similar;

• a modified scaling is also possible for the evolution with redshift (Eq. (B.1); see Appendix B for more details), although the small considered redshift range did not allow definitive conclusions.

We also compared our results to recent simulations (The300, Cui et al. 2018; MACSIS, Barnes et al. 2017) that are able to produce a sufficient number of massive clusters at redshifts 0.2 < z < 0.6. We found the following:

• MACSIS clusters are in good agreement with observations, both in shape and in normalisation (Fig. 7), where differences lower than ∼6% are measured. Slightly larger differences (offset ∼ − 9%) are observed in comparison with The300.

• Regarding the correlation with morphological indicators, simulations are able to reproduce the measured behaviour for observations. We verified this both through sub-samples of defined relaxed and disturbed clusters (Fig. 8) and, more quantitatively, through a correlation of the gas entropy with c and w as a function of radius (Fig. 9).

Finally, we tested the rescaling of the entropy profiles adopting masses derived from a weak lensing analysis and produced within the CHEX-MATE collaboration, in order to reduce covariance (as highlighted by Pratt et al. 2022 for densities) and to recover the ‘true’ dispersion. We considered for this study a subset of 16 HIGHMz clusters, for which weak-lensing masses were already available. We found that:

• the different statistical errors on weak lensing and Y_X masses have an impact on the results, if they are not properly taken into account in the analysis; and

• after propagating the statistical errors on mass measurements, there is an indication that the rescaling through Y_X masses leads to a systematic suppression (∼30%, although not statistically significant) of the intrinsic dispersion of the profiles beyond the core (Fig. 13).

Although we managed to derive several observational constraints on the properties of the gas entropy at the massive cluster scale, we could not provide definitive conclusions on a few other aspects. For example, both the study of the central entropy distribution and the rescaling using weak-lensing masses would require a larger cluster sample size. Furthermore, as discussed in Sect. 4.3, we emphasise that in comparing our data with other observational samples, we tried to minimise systematic errors that may result from a heterogeneous sample, for example by using the same code to compute median profiles and adopting the same mass rescalings. However, some residual hidden systematics may still be present, affecting our results to some degree. Further constraints on these and more aspects will be possible through the future study of the entropy profiles of the entire CHEX-MATE sample.

In conclusion, our results have set an observational baseline for the radial entropy distribution in the highest mass systems, where gravitational infall is expected to be the dominant entropy generation mechanism. Outside the core regions, the profiles converge rapidly, both in slope and normalisation, to the theoretical expectation from gravity-only numerical simulations. The results further confirm the findings of Pratt et al. (2010), in the sense that high central entropy is correlated with the morphological appearance of a cluster, and that low central entropy is correlated with the presence of a cool core. In this paper, for the first time, we have confirmed these same trends in recent cosmological numerical simulations of objects of a similar mass. If ICM morphology is correlated to dynamical state (e.g. Bauer et al. 2005; Poole et al. 2006; Jeltema et al. 2008; Chon et al. 2012; Rasia et al. 2013; Lovisari et al. 2017; Campitiello et al. 2022), these results therefore indicate that a high central entropy at these high masses is caused by gas mixing due to merging processes.

Comparison of the HIGHMz sample with other samples from the literature at a lower average mass reveals the presence of additional entropy, above that expected from gravitational generation alone. A modification of the expected scaling to take this into account reveals a mass dependence that is significantly less than self-similar. In other words, the ICM in lower-mass systems has a higher entropy than that expected from gravitational collapse alone. Outside the core regions, this entropy increase is radially dependent, reaching a larger radius in lower-mass systems. This is consistent with a scenario where radiative cooling of the gas, or additional energy injection from AGN activity or supernovae, raises the ICM in lower-mass systems to a higher adiabat.

An important avenue for further work will therefore be to disentangle the non-gravitational energy injection from the effect of gas mixing due to merging activity. We intend to pursue this with the full CHEX-MATE sample.

Acknowledgments

We thank the anonymous referee for useful comments, which have improved our paper. This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. The research has received funding from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158). 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”. G.R. acknowledges C. Grillo for useful discussions. G.R., M.R., I.B., H.B., S.D.G., F.D.L., S.E., F.G., S.G., L.L., P.M. and S.M. acknowledge the financial contribution from the contracts Prin-MUR 2022, supported by Next Generation EU (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). G.W.P acknowledges long-term support from CNES, the French space agency. K.U. acknowledges support from the National Science and Technology Council of Taiwan (grant NSTC 112-2112-M-001-027-MY3) and the Academia Sinica Investigator Award (grant AS-IA-112-M04). L.L. also acknowledges support from INAF mini grant 1.05.12.04.01. M.S. acknowledges financial contributions 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) and INAF Theory Grant 2023: Gravitational lensing detection of matter distribution at galaxy cluster boundaries and beyond (1.05.23.06.17). H.B., P.M., and F.D.L. acknowledge financial contribution from the contracts ASI-INAF Athena 2019-27-HH.0, “Attività di Studio per la comunità scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASI-INAF n. 2017-14- H.0), support 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 CN00000013 ‘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). H.B., P.M. and F.D.L. also acknowledge the financial contribution from the contract PrinMUR 2022 supported by Next Generation EU (n.20227RNLY3 The concordance cosmological model: stress-tests with galaxy clusters). Ma.G. acknowledges support from the ERC Consolidator Grant BlackHoleWeather (101086804). B.J.M acknowledges support from Science and Technology Facilities Council grant ST/Y002008/1. E.P. acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-22-CE31-0010 (project BATMAN). J.S. was supported by NASA Astrophysics Data Analysis Program (ADAP) Grant 80NSSC21K1571. We made use of the following software 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: The300 and MACSIS datasets

We report in Tables A.1 and A.2 additional information on the selected The300 and MACSIS clusters, respectively, which were used for comparison in Sect. 4.4. In particular, we specify their masses (M₅₀₀, and so the derived R₅₀₀), used in the rescaling of the profiles, together with morphological indicators such as c and w, computed as described in Campitiello et al. (2022) and Towler et al. (2023), respectively.

In Sect. 4.4, we showed that the scaled entropy profiles of HIGHMz clusters lie between those of The300 and MACSIS clusters, in terms of their normalisation, albeit with some small systematic differences at all radii. To better understand the origin of these discrepancies, and in general to further investigate the comparison with the simulated datasets, we broke down the contributions of MACSIS and The300 temperature and electron density profiles; these profiles are shown in Fig. A.1, in comparison with HIGHMz. In order to compute the median and the intrinsic scatter of the profiles, we followed the same procedure outlined in Sect. 4.1. The temperature profiles were scaled using T₅₀₀ computed from Eq. 2 of Arnaud et al. (2005), for an overdensity Δ = 500 and T > 3.5 keV; densities are scaled self-similarly by E(z)².

Fig. A.1. Comparison of temperature (left) and density (right) profiles of MACSIS (light blue) and The300 (violet) clusters, with respect to HIGHMz ones (black). Upper, central and lower panels show the same information as in Fig. 7.

We notice that, beyond the core, the median temperature profile of the The300 systematically underestimates the observed one, while densities are overestimated. A similar, but opposite, behaviour is found for MACSIS. Given the way temperature and density combine to yield the entropy, the resulting entropies are underestimated for The300 and overestimated for MACSIS, as already observed in Fig. 7. When looking at the central regions (≲0.1 R₅₀₀), the differences between simulations and observations become more evident. MACSIS clusters exhibit, in median, a stronger cooling and a more peaked central density profile than observed in HIGHMz clusters, as also found in previous works comparing simulations to observations (Barnes et al. 2017; Braspenning et al. 2024), while The300 clusters are better at reproducing the central temperature distribution of observations. Once again, these discrepancies are the origin of the differences noted in Fig. 7 for the entropy profiles.

As a final consideration, we recall that in Sect. 4.4 we considered mass-weighted entropy profiles for the simulated clusters, in comparison with HIGHMz ones. Mazzotta et al. (2004) discussed an alternative weighting method of the particles in simulations, called spectroscopic-like, which weights more dense and cold sub-structures, as it may happen for real observations if these regions are not masked properly in the analysis. However, since sub-structures are not masked in simulations, the use of spectroscopic-like profiles may bias the reproduced temperatures towards excessively low values. For this reason, in Sect. 4.4, we preferred to keep mass-weighted profiles as our reference. However, we present in Fig. A2 a direct comparison of mass-weighted to spectroscopic-like temperature profiles. In general, spectroscopic-like temperature profiles are lower than mass-weighted ones at all radii, particularly so in the cluster outskirts, where the presence of in-falling sub-structures causes a net increase of the ratio T_mw/T_sl. Another consequence of that is also the increase of the intrinsic scatter for spectroscopic-like profiles, which becomes consistent with observations at ∼R₅₀₀, both for MACSIS and The300.

Fig. A.2. Comparison of mass-weighted (Tmw) and spectroscopic-like (Tsl) temperature profiles, for MACSIS (light blue) and The300 (pink). Plotted are the ratios of the two measurements for each cluster. The two more marked lines are the median ratios for the two samples.

Appendix B: Residual dependence on redshift

We further explore the possibility of an alternative rescaling of the entropy profiles that goes beyond self-similarity (Eq. 1). In particular, in this section we investigate a possible residual dependence on both mass and redshift, using the clusters in the combined HIGHMz – REXCESS – ESZ sample. Similarly to Sect 5.3, we introduced two parameters, A_m and A_z, which quantify the departure from self-similar predictions for the mass and redshift dependencies, respectively, and built the modified entropy rescaling:

$$\begin{array}{c}\hfill \stackrel{\sim }{K_{500}}=M_{500}^{2/3+A_m}E{(z)}^{-2/3+A_z}.\end{array}$$(B.1)

We generated 500 random values of A_m and A_z, uniformly distributed in the ranges [ − 0.75, +0.25] and [ − 2.5, +1.0], respectively. At each iteration, we rescaled the entropy profiles using a different realisation of $\stackrel{\sim }{K_{500}}$ and fitted them in the radial range [0.3 − 0.8] R₅₀₀, using the model described in Eq. 8. We then studied how the measured scatter varies as a function of A_m and A_z; our best-fitting residual mass and redshift dependencies are the values of A_m and A_z that minimise the dispersion of the profiles. The results of this procedure are shown in Fig. B1.

Fig. B.1. Constraints to both residual mass (Am) and redshift (Az) dependencies, as defined in Eq. B.1. In the figure, we plot the distribution of the intrinsic scatter of the combined HIGHMz – REXCESS – ESZ sample, as a function of Am (top) and Az (bottom). In the upper panel, measured points are colour coded according Az, while in the lower one according to Am values. Vertical red lines mark the positions of the two minimum values, while black ones are predictions from the self-similar scenario (i.e. Am = Az = 0).

The distribution of the intrinsic scatter σ_int as a function of A_m and A_z describes a well defined surface in the 3D parameter space and the minimum of it is not achieved for values (A_m, A_z) = (0, 0), as it would be the case if entropy profiles scaled self-similarly. The projection of the 3D distribution along the A_m (top) and A_z (bottom) axes allows for a more precise identification of the two values for which the measured scatter is at its minimum. In particular, we find that the pair (A_m, A_z) = (−0.18 ± 0.05, −0.69 ± 0.30) is able to minimise the dispersion of the considered entropy profiles. The dependence of the intrinsic scatter on A_m is stronger, as suggested by the well defined minimum and by the evident curvature, while the one on A_z is less pronounced; this is likely simply due to the redshift range of the considered clusters (z < 0.6), which does not allow more robust constraints on this parameter. The measured value A_m that minimises the dispersion is in agreement with predictions by Ettori et al. (2023), who found A_m = −0.16 ± 0.01; however, they also derived a weaker residual dependence on redshift (A_z = −0.07 ± 0.02) than our measurement, although this is a 2σ effect. Similar results were found by Pratt et al. (2022) regarding the density profiles of ∼120 clusters up to z ∼ 1.1. If we assume that the departure from self-similar predictions comes from density profiles only, their measurements would translate into A_m ∼ −0.15 and A_z ∼ −0.06 for the ICM entropy. Finally, given the adopted colour coding, we notice that, at fixed A_z (top) and A_m (bottom), the distribution of σ_int describes well defined curves in the parameter space; for example, for A_z = 0 we obtain the result of Fig. 12.

Appendix C: ICM entropy from mean density and other mass rescalings

Although we considered the profiles obtained from median densities and scaled by M₅₀₀^Y_X as our reference, a combination of other entropy measurements (e.g. from mean densities) and mass rescalings (e.g. by $M_{500}^{\mathrm{HE}}$ and M₅₀₀^Y_SZ, with 1 − b = 0.8) were used throughout the paper for comparison with other samples. In this section we show the differences of these scaled profiles with respect to the reference ones (Fig. C1), with particular focus on those derived from azimuthal mean densities and scaled by M₅₀₀^Y_X (Fig. C2), which were used in comparison with REXCESS and ESZ samples in Sect. 4.3.

Fig. C.1. Different rescalings presented throughout the paper. In particular, we compare profiles derived with i) mean densities, scaled by M500YX (green triangles); ii) median densities, scaled by $M_{500}^{\mathrm{HE}}$ (red stars); iii) median densities, scaled by M500YSZ (1 − b = 0.8; orange squares) to the reference ones, i.e. obtained with median densities and scaled by M500YX (black dots).

Fig. C.2. Entropy profiles of HIGHMz clusters, derived from azimuthal mean densities. The description to the figure is the same as for for Fig. 4.

In general, the profiles shown in Fig. C2 have similar characteristics to those presented in Sect. 4.1, such as the observed large dispersion in the central regions, which reflects the morphological state of the clusters, and a smaller one towards the outskirts. Also here, the exclusion of temperature bins with S/B < 0.2 from the deprojection allows a regularisation of the profiles at all radii and alleviates the observed flattening at ∼R₅₀₀. Conversely, using mean densities, which have better resolution than median ones, we can explore further the core region of the clusters and derive the ICM entropy (assuming a constant core temperature) down to R ∼ 0.005 R₅₀₀^Y_X. A quantitative comparison between azimuthal median and mean scaled entropy profiles (Fig. C1, green triangles) shows that the latter are ≲5% lower at all radii, a small effect given both to higher densities and masses derived from the Y_X proxy, adopted in the scaling.

To compare with other samples, both observational and simulated ones, we also used two additional different versions of the scaled entropy profiles. In particular, we used profiles derived from median densities and scaled by $M_{500}^{\mathrm{HE}}$ in comparison with X-COP (Sect. 4.3) and M₅₀₀^Y_SZ masses, corrected to mimic a hydrostatic bias of 1 − b = 0.8, in comparison with MACSIS and The300 (Sect. 4.4). The medians of these scaled profiles are shown in Fig. C1 (red stars and orange squares, respectively), compared to the reference one. Using hydrostatic-equilibrium masses to compute K₅₀₀ (Eq. 1) leads to a small increase (≲3%) of the scaled entropy, not significant though; correcting Planck masses and subsequently using them in the rescaling makes the scaled entropy a ∼10% lower with respect to the reference profile.

All Tables

All Figures

	Fig. 1. Distribution of HIGHMz (black dots) clusters in the mass-redshift plane. X-COP (red squares), REXCESS (orange triangles) and ESZ (green stars) clusters are also plotted for comparison. Vertical and horizontal coloured lines mark the median redshift and mass, respectively, of the four samples.
In the text

Fig. 2. Properties of MACSIS (light blue) and The300 (violet) clusters, in comparison to HIGHMz (black). Top panel: distribution of the cluster masses. Bottom panel: distribution in the concentration-centroid shift plane. In both panels, coloured lines show the median value of the corresponding axis for the three samples. MACSIS and The300 clusters are all at redshifts z = 0.34 and z = 0.333, respectively.

In the text

	Fig. 3. Azimuthal median emission measure (left) and projected temperature (right) profiles of clusters in HIGHMz. The profiles are colour coded according to cluster masses (M500YSZ; Table 1). Temperature measurements with S/B < 0.2 are marked in black.
In the text

Fig. 4. Entropy profiles of HIGHMz clusters, derived using azimuthal median densities. The profiles are shown in physical units in the upper panels, while they are scaled using K500YX in the centre and bottom. In the upper and central panels, they are colour coded according to the mass, while in the bottom we show the dependence on the concentration (left) and centroid shift (right), as measured by Campitiello et al. (2022). The two labels ‘Without S/B cut’ and ‘S/B > 0.2’ indicate whether all the temperature bins were used in the deprojection or we limited to measurements with S/B > 0.2, respectively. In the central right panel, we also plot the median profile (solid black) and the measured intrinsic dispersion (grey shaded area). At small radii, dotted lines denote profiles derived assuming a constant core temperature on the density radial grid. Black dash-dotted lines in the central and lower panels are predictions from non-radiative simulations (Voit et al. 2005).

In the text

Fig. 5. Joint fit of the HIGHMz entropy profiles using the constant-plus-power-law parameterisation (Eqs. (4) and (5)). Top panel: best-fitting model (black), super-imposed to observational measurements (grey dots). Black dashed lines mark the intrinsic scatter, while red dash-dotted line is the prediction from non-radiative simulations (Voit et al. 2005). Central panel: ratio between our best-fitting model and predictions from non-radiative simulations. Grey shaded area is the associated statistical error. Bottom panel: local slope of the best-fitting analytical model (black), together with the associated statistical error (grey). In blue are the slopes obtained using fitting with piece-wise power laws (Eq. (6)). Red dash-dotted line is the canonical 1.1 slope.

In the text

Fig. 6. Comparison with other observational samples. On the left, we compare our results to X-COP (red squares), while on the right to REXCESS (yellow triangles) and ESZ (green stars) clusters. Top: median entropy profiles, with their intrinsic scatter. Centre: ratio to the HIGHMz median profile. The grey shaded area is the intrinsic scatter of the HIGHMz clusters, while the horizontal dotted lines indicate a 20% variation. Bottom: radial profiles of intrinsic scatter.

In the text

Fig. 7. Comparison with simulated datasets. From the top to the bottom: (i) measured profiles for MACSIS (light blue), The300 (violet), and HIGHMz, rescaled to mimic a 20% hydrostatic bias, as described in the text (black); (ii) median entropy profiles of the three samples, together with their intrinsic scatters and the red dash-dotted line showing the prediction from Voit et al. (2005); (iii) ratio to HIGHMz, with the grey shaded area showing the intrinsic dispersion of HIGHMz entropies and horizontal dotted lines marking a 20% discrepancy; (iv) radial profile of the intrinsic dispersion for the three samples.

In the text

Fig. 8. Comparison between relaxed and disturbed systems, for HIGHMz (left), The300 (centre), and MACSIS (right). Green is used for the most relaxed clusters, while yellow for most disturbed ones. Black is used for the median profile obtained using all the clusters in each sample. In the bottom panels, we show the ratios to the medians of the full samples, with the dotted lines marking a 20% discrepancy.

In the text

	Fig. 9. Pearson’s correlation coefficient between c (left) or w (right) and scaled entropy K/K500 as a function of cluster radius. Light blue, violet and black are used for MACSIS, The300, and HIGHMz, respectively. The grey shaded area is the statistical error associated with HIGHMz measurements.
In the text

Fig. 10. Relative contribution of potential density and temperature-related systematic uncertainties to the median entropy profile. We show the ratios between the medians $\stackrel{\sim }{K}$ of the scaled entropy profiles, reconstructed using: (i) median and mean densities, without cut for S/B > 0.2 (grey); (ii) mean densities, with and without cut for S/B > 0.2 (purple) and (iii) a combination of the two (red), as detailed in the text. Shaded areas are the 1σ statistical errors.

In the text

Fig. 11. Distribution of central entropies K0 (left) and outer slopes α (right) of the HIGHMz clusters. Plot bars are colour coded according to the median centroid shift w in each bin. Bottom panels are the cumulative distributions of the two histrograms. Blue and golden vertical lines mark the positions of the peaks identified with ACCEPT (Cavagnolo et al. 2009) and REXCESS (Pratt et al. 2010), respectively. The red dotted line on the right is the canonical slope value of 1.1.

In the text

Fig. 12. Constraints to the residual mass dependence Am in the entropy rescaling. Left panel: distribution of the intrinsic scatter of different samples as a function of Am, to identify the global (i.e. within [0.3 − 0.8] R500) Am correction (Sect. 5.3.1). Yellow is used for the combined X-COP – HIGHMz sample, brown for REXCESS – ESZ – HIGHMz, while blue for HIGHMz only. Vertical dashed lines and shaded areas mark the positions of the minima, together with their statistical errors. Right panel: radial dependence of the Am values which minimise the scatter of the combined REXCESS – ESZ – HIGHMz sample (Sect. 5.3.2), together with their statistical errors (black solid line and grey shaded area). Shown in brown (blue) is the global measurements using REXCESS – ESZ – HIGHMz (HIGHMz only) sample. In red are results from REXCESS only (Pratt et al. 2010). Both in the left and right panels, black dashed lines are the prediction from self-similar scenario, i.e. Am = 0.

In the text

Fig. 13. Rescaling using weak-lensing masses. From top to the bottom: (i) entropy profiles of the sub-sample of 16 clusters with weak-lensing masses available, scaled using $M_{500}^{\mathrm{WL}}$, where shaded grey lines show the impact of accounting for statistical errors on $M_{500}^{\mathrm{WL}}$ in the rescaling, as described in the text; (ii) median entropy profiles, obtained using M500YX (red dotted line) and $M_{500}^{\mathrm{WL}}$ (blue dotted line) in the rescaling; continuous lines are the measurements after including statistical errors on mass measurements; (iii) same as centre, but for the intrinsic dispersion; here, simple scalings without including errors on mass measurements are shown as shaded areas.

In the text

	Fig. A.1. Comparison of temperature (left) and density (right) profiles of MACSIS (light blue) and The300 (violet) clusters, with respect to HIGHMz ones (black). Upper, central and lower panels show the same information as in Fig. 7.
In the text

	Fig. A.2. Comparison of mass-weighted (Tmw) and spectroscopic-like (Tsl) temperature profiles, for MACSIS (light blue) and The300 (pink). Plotted are the ratios of the two measurements for each cluster. The two more marked lines are the median ratios for the two samples.
In the text

Fig. B.1. Constraints to both residual mass (Am) and redshift (Az) dependencies, as defined in Eq. B.1. In the figure, we plot the distribution of the intrinsic scatter of the combined HIGHMz – REXCESS – ESZ sample, as a function of Am (top) and Az (bottom). In the upper panel, measured points are colour coded according Az, while in the lower one according to Am values. Vertical red lines mark the positions of the two minimum values, while black ones are predictions from the self-similar scenario (i.e. Am = Az = 0).

In the text

	Fig. C.1. Different rescalings presented throughout the paper. In particular, we compare profiles derived with i) mean densities, scaled by M500YX (green triangles); ii) median densities, scaled by $M_{500}^{\mathrm{HE}}$ (red stars); iii) median densities, scaled by M500YSZ (1 − b = 0.8; orange squares) to the reference ones, i.e. obtained with median densities and scaled by M500YX (black dots).
In the text

	Fig. C.2. Entropy profiles of HIGHMz clusters, derived from azimuthal mean densities. The description to the figure is the same as for for Fig. 4.
In the text