© M. G. Campitiello et al. 2022

1. Introduction

Clusters of galaxies represent a common ground between astrophysics and cosmology. On the one hand, they are the most massive and virialised systems in the Universe and provide a unique opportunity to study processes related to structure formation on both large and small scales, such as cluster and galaxy scales, respectively. On the other hand, the description of their abundance and spatial distribution allows us to derive relevant information on the underlying cosmology, the gravitational processes, and the initial conditions that characterised our Universe. These two areas of research focus on objects that, throughout their lives, experience very different dynamical states, from being very relaxed to very disturbed due, for instance, to a major merger and other astrophysical processes (such as turbulence, feedback, or feeding flows). Relaxed systems are particularly suitable for deriving the cluster total mass, which is, together with the redshift, the most important cluster property used in cosmology. Given the absence of signs of mergers and turbulence, the mass of relaxed systems is derived assuming that both the intracluster medium (ICM) and the galaxies are in hydrostatic equilibrium (HE) within the binding cluster potential (see, e.g., Ettori et al. 2013; Pratt et al. 2019). The X-ray mass estimation turns out to be close to gravitational lensing mass (e.g., Meneghetti et al. 2010; Rasia et al. 2012) and consequently it is considered robust. Disturbed systems do not satisfy the HE assumption, and their mass estimations are characterised by larger uncertainties. However, astrophysical studies mainly focus on these systems, since phenomena such as turbulence are more prominent. For example, since a correlation between the X-ray emission and the presence of giant radio halos has been observed in merging systems, disturbed clusters are considered the perfect laboratories for the study of particle acceleration mechanisms resulting in non-thermal radio emission (e.g., Cassano et al. 2010; Mann & Ebeling 2012).

The identification of the most relaxed and disturbed systems is also essential for the understanding of the absolute scatter that characterises the scaling relations. This topic is one of the most important open issues in the study of clusters, since it is linked to the constraints used for cosmological models (e.g., Lima & Hu 2005). It was found that relaxed and disturbed objects lie in different regions of the scaling relations and that their dynamical state provides the major contribution to the scatter about the relations (e.g., Pratt et al. 2009; Lovisari et al. 2020). For example, Fabian et al. (1994) noticed, for the first time, that the offset of a cluster from the mean relation is linked to the presence of a cool core, which is a typical feature of relaxed systems, while Andrade-Santos et al. (2012) showed that disturbed systems can be used in the scaling relations when the level of sub-structures is known and parametrised so that their positions in the mass-observable planes can be corrected (Ventimiglia et al. 2008).

Finally, the dynamical state is also crucial in the interpretation of survey data because of its impact on the selection function. Merger events could influence the observable used to detect clusters, affecting their identification and selection and increasing (or decreasing) the number of objects observed with respect to what is expected from theoretical mass function. For example, by comparing Planck Sunyaev-Zeldovich effect (SZE, Sunyaev & Zeldovich 1972) clusters with X-ray-selected samples, it was found that the latter ones mainly tend to detect centrally peaked and more relaxed clusters. This behaviour is related to the different dependences of the Sunyaev-Zeldovich (SZ) signal and X-ray emission on the gas density. Since the X-ray emission scales with the square of the gas density, X-ray surveys tend to preferentially detect centrally peaked, more relaxed galaxy clusters (Eckert et al. 2011), which result to be more luminous at a given mass. On the contrary, the SZ signal is less sensitive to the central gas density, and simulations have shown that SZ-surveys are not strongly influenced by the dynamical state of the clusters (Motl et al. 2005). Hence, SZ samples are expected to provide a clean reconstruction of the underlying cluster population. This is a key property for those statistical cluster studies that aim to constrain cosmological models or to probe the physics of structure formation.

In this context, it is clear that the dynamical classification plays an important role when dealing with large samples of galaxy clusters, since it allows either to identify the most suitable set of systems to consider in the analysis or to estimate any systematic effect introduced by the relative fraction of relaxed and disturbed systems. However, obtaining an accurate characterisation of the dynamical state of galaxy clusters is very challenging because multi-wavelength information is required but only available for a few objects. To overcome this limitation it is possible to resort to the analysis of the distribution of the X-ray emission of galaxy clusters: all those processes that can alter the dynamical state of clusters, such as mergers, are indeed expected to leave traces in the ICM distribution, and thus in the X-ray images. Therefore, the identification of a proper method for the characterisation of the X-ray morphology has drawn the attention of the X-ray community over the past thirty years, favouring the development of many procedures. Initially, images were inspected by eye to detect and characterise sub-structures. For example, Jones & Forman (1992) distinguished clusters as follows: ‘single’, ‘double’, ‘primary with small secondary’, ‘complex’, ‘elliptical’ (according to the X-ray contours), ‘off-centre’ (either presenting a difference between the centres in optical and X-ray or showing an X-ray tail extended only in one sector off the X-ray peak), and ‘galaxy’ (when the main contribution to the X-ray emission is provided by the central galaxy). However, a classification based only on the presence of sub-structures may not include all those systems that have not recently interacted with merging massive systems, but despite this show traces of previous interactions in the form of either a strong elliptical shape (Buote & Tsai 1996; Pinkney et al. 1996; Plionis 2002) or a variation of their X-ray centroid (Mohr et al. 1995). Furthermore, the visual classification turns out to be both subjective, since different researchers may provide different classifications for individual clusters, and very time-consuming in the case of large samples of clusters such as those expected from future surveys. It then became necessary to define more robust indicators able to objectively quantify even small deviations from a perfectly regular and spherically symmetric emission. Among the commonly used parameters, we point out in particular the axial ratios (Mohr et al. 1993), the centroid shift (Mohr et al. 1995; O’Hara et al. 2006), the power ratios (Buote & Tsai 1995, 1996), the light concentration (Santos et al. 2008), and two parameters arising from morphological analysis of galaxies, that are, the asymmetry and the smoothness (Lotz et al. 2004). To make the classification more effective and less affected by bias and projection effects, these and a combination of other parameters were also used (e.g., Rasia et al. 2013; Mantz et al. 2015; Andrade-Santos et al. 2017; Lovisari et al. 2017; McDonald et al. 2017; Nurgaliev et al. 2017; Ghirardini et al. 2021).

In recent years, this approach has been used with the morphological analysis of simulated samples (Rasia et al. 2013; Barnes et al. 2018). In the case of hydrodynamical simulations, the cluster state is typically parametrised though a few indicators that are evaluated in 3D. These consider and quantify various conditions of a dynamically active object, such as the presence of well-identified sub-structures (indicating recent mergers), the displacement between the barycentre and the total density peak or the minimum of the potential well (indicating either past merger activities or a minor merger with a small impact parameter), and the ratio between the kinetic or thermal energy and the gravitational energy. These parameters allow us to obtain prior knowledge of the dynamical state of the considered systems (e.g., Cialone et al. 2018; Cui et al. 2018; Capalbo et al. 2021; De Luca et al. 2021). For this reason, simulations can be used to calibrate and check the robustness of the morphological parameters (e.g., Weißmann et al. 2013). From the comparison between observations and simulations it was found that, when used on single clusters, the morphological parameters should be used with extreme caution because they could be affected by substantial uncertainties (Böhringer et al. 2010). Jeltema et al. (2008), for example, found that less than one half of clusters that are classified as relaxed by the power ratios are truly relaxed, and 4%–10% of these relaxed clusters have very disturbed morphologies when viewed from other orientations. On the other hand, not only are observations limited by projection, resolution, and background, but also simulations are still unable to capture the complex multi-scale physics of the core, especially with regard to the AGN physics implementation (see Gaspari et al. 2020 for a review).

In this work, we performed a morphological analysis of the X-ray images obtained from XMM-Newton for the 118 clusters that constitute the CHEX-MATE (Cluster HEritage project with XMM-Newton – Mass Assembly and Thermodynamics at the Endpoint of structure formation) sample. The aim of the CHEX-MATE project is to set the stage for future X-rays missions by providing both an overview of the statistical properties of the underlying cluster population and an improvement of the analysis techniques developed to date. The morphological analysis itself reflects these two goals. The complete and homogeneous X-ray exposures of the CHEX-MATE objects allow us to derive, for the first time, a uniform characterisation of the X-ray morphology of the cluster’s underlying population, providing the entire community with an overview of the dynamical state of the CHEX-MATE clusters, which will be useful for the identification of the proper set of systems to use in specific analyses. Furthermore, our analysis aims to check the techniques developed up to now by testing the efficiency of the morphological parameters.

The paper is structured as follows. In Sect. 2, we present the sample and the dataset that we analysed. In Sect. 3, we describe the results of a preliminary visual classification of the CHEX-MATE sample, which we used as a reference for the interpretation of the results. In Sect. 4, we describe the set of morphological parameters that we estimated from the X-ray images. In Sect. 5, we report the results of the morphological analysis carried out on the X-ray observations. In Sect. 6, we investigate the robustness of the morphological parameters by assessing the systematic uncertainties that affect their measurements. In Sect. 7, we combine them to obtain a unique indicator of the grade of relaxation of clusters. In Sect. 8, we discuss our results, and we draw our conclusions in Sect. 9. Throughout the paper, if not otherwise stated, we assume a flat ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7 and H₀ = 70 km s⁻¹ Mpc⁻¹. The variables M_Δ and R_Δ are the total mass and radius corresponding to a total density contrast Δρ_c(z), where ρ_c(z) is the critical density of the Universe at the cluster redshift (for example, M₅₀₀ = (4π/3)500${\rho }_{\text{c}}(z)R_{500}^3$).

2. Dataset

2.1. The CHEX-MATE sample

The CHEX-MATE programme¹ is described in detail in CHEX-MATE Collaboration (2021). It is a 3 mega-second Multi-Year XMM-Newton Heritage Programme to obtain X-ray observations of a minimally biased, signal-to-noise-limited sample of 118 galaxy clusters detected by Planck through the SZE. The project has been developed to provide an accurate vision of the statistical properties of the underlying population of clusters, to measure how the gas properties are shaped by collapse into the dark matter halo, to uncover the origin of non-gravitational heating, and to resolve the major uncertainties in mass determination that limit the use of clusters for cosmological parameter estimation. To achieve these aims, a sample of 118 Planck clusters (Planck Collaboration VIII 2011; Planck Collaboration XXVII 2016; Planck Collaboration XXIX 2014) populating two different sub-samples was selected through their SZE signal (S/N > 6.5) according to the following criteria: Tier 1, consisting of 61 objects located at low redshift in the northern sky (0.05 < z < 0.2 and Dec > 0, with 2 × 10¹⁴ M_⊙< M₅₀₀ < 9 × 10¹⁴ M_⊙) and providing an unbiased view of the population at the most recent time; Tier 2, including the most massive systems to have formed thus far in the history of the Universe (z < 0.6 with M₅₀₀ > 7.25 × 10¹⁴ M_⊙). These two sub-samples have four clusters in common.

The XMM-Newton observations are characterised by an exposure time that ensures a S/N = 150 within R₅₀₀ in the [0.3–2.0] keV band. This condition has been requested to estimate the temperature profile at least up to R₅₀₀ (with a precision of ±15% in the region [0.8–1.2]R₅₀₀) and to obtain a measurement of both the mass derived from the Y_X mass proxy (Kravtsov et al. 2006), where Y_X = M_g, 500T_X, M_g, 500 is the mass of gas within R₅₀₀ and T_X is the spectroscopic temperature estimated in the [0.15–0.75]R₅₀₀ range) with ±2% of uncertainty and the mass derived from HE at R₅₀₀ with an ∼15–20% precision level. For more details on the sample, on the scientific goals, and on the strategy used to observe these systems in X-ray homogeneously and to follow them up in other wave-bands, we refer the reader to CHEX-MATE Collaboration (2021).

2.2. Preparation of the X-ray images

Images were produced using the pipeline developed during the XMM-Newton Cluster Outskirts Project (X-COP, Eckert et al. 2017; Ghirardini et al. 2019) and adopted by the CHEX-MATE collaboration. In particular, the XMM-Newton data were processed using the SAS software (version 16.1.0) and the extended source analysis software (ESAS) package (Snowden et al. 2008). Count images, exposure maps, and particle background maps are extracted in the narrow [0.7–1.2] keV band, where the ratio between the source and background emission is maximised and, consequently, the systematics related to the subtraction of the EPIC background are minimised (Ettori et al. 2010). A detailed description of the procedure adopted will be presented in Bartalucci et al. (in prep.), and a complete gallery of the images is shown in Fig. 6 in CHEX-MATE Collaboration (2021). According to the CHEX-MATE pipeline, point sources in our observations are identified using the SAS tool ewavelet in two bands (0.5 − 2 keV and 2 − 7 keV). Furthermore, we applied a filter to the LogN − LogS distribution as described in Ghirardini et al. (2019) to ensure a uniform level of the cosmic X-ray background (CXB) emission across the field of view. For the estimation of morphological parameters, we further inspected images by eye in order to identify residual point sources, which could affect our measurement. We masked identified point sources and filled the holes using an interpolation with the surrounding pixels. We point out that only point sources were masked, and thus sub-structures related to major or minor mergers are still present in the images, allowing a correct identification of the cluster’s dynamical state. The pixel size is 2.5 arcsec.

3. Initial visual classification

To study the link between the morphological parameters and the dynamical state of clusters, we first realised a visual classification of the objects of the CHEX-MATE sample. In particular, a group of seven X-ray astronomers inspected the images and rated the relaxation state of the clusters with a grade that ranges from 0 (most relaxed, i.e. circular X-ray isophotes and without sub-structures) to 2 (most disturbed, i.e. double or complex objects with clear sign of merging). Results were then averaged and objects with rounded values equal to 0 were classified as relaxed (R, 19 clusters), objects with rounded values equal to 2 as disturbed (D, 37 clusters), and objects with rounded values equal to 1 as mixed (M, 62). Examples of a relaxed, mixed, and disturbed system, as identified by the visual classification, are shown in Fig. 1, while in Fig. 2 we report the final dynamical state as a function of the redshift and the mass. It is possible to notice that the majority of the relaxed objects are located at redshift lower than z < 0.25 and have a mass higher than ∼4 × 10¹⁴ M_⊙. No particular trend is observed for the disturbed and mixed class. The dynamical state obtained from this first analysis and the related uncertainty are reported in Table A.1, ‘visual’ column.

Fig. 1. Example of relaxed (left), mixed (centre), and disturbed (right) systems identified by the visual classification. The horizontal and vertical dimensions of the boxes are equal to 2R500.

Fig. 2. Distribution in the M500 − z plane of the 118 CHEX-MATE clusters. The colour scale represents the dynamical state obtained from the visual classification. The shaded area indicate the Tier 1 and Tier 2 redshift ranges in grey and violet, respectively.

4. Morphological parameters

In this section, we introduce the methods for the sub-structure and morphology characterisation of the CHEX-MATE objects. To estimate the morphological parameters listed below, we considered a circular region within R₅₀₀ centred on the cluster X-ray peak. This choice was adopted to avoid the contamination of signatures related to accretion processes, which are expected at larger radii (e.g., Roncarelli et al. 2006). Hydrodynamical simulations (e.g., De Luca et al. 2021) and observations (Ghirardini et al. 2019; Eckert 2014) indeed show that within R₅₀₀ clusters are relatively relaxed, unless a merger event modifies the existing conditions. Furthermore, the CHEX-MATE observation strategy provides coverage of this area for all the clusters of the sample (CHEX-MATE Collaboration 2021). For these reasons, we considered R₅₀₀ the optimal assumption to obtain a complete view of the dynamical state of clusters. The clusters’ centres are set at the brightest pixel of a Gaussian-smoothed (σ ∼ 15 pixels), background-subtracted, exposure-corrected surface-brightness image. The value of R₅₀₀, was derived from the Planck PSZ2 masses (Planck Collaboration XXVII 2016). In our analysis, seven morphological parameters were taken into consideration: the light concentration, c, the centroid shift, w, two power ratios, the asymmetry, A, the smoothness, S, and the ellipticity, η. It results that A and S are not robust parameters, since they are strongly influenced by the signal-to-noise ratio (see Appendix B for the discussion). For this reason, we decided to cover their analysis in the appendix, excluding them from the main text. The analysis of ellipticity is also discussed in Appendix B, since its behaviour reproduces that of the quadrupole power ratio P₂ (a measure of the ellipticity of clusters; see the following sub-sections for more details).

4.1. Light concentration

The light concentration parameter (hereinafter, concentration), c, is defined as the ratio of the surface brightness (SB) inside two concentric apertures. It was introduced by Santos et al. (2008) to identify cool-core clusters at high redshift, using two apertures at 40 and 400 kpc, chosen to maximise the separation of concentration values between cool-core (CC) and non cool-core (NCC) clusters in the sample analysed there. However, the cluster volume enclosed in a fixed aperture evolves significantly with redshift, and this behaviour can affect the selection of relaxed and disturbed systems (Hallman & Jeltema 2011). Since our sample spans a large redshift range, we chose to define the concentration parameter as a function of the overdensity radii using the following apertures:

$$\begin{array}{c}\hfill c=\frac{\mathrm{NC}(r<0.15R_{500})}{\mathrm{NC}(r<R_{500})},\end{array}$$(1)

where NC is the number of counts corrected for the exposure map measured in the considered aperture. On the basis of the above definition, the concentration computed from images is not corrected for the point spread function (PSF). The effects of this choice are that more distant objects are characterised by systematically lower concentrations than the low-redshift clusters, because more photons located in the centre are spread out across larger regions. The implication of our choice is discussed detail in Sect. 8.4.3.

4.2. Centroid shift

The centroid shift parameter (Poole et al. 2006; Maughan et al. 2008), w, is defined as the standard deviation of the projected separation between the X-ray peak and the centroid of the X-ray surface brightness computed within N (=10 in our case) apertures of increasing radius:

$$\begin{array}{c}\hfill w=\frac{1}{R_{500}}[\frac{1}{N-1}{\displaystyle \sum _i{({\mathrm{\Delta }}_i-\overline{\mathrm{\Delta }})}^2{]}^{\frac{1}{2}},}\end{array}$$(2)

where Δ_i is the distance between the X-ray peak and the centroid of the ith aperture, and R₅₀₀ is the radius of the largest aperture. This parameter is useful for the characterisation of the dynamical state of clusters because it is sensitive to the presence of X-ray bright clumps and sub-structures, which can produce significant changes on the X-ray centroid.

4.3. Power ratios

The power ratio parameters were first introduced by Buote & Tsai (1995) and are based on the idea that the X-ray surface brightness of a cluster could be the representation of its projected mass distribution. They are computed as a multipole decomposition of the two-dimensional projected mass distribution inside a certain aperture; but, instead of the mass, the X-ray surface brightness is used. The m-order (m > 0) power ratio is defined as P_m/P₀, with:

$$\begin{array}{c}\hfill P_0={[a_0\ln (R_{500})]}^2,\end{array}$$(3)

and

$$\begin{array}{c}\hfill P_m=\frac{1}{2m^2{R}_{500}^{2m}}(a_m^2+b_m^2),\end{array}$$(4)

where a₀ is the total intensity within the aperture radius R₅₀₀ and the moments a_m and b_m are calculated by:

$$\begin{array}{c}\hfill a_m(R)={\displaystyle {\int }_{R<R_{500}}S(x)R^{m}cos(m\varphi ){\mathrm{d}}^{2}x}\end{array}$$(5)

and

$$\begin{array}{c}\hfill b_m(R)={\displaystyle {\int }_{R<R_{500}}S(x)R^{m}sin(m\varphi ){\mathrm{d}}^{2}x,}\end{array}$$(6)

where S(x) is the surface brightness at the position x = (R, ϕ). The quadrupole power P₂ quantifies the ellipticity of the clusters, P₃ informs us about bimodal distribution and is the most useful to identify asymmetries or presence of sub-structures, and P₄ is similar to P₃ but more sensitive to smaller scales (for this reason they are strongly correlated). In the following analysis, we focus on the ratios P₂/P₀ and P₃/P₀ (hereafter P₂₀ and P₃₀, respectively). This choice is due to the properties of P₃₀, which is one of the most unambiguous indicators of an asymmetric cluster structure (Jeltema et al. 2005, 2008; Cassano et al. 2010; Weißmann et al. 2013; Rasia et al. 2013; Lovisari et al. 2017; Cialone et al. 2018), and to our aim to investigate the power of an indicator of the ellipticity of clusters, such as P₂₀.

5. Morphological analysis of the CHEX-MATE sample

Starting from the dataset and the methods introduced above, we carried out a morphological analysis of the CHEX-MATE sample. For each cluster, we estimated the values of the four parameters presented in Sect. 4 and their associated statistical errors through a Monte Carlo 100 re-sampling of the counts per pixel of the original image according to their Poissonian error the technique is implemented in an IDL routine already used for the analyses presented in Cassano et al. (2010), Donahue et al. (2016), and Lovisari et al. (2017).

We then investigated the presence of correlations between pairs of these parameters by determining the Spearman rank correlation coefficient, r, a value that varies in the −1 ≤ r ≤ 1 range, assuming extreme values (−1 or 1) when each of the variables tested is a perfect monotone function of the other. The results of the analysis are reported in Fig. 3, left panel. All the parameter pairs are characterised by significant correlations (i.e. |r|> 0.5), with P₂₀ − P₃₀ and c − w showing the highest Spearman coefficient (r = 0.71 and r = −0.7, respectively). The corner plot also shows the parameter distributions of the relaxed and disturbed systems, as defined by the visual classification (see Sect. 3). We found that clusters with different dynamical states are located in distinct regions of the parameter-parameter planes; for example, in the c − w plot, relaxed clusters occupy the lower right region, while disturbed systems are placed in the upper left region. It is thus possible to identify thresholds above (or below) which clusters can be classified as relaxed (or disturbed) and vice versa, depending on the behaviour of the considered morphological parameter. As highlighted by the density plot of Fig. 3 (left panel), the concentration and the centroid shift are confirmed to be powerful identifiers of the relaxed and disturbed populations (Santos et al. 2008; Cassano et al. 2010; Hudson et al. 2010; Lovisari et al. 2017); the distribution of the values of these two classes of objects are indeed well separated. Also, the power ratios are able to distinguish between relaxed and disturbed systems. However, it is possible to observe an overlap of the distributions of the values of these parameters, thus suggesting that by defining a threshold for the identification of the dynamical state, some clusters may be erroneously classified as relaxed instead of disturbed and vice versa. Finally, for all four parameters, the systems classified as ’mixed’ represent an intermediate class between the relaxed and disturbed populations.

Fig. 3. Distribution of the morphological parameters for the CHEX-MATE (left) and simulated (right) samples. The r values reported in the boxes are the Spearman coefficients of the considered parameter pairs. Left: the different colours represent the dynamical state obtained from the visual classification (see Sect. 3). Right: the different colours represent the dynamical state obtained using the dynamical indicator χ (see Sect. 8.4.1). The background colours of the plots represent the level of agreement with the observations: green is for 1σ, yellow is for 1–3σ, and red is for a level of agreement over 3σ.

6. Systematics in the measurements of morphological parameters

In this section, we discuss how we tested the robustness of the morphological parameters by investigating whether their behaviour is influenced by the quality of the images used or by the assumptions or criteria adopted in our analysis. In particular, we investigated how the exposure time and the region considered for the estimation of the morphological parameters can alter our analysis.

6.1. Effects of the exposure time

The CHEX-MATE observations are characterised by their high quality. To understand the robustness of our analysis, we tested the stability of the morphological parameters when images with lower exposure times are used. To this aim, we considered a sub-sample of 20 CHEX-MATE clusters (5 Relaxed, 5 Disturbed, and 10 Mixed as defined by the visual classification) and we repeated the morphological analysis of the X-ray images now selected with two different exposure times, one of which halved the total exposure time, t_exp = 0.5 × t_exp, total, and the other with a minimal value of 5 ks. We then computed the dispersion Δ following this relation:

$$\begin{array}{c}\hfill \mathrm{\Delta }=|{log}_{10}({\mathcal{P}}_{\mathrm{red}})-{log}_{10}({\mathcal{P}}_{\mathrm{or}})|,\end{array}$$(7)

where 𝒫_red are the parameters computed using images with reduced exposure times and 𝒫_or are the parameters computed using the original images. The median values and the interquartile ranges (IQR) of Δ obtained from the two types of images are reported in Table 1 (first and second rows). In Fig. 4 (first row), we also plot the values of the morphological parameters computed from the new images (violet, 50% t_exp, total, and grey t_exp = 5 ks) as a function of the values computed using the original ones. It is possible to observe a low scatter between the two types of estimations.

Fig. 4. Analysis of the robustness of the morphological parameters. First row: comparison between the parameters estimated from the original images (x-axis) and the parameters estimated from images with halved texp (violet) and with texp = 5 ks (grey). Second row: comparison between the parameters estimated using R500 (x-axis) and the parameters computed in a region of r = 0.95R500 (violet) and r = 1.05R500 (grey).

6.2. Effects of the assumed R₅₀₀

As presented in Sect. 4, all the parameters are estimated inside a region of radius r = R₅₀₀ derived from the Planck PSZ2 masses. However, uncertainties that may affect the measure of the Planck masses could also influence the estimation of R₅₀₀. The understanding of how this could impact our analysis is crucial. To investigate this point, we considered the CHEX-MATE clusters and computed the morphological parameters inside two circular regions of radius r ∼ 1.05 ⋅ R₅₀₀ and 0.95 ⋅ R₅₀₀. Also in this case, we computed the dispersion Δ:

$$\begin{array}{c}\hfill \mathrm{\Delta }=|{log}_{10}({\mathcal{P}}_{R_{500}\pm 5\%})-{log}_{10}({\mathcal{P}}_{R_{500}})|,\end{array}$$(8)

where 𝒫_{R500 ± 5%} are the parameters computed using images with increased or decreased radius and 𝒫_R500 are the parameters computed using the original aperture. In Fig. 4 (second row), we show the scatter plot of the new parameters and the old ones, while in Table 1 (third and fourth rows) we report the medians and the interquartile ranges of Δ. We found that no significant difference is present when taking into account possible uncertainties related to the radius estimations.

In addition to this test, we also compared our estimations of the morphological parameters with the ones obtained using a radius region equal to 0.5 ⋅ R₅₀₀. The change of the radius in which to compute the morphological parameters could indeed affect the classification of the dynamical state. The results of this analysis are reported in Appendix C.

7. Construction of the M parameter

The test presented in the previous section highlights that the centroid shift, the power ratios, and the concentration could be considered powerful and robust morphological indicators. Therefore, we decided to combine the information included in these parameters in order to build a unique indicator of the grade of relaxation of a cluster. This new indicator allows us to establish a ranking of the dynamical state of the clusters of a sample, which can be used to to identify the population of the most relaxed and most disturbed objects. The definition of this new quantity is (see e.g., Rasia et al. 2013; Cialone et al. 2018):

$$\begin{array}{c}\hfill M={\displaystyle \sum }\frac{{log}_{10}({\mathcal{P}}^{{\alpha }_{\mathcal{P}}})-⟨{log}_{10}({\mathcal{P}}^{{\alpha }_{\mathcal{P}}})⟩}{{\sigma }_{{log}_{10}({\mathcal{P}}^{{\alpha }_{\mathcal{P}}})}},\end{array}$$(9)

where σ is the standard deviation of the considered parameter, 𝒫, and the term α_𝒫 is considered equal to −1 only in the case of the concentration, otherwise it is fixed as equal to +1. For each cluster, M represents the sum of the differences of the four parameters for the mean of their distributions, normalised by their standard deviations. The log-scale was introduced to take into account the different ranges of values covered by the four morphological parameters. According to this definition, relaxed systems are expected to show low values of M, while disturbed systems should be characterised by high values. Furthermore, it appears clear that the M parameter does not provide an ‘absolute’ grade of relaxation, but just a relative value based on the analysed sample: it ranks the clusters from the most relaxed to the most disturbed, given the distribution of the parameters of the sample.

We estimated M for the objects of the CHEX-MATE sample, and we identified the most relaxed and disturbed systems of the sample using as reference the visual classification presented in Sect. 3. In particular, using the fractions of relaxed and disturbed systems obtained by that classification, we verified whether the first 19 (or 37) objects with the lowest (or highest) values of M are effectively classified as relaxed (or as disturbed) by the visual classification. If that is the case, we refer to correct detection, C, while if the objects are classified as relaxed (or disturbed) by M and as disturbed (or relaxed) by the visual classification we refer to wrong detection, W. We found M to be highly efficient in identifying the most relaxed and disturbed systems, with correct detection equal to C_R = 79% and C_D = 70% and no wrong detection (W_R and W_D equal to zero). This means that the incorrect detections are mostly due to the mixed population. A comparison between the M and visual classifications is shown in Fig. 5, left panel. In Appendix D, we report the results obtained from the principal component analysis (PCA), which is another method used to combine the information included in the four morphological parameters.

Fig. 5. Distribution of the M parameter. Left: comparison between the CHEX-MATE dynamical state obtained from the visual classification and from the M parameter. The colours represent the dynamical classification obtained on the basis of the M parameter (blue for relaxed, red for disturbed, and green for mixed clusters). Right: distribution of the simulated sample in the χ − M space. The colours represent the relaxed (blue), disturbed (red), and mixed (green) class defined on the basis of χ (see Sect. 8.4.1).

8. Discussion

8.1. Distribution of the morphological parameters

In this sub-section, we describe our investigation into the presence or absence of bimodality in our parameter distributions. In particular, we performed a maximum likelihood fit of our unbinned data using the normal mixture model and a commonly used positively skewed function: the Weibull distribution.

For this part of the analysis, we used the MCLUST (Fraley & Raftery 2002; Scrucca et al. 2016) package and the FITDISTR and FITDIST function of the MASS (Venables & Ripley 2013) and fitdistrplus (Delignette-Muller & Dutang 2015) packages in the software environment R, version 3.6.3 (R Core Team 2021). MCLUST allowed us to perform a cluster analysis (or ’clustering’), which consists of grouping together a set (or cluster) of similar objects in a dataset. The objects of each cluster are comparatively more similar to objects of that group than those of the other clusters. Clustering is a main task of exploratory data analysis (EDA) and provides us the opportunity to see what the data can tell us beyond the formal modelling or the a priori hypothesis. The MCLUST method realises a maximum likelihood fit assuming that one to nine normal components are present in the data. Also the function FITDISTR performs a maximum likelihood fit of the data using some probability distribution functions, either calculated using analytic formulae (as, e.g., in the log-normal case) or computed by optimising the likelihood. This method was used to fit our data with the Weibull distribution.

To compare and select the most appropriate models for the description of the empirical dataset, we considered the Bayesian information criterion (BIC, Schwarz 1978), defined as BIC = 2ln(L)−klog(n), where L is the likelihood, k is the number of parameters of the model, and n is the number of data points; klog(n) is the penalty term that compensates the difference in likelihood due to an increase in the number of fitting parameters. The best model is the one that maximises the BIC. For the interpretation of the differences between the BIC values obtained with different models, we considered the following commonly adopted thresholds: a BIC difference of 0–2 is a weak confirmation, 2–6 a positive confirmation, 6–10 a strong confirmation, and > 10 a very strong confirmation of the model with the greater BIC value (Kass & Raftery 1995; Raftery 1995).

Since the Weibull function is not defined for negative values, we applied a normalisation to the morphological parameters. In particular, we multiplied c, w, P_20, and P_30, respectively, by 10³, 10⁴, 10⁹, and 10¹⁰, and we added 10 to the M values. This way, we were able to consider their log distributions. The results of our analysis are shown in Fig. 6 and Table 2. For all five parameters, the clustering analysis revealed the presence of a single component. For this reason, the BIC values reported in Table 2 are related to a single Gaussian component model. We concluded that there are no parameters showing signs of bimodality and, therefore, that the cluster population cannot be easily divided in two populations. A similar result was already observed for c. P₂₀ and P₃₀ in a recent analysis on the eFEDS sample (Ghirardini et al. 2021). Given the BIC values obtained, it is not possible to unambiguously identify the model that best fits our data. In particular, for w, P_20, and P_30, the discrepancy between the BIC values is ∼2 and represents only a weak confirmation for the model with the highest BIC value (i.e. Weibull for w, P_30, and M, and Gaussian for P₂₀). Concerning the concentration, the discrepancy between the BIC values is 5 and could be considered a positive confirmation for the single Gaussian component model. We thus conclude that the distribution of the concentration is log-normal. Similar behaviour was already found both in the Planck-selected sample analysed in Rossetti et al. (2017), where a log-normal best fit for c contrasted the bimodal best-fit distribution of the X-ray-selected sample ME-MACS and the eFEDS sample analysed in Ghirardini et al. (2021). In this latter analysis, the distributions of P₂₀ and P₃₀ were also taken into consideration. As in our case, it was not possible to unambiguously identify the best-fit model for these parameters, but skewed and log-normal distributions were preferred.

Fig. 6. Distribution of c, w, P20, P30, and M. The solid line represents the model obtained from the clustering analysis, while the dashed one represents the Weibull function. The parameters c, w, P20, and P30 were multiplied by 103, 104, 109, and 1010, respectively. For the M parameter, we added 10 to the original values instead.

Analogous behaviour was observed by repeating this analysis on Tier 1 and Tier 2, separately. The results are shown in Table 2. Also in this case, the BIC values reported for the Gaussian mixture models are related to a single Gaussian component function, which proves to be the best one in reproducing the distributions of the five parameters. Also, in this case we found no sign of bimodality either in the uniform low-z sub-sample or in the uniform high-mass sub-sample. However, since the BIC values are very similar, it is not possible to derive which is the best-fit model for our data.

Moreover, we investigated the distributions of the morphological parameters with respect to two intrinsic properties of clusters: the mass and the redshift. As one can see from Fig. 7, no specific trend is observed, and the low values of the Spearman coefficients highlight the absence of correlations. This result is in agreement with previous findings by Böhringer et al. (2010), Mantz et al. (2015), Lovisari et al. (2017), and Rossetti et al. (2017).

Fig. 7. Mass and redshift distributions of the morphological parameters. The different colours represent the Tier 1 (grey) and Tier 2 (violet) objects. In the boxes is reported the Spearman coefficient, r.

8.2. Final CHEX-MATE classification

In Fig. 5 (left panel), we present a comparison between the dynamical state derived from the visual classification (y-axis) and the one derived on the basis of the M values (colour scale). As already highlighted by the quantities C_R and C_D in Sect. 7, some clusters may be classified as mixed instead of relaxed or disturbed by one of the two classifications. In order to investigate the differences between the two classifications, we inspected the objects for which an agreement was not obtained by eye. For example, we verified that some objects classified as relaxed by M show not only a centrally peaked emission but also some sub-structures in their outskirts. Due to this feature, these systems were classified as mixed by the visual inspection. In order to obtain a classification of the sample that is as accurate as possible, we decided to only define the clusters for which the two classifications provide the same results as relaxed (or disturbed). Therefore, our analysis identified 15 relaxed and 25 disturbed clusters. To have an overview of the dynamical state of the entire sample, we realised a continuous classification using the following criteria. The first 15 and the last 25 objects of the continuous classification are, respectively, the most relaxed and disturbed systems identified by both the visual and M classification. All the other systems are classified as mixed. Inside these three populations, clusters are ranked on the basis of M (from the lowest to the highest values). The final rank is presented in Table A.1.

8.3. Comparison with other observed samples

The characterisation of the dynamical state of cluster samples has been extensively explored. However, it is not easy to compare the results obtained from different analyses. The definition of the morphological parameters is often different and depends on the goals or the limitations of the analysis. In addition to this, we also have to consider that selection effects may affect the definition of cluster samples. In particular, many studies have highlighted the presence of discrepancies between SZ and X-ray samples; the latter typically show a higher fraction of cool core systems (i.e. systems with a centrally peaked emission usually defined as relaxed), while the SZ-selected clusters are characterised by a higher fraction of sub-structures than the X-ray selected systems (Rossetti et al. 2017). However, not all the analyses converge towards this result. For example, using the photon asymmetry and the centroid shift parameters, Nurgaliev et al. (2017) found no significant statistical difference between the X-ray morphology of X-ray- and SZ-selected samples, suggesting that the two are probing similar populations of clusters. A comparison between the X-ray morphology of the CHEX-MATE clusters and other samples analysed in literature may help us understand how the CHEX-MATE sample is able to represent the underlying cluster population.

To this aim, we compared CHEX-MATE with three X-ray-selected samples and one SZ-selected sample. The first of these is an extended version of the MACS sample described in Mann & Ebeling (2012; ME-MACS hereafter). As MACS, this sample is built from the RASS Bright Source Catalogue by Voges et al. (1999), with a flux limit of f_RASS[0.1–2.4] KeV > 1 × 10⁻¹² erg cm⁻² s⁻¹. In contrast to MACS, which is composed of the most distant systems (z > 0.3), ME-MACS also includes lower redshift clusters (z > 0.15) and has an additional luminosity cut of L_RASS[0.1–2.4 keV] > 5 × 10⁴⁴ erg s⁻¹. These features make the ME-MACS sample a purely X-ray-selected sample, based on a flux-limited survey whose distribution is similar to the CHEX-MATE redshift distribution. Secondly, we have the The REpresentative XMM-Newton Cluster Structure Survey sample (REXCESS, Böhringer et al. 2007), which is a representative and statistically unbiased sub-sample of 33 galaxy clusters extracted from the REFLEX cluster catalogue with a rigorous selection in the luminosity-redshift space (see details in Böhringer et al. 2007). Thirdly is the high-z part of the ROSAT PSPC 400 deg² cluster survey (Burenin et al. 2007), abbreviated hereafter as 400d, for our X-ray-selected sample. This sample consists of 36 clusters in the redshift range of 0.35 < z < 0.9 and the mass range of 10¹⁴ M_⊙ < M₅₀₀ < 5  ×  10¹⁴ M_⊙. Finally, we have the 2500 deg² SPT survey of Bleem et al. (2015). This sample is composed of 90 clusters which are among the most massive of the SPT-selected clusters. The systems span the redshift range of 0.25 < z < 1.2 and the mass range range 2  ×  10¹⁴ M_⊙ < M₅₀₀ < 2  ×  10¹⁵M_⊙.

The parameters used for this comparison are the concentration estimated by Rossetti et al. (2017) for the ME-MACS sample, the centroid shift obtained by Böhringer et al. (2010) for the REXCESS sample, and the centroid shift obtained by Nurgaliev et al. (2017) for the 400d and SPT samples. Using the procedure adopted in our analysis, we recomputed c and w following the definition used in those studies. In particular, for the concentration we adopted two radius apertures equal to 40 and 400 kpc. For the centroid shift of Böhringer et al. (2010), we used a region included from [0.1–1]R₅₀₀. For the centroid shifts of Nurgaliev et al. (2017), we used the definition of Böhringer et al. (2010), and we scaled the values by a factor of 1.5. This normalisation is required because the definition of w presented in Böhringer et al. (2010) (which is the same adopted in this paper) differs from the one of Nurgaliev et al. (2017). According to this latter, w is the squared difference between the position of the centroid and the average position of the centroid. The value of the normalisation (1.5) was already determined in Nurgaliev et al. (2017). Results are shown in Fig. 8. We found that CHEX-MATE has a higher fraction of objects with low concentration compared to ME-MACS and a lower fraction of objects with lower values of the centroid shift. These results are confirmed by the KS test, which show a p-value lower than 0.01 in both cases. Concerning 400d, we found that, despite it being an X-ray-selected sample, it shows behaviour similar to that of the SZ-selected samples (CHEX-MATE and SPT). The differences observed between this sample and other X-ray-selected samples have already been debated in the literature (e.g., Santos et al. 2010; Ma et al. 2011; Mantz et al. 2015; Rossetti et al. 2017) and are related to the detection procedures adopted, which made 400d rather unique among X-ray samples.

Fig. 8. Left: comparison between the concentration values of the CHEX-MATE and ME-MACS sample. Right: comparison between the centroid shift values of the CHEX-MATE, REXCESS, 400d and SPT samples.

8.4. Comparison with simulations

In this sub-section, we present the result of the comparison between the CHEX-MATE observations and the simulations. The simulated sample at our disposal is provided by THE THREE HUNDRED² collaboration (Cui et al. 2018) and is composed of 1564 objects spanning a wide range of redshift (0 < z < 0.59) and masses (M₅₀₀ > 1.1 × 10¹⁴ h⁻¹ M_⊙). For each object, three images related to three different orientations are provided, thus allowing us to build a sample of 4692 (1564 × 3) different maps.

8.4.1. Dynamical state of the simulations

The additional information provided by the analysis of smoothed-particle hydrodynamics (SPH) simulations resides in the knowledge of the physical properties of each particle. In order to determine the dynamical state of simulated clusters, it is thus possible to take advantage of this by referring to quantities computed in 3D that would be unreachable with the current observational techniques. In this analysis, we characterised the grade of relaxation of the simulated sample by means of the following indicators. The first one is the mass fraction of all sub-halos in the cluster, f_s, where the sub-halos are identified with the Amiga Halo Finder, AHF3 (Knollmann & Knebe 2009), whenever the structure has at least 20 particles. This parameter is defined as

$$\begin{array}{c}\hfill f_{\mathrm{s}}=\frac{{\mathrm{\Sigma }}_i{M}_i}{M_{500}},\end{array}$$(10)

where M₅₀₀ is the mass of the cluster enclosed in R < R₅₀₀ and M_i is the mass of the sub-halos in the same volume. The second one is the offset of the centre of mass, Δ_r, defined as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_r=\frac{|r_{\mathrm{cm}}-r_{\mathrm{c}}|}{R_{500}},\end{array}$$(11)

where r_cm is the centre-of-mass position of the cluster and r_c is the theoretical centre of the cluster, identified as the position of the highest density peak. Finally we considered the virial ratio, η, based on the virial theorem and defined as

$$\begin{array}{c}\hfill \eta =\frac{2T-E_{\mathrm{s}}}{|W|},\end{array}$$(12)

where T is the total kinetic energy, E_s is the surface pressure energy from both collisionless and gas particles, and W is the total potential energy (see Klypin et al. 2016, Cui et al. 2017, and John et al. 2019 for more details).

To obtain a continuous classification of the dynamical state of clusters, it is possible to combine these indicators using the following relation (Haggar et al. 2020):

$$\begin{array}{c}\hfill \chi =\sqrt{\frac{3}{(\frac{{\mathrm{\Delta }}_r}{0.04}{)}^2+(\frac{f_S}{0.1}{)}^2+(\frac{|1-\eta |}{0.1}{)}^2}}.\end{array}$$(13)

For a relaxed cluster, Δ_r and f_s are expected to be minimal, and η→ 1 (Cui et al. 2017). Therefore, they are expected to show high values of χ. Unfortunately, in the literature there is not a unique selection of the thresholds to use to segregate relaxed and disturbed clusters (see also Cui et al. 2017). The variety of choices made by different authors is partially explained either by the fact that different kinds of simulations were taken into account (e.g., dark matter versus hydrodynamical runs with different treatments for the baryon physics) or by the fact that different volumes (e.g., within R₅₀₀ or R₂₀₀) were used to estimate the dynamical state. By considering the most external regions, more sub-structures that are still in the process of merging may indeed be included, and the cluster result may be less virialised. For this reason, we decided to use the dynamical information provided by simulations to investigate the behaviour of the 12.7 % most relaxed and 21.1 % most disturbed objects of the sample (i.e. the 12.7% of objects with the highest and the 21.1 % of objects with the lowest values of χ). These percentages represent the fraction of relaxed and disturbed objects of the CHEX-MATE sample identified in Sect. 8.2.

8.4.2. Morphological analysis of the simulated sample

The X-ray images associated with each cluster are produced with Smac (Dolag et al. 2005) without the inclusion of any background component and are filtered in the 0.7–1.2 keV energy band. The pixel size of each map is 4 kpc and is fixed at all redshifts. In order to reproduce the XMM-Newton observations, we smoothed the simulated images with a Gaussian function of σ = 6 arcsec, which represents the FWHM of the XMM-Newton PSF³, and we binned them using the same scale of the observations (i.e. 1 pixel = 2.5 arcsec). Starting from these images, we then computed the morphological parameters described in Sect. 4 for each simulated cluster.

To properly compare the distribution of the morphological parameters obtained from the simulated images with the observed ones, we randomly extracted 10⁴ sub-samples, each consisting of 118 systems, from the 4692 X-ray observations of the 1564 simulated clusters. Each sub-sample was built with the aim of reproducing the distribution in mass and redshift of the CHEX-MATE objects and avoiding the selection of the same cluster in more snapshots of the simulations. In particular, we only considered the simulated objects reproducing the properties of Tier 1 and Tier 2, and we extracted a number of clusters corresponding to the number of the CHEX-MATE objects located at redshift $\frac{z_i-z_{i-1}}{2}<z<\frac{z_i+z_{i-1}}{2}$ from each snapshot, i. By computing the first, second, and third quartiles of the distribution of the morphological parameters for each extraction, we obtained the mean values of these quantities and their associated uncertainties (namely their standard deviation). In Table 3, we report both a comparison of this result with the one arising from the observations and also the p value of the KS test derived to compare the distributions of the observed and simulated morphological parameters (in this case the value and the uncertainties reported are also, respectively, the mean and the standard deviation of the values obtained for each randomly extracted sub-sample). We noticed a good agreement between the observed and simulated P₂₀, P₃₀, and w (high p-values of the KS test), while c shows borderline behaviour that is further discussed in Sect. 8.4.3. We also computed the statistical relative error associated with each parameter (defined as the ratio between the standard deviation and the mean of the the values obtained for each randomly extracted sub-sample), finding values of 5%, 10%, 25%, and 38% for c, w, P_20, and P_30, respectively.

For each randomly extracted sub-sample, we also evaluated the grade of correlation between these four parameters by estimating the Spearman coefficients (see Table 4). The results of this analysis are shown in the corner plot of Fig. 3, right panel. In particular, the couple of parameters for which the observed Spearman coefficient is within one standard deviation (1σ) from the mean of the simulated Spearman coefficients are represented in green, the Spearman coefficients between 1 and 3σ are shown in yellow, and the correlations over 3σ are given in red. This comparison highlights that all the couples show correlations included in the 3σ interval, with the exception of P₂₀–P_30, whose trend does not reproduce the one arising from the observations. This behaviour could be related to the fact that in simulations the signal is strong up to R₅₀₀.

As was done for the observations, in this case we also combined the four parameters to build the parameter M (see Sect. 7 for the definition). In Fig. 5 (right panel), we compare the distribution of this quantity with the distribution of the dynamical indicator χ. As it is possible to observe, a correlation between the two quantities is present (Spearman coefficient, r = 0.5). Furthermore, it is possible to observe that the majority of the objects classified as relaxed (or disturbed) with the definition presented in Sect. 8.4.1 are characterised by M < 0 (or M > 0), which is in agreement with the result obtained for the CHEX-MATE sample.

8.4.3. Discussion on the concentration

As presented in Sect. 8.4.2, some discrepancies are present between the observed and simulated distributions of the concentration. To investigate the nature of these differences, we first checked the validity of our procedure by comparing our result with the concentration estimated from surface brightness profiles by Bartalucci et al. (in prep.). In this latter analysis, the SB profiles were centred on the X-ray peak and then subtracted of the background and corrected for the vignetting. The concentration was then computed within two apertures of radii 0.1R₅₀₀ and R₅₀₀ using the following relation:

$$\begin{array}{c}\hfill C_{\mathit{SB}}=\frac{{\int }_0^{0.1R_{500}}{S}_X(r)r\mathrm{d}r}{{\int }_0^{R_{500}}{S}_X(r)r\mathrm{d}r}.\end{array}$$(14)

Eventually, the PSF is taken into account using the King function (e.g., Read et al. 2011):

$$\begin{array}{c}\hfill B(r)=\frac{A}{{[1+{(r/R_0)}^2]}^{\alpha }},\end{array}$$(15)

where R₀ = 0.088 arcmin is the core radius and α = 1.59 is the index. To be consistent with this analysis based on the profiles, we recomputed the values of the concentration using the same apertures. The results of this comparison is shown in Table 5. It is possible to observe that a good agreement is present between the two procedures, with the medians and the first and third quartiles being very similar in both cases. The lack of an influence of the PSF on the estimation of c (which is actually observed in literature) is probably related to the choice of the radii in which c is computed. As presented in Fig. 9, 0.1R₅₀₀ is larger than the dimension of the on-axis PSF for all the CHEX-MATE clusters, and consequently the photons located in the centre are not spread across larger regions. From this check, we verified that the origin of the discrepancy between the observed and simulated sample is not related to the method adopted for the estimation of c.

Fig. 9. Comparison between the dimension of the inner aperture used for the computation of c (i.e., 0.1R500 and the dimension of the XMM-Newton PSF; violet dashed line).

We then focused on the simulated sample and we noticed that many systems with high c values show a ring emission in their central region (i.e. r < 40 kpc; for an example, see Fig. E.1). This effect is related to the presence of particles with high SPH density in the central regions of the simulated clusters, due to the action of the isotropic feedback from AGNs. In these simulations, the AGN feedback is modelled as a single bubble that expands and compresses gas. Therefore, we re-computed c for both the simulated and observed sample excluding the central region of the cluster (i.e. r < 0.15R₅₀₀) and adopting 0.3 ⋅ R₅₀₀ and R₅₀₀ as inner and outer apertures. The results obtained are presented in Table 3 (last row). It is possible to observe that a difference between the observed and simulated distributions is still present. Therefore, we conclude that the discrepancies between the CHEX-MATE X-ray images and simulated maps are not limited to the description of the core. A more accurate analysis of the simulated and observed SB profiles is necessary to better investigate the behaviour of the c distributions.

8.4.4. Effects of the orientation

Finally, simulations are useful to understanding how the orientation of a cluster can influence the estimation of the morphological parameters and consequently the classification of its dynamical state. For each of the 1564 simulated clusters, we had at our disposal three images related to three different orientations of the object (namely, the orientations along the x, y, or z axes⁴). Therefore, we computed the morphological parameters for each orientation, and we estimated the dispersion Δ between the values obtained for the x projection and the values estimated along the other projections following the relation

$$\begin{array}{c}\hfill \mathrm{\Delta }=|{log}_{10}({\mathcal{P}}_x)-{log}_{10}({\mathcal{P}}_{y,z})|,\end{array}$$(16)

with 𝒫 representing a general parameter and the subscript representing the orientation considered (x, y, and z). The mean distribution of the dispersion along the y and z axis is shown in Fig. 10, first row. The median values of Δ obtained are 0.03, 0.17, 0.45, and 0.56 for c, w, P₂₀, and P₃₀, respectively.

Fig. 10. Comparison between parameters computed along the x direction (x-axis in the plot) and the y and z directions (grey and violet, respectively).

8.5. Robustness and efficiency of the parameters

Using the median value of the parameter Δ obtained from the investigation of the effects of the exposure time (see Sect. 6.1), radius (see Sect. 6.2), and orientation (see Sect. 8.4.4) on the estimates of the morphological parameters, we computed the associated relative error:

$$\begin{array}{c}\hfill ϵ={10}^{\mathrm{\Delta }}-1.\end{array}$$(17)

The results obtained are reported in Table 6. We noticed that the concentration is the most stable parameter, and it is also the indicator that shows the lowest uncertainties related to orientation effects and to the exposure time of the observations. Concerning the other parameters, it seems that the orientation is the major source of uncertainty, especially for the power ratios, which show percentage relative errors of 180% and 263% for P₂₀ and P₃₀, respectively. However, we have to take into account the fact that these two parameters cover a wide range of values. We found that for more than 93 % of the clusters the variation of the power ratios with the orientation is less than 10% of the total observed range of values. This result is similar to the one obtained by Jeltema et al. (2008).

The uncertainties related to the estimation of w, P_20, and P₃₀ increase with the decreasing exposure time. High values of the relative percentage errors are indeed obtained using images with t_exp equal to 5 ks. Given these results, we concluded that under the conditions of our analysis (i.e. high-quality observations and low uncertainties on the R₅₀₀ estimation), our four parameters could be considered stable.

For each morphological parameter, we assessed the threshold values above or below which the relaxed and disturbed systems lie (defined from the final classification, see Sect. 8.2). We found that the relaxed population is characterised by values of c > 0.49, w < 0.006, P₂₀ < 1.0 × 10⁻⁶, and P₃₀ < 0.4 × 10⁻⁷, while the disturbed population shows values of c < 0.19, w > 0.01, P₂₀ > 4.0 × 10⁻⁶ and P₃₀ > 0.5 × 10⁻⁷.

9. Summary and conclusions

In this paper, we present the morphological analysis performed on the X-ray images of the 118 clusters of the CHEX-MATE sample, and we provide a classification of their dynamical state. To achieve this aim, we investigated the behaviour of four morphological parameters commonly used in literature to identify the most relaxed and most disturbed systems. These parameters are the concentration, c, the centroid shift, w, and the power ratios P₂₀ and P₃₀. To verify the ability of these indicators in reproducing the classification realised by trained eyes, we first realised a visual classification of the sample. In particular, seven astronomers inspected images by eye, assigning a grade from 0 (most relaxed) to 2 (most disturbed) to each object. By averaging these seven classifications, we found that 19 clusters are defined as more relaxed systems and 37 as more disturbed. All the other systems (62) are not clearly classifiable and are thus defined as mixed. Using this classification as a reference for our analysis, our conclusions are the following:

• The distributions of all four morphological parameters do not show signs of bimodality (see Sect. 8.1).

• The median (first-third quartile) of the distributions of the morphological parameters are, respectively, 0.29 (0.21–0.43) for c; 0.011 (0.005–0.024) for w; 2.6 × 10⁻⁶ (1.0–7.5 × 10⁻⁶) for P_20, and 1.3 × 10⁻⁷ (0.3–5.8 × 10⁻⁷) for P₃₀ (see e.g., Table 3).

• The four morphological parameters are characterised by strong correlations (i.e. Spearman coefficient r > 0.5, see Sect. 5 and Fig. 3, left panel), with the tails of their distributions that are populated by the most relaxed and most disturbed systems identified by the visual classification.

• These four morphological parameters are robust enough to be not biased under the conditions of our analysis. They have proved to be not influenced by the exposure time of the observations and by uncertainties related to the definition of the region where they are computed. The systematic relative errors due to the exposure time are 0.7%, 7.2 %, 12%, and 32% for c, w, P_20, and P₃₀, respectively, when the exposure time of the observations is halved (while it is 0.9%, 15%, 26%, and 58% for c, w, P_20, and P₃₀, when the exposure time of the observations is fixed at 5 ks; see Sect. 6.1). Concerning the assumption of R₅₀₀, we found systematic relative errors of ∼5%, ∼5%, ∼15%, and ∼32% for c, w, P_20, and P₃₀, respectively (see Sect. 6.2).

• These values are, on average, lower than (or comparable with) the mean statistical relative errors of 5%, 10%, 25%, and 38% for c, w, P_20, and P_30, respectively (see Sect. 8.4.2), confirming in particular how the concentration is less prone to systematic effects.

• We combined the parameters in a single quantity, M, that is able to assign to each object of the sample a grade of relaxation. We also noticed that this new parameter does not show signs of bimodality. We then realised a continuous classification of the dynamical state of the sample based on the values of the M parameter (see Sect. 7).

• By comparing the visual classification and the classification based on the M parameter, we finally classified the objects for which the two classifications are in agreement as relaxed and as disturbed. Using this criteria, we found that 15 systems (12.7 %) are relaxed and 25 (21.1 %) are disturbed (see Sect. 8.2).

• According to the classification described in the above point, we found that the relaxed population is characterised by values of c > 0.49, w < 0.006, P₂₀ < 1.0 × 10⁻⁶, and P₃₀ < 0.4 × 10⁻⁷, while the disturbed population shows values of c < 0.19, w > 0.01, P₂₀ > 4.0 × 10⁻⁶, and P₃₀ > 0.5 × 10⁻⁷ (see Sect. 8.5).

• From the comparison of our analysis with previous works related to X-ray-selected samples, we found that the CHEX-MATE objects tend to be more dynamically disturbed (i.e. high centroid shift and low concentration) than the X-ray-selected samples, in agreement with what has been obtained by other recent studies (see Sect. 8.3).

• By repeating our analysis on a simulated sample provided by THE THREE HUNDRED collaboration, we found that a good agreement is present for what concerns the level of correlation between the parameter pairs (see e.g., Fig. 3, right panel). However, the comparison between the observed and simulated distributions of the morphological parameters highlighted the presence of a discrepancy: the simulated values of the concentration are on average higher than the observed ones (see Sect. 8.4.2).

• By investigating the behaviour of c, we found that simulations are characterised by a distribution with higher values due to the non-negligible presence of particles with high SPH density produced from the action of the isotropic thermal AGN feedback in the central regions of the simulated objects. This behaviour indicates that a more realistic and higher resolution feedback model is required (such as a circulating mechanical anisotropic AGN feedback; see Sect. 8.4.3). However, the discrepancy between the observed and disturbed distribution remains even if the core region is masked when computing the morphological parameters. Therefore, further analyses are necessary to investigate the properties of the observed and simulated SB profiles.

• Thanks to simulations, we investigated the systematic relative error associated with the orientation of the cluster, finding values of 7%, 48%, 180%, and 263% for c, w, P_20, and P₃₀ (see Sect. 8.4.2), which represent the dominant component in their error budget, as expected for quantities based on the projected distribution of X-ray counts in the plane of the sky.

This study of the X-ray morphological properties of the CHEX-MATE objects provides the parameters that will be used in forthcoming analyses to assess the role of the dynamical state in the reconstruction and characterisation of their intrinsic physical quantities, from the thermodynamic profiles to their distribution in the scaling relations.

Acknowledgments

S.E., L.L., I.B., M.R., F.G., S.G. and S.M. 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), INAF mainstream project 1.05.01.86.10, and funding from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158). G.W.P. acknowledges support from CNES, the French space angency.

References

Appendix A: Parameter values

All the parameter values used in this paper and calculated within R₅₀₀ are listed in Table A.1. The objects marked with an asterisk are those for which the estimation of the morphological parameters may not be accurate. These systems are: G048.10+57.16, for which the emission inside R₅₀₀ is not fully covered if we choose the X-ray peak as centre; G283.91+73.97, for which the characterisation of the background is complex, since it is located behind Virgo; G028.63+50.15, whose X-ray emission is very complex and may not be described accurately; G107.10+65.32 and G124.20-36.48, which show a disturbed morphology due to the ongoing merger between two sub-structures of similar dimensions. Since a double X-ray peak is present, and since it is not possible to identify the principal cluster, the parameters could be influenced by the choice of the X-ray peak.

Appendix B: Other morphological parameters

B.1. Ellipticity

The ellipticity parameter is able to describe the shape of the surface brightness distribution and could provide a link to the dynamical state. Indeed, relaxed systems are expected to be rounder than the disturbed clusters. The ellipticity parameter is defined as

$$\begin{array}{c}\hfill \eta =1-\frac{b}{a},\end{array}$$(B.1)

where a and b are, respectively, the major and the minor semi-axes. Here, for simplicity, we only consider the axial ratio (i.e. b/a), and consequently η takes values from 0 (elliptical shape) to 1 (circular shape). For the estimation of a and b, we diagonalised the 2×2 inertia tensor T, defined as:

$$\begin{array}{c}\hfill T_{k,l}={\displaystyle \sum _{i=1}^n(r_k{r}_l)S{B}_i}\end{array}$$(B.2)

where k, l = 1,  2, r₁ = x, r₂ = y, SB_i is the surface brightness computed in the i-th pixel, and n is the number of pixels covering the cluster emission. The origin of the coordinate axis is the X-ray peak of the cluster. Denoting the eigenvalues of I as λ_1, 2 (so that $|T-1{\lambda }_{1,2}^2|=0$), we defined a as λ_major and b as λ_minor.

Observing the correlation of η with the other parameters, we noticed a strong correlation from η – P₂₀ (r = -0.81), which highlights that these two parameters represent the same quantity, that is the cluster ellipticity. For this reason, we decided to only consider P₂₀ for the rest of the analysis and to report the correlations of η with the other parameters in this appendix. For the sake of completeness, in Fig. B.1 the distribution of η for the observed and simulated samples is reported.

Fig. B.1. Distribution of η for the observed (violet solid line) and simulated (black dashed line) samples. The vertical lines represent the medians of the two distributions.

B.2. Asymmetry and smoothness

In this sub-section, we present the analysis realised using the asymmetry and the smoothness, which are two parameters originally used to study the morphology of galaxies. The asymmetry parameter, A, is a measure of how the light distribution differs from a spherically symmetric distribution and is computed by subtracting an image rotated by 180° from the original image (Conselice et al. 2000; Lotz et al. 2004):

$$\begin{array}{c}\hfill A=\frac{{\sum }_{i,j}|\text{I}(i,j)-{\text{I}}_{180}(i,j)|}{{\sum }_{i,j}|\text{I}(i,j)|},\end{array}$$(B.3)

where I is the cluster’s image, I₁₈₀ is the image rotated by 180° around the pixel corresponding to the X-ray peak.

The smoothness parameter, S, is obtained by subtracting a Gaussian-smoothed image, I_s, from the original one, I (see also Lotz et al. 2004; Conselice et al. 2003):

$$\begin{array}{c}\hfill S=\frac{{\sum }_{i,j}|\text{I}(i,j)-{\text{I}}_{\text{s}}(i,j)|}{{\sum }_{i,j}|\text{I}(i,j)|},\end{array}$$(B.4)

$$\begin{array}{c}\hfill S_{\text{bkg}}=\frac{{\sum }_{i,j}|\text{B}(i,j)-{\text{B}}_{\text{s}}(i,j)|}{{\sum }_{i,j}|\text{I}(i,j)|}.\end{array}$$(B.5)

In order to compare S values of clusters spanning a wide range of redshifts, we decided to suppress inhomogeneities smaller than the minimal scale resolved in high-redshift images, that is ∼ 100 kpc. Thanks to this condition, S is sensitive to the same sub-clusters scales, independently of the redshift of the considered objects.

We then computed these two parameters for the CHEX-MATE and simulated sample, and we report the results in Table B.1. In particular, concerning the simulated sample, the values reported are obtained with the same method adopted in Sect. 8.4.2 (i.e. building 10⁴ randomly extracted sub-samples, computing the median value of the Spearman coefficient and the corresponding standard deviation). It is possible to observe that the observed median and Spearman coefficient differ from the simulated one. Furthermore, a very strong correlation is observed from A – S, with r = 0.94. The discrepancies between the observed and simulated distributions of the asymmetry and smoothness could be an indications that these two parameters are strongly influenced by the signal-to-noise ratio (S/N), which is equal to S/N ∼ 150 for the observations and S/N ∼ ∞ for the simulations. To test this hypothesis, we considered a sub-sample of 20 CHEX-MATE clusters, and we repeated the morphological analysis on images with both halved t_exp and t_exp= 5 ks. In Figure B.2, we show the comparison between the distribution of A and S computed using the original images (solid line) and the distributions obtained using images with a reduced exposure time (dotted line for ∼ 50 % t_ext and dashed line for t ∼ 5 ks). As it is possible to notice, these two parameters show a shift in their distributions, which increase as the exposure time decreases.

Fig. B.2. Comparison between the distribution of A and S computed using images with the original exposure time (solid line), with halved exposure time (dotted line) or with an exposure time of 5 ks. The bottom plots show the results obtained when no threshold is applied, the upper plots instead show the results obtained when only pixels with value higher than 0 are considered.

Before drawing any conclusions, we make the following consideration: originally, the asymmetry and the smoothness were defined to classify the morphology of galaxies in the optical band, where images are characterised by a larger number of photon counts. X-ray images, instead, are characterised by a sparse emission. When computing A (or S) for each pixel, the absolute value of the difference between the original and the rotated (or smoothed) image is realised. The local difference between pixels with and without counts could thus result in high values of the asymmetry (or smoothness), even if the overall shape of the emission is still symmetric (as already observed by Nurgaliev et al. 2017). The difference observed is not due to a real decrease of the flux, but due to the sparse distribution of the emission. This consideration could explain why the Spearman coefficient of the A and S pair is so high: they do not quantify the shape of the X-ray emission, but the number of pixels with no emission. We thus decided to recompute these two parameters by only considering the pixels with values higher than 0, and we report the new results in Table B.1. The Spearman coefficient obtained for the A – S pair is consistent with what was obtained by Parekh et al. (2015); however, discrepancies are still present between the observed and simulated S. Given these results, we prefer not to include this parameter in our analysis.

Appendix C: Morphological parameters at 0.5 R₅₀₀

In this section, we report the result of the morphological analysis realised using a region with a radius of 0.5 R₅₀₀. First of all, we compared the values of the morphological parameters obtained using these new regions, with the values estimated inside R₅₀₀. The results are shown in Fig. C.1. It is possible to observe that strong correlations (r>0.5) are obtained for all the parameters. We then realised a cluster analysis as the one realised in Sect. 8.1, and we present the results in Table C.1 and in Fig. C.2. The Gaussian mixture model reveals the presence of a single component for all the parameters, with the exception of P₂₀. For this latter indicator, two components are identified, suggesting that a bimodality may characterise its distribution. However, it is possible to observe that the Weibull function has a highest BIC value. In particular, since the discrepancy between the Weibull and Gaussian BIC values is ∼ 6, it is possible to assess that the Weibull model is favoured (see Sect. 8.1 for more details). Therefore, we conclude that the four parameters computed within 0.5 R₅₀₀ do not show signs of bimodality. Given the BIC values, we found that the single Gaussian component model is favoured in describing the distribution of the concentration. For the other parameter, it is not possible to unambiguously identify the best-fit model.

Fig. C.1. Correlations between the morphological parameters estimated inside R500 (x-axis) and 0.5 ⋅ R500 (y-axis).

Fig. C.2. Distribution of c, w, P20, and P30 computed inside 0.5 ⋅R500. The solid line represents the model obtained from the clustering analysis, while the dashed line represents the Weibull function. The parameters c, w, P20, and P30 were multiplied by 103, 104, 109, and 1010, respectively.

Appendix D: Principal component analysis

Some parameters may show good correlations because they essentially provide the same type of information. In order to recognise redundancies and determine the minimal number of dimensions able to describe the properties of the observed clusters, we applied the so-called principal component analysis (PCA; Pearson 1901; Hotelling 1933; Jolliffe et al. 2011) to our dataset. This procedure is used to represent an original set of n₀ mutually correlated random variables (in this case, our four morphological parameters) with a smaller set (n < n₀) of independent hypothetical variables, and it allows us to reduce the number of dimensions able to describe a dataset, without much information loss.

Before starting this analysis, we realised a standardisation of the dataset, which consists of re-scaling the values of each parameter in order to obtain a standard normal distribution, with a mean of zero and a standard deviation of one. We then computed the eigenvectors and eigenvalues of the covariance matrix; the eigenvectors represent the axis of greatest variance, while the eigenvalues associated with the eigenvectors indicate their magnitude (and thus the variance in the direction of the eigenvector). By ranking the eigenvectors in order of their eigenvalues (highest to lowest), it is possible to find the principal components in order of significance.

The results obtained using P₂₀, P₃₀, c, and w as input parameters are reported in Table D.1 and shown in Fig. D.1, left panel. We noticed that two components are sufficient to explain ∼ 90 % of the variance and that the relaxed, disturbed, and mixed systems (as defined by the visual classification, see Sect. 3) occupy different regions in the plot. For completeness, we also tested this procedure using only c and w as input parameters, which are considered in literature among the most powerful parameters for the identification of the dynamical state of clusters. We found that the so-called principal component 1, which is essentially a weighted combination of c and w, is able to clearly identify the three populations of objects (Fig. D.1). The behaviour observed confirms that the combination of the concentration and of the centroid shift is particularly suitable for the identification of the dynamical state of clusters, as already observed by Lovisari et al. (2017).

Fig. D.1. Results obtained by applying PCA to observations and using c, w,P20, and P30 (left panel), or only c, w (right panel), as input parameters. The colour scale represents the grade of relaxation of a system (with those rated 0 being the most relaxed systems, and those rated 3 being the most disturbed ones) defined by the visual classification (see Sect. 3).

Appendix E: Ring emission in simulations

In Fig. E.1, the ring emission observed in the central regions (i.e. r< 40 kpc) of the simulated clusters is shown. As presented in Sect. 8.4.3, this effect is related to the isotropic model used for the description of the AGN feedback and has an impact on the estimation of the concentrations of the simulated sample.

Fig. E.1. Example of ring emission observed in the central regions of simulated clusters. The green circle represents a circular region of radius r = 40 kpc.

All Tables

All Figures

	Fig. 1. Example of relaxed (left), mixed (centre), and disturbed (right) systems identified by the visual classification. The horizontal and vertical dimensions of the boxes are equal to 2R500.
In the text

	Fig. 2. Distribution in the M500 − z plane of the 118 CHEX-MATE clusters. The colour scale represents the dynamical state obtained from the visual classification. The shaded area indicate the Tier 1 and Tier 2 redshift ranges in grey and violet, respectively.
In the text

Fig. 3. Distribution of the morphological parameters for the CHEX-MATE (left) and simulated (right) samples. The r values reported in the boxes are the Spearman coefficients of the considered parameter pairs. Left: the different colours represent the dynamical state obtained from the visual classification (see Sect. 3). Right: the different colours represent the dynamical state obtained using the dynamical indicator χ (see Sect. 8.4.1). The background colours of the plots represent the level of agreement with the observations: green is for 1σ, yellow is for 1–3σ, and red is for a level of agreement over 3σ.

In the text

Fig. 4. Analysis of the robustness of the morphological parameters. First row: comparison between the parameters estimated from the original images (x-axis) and the parameters estimated from images with halved texp (violet) and with texp = 5 ks (grey). Second row: comparison between the parameters estimated using R500 (x-axis) and the parameters computed in a region of r = 0.95R500 (violet) and r = 1.05R500 (grey).

In the text

Fig. 5. Distribution of the M parameter. Left: comparison between the CHEX-MATE dynamical state obtained from the visual classification and from the M parameter. The colours represent the dynamical classification obtained on the basis of the M parameter (blue for relaxed, red for disturbed, and green for mixed clusters). Right: distribution of the simulated sample in the χ − M space. The colours represent the relaxed (blue), disturbed (red), and mixed (green) class defined on the basis of χ (see Sect. 8.4.1).

In the text

	Fig. 6. Distribution of c, w, P20, P30, and M. The solid line represents the model obtained from the clustering analysis, while the dashed one represents the Weibull function. The parameters c, w, P20, and P30 were multiplied by 103, 104, 109, and 1010, respectively. For the M parameter, we added 10 to the original values instead.
In the text

	Fig. 7. Mass and redshift distributions of the morphological parameters. The different colours represent the Tier 1 (grey) and Tier 2 (violet) objects. In the boxes is reported the Spearman coefficient, r.
In the text

	Fig. 8. Left: comparison between the concentration values of the CHEX-MATE and ME-MACS sample. Right: comparison between the centroid shift values of the CHEX-MATE, REXCESS, 400d and SPT samples.
In the text

	Fig. 9. Comparison between the dimension of the inner aperture used for the computation of c (i.e., 0.1R500 and the dimension of the XMM-Newton PSF; violet dashed line).
In the text

	Fig. 10. Comparison between parameters computed along the x direction (x-axis in the plot) and the y and z directions (grey and violet, respectively).
In the text

	Fig. B.1. Distribution of η for the observed (violet solid line) and simulated (black dashed line) samples. The vertical lines represent the medians of the two distributions.
In the text

	Fig. B.2. Comparison between the distribution of A and S computed using images with the original exposure time (solid line), with halved exposure time (dotted line) or with an exposure time of 5 ks. The bottom plots show the results obtained when no threshold is applied, the upper plots instead show the results obtained when only pixels with value higher than 0 are considered.
In the text

	Fig. C.1. Correlations between the morphological parameters estimated inside R500 (x-axis) and 0.5 ⋅ R500 (y-axis).
In the text

	Fig. C.2. Distribution of c, w, P20, and P30 computed inside 0.5 ⋅R500. The solid line represents the model obtained from the clustering analysis, while the dashed line represents the Weibull function. The parameters c, w, P20, and P30 were multiplied by 103, 104, 109, and 1010, respectively.
In the text

	Fig. D.1. Results obtained by applying PCA to observations and using c, w,P20, and P30 (left panel), or only c, w (right panel), as input parameters. The colour scale represents the grade of relaxation of a system (with those rated 0 being the most relaxed systems, and those rated 3 being the most disturbed ones) defined by the visual classification (see Sect. 3).
In the text

	Fig. E.1. Example of ring emission observed in the central regions of simulated clusters. The green circle represents a circular region of radius r = 40 kpc.
In the text