Subscriber Authentication Point
Free Access
 Issue A&A Volume 634, February 2020 A113 21 Cosmology (including clusters of galaxies) https://doi.org/10.1051/0004-6361/201936742 19 February 2020

## 1. Introduction

Clusters of galaxies are the endpoint of the process of cosmic structure formation. As such, they are optimal tracers of the growth of structures and useful tools for estimating cosmological parameters, such as those measuring the amount of matter, ΩM, and dark energy, ΩΛ, and the normalization of the power spectrum of density fluctuations, σ8 (see reviews by Voit 2005; Allen et al. 2011; Kravtsov & Borgani 2012). Cosmological studies based on galaxy clusters rely on the measurement of (i) the baryon fraction (Allen et al. 2008; Ettori et al. 2009; Mantz et al. 2014), and (ii) the evolution of the cluster mass function, or the number density of clusters per unit mass and redshift interval (Vikhlinin et al. 2009; Planck Collaboration XX 2014; Costanzi et al. 2019; Bocquet et al. 2019). The key quantity entering into both these techniques is the cluster mass. It is therefore crucial to refine the estimate of the total gravitating mass in clusters of galaxies as precisely as possible by constraining in great detail the sources of bias or by providing reliable corrections (see review by Pratt et al. 2019).

Indeed, when the Planck collaboration pointed out that σ8 derived from the cosmic microwave background anisotropies was inconsistent with the σ8 obtained from cluster number counts (Planck Collaboration XX 2014), the so-called mass bias, defined as an underestimation of the total mass, was immediately suspected to be the root of the problem. Even though in the Planck-specific comparison other systematic errors might have played a role in increasing the discrepancy between the two σ8 values, there is recent observational evidence that the masses estimated under the assumption of hydrostatic equilibrium (HE), such as those derived from the X-ray analysis used in the Planck analysis, are actually underestimated with respect to masses measured using other methods (see reviews by Ettori et al. 2013; Pratt et al. 2019, and references therein).

In this context, cosmological hydrodynamical simulations have been insightful tools for understanding the nature and origin of HE mass bias (Evrard 1990). Advanced numerical models providing realistic populations of simulated clusters allow the precise quantification of the underestimation of the total mass. Such bias can be connected to the intrinsic properties of the simulated objects, such as their dynamical state, or to other quantities mimicking those derived from X-ray, optical, or millimeter observations. Since the earliest simulation studies, numerical resolution has largely increased, and progressively refined descriptions of star formation, feedback in metals, and energy from stars and from active galactic nuclei (AGNs) have vastly improved the level of realism and reliability of such simulations. At the same time, the techniques used to analyze the simulations and compare them to observational data, including mock images, have further enhanced their predictive power. Still, the level of the above mass bias for simulated clusters considered “dynamically relaxed” has always been consistently found to be around 10–20% (Rasia et al. 2006, 2012; Nagai et al. 2007a; Jeltema et al. 2008; Piffaretti & Valdarnini 2008; Lau et al. 2009; Meneghetti et al. 2010; Nelson et al. 2012; Battaglia et al. 2012; Shi et al. 2015; Biffi et al. 2016; Vazza et al. 2016, 2018; Henson et al. 2017; Barnes et al. 2017a; Cialone et al. 2018; Angelinelli et al. 2019). Thanks to simulations, we understand that the main sources of the HE bias are the residual, nonthermalized gas velocities in the form of both bulk motion and turbulence, and intra-cluster medium (ICM) inhomogeneities.

The existence of gas velocities compromises the assumption of HE: the gravitational force is not completely in equilibrium with the hydrodynamical pressure force as some nonthermal pressure support in the form of residual gas motions is still present in the gas (Rasia et al. 2004; Lau et al. 2009; Vazza et al. 2009; Fang et al. 2009; Suto et al. 2013; Shi et al. 2015; Biffi et al. 2016). Residual gas motions are expected because clusters are recently assembled systems and lie at the intersection of cosmic filaments that define the preferential directions along which there is continuous mass accretion in the form of both diffuse gas and over-dense clumps (Vazza et al. 2013; Zinger et al. 2018). Moreover, indirect hints of the presence of residual kinetic energy associated with random gas motions in the ICM comes from the observational evidence of diffuse radio emission connected with clusters undergoing mergers (and likely powered by the dissipation of turbulent motions via Fermi-like mechanisms; see van Weeren et al. 2019 for a recent review) and from the observed correlation between X-ray surface brightness fluctuations and radio power Eckert et al. 2017) which suggests a dynamical link between perturbed X-ray morphologies and residual gas motions. Unfortunately, for the direct ICM velocity measurements we need to wait for next-generation X-ray spectrometers (Biffi et al. 2013; Roncarelli et al. 2018; ZuHone et al. 2018; Simionescu et al. 2019; Cucchetti et al. 2019; Clerc et al. 2019) which might obtain such measurements in the external regions of (a few) clusters, such as those obtained by Hitomi for the core of the Perseus cluster (Hitomi Collaboration 2016).

In addition, in the presence of unresolved small, cold, and dense clumps, the gas density can be boosted towards higher values with respect to the smoother ICM distribution (Nagai & Lau 2011; Vazza et al. 2013; Zhuravleva et al. 2013; Roncarelli et al. 2013; Planelles et al. 2017; Walker et al. 2019). If the clumps are not in pressure equilibrium with the ambient medium, the pressure signal as derived from the Compton parameter measured by SZ is either boosted or hampered depending on whether the structures are at a hyperbaric or a hypobaric level (Battaglia et al. 2015; Khatri & Gaspari 2016; Planelles et al. 2017; Ruppin et al. 2018). Since clumps are more prominent in the outskirts of clusters, they lead to an apparent decrease in the slope of the gas profiles and, as a consequence, to an underestimate of the derived HE mass (which is proportional to the derivative of the gas density or pressure profile). Moreover, inhomogeneities in the temperature structure might cause an additional bias if the multi-temperature nature of the ICM is not recognized from the residuals of the spectroscopic fitting analysis. Indeed, X-ray CCD detectors onboard Chandra, XMM-Newton, and Suzaku have responses that tend to emphasize colder gas components in a thermally complex medium (Gardini et al. 2004; Mazzotta et al. 2004; Vikhlinin et al. 2006). A bias in the temperature measurement at a fixed radius automatically translates into a mass bias at the same radius. Simulations show that thermal inhomogeneities in the ICM increase with radius (Rasia et al. 2014), and so does the corresponding mass bias. However, measuring the temperature fluctuations from observations is extremely challenging (Frank et al. 2013). At the same time, direct measurement of clumpiness from observations requires a spatial resolution allowing to discern appropriately the X-ray or SZ signal on kiloparsec scales associated with exquisite spectral-imaging capabilities. These two requirements are not met by existing SZ telescopes, while with current X-ray instruments this level of clumpiness can be at least indirectly estimated (Walker et al. 2012, 2019; Morandi et al. 2013, 2017; Urban et al. 2014; Morandi & Cui 2014; Eckert et al. 2015; Ghirardini et al. 2018) and positively compares with results from simulations (Eckert et al. 2015).

Since quantifying the level of gas inhomogeneities is feasible with current X-ray instruments, in this paper we investigate how these estimates can be used to statistically correct the mass bias at R500. We also investigate an approach proposed almost a decade ago (e.g., Ameglio et al. 2009), and recently adopted in observational samples (e.g., Ettori et al. 2019), to directly include the pressure, which can be derived from SZ observations out to large radii, within the HE mass equation under the assumption that clump pressure is close to the isobaric level of the ICM. These goals are pursued with a detailed investigation of the simulated clusters of the Three Hundred Project (Cui et al. 2018), analyzed at z  =  0. The sample employed includes a large number of massive clusters (more than three hundred), comparable only to the MACSIS project (Barnes et al. 2017b). The large sample size allows us to adopt a conservative approach and discard all objects that have significant interactions and/or have an extremely complex morphology. The methods employed to extract the quantities of interest are simplified with respect to earlier work where mock Chandra or XMM-Newton event files were produced (Rasia et al. 2006, 2012; Nagai et al. 2007b; Meneghetti et al. 2010). In fact, previous analyses already demonstrated that gas density and temperature profiles can be recovered relatively accurately from mock event files (a more detailed discussion of this will be presented in Sect. 4). In this paper we simply build X-ray surface-brightness maps that are not convolved with any instrument response. We remark that the quantities considered have already been derived in X-ray studies from observations with a modest exposure time.

This paper is structured as follows: we present the simulations in Sect. 2. The subsamples based on the morphological classification are introduced in Sect. 3, and in Sect. 4 we describe in detail how we derive all quantities from the simulated clusters, through either the 2D analysis of the X-ray maps or the 3D intrinsic analysis. We discuss the results on the gas inhomogeneities from the maps and cluster clumpiness factor in Sects. 5 and 6, where we also highlight their mutual relation. The mass bias profiles and mass bias distribution at R500 are investigated in Sect. 7, where we also consider the dependence of the mass bias on all quantities linked to the gas inhomogeneities and on the asymptotic external slope of the gas density and pressure profiles. Finally, conclusions are drawn in Sect. 8.

## 2. Simulations

The hydrodynamical simulated clusters used in this work are part of the Three Hundred Project1 introduced in Cui et al. (2018) and analyzed in Wang et al. (2018), Mostoghiu et al. (2019), and Arthur et al. (2019) to respectively study galaxy properties in a rich environment, the evolution of the density profile, and the effect of ram pressure stripping on the gas content of halos and substructures.

These simulations are based on a set of 324 Lagrangian regions centered on as many galaxy clusters, which have been previously selected as the most massive within the parent MultiDark simulation (Klypin et al. 2016)2 and precisely the MultiDarkPlanck2 box. This dark-matter(DM)-only simulation consists of a periodic cube of 1.5 Gpc in comoving length containing 38403 DM particles. As the name suggests, this simulation assumes the best-fitting cosmological parameters from the Planck Collaboration XIII (2016): h = 0.6777 for the reduced Hubble parameter, n = 0.96 for the primordial spectral index, σ8 = 0.8228 for the amplitude of the mass density fluctuations in a sphere of 8 h−1 Mpc comoving radius, and ΩΛ = 0.692885, Ωm = 0.307115, and Ωb = 0.048206 for the density parameters of dark energy, matter, and baryonic matter respectively.

The Lagrangian regions to be re-simulated at higher resolution are identified at z = 0 as the volume centered on the selected massive clusters (all with virial masses3 greater than 1.2 × 1015M) and extending for a radius of 22 Mpc. With respect to the MACSIS sample (Barnes et al. 2017b), the Three Hundred Project set has the advantage of being a volume-limited mass-complete sample for all objects with M500 >  6.5 × 1014M. Nevertheless, in this paper we extend the sample to smaller mass objects to study a possible mass dependence of the results (Appendix A).

Initial conditions are generated at the initial redshift z = 120 with the GINNUNGAGAP4 code by refining the mass resolution in the central region and degrading it in the outer part with multiple levels of mass refinement. This step allows us to keep the information on the large-scale tidal fields without raising the computational cost. The high-resolution dark-matter particle mass is equal to mDM = 1.9 × 109M, while the initial gas mass is equal to 3.5 × 108M. The gas softening of the simulations is fixed to be 15 kpc in comoving units for z >  0.6 and 9.6 kpc in physical units afterwards. The minimum value of the SPH smoothing length allowed is one-thousandth of the gas softening, but it is de facto around 1 kpc.

Since the Three Hundred Project was born as a comparison project, the 324 regions are re-simulated with three hydrodynamical codes: GADGET-X (Rasia et al. 2015), GADGETMUSIC (Sembolini et al. 2013), and GIZMO-SIMBA (Davé et al. 2019). It is beyond the scope of this paper to compare different codes; therefore in this work we refer only to the GADGET-X sample.

The Tree–Particle–Mesh gravity solver of GADGET-X corresponds to that of the GADGET3 code, which is an updated and more efficient version of the GADGET2 code (Springel 2005). GADGET-X includes an improved SPH scheme with artificial thermal diffusion, time-dependent artificial viscosity, high-order Wendland C4 interpolating kernel and wake-up scheme as described in Beck et al. (2016). Radiative gas-cooling depends on the metallicity as in Wiersma et al. (2009). The star formation and thermal feedback from supernovae closely follow the original prescription by Springel & Hernquist (2003) and are connected to a detailed chemical evolution and enrichment model as in Tornatore et al. (2007). More details on the chemical enrichment model are presented in Biffi et al. (2017, 2018) and Truong et al. (2019). Finally, the gas accretion onto super-massive black holes powers AGN feedback following the model by Steinborn et al. (2015), which considers both hot and cold accretion (see also, Churazov et al. 2005; Gaspari et al. 2018). The impact of these physical processes on the ICM properties has been discussed in relation to observed quantities in several papers. Relevant for this work, it is worth mentioning that previous simulated samples carried out with this code were shown to broadly agree with observed gas density and entropy profiles (Rasia et al. 2015), pressure profiles and ICM clumpiness (Planelles et al. 2017), and global ICM quantities (Truong et al. 2018). Li et al. (in prep.) find that the gas density and temperature profiles of the Three Hundred sample are in good agreement with the observational results of Ghirardini et al. (2019) at around R500.

### 2.1. Generation of maps

For the current analysis, in each Lagrangian region we consider the most-massive clusters at z = 0 without any low-resolution particles5 within R100, which is close to the virial radius for our cosmology. Whenever the main object had at least one low-resolution particle within R100 we considered the second most massive cluster. Seven regions have no available objects with M500 >  3 × 1013M, which is the mass limit that we impose. For each cluster, we produced with the code Smac6 (Dolag et al. 2005) three X-ray surface brightness maps in the soft-energy band, [0.5–2] keV, along three orthogonal lines of sight. The map is created by summing over the contributions of the gas particles that are in the hot phase – meaning with density lower than the density threshold for star formation7 – and emit in the X-ray band – meaning with temperature above 106 K. Each particle emissivity is weighted by a spline kernel of width equal to the gas particle smoothing length. The center of each map coincides with our definition of the system theoretical center: the position of the minimum of the potential well. From this position, all objects have their R500 radii inscribed in the map8, the side of which is fixed equal to 4 Mpc and is divided into 1024 pixels on a side, leading to a physical resolution of 3.9 kpc per pixel. The integration length along the line of sight is equal to 10 Mpc (5 Mpc from each side with respect to the center, corresponding to about three to four times R500). As an example, three maps are shown in Fig. 1, one for each of the first three morphological classes (very regular, regular, intermediate/irregular) described in Sect. 3. Fig. 1. X-ray maps of three clusters, which are representative of the first three classes described in Sect. 3.1: Very Regular (left panel), Regular (central panel), and Intermediate/Irregular (right panel). The three objects have a similar mass, M500 ≈ 1 × 1015 M⊙, corresponding to R500 of about 1.5 Mpc (shown as a white circle in the images). The black lines are the log-spaced iso-flux contours smoothed over 4 pixels.

### 2.2. Center definitions and ICM ellipticity

As anticipated, our reference theoretical center for the analysis on the simulated clusters is the position of the minimum of the potential well. This center is used to compute the gas profiles. However, we consider two additional centers identified from the X-ray surface brightness maps: the emission peak and the center of the iso-flux contours. Both are commonly adopted in X-ray analysis since the former can be easily identified and the latter is generally considered to be close to the center of mass.

The first center corresponds to the pixel of maximum flux and is therefore dubbed “MF”. At first, its identification is performed automatically. However, whenever its distance from the minimum of the potential, DMF, exceeds 0.4 × R500, we proceed to visually inspect the individual map. During this process, we verify whether the X-ray peak is associated with a denser companion rather than to the main cluster. In this case, we mask the secondary object by excluding the pixels associated with it, and recompute the maximum. In the large majority of the other cases, MF is located at a distance of a few pixels from the map center. For this, we express DMF in pixels rather than in units of R500 which could be misleading because one single pixel corresponds to a different fraction of the cluster radii. Considering all the maps, the median value of this distance is equal to 18 kpc. Rossetti et al. (2016) looked at a similar estimator in a large sample of massive clusters selected from the Planck catalogue. They found that the observed distance between the X-ray peak (i.e., MF) and the center of the brightest cluster galaxy (most of the time coincident with the minimum of the potential well; see Cui et al. 2016) was equal to 21.5 kpc, which is very close to the median value found in our sample.

The second center is the center of the ellipse that best describes the iso-flux contours of the images drawn at around 0.8 R500 and it is dubbed “CE”. In practice, we follow this procedure: we compute the mean of the flux of all pixels at that distance from the map center; then, we use the mean flux value as a threshold to separate two regions with pixels below and above this limit; finally, we recognize as a contour the border between the two regions. Whenever multiple contours are identified within the maps, we consider the longest one. We expect that the ellipse center better approximates the center of the mass distribution on large scales. Since the distance from the minimum of the potential well, DCE, has a broader distribution than for DMF, we measure it in units of R500.

In addition to CE and DCE, we also save the value of the ellipticity of the best-fitting ellipse, which is defined as ε = (a − b)/a, where a is the major axis and b the minor one.

In this preparatory phase, we discarded 29 clusters (and their associated 87 maps) because we encountered difficulties in the determination of their best-fitting ellipse. Visual inspection of these maps confirms that these objects are indeed associated with systems that are either interacting with other massive clusters or have extremely irregular flux maps. Their peculiar morphology prevents any rigorous identification of either centers. Moreover, for a large number of them we recognize that the minimum of the potential is associated with neither object but is located in between the interacting systems. In these cases, no profiles, either 2D or 3D, are informative; rather they will most likely introduce an uncontrollable bias in the results9.

All the other images (precisely, 864 maps corresponding to 288 clusters) and their measured DMF, DCE, and ε are used to build the cluster subsamples as outlined below. Observational works in the literature have used different morphological parameters to classify the regularity of the X-ray appearance of clusters (e.g., Rasia et al. 2013; Mantz et al. 2015; Lovisari et al. 2017) and their connection with mass bias has already been studied (e.g., Jeltema et al. 2008; Piffaretti & Valdarnini 2008; Rasia et al. 2012; Cialone et al. 2018). Here we focus on a more easily derived quantity, the ICM ellipticity, which has been shown to correlate well with the mass-accretion history (Chen et al. 2019) and the overall dynamical state of the clusters (Laganá et al. 2019). Furthermore, Mantz et al. (2015), studying 350 clusters observed in X-ray, showed that the ellipticity is a good proxy to select either very relaxed or, alternatively, very dynamically active clusters.

## 3. Cluster classification

The main goal of the paper is to connect measurements derived from X-ray or SZ observations to the intrinsic estimate of the mass bias. Therefore, the primary classification of the paper is based on the 2D analysis of the X-ray surface brightness maps and uses the parameters introduced before. Since we produced three images for each cluster, a system can be part of more than one 2D subsample. Our study also involves quantities that are measured directly from the simulated clusters in 3D, such as the intrinsic clumpiness, gas density, temperature, and pressure profiles, and the bias of the hydrostatic-equilibrium mass. For this reason, we also introduce a 3D classification.

### 3.1. Two-dimensional classification

We use the computed values of DCE, DMF, and ε from each image to sort the cluster maps and broadly divide them into five classes: very regular (VR), regular (R), intermediate/irregular (IR), very irregular (VI), and extremely irregular (EI). We stress that the classification aims at distinguishing the regular systems and very irregular systems from the bulk of the cluster population.

The division is based on the two parameters linked to the best-fitting ellipse of the external iso-flux contour, because after looking at some images we recognize that both DCE and ε are sensitive to even minor mergers (in agreement with similar indications from observational analyses; see, e.g., Lopes et al. 2018). The parameters are shown in Fig. 2, color-coded with the classes presented below. In contrast, the third parameter, DMF, does not always reflect the irregular morphology or the disturbed state of the ICM, especially if only the cluster outskirts clearly manifest a nonrelaxed status. In other words, for the purpose of the classification, the parameter DMF is not effective at separating regular from intermediate/irregular objects. Nonetheless, it remains very useful for immediately identifying very/extremely irregular systems. In the following, we therefore impose limits on DMF as a secondary condition. Fig. 2. Distribution of DCE (in units of R500) and ε for all maps in the categories: VR (blue triangles), R (red pentagons), IR (green squares), and VI (magenta circles). Top and right panels: abundance of DCE and ε of each class.

After a first automatic classification made on the basis of the parameters as detailed below, we visually inspect all 864 maps to derive the final classification:

Very Regular (VR). 25 maps. We pre-select as part of this class all objects with DCE ≤ 0.052 R500 and ε ≤ 0.1. These threshold values correspond to about the 20th percentiles of the corresponding distributions. In addition, we further impose that DMF ≤ 5 pixels (less than 20 kpc), even though DMF is within 1 or 2 pixels for most of the pre-selected objects. We visually check all maps and keep in this class only those with regular shapes and without substructures. We also add one map to this class even if its DCE is larger than the imposed limit because it has the roundest iso-flux contours (ε = 0.012). We also include three other maps that appear very regular despite having a slightly larger ellipticity (ε >  0.1) with respect to the rest. One cluster has all three projections classified as VR, and four objects have two projections in the VR class and the third in the (adjacent) regular class. Indeed, imposing tight limits on DCE and ε leads us to select objects that are regular in more than one projection and therefore likely to be in a truly relaxed dynamical state. In this class, the median values of the ellipse-center distance and of the ellipticity are: DCE = 0.035 R500 and ε = 0.073.

Regular (R). 102 maps. The second class again includes regular images, but the iso-flux contours are allowed to have a small mis-centering (DCE ≤ 0.12 R500) and to be more elliptical (ε ≤ 0.15, with three exceptions at ε ∼ 0.22 but DCE ≤ 0.05 R500). In addition, small substructures can be present within R500 and the condition on DMF is more relaxed (DMF <  15 pixels – about 60 kpc – even though DMF is within 5 pixels for the majority of the clusters). The following values correspond to the median CE distance and ellipticity: DCE = 0.057 R500 and ε = 0.116.

Intermediate/Irregular (IR). 424 maps. Their CE center can have a non-negligible offset with respect to the minimum of the potential well (but still DCE <  0.2 R500) and the axes of the best-fitting ellipse are characterized by a ratio that more strongly departs from spherical symmetry (ϵ <  0.4 for 95% of the IR clusters and ϵ <  0.5 for 98.5% of them). Some maps with parameters within the limits of the previous two classes are classified as IR for their asymmetric emission or the presence of some larger substructures. No limits on DMF are imposed. Another 13 maps were originally assigned to this class but then removed from the analysis because their 2D profiles centered either in CE or MF do not reach R500. We recall that in this class there might be objects that other authors could classify as “intermediate regular”. For example, the rightmost map in Fig. 1 still shows almost regular iso-flux contours at about R500 without major substructures. The median values of distance and ellipticity for the irregular class are: DCE = 0.107 R500, and ε = 0.211.

Very Irregular (VI). 118 maps. All these maps have DCE between 0.2 and 0.5 R500, implying that either the gas distribution is extremely disturbed at larger distances or they have significant substructures that impact the ellipse fitting at about R500. Also in this case, we do not consider any threshold on DMF. The median value of the ellipse center distance is significantly larger, DCE = 0.257, while the median value of the ellipticity is ε = 0.268 (95% of the VI objects have ϵ <  0.5). The matching between the 2D quantities extracted from the surface brightness maps and the 3D cluster profiles needs extra care because of possible mis-centering between the two sets of information. We therefore exclude this class from the main paper.

Extremely Irregular (EI). 195 maps. These appear extremely disturbed and to strongly interact either with another cluster of similar mass or with multiple groups. The distance between the maximum of the X-ray flux and the minimum of the potential well can be significant. Identifying the center, and thus extracting the profiles in both 2D and 3D, is challenging for most of these objects. Since these features are recognizable in more than one projection, we discard the three maps of all the 65 EI objects from the analysis.

We note that in the paper by Mantz et al. (2015) ellipticity values around 0.2 ± 0.1 are linked to clusters with a mixture of dynamical states, while clusters with ε <  0.12 or with ε >  0.3 are almost exclusively relaxed or unrelaxed. These limits closely match with the ε thresholds used in our classifications.

The rest of the analysis is focused on the approximately 550 maps of the first three classes. One-quarter of them show regular morphology, which means that the images are classified as either VR or R. The analysis on the VI clusters is excluded from the main text but their results are reported in Appendix B since they might be useful for comparisons with observational samples which include very disturbed systems.

### 3.2. Three-dimensional classification

In Sects. 6 and 7, we show results from the 3D analysis of the simulated clusters. In these sections, whenever we compare a 3D quantity with a 2D measurement we use the 2D classification. However, when we compare 3D quantities amongst themselves, such as clumpiness factor and mass bias, the classification based on the maps is less appropriate because the same cluster might belong to different 2D classes depending on projection. Therefore, whenever we rely on 3D properties, we refer to two simplified classes: the 3D regular, R3D, and the 3D irregular, IR3D, clusters. In the former class we include all objects that have at least two projected images classified as R or VR; in the latter, we instead require that at least two projected images are in the IR class. We do not consider any cluster that has two projections in the VI or in the EI classes. The classes R3D and IR3D contain 30 and 150 clusters, respectively.

Basic properties of the cluster subsamples corresponding to the studied 2D and 3D classes are summarized in Table 1. For completeness, we report that the mass range of the 29 objects for which the ellipse fitting cannot be performed is [0.5 − 14.2]  ×  1014M, while this is [0.4 − 26.0]  ×  1014M for the EI sample, where the second most massive clusters has M500 = 1.6 × 1015M. The mass coverage is therefore very similar to the other subsamples whose mass distribution is shown in Fig. 3. The top panel refers to the 2D classification, the bottom to the 3D classification. Each histogram is normalized by the number of objects of each class. Ninety percent of the mass distributions in the IR and (VR+R) classes have values in the range 7 × 1013M <  M500 <  1.4 × 1015M. These clusters are the main objects in their respective Lagrangian regions. Those with M500 <  3 × 1014M are instead the second most massive clusters of their Lagrangian regions, introduced in the sample because of contaminated particles within the primary object. We conclude that the minimum, median, and maximum mass values of the samples (Table 1) and the overall distribution of masses (Fig. 3) are similar, or in other words there is no particular selection mass bias.

Table 1.

Basic properties of the subsamples for the 2D analysis (upper part) and the 3D analysis (lower part). Fig. 3. Distribution of M500 for the 2D classification (top panel) and the 3D classification (bottom panel).

## 4. Methods

This section describes how we derive all the ICM quantities. It is important to note that all of them are obtained as radial profiles (either in 2D in Sect. 4.1 or 3D in Sects. 4.2 and 4.3). The profiles are also used to evaluate a particular quantity at R500. Indeed, we prefer to consider the quantity mean value computed by averaging it over the radial bins from 0.9 to 1.1 R500, instead of using the interpolation. This procedure limits the dependence of our results on the precise radial binning adopted.

### 4.1. Two-dimensional gas inhomogeneity

Starting from the X-ray flux maps, we extract different indicators of the gas inhomogeneities:

ε: the ellipticity of the best-fitting ellipse to the external iso-flux contour (Sect. 2.2);

σA: the azimuthal scatter of the X-ray surface brightness profiles (various ways of calculating this quantity are described below);

MM: the ratio between the mean and median of the X-ray surface brightness profiles (this is a byproduct of the σA measurements) minus one.

We note that we use the information from the entire map without removing any substructures before the analysis and we never apply any extrapolation to extract the profiles. For this reason, 15 maps (1 in the R class and 14 in the IR class) are not considered when discussing the properties of the quantities evaluated at R500. Indeed, in these few cases the projected annulus around R500 computed from the X-ray center (not coincident with the center of the map) is not entirely contained in the map. To measure MM and σA we center on the X-ray peak, MF, and divide the image into 12 sectors. We then derive the corresponding 12 surface brightness profiles in radial bins spanning from 0.4 R500 to 1.2 R500 and linearly equally spaced with a distance equal to 5% of R500 (see, e.g., the top panels of Fig. 4 for each of the maps shown in Fig. 1). From the 12 surface brightness profiles, we compute the median and mean (solid thick line in Fig. 4) surface brightness profiles and extract: the ratio between the median and the mean value of the surface brightness computed in each radial bin, similarly to Eckert et al. (2015): (1) Fig. 4. Upper panels: surface brightness profiles (gray lines) extracted from the 12 sectors centered on MF of the three images of Fig. 1, representative of the classes: VR (left panel), R (central panel), and IR (right panel). The thick solid lines show the mean profile. The median profile is not shown for clarity. Lower panels: azimuthal scatter profiles derived as in Eq. (2) and computed with respect to the median profile and centered on MF, , for clarity simply referred to as σA in the plot.

and the azimuthal scatter, defined in Vazza et al. (2011) and Roncarelli et al. (2013): (2)

where the reference profile in the formula, ⟨X(r)⟩, is either the mean or the median.

We then repeat the entire procedure by centering the sectors on CE (Sect. 2.2) rather than on MF.

In conclusion, for each map we have two values of the difference between the mean and median, MMCE and MMMF, and four different versions of the azimuthal scatter profile, , , , and (examples of the latter are shown in the bottom panels of Fig. 4). We compare their ability to capture gas inhomogeneities in Sect. 5; we investigate their reliability as a proxy for the intrinsic 3D clumpiness level in Sect. 6; and finally we relate them to the HE mass bias at R500 in Sect. 7. For the latter-mentioned study, we include ε as well as another variation of the azimuthal scatter: σA, R. This quantity is defined as the mean value of the azimuthal scatter averaged over the entire radial range (from 0.4 to 1.2 R500). We consider it in relation to the HE mass bias at R500 because the equilibrium assumption at that radius can be broken by a clump that already moved away by generating motion in the ICM.

### 4.2. Clumpiness factor

As stated in the introduction, the level of clumpiness can only be indirectly estimated from X-ray observations, whereas the clumpiness factor can be precisely measured in simulations by adopting its definition: (3)

where ρ is the gas density and the brackets, ⟨⟩, indicate the average taken over the region of interest (Mathiesen et al. 1999). More precisely in our analysis, based on SPH simulated clusters, we compute the clumpiness factor by adopting the following formula discussed in Battaglia et al. (2015) and Planelles et al. (2017): (4)

where mi and ρi are the mass and density of the ith gas particle. The sum is extended over all the gas particles used for the observational-oriented quantities, in other words not star-forming and with temperature above 106 K in order to consider only X-ray-emitting gas. We extract the clumpiness profiles by adopting the same radial range and binning as for the azimuthal scatter even though the shells are now spherical and centered on the minimum of the potential well.

Past studies have stressed the importance of additionally computing the residual clumpiness (Roncarelli et al. 2013; Vazza et al. 2013; Zhuravleva et al. 2013; Khedekar et al. 2013). For example, in Roncarelli et al. (2013) the residual clumpiness is derived after excluding in Eq. (4) the densest particles of the radial shell, defined as the particles that account for 1% of the total volume of the shell. This work emphasizes how the residual clumpiness rather than the clumpiness factor is a more appropriate proxy for large-scale inhomogeneities.

However, in the set of simulations analyzed here the difference between clumpiness and residual clumpiness is not as evident as in previous analyses. To reach this conclusion, we computed the residual clumpiness following Roncarelli et al. (2013): we ordered all particles by their density, ρi, and removed the densest ones until Σimi/ρi = Σi(Vi) = 0.01 × Vshell, where Vi and Vshell are the volumes associated with the i-th particle and the shell, respectively. For the three representative cases already shown in Fig. 4, we plot both the clumpiness and the residual clumpiness10 in Fig. 5. The small offset between the two sets of curves should be compared with Fig. 4 of Roncarelli et al. (2013), which shows that clumpiness and residual clumpiness can differ by one order of magnitude within R500. Among the irregular (IR3D) sample, where we expect the largest difference between the two clumpiness profiles, we find that 95% of the objects have a maximum difference lower than a factor of 1.5. The artificial conduction introduced in our code indeed leads to better mixing of the medium and consequently to a net reduction in the number of clumps (Biffi & Valdarnini 2015; Planelles et al. 2017). In addition, the cores of the main halos and of the substructures are smoother. The difference in the clumping level of SPH clusters and adaptive-mesh-refinement (AMR) objects shown in Rasia et al. (2014) are now almost completely erased for nonradiative runs, and cosmological simulations that include AGN feedback such as the one investigated here. In the rest of this paper we therefore focus on correlations with respect to the clumpiness measurement to emphasize the signal of inhomogeneities on all scales. Fig. 5. Upper panels: clumpiness (solid line) and residual clumpiness (dashed line). Lower panels: profiles of the hydrostatic-equilibrium mass bias as in Eq. (12). The quantities (see footnote 8) are measured in 3D for the clusters whose maps are representative of the 2D classification of Fig. 1.

### 4.3. Hydrostatic-equilibrium mass bias

From the 3D distributions of the gas particles we compute the gas density, temperature, and pressure radial profiles. The profiles are computed by centering on the minimum of the potential using radial bins that are logarithmically equally spaced, with the external radius of each shell fixed to be 1.1 times the inner radius. The binned gas profiles are not directly used in the hydrostatic mass equation. Instead, we search for best-fitting analytic formulae that could appropriately reproduce each profile. The formulae adopted are taken from observational work. In this way, we follow the X-ray procedure more closely and, at the same time, we prevent singularities in the derivatives of the gas profiles, which can emerge near to negative or positive spikes in the gas density (see e.g., the noisy profiles of Biffi et al. 2016 or Cialone et al. 2018 who directly use the intrinsic numerical gas profiles in the HE mass equation). The data points are always fitted over the radial range between 0.4 R500 and 1.2 R500 to search for the best constraints on the profiles around R500.

#### 4.3.1. Gas density

Rasia et al. (2006) and Nagai et al. (2007b) proved the reliability of the X-ray reconstruction of the gas density profiles. Both works analyzed simulated clusters and produced mock X-ray images including Chandra ACIS-S and ACIS-I responses. The quantities obtained from the X-ray analysis, such as the surface brightness profiles and the projected and de-projected gas density profiles, were found to agree with the input simulated data set (see also Meneghetti et al. 2010 who tested different X-ray procedures and Avestruz et al. 2014 who extended the X-ray comparison to large radii). In light of these previous tests, which were also based on different exposure times, we decided to pursue a straightforward analysis of the gas density profiles of the simulated objects rather than a more complicated analysis of mock images (albeit, see Henson et al. 2017).

The gas density, ρ, is computed as the total gas mass in the spherical shell divided by the shell volume, ρ = Σmi/Vshell. Each gas density profile is fitted by the (simplified) parametric formula by Vikhlinin et al. (2006): (5)

where ρ0, r0, rs, β, and ϵ are free parameters. With respect to the original formula proposed by Vikhlinin et al. (2006), we impose, as often done, that the parameter γ be equal to 3, and we avoid the second beta model that describes the inner core because we are only interested in obtaining a precise analytic fit of the gas density slope around R500 and the radial range investigated excludes the central 40% of R500.

In Sect. 7, we refer also to the asymptotic external slope (for r ≫ r0 and r ≫ rs) of the analytic density profile to correct for the HE mass bias. Accordingly to the adopted formula, this is given by 𝒟 = 3β + ϵ/2. More rigorously, when the formula was introduced, the second term was included to improve the description of the external slope which observations suggested could not always be represented by the simple beta model (the first term in the formula). Therefore the scaling radius of the second term, rs, should be larger than the core radius of the beta-model, r0. However, since we did not impose any particular condition on the relative values of the two scale radii the second term of the expression does not necessarily represent the trend of the density profile in the outskirts. To confirm that 𝒟 is a good representation of the external density slope, we consider the density profiles obtained from the best analytic fits and calculate their derivative at a large distance, precisely at 100 times their maximum scale radius (either r0 or rm). We find that on average there is no difference between the resulting derivative and 𝒟 and that the maximum deviation between the two is of the order of a few thousand. This quick test validates the use of 𝒟 as a reference for the asymptotic external density slope.

#### 4.3.2. Temperature

It is significantly more complicated to test whether the temperature computed in simulated clusters reflects the X-ray temperature from spectroscopic analysis. From a numerical point of view, the temperature is measured as a weighted average over an ensemble of gas elements. It is now well established that using the X-ray emission to weight the temperature leads to biased results (Gardini et al. 2004) and that another definition should be used to reproduce Chandra or XMM-Newton measurements: the spectroscopic-like temperature (Mazzotta et al. 2004). Nonetheless, the spectroscopic-like temperature reproduces the projected temperature obtained directly from the spectra, which is not the temperature that enters the equation of hydrostatic equilibrium used to derive the X-ray mass. Instead, X-ray observers deproject the projected temperature profile using the gas density obtained from the imaging as input for the spectroscopic-like weighting. The final X-ray deprojected temperature profile, already deconvolved from the instrumental response, can then be considered as the “true” (un-weighted) temperature profile (e.g., see the review by Ettori et al. 2013). The accuracy and precision of this procedure depend on the ability to correctly treat the background and on the exposure times, since at least 1000 counts should be collected to measure the temperature. Having a large exposure time allows reconstruction of the profiles in finer radial bins. This improves the deprojection technique and reduces possible biases due to the co-existence of multi-temperature gas components. Another source of complication is that the temperature distribution in simulated clusters depends not only on the ICM physics included in the simulations, such as thermal conduction, viscosity, and sources of feedback, but also on details of the hydrodynamical methods employed (Vazza et al. 2011; Rasia et al. 2014; Sembolini et al. 2016; Richardson et al. 2016; Cui et al. 2018; Power et al. 2019; Huang et al. 2019).

Owing to all these problems, we decided to follow a theoretical approach and thus to use the mass-weighted temperature: T = Σ(mi × Ti)/Σmi, where the i - th gas particle has a temperature greater than the lowest energy band of current X-ray telescopes, 0.3 keV (≈3.5 × 106 K). Indeed, the mass-weighted temperature is the one to be considered for the derivation of the cluster mass under the assumption of hydrostatic equilibrium11. Furthermore, in some observational analyses (e.g., Vikhlinin et al. 2006) the masses are similarly derived by weighting the temperatures by the gas mass.

The temperature profiles are fitted by the functional form (6)

where T0, r0, αrmT, and βT are free parameters. With respect to the original formula, introduced by Vikhlinin et al. (2006), we neglect the extra term describing the temperature drop in the core region for the same reasons listed at the end of the previous section.

#### 4.3.3. Pressure

The pressure profile is measured starting from the pressure of the individual gas particles. These profiles are very similar to those obtained by multiplying the gas density and mass-weighted temperature profiles.

The pressure profile is described by a generalized Navarro-Frenk-White model (Nagai et al. 2007a; Arnaud et al. 2010): (7)

where P0, rp, αp, and βp are free parameters, and the internal slope, γp, is fixed equal to 0.31 as in Planelles et al. (2017). The asymptotic external slope (for r ≫ rp) of the analytic pressure profile is given by βp and, similarly to 𝒟, it is used in Sect. 7.

#### 4.3.4. Goodness of the fit

On top of the best-fitting parameters, we also save the normalized root-mean-square (NRMS) value as a measure of the goodness of the fit: (8)

where the sum is extended over all radial bins.

We opt for normalizing the residuals to obtain comparable values of the three gas profiles (density, temperature, or pressure) from the fitting procedures. We generate the three respective distributions of the NRMS to identify the clusters poorly described by their best-fitting curves (Sect. 7). Namely, these are the objects that belong to the highest quintile in any of the three NRMS distributions.

#### 4.3.5. Hydrostatic mass equations

For consistency with the measurements of gas inhomogeneities, the analytic profiles are then computed adopting the same radial binning as that in Sect. 4.1 and then folded into three versions of the hydrostatic equilibrium equation to derive an estimate of the total mass: (9)

This expression has been used to exploit the advantages of both SZ and X-ray signals in providing the pressure and the gas density, respectively, with good accuracy at large distances from the center (see an early study by Ameglio et al. 2009 or the recent works by Eckert et al. 2019; Ettori et al. 2019); (10)

This latter equation is typically used in X-ray analyses where the gas density is derived from the imaging and the temperature from spectroscopy (see review by Pratt et al. 2019, and references therein); (11)

this latest hybrid version helps to separate the influence on the mass-bias calculation of the linear multiplicative factor and the term with the derivatives sum.

In all equations, kb is the Boltzmann constant and A = 1/(Gμmp) = 3.7 × 1013M keV−1, where G, mp, and μ are the gravitational constant, the proton mass, and the mean molecular weight, equal to 0.59 in our simulations. To consider the same multiplicative factor, A, in all expressions, the pressure in Eq. (9) is computed from the gas mass density rather than the electron number density.

#### 4.3.6. Hydrostatic mass bias

The hydrostatic-equilibrium mass is a locally defined quantity because all gas profiles and their derivatives are measured or computed at a precise radius. The bias between the HE mass and the true mass can therefore be evaluated within each radial bin. We define the bias parameter as (12)

for each of the three versions of MHE. The bias, bHE, is zero when the HE mass coincides with the true mass, while it is negative or positive for overestimated or underestimated values of MHE.

For the representative clusters shown in Fig. 1, we show the corresponding profiles of (1 − bHE, X) in the bottom panels of Fig. 5.

## 5. Results: 2D gas inhomogeneity

We begin the presentation of our results with a comparison among the profiles obtained by the various estimators of the 2D gas inhomogeneities presented in Sect. 4.1. We focus here on the trend of these estimators over the entire radial range and therefore we do not consider the X-ray ellipticity, ε, which was computed at a fixed distance.

### 5.1. Scatter with respect to the median or mean

The two estimates of the azimuthal scatter σA, computed with respect to the mean or to the median surface brightness profiles, centered on MF are compared in Fig. 6 as the ratio between the two options. To prepare this plot, we first calculated the ratio between each individual pair of σA and then, in each radial bin, we computed the median of the ratio distribution. The median is shown as a solid line and the distribution between the 16th and the 84th percentiles as a shaded area. Fig. 6. Ratio between the azimuthal scatters computed in Eq. (2) with respect to the median, , and to the mean, and centered on MF. The solid lines represent the median of the ratios in each radial bin over the VR (blue), R (red), and the IR (green) subsamples. The shaded regions comprise the distribution between the 16th and the 84th percentiles.

For the VR class, the profiles are between a few and five percent higher than the profiles over the entire radial range; for the R class, they are ten percent higher. Three quarters of the VR (R) objects have ratios smaller than 12 (18) percent at R500. The two choices of the azimuthal scatter therefore provide similar results for the regular classes. We verify that this result holds independently of the chosen center (MF or CE). In the IR class, the scatter measured with respect to the median surface brightness profile presents larger fluctuations than the scatter measured with respect to the mean surface brightness profiles. Indeed, at all radii, there is a difference between the two σA of about 20–30% using MF as center (green curve in Fig. 6) and 25–40% using CE as center.

The median surface brightness profile is always smaller than the mean one because it is less affected by (and thus more stable against) the presence of substructures (e.g., Zhuravleva et al. 2013). The azimuthal scatter computed with respect to the median will therefore enhance the effect of gas inhomogeneities, including not only substructures but also overall large-scale irregularities. The qualitative results found in this section are in line with what is presented in previous studies (Zhuravleva et al. 2013; Khedekar et al. 2013) based on different simulations. Here, it is important to stress the quantitative evaluation of this effect, since the simulations analyzed are characterized by a higher level of mixing with respect to several previous analyses: using an azimuthal scatter computed with respect to the median enhances the imprint of inhomogeneities at R500 by at least 30% for half of the IR objects and 10% for half of the R clusters.

### 5.2. Azimuthal scatter and MM

The findings of the previous section are clearly connected to the parameter MM (Eq. (1)). Indeed, the median profiles of the MM parameters (not shown) have the same trends as the solid lines of Fig. 6. In this figure, the small difference found between the two scatter estimates in the regular classes (VR and R) essentially reflects the similarity between the mean and median surface brightness profiles. For these objects, MM always shows little deviation from zero; for example, at R500 its median value is about 0.05. This finding reflects the fact that the objects in these classes are characterized by an homogeneous and symmetrical X-ray distribution around MF, which by the definition of the VR and R classes is close to the minimum of the potential well and to the center of the best-fitting ellipse.

The increase in the scatter in the IR class results from a larger offset between the median and the mean of the surface-brightness profiles over the 12 sectors. The latter is on average 15–20% higher than the former at all radii, but MM can reach a value of 1 in about 10% of the IR objects at R500, implying that the mean is twice as high as the median, with clear consequences for the two derived azimuthal scatters. That said, at that radius the median value of MM is much lower and equal to 0.22 and 80% of the objects have MM <  0.6.

### 5.3. Effects of centering

We proceed to assess the impact of the choice of center (MF vs. CE) for the 12 sectors. Based on the results of Sect. 5.1, we consider the azimuthal scatter computed with respect to the median profile in order to be more sensitive to the presence of inhomogeneities. The median ratio is computed following the above procedure and is shown in Fig. 7. Fig. 7. Ratio between the azimuthal scatters centered on MF, , and centered on CE, and computed in Eq. (2) with respect to the median. The color code and meaning of the shaded area and solid lines are the same as in Fig. 6.

The classes VR, R, and IR show similar behaviors at all radii. The majority of the objects in each class have a smaller scatter when the sectors are centered on MF. This result is expected for the innermost region since the X–ray emission in the central part is supposedly smoothly distributed around its maximum. Here, in fact, tends to be 20–30% smaller than the scatter computed with the sectors centered on CE for all classes and reaches a difference equal to or greater than a factor of two for one-quarter of the IR systems. On the other hand, for radii between 0.8 and 1.2 R500, the median of the ratios is almost constant and approaches unity with a small deviation of about 5–8%. Even though the overall difference between the two scatters at R500 is very small, the large majority of the systems (70 and 60% in the VR class and in the IR class) have a ratio below 1, while the naive expectation was that the difference between the two centers would disappear in the external regions. The highest discrepancies (those with ) seem to be caused by massive and extended substructures that not only distort the ellipse (and thus its CE centre), but also increase at about R500.

It is important to stress that while the three classes show similar trends for the ratio in Fig. 7, the median azimuthal scatter of the three classes is rather different. Indeed, σA of the IR class is typically twice as high as that of the most relaxed objects (Fig. 8 and following section). As a consequence, 55% of the IR clusters have while only two objects satisfy this condition in the VR class. More extremely, one-quarter of the IR objects have . Fig. 8. Azimuthal scatter computed as in Eq. (2) with respect to the median and centered in MF, . The color code and the meaning of the shaded area and solid lines are the same as in Fig. 6.

Taking all of this into consideration as well as the uncertainties associated with the automatic determination of the best-fitting ellipse and especially for the vicinity of the MF center to the minimum of the potential well, we chose MF as the cluster center when measuring both σA and MM.

### 5.4. Azimuthal scatter for the VR, R, and IR classes

From the previous analysis, we establish that the best choice to compute the azimuthal scatter is with respect to the median and centered on MF. Hereafter, we use the symbol σA to refer to . The median behavior of the azimuthal scatter profiles is shown for the three classes in Fig. 8. The azimuthal scatter grows from 0.2–0.3 to 0.5–0.7 from the VRR classes to the IR class, consistent with Vazza et al. (2011) and Roncarelli et al. (2013).

For the regular systems (VR and R) not only are the profiles of the median values in each radial bin flat but also the dispersion around these values is small, implying that the individual σA profiles show little spread without significant bumps. Vice versa, the IR class is characterized by a significant scatter which increases with radius. Several profiles present spikes at different radii making the distribution highly skewed in all radial bins.

A high percentage of images in the IR class at a certain point have σA >  1. This extreme condition can be verified when there is a flux significantly higher than the median behavior in one or more sectors or when numerous sectors have simultaneously higher and lower emission than the median value. These situations reflect the presence of one or more bright substructures or an asymmetric distribution of the ICM characterized by a pronounced ellipticity. To investigate which has the greatest impact on σA we studied the objects with high values of ellipticity in more detail. We selected the two most elliptical images in the VR class (with ε >  0.1) and the 20 most elliptical images in the R class (with 0.14 <  ε <  0.23). We inspected the maps to make sure that there are no substructures (or even small clumps) close to R500. The maximum value of σA at R500 for all these maps is equal to only 0.6. We also searched among the IR images and found one cluster selected with high ellipticity, ε = 0.34, but without major substructures. Even in this case, the azimuthal scatter at R500 is limited to σA = 0.73. We thus conclude that an azimuthal scatter higher than ≈0.8 is mostly caused by the presence of substructures rather than by elongated X-ray contours.

## 6. Results: 3D clumpiness factor

### 6.1. Clumpiness for the R3D and the IR3D classes

We reiterate that we opt to investigate the 3D clumpiness rather than the residual clumpiness to enhance any signal from gas inhomogeneities and/or small-scale irregularities, and that the clumpiness profiles are always centered on the minimum of the potential well and evaluated in 3D. The median values of the clumpiness factor for each radial bin are shown in Fig. 9 as a solid line. The shaded area includes the distribution between the 16th and 84th percentiles. Since we are not considering any 2D quantity, we are presenting the clumpiness factor profile by dividing the clusters according to the 3D classification into R3D and IR3D. Fig. 9. Clumpiness profiles for the R3D class in red and for the IR3D class in olive green. The solid line refers to the median profile and the shaded area shows the values between the 16th and 84th percentiles of the 𝒞 distribution.

Similar to the azimuthal scatter, the median values of the clumpiness factor profiles within the regular class are flat and have low values (𝒞  ≈  1.05 − 1.08) and very low dispersion. On the other hand, the irregular class shows proof of a slight increase in the clumpiness factor profiles towards the largest radii. In reality, not only does the median value of the distribution grow from about 1.1 to about 1.2 over the considered radial range but the overall distribution of IR3D objects also shifts to higher clumpiness values at farther distances. This trend is consistent with all other studies based on simulations (Planelles et al. 2017; Battaglia et al. 2015) and observations (Eckert et al. 2015), which show how the clumpiness profile gently increases out to R500 (Nagai & Lau 2011; Zhuravleva et al. 2013; Vazza et al. 2013; Roncarelli et al. 2013; Khedekar et al. 2013; Morandi et al. 2013).

### 6.2. Clumpiness factor and azimuthal scatter

The relation between the clumpiness factor and the azimuthal scatter is shown in Fig. 10 for the values of the two quantities computed in each radial bin and for the IR systems. The figure zooms into the part of the plane (𝒞 <  1.5 and σA <  2) where most of the points are located; indeed the median values of the two quantities are equal to 𝒞 = 1.14 and σ = 0.5. For the VR and R classes not considered in the plot, three-quarters of their points are in the bottom-left corner (𝒞 <  1.1 and σA <  0.4). If more than a projection is part of the sample, the same objects are counted multiple times with different values of σA(r) and a single measurement of 𝒞(r). This visualization highlights the link between 𝒞 and σA. The Spearman correlation coefficient (calculated with the IDL routine R_CORRELATE ) is relatively strong corr(𝒞, σA) = 0.60 with null probability of consistency with zero. The correlation is evaluated in about 9,000 points from all maps and using the information in all radial bins. We have verified that this value does not vary when we refer to the residual clumpiness instead of the clumpiness, or when we consider σA computed with respect to the mean and/or centered in CE. Using all points from all classes, we search for the best linear relation between clumpiness factor and azimuthal scatter by employing an outlier-resistant two-variable linear regression routine in IDL (ROBUST_LINEFIT performed with the bisector method). The best-fitting procedure returns the relation: 𝒞 = 1.01 + 0.22 × σA. Fig. 10. Distribution of clumpiness vs. azimuthal scatter for the IR systems and for all radial bins. Only points with 𝒞 <  1.5 and σA <  2 are considered. For clarity, the distribution is shown in small bins of the two quantities. The colors indicate the number of points per bin. The color scale is saturated at ten.

Roncarelli et al. (2013) described the clumpiness factor as a function of both the azimuthal scatter and the radius: (13)

with σ0 and r0 approximately equal to 16 and 6 × R200 when they extract the σA values from the surface brightness maps produced within the same energy band used in this paper ([0.5 − 2] keV). We fit the same relation to our data sets but do not detect any actual need to include the dependence on the radial distance. To confirm this result, we restrict the fitting procedure to σA and 𝒞 computed in three different regions: the first with R <  0.6 R500), the second with 0.6 <  R/R500 <  0.8, and the third with R >  0.8 R500. We retrieve the values of the intercepts and the slopes and find that they are always consistent with each other. This proves that in our simulations the explicit dependence of clumping on the radius, as in Eq. (13), is not required; this most likely depends on the radial range investigated because we focus within R500, where the clumpiness factor is still reduced.

Before investigating the region around R500 in more detail, which is the one we focus on when discussing the HE mass bias, we briefly examine the possible origins of the scattering of the points over the plane shown in Fig. 10. For simplicity, we consider two classes of outliers deviating from the diagonal; each includes 31 points. These are less than 0.5% of the total number of points but they represent the most extreme situations. Outliers in the bottom-right part of Fig. 10 with σA >  1.72 (90th percentile of the σA distribution) and 𝒞 <  1.07 (20th percentile of the 𝒞 distribution) and those to the top-left with σA <  0.3 (20th percentile) and 𝒞 >  1.4 (90th percentile). Within the first outlier class, with large σA but low 𝒞, 22 of the points have R >  0.8 × R500. Most of them are associated with images with projected substructures. These increase σA, being present in the 2D map, but they lie outside the sphere used to compute the clumpiness in 3D, and are therefore present only in one or two projections. The other class of outliers, with low σA and high 𝒞, is linked to the presence of inhomogeneities that cannot be easily identified in the images because they are aligned with the cluster core that dominates the emission. This situation is present at all radii, near and far from the cluster center.

In Fig. 11 we show the relation between σA and 𝒞 at R500. In the top panel, each map is represented by a single point, while each cluster produces three points, all with the same clumpiness value but different σA. In the bottom panel, the azimuthal scatter is instead computed for each cluster as the mean value of the σA of each projection: ⟨σA⟩. Fig. 11. Clumpiness vs. azimuthal scatter computed around R500. Top panel: each σA refers to a separate map. The color code reflects the 2D classification: VR in blue, R in red, and IR in green. Bottom panel: each cluster is represented only once and the azimuthal scatter is averaged over all its considered projections, ⟨σA⟩. The R3D objects are shown in brown, and the IR3D clusters in olive green.

The distribution of points in the top panel resembles that of Fig. 10. The correlation coefficient is similar, being 0.56, and the parameters of the linear fit are identical12. The scatter on the relation is still significant, but drastically reduces when we average the three scatters for each cluster, ⟨σA⟩. The bottom-right outliers (high σA and low 𝒞) are now sparse. Indeed, the average scatter ⟨σA⟩ is reduced because in at least one line of sight the substructure is correctly identified as being external to the cluster and thus does not have an influence on the value of the azimuthal scatter. The top-left outliers (low σA and high 𝒞) have almost disappeared. These were related to objects with substructures aligned with the cluster center and therefore more likely to have a low value of σA only in one projection (the one with the perfect alignment). By reducing the effect of both classes of outliers, the correlation between clumpiness and scatter in 3D is even stronger: corr(𝒞(R500),⟨σA(R500)⟩) = 0.65

From our results we therefore conclude that an azimuthal scatter with a value significantly greater than 1 is a strong indication of substructures either projected or in 3D. On the other hand, substructures can also be masked by the core emission in the case of close alignment along the line of sight. Since this case is difficult to pick up, one might use statistical considerations: among all objects with σA <  0.5 = median(σA), the incidence of 𝒞 >  1.35 is around 5% and the incidence of 𝒞 >  1.2 is around 10%.

To conclude, we checked the correlation between the clumpiness factor at R500 with the other 2D estimators of the gas inhomogeneities: the MM parameter and the ellipticity, ε. In both cases we find a weaker correlation: corr(𝒞, MM) = 0.47 and corr(𝒞, ε) = 0.37. This is not surprising because these estimators are thought to describe the large-scale inhomogeneity rather than the distribution of individual small clumps.

## 7. Results: hydrostatic mass bias

In Sect. 4 we presented different expressions to obtain the mass under the assumption of hydrostatic equilibrium: (i) the X-ray mass, MHE, X, (ii) the SZ/X-ray mass, MHE, SZ, and (iii) what we called the hybrid estimator, MHE, T. The latter is useful for understanding the relative weight of all the factors entering in the HE mass equation, although this is inconvenient to derive from an observational point of view because of the difficulties in obtaining precise temperature measurements in small radial bins in the cluster outskirts. These three estimators (Eqs. (9)–(11)) refer to physically equivalent quantities, but operationally they might lead to dissimilar results because clumps affect the distinct thermo-dynamical quantities differently (pressure vs. gas density vs. temperature; see e.g., Ruppin et al. 2018).

In Sect. 7.1, we compare the three estimates of (1 − bHE) and in Sect. 7.2 we relate the bias to the clumpiness level. For these parts, we consider only the 3D analysis of the simulated sample, therefore each cluster is counted only once and we use the 3D classification. Further below, we attempt to correct for the mass bias using information from the X-ray images and from the gas fitting procedure (Sect. 7.3). It is worth stressing that a solution is effective not only when the median of all corrected biases is close to one (implying that the HE mass is identical to the true one) but also when both scatter and skewness of the bias distribution are reduced. Otherwise, any proposed solution is equivalent to simply adding a constant equal to the median of the bias values to all HE masses. We summarize all results related to the mass bias measured at R500 in Table 2 for both the 3D (top panel) and 2D (bottom panel) subsamples. In Table 3, we instead report the Spearman correlation coefficient between the mass bias and all other investigated quantities.

Table 2.

Summary of HE mass bias result at R500.

Table 3.

Spearman rank correlation coefficient and the number of standard deviations from the null-hypothesis expected value (in parenthesis) computed between the mass biases and different quantities listed in the first column and obtained in the 2D and 3D analyses at R500.

### 7.1. The mass bias

We show in Fig. 12 the median HE mass bias profiles, (1 − bHE), as expressed in Eq. (12) for the regular (brown) and irregular (olive green) subsamples defined in Sect. 3.2. From the left to the right panel, MHE is given by Eq. (10) for bHE, X, Eq. (9) for bHE, SZ and Eq. (11) for bHE, T. For all expressions and classes, the median bias is increasingly departing from one as the radius grows: we find that the total mass at 0.5 R500 (which is approximately R2500) and at R500 is underestimated by 5–10% and 10–15% respectively. The irregular systems tend to have a higher bias by between a few and 5% and they have a much wider spread. Indeed, the area between the 16th and 84 percentiles of the bias distribution of the IR3D subsample exceeds the respective percentiles of the R3D distribution. These findings are common to most studies based on the direct analysis of simulated samples. At R500 the shaded area is above the value (1 − b) = 0.80 in all panels. In the entire sample (R3D plus IR3D), we find that only four clusters (less than 2.5%) have (1 − bHE, X) < 0.70, in conflict with the mass bias required to solve the discrepancy on the cosmological parameters derived from cluster number counts and cosmic-microwave-background power spectrum (see Salvati et al. 2019, and their discussion). Fig. 12. Median profile of the HE mass bias: (1 − bHE, X), (1 − bHE, SZ), and (1 − bHE, T) from the left to the right panel. The color code and the meaning of the shaded area and solid lines are the same as in Fig. 9.

From Fig. 12, we could conclude that the three bias measurements are, to a first approximation, all very similar. Looking more carefully however we notice subtle differences which can help to better understand the contribution of each term in the HE mass equation. Indeed, comparing bHE, X and bHE, T allows a better understanding of the impact on the mass bias of the derivatives and specifically of the pressure derivative versus the sum of the gas density and temperature derivatives. The first thing to notice from the figure is that the derivative of the pressure profiles plays a decisive role in increasing the scatter of the distribution especially at R500. A clump manifests its impact more strongly on the derivative of the pressure profile rather than on the sum of the derivatives of the gas density profile and the temperature profile. Indeed, while for the entire sample bHE, X(R500) has the lowest standard deviation, σ(bHE, X) = 0.09, the other two biases have σ(bHE, SZ) = 0.11 and σ(bHE, T) = 0.12 (see the first row of Table 2 for the first two bias expressions). It also appears that the distribution of (1 − bHE, X) is quite symmetric with respect to the median behavior in the entire radial range. Specifically, at R500, bHE, X has a low value of skewness (0.20). Alternatively, the other two biases, bHE, SZ and bHE, T, at the same radius have skewness values respectively equal to 0.73 and 0.90, indicating a (small) predominance of the tail towards and beyond the zero bias (or (1 − b) = 1) over the other tail. For completeness, we report that the distributions of the three biases around R500 are all characterized by a kurtosis parameter, measuring the ratio between the peak and the tails, as expected in a normal distribution.

The differences between bHE, SZ and bHE, T highlight the influence of the multiplicative factor: temperature, rather than the ratio between pressure and gas density. The SZ HE mass bias, (1 − bHE, SZ), shows a general shift towards the value of 1, which is more evident for the IR3D systems and for the central regions. This is consistent with the expectation that the temperature can induce an extra bias in the mass determination in the presence of gas of multiple temperatures. The excess of positive bias, or overestimation of the true mass, for (1 − bHE, SZ) is present at all radii, including R500. There, the number of clusters with (1 − bHE, SZ) > 1 is 50% higher than those with (1 − bHE, X) > 1. At the same radius, but on the other side of the bias distribution, we find that about 10% of the systems have (1 − bHE, T) < 0.75 while only 6% show the same amount of bias using the SZ formulation.

As a summary, we can conclude that the X-ray and the SZ mass biases are substantially providing the same answer but the derivative of the pressure in bHE, SZ can induce a higher scatter and using P/ρ instead of T reduces the overall bias. Even though the last characteristic is clearly desirable, the larger scatter and the higher skewness value make MHE, SZ less appealing because modeling its final distribution is more difficult.

### 7.2. Bias and 3D clumpiness factor

In Fig. 13 we relate the X-ray-mass and the SZ-mass bias to the clumpiness factor. All quantities are computed in 3D and around R500. In the plot we draw with empty circles the points associated with poor analytic fits of the gas profiles (see Sect. 4.3.4 for details). Hereafter, we refer to the remaining clusters (filled circles) as “well-fitted” clusters. Fig. 13. Mass bias as a function of clumpiness. All quantities are considered at R500. Empty and filled circles refer to clusters whose thermodynamical profiles are either poorly or well fitted with the assumed analytic function, respectively (Sect. 4.3.4). The R3D objects are show in brown and the IR3D clusters in olive green.

The most interesting feature is that most of the clusters with positive bias, (1 − bHE, X) ≥ 1 or (1 − bHE, SZ) ≥ 1, are associated with poorly fitted systems and that the only well-fitted system with high positive bias, (1 − bHE, X) ≥ 1.10 and (1 − bHE, SZ) > 1.20, is characterized by a high clumpiness value, 𝒞 >  1.3. The poorly fitted objects cover all values of the bias; indeed, their median bias is similar to that of the well-fitted systems, even though the former have a much higher dispersion, as we can see by comparing the results of the well-fitted clusters with those of the entire sample reported in Table 2.

If we divide the objects according to their clumpiness level, we find that the median bias for less clumped systems is closer to 1 than the mass bias of the most clumped systems13. Specifically, the median of both biases moves from −9% with a dispersion of 6% in the case of 𝒞 <  1.1 (third row in Table 2) to −15% with a dispersion of 11–13% in the case of 𝒞 >  1.2. As expected from Fig. 9, the objects that were classified regular in 3D present a low clumpiness level and have a median bias value which is approximately 10–15% lower than the irregular systems (fourth and fifth row of Table 2). These results are oriented towards the expected conclusions: most regular systems (and generically less clumpy objects) are modestly biased and their distribution has a 20–25% lower dispersion. Most of the poorly fitted systems are classified as irregular.

From these results one expects that the bias and the clumpiness values should be tightly related, but instead the two quantities exhibit a low level of correlation (Table 3 for the quantities measured at R500). The Spearman correlation coefficient is indeed corr(1 − bHE, 𝒞)≲0.20 for the entire sample and corr(1 − bHE, 𝒞)≲0.30 for the well-fitted subsample. Independent simulations have indeed shown that even after removing the 50% densest cells at all cluster radii, the ratio of non-thermal to thermal pressure support is nearly unchanged, suggesting that the HE bias is not simply related to high clumping factor values (Angelinelli et al. 2019). It seems therefore difficult to correct the bias using information exclusively from the clumpiness. Nevertheless, we still attempt to extract a correction to (1 − bHE) by looking for a linear relation between (1 − bHE) and 𝒞. By subtracting the best-fitting line from the individual bias, we succeed in obtaining a median value for the corrected bias very close to one, but we do not find any gain in the standard deviation of the distribution of the corrected biases. Since the bias has a weak correlation with the clumpiness, we search for another parameter among those investigated that could improve the result when combined with the clumpiness.

The only promising parameter that we find is the asymptotic external slope of the gas profiles (see Table 3). Precisely, in Fig. 14 we relate (1 − bHE, SZ) to the value of βP (the asymptotic slope of the pressure profile as in Eq. (7)) and (1 − bHE, X) to the value of 𝒟 = 3β + ϵ/2 (asymptotic slope of the analytic gas density profile of Eq. (5)). As in Fig. 13, we divide the clusters into poorly and well-fitted objects (empty and filled circles) and we further divide the last class into three bins of clumpiness: the 25% with the lowest clumpiness are shown in magenta, the 25% with the highest clumpiness are in cyan, and those in between are plotted in black. The vertical lines indicate the value of the slopes that contain three-quarters of all well-fitted clusters and are equal to 5.5 for the gas density and 6 for the pressure slopes. The density and pressure asymptotic slopes are rarely below a value of 2 and 3, respectively. Fig. 14. Mass biases (1 − bHE, X) and (1 − bHE, SZ) are shown as a function of the slope of the gas density, 𝒟, and of the pressure profile, βP. Empty and filled circles have the same meaning as in Fig. 13. Magenta points are well-fitted clusters with the lowest clumpiness values, the cyan points show those with the highest clumpiness values. The vertical lines represent the value 𝒟 = 5.5 and β = 6 which are approximately the 75th percentiles of the values of the two slopes. All quantities are measured in 3D so each cluster appears only once.

The large majority of clusters with slopes lower than these values or higher than the 75th percentile (shown by the vertical lines) are typically characterized either by high clumpiness values (cyan points) or by high NRMS linked to the fit of the gas profile (empty points). These outliers are also responsible for increasing the dispersion of the bias distributions.

It is remarkable how the SZ bias for βP <  6 has a clear separation between low (magenta) and high (cyan) clumped objects, and overall, both biases have a medium level of correlation with the asymptotic slopes (Table 3): corr(1 − bHE, SZ, βP) ≈ corr(bHE, X, 𝒟)≈0.35 − 0.45. In light of these results, we try to correct the mass bias by using the best-fitting plane of bias, slope, and clumpiness values. The fitting procedure to determine the parameters of the plane was restricted only to the well-fitted clusters. However, we find that correcting all data points (both well and poorly fitted) following this expression for (1 − bHE, SZ): (14)

and equivalently for (1 − bHE, X): (15)

decreases the standard deviations by about 10% (5% for bHE, X). In addition, the skewness of the SZ bias drops from 0.73 to 0.18 indicating that both mass biases now have Gaussian distributions (see fifth row in Table 2).

### 7.3. Bias and 2D gas inhomogeneities

In the previous section we studied the HE bias in 3D and attempted to connect it to 3D properties of the ICM. In this section, we follow a more observation-oriented approach by relating the X-ray and SZ mass bias to all quantities defined from the X-ray maps as listed at the beginning of Sect. 4.1. None of these 2D quantities show a strong correlation with the two biases (Table 3). Only mild correlations for the well-fitted clusters and all 2D proxies (σA, σA, R, ε, and MM) are found with respect to (1 − bHE, X) with values in the range −0.26 <  corr <  −0.29. For the entire sample, including poorly fitted objects, the correlation is even weaker. The correlation values are also lower with respect to (1 − bHE, SZ), for which the only notable correlation value is −0.24 in relation to σA, R and for the well-fitted clusters. These numbers are similar to those reported in other studies that search for a connection between morphological parameters and HE bias (e.g. Jeltema et al. 2008; Rasia et al. 2012).

This lack of correlation is not surprising given that already for the 3D clumpiness the correlation is rather weak. As a consequence, any correction to either (1 − bHE, X) or (1 − bHE, SZ) based on a best linear fit between the bias and any of the 2D quantities is not expected to substantially change the values of the dispersion of the bias distributions or their skewness values, similar to the findings of the previous section. Nevertheless, also regarding the 2D analysis, we notice that if we compare the mass-bias statistics of the subsample with the lowest azimuthal scatter, σA <  0.4, to that of the subsample with the highest scatter12, σA >  1, we find a clear trend since the median biases move from about −9% with a dispersion around 8% (second row in the second panel of Table 2) to −13% with a dispersion of 12%. We therefore expect that these estimators can provide some level of improvement in the correction.

We then proceeded to include the information of the azimuthal scatter and of the asymptotic slope of the gas density, 𝒟, for the X-ray mass bias: (16)

At the same time, we correct the SZ mass bias by invoking the asymptotic slope of the pressure profile, βP: (17)

The impact of this correction on the mass bias distribution is listed in the second panel of Table 2 where we also report on the similar gains achieved by substituting the last term in Eqs. (16) and (17) (i.e., 0.0075 × σA) with factors that depend on the other estimators of ICM gas inhomogeneity:

−0.04 + 0.20 × ε (for the ellipticity);

−0.02 + 0.12 × MM (for the MM value);

−0.01 + 0.025 × σA, R (for the azimuthal scatter averaged over the entire radial range).

In all cases, the largest correction for the bias comes from the slopes of the gas profiles, 𝒟 and βP, that on average contribute 5 and 9%, respectively. The azimuthal scatter, ellipticity, MM, and the azimuthal scatter averaged over the entire radial range account for an extra few percent. However, it is thanks to their inclusion that we can efficiently correct the 2D samples in their entirety, including irregular objects or, generically, those with evidence of substructures (σA >  1). The (albeit small) contribution of the 2D gas inhomogeneity estimators is therefore essential to extend the mass determination to a mixed sample of objects (see also Appendix B for the VI class).

The above equations reduce the X-ray and SZ bias scatter by 10–15% and the skewness value of the SZ bias by almost a factor of three, from 0.96 to 0.34 (Table 2). The changes to the distributions can be appreciated in Fig. 16: the improvement on the skewness values clearly changes the red histograms, making them more closely resemble a Gaussian distribution. Fig. 15. Mass bias as a function of azimuthal scatter averaged over the radial range explored, σA, R. Empty and filled circles have the same meaning as in Fig. 13. Magenta and cyan points highlight the well-fitted clusters with the lowest and highest asymptotic slopes, respectively, of the gas density profile (top panels) and of the pressure (bottom panels). Each point represents a map. Fig. 16. Distribution of the mass biases, (1 − bHE, X), on the left, and (1 − bHE, SZ), on the right, before (top panels) and after (bottom panels) the corrections expressed respectively in Eqs. (16) and (17). The empty histograms show the overall distribution of all the 175 clusters, the filled histograms are restricted to the 97 well-fitted objects. The parameters characterizing the histograms are reported in Table 2.

## 8. Conclusions

This work characterizes the inhomogeneities present in the ICM as measurable from X-ray observations, and links them to the 3D clumpiness level and to the bias of the total cluster mass derived under the assumption of hydrostatic equilibrium. We analyze an extended set of simulated galaxy clusters taken at z = 0 from the Three Hundred Project (Cui et al. 2018). The simulations were performed with the GADGET-X code, which includes an improved formulation of SPH, with respect to GADGET-2 (Springel 2005), and thus promotes the mixing of gas phases with different entropy levels (Beck et al. 2016). The runs incorporate stellar feedback in kinetic form and thermal AGN feedback generated by gas accretion onto super-massive black holes (see Rasia et al. 2015). X-ray images in the soft ([0.5–2] keV) band are produced and processed to extract 2D measurements of gas inhomogeneities. We consider two centers associated with each map: the maximum of the flux (MF) and the center of the ellipse (CE) that best fits the iso-flux contour at around 0.8 R500. The distance of both centers from the theoretical center, the minimum of the potential well, and the ellipticity of the best-fit ellipse are used to provide a first separation of the clusters into different morphological classes: very-regular, VR, regular, R, and intermediate-irregular, IR, objects. From the surface-brightness maps, we compute the azimuthal scatter, σA, over 12 sectors. Two measurements of this scatter are carried out with respect to the mean and the median of the surface brightness profiles. We further consider the ratio between the mean and median, MM, as another estimate of the gas inhomogeneities. Our findings can be summarized as follows.

We compare four estimates of the azimuthal scatter that combine (1) the two options for the sectors center and (2) the two choices for the reference profile (mean or median) used to compute the scatter as expressed in Eq. (2): and . We conclude that the most useful option is because the median boosts the signal revealing the presence of gas inhomogeneities and MF is closer to the theoretical center and is easily determined.

The 2D azimuthal scatter grows from 0.2–0.3 for the regular systems to 0.5–0.7 for the intermediate/irregular clusters. On average, the scatter profile does not vary over the radial range investigated ([0.4 − 1.2] R500), although the IR class presents a high dispersion and skewness, indicating that individual objects at fixed radii can have high scatter values. In our sample of about 540 maps, we found that a scatter higher than unity is a strong indication of the presence of one or more substructures and cannot be ascribed only to an elongated gas distribution. Indeed, even the most elliptical cluster among the substructure-free objects in our IR subsample has σA = 0.73.

The 3D clumpiness is lower for regular objects, 𝒞 ≈ 1.05 − 1.08, than for irregular ones, 𝒞 ≈ 1.1 − 1.2, which are also characterized by a much higher scatter. The clumpiness is closely linked to the azimuthal scatter with a correlation coefficient of about 0.6. However, some outliers are present in the overall distribution. We studied two most extreme cases which include objects in which a clump is aligned with the cluster core (leading to high values of clumpiness and low values of scatter) and objects in which substructures have a small projected distance from the map center although they are external to the cluster in 3D. The latter situation is not very frequent since among all maps with σA <  0.5 only 5% have 𝒞 >  1.35.

We consider three expressions for the hydrostatic mass: two used in X-ray and SZ observational works, MHE, X (Eq. (10)) and MHE, SZ (Eq. (9)) respectively, and a third formulation, MHE, T (Eq. (11)), which helps us to separately evaluate the impact of each term of the expression. All three estimates introduce a similar bias which spans from 5-10% around R2500 to 10–15% around R500. Regular clusters (VR or R), less clumped clusters (𝒞 <  1.1), objects with low azimuthal-scatter (σA <  0.4), and systems with well-behaved gas profiles tend to have a slightly lower mass bias and a reduced scatter. When using the SZ formulation for the HE mass, we find a broadening of the bias distribution such that the tail corresponding to lowest bias reaches the value of (1 − b) = 1 and even greater values, (1 − b) > 1. For cosmological studies, using the pressure profile is less suitable than using the combination of the gas density and temperature profiles, because despite the lower average bias, the distribution of (1 − bHE, SZ) is broader and more skewed.

Even though the average mass bias of less clumped systems or of the objects with smaller azimuthal scatter is lower than the average over the entire sample, the mass bias is not strongly correlated with either parameter. Therefore, they cannot be efficiently used to reduce the dispersion of the mass bias of the entire sample. Adding extra information such as the external slope of either the gas density or the pressure profile diminishes the dispersion by 10% and reduces the skewness of the bias distribution. This essentially improves the accuracy of a Gaussian model to describe the distribution of the mass bias and possibly correct for it. These corrections are suitable not only for regular objects but also for irregular systems and therefore should be sought to exploit large sample sizes.

Modern state-of-the-art numerical models such as those analyzed here produce a realistic description of the gas properties and of the clumpy structure of the ICM. The large number of massive clusters analyzed in this paper allows us to robustly determine that, although the X-ray images cannot provide a compelling correction to the mass bias, it is still possible to use information from X-ray and SZ data to obtain a mass bias distribution which can be appropriately modeled by a Gaussian distribution. This result will be useful for measurements which combine X-ray observations (such as those from XMM-Newton, Chandra, or e-Rosita) and SZ measurements (such as those from Planck, Bolocam, NIKA–2, MUSTANG–2,SPT, or ACT) as recently investigated in Shitanishi et al. (2018) by combining Chandra and Bolocam, in the X-COP collaboration (Ghirardini et al. 2019) using XMM-Newton and Planck data, and in Ruppin et al. (2017) with the combined analysis of NIKA, Planck, and XMM-Newton data. In cases of deep data, more sophisticated algorithms can be used to distinguish between clumpy and linear inhomogeneities (Bourdin et al. 2015; Vafaei Sadr et al. 2018) and to provide indications of inhomogeneities present in SZ maps (Baldi et al. 2019). As a caveat, we note that in this paper we estimate the mass bias by looking directly at the 3D profiles as derived from the simulated clusters since it has been demonstrated that the deprojection of X-ray profiles works in a satisfactory manner. It remains to be proven that the same applies to the pressure profiles as derived from SZ observations, which are typically characterized by a poorer spatial resolution (especially concerning Planck’s observations). This remains to be explored further in a future study.

2

The MultiDark simulations are publicly available at the https://www.cosmosim.org database.

3

We refer to the mass MΔ as the mass of the sphere of radius RΔ with density Δ times the critical density of the universe at that redshift. Virial masses are defined for Δ = 98 following Bryan & Norman (1998). In the rest of the paper we mostly consider Δ = 500.

5

There are multiple levels of low-resolution dark-matter particles, each level having its particle mass ∼10 times more massive than its inner level. For our purpose, we consider all low-resolution particles as contaminant.

7

The density of star formation is approximately equal to 1.95 × 10−25 g cm−3 or 2.88 × 106 kpc−3M.

8

With the exception of the three most massive clusters, which in any case are classified as very irregular in Sect. 3.1 and are therefore not part of the sample analyzed in the main paper.

9

We note that, in a similar way, objects with a very perturbed X-ray morphology (e.g., with ongoing disruptive mergers and/or with large deviations from spherical symmetries) are also discarded from observational analysis of the hydrostatic mass bias in galaxy clusters, or are subject to ad-hoc analysis procedures (e.g. Ghirardini et al. 2018).

10

Here we opt to show the clumpiness (a 3D quantity) for the clusters shown in Figs. 1 and 4, chosen to represent the 2D classifications of VR, R, and IR systems. According to the 3D classification, the first two objects are part of the R3D class, while the latter is part of the IR3D class.

11

To ease the comparison with other numerical works that use the spectroscopic-like temperature in the HE mass derivation rather than the mass-weighted temperature, we computed the ratio of these two temperatures at R500. On average an offset of 10% is found leading to a HE mass bias of 20%. The mismatch between the two temperatures and its impact on the HE mass bias confirms previous results, starting from Rasia et al. (2006), where it was discussed for the first time, to the recent paper by Pearce et al. (2019).

12

We note that the correlation coefficient between clumpiness and azimuthal scatter is not biased in a particular way from the mis-centering between MF and CE. Nevertheless, we notice that when we select objects with DMF <  2 pixels the relation becomes steeper: 𝒞 = 1.0 + 0.3 × σA.

13

Both subsamples include about one third of the points of the entire sample.

## Acknowledgments

We thank the Referee for comments and suggestions that improved the presentation of the results, Frazer Pearce, Franco Vazza for useful discussions and insightful comments, and Volker Springel for making the GADGET3 Code available. S.A. thanks the Observatory and the University of Trieste for warm hospitality during his visit. This work is part of The Three Hundred collaboration. The simulations used in this paper have been performed in the MareNostrum Supercomputer at the Barcelona Supercomputing Center, thanks to CPU time granted by the Red Española de Supercomputacíon. As part of the Three Hundred project, his work has received financial support from the European Union’s Horizon 2020 Research and Innovation programme under the Marie Sklodowskaw-Curie grant agreement number 734374, the LACEGAL project. We further acknowledge financial contribution from the contracts/grants: AYA2015-63810-P (W.C. and G.Y.); ASI 2015-046-R.0 (S.E.); ASI-INAF n.2017-14-H.0 (E.R., S.B., S.E.); PGC2018-094975-B-C21 (G.Y.); PRIN-MIUR 2015W7KAWC (S.B.); INFN INDARK grant (E.R., S.B.); the European Research Council under grant number 670193 (W.C.); DFG Cluster of Excellence “Origins” (K.D.); DFG project BI 1699/2-1 (V.B.); the sabbatical program for PhD student of Iran Ministry of Science, Research and Technology (S.A.). Part of this work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number Ts 17/2–1. As requested by the Three Hundred policy, we state the contribution by each author to this paper: S. A.analysed the data and produced the plots; E.R. developed and supervised the project; E.R. wrote the paper with support from V.B.; E.R., V.B., S.B., contributed to the interpretation of the results; M.D.P, S.E., F.P., S.M.S.M provided critical feedback; W.C., K.D., G.M., G.Y. contributed to the simulated data set. The author list is in alphabetic order with the exclusion of the first four authors.

## References

1. Allen, S. W., Rapetti, D. A., Schmidt, R. W., et al. 2008, MNRAS, 383, 879 [NASA ADS] [CrossRef] [Google Scholar]
2. Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409 [NASA ADS] [CrossRef] [Google Scholar]
3. Ameglio, S., Borgani, S., Pierpaoli, E., et al. 2009, MNRAS, 394, 479 [NASA ADS] [CrossRef] [Google Scholar]
4. Angelinelli, M., Vazza, F., Giocoli, C., et al. 2019, MNRAS, submitted [arXiv:1905.04896] [Google Scholar]
5. Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
6. Arthur, J., Pearce, F. R., Gray, M. E., et al. 2019, MNRAS, 484, 3968 [NASA ADS] [CrossRef] [Google Scholar]
7. Avestruz, C., Lau, E. T., Nagai, D., & Vikhlinin, A. 2014, ApJ, 791, 117 [NASA ADS] [CrossRef] [Google Scholar]
8. Baldi, A. S., Bourdin, H., Mazzotta, P., et al. 2019, A&A, 630, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
9. Barnes, D. J., Kay, S. T., Bahé, Y. M., et al. 2017a, MNRAS, 471, 1088 [NASA ADS] [CrossRef] [Google Scholar]
10. Barnes, D. J., Kay, S. T., Henson, M. A., et al. 2017b, MNRAS, 465, 213 [NASA ADS] [CrossRef] [Google Scholar]
11. Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012, ApJ, 758, 74 [NASA ADS] [CrossRef] [Google Scholar]
12. Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2015, ApJ, 806, 43 [NASA ADS] [CrossRef] [Google Scholar]
13. Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110 [Google Scholar]
14. Biffi, V., & Valdarnini, R. 2015, MNRAS, 446, 2802 [NASA ADS] [CrossRef] [Google Scholar]
15. Biffi, V., Dolag, K., & Böhringer, H. 2013, MNRAS, 428, 1395 [NASA ADS] [CrossRef] [Google Scholar]
16. Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112 [NASA ADS] [CrossRef] [Google Scholar]
17. Biffi, V., Planelles, S., Borgani, S., et al. 2017, MNRAS, 468, 531 [NASA ADS] [CrossRef] [Google Scholar]
18. Biffi, V., Planelles, S., Borgani, S., et al. 2018, MNRAS, 476, 2689 [NASA ADS] [CrossRef] [Google Scholar]
19. Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55 [NASA ADS] [CrossRef] [Google Scholar]
20. Bourdin, H., Mazzotta, P., & Rasia, E. 2015, ApJ, 815, 92 [NASA ADS] [CrossRef] [Google Scholar]
21. Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]
22. Chen, H., Avestruz, C., Kravtsov, A. V., Lau, E. T., & Nagai, D. 2019, MNRAS, 490, 2380 [Google Scholar]
23. Churazov, E., Sazonov, S., Sunyaev, R., et al. 2005, MNRAS, 363, L91 [NASA ADS] [Google Scholar]
24. Cialone, G., De Petris, M., Sembolini, F., et al. 2018, MNRAS, 477, 139 [NASA ADS] [CrossRef] [Google Scholar]
25. Clerc, N., Cucchetti, E., Pointecouteau, E., & Peille, P. 2019, A&A, 629, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
26. Costanzi, M., Rozo, E., Simet, M., et al. 2019, MNRAS, 488, 4779 [NASA ADS] [CrossRef] [Google Scholar]
27. Cucchetti, E., Clerc, N., Pointecouteau, E., Peille, P., & Pajot, F. 2019, A&A, 629, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
28. Cui, W., Power, C., Biffi, V., et al. 2016, MNRAS, 456, 2566 [NASA ADS] [CrossRef] [Google Scholar]
29. Cui, W., Knebe, A., Yepes, G., et al. 2018, MNRAS, 480, 2898 [NASA ADS] [CrossRef] [Google Scholar]
30. Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827 [NASA ADS] [CrossRef] [Google Scholar]
31. Dolag, K., Hansen, F. K., Roncarelli, M., & Moscardini, L. 2005, MNRAS, 363, 29 [NASA ADS] [CrossRef] [Google Scholar]
32. Eckert, D., Roncarelli, M., Ettori, S., et al. 2015, MNRAS, 447, 2198 [NASA ADS] [CrossRef] [Google Scholar]
33. Eckert, D., Gaspari, M., Vazza, F., et al. 2017, ApJ, 843, L29 [NASA ADS] [CrossRef] [Google Scholar]
34. Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
35. Ettori, S., Morandi, A., Tozzi, P., et al. 2009, A&A, 501, 61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
36. Ettori, S., Donnarumma, A., Pointecouteau, E., et al. 2013, Space Sci. Rev., 177, 119 [NASA ADS] [CrossRef] [Google Scholar]
37. Ettori, S., Ghirardini, V., Eckert, D., et al. 2019, A&A, 621, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
38. Evrard, A. E. 1990, ApJ, 363, 349 [NASA ADS] [CrossRef] [Google Scholar]
39. Fang, T., Humphrey, P., & Buote, D. 2009, ApJ, 691, 1648 [NASA ADS] [CrossRef] [Google Scholar]
40. Frank, K. A., Peterson, J. R., Andersson, K., Fabian, A. C., & Sanders, J. S. 2013, ApJ, 764, 46 [NASA ADS] [CrossRef] [Google Scholar]
41. Gardini, A., Rasia, E., Mazzotta, P., et al. 2004, MNRAS, 351, 505 [NASA ADS] [CrossRef] [Google Scholar]
42. Gaspari, M., McDonald, M., Hamer, S. L., et al. 2018, ApJ, 854, 167 [NASA ADS] [CrossRef] [Google Scholar]
43. Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
44. Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
45. Henson, M. A., Barnes, D. J., Kay, S. T., McCarthy, I. G., & Schaye, J. 2017, MNRAS, 465, 3361 [NASA ADS] [CrossRef] [Google Scholar]
46. Hitomi Collaboration (Aharonian, F., et al.) 2016, Nature, 535, 117 [Google Scholar]
47. Huang, S., Katz, N., Davé, R., et al. 2019, MNRAS, 484, 2021 [NASA ADS] [CrossRef] [Google Scholar]
48. Jeltema, T. E., Hallman, E. J., Burns, J. O., & Motl, P. M. 2008, ApJ, 681, 167 [Google Scholar]
49. Khatri, R., & Gaspari, M. 2016, MNRAS, 463, 655 [NASA ADS] [CrossRef] [Google Scholar]
50. Khedekar, S., Churazov, E., Kravtsov, A., et al. 2013, MNRAS, 431, 954 [NASA ADS] [CrossRef] [Google Scholar]
51. Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340 [NASA ADS] [CrossRef] [Google Scholar]
52. Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353 [NASA ADS] [CrossRef] [Google Scholar]
53. Laganá, T. F., Durret, F., & Lopes, P. A. A. 2019, MNRAS, 484, 2807 [NASA ADS] [CrossRef] [Google Scholar]
54. Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129 [NASA ADS] [CrossRef] [Google Scholar]
55. Lopes, P. A. A., Trevisan, M., Laganá, T. F., et al. 2018, MNRAS, 478, 5473 [NASA ADS] [CrossRef] [Google Scholar]
56. Lovisari, L., Forman, W. R., Jones, C., et al. 2017, ApJ, 846, 51 [NASA ADS] [CrossRef] [Google Scholar]
57. Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2014, MNRAS, 440, 2077 [NASA ADS] [CrossRef] [Google Scholar]
58. Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2015, MNRAS, 449, 199 [NASA ADS] [CrossRef] [Google Scholar]
59. Mathiesen, B., Evrard, A. E., & Mohr, J. J. 1999, ApJ, 520, L21 [NASA ADS] [CrossRef] [Google Scholar]
60. Mazzotta, P., Rasia, E., Moscardini, L., & Tormen, G. 2004, MNRAS, 354, 10 [NASA ADS] [CrossRef] [Google Scholar]
61. Meneghetti, M., Rasia, E., Merten, J., et al. 2010, A&A, 514, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
62. Morandi, A., & Cui, W. 2014, MNRAS, 437, 1909 [NASA ADS] [CrossRef] [Google Scholar]
63. Morandi, A., Nagai, D., & Cui, W. 2013, MNRAS, 436, 1123 [NASA ADS] [CrossRef] [Google Scholar]
64. Morandi, A., Sun, M., Mulchaey, J., Nagai, D., & Bonamente, M. 2017, MNRAS, 469, 2423 [NASA ADS] [CrossRef] [Google Scholar]
65. Mostoghiu, R., Knebe, A., Cui, W., et al. 2019, MNRAS, 483, 3390 [CrossRef] [Google Scholar]
66. Nagai, D., & Lau, E. T. 2011, ApJ, 731, L10 [NASA ADS] [CrossRef] [Google Scholar]
67. Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007a, ApJ, 668, 1 [NASA ADS] [CrossRef] [Google Scholar]
68. Nagai, D., Vikhlinin, A., & Kravtsov, A. V. 2007b, ApJ, 655, 98 [NASA ADS] [CrossRef] [Google Scholar]
69. Nelson, K., Rudd, D. H., Shaw, L., & Nagai, D. 2012, ApJ, 751, 121 [NASA ADS] [CrossRef] [Google Scholar]
70. Pearce, F. A., Kay, S. T., Barnes, D. J., Bower, R. G., & Schaller, M. 2019, MNRAS, 2606 [Google Scholar]
71. Piffaretti, R., & Valdarnini, R. 2008, A&A, 491, 71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
72. Planck Collaboration XX. 2014, A&A, 571, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
73. Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
74. Planelles, S., Fabjan, D., Borgani, S., et al. 2017, MNRAS, 467, 3827 [NASA ADS] [CrossRef] [Google Scholar]
75. Power, C., Elahi, P. J., Welker, C., et al. 2019, MNRAS, 2758 [Google Scholar]
76. Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25 [NASA ADS] [CrossRef] [Google Scholar]
77. Rasia, E., Tormen, G., & Moscardini, L. 2004, MNRAS, 351, 237 [NASA ADS] [CrossRef] [Google Scholar]
78. Rasia, E., Ettori, S., Moscardini, L., et al. 2006, MNRAS, 369, 2013 [NASA ADS] [CrossRef] [Google Scholar]
79. Rasia, E., Meneghetti, M., Martino, R., et al. 2012, New J. Phys., 14, 055018 [NASA ADS] [CrossRef] [Google Scholar]
80. Rasia, E., Meneghetti, M., & Ettori, S. 2013, Astron. Rev., 8, 40 [Google Scholar]
81. Rasia, E., Lau, E. T., Borgani, S., et al. 2014, ApJ, 791, 96 [NASA ADS] [CrossRef] [Google Scholar]
82. Rasia, E., Borgani, S., Murante, G., et al. 2015, ApJ, 813, L17 [NASA ADS] [CrossRef] [Google Scholar]
83. Richardson, M. L. A., Scannapieco, E., Devriendt, J., et al. 2016, ApJ, 825, 83 [NASA ADS] [CrossRef] [Google Scholar]
84. Roncarelli, M., Ettori, S., Borgani, S., et al. 2013, MNRAS, 432, 3030 [NASA ADS] [CrossRef] [Google Scholar]
85. Roncarelli, M., Gaspari, M., Ettori, S., et al. 2018, A&A, 618, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
86. Rossetti, M., Gastaldello, F., Ferioli, G., et al. 2016, MNRAS, 457, 4515 [Google Scholar]
87. Ruppin, F., Adam, R., Comis, B., et al. 2017, A&A, 597, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
88. Ruppin, F., Mayet, F., Pratt, G. W., et al. 2018, A&A, 615, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
89. Salvati, L., Douspis, M., Ritz, A., Aghanim, N., & Babul, A. 2019, A&A, 626, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
90. Sembolini, F., Yepes, G., De Petris, M., et al. 2013, MNRAS, 429, 323 [NASA ADS] [CrossRef] [Google Scholar]
91. Sembolini, F., Elahi, P. J., Pearce, F. R., et al. 2016, MNRAS, 459, 2973 [NASA ADS] [CrossRef] [Google Scholar]
92. Shi, X., Komatsu, E., Nelson, K., & Nagai, D. 2015, MNRAS, 448, 1020 [NASA ADS] [CrossRef] [Google Scholar]
93. Shitanishi, J. A., Pierpaoli, E., Sayers, J., et al. 2018, MNRAS, 481, 749 [NASA ADS] [CrossRef] [Google Scholar]
94. Simionescu, A., ZuHone, J., Zhuravleva, I., et al. 2019, Space Sci. Rev., 215, 24 [NASA ADS] [CrossRef] [Google Scholar]
95. Springel, V. 2005, MNRAS, 364, 1105 [NASA ADS] [CrossRef] [Google Scholar]
96. Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 [Google Scholar]
97. Steinborn, L. K., Dolag, K., Hirschmann, M., Prieto, M. A., & Remus, R.-S. 2015, MNRAS, 448, 1504 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
98. Suto, D., Kawahara, H., Kitayama, T., et al. 2013, ApJ, 767, 79 [NASA ADS] [CrossRef] [Google Scholar]
99. Tornatore, L., Borgani, S., Dolag, K., & Matteucci, F. 2007, MNRAS, 382, 1050 [NASA ADS] [CrossRef] [Google Scholar]
100. Truong, N., Rasia, E., Mazzotta, P., et al. 2018, MNRAS, 474, 4089 [NASA ADS] [CrossRef] [Google Scholar]
101. Truong, N., Rasia, E., Biffi, V., et al. 2019, MNRAS, 484, 2896 [NASA ADS] [CrossRef] [Google Scholar]
102. Urban, O., Simionescu, A., Werner, N., et al. 2014, MNRAS, 437, 3939 [NASA ADS] [CrossRef] [Google Scholar]
103. Vafaei Sadr, A., Movahed, S. M. S., Farhang, M., Ringeval, C., & Bouchet, F. R. 2018, MNRAS, 475, 1010 [NASA ADS] [CrossRef] [Google Scholar]
104. van Weeren, R. J., de Gasperin, F., Akamatsu, H., et al. 2019, Space Sci. Rev., 215, 16 [NASA ADS] [CrossRef] [Google Scholar]
105. Vazza, F., Brunetti, G., Kritsuk, A., et al. 2009, A&A, 504, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
106. Vazza, F., Roncarelli, M., Ettori, S., & Dolag, K. 2011, MNRAS, 413, 2305 [NASA ADS] [CrossRef] [Google Scholar]
107. Vazza, F., Eckert, D., Simionescu, A., Brüggen, M., & Ettori, S. 2013, MNRAS, 429, 799 [NASA ADS] [CrossRef] [Google Scholar]
108. Vazza, F., Wittor, D., Brüggen, M., & Gheller, C. 2016, Galaxies, 4, 60 [NASA ADS] [CrossRef] [Google Scholar]
109. Vazza, F., Angelinelli, M., Jones, T. W., et al. 2018, MNRAS, 481, L120 [NASA ADS] [CrossRef] [Google Scholar]
110. Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691 [NASA ADS] [CrossRef] [Google Scholar]
111. Vikhlinin, A., Kravtsov, A. V., Burenin, R. A., et al. 2009, ApJ, 692, 1060 [NASA ADS] [CrossRef] [Google Scholar]
112. Voit, G. M. 2005, Rev. Mod. Phys., 77, 207 [NASA ADS] [CrossRef] [Google Scholar]
113. Walker, S. A., Fabian, A. C., Sanders, J. S., & George, M. R. 2012, MNRAS, 424, 1826 [NASA ADS] [CrossRef] [Google Scholar]
114. Walker, S., Simionescu, A., Nagai, D., et al. 2019, Space Sci. Rev., 215, 7 [NASA ADS] [CrossRef] [Google Scholar]
115. Wang, Y., Pearce, F., Knebe, A., et al. 2018, ApJ, 868, 130 [NASA ADS] [CrossRef] [Google Scholar]
116. Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99 [Google Scholar]
117. Zhuravleva, I., Churazov, E., Kravtsov, A., et al. 2013, MNRAS, 428, 3274 [NASA ADS] [CrossRef] [Google Scholar]
118. Zinger, E., Dekel, A., Birnboim, Y., et al. 2018, MNRAS, 476, 56 [NASA ADS] [CrossRef] [Google Scholar]
119. ZuHone, J. A., Miller, E. D., Bulbul, E., & Zhuravleva, I. 2018, ApJ, 853, 180 [NASA ADS] [CrossRef] [Google Scholar]

## Appendix A: Dependence of the mass bias from the cluster mass

The sample studied in the paper is mostly composed of massive objects as suggested by the median value of M500 reported in Table 1 for all the subsamples investigated. Restricting the distribution of masses (shown in Fig. 3) to the VR, R, and IR subsamples, we find that 85% of the systems have M500 >  6 × 1014M. Half of the remaining clusters has a mass below 1014M and the others have a mass between 1 and 2 × 1014M. Owing to this mass distribution, the overall sample is not particularly suited to investigating the dependence of the mass bias on the cluster mass, however we can still discuss here the general trend of our clusters.

In Fig. A.1, we show the distribution of the X-ray and SZ mass bias, (1 − bHE), as function of the total mass of the clusters. Since both measurements are derived in 3D, each point refers to a single cluster rather than to a map. The correlation coefficients of both biases are below 0.15 and are consistent with zero within 2σ. Considering only the well-fitted objects (filled points in the plot) the correlation is reduced to values even lower than 0.10. Looking at the median bias values computed in different mass bins, we find that the bias shifts between 0.88 and 0.92 without any special trend. We therefore conclude that in our sample we do not find any dependence of the mass bias on the cluster mass. This result might be surprising because the expectation is that the least massive systems have a reduced bias because they are typically formed earlier and have more time to relax. Our result could be influenced by the poor representation of small groups in our sample. That said, among the most massive objects (between 5 and 12 × 1014M), which are a complete sample and have a mass range similar to the Planck clusters, we can robustly say that there is no trend between (1 − bHE) and M500, the correlation being below 0.10 for both biases. Fig. A.1. Hydrostatic mass bias vs. the cluster mass at R500. Top panel: X-ray mass bias, (1 − bHE, X). Bottom panel: SZ mass bias, (1 − bHE, SZ). Similarly to Fig. 14 of the main text, well-fitted clusters are represented by filled points and poorly fitted objects by empty points.

## Appendix B: The “very irregular” class

In the analysis presented in the paper, we exclude both the very-irregular (VI) and extremely-irregular (EI) objects. The former class is excluded because the center of the X-ray analysis, used to derive azimuthal scatter, ellipticity, and the MM parameter, could be significantly different from the theoretical center used to derive the quantities linked to the mass bias. The latter class is even more complicated because not only could the theoretical center be at a large distance from the gas centroid, but also because the respective X-ray emission is strongly asymmetric. For these reasons, as mentioned in the main text, we consider the extension of the analysis to the extremely-irregular objects to be fruitless. However, it can be useful to have indications of the mass bias of the VI class and on its correlation with the gas inhomogeneity parameters studied in the paper. In this section, all quantities are measured at R500.

In Table B.1, we report median value, standard deviation, and skewness of the hydrostatic mass bias distribution of the very-irregular clusters. The average and standard deviation are 50% higher in this class with respect to the first three less-disturbed classes (VR, R, and IR, see Table 2). Specifically, the mean of the mass bias grows from 10 to 15% and the standard deviation increases from about 0.10 to 0.15.

Table B.1.

Summary of results on the HE mass bias at R500.

Including the very-irregular objects in the main sample, the correlation between the parameters that describe the gas inhomogeneity and the mass bias are stronger, as we can see by comparing the values in Table B.2 with those in Table 3. By applying the same corrections proposed in the paper, we reduce the bias and generate a more symmetric distribution (meaning that the skewness parameter strongly decreases). The gain on the dispersion is limited to less than 10% and the resulting dispersion is always higher than the maximum dispersion found in the other three classes (see Table 3). We therefore conclude that there is no net gain in adding objects that clearly appear morphologically disturbed, such as those with large distance between the center of the brightest cluster galaxy and the X-ray peak, because they are not in equilibrium. For completeness, we report that among the VI clusters, 18 objects have (1 − b) < 0.70, 3 of which have the bias actually lower than 0.60. In the entire sample of VR, R, IR, VI only 3% of objects have an HE mass lower than the true mass by more than 30%.

Table B.2.

Spearman rank correlation coefficient and the number of standard deviations from the null-hypothesis expected value (in parenthesis) similarly to Table 3 but now including also the VI class.

## Appendix C: The HE mass bias at R2500

In the VR, R, and IR classes, the R2500 radius is on average equal to 0.45 R500 and thus it is within the radial range investigated. In this section, we proceed to compute the mass bias and the gas inhomogeneity parameters, such as clumpiness, residual clumpiness, azimuthal scatter, and MM at R2500, and we further link the mass bias at the same radius with the ellipticity value of the X-ray iso-flux contours drawn at 0.8 R500, the radial average of the azimuthal scatter, and the asymptotic slope of the analytic profiles that best describe the gas density and pressure profiles.

For this analysis, we discard all maps corresponding to cluster whose R2500 radius is smaller than our innermost radial bin, 0.4 R500. This sample reduction led to a total of 467 maps.

On average we find that (1 − bHE, X, 2500)∼(1 − bHE, SZ, 2500) = 0.92, with respective standard deviations of 0.10 and 0.13. These values should be compared with the first line of the 2D sample (“all”) in the second panel of Table 2. The derived X-ray and SZ mass are thus slightly closer to the true mass at R2500 than at R500 as suggested also from Fig. 12 where we notice an increase in the departure from the true mass at increasing radius. At R2500 both X-ray and SZ mass bias distribution have a reduced level of skewness (equal to −0.06 and −0.50) indicating that the two distributions are already quite symmetric in contrast with the previous result Sk(bHE, SZ) = 0.96.

The correlation coefficients between mass bias and measurements of the gas inhomogeneities at R2500 are reported in Table C.1. The values listed here should be compared with those in the second section of Table 3. The correlation is in general weaker. The highest variation relates to the correlation between the mass bias and the asymptotic slope. Previously, bHE, X(R500) exhibits the strongest correlation with 𝒟, while at R2500 the two quantity are completely uncorrelated.

Table C.1.

Spearman rank correlation coefficient and the number of standard deviations from the null-hypothesis expected value (in parenthesis) similarly to Table 3.

## All Tables

Table 1.

Basic properties of the subsamples for the 2D analysis (upper part) and the 3D analysis (lower part).

Table 2.

Summary of HE mass bias result at R500.

Table 3.

Spearman rank correlation coefficient and the number of standard deviations from the null-hypothesis expected value (in parenthesis) computed between the mass biases and different quantities listed in the first column and obtained in the 2D and 3D analyses at R500.

Table B.1.

Summary of results on the HE mass bias at R500.

Table B.2.

Spearman rank correlation coefficient and the number of standard deviations from the null-hypothesis expected value (in parenthesis) similarly to Table 3 but now including also the VI class.

Table C.1.

Spearman rank correlation coefficient and the number of standard deviations from the null-hypothesis expected value (in parenthesis) similarly to Table 3.

## All Figures Fig. 1. X-ray maps of three clusters, which are representative of the first three classes described in Sect. 3.1: Very Regular (left panel), Regular (central panel), and Intermediate/Irregular (right panel). The three objects have a similar mass, M500 ≈ 1 × 1015 M⊙, corresponding to R500 of about 1.5 Mpc (shown as a white circle in the images). The black lines are the log-spaced iso-flux contours smoothed over 4 pixels. In the text Fig. 2. Distribution of DCE (in units of R500) and ε for all maps in the categories: VR (blue triangles), R (red pentagons), IR (green squares), and VI (magenta circles). Top and right panels: abundance of DCE and ε of each class. In the text Fig. 3. Distribution of M500 for the 2D classification (top panel) and the 3D classification (bottom panel). In the text Fig. 4. Upper panels: surface brightness profiles (gray lines) extracted from the 12 sectors centered on MF of the three images of Fig. 1, representative of the classes: VR (left panel), R (central panel), and IR (right panel). The thick solid lines show the mean profile. The median profile is not shown for clarity. Lower panels: azimuthal scatter profiles derived as in Eq. (2) and computed with respect to the median profile and centered on MF, , for clarity simply referred to as σA in the plot. In the text Fig. 5. Upper panels: clumpiness (solid line) and residual clumpiness (dashed line). Lower panels: profiles of the hydrostatic-equilibrium mass bias as in Eq. (12). The quantities (see footnote 8) are measured in 3D for the clusters whose maps are representative of the 2D classification of Fig. 1. In the text Fig. 6. Ratio between the azimuthal scatters computed in Eq. (2) with respect to the median, , and to the mean, and centered on MF. The solid lines represent the median of the ratios in each radial bin over the VR (blue), R (red), and the IR (green) subsamples. The shaded regions comprise the distribution between the 16th and the 84th percentiles. In the text Fig. 7. Ratio between the azimuthal scatters centered on MF, , and centered on CE, and computed in Eq. (2) with respect to the median. The color code and meaning of the shaded area and solid lines are the same as in Fig. 6. In the text Fig. 8. Azimuthal scatter computed as in Eq. (2) with respect to the median and centered in MF, . The color code and the meaning of the shaded area and solid lines are the same as in Fig. 6. In the text Fig. 9. Clumpiness profiles for the R3D class in red and for the IR3D class in olive green. The solid line refers to the median profile and the shaded area shows the values between the 16th and 84th percentiles of the 𝒞 distribution. In the text Fig. 10. Distribution of clumpiness vs. azimuthal scatter for the IR systems and for all radial bins. Only points with 𝒞 <  1.5 and σA <  2 are considered. For clarity, the distribution is shown in small bins of the two quantities. The colors indicate the number of points per bin. The color scale is saturated at ten. In the text Fig. 11. Clumpiness vs. azimuthal scatter computed around R500. Top panel: each σA refers to a separate map. The color code reflects the 2D classification: VR in blue, R in red, and IR in green. Bottom panel: each cluster is represented only once and the azimuthal scatter is averaged over all its considered projections, ⟨σA⟩. The R3D objects are shown in brown, and the IR3D clusters in olive green. In the text Fig. 12. Median profile of the HE mass bias: (1 − bHE, X), (1 − bHE, SZ), and (1 − bHE, T) from the left to the right panel. The color code and the meaning of the shaded area and solid lines are the same as in Fig. 9. In the text Fig. 13. Mass bias as a function of clumpiness. All quantities are considered at R500. Empty and filled circles refer to clusters whose thermodynamical profiles are either poorly or well fitted with the assumed analytic function, respectively (Sect. 4.3.4). The R3D objects are show in brown and the IR3D clusters in olive green. In the text Fig. 14. Mass biases (1 − bHE, X) and (1 − bHE, SZ) are shown as a function of the slope of the gas density, 𝒟, and of the pressure profile, βP. Empty and filled circles have the same meaning as in Fig. 13. Magenta points are well-fitted clusters with the lowest clumpiness values, the cyan points show those with the highest clumpiness values. The vertical lines represent the value 𝒟 = 5.5 and β = 6 which are approximately the 75th percentiles of the values of the two slopes. All quantities are measured in 3D so each cluster appears only once. In the text Fig. 15. Mass bias as a function of azimuthal scatter averaged over the radial range explored, σA, R. Empty and filled circles have the same meaning as in Fig. 13. Magenta and cyan points highlight the well-fitted clusters with the lowest and highest asymptotic slopes, respectively, of the gas density profile (top panels) and of the pressure (bottom panels). Each point represents a map. In the text Fig. 16. Distribution of the mass biases, (1 − bHE, X), on the left, and (1 − bHE, SZ), on the right, before (top panels) and after (bottom panels) the corrections expressed respectively in Eqs. (16) and (17). The empty histograms show the overall distribution of all the 175 clusters, the filled histograms are restricted to the 97 well-fitted objects. The parameters characterizing the histograms are reported in Table 2. In the text Fig. A.1. Hydrostatic mass bias vs. the cluster mass at R500. Top panel: X-ray mass bias, (1 − bHE, X). Bottom panel: SZ mass bias, (1 − bHE, SZ). Similarly to Fig. 14 of the main text, well-fitted clusters are represented by filled points and poorly fitted objects by empty points. In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
{"ft_lang_label":"__label__en","ft_lang_prob":0.8870734,"math_prob":0.9190379,"size":125926,"snap":"2022-27-2022-33","text_gpt3_token_len":32220,"char_repetition_ratio":0.18258128,"word_repetition_ratio":0.14120752,"special_character_ratio":0.2551578,"punctuation_ratio":0.16396059,"nsfw_num_words":6,"has_unicode_error":false,"math_prob_llama3":0.96824044,"pos_list":[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158],"im_url_duplicate_count":[null,6,null,6,null,6,null,3,null,6,null,6,null,3,null,3,null,3,null,3,null,3,null,3,null,3,null,6,null,3,null,3,null,3,null,3,null,3,null,3,null,3,null,3,null,6,null,6,null,6,null,3,null,3,null,3,null,6,null,6,null,6,null,3,null,3,null,3,null,3,null,3,null,6,null,6,null,3,null,6,null,6,null,3,null,6,null,6,null,6,null,6,null,3,null,3,null,3,null,3,null,6,null,6,null,3,null,3,null,3,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null,6,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2022-07-07T11:16:45Z\",\"WARC-Record-ID\":\"<urn:uuid:c7ee0089-c662-4b9e-9952-69fcf42dbe70>\",\"Content-Length\":\"413790\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:a4ff137a-a8eb-4225-a8a0-f5da1536b258>\",\"WARC-Concurrent-To\":\"<urn:uuid:2116b368-b9d6-4f8e-8d91-66de86a2f2ed>\",\"WARC-IP-Address\":\"167.114.155.65\",\"WARC-Target-URI\":\"https://www.aanda.org/articles/aa/full_html/2020/02/aa36742-19/aa36742-19.html\",\"WARC-Payload-Digest\":\"sha1:VWT3CWJXEHHEDI2NZBCVI76D4SIVANPV\",\"WARC-Block-Digest\":\"sha1:2LRELTMSCWLVGHETGRMDLE3HPL54ZGMM\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2022/CC-MAIN-2022-27/CC-MAIN-2022-27_segments_1656104690785.95_warc_CC-MAIN-20220707093848-20220707123848-00522.warc.gz\"}"}