Spatial distribution of dark matter in and around galaxy clusters traced by galaxies, gas and intracluster stars in a simulated universe

To understand how well galaxies, gas and intracluster stars trace dark matter in and around galaxy clusters, we use the IllustrisTNG cosmological hydrodynamical simulation and compare the spatial distribution of dark matter with those of baryonic components in clusters. To quantify the global morphology of the density distribution of each component in clusters, we fit an ellipse to the density contour of each component and derive shape parameters at different radii. We find that ellipticity of dark matter is better correlated with that of galaxy mass-weighted number density, rather than with that of galaxy number density or galaxy velocity dispersion. We thus use the galaxy mass-weighted number density map as a representative of the galaxy maps. Among three different density maps from galaxies, gas, and intracluster stars, the ellipticity of dark matter is best reproduced by that of the galaxy map over the entire radii. The 'virialized' galaxy clusters show a better correlation of spatial distribution between dark matter and other components than the 'unvirialized' clusters, suggesting that it requires some time for each component to follow the spatial distribution of dark matter after merging events. Our results demonstrate that galaxies are still good tracers of dark matter distribution even in the non-linear regime corresponding to the scales in and around galaxy clusters, being consistent with the case where galaxies trace well the matter distribution in cosmologically large scales.


INTRODUCTION
Galaxy clusters are important probes in studying cosmology and structure formation (Allen et al. 2011). In particular, mapping the mass distribution in and around them provides valuable information that can be used for tions for testing structure formation models. The large masses and three baryonic ingredients 1 .
The total mass of clusters as a function of redshift can offer useful constraints on cosmological parameters (Vikhlinin et al. 2009;Bocquet et al. 2019). The radial profile of cluster masses derived from observations can not only be directly compared with the (universal) profile of dark matter halos in simulations (e.g., Bullock et al. 2001), but also be used for comparisons of mass measurements from different tracers Rines et al. 2013).
The two-dimensional maps of cluster mass distribution also provide strong constraints on the mass growth of clusters and the formation of galaxies around them: e.g., existence of dark structures without galaxies (Clowe et al. 2012;Jee et al. 2014), and subhalo mass function (Okabe et al. 2014). However, such maps can highly depend on the systematics of the tracers used for the map construction, so it is very important to cross check the mass measurements from different methods (Rines et al. 2016;Sohn et al. 2018). In this regard, Geller et al. (2014) showed a seminal example of such a comparison between weak lensing and redshift surveys for Abell 383 (see their Figure 8). Their results indicate that the weak-lensing cluster masses could be affected by foreground/background structures even within the virial radii of clusters (see also Hoekstra 2001;Dodelson 2004). This comparison demonstrates the importance of dense redshift surveys to understand the mass distribution in weak-lensing maps. Hwang et al. (2014) extended such comparisons to nine relaxed clusters for more quantitative analysis using cross-correlation technique between the weak-lensing convergence map and the galaxy number density map. This line of study was further extended to a merging cluster that has both weak lensing and redshift survey data exist (Liu et al. 2018); many substructures with multiple merging processes in merging clusters make the comparison of mass distributions from different components more interesting as in the Bullet cluster (Clowe et al. 2004). In these analyses, galaxies were used as test particles with the assumption that they trace dark matter well. This issue has been 1 The total mass of galaxy clusters is 10 14 M ; the contributions of dark matter, galaxies, intracluster objects and hot gas are > 80%, 1 − 2%, < 1% and 5 − 15%, respectively (Böhringer 2002;Rudick et al. 2011). (i.e. galaxies, intracluter objects, and hot gas with temperature Tgas 10 7 K) of galaxy clusters provide various tools to trace the underlying dark matter distribution: e.g., galaxy redshift surveys for dynamical analysis Hwang et al. 2014;Song et al. 2017) optical/near-infrared imaging surveys for weak-lensing analysis (Okabe et al. 2010;Umetsu 2020) and X-ray and microwave surveys for the analysis of hot gas (Ettori et al. 2013;Planck Collaboration et al. 2013). tested since Kaiser (1984), and the galaxies as biased tracers of dark matter is reviewed thoroughly in Desjacques et al. (2018).
Here we use the cosmological hydrodynamic simulations where we know well the distributions of dark matter and other components including galaxies, gas and intracluster stars in clusters. Then the comparison of spatial distribution of dark matter with those of other components can provide an important test for the reliability of the use of galaxies to infer the dark matter distribution in galaxy clusters (i.e. non-linear regime). Among many ways to analyze the spatial distribution of each component in clusters, we focus on the comparison of their global morphology; this includes ellipticity, position angle, and centroid offset of the projected density distribution of each component derived from the ellipse fitting method (see Section 2.3 for details). We thus adopt the method fitting the contours of mass or number density of each component with ellipse to derive shape parameters at different radii. We note that the ellipse fitting may miss some possible small-scale differences among different components, but we think that our approach can provide quantitative and straightforward comparisons among different tracers; for example, see Montes & Trujillo (2019) for non-parametric methods.
On the other hand, thanks to recent wide-field, deep images available for galaxy clusters, the intracluster objects other than galaxies including globular clusters (e.g. West et al. 1995;Lee et al. 2010;Ko et al. 2017), planetary nebulae (e.g. Feldmeier et al. 2004;Longobardi et al. 2018), and stars (e.g. Mihos et al. 2005;Ko & Jee 2018) have become valuable probes of dark matter in clusters. In particular, whether the intra-cluster light (ICL) could be an excellent tracers of underlying dark matter distribution or not has been seriously examined (Montes & Trujillo 2019;Alonso Asensio et al. 2020;Sampaio-Santos et al. 2020) in both observations and simulations. A study focusing on the detailed connection between ICL and dark matter using various methods will be a topic of separate papers in our group (Yoo et al., accepted). Here we focus on the large-scale spatial distribution of galaxies along with their kinematics in and around clusters in comparison with other components including gas and intracluster stars.
Different mass accretion and/or merging histories of galaxy clusters can affect not only the current dynamical status of galaxy clusters, but also the spatial distribution of dark matter and other components therein (Mohr et al. 1993;Buote & Tsai 1995;Bershady et al. 2000;Lotz et al. 2004;Parekh et al. 2015;Cui et al. 2017;Kim et al. 2022). To investigate how the different dynamical states are imprinted in the large-scale spatial distribution of each component around galaxy clusters, we divide the galaxy clusters into two sub-samples as 'unvirialized and 'virialized, and compare the results between the two.
We describe the simulation data and the method for measuring the shape parameters we adopt in Section 2, and compare one-and two-dimensional distributions of dark matter with those of galaxies (including two-dimensional number density and velocity dispersion maps), gas and intracluster stars in Section 3. We discuss and summarize the results in Sections 4 and 5, respectively. Throughout, we adopt cosmological parameters as in the TNG simulation: Ω Λ,0 = 0.6911, Ω m,0 = 0.3089, Ω b,0 = 0.0486, σ 8 = 0.8159 and n s = 0.9667 and h = 0.6774 (Planck Collaboration et al. 2016). All quoted errors in measured quantities are 1σ, and spatial quantities and coordinates are expressed in comoving units.

Data
We use the publicly released data from the TNG300 (Nelson et al. 2019;Pillepich et al. 2018a;Springel et al. 2018), which is one of a suite of cosmological hydrodynamic simulations of the IllustrisTNG (Pillepich et al. 2018b): the successor of the Illustris simulation (Vogelsberger et al. 2014). The IllustrisTNG is composed of three different volumes whose one side length is roughly 300, 100 and 50 Mpc: TNG300, TNG100, and TNG50, respectively. Among the TNG series, the TNG300 has the largest simulation box and thus enables us to obtain the most reliable statistics on galaxy clusters with the least cosmic variance even though the mass (spatial) resolution of dark matter particle is limited most as 5.9 × 10 7 M (1.48 kpc); the resolutions of dark matter for TNG100 and TNG50 are 7.5×10 6 M (0.74 kpc) and 4.5 × 10 5 M (0.29 kpc), respectively. Because we are interested in a comparison of the global morphology of the clusters matter distribution with that of each different component, rather than in the detailed structures that are highly affected by a choice of the resolution, we decide to use the TNG300 that gives best statistics for this study. Note that projected mass density maps generated for this study (see Section 2.2) are smoothed by a Gaussian kernel with a value of 90 − 240 kpc, and thus the 1.48 kpc resolution of the TNG300 is much smaller than the resolution that we deal with in this study.
We use the properties of the FoF (friends-of-friends) halos listed in the group catalog, which are identified with the FoF algorithm with a linking length b = 0.2 at z = 0. To choose the galaxy cluster analogs from the catalog, we select massive FoF halos whose enclosed mass in a sphere of R vir or R 200 (the radius whose mean density is 200 times the mean density of the universe) is larger than 10 14 M : the lower-mass limit typically used for galaxy clusters 2 . The number of selected halos considered as galaxy clusters is 426 in total, whose mass and size is ranging 10 14 < M vir / M < 2.9 × 10 15 and 0.9 < R vir /Mpc < 2.4, respectively. We use this sample for following analyses. Here, we consider halos and subhalos in simulations as galaxy clusters and cluster galaxies in observations, respectively.

Projected Mass Density Maps
Using the snapshot at z = 0, we make projected mass density maps of dark matter, gas 3 , and star for each galaxy cluster. To do this, we first extract all particles/cells belonging to each cluster and then project them to three different directions of x-y, y-z, and x-z planes. The projected distribution of the particles/cells is counted in a 2-dimensional array of 900 × 900 pixels 2 using a count-in-cell (CIC) subroutine built in the IDL astronomy user's library (see the 1st columns of Figure 1 and 2 ). Each dimension of the maps covers ±3R vir from the cluster center, where the gravitational potential energy is at the minimum, and thus 15 pixels correspond to 0.1R vir .
Before making the smoothed image of the projected mass density, the local stellar density spikes originated from individual galaxies are removed because stars inside the galaxies tend to reflect the mass distribution of host galaxies rather than that of galaxy clusters; we are interested in global distribution of the diffuse intracluster stars that follows the cluster potential. Note that the individual galaxies are treated as a separated component. Using the subhalo catalog containing galaxy properties measured with the SUBFIND algorithm (Springel et al. 2001), the pixels out to 2R h from each galaxy center are masked, where R h is the radius containing half of each galaxy mass (see the 2nd column of Figure 2). Very small-sized galaxies whose 2R h is smaller than a half of pixel size are omitted from masking and treated as the diffuse intracluster stars. Note that we also apply no mask to the brightest galaxy in the cluster center (i.e. BCG) because the mask for the BCG can remove most of the central region of the cluster.   . Spatial distribution of galaxies in the example cluster. This is color-coded by galaxy mass (left) and by velocity difference from the cluster mean motion (right). All of Figure 1-5 show the same object and are xy projected. Because the brightest cluster galaxy is extremely massive (10 15 M ), the upper limit in the left color is determined by the mass of 2nd BCG.  We first use these projected mass density maps to calculate the one-dimensional density profile for each component (see Section 3.1). We then smooth the projected mass density maps for the comparison of twodimensional distributions (see Section 3.2). To do that, we apply a Gaussian kernel with the standard deviation of 15 pixels to the maps using a GAUSS SMOOTH/IDL routine (see the last columns of Figure 1 and 2). The masked regions for the stellar mass density map are treated as non-existent pixels when pixel values of the unmasked regions are updated by smoothing, but their pixel values can be updated by those of the unmasked regions. Most of masked regions are smoothed out by nearby unmasked regions. However, a large mask (radius> 3σ) produces some defects including image distortion because the routine is set to update a pixel value using the surrounding pixels inside the radius of 3σ. Thus, if the mask size is too large (i.e. radius> 3σ), it even leaves a hole in the image because all pixel values cannot be updated with a real number. We apply the empirical criteria as explained in Section 2.3 to reduce the large mask effect on the ellipse fitting results. Note that the smoothed stellar mass density map after masking represents the mass density of diffuse intracluster stars including the BCG.
To examine in detail how galaxies can be good tracers of dark matter, we generate three different galaxy maps: number density, mass-weighted number density and velocity dispersion (see Figure 3 and 4). For this, we extract position, velocity, and stellar mass of the cluster member galaxies from the subhalo catalog. Here, the galaxy mass reaches down to ∼ 1.2 × 10 9 M , which is the minimum mass (20 dark matter particles) to be identified as a galaxy with the SUBFIND algorithm. The galaxy number density, mass-weighted number density and velocity dispersion (galaxy N, galaxy M, and galaxy σ v , hereafter) are calculated in a 2D array of 300 × 300 pixels 2 , where each dimension covers ±3R vir from the cluster center. Because the number of galaxies is much smaller than those of dark matter, gas, and star parti-cles, we decrease the pixel number by three times compared to those of the other density maps with 900 × 900 pixels 2 to enhance the statistical significance on galaxy counting. We then smoothed the 2D array with σ = 5 pixels, which is comparable to that of density maps for dark matter, gas, and intracluster stars (900 × 900 pixels 2 with σ = 15 pixels). The member galaxies we use as test particles are distributed at discrete positions and have a wide mass range of 10 9 −10 12 M . Thus a simple linear mass-weight scheme in constructing galaxy density maps can leave many clumpy structures. To avoid this issue, we rescale the mass weight so that the member galaxies with minimum and maximum masses have the mass weights of one and 10, respectively (except the BCG). Here, galaxy σ v is a measure of velocity dispersion of galaxies within the pixel having at least three galaxies. To have a fair comparison of spatial distribution among different components and to allocate the same numerical error that arises from the ellipse fitting, we rescale the density maps of dark matter, gas and intra-cluster stars to be 300 × 300 pixels 2 , which is the same as for the galaxy density maps.

Ellipse Fitting
To parameterize the spatial distribution of each component, we perform ellipse fitting on the six different projected maps of dark matter, gas, ICS, galaxy N, galaxy M, and galaxy σ v that are obtained in the previous section. We apply the ellipse fitting to each contour determined at three different radii, R/R vir = 0.5, 1 and 2; these radii were chosen to study the structures of inner, intermediate and outer regions, respectively. We use the radii normalized by R vir to have a fair comparison among clusters. We first estimate the density level of the map at each radius from the average pixel value within the annulus with a width of 0.1R vir . We then extract the contour line at each level with the CON-TOUR/IDL routine and convert the contour line into pixel points for the input of the ellipse fitting. If there are several peaks above the density level, the contour line may consist of several loops with different sizes (e.g. three yellow loops in the top-left panel of Figure 5). In that case, we only adopt the longest loop that best reflects the primary structure. To prevent that the distortion caused by the large masks affects the fitting, the pixel points with the distance from large galaxies (R h > 0.1R vir ) less than 2R h are excluded from the input for the intracluster star density map. We use the MP-FIT/IDL package (Markwardt 2009) to derive the shape parameters of the ellipse: centroid, major axis length, ellipticity and position angle. The best-fitting results of an example cluster are shown in Figure 5.
Samples for the comparison are constrained to satisfy the following conditions. First, massive galaxies with R h > 400kpc should not be located within 2.5R vir from the cluster center to avoid the effect of large masks on the ellipse fitting. There are 14 clusters excluded from the original sample of 426 clusters. They are mostly in a process of cluster (or group) merging and thus have a massive galaxy with R h > 400kpc, which was the BCG (or BGG) of the recently merged cluster (or group). Second, the number fraction of the non-zero pixel used for estimating the level of the map at a given radius should be larger than 30%. Last, the number fraction of the valid pixels along the contour line for a given level should be larger than 80% to obtain reliable ellipse fitting results. Note that the second and third criteria are applied to individual density maps at different radii, while the first criterion to individual clusters. A total number of reliable samples after these selection criteria decreases to 14,973 from 19,170 (426 clusters × 3 projections × 5 components × 3 radii).

1D Radial profiles
We first compare projected radial density profiles among different components. To calculate azimuthally averaged radial density profiles of dark matter, gas and intracluster stars, we use the projected mass density maps before smoothing: the 1st column of Figure 1 (dark matter, gas) and the 2nd panel of Figure 2 (intracluster stars). The galaxy number density profiles with and without mass-weight are calculated from the 1st row of Figure 4. Figure 6 shows the radial density profile of the 426 clusters for each component: dark matter, galaxy N (i.e. number density without mass weight), galaxy M (i.e. number density with mass weight), gas, and intracluster stars, as well as their averaged density profiles. Since the dark matter particle mass is fixed as 5.9 × 10 7 M , the dark matter mass density profile can be compared with the galaxy number density profiles even though their absolute value can not be directly compared to each other. dark matter 2.0 ± 3.2 4.0 ± 1.4 1.0 ± 0.2 0.9 ± 0.8 galaxy N 3.4 ± 1.5 4.2 ± 2.5 0.3 ± 0.1 0.9 ± 1.4 galaxy M 1.9 ± 2.3 3.1 ± 1.8 0.9 ± 0.6 0.8 ± 2.3 gas 3.5 ± 4.8 4.0 ± 1.1 0.6 ± 0.1 0.7 ± 0.3 intracluster stars 2.7 ± 4.3 5.9 ± 2.0 2.1 ± 0.3 0.8 ± 0.8 Figure 6. Projected density profiles of dark matter, galaxy N, galaxy M, gas, and intracluster stars. Individual density profiles of 426 clusters with three different projections are over-plotted with the gray lines, while their averaged profiles are shown together as the colored lines. The projected mass density profiles for dark matter, gas, and intracluster stars are in a unit of M kpc −2 . The projected galaxy number density profiles, galaxy N and M, are scaled up by a factor of 10 8 to be shown in the same y-axis range of the other profiles. Figure 7. Averaged density profiles and their fitted lines with the double-power law. The overall density profiles are well described with the double-power law whose breaks are commonly in 0.7 − 1.0Rvir. Thus, the density profiles can be approximated as two separate power-laws with slopes of γ and β for R/Rvir < 1 and R/Rvir > 1, respectively.
We fit the averaged density profiles with a doublepower law of where r 0 is a scale radius, ρ 0 is a central density. Depending on r/r 0 , ρ can be approximated as follow: where α controls the sharpness of the break between the β and γ regimes. The double power-law function has been frequently used for describing the density profile of the dark matter structures formed in a cosmological context (Hernquist 1990;Zhao 1996;Navarro et al. 1996). It should be noted that the density profiles in  Figure 7 shows the averaged density profiles and their fitting results, which are summarized in Table 1. We adopt the standard deviation of the 426 profiles for each component as 1-sigma uncertainty for the Levenberg-Marquardt least-square fitting process using MPFIT-FUN/IDL routine (Markwardt 2009). Since the fitted profiles for the individual components have similar breaks around r 0 = 0.7 − 1.0R vir , they can be approximated by a broken power-law bent at similar r 0 . Interestingly, the Galaxy M profile shows the least conspicuous break. Its inferred that is from different radial profiles of the high-and the low-mass galaxies inside the cluster region. Because the higher (or lower) massweight is more frequently given for galaxies locating at the center (outskirt) region, the slope for the Galaxy M profile for R/R vir < 1 can be enhanced compared to that of the Galaxy N . In the inner region of R/R vir < 1, the dark matter density profile can be approximately described as a power-law of γ = 1.0 ± 0.2; the galaxy M gives γ = 0.9 ± 0.6, the most similar to that of dark matter. The dark matter profile in the outer region of R/R vir > 1 can be described as a power-law with β = 4.0 ± 1.4 despite the scatters in the outskirts of the halo. The β for the dark matter profile shows no significant difference from that of the galaxy M, β = 3.1 ± 1.8. The gas and the galaxy N components, respectively, give β of 4.0±1.1 and 4.2±2.5, the most similar to that of the dark matter. Similar to dark matter, the gas shows a significant departure from the fit at the end of the profile. Considering the fitted slopes and their uncertainties, the galaxy M profile described by a double-power law appears most consistent with that of dark matter over the entire radial range. The profile for Galaxy M shows the least conspicuous break. We do not have comprehensive explanation for this, but can at least comment on the comparison of the profile between galaxy M and N. The inner slope of the profile for galaxy M (i.e. γ) is slightly larger than that for galaxy N, which may be because of the amplified mass weight for the BCG and/or of the high-mass galaxies segregated into the inner region. On the other hand, the difference in the outer slope of the two profiles is not statistically significant given the error of each profile.

2D Spatial Distributions
We compare the global morphology of the projected density distributions using the ellipse fitting results of ellipticity, position angle, and centroid offset. Before comparing the ellipse fitting results of three different components (i.e. galaxies, gas, and intracluster stars) with that of the dark matter, we first choose the representative galaxy map from galaxy N, galaxy M, and galaxy σ v , whose ellipticity is best correlated with that of dark matter. Figure 8 shows the correlations of the measured ellipticities of the maps between dark matter and galaxies (i.e. galaxy N, the galaxy M, and the galaxy σ v ) at three different radii of 0.5, 1, and 2R vir . We estimate the significance of the correlation using the Pearson correlation coefficient ρ, standard deviation from one-toone correlation σ, and fitted slope µ, which are shown in each panel. The ellipticity of each galaxy map at R = 2R vir shows the strongest correlation with that of dark matter with the highest ρ and the lowest σ values. Particularly, the ellipticities of the galaxy N and the galaxy M show strong correlations with that of dark matter with ρ = 0.92±0.03 and 0.93±0.03, respectively. The correlation for the inner region of R ≤ 1R vir is getting weakened for all the cases not only by an intrinsic scatter coming from the different distribution between dark matter and galaxies (e.g. astrophysical effects including dynamical friction or baryonic feedback, Massey et al. 2015;Springel et al. 2018), but also by a numerical error coming from the small number of pixels used for the ellipse fitting. The ellipticity measured by the galaxy N and the galaxy M still show strong correlations with that of dark matter even at R = 1R vir with ρ = 0.79 ± 0.05 and 0.81 ± 0.05, while moderate correlations at R = 0.5R vir as ρ = 0.49 ± 0.09 and 0.51 ± 0.08, respectively. In summary, the case of galaxy M shows the highest ρ values at all radii among the three galaxy maps. The galaxy σ v case differs from other cases significantly at R = 2R vir . Hereafter, we use the galaxy M as a representative of the galaxy maps when we compare the spatial distribution of three different components of galaxies, gas, and intracluster stars with that of dark matter in and around galaxy clusters.

Comparison among Different Components
The ellipticity of dark matter at three different radii of 0.5, 1, and 2R vir is compared with those of the galaxies, gas, and the intracluster stars in Figure 9. At R = 2R vir , the ellipticity of the gas component shows the most prominent correlation with that of dark matter as ρ = 0.96 ± 0.01, which is significantly larger than that of intracluster stars (ρ = 0.74 ± 0.06) and slightly larger than that of galaxies (ρ = 0.93 ± 0.03). In the inner regions of R ≤ 1R vir , the ellipticity measured for the gas component is getting biased to a lower value than that of dark matter due to 'virialization; kinetic energy of infalling gas is transformed to thermal energy (to be discussed in Section 4). Because gas distribution is getting circularized in the inner region, the fitted linear-slope for the ellipticity correlation is significantly reduced to µ = 0.43 ± 0.03. Compared to the galaxy and the gas components, the ellipticity of the intracluster stars shows a weaker correlation with that of dark matter in general, except for the inner region. However, the difference of ρ values among different components at the inner region is negligible, less than 1-sigma error. The radial difference of ρ for the intracluster stars is relatively less prominent than those of other components. This seems to be because the BCG is not excluded, which makes a slightly better ρ value in the inner region. In summary, among the maps based on different components, the ellipticity of the galaxy M map shows the best correspondence to that of the dark matter map.
We now examine the difference of position angle (shortly, δPA) of the measured ellipse between each component and dark matter (see Figure 10). Because the PAs of different components at different radii do not show a systematic difference to that of dark matter, we show only the distribution of the scatter, so-called δPA, rather than the PA correlation of different components. Here, the δPA value closer to zero means that the ellipse's major axis of a component is more aligned to that of dark matter. The δPA distributions peak around zero (−0.5 • < δPA < 0.5 • ) with different standard deviation σ δPA for different components and radii. The position angles of all the components show better alignments with those of dark matter at the outer region where the ellipse fit is less affected by numerical errors arising from the ellipse fitting. At R = 2R vir , galaxies show the smallest σ δPA , but the difference from the other components is not so significant. We have perform the Kolmogorov-Smirnov test (K-S test) to examine how different the δPA distributions of galaxy M , gas, and ICS at R = 2R vir are. The K-S test for the δPA distribution of galaxy M and gas gives p-value of 0.45, indicating the null hypothesis that the two distributions are drawn from the same population cannot be rejected. The pvalues from the K-S test between galaxy M and ICs, and between gas and ICS are both < 0.001, suggesting a significant difference between the two samples. Figure 11 shows how the centroid of the ellipse measured for each component coincides with that of dark matter, which is measured by an offset of the centroid for each component from that of dark matter. Because the area of the annulus centered by the centroid of dark matter varies with radius, the count for each annulus changes with radius. Thus, to remove the difference introduced by this in the centroid-offset histogram, we normalize them by the area of the annulus. The normalized centroid offset distribution is highly peaked around zero, and the probability decreases exponentially with the centroid offset. The deviation of the normalized histogram from zero (σ offset ) is measured to quantify the centroid consistency between each component and dark matter. The σ offset for each component shows an insignificant change and/or trend for different radii. For a given radius, the gas component shows the lowest σ offset , which is different from other components by 2-sigma at least (to be discussed in Section 4).

Dependence on Dynamical Status of Galaxy Clusters
To investigate how different dynamical states of galaxy clusters affect the correlation of global morphology between dark matter and other components (e.g. ellipticity), we divide galaxy clusters into two sub-samples using the virial ratio, which can be derived as where K is a total kinetic energy, U is a total potential energy. Here, E s is the energy from surface pressure at boundary, which is introduced to correct an assumption that the system is in complete isolation (Chandrasekhar 1961). Because a halo is embedded in a cosmological density field, particles close to the halo are continuously accreted across the halo boundary. To account for these particles giving pressure on the halo, the conventional virial theorem should be corrected by the additional term, so-called the surface pressure. We follow an ap-  proximated E s formula of Shaw et al. (2006) and Cui et al. (2017) for collisionless and collisional components, respectively. Figure 12 shows the virial ratio distribution of 426 clusters. Despite the surface pressure correction, the Q value peaks at 1.05, which is slightly higher than 1. Although we could not completely understand the cause of this difference, our Q values are still good enough to construct subsamples for the comparison: i.e. virialized vs. unvirialized. We choose the 200 closest clusters to the peak Q value and define them as the 'virialized' sub-sample, while the 200 farthest clusters as the 'unvirialized' sub-sample 4 . Figure 13 shows how the ellipticity correlation shown in Figure 8 is different for the two sub-samples. The 'virialized' sub-sample shows a marginally higher ρ value than the 'unvirialized' sub-sample in all panels, but the difference is not so significant in most cases: only the outer region of R = 2R vir shows a bit larger difference than the 1-sigma error. To statistically quantify the similarity of the ellipticity correlation between the 'virial- 4 We perform a similar analysis for two 'unvirialized sub-samples: 86 clusters with 'Q < Q peak ' and 114 clusters with 'Q > Q peak ', where Q peak is the peak of Q distribution. Main results for the two 'unvirialized sub-samples do not show any significant difference. A detailed discussion on the implication of the unvirialized clusters with different Q values is beyond the scope of this paper. ized and 'unvirialized sub-samples, the two-dimensional KS-test is performed in each panel. The ellipticity distribution between gas and dark matter components at R = 1 and 2R vir shows a significantly low p-value as 0.02 and 0, respectively, which indicates the non-negligible difference between the two sub-samples. Note that the ρ values can be used to compare how the ellipticity correlation from the two sub-samples are alike to each other, while the KS p-value is for the similarity of the ellipticity distributions between the two sub-samples, regardless of the ellipticity correlation. Similar to the case of ellipticity, we examine how the two sub-samples show a difference for the δPA and the centroid offset distributions (see Figures 14 and 15). We find no statistically significant difference of the δPA distribution between the two sub-samples not only for all the components but also for all radii. In case of the centroid offset distribution, σ offset of the 'virialized' subsample at R = 2R vir is significantly lower than that of the 'unvirialized' sub-sample for all the components. It indicates that the centroid of ellipse measured at the outer region by each component is closer to that of dark matter for the 'virialied' sub-sample, compared to the 'unvirialized' one.

DISCUSSION
It has been well known that galaxies are good tracers of matter field in cosmologically large scales, which is still in linear regime (Kaiser 1984;Desjacques et al. 2018). However, non-linear evolution in smaller scales can introduce a scale-dependent bias, which can not be easily predicted from a theory. Using the IllustrisTNG cosmological hydrodynamic simulations, Springel et al. (2018) have studied the non-linear galaxy bias over a wide range of scales. In this study, we use the Il-lustrisTNG simulation and quantitatively compare the global spatial distribution of dark matter not only with that of the galaxies but also with those of gas and intracluster stars to better understand how well galaxies and other components (i.e. gas, intracluster stars) trace dark matter in and around clusters.
We use the ellipse fitting method to reflect the global morphology of the galaxy clusters, and compare the ellipticity of the dark matter density map with that of three different galaxy maps of galaxy N, galaxy M, and galaxy σ v , which have been widely used in observations to infer underlying matter distribution in cluster region Zahid et al. 2016;Liu et al. 2018;Sohn et al. 2020). We find that the ellipticity of dark matter is better correlated with that of galaxy N or M rather than galaxy σ v over the entire radii. It is inferred that a small number of galaxies used for measuring σ v may not show the underlying dark matter distribution well. The measured ellipticity of the galaxy N and M and their correlation with that of dark matter are similar to each other (see Figure 5 and 8). Between the galaxy N and M, the ellipticity measured from galaxy M (massweighted galaxy number density) shows a slightly better correlation with that of dark matter even though the difference is not significant. This is consistent with the expectation from the studies of large-scale structures focusing on the reconstruction of dark matter distribution with the halo/galaxy-mass based techniques (e.g., Wang et al. 2009, see also Hong et al. 2021 for the importance of peculiar velocity information in the reconstruction).
We then use the similar ellipse fitting method to compare the global morphology of different components in galaxy clusters (i.e. galaxies, gas, intracluster stars). The three components show the strongest ellipticity correlation with dark matter at the outer region, and their correlations are getting weakened for the inner region. The larger scatter in the inner region seems to originate from the following two factors: 1) an intrinsic error coming from a fact that each test particle of galaxies, gas and intracluster stars tends to fail to trace the dark matter distribution at the inner region where the non-linear evolution takes place intensively (e.g. astrophysical effects including dynamical friction or baryonic feedback, Massey et al. 2015;Springel et al. 2018), and 2) a numerical error originated from the smaller number of pixels for the inner contours, which results in the larger uncertainty on the ellipse fitting. To avoid the numerical error, a non-parametric comparison of contours without the fitting, such as the cross-correlation technique  In the outer region, the ellipticity of gas distribution shows the strongest correlation with that of dark matter even though the difference from the other components is not significant. The higher ρ value of the gas component than that of the others could be because the gas component has a higher mass-resolution than other components so that the gas can describe the density distribution in a continuous way, similar to the case of dark matter. The gas ellipticity at the inner region, however, is getting biased to a lower value than that of dark mat-ter, which is also shown as circular shape at Figure 5. It is because infalling gas particles to the galaxy cluster are heated by a virial shock, and thus their kinetic energy is transformed to thermal energy: so-called 'virialization'. The biased gas ellipticity shows the effect of the gas virialization only for the inner and middle regions (0.5 and 1 R vir ), while the effect has been believed to begin to occur at much further distance around 2-3 R vir (Zinger et al. 2018). Although the ellipticity is biased to the lower value at the inner region, the gas component shows the smallest centroid offset with respect to dark matter's over the entire radii. It is inferred that gas particle can quickly adapt to the centrally concentrated mass dis-  tribution by collisions between themselves with much shorter timescale than those of other components. Separation between gas and dark matter of merging galaxy clusters is reproduced in the recent simulations of Moura et al. (2021), which is similarly shown in Figure 14; the centroid of each component from the 'unvirialized subsample is less aligned to that of dark matter. See the following paragraphs for details.
The diffuse intracluster light, which is not bounded to individual galaxies but to the cluster potential, is expected to trace the global matter distribution of clusters (Montes & Trujillo 2019;Alonso Asensio et al. 2020;Sampaio-Santos et al. 2020). The intracluster stars in this study include the intracluster light and the BCG by excluding individual galaxies that are masked from the original stellar density map. However, the ellipticity of the intracluster star shows a weak correlation with that of dark matter compared to the other components in general. It is inferred that 1) the low-mass galaxies, not found with the SUBFIND algorithm (< 1.2 × 10 9 M ) or whose 2R h are smaller than a half of pixel size, are missing from the masking and thus lead to clumpy structures especially at the outer region where the low-mass galaxies are more abundant, and 2) the detailed tidal features of intracluster stars, which can trace interaction between galaxies, result in scatters on the ellipticity correlation since global spatial distribution is described by an ellipse only.
Compared to the 'unvirialized sub-sample, the 'virialized galaxy clusters show a better correlation of spatial distribution between dark matter and other components at all radii even though the difference in the correlation between the sub-samples is not significant; the marginally higher ρ value for the 'virialized subsample can mean that each component follows better the spatial distribution of dark matter some time after merging events. The two-dimensional KS-tests for the ellipticity distributions of dark matter and gas show that there is a non-negligible difference between the 'virialized and 'unvirialized sub-samples only at R = 1 and 2R vir . Since the 'unvirialized' features can remain in the smaller-scale structures that may be missed by the the ellipse fitting, a different comparison focusing on detailed features would be necessary to better understand the difference between the sub-samples.
Interestingly, the fitted linear-slope (µ) of ellitpicity correlation between dark matter and gas of the 'virialized' sub-sample at R = 0.5 and 1R vir shows a slightly lower value than that of the 'unvirialized' sub-sample. It means that the gas ellipticity of the 'virialized' galaxy clusters tends to be more biased to a lower value than the ellipticity of dark matter, compared to that of the 'unvirialized' galaxy clusters. Because the 'virialized' clusters tend to be growing in a continuous way by a smooth accretion, while the 'unvirialized' clusters in a radical way by quite recent mergers, the smaller µ value for the 'virialized' sub-sample can be circumstantial evidence that gas residing longer inside the galaxy clusters is heated more by a virial shock and thus circularized more.
Unlikely to the insignificant difference of the σ δPA between 'virialized and 'unvirialized sub-samples for all components at all radii (see Figure 14), the σ offset at R = 2R vir shows a significant difference (see Figure 15), where the centroid of each component from the unvirialized sub-sample is less aligned to that of dark matter. This result is consistent with the observations showing that X-ray peak and the BCG tend to be separated more for non-cool-core clusters, which are believed to be associated with recent major merging events and thus tend to be unvirialized more (Katayama et al. 2003;Sanderson et al. 2009;Hudson et al. 2010;Hoffer et al. 2012;Kim et al. 2017). At R = 0.5 and 1R vir , on the other hand, the σ offset between the two sub-samples does not show a significant difference. The non-parametric comparison (i.e. the Modified Hausdorff distance method, Montes & Trujillo 2019) should be applied to check whether the insignificant difference is originated from an intrinsic error coming from a fact that spatial distribution at the inner and middle regions is less sensitive to the recent merging events.

SUMMARY
We use the IllustrisTNG cosmological hydrodynamical simulation to understand how well dark matter distribution in and around galaxy clusters is traced by galaxies, gas, and intracluster stars. Using the ellipse fitting method, we quantify the global morphology of dark matter distribution of galaxy clusters and compare the result (e.g. ellipticity) with those of galaxies, gas, and intracluster stars. Our primary results are as follows: (1) Among three different galaxy maps of galaxy N (number density), galaxy M (mass-weighted number density), and galaxy σ v (velocity dispersion), ellipticity of the galaxy M shows the best correlation with that of the dark matter map over entire radii.
(2) The ellipticities of three different density maps from galaxies (galaxy M ), gas, and intracluster stars are compared with that of dark matter. We find that ellipticity of the dark matter map is best reproduced by that of the galaxy map.
(3) The ellipticity of the gas component shows the most prominent correlation with that of dark matter at the outskirt of galaxy clusters (R = 2R vir ). However, in the inner region of R ≤ 1R vir , it is significantly biased to a lower value than that of dark matter due to 'virialization.
(4) We have examined how position angle and centroid of the measured ellipse for each component coincides with that of dark matter: δPA and centroid offset, respectively. We find that δPA and centroid offset distributions do not show an insignificant difference and/or trend for different components and radii.
(5) The 'virialized galaxy clusters show a better correlation of spatial distribution between dark matter and other components than the 'unvirialized clusters. This can suggest that it requires some time for each component to follow the spatial distribution of dark matter after merging events.
This study demonstrates that galaxies are still good tracer of dark matter distribution in and around galaxy clusters, similar to the case where galaxies well trace the matter distribution in cosmologically large scales. The galaxy density map with mass weight works best as expected from observations that total stellar masses of cluster member galaxies could be a good proxy for cluster masses (e.g. Pereira et al. 2018;Palmese et al. 2020) and that the mass-weighted galaxy number density maps of clusters show good correspondence to the weak-lensing maps (e.g. Fig. 5 of Hwang et al. 2014). It would be interesting to make a comparison as in this study, but focusing on the detailed features that might have been missed in our method (e.g. Yoo et al., in preparation). Also, this study will be very helpful for understanding the results of the mass distribution in and around galaxy clusters from wide-field surveys with various tracers, which include DESI (galaxies; DESI Collaboration et al. 2016), eROSITA (hot gas; Predehl et al. 2021), and the Vera C. Rubin Observatory (dark matter and intracluster stars; LSST Science Collaboration et al. 2009).