Issue 
A&A
Volume 634, February 2020



Article Number  A113  
Number of page(s)  21  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201936742  
Published online  19 February 2020 
The Three Hundred Project: Correcting for the hydrostaticequilibrium mass bias in Xray and SZ surveys
^{1}
Department of Physics, Shahid Beheshti University, G.C., Evin, Tehran 19839, Iran
email: s_ansarifard@sbu.ac.ir
^{2}
Dipartimento di Fisica, Sezione di Astronomia, Università di Trieste, via Tiepolo 11, 34143 Trieste, Italy
^{3}
INAF–Osservatorio Astronomico Trieste, via Tiepolo 11, 34123 Trieste, Italy
email: elena.rasia@inaf.it
^{4}
Institute of Fundamental Physics of the Universe, via Beirut 2, 34151 Grignano, Trieste, Italy
^{5}
HarvardSmithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA
^{6}
UniversitätsSternwarte München, Scheinerstr. 1, 81679 München, Germany
^{7}
INFN, Instituto Nazionale di Fisica Nucleare, Trieste, Italy
^{8}
Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK
^{9}
Dipartimento di Fisica, Sapienza Università di Roma, p.le Aldo Moro 5, 00185 Rome, Italy
^{10}
Max Planck Institut for Astrophysics, 85748 Garching, Germany
^{11}
INAF, Osservatorio di Astrofisica e Scienza dello Spazio, via Pietro Gobetti 93/3, 40129 Bologna, Italy
^{12}
INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy
^{13}
Departamento de Física Teórica M8, Universidad Autónoma de Madrid, Cantoblanco 28049, Madrid, Spain
^{14}
Centro de Investigación Avanzada en Física Fundamental (CIAFF), Facultad de Ciencias, Universidad Autónoma de Madrid, Cantoblanco 28049, Madrid, Spain
Received:
19
September
2019
Accepted:
18
November
2019
Accurate and precise measurement of the masses of galaxy clusters is key to deriving robust constraints on cosmological parameters. However, increasing evidence from observations confirms that Xray masses obtained under the assumption of hydrostatic equilibrium might be underestimated, as previously predicted by cosmological simulations. We analyze more than 300 simulated massive clusters from the Three Hundred Project, and investigate the connection between mass bias and several diagnostics extracted from synthetic Xray images of these simulated clusters. We find that the azimuthal scatter measured in 12 sectors of the Xray flux maps is a statistically significant indication of the presence of an intrinsic (i.e., 3D) clumpy gas distribution. We verify that a robust correction to the hydrostatic mass bias can be inferred when estimates of the gas inhomogeneity from Xray maps (such as the azimuthal scatter or the gas ellipticity) are combined with the asymptotic external slope of the gas density or pressure profiles, which can be respectively derived from Xray and millimeter (SunyaevZeldovich effect) observations. We also obtain that mass measurements based on either gas density and temperature or gas density and pressure result in similar distributions of the mass bias. In both cases, we provide corrections that help reduce both the dispersion and skewness of the mass bias distribution. These are effective even when irregular clusters are included leading to interesting implications for the modeling and correction of hydrostatic mass bias in cosmological analyses of current and future Xray and SZ cluster surveys.
Key words: galaxies: clusters: general / galaxies: clusters: intracluster medium / Xrays: galaxies: clusters / largescale structure of Universe / methods: numerical
© ESO 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 socalled mass bias, defined as an underestimation of the total mass, was immediately suspected to be the root of the problem. Even though in the Planckspecific 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 Xray 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 Xray, 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 intracluster 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 overdense 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 Fermilike mechanisms; see van Weeren et al. 2019 for a recent review) and from the observed correlation between Xray surface brightness fluctuations and radio power Eckert et al. 2017) which suggests a dynamical link between perturbed Xray morphologies and residual gas motions. Unfortunately, for the direct ICM velocity measurements we need to wait for nextgeneration Xray 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 multitemperature nature of the ICM is not recognized from the residuals of the spectroscopic fitting analysis. Indeed, Xray CCD detectors onboard Chandra, XMMNewton, 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 Xray or SZ signal on kiloparsec scales associated with exquisite spectralimaging capabilities. These two requirements are not met by existing SZ telescopes, while with current Xray 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 Xray instruments, in this paper we investigate how these estimates can be used to statistically correct the mass bias at R_{500}. 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 XMMNewton 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 Xray surfacebrightness maps that are not convolved with any instrument response. We remark that the quantities considered have already been derived in Xray 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 Xray 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 R_{500} 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 Project^{1} 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 darkmatter(DM)only simulation consists of a periodic cube of 1.5 Gpc in comoving length containing 3840^{3} DM particles. As the name suggests, this simulation assumes the bestfitting 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 resimulated at higher resolution are identified at z = 0 as the volume centered on the selected massive clusters (all with virial masses^{3} greater than 1.2 × 10^{15} M_{⊙}) 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 volumelimited masscomplete sample for all objects with M_{500} > 6.5 × 10^{14} M_{⊙}. 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 GINNUNGAGAP^{4} 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 largescale tidal fields without raising the computational cost. The highresolution darkmatter particle mass is equal to m_{DM} = 1.9 × 10^{9} M_{⊙}, while the initial gas mass is equal to 3.5 × 10^{8} M_{⊙}. 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 onethousandth 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 resimulated with three hydrodynamical codes: GADGETX (Rasia et al. 2015), GADGETMUSIC (Sembolini et al. 2013), and GIZMOSIMBA (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 GADGETX sample.
The Tree–Particle–Mesh gravity solver of GADGETX corresponds to that of the GADGET3 code, which is an updated and more efficient version of the GADGET2 code (Springel 2005). GADGETX includes an improved SPH scheme with artificial thermal diffusion, timedependent artificial viscosity, highorder Wendland C4 interpolating kernel and wakeup scheme as described in Beck et al. (2016). Radiative gascooling 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 supermassive 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 R_{500}.
2.1. Generation of maps
For the current analysis, in each Lagrangian region we consider the mostmassive clusters at z = 0 without any lowresolution particles^{5} within R_{100}, which is close to the virial radius for our cosmology. Whenever the main object had at least one lowresolution particle within R_{100} we considered the second most massive cluster. Seven regions have no available objects with M_{500} > 3 × 10^{13} M_{⊙}, which is the mass limit that we impose. For each cluster, we produced with the code Smac^{6} (Dolag et al. 2005) three Xray surface brightness maps in the softenergy 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 formation^{7} – and emit in the Xray band – meaning with temperature above 10^{6} 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 R_{500} radii inscribed in the map^{8}, 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 R_{500}). 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.
Xray 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, M_{500} ≈ 1 × 10^{15} M_{⊙}, corresponding to R_{500} of about 1.5 Mpc (shown as a white circle in the images). The black lines are the logspaced isoflux 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 Xray surface brightness maps: the emission peak and the center of the isoflux contours. Both are commonly adopted in Xray 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, D_{MF}, exceeds 0.4 × R_{500}, we proceed to visually inspect the individual map. During this process, we verify whether the Xray 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 D_{MF} in pixels rather than in units of R_{500} 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 Xray 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 isoflux contours of the images drawn at around 0.8 R_{500} 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, D_{CE}, has a broader distribution than for D_{MF}, we measure it in units of R_{500}.
In addition to CE and D_{CE}, we also save the value of the ellipticity of the bestfitting 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 bestfitting 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 results^{9}.
All the other images (precisely, 864 maps corresponding to 288 clusters) and their measured D_{MF}, D_{CE}, 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 Xray 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 massaccretion 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 Xray, 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 Xray 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 Xray 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 hydrostaticequilibrium mass. For this reason, we also introduce a 3D classification.
3.1. Twodimensional classification
We use the computed values of D_{CE}, D_{MF}, 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 bestfitting ellipse of the external isoflux contour, because after looking at some images we recognize that both D_{CE} 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, colorcoded with the classes presented below. In contrast, the third parameter, D_{MF}, 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 D_{MF} 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 D_{MF} as a secondary condition.
Fig. 2.
Distribution of D_{CE} (in units of R_{500}) 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 D_{CE} 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 preselect as part of this class all objects with D_{CE} ≤ 0.052 R_{500} and ε ≤ 0.1. These threshold values correspond to about the 20th percentiles of the corresponding distributions. In addition, we further impose that D_{MF} ≤ 5 pixels (less than 20 kpc), even though D_{MF} is within 1 or 2 pixels for most of the preselected 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 D_{CE} is larger than the imposed limit because it has the roundest isoflux 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 D_{CE} 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 ellipsecenter distance and of the ellipticity are: D_{CE} = 0.035 R_{500} and ε = 0.073.
Regular (R). 102 maps. The second class again includes regular images, but the isoflux contours are allowed to have a small miscentering (D_{CE} ≤ 0.12 R_{500}) and to be more elliptical (ε ≤ 0.15, with three exceptions at ε ∼ 0.22 but D_{CE} ≤ 0.05 R_{500}). In addition, small substructures can be present within R_{500} and the condition on D_{MF} is more relaxed (D_{MF} < 15 pixels – about 60 kpc – even though D_{MF} is within 5 pixels for the majority of the clusters). The following values correspond to the median CE distance and ellipticity: D_{CE} = 0.057 R_{500} and ε = 0.116.
Intermediate/Irregular (IR). 424 maps. Their CE center can have a nonnegligible offset with respect to the minimum of the potential well (but still D_{CE} < 0.2 R_{500}) and the axes of the bestfitting 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 D_{MF} 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 R_{500}. 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 isoflux contours at about R_{500} without major substructures. The median values of distance and ellipticity for the irregular class are: D_{CE} = 0.107 R_{500}, and ε = 0.211.
Very Irregular (VI). 118 maps. All these maps have D_{CE} between 0.2 and 0.5 R_{500}, implying that either the gas distribution is extremely disturbed at larger distances or they have significant substructures that impact the ellipse fitting at about R_{500}. Also in this case, we do not consider any threshold on D_{MF}. The median value of the ellipse center distance is significantly larger, D_{CE} = 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 miscentering 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 Xray 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. Onequarter 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. Threedimensional 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, R_{3D}, and the 3D irregular, IR_{3D}, 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 R_{3D} and IR_{3D} 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] × 10^{14} M_{⊙}, while this is [0.4 − 26.0] × 10^{14} M_{⊙} for the EI sample, where the second most massive clusters has M_{500} = 1.6 × 10^{15} M_{⊙}. 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 × 10^{13} M_{⊙} < M_{500} < 1.4 × 10^{15} M_{⊙}. These clusters are the main objects in their respective Lagrangian regions. Those with M_{500} < 3 × 10^{14} M_{⊙} 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.
Basic properties of the subsamples for the 2D analysis (upper part) and the 3D analysis (lower part).
Fig. 3.
Distribution of M_{500} 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 R_{500}. Indeed, we prefer to consider the quantity mean value computed by averaging it over the radial bins from 0.9 to 1.1 R_{500}, instead of using the interpolation. This procedure limits the dependence of our results on the precise radial binning adopted.
4.1. Twodimensional gas inhomogeneity
Starting from the Xray flux maps, we extract different indicators of the gas inhomogeneities:
ε: the ellipticity of the bestfitting ellipse to the external isoflux contour (Sect. 2.2);
σ_{A}: the azimuthal scatter of the Xray surface brightness profiles (various ways of calculating this quantity are described below);
MM: the ratio between the mean and median of the Xray 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 R_{500}. Indeed, in these few cases the projected annulus around R_{500} computed from the Xray 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 Xray 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 R_{500} to 1.2 R_{500} and linearly equally spaced with a distance equal to 5% of R_{500} (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):
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):
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, MM^{CE} and MM^{MF}, 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 R_{500} in Sect. 7. For the lattermentioned 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 R_{500}). We consider it in relation to the HE mass bias at R_{500} 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 Xray observations, whereas the clumpiness factor can be precisely measured in simulations by adopting its definition:
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):
where m_{i} and ρ_{i} are the mass and density of the ith gas particle. The sum is extended over all the gas particles used for the observationaloriented quantities, in other words not starforming and with temperature above 10^{6} K in order to consider only Xrayemitting 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 largescale 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 Σ_{i}m_{i}/ρ_{i} = Σ_{i}(V_{i}) = 0.01 × V_{shell}, where V_{i} and V_{shell} are the volumes associated with the ith particle and the shell, respectively. For the three representative cases already shown in Fig. 4, we plot both the clumpiness and the residual clumpiness^{10} 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 R_{500}. Among the irregular (IR_{3D}) 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 adaptivemeshrefinement (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 hydrostaticequilibrium 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. Hydrostaticequilibrium 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 bestfitting analytic formulae that could appropriately reproduce each profile. The formulae adopted are taken from observational work. In this way, we follow the Xray 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 R_{500} and 1.2 R_{500} to search for the best constraints on the profiles around R_{500}.
4.3.1. Gas density
Rasia et al. (2006) and Nagai et al. (2007b) proved the reliability of the Xray reconstruction of the gas density profiles. Both works analyzed simulated clusters and produced mock Xray images including Chandra ACISS and ACISI responses. The quantities obtained from the Xray analysis, such as the surface brightness profiles and the projected and deprojected gas density profiles, were found to agree with the input simulated data set (see also Meneghetti et al. 2010 who tested different Xray procedures and Avestruz et al. 2014 who extended the Xray 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, ρ = Σm_{i}/V_{shell}. Each gas density profile is fitted by the (simplified) parametric formula by Vikhlinin et al. (2006):
where ρ_{0}, r_{0}, r_{s}, β, 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 R_{500} and the radial range investigated excludes the central 40% of R_{500}.
In Sect. 7, we refer also to the asymptotic external slope (for r ≫ r_{0} and r ≫ r_{s}) 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, r_{s}, should be larger than the core radius of the betamodel, r_{0}. 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 r_{0} or r_{m}). 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 Xray 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 Xray emission to weight the temperature leads to biased results (Gardini et al. 2004) and that another definition should be used to reproduce Chandra or XMMNewton measurements: the spectroscopiclike temperature (Mazzotta et al. 2004). Nonetheless, the spectroscopiclike 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 Xray mass. Instead, Xray observers deproject the projected temperature profile using the gas density obtained from the imaging as input for the spectroscopiclike weighting. The final Xray deprojected temperature profile, already deconvolved from the instrumental response, can then be considered as the “true” (unweighted) 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 coexistence of multitemperature 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 massweighted temperature: T = Σ(m_{i} × T_{i})/Σm_{i}, where the i  th gas particle has a temperature greater than the lowest energy band of current Xray telescopes, 0.3 keV (≈3.5 × 10^{6} K). Indeed, the massweighted temperature is the one to be considered for the derivation of the cluster mass under the assumption of hydrostatic equilibrium^{11}. 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
where T_{0}, r_{0}, α_{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 massweighted temperature profiles.
The pressure profile is described by a generalized NavarroFrenkWhite model (Nagai et al. 2007a; Arnaud et al. 2010):
where P_{0}, r_{p}, α_{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 ≫ r_{p}) 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 bestfitting parameters, we also save the normalized rootmeansquare (NRMS) value as a measure of the goodness of the fit:
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 bestfitting 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:
This expression has been used to exploit the advantages of both SZ and Xray 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);
This latter equation is typically used in Xray 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);
this latest hybrid version helps to separate the influence on the massbias calculation of the linear multiplicative factor and the term with the derivatives sum.
In all equations, k_{b} is the Boltzmann constant and A = 1/(Gμm_{p}) = 3.7 × 10^{13} M_{⊙} keV^{−1}, where G, m_{p}, 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 hydrostaticequilibrium 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
for each of the three versions of M_{HE}. The bias, b_{HE}, is zero when the HE mass coincides with the true mass, while it is negative or positive for overestimated or underestimated values of M_{HE}.
For the representative clusters shown in Fig. 1, we show the corresponding profiles of (1 − b_{HE, 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 Xray 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 R_{500}. 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 largescale 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 R_{500} 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 R_{500} 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 Xray 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 bestfitting ellipse.
The increase in the scatter in the IR class results from a larger offset between the median and the mean of the surfacebrightness 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 R_{500}, 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 onequarter of the IR systems. On the other hand, for radii between 0.8 and 1.2 R_{500}, 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 R_{500} 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 R_{500}.
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, onequarter 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 bestfitting 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 VR–R 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 R_{500}. The maximum value of σ_{A} at R_{500} 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 R_{500} 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 Xray contours.
6. Results: 3D clumpiness factor
6.1. Clumpiness for the R_{3D} and the IR_{3D} 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 smallscale 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 R_{3D} and IR_{3D}.
Fig. 9.
Clumpiness profiles for the R_{3D} class in red and for the IR_{3D} 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 IR_{3D} 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 R_{500} (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, threequarters of their points are in the bottomleft 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 outlierresistant twovariable linear regression routine in IDL (ROBUST_LINEFIT performed with the bisector method). The bestfitting 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:
with σ_{0} and r_{0} approximately equal to 16 and 6 × R_{200} 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 R_{500}), the second with 0.6 < R/R_{500} < 0.8, and the third with R > 0.8 R_{500}. 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 R_{500}, where the clumpiness factor is still reduced.
Before investigating the region around R_{500} 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 bottomright 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 topleft 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 × R_{500}. 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 R_{500}. 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 R_{500}. 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 R_{3D} objects are shown in brown, and the IR_{3D} 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 identical^{12}. The scatter on the relation is still significant, but drastically reduces when we average the three scatters for each cluster, ⟨σ_{A}⟩. The bottomright 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 topleft 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(𝒞(R_{500}),⟨σ_{A}(R_{500})⟩) = 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 R_{500} 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 largescale 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 Xray mass, M_{HE, X}, (ii) the SZ/Xray mass, M_{HE, SZ}, and (iii) what we called the hybrid estimator, M_{HE, 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 thermodynamical 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 − b_{HE}) 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 Xray 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 R_{500} 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.
Summary of HE mass bias result at R_{500}.
Spearman rank correlation coefficient and the number of standard deviations from the nullhypothesis 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 R_{500}.
7.1. The mass bias
We show in Fig. 12 the median HE mass bias profiles, (1 − b_{HE}), 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, M_{HE} is given by Eq. (10) for b_{HE, X}, Eq. (9) for b_{HE, SZ} and Eq. (11) for b_{HE, 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 R_{500} (which is approximately R_{2500}) and at R_{500} 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 IR_{3D} subsample exceeds the respective percentiles of the R_{3D} distribution. These findings are common to most studies based on the direct analysis of simulated samples. At R_{500} the shaded area is above the value (1 − b) = 0.80 in all panels. In the entire sample (R_{3D} plus IR_{3D}), we find that only four clusters (less than 2.5%) have (1 − b_{HE, X}) < 0.70, in conflict with the mass bias required to solve the discrepancy on the cosmological parameters derived from cluster number counts and cosmicmicrowavebackground power spectrum (see Salvati et al. 2019, and their discussion).
Fig. 12.
Median profile of the HE mass bias: (1 − b_{HE, X}), (1 − b_{HE, SZ}), and (1 − b_{HE, 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 b_{HE, X} and b_{HE, 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 R_{500}. 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 b_{HE, X}(R_{500}) has the lowest standard deviation, σ(b_{HE, X}) = 0.09, the other two biases have σ(b_{HE, SZ}) = 0.11 and σ(b_{HE, T}) = 0.12 (see the first row of Table 2 for the first two bias expressions). It also appears that the distribution of (1 − b_{HE, X}) is quite symmetric with respect to the median behavior in the entire radial range. Specifically, at R_{500}, b_{HE, X} has a low value of skewness (0.20). Alternatively, the other two biases, b_{HE, SZ} and b_{HE, 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 R_{500} 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 b_{HE, SZ} and b_{HE, T} highlight the influence of the multiplicative factor: temperature, rather than the ratio between pressure and gas density. The SZ HE mass bias, (1 − b_{HE, SZ}), shows a general shift towards the value of 1, which is more evident for the IR_{3D} 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 − b_{HE, SZ}) is present at all radii, including R_{500}. There, the number of clusters with (1 − b_{HE, SZ}) > 1 is 50% higher than those with (1 − b_{HE, 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 − b_{HE, T}) < 0.75 while only 6% show the same amount of bias using the SZ formulation.
As a summary, we can conclude that the Xray and the SZ mass biases are substantially providing the same answer but the derivative of the pressure in b_{HE, 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 M_{HE, SZ} less appealing because modeling its final distribution is more difficult.
7.2. Bias and 3D clumpiness factor
In Fig. 13 we relate the Xraymass and the SZmass bias to the clumpiness factor. All quantities are computed in 3D and around R_{500}. 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 “wellfitted” clusters.
Fig. 13.
Mass bias as a function of clumpiness. All quantities are considered at R_{500}. 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 R_{3D} objects are show in brown and the IR_{3D} clusters in olive green. 
The most interesting feature is that most of the clusters with positive bias, (1 − b_{HE, X}) ≥ 1 or (1 − b_{HE, SZ}) ≥ 1, are associated with poorly fitted systems and that the only wellfitted system with high positive bias, (1 − b_{HE, X}) ≥ 1.10 and (1 − b_{HE, 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 wellfitted systems, even though the former have a much higher dispersion, as we can see by comparing the results of the wellfitted 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 systems^{13}. 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 R_{500}). The Spearman correlation coefficient is indeed corr(1 − b_{HE}, 𝒞)≲0.20 for the entire sample and corr(1 − b_{HE}, 𝒞)≲0.30 for the wellfitted subsample. Independent simulations have indeed shown that even after removing the 50% densest cells at all cluster radii, the ratio of nonthermal 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 − b_{HE}) by looking for a linear relation between (1 − b_{HE}) and 𝒞. By subtracting the bestfitting 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 − b_{HE, SZ}) to the value of β_{P} (the asymptotic slope of the pressure profile as in Eq. (7)) and (1 − b_{HE, 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 wellfitted 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 threequarters of all wellfitted 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 − b_{HE, X}) and (1 − b_{HE, 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 wellfitted 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 75^{th} 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 − b_{HE, SZ}, β_{P}) ≈ corr(b_{HE, X}, 𝒟)≈0.35 − 0.45. In light of these results, we try to correct the mass bias by using the bestfitting plane of bias, slope, and clumpiness values. The fitting procedure to determine the parameters of the plane was restricted only to the wellfitted clusters. However, we find that correcting all data points (both well and poorly fitted) following this expression for (1 − b_{HE, SZ}):
and equivalently for (1 − b_{HE, X}):
decreases the standard deviations by about 10% (5% for b_{HE, 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 observationoriented approach by relating the Xray and SZ mass bias to all quantities defined from the Xray 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 wellfitted clusters and all 2D proxies (σ_{A}, σ_{A, R}, ε, and MM) are found with respect to (1 − b_{HE, 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 − b_{HE, SZ}), for which the only notable correlation value is −0.24 in relation to σ_{A, R} and for the wellfitted 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 − b_{HE, X}) or (1 − b_{HE, 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 massbias statistics of the subsample with the lowest azimuthal scatter, σ_{A} < 0.4, to that of the subsample with the highest scatter^{12}, σ_{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 Xray mass bias:
At the same time, we correct the SZ mass bias by invoking the asymptotic slope of the pressure profile, β_{P}:
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 Xray 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 wellfitted 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 − b_{HE, X}), on the left, and (1 − b_{HE, 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 wellfitted 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 Xray 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 GADGETX code, which includes an improved formulation of SPH, with respect to GADGET2 (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 supermassive black holes (see Rasia et al. 2015). Xray 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 isoflux contour at around 0.8 R_{500}. The distance of both centers from the theoretical center, the minimum of the potential well, and the ellipticity of the bestfit ellipse are used to provide a first separation of the clusters into different morphological classes: veryregular, VR, regular, R, and intermediateirregular, IR, objects. From the surfacebrightness 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] R_{500}), 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 substructurefree 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 Xray and SZ observational works, M_{HE, X} (Eq. (10)) and M_{HE, SZ} (Eq. (9)) respectively, and a third formulation, M_{HE, 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 510% around R_{2500} to 10–15% around R_{500}. Regular clusters (VR or R), less clumped clusters (𝒞 < 1.1), objects with low azimuthalscatter (σ_{A} < 0.4), and systems with wellbehaved 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 − b_{HE, 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 stateoftheart 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 Xray images cannot provide a compelling correction to the mass bias, it is still possible to use information from Xray 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 Xray observations (such as those from XMMNewton, Chandra, or eRosita) 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 XCOP collaboration (Ghirardini et al. 2019) using XMMNewton and Planck data, and in Ruppin et al. (2017) with the combined analysis of NIKA, Planck, and XMMNewton 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 Xray 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.
The MultiDark simulations are publicly available at the https://www.cosmosim.org database.
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.
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.
We note that, in a similar way, objects with a very perturbed Xray 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 adhoc analysis procedures (e.g. Ghirardini et al. 2018).
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 R_{3D} class, while the latter is part of the IR_{3D} class.
To ease the comparison with other numerical works that use the spectroscopiclike temperature in the HE mass derivation rather than the massweighted temperature, we computed the ratio of these two temperatures at R_{500}. 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).
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 SklodowskawCurie grant agreement number 734374, the LACEGAL project. We further acknowledge financial contribution from the contracts/grants: AYA201563810P (W.C. and G.Y.); ASI 2015046R.0 (S.E.); ASIINAF n.201714H.0 (E.R., S.B., S.E.); PGC2018094975BC21 (G.Y.); PRINMIUR 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/21 (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
 Allen, S. W., Rapetti, D. A., Schmidt, R. W., et al. 2008, MNRAS, 383, 879 [NASA ADS] [CrossRef] [Google Scholar]
 Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409 [NASA ADS] [CrossRef] [Google Scholar]
 Ameglio, S., Borgani, S., Pierpaoli, E., et al. 2009, MNRAS, 394, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Angelinelli, M., Vazza, F., Giocoli, C., et al. 2019, MNRAS, submitted [arXiv:1905.04896] [Google Scholar]
 Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arthur, J., Pearce, F. R., Gray, M. E., et al. 2019, MNRAS, 484, 3968 [NASA ADS] [CrossRef] [Google Scholar]
 Avestruz, C., Lau, E. T., Nagai, D., & Vikhlinin, A. 2014, ApJ, 791, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Baldi, A. S., Bourdin, H., Mazzotta, P., et al. 2019, A&A, 630, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barnes, D. J., Kay, S. T., Bahé, Y. M., et al. 2017a, MNRAS, 471, 1088 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, D. J., Kay, S. T., Henson, M. A., et al. 2017b, MNRAS, 465, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012, ApJ, 758, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2015, ApJ, 806, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110 [Google Scholar]
 Biffi, V., & Valdarnini, R. 2015, MNRAS, 446, 2802 [NASA ADS] [CrossRef] [Google Scholar]
 Biffi, V., Dolag, K., & Böhringer, H. 2013, MNRAS, 428, 1395 [NASA ADS] [CrossRef] [Google Scholar]
 Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Biffi, V., Planelles, S., Borgani, S., et al. 2017, MNRAS, 468, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Biffi, V., Planelles, S., Borgani, S., et al. 2018, MNRAS, 476, 2689 [NASA ADS] [CrossRef] [Google Scholar]
 Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Bourdin, H., Mazzotta, P., & Rasia, E. 2015, ApJ, 815, 92 [NASA ADS] [CrossRef] [Google Scholar]
 Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, H., Avestruz, C., Kravtsov, A. V., Lau, E. T., & Nagai, D. 2019, MNRAS, 490, 2380 [Google Scholar]
 Churazov, E., Sazonov, S., Sunyaev, R., et al. 2005, MNRAS, 363, L91 [NASA ADS] [Google Scholar]
 Cialone, G., De Petris, M., Sembolini, F., et al. 2018, MNRAS, 477, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Clerc, N., Cucchetti, E., Pointecouteau, E., & Peille, P. 2019, A&A, 629, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Costanzi, M., Rozo, E., Simet, M., et al. 2019, MNRAS, 488, 4779 [NASA ADS] [CrossRef] [Google Scholar]
 Cucchetti, E., Clerc, N., Pointecouteau, E., Peille, P., & Pajot, F. 2019, A&A, 629, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cui, W., Power, C., Biffi, V., et al. 2016, MNRAS, 456, 2566 [NASA ADS] [CrossRef] [Google Scholar]
 Cui, W., Knebe, A., Yepes, G., et al. 2018, MNRAS, 480, 2898 [NASA ADS] [CrossRef] [Google Scholar]
 Davé, R., AnglésAlcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827 [NASA ADS] [CrossRef] [Google Scholar]
 Dolag, K., Hansen, F. K., Roncarelli, M., & Moscardini, L. 2005, MNRAS, 363, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Eckert, D., Roncarelli, M., Ettori, S., et al. 2015, MNRAS, 447, 2198 [NASA ADS] [CrossRef] [Google Scholar]
 Eckert, D., Gaspari, M., Vazza, F., et al. 2017, ApJ, 843, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ettori, S., Morandi, A., Tozzi, P., et al. 2009, A&A, 501, 61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ettori, S., Donnarumma, A., Pointecouteau, E., et al. 2013, Space Sci. Rev., 177, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Ettori, S., Ghirardini, V., Eckert, D., et al. 2019, A&A, 621, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Evrard, A. E. 1990, ApJ, 363, 349 [NASA ADS] [CrossRef] [Google Scholar]
 Fang, T., Humphrey, P., & Buote, D. 2009, ApJ, 691, 1648 [NASA ADS] [CrossRef] [Google Scholar]
 Frank, K. A., Peterson, J. R., Andersson, K., Fabian, A. C., & Sanders, J. S. 2013, ApJ, 764, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Gardini, A., Rasia, E., Mazzotta, P., et al. 2004, MNRAS, 351, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Gaspari, M., McDonald, M., Hamer, S. L., et al. 2018, ApJ, 854, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Henson, M. A., Barnes, D. J., Kay, S. T., McCarthy, I. G., & Schaye, J. 2017, MNRAS, 465, 3361 [NASA ADS] [CrossRef] [Google Scholar]
 Hitomi Collaboration (Aharonian, F., et al.) 2016, Nature, 535, 117 [Google Scholar]
 Huang, S., Katz, N., Davé, R., et al. 2019, MNRAS, 484, 2021 [NASA ADS] [CrossRef] [Google Scholar]
 Jeltema, T. E., Hallman, E. J., Burns, J. O., & Motl, P. M. 2008, ApJ, 681, 167 [Google Scholar]
 Khatri, R., & Gaspari, M. 2016, MNRAS, 463, 655 [NASA ADS] [CrossRef] [Google Scholar]
 Khedekar, S., Churazov, E., Kravtsov, A., et al. 2013, MNRAS, 431, 954 [NASA ADS] [CrossRef] [Google Scholar]
 Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340 [NASA ADS] [CrossRef] [Google Scholar]
 Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353 [NASA ADS] [CrossRef] [Google Scholar]
 Laganá, T. F., Durret, F., & Lopes, P. A. A. 2019, MNRAS, 484, 2807 [NASA ADS] [CrossRef] [Google Scholar]
 Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129 [NASA ADS] [CrossRef] [Google Scholar]
 Lopes, P. A. A., Trevisan, M., Laganá, T. F., et al. 2018, MNRAS, 478, 5473 [NASA ADS] [CrossRef] [Google Scholar]
 Lovisari, L., Forman, W. R., Jones, C., et al. 2017, ApJ, 846, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2014, MNRAS, 440, 2077 [NASA ADS] [CrossRef] [Google Scholar]
 Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2015, MNRAS, 449, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Mathiesen, B., Evrard, A. E., & Mohr, J. J. 1999, ApJ, 520, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Mazzotta, P., Rasia, E., Moscardini, L., & Tormen, G. 2004, MNRAS, 354, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Meneghetti, M., Rasia, E., Merten, J., et al. 2010, A&A, 514, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morandi, A., & Cui, W. 2014, MNRAS, 437, 1909 [NASA ADS] [CrossRef] [Google Scholar]
 Morandi, A., Nagai, D., & Cui, W. 2013, MNRAS, 436, 1123 [NASA ADS] [CrossRef] [Google Scholar]
 Morandi, A., Sun, M., Mulchaey, J., Nagai, D., & Bonamente, M. 2017, MNRAS, 469, 2423 [NASA ADS] [CrossRef] [Google Scholar]
 Mostoghiu, R., Knebe, A., Cui, W., et al. 2019, MNRAS, 483, 3390 [CrossRef] [Google Scholar]
 Nagai, D., & Lau, E. T. 2011, ApJ, 731, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007a, ApJ, 668, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Nagai, D., Vikhlinin, A., & Kravtsov, A. V. 2007b, ApJ, 655, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, K., Rudd, D. H., Shaw, L., & Nagai, D. 2012, ApJ, 751, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Pearce, F. A., Kay, S. T., Barnes, D. J., Bower, R. G., & Schaller, M. 2019, MNRAS, 2606 [Google Scholar]
 Piffaretti, R., & Valdarnini, R. 2008, A&A, 491, 71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XX. 2014, A&A, 571, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planelles, S., Fabjan, D., Borgani, S., et al. 2017, MNRAS, 467, 3827 [NASA ADS] [CrossRef] [Google Scholar]
 Power, C., Elahi, P. J., Welker, C., et al. 2019, MNRAS, 2758 [Google Scholar]
 Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Rasia, E., Tormen, G., & Moscardini, L. 2004, MNRAS, 351, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Rasia, E., Ettori, S., Moscardini, L., et al. 2006, MNRAS, 369, 2013 [NASA ADS] [CrossRef] [Google Scholar]
 Rasia, E., Meneghetti, M., Martino, R., et al. 2012, New J. Phys., 14, 055018 [NASA ADS] [CrossRef] [Google Scholar]
 Rasia, E., Meneghetti, M., & Ettori, S. 2013, Astron. Rev., 8, 40 [Google Scholar]
 Rasia, E., Lau, E. T., Borgani, S., et al. 2014, ApJ, 791, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Rasia, E., Borgani, S., Murante, G., et al. 2015, ApJ, 813, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Richardson, M. L. A., Scannapieco, E., Devriendt, J., et al. 2016, ApJ, 825, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Roncarelli, M., Ettori, S., Borgani, S., et al. 2013, MNRAS, 432, 3030 [NASA ADS] [CrossRef] [Google Scholar]
 Roncarelli, M., Gaspari, M., Ettori, S., et al. 2018, A&A, 618, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rossetti, M., Gastaldello, F., Ferioli, G., et al. 2016, MNRAS, 457, 4515 [Google Scholar]
 Ruppin, F., Adam, R., Comis, B., et al. 2017, A&A, 597, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ruppin, F., Mayet, F., Pratt, G. W., et al. 2018, A&A, 615, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salvati, L., Douspis, M., Ritz, A., Aghanim, N., & Babul, A. 2019, A&A, 626, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sembolini, F., Yepes, G., De Petris, M., et al. 2013, MNRAS, 429, 323 [NASA ADS] [CrossRef] [Google Scholar]
 Sembolini, F., Elahi, P. J., Pearce, F. R., et al. 2016, MNRAS, 459, 2973 [NASA ADS] [CrossRef] [Google Scholar]
 Shi, X., Komatsu, E., Nelson, K., & Nagai, D. 2015, MNRAS, 448, 1020 [NASA ADS] [CrossRef] [Google Scholar]
 Shitanishi, J. A., Pierpaoli, E., Sayers, J., et al. 2018, MNRAS, 481, 749 [NASA ADS] [CrossRef] [Google Scholar]
 Simionescu, A., ZuHone, J., Zhuravleva, I., et al. 2019, Space Sci. Rev., 215, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 [Google Scholar]
 Steinborn, L. K., Dolag, K., Hirschmann, M., Prieto, M. A., & Remus, R.S. 2015, MNRAS, 448, 1504 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suto, D., Kawahara, H., Kitayama, T., et al. 2013, ApJ, 767, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Tornatore, L., Borgani, S., Dolag, K., & Matteucci, F. 2007, MNRAS, 382, 1050 [NASA ADS] [CrossRef] [Google Scholar]
 Truong, N., Rasia, E., Mazzotta, P., et al. 2018, MNRAS, 474, 4089 [NASA ADS] [CrossRef] [Google Scholar]
 Truong, N., Rasia, E., Biffi, V., et al. 2019, MNRAS, 484, 2896 [NASA ADS] [CrossRef] [Google Scholar]
 Urban, O., Simionescu, A., Werner, N., et al. 2014, MNRAS, 437, 3939 [NASA ADS] [CrossRef] [Google Scholar]
 Vafaei Sadr, A., Movahed, S. M. S., Farhang, M., Ringeval, C., & Bouchet, F. R. 2018, MNRAS, 475, 1010 [NASA ADS] [CrossRef] [Google Scholar]
 van Weeren, R. J., de Gasperin, F., Akamatsu, H., et al. 2019, Space Sci. Rev., 215, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Vazza, F., Brunetti, G., Kritsuk, A., et al. 2009, A&A, 504, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vazza, F., Roncarelli, M., Ettori, S., & Dolag, K. 2011, MNRAS, 413, 2305 [NASA ADS] [CrossRef] [Google Scholar]
 Vazza, F., Eckert, D., Simionescu, A., Brüggen, M., & Ettori, S. 2013, MNRAS, 429, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Vazza, F., Wittor, D., Brüggen, M., & Gheller, C. 2016, Galaxies, 4, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Vazza, F., Angelinelli, M., Jones, T. W., et al. 2018, MNRAS, 481, L120 [NASA ADS] [CrossRef] [Google Scholar]
 Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691 [NASA ADS] [CrossRef] [Google Scholar]
 Vikhlinin, A., Kravtsov, A. V., Burenin, R. A., et al. 2009, ApJ, 692, 1060 [NASA ADS] [CrossRef] [Google Scholar]
 Voit, G. M. 2005, Rev. Mod. Phys., 77, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Walker, S. A., Fabian, A. C., Sanders, J. S., & George, M. R. 2012, MNRAS, 424, 1826 [NASA ADS] [CrossRef] [Google Scholar]
 Walker, S., Simionescu, A., Nagai, D., et al. 2019, Space Sci. Rev., 215, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y., Pearce, F., Knebe, A., et al. 2018, ApJ, 868, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99 [Google Scholar]
 Zhuravleva, I., Churazov, E., Kravtsov, A., et al. 2013, MNRAS, 428, 3274 [NASA ADS] [CrossRef] [Google Scholar]
 Zinger, E., Dekel, A., Birnboim, Y., et al. 2018, MNRAS, 476, 56 [NASA ADS] [CrossRef] [Google Scholar]
 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 M_{500} 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 M_{500} > 6 × 10^{14} M_{⊙}. Half of the remaining clusters has a mass below 10^{14} M_{⊙} and the others have a mass between 1 and 2 × 10^{14} M_{⊙}. 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 Xray and SZ mass bias, (1 − b_{HE}), 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 wellfitted 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 × 10^{14} M_{⊙}), 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 − b_{HE}) and M_{500}, the correlation being below 0.10 for both biases.
Fig. A.1.
Hydrostatic mass bias vs. the cluster mass at R_{500}. Top panel: Xray mass bias, (1 − b_{HE, X}). Bottom panel: SZ mass bias, (1 − b_{HE, SZ}). Similarly to Fig. 14 of the main text, wellfitted 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 veryirregular (VI) and extremelyirregular (EI) objects. The former class is excluded because the center of the Xray 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 Xray emission is strongly asymmetric. For these reasons, as mentioned in the main text, we consider the extension of the analysis to the extremelyirregular 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 R_{500}.
In Table B.1, we report median value, standard deviation, and skewness of the hydrostatic mass bias distribution of the veryirregular clusters. The average and standard deviation are 50% higher in this class with respect to the first three lessdisturbed 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.
Summary of results on the HE mass bias at R_{500}.
Including the veryirregular 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 Xray 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%.
Spearman rank correlation coefficient and the number of standard deviations from the nullhypothesis expected value (in parenthesis) similarly to Table 3 but now including also the VI class.
Appendix C: The HE mass bias at R_{2500}
In the VR, R, and IR classes, the R_{2500} radius is on average equal to 0.45 R_{500} 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 R_{2500}, and we further link the mass bias at the same radius with the ellipticity value of the Xray isoflux contours drawn at 0.8 R_{500}, 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 R_{2500} radius is smaller than our innermost radial bin, 0.4 R_{500}. This sample reduction led to a total of 467 maps.
On average we find that (1 − b_{HE, X, 2500})∼(1 − b_{HE, 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 Xray and SZ mass are thus slightly closer to the true mass at R_{2500} than at R_{500} as suggested also from Fig. 12 where we notice an increase in the departure from the true mass at increasing radius. At R_{2500} both Xray 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(b_{HE, SZ}) = 0.96.
The correlation coefficients between mass bias and measurements of the gas inhomogeneities at R_{2500} 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, b_{HE, X}(R_{500}) exhibits the strongest correlation with 𝒟, while at R_{2500} the two quantity are completely uncorrelated.
Spearman rank correlation coefficient and the number of standard deviations from the nullhypothesis expected value (in parenthesis) similarly to Table 3.
All Tables
Basic properties of the subsamples for the 2D analysis (upper part) and the 3D analysis (lower part).
Spearman rank correlation coefficient and the number of standard deviations from the nullhypothesis 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 R_{500}.
Spearman rank correlation coefficient and the number of standard deviations from the nullhypothesis expected value (in parenthesis) similarly to Table 3 but now including also the VI class.
Spearman rank correlation coefficient and the number of standard deviations from the nullhypothesis expected value (in parenthesis) similarly to Table 3.
All Figures
Fig. 1.
Xray 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, M_{500} ≈ 1 × 10^{15} M_{⊙}, corresponding to R_{500} of about 1.5 Mpc (shown as a white circle in the images). The black lines are the logspaced isoflux contours smoothed over 4 pixels. 

In the text 
Fig. 2.
Distribution of D_{CE} (in units of R_{500}) 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 D_{CE} and ε of each class. 

In the text 
Fig. 3.
Distribution of M_{500} 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 hydrostaticequilibrium 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 R_{3D} class in red and for the IR_{3D} 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 R_{500}. 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 R_{3D} objects are shown in brown, and the IR_{3D} clusters in olive green. 

In the text 
Fig. 12.
Median profile of the HE mass bias: (1 − b_{HE, X}), (1 − b_{HE, SZ}), and (1 − b_{HE, 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 R_{500}. 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 R_{3D} objects are show in brown and the IR_{3D} clusters in olive green. 

In the text 
Fig. 14.
Mass biases (1 − b_{HE, X}) and (1 − b_{HE, 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 wellfitted 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 75^{th} 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 wellfitted 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 − b_{HE, X}), on the left, and (1 − b_{HE, 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 wellfitted 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 R_{500}. Top panel: Xray mass bias, (1 − b_{HE, X}). Bottom panel: SZ mass bias, (1 − b_{HE, SZ}). Similarly to Fig. 14 of the main text, wellfitted 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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.