The Most Massive galaxy Clusters (M2C) across cosmic time: link between radial total mass distribution and dynamical state

# The Most Massive galaxy Clusters (M2C) across cosmic time: link between radial total mass distribution and dynamical state

I. Bartalucci AIM, CEA, CNRS, UniversitÃ© Paris-Saclay, UniversitÃ© Paris Diderot, Sorbonne Paris CitÃ©, F-91191 Gif-sur-Yvette, France    M. Arnaud AIM, CEA, CNRS, UniversitÃ© Paris-Saclay, UniversitÃ© Paris Diderot, Sorbonne Paris CitÃ©, F-91191 Gif-sur-Yvette, France    G. W. Pratt AIM, CEA, CNRS, UniversitÃ© Paris-Saclay, UniversitÃ© Paris Diderot, Sorbonne Paris CitÃ©, F-91191 Gif-sur-Yvette, France    J. DÃ©moclès AIM, CEA, CNRS, UniversitÃ© Paris-Saclay, UniversitÃ© Paris Diderot, Sorbonne Paris CitÃ©, F-91191 Gif-sur-Yvette, France    L. Lovisari Center for Astrophysics — Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA
Submitted to A&A on 29/05/2019
###### Key Words.:
intracluster medium – X-rays: galaxies: clusters – Dark matter

We study the dynamical state and the integrated total mass profiles of massive () SZ-selected clusters at . The sample is built from the Planck catalogue, with the addition of four SPT clusters at . Using XMM-Newton imaging observations, we characterise the dynamical state with the centroid shift , the concentration , and their combination, , which simultaneously probes the core and the large scale gas morphology. Using spatially-resolved spectroscopy and assuming hydrostatic equilibrium, we derive the total integrated mass profiles. The mass profile shape is quantified by the sparsity, the ratio of to , the masses at density contrast 500 and 2500, respectively. We study the correlations between the various parameters and their dependence on redshift. We confirm that SZ-selected samples, thought to reflect most closely the underlying cluster population, are dominated by disturbed and non-cool core objects at all redshifts. There is no significant evolution or mass dependence of either the cool core fraction or the centroid shift parameter. The parameter evolves slightly with , having a correlation coefficient of and a null hypothesis -value of . In the high mass regime considered here, the sparsity evolves minimally with redshift, increasing by between and , an effect significant at less than . In contrast, the dependence of the sparsity on dynamical state is much stronger, increasing by a factor of from the 1/3 most relaxed to the most disturbed objects, an effect significant at more than . This is the first observational evidence that the shape of the integrated total mass profile in massive clusters is principally governed by the dynamical state, and is only mildly dependent on redshift. We discuss the consequences for the comparison between observations and theoretical predictions.

## 1 Introduction

The shape of the dark matter profile in galaxy clusters is a sensitive test of the nature of dark matter and of the theoretical scenario of structure formation. In the standard framework, cosmological structures form hierarchically from initial density fluctuations that grow under the influence of gravity. In a cold dark matter (CDM) Universe, the dark matter (DM) collapse is scale-free and one expects objects to form with similar internal structure. The DM shape is expected to depend on the halo assembly history, which is a function of the redshift, the total mass, and the underlying cosmology (e.g. Dolag et al. 2004; Kravtsov & Borgani 2012). A certain scatter of the shape and break of strict self-similarity are expected, reflecting the detailed formation history of each halo.

Clusters of galaxies are ideal targets to test the above scenario: the dark matter is the dominant component by far except in the very centre. Furthermore, complementary techniques can measure the total mass density profile, e.g., galaxy velocities; strong gravitational lensing in the centre and weak lensing at large scale; X-ray estimates using the intra-cluster medium (ICM) density and temperature profiles and the hydrostatic equilibrium (HE) equation (see Pratt et al. 2019, for a review).

Numerical simulations indeed predict that the cold dark matter density profiles of virialised objects follows a ’universal’ form. Well-known parametric models for the DM density profiles include those proposed by Navarro et al. (1997, herafter, NFW) and the Einasto profile (Einasto 1965), which is currently considered to be a more accurate description of the profiles in state-of-the-art simulations (Navarro et al. 2004; Wu et al. 2013).

A fundamental parameter of these parametric models is their concentration, which in general terms describes the relative distribution in the core as compared to the outer region (Klypin et al. 2016). It is characterised in the NFW and Einasto111In this model, the actual shape also depends on a second parameter. models by the ratio c222 is defined as the radius enclosing times the critical matter density at the cluster redshift. is the corresponding mass., where is the radius at which the logarithmic density slope is equal to . The relation between the concentration and the total mass (hereafter c-) has been very widely used as an indicator of the dark matter profile shape in cosmological simulations (e.g. Diemer & Kravtsov 2015, and references therein). However, it has been shown that the c- relation of the most relaxed haloes is different to that derived for the full population (e.g. Neto et al. 2007; Bhattacharya et al. 2013). In fact, the capability of these parametric models to reproduce the DM profile (i.e. the goodness of the fit) depends on the dynamical state of the halo and on the formation time (Jing 2000; Wu et al. 2013, for the NFW case). Results on the dependence of the DM profile shape on mass and redshift based on c- relations may thus be ambiguous, and depend on the sample selection (Klypin et al. 2016; Balmès et al. 2014). For this reason, recent works based on simulations have often focused on relaxed haloes (e.g. Dutton & Macciò 2014; Ludlow et al. 2014; Correa et al. 2015).

However, comparison with cluster observations cannot rigorously be made, as there is a continuous distribution of dynamical states, and the definition of what is a relaxed object is therefore somewhat arbitrary. Moreover, in the absence of fully realistic simulations including the complex baryonic physics, it is nearly impossible to define common criteria that quantify the dynamical state, in a consistent manner, both in simulations and observations. Ideally we would compare the full cluster population from numerical simulations to observed samples chosen to reflect as closely as possible the true underlying population. The advent of Sunyaev-Zeldovich (SZ)-selected cluster catalogues offers a unique opportunity to build such observational samples. Surveys such as those from the Atacama Cosmology Telescope (ACT, Marriage et al. 2011), the South Pole Telescope (SPT, Reichardt et al. (2013) and Bleem et al. 2015) and the Planck Surveyor (Planck Collaboration VIII 2011, Planck Collaboration XXIX 2014, and Planck Collaboration XXVII 2016) have provided SZ-selected cluster samples up to . The magnitude of the SZ effect being closely linked to the underlying mass with small scatter (da Silva et al. 2004), these are thought to be as near as possible to being mass selected, and as such unbiased.

The observational study of the DM profile shape using such samples requires investigation of the dependence on fundamental cluster quantities such as mass, dynamical state, and redshift. X-ray observations are an excellent tool with which to undertake such studies – the ICM morphology can be used to infer the dynamical state, while the total mass profile can be derived by applying the hydrostatic equilibrium (HE) equation to spatially-resolved density and temperature profiles. While this method yields the highest statistical precision on individual profiles over a wide radial range and up to high (Amodeo et al. 2016; Bartalucci et al. 2018), it has the drawback of a systematic uncertainty due to any departure of the gas from HE, which must be taken account.

Here we apply the sparsity parameter, introduced by Balmès et al. (2014), to quantify the shape of the total mass profiles derived from the X-ray observations. The sparsity is defined as the ratio of the integrated mass at two overdensities. This non-parametric quantity is capable of efficiently characterising the profile shape, as long as the two overdensities are separated enough to probe the shape of the mass profile (Balmès et al. 2014; Corasaniti et al. 2018). Formally the DM profiles can be derived by subtracting the gas and galaxy distribution from the total. However in the following we focus on the total distribution, in view the negligible impact of the baryonic component on the total profile shape for the halo masses and density contrasts under consideration (e.g. Velliscig et al. 2014, Fig. 2).

In this work, we present the dynamical properties and the individual spatially resolved radial total mass profiles of a sample of 75 SZ-selected massive clusters in the redshift range with . We discuss the dispersion and evolution of the total mass profile shape, and demonstrate the link between the diversity in profile shape and the underlying dynamical state. In Sect. 2 we present the sample; in Sects. 3 and 4 we describe the methodology used to derive the HE mass profiles and the morphological parameters of each cluster, respectively; in Sect. 5 we discuss the morphological properties of the sample and its evolution; in Sect. 6 we investigate the dependence of the profile shape on mass, redshift, and dynamical state using the sparsity; finally, in Sects. 7 and Sect. 8 we discuss our results and present our conclusions.

We adopt a flat -cold dark matter cosmology with , , km Mpc s, and throughout. Uncertainties are given at the 68 % confidence level (). All fits were performed via minimisation.

## 2 The sample

### 2.1 The high-z SZ-selected sample

Our initial sample is built from the clusters with spectroscopic in the sample of the XMM-Newton Large Programme (LP) ID069366 (with re-observation of flared targets in ID72378), consisting of clusters detected at high with Planck, and confirmed by autumn 2011 to be at . The exposure times of these observations were optimised to allow the determination of spatially-resolved radial total mass profiles at least up to . To extend the redshift coverage to we included the five most massive clusters detected by SPT and Planck in the redshift range . Using deep XMM-Newton and Chandra observations, Bartalucci et al. (2017) and Bartalucci et al. (2018) examined the X-ray properties of these objects and determined ICM and HE total mass profiles up to . The resulting observation properties are detailed in Table 1 of Bartalucci et al. (2017).

Combining the above, we obtain the full ‘High-z’ sample of 33 objects, with , shown with orange open squares in Fig. 1. The observation details of this sample are detailed in Table 2.

### 2.2 The low-z SZ selected sample

Any evolution study requires a local reference sample with similar selection and quality criteria to act as as an ‘anchor’ to compare to high redshifts. The ongoing (AO17) XMM-Newton heritage program ‘Witnessing the culmination of structure formation in the Universe’ (PIs M. Arnaud and S. Ettori), based on the final Planck PSZ2 catalogue, will serve this function in the future (observations are to be completed by 2021). In the interim, for the present study, we use published XMM-Newton follow-up of Planck-selected local clusters taken from the ESZ sample.

The early SZ (ESZ, Planck Collaboration VIII 2011) catalogue represents the first release derived from the Planck all-sky SZ survey, containing clusters mostly at . There are a number of studies in the literature describing the X-ray properties of this sample. In particular, Lovisari et al. (2017) characterised the global properties and the morphological state of the ESZ clusters covered by XMM-Newton observations. We use their results on the morphological properties of the ESZ clusters at for which is within the field of view, and the relative error on the morphological parameters is less than . This is about of the corresponding parent ESZ sub-sample defined with the same and size criteria, thus we do not expect any major bias due to the incomplete XMM-Newton coverage. These objects cover one decade in mass, , and are shown with blue open squares in Fig. 1. Henceforth we refer to this sample as the ‘Low-z’ sample.

The spatially-resolved thermodynamic properties of an ESZ sub-sample were analysed by the Planck Collaboration, who combined X-ray and SZ data to calibrate the local scaling relations (Planck Collaboration XI 2011) and measure the pressure profiles (Planck Collaboration V 2013). We use the published thermodynamic profiles of clusters for which is within the field of view (the sub-sample ”A” defined in Sect. 3.1 of Planck Collaboration XI 2011) to derive the total mass profiles. These clusters are shown as filled blue squares in Fig. 1. We refer to this sample as the ”Low-z PXI” sample henceforth. Its representativity with respect to the full Low-z sample is excellent, as discussed in Appendix B.

### 2.3 Data preparation

The observations used in this work were taken using the European Photon Imaging Camera (EPIC, Turner et al. 2001 and Strüder et al. 2001) instrument on board the XMM-Newton satellite. This instrument is composed of three CCD arrays, namely MOS1, MOS2, and PN, which simultaneously observe the target. Datasets were reprocessed using the Science Analysis System333cosmos.esa.int/web/xmm-newton (SAS) pipeline ver. and calibration files as available in December . Event files with the this calibration applied were produced using the emchain and epchain tools.

The reduced datasets were filtered in the standard fashion. Events for which the keyword PATTERN is and for MOS, and PN cameras, respectively, were filtered out from the analysis. Flares were removed by extracting a light curve, and removing from the analysis the time intervals where the count rate exceeded times the mean value. We created the exposure map for each camera using the SAS tool eexpmap. We merged multiple observation datasets of the same object, if available. We report in Table 2 the effective exposure times after all these procedures. We corrected vignetting following the weighting scheme detailed in Arnaud et al. (2001). The weight for each event was computed by running the SAS evigwieght tool on the filtered observation datasets.

Point sources were identified using the Multi-resolution wavelet software (Starck et al. 1998) on the exposure-corrected , keV and keV images. We inspected each resulting list by eye to check for false detections and missed sources. We defined a circular region around each detected point sourceand excised these from the subsequent analysis.

We also defined regions encircling obvious sub-structures. Identified by eye, these regions were considered in the morphological analysis because the parameters we used to characterise the morphological state of a cluster (Sect. 4) are sensitive to the presence of any such sub-structures. However, these regions were excluded in the radial profile analysis detailed in Sect. Sect. 3.

X-ray observations are affected by instrumental and sky backgrounds. The first is due to the interaction of the instrument with energetic particles. The second is caused by Galactic thermal emission and the superimposed emission of all the unresolved point sources, namely the cosmic X-ray background (Lumb et al. 2002; Kuntz & Snowden 2000). These components were estimated differently for the radial profile 1D analysis and for the morphological 2D analysis, as described in Sects. 3 and 4, respectively.

### 3.1 Instrumental background estimation

We evaluated the instrumental background for the radial profile analysis following the procedures described in Pratt et al. (2010). Briefly, observations taken with the filter wheel in CLOSED position were renormalised to the source observation count rate in the , keV band for the EMOS and PN cameras, respectively. We then projected these event lists in sky coordinates to match our observations. We applied the same point source masking and vignetting correction to the CLOSED event lists as for the source data. We also produced event lists to estimate the out-of-time (OOT) events using the SAS-epchain tool.

### 3.2 Density and 3D temperature profiles

To determine the radial profiles of density and temperature of the ICM we followed the same procedures and settings detailed in Sect. 3.2 and Sect. 3.3 of Bartalucci et al. (2017). Briefly, we firstly determined the X-ray peak by identifying the peak of the emission measured in count-rate images in the [0.3-2.5] keV band smoothed using a Gaussian kernel of width between 3 and 5 pixels of width. We extracted the vignetted-corrected and background subtracted surface brightness profiles, from concentric annuli centred on the X-ray peak from both source and background event lists. These profiles were used to derive the radial density profiles, , employing the deprojection and PSF correction with regularisation technique described in Croston et al. (2006).

We obtained the deprojected temperature profiles by performing the spectral analysis described in detail in Pratt et al. (2010) and Sect. of Bartalucci et al. (2017). The background-subtracted spectrum of a region free of cluster emission was fitted with two unabsorbed MeKaL thermal models plus an absorbed power-law with fixed slope of . The resulting best-fitting model, renormalised by the ratio of the extraction areas, was then added as an extra component in each annular fit. The cluster emission was modelled by an absorbed MeKaL model with given in Table 2, using the absorption cross sections from Morrison & McCammon (1983). Spectral fitting was performed using XSPEC ver. 12.8.2. The deprojected 3D temperature profile, , was derived from the projected profile using the ”Non parametric-like” technique described in Sect. 2.3.2 of Bartalucci et al. (2018).

For the Low-z PXI sample, we used the temperature and density profiles published in Planck Collaboration V (2013) and Planck Collaboration XI (2011), which were derived employing identical methods to those used in this work.

### 3.3 Global properties

We determined the mass at density contrast , , and corresponding radius, from the mass proxy . This was computed iteratively from the  relation, calibrated by Arnaud et al. (2010) using HE mass estimates of local relaxed clusters. We assumed that the  relation obeys self-similar evolution. is defined as the product of the temperature measured in the region and the gas mass within (Kravtsov et al. 2006), the gas mass profiles being computed from the density profiles. The mass for each cluster and associated errors are reported in Table 2.

We used the published by Lovisari et al. (2017) for the Low-z sample, and the computed by us for the Low-z PXI and High-z samples. We investigated the coherence between these measurements by comparing the masses for the clusters in common of the Low-z PXI. The excellent agreement between the two is discussed in Appendix A, and the comparison is shown the right panel of Fig. 7.

### 3.4 Derivation of the total mass profiles

Under the assumption that the ICM is in hydrostatic equilibrium in the gravitational potential well the relation between the total halo mass within radius and the ICM thermodynamic properties is:

 M(≤R)=−kT(R)RGμmH[dlnne(R)dlnR+dlnT(R)dlnR], (1)

where G is the gravitational constant, is the hydrogen atom mass, and is the mean molecular weight in a.m.u. We used the deprojected density and temperature profiles and the relation in Eq. 1 to derive the mass profiles, applying the ‘forward non-parametric-like technique’ detailed in Sect. 2.4.1 of Bartalucci et al. (2018). The ‘forward’ in the technique name is due to the fact that we started our analysis from the surface brightness and projected temperature observables to derive the mass profiles at the end. The ‘non-parametric-like’ refers to the fact that we used a deprojection technique to derive the density and temperature profiles, i.e. without using parametric models.

The mass profiles of the Low-z PXI and High-z samples, and their associated uncertainties, are shown in Fig.2 (those of the five highest-redshift clusters are reproduced as published in Bartalucci et al. 2018). The profiles of the High-z clusters are mapped at least up to , with objects measured up to , thus not requiring extrapolation to compute the HE mass at a density contrast of , . The median statistical error is about . The quality of the Low-z PXI mass profiles, based on archival data, is lower and less homogeneous. Only clusters have temperature profiles extending to , with a maximum radius between and . Bartalucci et al. (2018) have shown that, while the HE mass is very robust when the temperature profiles extend up to , the mass is very sensitive to the mass estimation method when extrapolation is required. This is particularly the case for irregular clusters. To minimise systematic errors, we used rather than in the following for all scaling with mass, and for the computation of the sparsity (see below). The possible impact of this choice on our results is discussed in Sect. 7.

### 3.5 Sparsity

The sparsity, , was introduced by Balmès et al. (2014) to quantify the shape of the dark matter profile. It is defined as the ratio of the integrated mass at two over-densities, and has the advantage of being non-parametric, as there is no a priori assumption on the form of the profile. It therefore represents a useful measure when dealing with a population of objects with a wide variety of dynamical states. Another non-parametric approach, advocated by Klypin et al. (2016), considers the maximum circular velocity, which is linked to the traditional NFW or Einasto concentration. As the velocity is basically the square root of , it can also be derived from observations. In practice however, measuring such a maximum is much more difficult than measuring a ratio of integrated masses.

In the following, we thus concentrate on the sparsity to investigate the shape of the mass profiles. It is defined as:

 SΔ1,Δ2≡MΔ1MΔ2, (2)

where is the mass corresponding to the density contrast and with . We recall that , is the mass enclosed within , such that:

 M(

Balmès et al. (2014) argue that the general properties of the sparsity do not depend on the choice of as long as the halo is well defined (i.e. is not too small), and that the interaction between dark and baryonic matter in the central region can be neglected (i.e. is not too large). We use and , the latter choice further discussed in Sect. 6.2.

## 4 Morphological analysis

### 4.1 Centroid shift

We produced count images for each camera in the soft band, keV, binned using pixel size, on which we excised and refilled the masked regions where point sources were detected using the Chandra interactive analysis of observation (CIAO, Fruscione et al. 2006) dmfilth tool. Sub-structures were not masked for this analysis. We estimated the background following a similar approach to that of Böhringer et al. (2010). We computed the background map for each camera by fitting the refilled count images using a linear combination of the vignetted and unvignetted exposure maps to account for the instrumental and sky background, respectively. We removed the cluster emission by masking a circular region within and centred on the X-ray peak. Exposure maps, background and count images of MOS1, MOS2, and PN were combined, weighting by the ratio of integrated surface brightness profile of each individual camera to that of the combined profile. The combined count images were then background subtracted and exposure corrected. We produced 100 realisations of the count rate maps by applying the same procedure to Poisson realisations of the count maps.

The centroid shift parameter, , has been introduced by Mohr et al. (1993) as a proxy to characterise the dynamical state of a cluster. The centroid within an aperture is defined as

 (xc,yc)≡1Ni∑knk(xk;yk), (4)

where is the total number of counts per second within the ith aperture and is the count rate in the pixel of coordinates . We computed the mean deviation of the centroid from the X-ray peak by measuring the displacement within apertures using the definition of Böhringer et al. (2010) :

 ⟨w⟩=(1N−110∑i=0(Δi−⟨Δ⟩)2)1/21RYX500, (5)

where is the projected distance between the X-ray peak and the centroid computed within the i th concentric annulus, each one being wide. Uncertainties on the centroid shift were estimated by measuring  on the 100 Poisson count map realisations, and taking the values within of the median.

Maughan et al. (2008) measured  on a sample of clusters observed by Chandra, excluding the inner  kpc to make the parameter less sensitive to very bright cores. The PSF of XMM-Newton is larger for all the clusters considered in this work and, for this reason, we did not excise the core from the analysis. The good agreement between the  values derived at by Chandra and XMM-Newton  shown by Bartalucci et al. (2017), indicates that the XMM-Newton PSF is not an issue. Nurgaliev et al. (2013) have shown that the centroid shift can be biased high in the case of observations with a low number () of counts. In our sample the minimum number of counts in the  keV band we used to measure  is . Furthermore, all clusters are in the high SN regime, the lowest SN in our sample being .

Following Pratt et al. (2009), we initially classify an object as ‘morphologically disturbed’ if and ‘morphologically regular’ if . The results of the centroid shift characterisation for the Low-z and High-z samples are shown in the top and bottom left panels of Fig. 3, and the corresponding values are given in Table 2.

### 4.2 Surface brightness concentration CSB

The ratio of the surface brightness profile within two concentric apertures, hereafter the , was introduced by Santos et al. (2008) to quantify the concentration of cluster X-ray emission. We computed the  using the following definition:

 CSB=∫0.1×RYX5000SX(r)dr∫0.5×RYX5000SX(r)dr, (6)

where the error was computed using a Monte-Carlo procedure on 100 Gaussian realisations of the surface brightness profiles and taking the value around the median. The  parameter is a robust X-ray measurement, as it relies on the extraction of surface brightness profiles only and is not model-dependent. We nonetheless corrected for the XMM-Newton PSF, in view of the high and small angular size of the high sample.

Santos et al. (2010) demonstrated that, for objects at high redshift, the emission within two apertures requires a different -correction due to the presence of a cool core i.e. cool cores will have typically softer spectrum than the surrounding regions. This correction is potentially important for this work, as we are comparing  in a wide redshift range. Santos et al. (2010) proposed a correction for this effect which requires spatially resolved temperature profiles. At the median redshift of the Low-z sample, the correction is negligible (), but it can be up to at the highest redshifts. We thus did not make this correction for the Low-z sample. For the High-z sample, we applied the -correction to the  of all the clusters as if they were observed at the median redshift of the Low-z sample, z=0.19. Henceforth, all the  values shown and used in this work are -corrected in this manner. The results of the  analysis are reported in the top and bottom central panels of Fig. 3 and in Table 2.

The  allows the identification of cool-core (hereafter CC) clusters, the parameter being tightly correlated with the cooling time (e.g. Croston et al. 2008; Santos et al. 2010; Pascut & Ponman 2015). From the correlation between central density and cooling time in the REXCESS sample, Pratt et al. (2009) defined a central density of as a threshold which segregates CC and non-CC clusters (their Fig. 2). We computed the central density for the Low-z PXI and High-z objects and used the correlation between this quantity and the  to translate this density threshold in terms of the , finding that CC clusters have for this classification scheme.

The  value can also be used as an indicator of the relaxation state of the cluster. A high concentration is an indication that the core has not been disturbed by recent merger events. The corresponding threshold defined to distinguish relaxed clusters, e.g. from the anti-correlation observed between  and  and/or visual inspection (Cassano et al. 2010; Lovisari et al. 2017) differs from that used to define the CC/NCC segregation.

### 4.3 Combined dynamical indicator, M

The combination of certain morphological parameters has been shown to identify the most disturbed and relaxed clusters (Cassano et al. 2010; Rasia et al. 2013; Lovisari et al. 2017; Cialone et al. 2018). We take advantage of the observed anti-correlation between the centroid shift and the  to compute the parameter introduced by Rasia et al. (2013) and use it as an additional dynamical indicator. It is defined as

 (7)

where and are the median values of the  and centroid shift, respectively, and and the first or third quartile depending on whether the parameter value is larger or smaller than the median, respectively. is thus an indicator that combines the large scale (i.e. centroid shift) and the core (i.e concentration) properties. According to this definition, clusters which are characterised by the presence of a cool core and are morphologically regular will have and very disturbed objects with a very diffuse core will have . Henceforth, we refer to the former and latter objects as ‘relaxed’ and ‘disturbed’, respectively. The choice of this dual classification is arbitrary and does not correspond to a strict segregation between two types of objects. The distribution of both for Low-z and High-z samples is continuous, as shown in the top right panel of Fig. 3. The results of the characterisation are reported in the top and bottom right panels of Fig. 3 and in Table 2. Note that the numerical value of that we use is smaller than the threshold used in Sect. 4.2 to define CC clusters, and is closer to the value chosen by Lovisari et al. (2017), after visual classification of clusters, when defining a similar M parameter 555 which is close to , used by Lovisari et al. (2017), after correction for the different aperture definition (their Fig. C1) .

### 4.4 Consistency of the morphological characterisation

We used the results of the morphological analysis published in Lovisari et al. (2017) to characterise the Low-z cluster morphological properties, using their  and  values to derive the parameter. We obtained the morphological parameters of the High-z sample using a different pipeline and different analysis settings. To avoid potentially biased conclusions on morphological evolution or the dependence of the mass profiles on morphology for the full sample, it is necessary to check the consistency between our morphological analysis and that of Lovisari et al. (2017). We thus compared the morphological parameters derived from each pipeline independently for the common clusters of the Low-z PXI sample. The agreement is excellent, as shown in the left and central panel of Fig. 7. Full details of the comparison are discussed in Appendix A.

## 5 Morphology and dynamical state

### 5.1 Sample characterisation and comparison

The results of the Low-z and High-z morphological characterisation are shown in Fig. 3, where the top and bottom of each panel show the normalised and cumulative distribution for each parameter, respectively. We also show the fraction of objects above the fiducial thresholds discussed in Sect. 4.1 in each panel. Errors are computed from 1000 Monte-Carlo bootstrap resamples. They are dominated by the number of clusters, the individual uncertainties being much smaller than the intrinsic dispersion.

The top left panel shows that both samples contain a majority of disturbed () objects. The shape of the distributions differs, the High-z sample having a prominent peak of objects around . This is reflected in the cumulative distribution, where there is an excess of High-z objects at . However, the fraction of disturbed objects in the Low-z sample, , and in the High-z sample, , are nearly identical and are consistent within the uncertainties. We investigated if the two samples are representative of the same population by performing the Kolmogorov-Smirnov (KS) test. We determine the p-value of the null hypothesis, i.e. that the two samples are drawn from the same distribution. The -value is , indicating no significant difference between the Low-z and High-z samples.

We obtained similar results studying the distribution of CC clusters using the  parameter, as shown in the central panels of Fig. 3. The fraction of CC clusters is low. The High-z sample has a slightly lower fraction of CC objects, , as compared to the Low-z sample, , but the difference is not significant. Consistently, the KS test yields a high -value of .

Interestingly, the distributions shown in the right panels of Fig. 3 suggest some evolution. While the distributions have qualitatively similar shapes, the High-z sample has a peak which is clearly shifted towards disturbed objects. Furthermore, the fraction of relaxed clusters in the Low-z sample, , is higher than that of the High-z sample, , a effect. This evolution can be seen in the cumulative distribution as a systematic over-abundance of disturbed objects in the High-z sample as compared to the Low-z PXI. The KS test yields a smaller -value of , however not small enough to reject the null-hypothesis.

### 5.2 Mass dependence and redshift evolution

For each morphological parameter, we further quantified the relation with mass and redshift by computing the Spearman’s rank (SR) coefficient and the null hypothesis -value, the probability that the observed coefficient is obtained by chance if the two parameters are completely independent. We also considered the sum square difference of ranks, , and the number of standard deviations by which deviates from its null-hypothesis expected value, . As in the previous section, we performed bootstrap resamples to estimate these values and their errors. The top and bottom panels of Fig. 4 show each parameter as a function of redshift and mass, respectively. The centroid shift,  and parameters are shown on the left, middle and right panels respectively, with the SR coefficient and corresponding -value indicated in the top left of each plot.

The only parameter for which there is a correlation with is the combined parameter, for which . The correlation is not very significant, with a null hypothesis -value of , and a standard deviation on the null hypothesis of . The Kendall test gives consistent results. This weak correlation of with comes from the amplification of the positive (but not significant) trend in  versus , while there is no correlation between  and redshift.

In summary, consistent with the trend observed in Sect. 5.1 above, there is weak evidence for clusters at higher redshift to be slightly more disturbed. On the other hand, there is no evidence for any trend with mass, the SR coefficient for all parameters being consistent with zero and the corresponding -value in the range . This is in agreement with the mass independence of the dynamical state, found by Böhringer et al. (2010) and Lovisari et al. (2017) for local clusters. We must note, however, that the mass range is limited in the present sample.

## 6 Total mass profile shape

The individual scaled HE integrated total mass profiles are shown in Fig. 5. In the left and central panels, the profiles are colour coded according to their morphological state according to the parameter, in blue for the relaxed () and in red for the disturbed () clusters.

There is a clear difference between the two populations in the range. As compared to the relaxed clusters, the disturbed objects have a shallower mass distribution on average (i.e. there is less mass within the central region), and these profiles show a larger dispersion. This effect is even more evident in the central panel of Fig. 5, where the mass profiles are plotted as a function of the total density contrast, . As is proportional to the mean total density within a sphere of radius (Eq. 3), it decreases with radius, more or less rapidly depending on the steepness of the density profile. For density profiles that are very peaked towards the centre, decreases rapidly, or equivalently slowly increases with decreasing . On the other hand, a flat profile within the core would correspond to a constant mean density, thus quasi–constant , and a very steep variation of with in the central region, with a maximum value of corresponding to that of the core. This is what is observed for the relaxed and disturbed clusters, respectively. In summary, the dynamical state of the cluster clearly has a very strong impact on the shape of the total mass profile.

The right panel of Fig. 5 shows the dispersion envelopes of the Low-z PXI and High-z samples with light and dark green, respectively. The black solid and dotted lines represent the median profiles for the Low-z PXI and High-z samples, respectively. The envelopes of the two samples are consistent. However, the median HE mass profile of the High-z sample is slightly lower and shallower than that of the Low-z PXI sample, being lower by and at and 0.3, respectively. The HE radial mass profiles thus show a hint of evolution, but this is not statistically significant considering the dispersion of the profiles. Interestingly, the dispersions of the two samples are similar. We recall that the samples have similar number of objects, the Low-z PXI and High-z having 42 and 33 objects, and the data quality ensures that HE mass profiles are computed for each object. In view of our finding above that the morphology evolution is negligible, the absence of evolution of the HE mass profile dispersion is a natural consequence of their shape being driven by the dynamical state.

### 6.2 The shape of the mass profiles

We further quantified the evolution and the impact of dynamical state on the shape of the profile using the sparsity, . We chose , which is large enough to encompass the dependence of the mass profile shape as a function of the dynamical state and is reached by all haloes, as shown in the middle panel of Fig. 5.

Figure 6 shows the sparsity as a function of the redshift in the top-left panel, and the three morphological parameters, ,  and in the top-right, bottom-left and bottom-right panels, respectively. The statistical errors on the sparsity are not negligible as compared to the intrinsic scatter, and the simple correlation tests used in the previous section cannot be applied. We thus computed the mean sparsity in bins. We defined three bins in redshift and for each parameter, the width of each bin being defined to have roughly the same number of objects in each bin. The logarithmic–mean of the sparsity was computed in each bin, weighting each value by the quadratic sum of the statistical error and the intrinsic scatter. This scatter and the weighted mean were estimated simultaneously by iteration. The mean and intrinsic scatter, together with their errors, were computed using 1000 bootstrap resamples. The results are reported in Table 1 and as filled blue rectangles in each panel of Fig 6. As there is a large scatter with the presence of strong outliers, we also computed the mean within the same bins excluding the outliers. The results are shown with the open blue rectangles.

There seems to be a slight, although not very significant, evolution of the sparsity with . Higher redshift () clusters have a slightly larger value of  by , i.e. these profiles are less concentrated, which is not consistent with the low redshift () clusters at . We found the same behaviour when excluding the outliers.

There is a much stronger, and significant, variation of the sparsity with dynamical state: morphologically disturbed (high ), non-CC (low ) and disturbed (low ) clusters have larger sparsity. For all parameters considered, the sparsity of the first and third bins are not consistent at more than and their difference is of the order of . At the same time, the intrinsic scatter increases significantly, reaching dex for the most disturbed or least concentrated objects. For example, the sparsity increases by at a significance level of between and , i.e. between the most relaxed and the most disturbed as defined by this parameter. Only an upper limit on the intrinsic scatter can be estimated from the former ( dex), while the intrinsic scatter for the latter reaches dex. There are strong outliers at high sparsity. Excluding the outliers yields the same qualitative results, although the variation between bins is weaker. This suggests a sparsity distribution skewed towards high values, with a skewness increasing with departure from dynamical relaxation.

## 7 Discussion

### 7.1 Dynamical state

The first result of our study is that SZ-selected samples are dominated at all redshifts by disturbed and non-CC objects. Recent observational work on local clusters has converged to similar results (e.g. Lovisari et al. 2017, Rossetti et al. 2017, Andrade-Santos et al. 2017, and Lopes et al. 2018). In particular, these works highlight the higher fraction of disturbed objects or the lower fraction of CC objects in SZ-selected samples as compared to X–ray selected samples, a fact interpreted to be due to preferential detection of relaxed, or more concentrated, clusters in X–ray surveys.

In contrast, there is little consensus in the literature concerning the evolution of the dynamical state, as determined from various morphological parameters. Studies of the evolution up to of the centroid shift and/or power ratios of X–ray selected clusters indicate a larger fraction of disturbed clusters at high (Maughan et al. 2008; Jeltema et al. 2005; Weißmann et al. 2013). However, the latter study is also consistent with no evolution, and Nurgaliev et al. (2017) did not find any significant evolution in the redshift range for the 400d X–ray selected clusters.

These somewhat contradictory results may simply be due to selection effects: X-ray detectability is clearly not independent of the cluster morphology. More peaked clusters (usually relaxed) are not only more luminous at given mass, but also easier to detect at given X–ray luminosity. Such effects are particularly important in flux-limited surveys, as shown by Chon & Böhringer (2017) in the context of a volume limited X–ray survey. It is thus difficult to disentangle selection effects and/or and mass dependence (see also the discussion in Mantz et al. 2015).

SZ detection does not suffer from these limitations, and an SZ-selected sample is expected to be close to mass-limited. The morphological evolution of SPT clusters was studied recently by Nurgaliev et al. (2017), using their newly-introduced parameter. They did not find a significant difference between the redshift ranges and . This absence of significant evolution was also observed by McDonald et al. (2017) up to , also for an SPT sample. The comprehensive study of classification criteria for the most relaxed clusters by Mantz et al. (2015) indicates that the fraction of relaxed clusters in the SPT and Planck samples is consistent with being constant with redshift.

In the present study, we extend the morphological analysis of the full SZ-selected population of high mass clusters, from very local systems up to , applying a consistent sample construction and analysis strategy over the full range. The high quality of our data allows us to investigate core (i.e., ) and bulk (i.e., ) properties at high precision. There are no significant trends either with or mass in these parameters individually. However, we find a significant evolution with of the dynamical indicator, which combines these large scale and core parameters, with a null hypothesis -value of .

It has been asserted that, naively, the fraction of disturbed clusters should increase with and mass in a hierarchical formation model (e.g. Böhringer et al. 2010; McDonald et al. 2017). To explain the observed absence of evolution in the SPT sample, McDonald et al. (2017) proposed a simple model, combining the merger rate from the simulations of Fakhouri et al. (2010) and a fiducial relaxation time of the hot gas equal to the crossing time. In fact, the link between cluster formation history and morphological state as observed in X-ray, as a function of and mass, is very complex. This first depends on relating the individual mass assembly history to dynamical state (Power et al. 2012, e.g.), and then the dynamical state to morphological indicators (e.g. Cui et al. 2017). Individual cluster history is never observed directly and has to be translated into to the ensemble properties of cluster samples at different (e.g. see Mostoghiu et al. 2019). A further complication is the relation between the gas dynamical history and that of the underlying dark matter, and how X–ray morphological observables relate to the gas dynamical state. To our knowledge the only theoretical prediction of the evolution of observed ICM morphological parameters is that of Jeltema et al. (2008). They claim a significant evolution of the mean  with , although a comparison of their results with our data in Fig. 4 shows that any such evolution is mild, is much smaller than the dispersion, and is fully in agreement with our results.

### 7.2 Total mass profile

Extending the pilot study of Bartalucci et al. (2018) with a fully SZ-selected sample, the overall picture emerging from the present work is that the shape of the dark matter profiles is affected both by evolution and by the dynamical state (Table 1). The evolution effect is mild, increasing the sparsity of objects by from to . In contrast, the dependence is much stronger, with a sparsity increasing by with decreasing , i.e. from the most relaxed to the most disturbed objects. This variation of profile shape with dynamical state is likely a fundamental property, rather than a secondary consequence of the mutual variation of the  and parameters with , which are both less significant. A multi-component analysis, requiring a larger sample, is needed to firmly assess this point.

An obvious question is whether the observed dependence of the sparsity on dynamical state is an artefact due to systematics in the X-ray mass estimate. As discussed in detail by Corasaniti et al. (2018) the sparsity derived from HE mass estimates is essentially bias-free (less than ). As it is a mass ratio, the sparsity is only sensitive to the radial dependence of any bias, which is usually small between the density contrasts under consideration. For instance the mean bias from HE mass estimate is of from the simulations of Biffi et al. (2016). Generally, although the exact radial dependence of the HE bias will differ from object to object, we expect it to increase towards larger radii (smaller ) so that sparsity  measured from HE profiles would be biased low as compared to the true value, especially for the most disturbed objects. This effect is the opposite to the observed increase of sparsity for increasingly perturbed systems. As detailed in Sect. 3.5 we are not directly using the HE mass at but its proxy, . For the clusters for which no extrapolation is needed, the differences between these measures are the order of , with no systematic trend with dynamical state. Thus the measured sparsity may be slightly lower than if we had used the HE mass, but the clear trend of the sparsity with dynamical state would not be changed.

The trend of sparsity with dynamical state indicates that morphologically disturbed objects are less concentrated than relaxed objects. There is also evidence for increased scatter. This is qualitatively in agreement with the difference in the relations for relaxed versus disturbed objects seen in numerical simulations (e.g. Bhattacharya et al. 2013, Fig. 1 and 4). In Le Brun et al. (2018) we performed a first investigation using numerical simulations tailored to cover the high–mass high–z range considered here. In these simulations, the sparsity of the 25 most massive clusters at all redshifts shows a correlation with the DM dynamical indicator r with a -value of (their Fig. 3). There is also a very significant increase in the dispersion with r. We will revisit the link between dynamical state and DM sparsity for the full simulated sample in a forthcoming paper.

## 8 Conclusions

We presented new XMM-Newton observations of a Planck SZ selected sample of 28 massive clusters in the redshift range . These were combined with the sample of Bartalucci et al. (2018) at , for a total High-z sample of 33 objects at masses . We characterised the dynamical state with the centroid shift , the concentration , and the combination of the two parameters, , which simultaneously probes the large scale and core morphology. The shape of the total mass profile, derived from the hydrostatic equilibrium equation, was quantified using the sparsity, the ratio of to , the masses at density contrast 500 and 2500, respectively. This parameter, introduced by Balmès et al. (2014) offers a non-parametric measure of the shape which is thought to be relatively insensitive to HE bias (Corasaniti et al. 2018) .

We first combined the High-z morphology measurements with that of the ESZ clusters at in Lovisari et al. (2017), for a total sample of 151 objects.

• We confirmed that SZ-selected samples, thought to best reflect the underlying cluster population, are dominated by disturbed () and non-CC () objects, at all redshift.

• There is no significant evolution or mass dependence of the fraction of cool core or of the centroid shift parameter. The only parameter for which there is a significant correlation with is the combined parameter, for which and a null hypothesis -value of .

We then combined the mass measurements obtained for our new data with those from a sub-sample of 42 ESZ objects, with spatially resolved ICM profiles presented in Planck Collaboration V (2013), and which we confirmed is representative of the full Low-z sample. The total sample of 75 objects covers the redshift range and mass range . We found that :

• The median scaled mass profile differs by less than and at and 0.3, respectively, between the Low-z PXI and High-z samples, with no difference in the dispersion. The evolution of the sparsity with is mild: it increases by only between and , an effect significant at less than .

• When expressed in term of scaled mass profiles, there is a clear difference between relaxed and disturbed objects. The latter have a less concentrated mass distribution, on average, and their scaled profiles show a much larger dispersion.

• Consequently there is a clear dependence of the sparsity on the dynamical state. When expressed as a function of the parameter, the sparsity increases by a from the 1/3 most relaxed to the most disturbed objects, an effect significant at more than level. We discussed that the HE bias will not change significantly this result.

The main result of this work is that the radial mass distribution is chiefly governed by the dynamical state of the cluster and only mildly dependent on redshift. This has important consequences:

• A coherent sample selection at all is key. For instance, one cannot compare the relation calibrated at low on the most relaxed X-ray selected clusters, to that of a SZ selected sample at high . Ideally, one should consider complete samples, representative of the full true underlying cluster population i.e. mass selected sample. This is even more critical when comparing theory and observation, in view of the difficulty of defining coherent dynamical indicators between the two.

• To test theoretical predictions, it is insufficient to simply compare median or stacked properties at each . The dispersion is a critical quantity, as is the profile distribution, in view of likely departure from log-normality. This requires the measurement of individual profiles at high statistical precision. Presently, X–ray observations currently provide the best way to obtain the required statistical precision.

• In view of this, observational and theoretical efforts to understand the HE bias and its radial dependence are all the more important.

In a forthcoming paper, we will extend our study of the dependence of the sparsity on the dynamical state using an extension of the dark matter simulations presented in Le Brun et al. (2018) to a larger sample of clusters at . On a longer time scale, the link between the true dark matter distribution and its dynamical state and the X-ray observables will need to be better understood. This includes the link between ICM morphological proxies and the true underlying dynamical state, and the potential critical issue of the radial variation of the HE bias. We will address these issues with dedicated simulations.

On the observational side, we will investigate the HE bias by comparing our results to weak lensing mass measurements for a sub-sample of the present data set. The observed lack of significant evolution needs to be tested with a larger sample, particularly at , where the present sample is limited, with data of the same or better quality. The fundamental link uncovered by the present paper between the mass profile and dynamical state will also need to be consolidated with better low data. Our current low–z sample relies on archival data, of uneven quality, and is not a complete sample. The building of a new local reference SZ-selected sample, with high quality ICM thermodynamic and HE mass profiles, will be one of the main outcomes of the AO17 XMM-Newton heritage program ‘Witnessing the culmination of structure formation in the Universe’, and will provide the necessary inputs.

###### Acknowledgements.
The results reported in this article are based on data obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This work was supported by CNES. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP72007-2013) ERC grant agreement no 340519. LL acknowledges support from NASA through contracts 80NSSCK0582 and 80NSSC19K0116.

## References

• Amodeo et al. (2016) Amodeo, S., Ettori, S., Capasso, R., & Sereno, M. 2016, A&A, 590, A126
• Andrade-Santos et al. (2017) Andrade-Santos, F., Jones, C., Forman, W. R., et al. 2017, ApJ, 843, 76
• Arnaud et al. (2001) Arnaud, M., Neumann, D. M., Aghanim, N., et al. 2001, A&A, 365, L80
• Arnaud et al. (2010) Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92
• Balmès et al. (2014) Balmès, I., Rasera, Y., Corasaniti, P.-S., & Alimi, J.-M. 2014, MNRAS, 437, 2328
• Bartalucci et al. (2017) Bartalucci, I., Arnaud, M., Pratt, G. W., et al. 2017, A&A, 598, A61
• Bartalucci et al. (2018) Bartalucci, I., Arnaud, M., Pratt, G. W., & Le Brun, A. M. C. 2018, A&A, 617, A64
• Bhattacharya et al. (2013) Bhattacharya, S., Habib, S., Heitmann, K., & Vikhlinin, A. 2013, ApJ, 766, 32
• Biffi et al. (2016) Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112
• Bleem et al. (2015) Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27
• Böhringer et al. (2010) Böhringer, H., Pratt, G. W., Arnaud, M., et al. 2010, A&A, 514, A32
• Cassano et al. (2010) Cassano, R., Ettori, S., Giacintucci, S., et al. 2010, ApJ, 721, L82
• Chon & Böhringer (2017) Chon, G. & Böhringer, H. 2017, A&A, 606, L4
• Cialone et al. (2018) Cialone, G., De Petris, M., Sembolini, F., et al. 2018, MNRAS, 477, 139
• Corasaniti et al. (2018) Corasaniti, P. S., Ettori, S., Rasera, Y., et al. 2018, ApJ, 862, 40
• Correa et al. (2015) Correa, C. A., Wyithe, J. S. B., Schaye, J., & Duffy, A. R. 2015, MNRAS, 452, 1217
• Croston et al. (2006) Croston, J. H., Arnaud, M., Pointecouteau, E., & Pratt, G. W. 2006, A&A, 459, 1007
• Croston et al. (2008) Croston, J. H., Pratt, G. W., Böhringer, H., et al. 2008, A&A, 487, 431
• Cui et al. (2017) Cui, W., Power, C., Borgani, S., et al. 2017, MNRAS, 464, 2502
• da Silva et al. (2004) da Silva, A. C., Kay, S. T., Liddle, A. R., & Thomas, P. A. 2004, MNRAS, 348, 1401
• Diemer & Kravtsov (2015) Diemer, B. & Kravtsov, A. V. 2015, ApJ, 799, 108
• Dolag et al. (2004) Dolag, K., Bartelmann, M., Perrotta, F., et al. 2004, A&A, 416, 853
• Dutton & Macciò (2014) Dutton, A. A. & Macciò, A. V. 2014, MNRAS, 441, 3359
• Einasto (1965) Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87
• Fakhouri et al. (2010) Fakhouri, O., Ma, C.-P., & Boylan-Kolchin, M. 2010, MNRAS, 406, 2267
• Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Proc. SPIE, Vol. 6270, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 62701V
• Hasselfield et al. (2013) Hasselfield, M., Hilton, M., Marriage, T. A., et al. 2013, J. Cosmology Astropart. Phys., 7, 008
• Jeltema et al. (2005) Jeltema, T. E., Canizares, C. R., Bautz, M. W., & Buote, D. A. 2005, ApJ, 624, 606
• Jeltema et al. (2008) Jeltema, T. E., Hallman, E. J., Burns, J. O., & Motl, P. M. 2008, ApJ, 681, 167
• Jing (2000) Jing, Y. P. 2000, ApJ, 535, 30
• Kalberla et al. (2005) Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775
• Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340
• Kravtsov & Borgani (2012) Kravtsov, A. V. & Borgani, S. 2012, ARA&A, 50, 353
• Kravtsov et al. (2006) Kravtsov, A. V., Vikhlinin, A., & Nagai, D. 2006, ApJ, 650, 128
• Kuntz & Snowden (2000) Kuntz, K. D. & Snowden, S. L. 2000, ApJ, 543, 195
• Le Brun et al. (2018) Le Brun, A. M. C., Arnaud, M., Pratt, G. W., & Teyssier, R. 2018, MNRAS, 473, L69
• Lopes et al. (2018) Lopes, P. A. A., Trevisan, M., Laganá, T. F., et al. 2018, MNRAS, 478, 5473
• Lovisari et al. (2017) Lovisari, L., Forman, W. R., Jones, C., et al. 2017, ApJ, 846, 51
• Ludlow et al. (2014) Ludlow, A. D., Navarro, J. F., Angulo, R. E., et al. 2014, MNRAS, 441, 378
• Lumb et al. (2002) Lumb, D. H., Warwick, R. S., Page, M., & De Luca, A. 2002, A&A, 389, 93
• Mantz et al. (2015) Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2015, MNRAS, 449, 199
• Marriage et al. (2011) Marriage, T. A., Acquaviva, V., Ade, P. A. R., et al. 2011, ApJ, 737, 61
• Maughan et al. (2008) Maughan, B. J., Jones, C., Forman, W., & Van Speybroeck, L. 2008, ApJS, 174, 117
• McDonald et al. (2017) McDonald, M., Allen, S. W., Bayliss, M., et al. 2017, ApJ, 843, 28
• Mohr et al. (1993) Mohr, J. J., Fabricant, D. G., & Geller, M. J. 1993, ApJ, 413, 492
• Morrison & McCammon (1983) Morrison, R. & McCammon, D. 1983, ApJ, 270, 119
• Mostoghiu et al. (2019) Mostoghiu, R., Knebe, A., Cui, W., et al. 2019, MNRAS, 483, 3390
• Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
• Navarro et al. (2004) Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039
• Neto et al. (2007) Neto, A. F., Gao, L., Bett, P., et al. 2007, MNRAS, 381, 1450
• Nurgaliev et al. (2017) Nurgaliev, D., McDonald, M., Benson, B. A., et al. 2017, ApJ, 841, 5
• Nurgaliev et al. (2013) Nurgaliev, D., McDonald, M., Benson, B. A., et al. 2013, ApJ, 779, 112
• Pascut & Ponman (2015) Pascut, A. & Ponman, T. J. 2015, MNRAS, 447, 3723
• Planck Collaboration V (2013) Planck Collaboration V. 2013, A&A, 550, A131
• Planck Collaboration VIII (2011) Planck Collaboration VIII. 2011, A&A, 536, A8
• Planck Collaboration XI (2011) Planck Collaboration XI. 2011, A&A, 536, A11
• Planck Collaboration XXIX (2014) Planck Collaboration XXIX. 2014, A&A, 571, A29
• Planck Collaboration XXVII (2016) Planck Collaboration XXVII. 2016, A&A, 594, A27
• Power et al. (2012) Power, C., Knebe, A., & Knollmann, S. R. 2012, MNRAS, 419, 1576
• Pratt et al. (2019) Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25
• Pratt et al. (2010) Pratt, G. W., Arnaud, M., Piffaretti, R., et al. 2010, A&A, 511, A85
• Pratt et al. (2009) Pratt, G. W., Croston, J. H., Arnaud, M., & Böhringer, H. 2009, A&A, 498, 361
• Rasia et al. (2013) Rasia, E., Meneghetti, M., & Ettori, S. 2013, The Astronomical Review, 8, 40
• Reichardt et al. (2013) Reichardt, C. L., Stalder, B., Bleem, L. E., et al. 2013, ApJ, 763, 127
• Rossetti et al. (2017) Rossetti, M., Gastaldello, F., Eckert, D., et al. 2017, MNRAS, 468, 1917
• Santos et al. (2008) Santos, J. S., Rosati, P., Tozzi, P., et al. 2008, A&A, 483, 35
• Santos et al. (2010) Santos, J. S., Tozzi, P., Rosati, P., & Böhringer, H. 2010, A&A, 521, A64
• Starck et al. (1998) Starck, J.-L., Murtagh, F., & Bijaoui, A. 1998, Image Processing and Data Analysis: The Multiscale Approach (New York, NY, USA: Cambridge University Press)
• Strüder et al. (2001) Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18
• Turner et al. (2001) Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27
• Velliscig et al. (2014) Velliscig, M., van Daalen, M. P., Schaye, J., et al. 2014, MNRAS, 442, 2641
• Weißmann et al. (2013) Weißmann, A., Böhringer, H., & Chon, G. 2013, A&A, 555, A147
• Wu et al. (2013) Wu, H.-Y., Hahn, O., Wechsler, R. H., Mao, Y.-Y., & Behroozi, P. S. 2013, ApJ, 763, 70

## Appendix A Low-z vs Low-z PXI characterisation

In this work we used the results of Lovisari et al. (2017) to characterise the morphological properties of the Low-z sample, using it as anchor for the local universe properties. For this reason, it is mandatory that the morphological parameters we derived for the High-z sample are coherent with the values computed from Lovisari et al. (2017). We derived the centroid shift, , and the for the Low-z PXI objects we have in common with Lovisari et al. (2017). The comparison between these values are shown in Fig. 7, denoting our and Lovisari et al. (2017) values with ”This work” and ”L+17” labels, respectively. The measurements of the centroid shifts shown in the left panel are in good agreement, with a ratio of and a standard deviation of . The strongest outlier is MACSJ2243.3-0935, shown with a red point. The difference is probably caused by the different choice of the X-ray peak (distant by ) amplified by the high ellipticity of this object. The concentration parameter  and the masses, shown in the central and right panels, respectively, are in excellent agreement. The median ratio for both quantities is excellent and the dispersion is remarkably small, the two analysis being performed with different pipelines. This comparison shows that the measurements which require more statistic, as the centroid shift, are more sensitive to analysis parameters such as the choice of the centre, the exclusion of point sources. Integrated quantities as the  are more robust.

## Appendix B Representativity of the Low-z PXI

We derived the spatially resolved HE mass profiles for a sub-sample of 42 clusters of the Low-z sample, namely the Low-z PXI. This sub-sample covers the same Low-z redshift range and comprises clusters with . We investigated if this sub-sample is representative of the Low-z population in terms of morphological status.

The Low-z and the Low-z PXI  distributions are shown in the left panels of Fig. 8 with green and blue polygons, respectively. In particular, the normalised and cumulative distributions are shown in the left and bottom panels, respectively. The two samples present qualitatively the same distribution and similar fraction of disturbed cluster, the Low-z and the Low-z PXI yielding a fraction of disturbed clusters equal to and to , respectively. This results is confirmed by the high value of the KS test probability, . We found similar results using the . The cluster distributions as a function of this parameter are shown in the right top and bottom panels of Fig. 8. The two samples have the same fraction of CC objects, , and the KS test probability value, , confirms that they are representative of the same population.