Comparison of spatial distributions of Intracluster light and Dark Matter

In a galaxy cluster, the relative spatial distributions of dark matter, member galaxies, gas, and intracluster light (ICL) may connote their mutual interactions over the cluster evolution. However, it is a challenging problem to provide a quantitative measure for the shape matching between two multi-dimensional scalar distributions. We present a novel methodology, named the {\em Weighted Overlap Coefficient (WOC)}, to quantify the similarity of 2-dimensional spatial distributions. We compare the WOC with a standard method known as the Modified Hausdorff Distance (MHD). We find that our method is robust, and performs well even with the existence of multiple sub-structures. We apply our methodology to search for a visible component whose spatial distribution resembled with that of dark matter. If such a component could be found to trace the dark matter distribution with high fidelity for more relaxed galaxy clusters, then the similarity of the distributions could also be used as a dynamical stage estimator of the cluster. We apply the method to six galaxy clusters at different dynamical stages simulated within the GRT simulation, which is an N-body simulation using the galaxy replacement technique. Among the various components (stellar particles, galaxies, ICL), the ICL+ brightest cluster galaxy (BCG) component most faithfully trace the dark matter distribution. Among the sample galaxy clusters, the relaxed clusters show stronger similarity in the spatial distribution of the dark matter and ICL+BCG than the dynamically young clusters. While the MHD results show weaker trend with the dynamical stages.


INTRODUCTION
One of the most promising routes to understanding how a cluster assembles is to observe cluster quantities which are sensitive to the evolutionary state of the cluster. Recent deep observations of nearby clusters show distinct diffuse intracluster light (ICL), which is light from stars not gravitationally bound to any individual cluster galaxy (Feldmeier et al. 2002;Lin & Mohr 2004;Gonzalez et al. 2005;Zibetti et al. 2005;Mihos et al. 2005Mihos et al. , 2017Contini et al. 2014Contini et al. , 2019DeMaio et al. 2018;Ko & Jee 2018;Jiménez-Teja et al. 2019;Kluge et al. 2020;Furnell et al. 2021;Yoo et al. 2021).
With regard to galaxy cluster evolution, one may raise a question: how does the spatial distribution of dark matter vary with redshift, mass, and dynamical stage of the galaxy cluster? The gravitational weak lensing is widely applied to investigating the dark-matter distribution in clusters. Although it is powerful, the weak lensing method also has 2 Yoo et al.
several drawbacks such as the difficulty of identifying background galaxies and accurately measuring their shapes in the presence of the telescope point spread function. Therefore, it would be worthwhile to search for other quantities that can recover the dark matter distribution, which could provide complementary information to enhance our understanding.
Unlike the hot gas component, the ICL is collisionless, and at the same time gravitationally bound to the global potential of the galaxy cluster (not to individual galaxies). In the current ΛCDM paradigm, dark matter halos are built up by hierarchical assembly of smaller halos, and the ICL grows through the assembly history of galaxy clusters (Mihos et al. 2005;Mihos 2016). Accordingly, ICL is expected to map the spatial distribution of dark matter, and several previous studies have shown it does (Montes & Trujillo 2019;Alonso Asensio et al. 2020). Moreover, the ICL could be used as an effective tool to measure the evolutionary stage of galaxy clusters by assuming that the more evolved clusters would contain a richer ICL component, and that its spatial distribution would be well aligned with the distribution of dark matter. Making a rigorous statistical comparison of the spatial distributions of ICL and dark matter would be a crucial step in the study of ICL as a luminous tracer for dark matter, and as an evolutionary stage estimator for galaxy clusters.
The one-dimensional radial profiles of ICL and dark matter were compared in Sampaio-Santos et al. (2021). However, the azimuthally-averaged radial profile did not fully represent the actual spatial distribution in 2-or 3-dimensions and thus consequently lost some valuable information. Previous studies of the two-dimensional spatial distribution of ICL and dark matter Montes & Trujillo (2019); Alonso Asensio et al. (2020) used the Modified Hausdorff Distance (MHD) method (Dubuisson & Jain 1994), which quantifies the difference between two nearby contours using the distances between points on the contours. The cross-correlation technique was adopted by Hwang et al. (2014) to compare the dark matter maps from weak lensing analysis and galaxy cluster redshift surveys. Outside of the astrophysics field, the environmental niche overlap (Schoener 1968;Warren et al. 2008) is used in ecology, which shows the co-existence (overlap of resource use) of different species. The Earth mover's distance (Kantorovich 1960), also referred to as the Wasserstein metric, has gained recent popularity as a loss measure in certain machine learning applications involving image generation or classification.
In this particular statistical analysis, we compare two spatial distributions while the specific relation between their signal strengths is either not known a priori or is not being sought. In this work, we suggest that under such a specific requirement the aforementioned methods could suffer from biases and may not be the most appropriate choice. We thus introduce a new method for comparing the morphology of fields, and compare it with one of the previous methods, the MHD method.
Although our method may be generalizable to a range of problems here we will restrict our analysis to the quantification and comparison of simulated massive galaxy clusters. Furthermore we will investigate the ICL versus dark matter spatial distribution comparison and galaxy cluster dynamics. A similar comparison from our group, but focusing more on the large-scale spatial distribution in and around clusters based on a parametric method is presented in Shin et al. (in preparation).
In §2 we define the simulated data we will use. In §3 we introduce our new methodology and discuss its practical application. Our analysis on our method is presented in §4, where we directly compare our new method with the MHD method. We discuss about the results and possible applications in §5. We summarize in §6.

SIMULATION DATA
We applied our method to six galaxy clusters lying at different dynamical stages selected from Galaxy Replacement Technique (GRT, Chun et al. 2022) simulation, which is run with the Gadget3 code (Springel 2005) for the study of cosmological gravitational evolution of galaxy clusters. First, we performed the N-cluster Run, a low-resolution dark matter only simulation with 512 3 particles of a uniform mass (10 9 h −1 M ) in a cubic box of a side length, 120 h −1 Mpc. Then, we resimulate the N-cluster Run but in a higher-resolution mode (10 h −1 pc), placing a stellar disk with high-resolution star particles according to the exponential disk galaxy model (Dutton 2009). We replace only the dark matter halos that fall into the galaxy cluster and contribute to its growth. When a dark matter halo in the low-resolution simulation reaches its maximum mass as a central halo, and if the mass is greater than 10 11 h −1 M , the dark matter halo is replaced by a high-resolution dark matter halo and high-resolution stellar disk. Since the replacement happened before the halo enters the virial radius of the host cluster, we assume it does not suffer from significant tidal interactions and does not contribute to the BCG+ICL. As these high-resolution halos fall into the clusters and are tidally stripped, the BCG and ICL form naturally. For the GRT simulation, we use the post-Planck cosmological model of Ω m = 0.3, Ω Λ = 0.7, Ω b = 0.047, and h = 0.684. Please refer to Chun et al. (2022)  4 5 6 7 log 10 (M /kpc 2 ) 4 5 6 7 log 10 (M /kpc 2 ) 4 5 6 7 log 10 (M /kpc 2 ) Figure 1. Galaxy clusters simulated in the GRT simulations. Projected images in 2.5 × 2.5 h −2 Mpc 2 boxes for the dark matter (first column), all stellar particles (second column), member galaxies excluding BCG (third column) and BCG+ICL (fourth column) are shown. Throughout this work, the spatial distribution comparison between dark matter and galaxies is conducted using all member galaxies including the BCG, however we show satellite galaxies without the BCG as an aid the visualize the relative distributions of the galaxies and BCG+ICL. The mass growth histories of each cluster are plotted in the fifth column. The major merging events (mass ratio smaller than 1:5) are marked as red triangles. detailed description of the resimulation. The GRT allows us to trace the tidal interaction histories of stellar particles in galaxies and galaxy clusters. The high resolution of GRT simulations (low mass of star particles of 5.4×10 4 h −1 M ) allows us to resolve the low surface brightness features down to µ V < 32 mag arcsec −2 . This surface brightness limit is approximately three mag arcsec −2 lower than in Illustris TNG100 (Nelson et al. 2018). The brightest cluster galaxy (BCG), other member galaxies, and ICL are defined among the stellar particles using the six-dimensional phase-space friends-of-friends halo finding algorithm, ROCKSTAR (Behroozi et al. 2013). We define the stellar particles which belong to the BCG and other member galaxies. The ICL is then defined as stellar particles which belong to the galaxy cluster but not to any member galaxies. We show the projected maps of each component in Figure 1. From the first to the fourth column, maps of the dark matter and each component are shown.
In the fifth column, the mass growth histories are shown, with the indication of the major merging events. Considering the formation time (z m/2 ; redshift at which 50% of the cluster mass is obtained) and the redshift when the last major merger occurred (z LM M ), three dynamically young clusters (cluster 4, 5, and 6) and three dynamically relaxed fossil clusters (cluster 1, 2, and 3) were defined. Clusters 4, 5, and 6 have relatively recent major merging events. Specifically, cluster 6 seems to currently be undergoing major merging. Meanwhile, clusters 1, 2, and 3 had their last major merger events at much earlier times. Their formation times (z m/2 , see Table 1) were also relatively early. We projected the simulated galaxy clusters in three different viewing angles (i.e., the x − y, x − z, and y − z planes) to avoid biases arising from any one particular projection, augmenting our data sample. The projected images were cropped to a scale of 2.5 × 2.5 h −2 Mpc 2 , which were then pixelized to 1024 × 1024 pixels. Thus, the pixel scale for the images is 2.44 h −1 kpc/pixel.

WOC method
We introduce a novel method to quantify the similarity of two 2-dimensional spatial distributions. The similarity is quantified by the Weighted Overlap Coefficient (WOC) whose functional form is given in equation 1. In this, we calculate the fraction of overlapping area between two distributions at various density threshold levels (see Figure 2). Using a weighted sum of those fractions, we normalize the similarity measure between 0 and 1.
The calculation proceeds as follows: 1. Prepare two-dimensional maps of two components smoothed on the same angular or spatial scale.
2. Determine the signal center in the reference map (for example, the peak position of the dark matter density map) and measure the one-dimensional radial profile. Find the signal levels at n equally-spacing radii (for example, if n = 3, r = 100, 200, and 300 kpc).
3. Find the enclosed area at each level of the isolevel contour (i.e. pixels with signal values larger than the input signal levels). When there are multiple disconnected contours in the map or multiple isolated areas, simply add up all the areas.
4. In the comparison map (for example ICL), find contour levels to have the same enclosed area size as those found in the reference map. We do this by starting from the peak signal and lowering the pixel value threshold in small steps until the area constraint is satisfied. This step is needed to quantify the overlapping area as a percentage.
2D compare 5 Figure 2. WOC method idea. We measure the overlap area of contours between two distributions at various density threshold levels. We calculate the fraction of the overlapping area and give weight using the signal strength to quantify their overall similarity as a normalized number between 0 and 1. See more details in Section 3.1.

5.
Measure the amount of overlapping areas between the reference and comparison maps at each signal level and calculate the overlapping fraction (see Figure 3).
6. Parametrize the similarity of the two component maps with n overlapping fractions multiplied by weightings as (1)  Finally the similarity of the two spatial distributions is quantified as the WOC, which is a single number between 0 and 1. For two maps I 1 and I 2 , the WOC is defined as where a i is the overlap area (%) for the i th level, w A,i is the area ratio of the i th contour compared to the largest (n th ) contour, w ρ1,i is the density level at the i th contour on the Map 1 and w ρ1,i is the density level at the i th contour on the Map 2. Throughout this weighting system, we take account of both the area of each contour, which are used to measure the overlapping region, and the strength of the signal of the selected contours of each map (i.e., dark matter and ICL). Among the contours on various levels, the contour for the stronger signal level would occupy smaller area than that of the weaker signal. Moreover, if the distribution is rough, the contour of the stronger signal would be smaller than the smooth distribution situation at the same level. If the sharp peaks of two distributions match well with each other, even though the chance of it would be smaller due to the small area, then we could regard those cases as the location of two distributions coincide. Therefore, using the weight w A (area ratio), we give more weight to smaller contours. The contour areas of two distributions for the same level to compare are set to be same, thus w A would represent the signal strength of both distributions. The signal strengths of the individual distributions are taken into account . WOC method measuring process for cluster 1, simulated in GRT simulations. All panels have a box size of 1.5 × 1.5 h −2 Mpc 2 . The dark matter distribution in cluster 1 (first row) was compared with all stellar particles (second row; cluster galaxies including BCG+ICL), member galaxies (third row; cluster galaxies including BCG), and BCG+ICL (fourth row). For three lower panels, the distributions of dark matter and the other compared components are marked with blue and orange colors, respectively. We measured the overlap area for three levels, in the order from highest to lowest, as the first (left), second (middle), and third (right) levels, and gave them weight to calculate the WOC. . The intermediary steps used for measurement with the WOC (third row) and MHD (fourth and fifth row) methods, for cluster 6 (third and fourth row) and cluster 1 (fifth row). The raw data images of dark matter and BCG+ICL of cluster 6 are plotted in the first and second row, respectively. All of the panels have a box size of 1.5 × 1.5 h −2 Mpc 2 . For both methods, the distributions of dark matter and BCG+ICL are marked with blue and orange colors, respectively. The three levels of dark matter and BCG+ICL distributions are plotted in the left (first level), middle (second level), and right (third level) panels. In the MHD method, when there were several contours for one level (marked by the dotted line), the largest contour (solid line) was selected to compare with the contour of the other component. more directly/ separately, via their density value w ρ1 and w ρ2 . If the distribution maps are density fields, w ρ1 and w ρ2 would be simply the pixel value (before normalization in equation 2) at the contour to compare on the Map 1 and Map 2, respectively.
Note that, we are not interested in the relative signal strengths of each component or the exact shape of their profiles, rather we focus on the spatial correspondence of the two components. One of the main advantages of this method is that it performs well when comparing two distributions when either or both contain disconnected regions. Furthermore, in the WOC method we do not need to compute the individual contours, and thus it suffers less bias even when working with masked maps.

MHD method
Previous studies (Montes & Trujillo 2019;Alonso Asensio et al. 2020) of the two-dimensional spatial distribution of ICL and dark matter have used the MHD method.
The MHD defined by Dubuisson & Jain (1994): where The two sets of points, X = {x 1 , x 2 , ..., x Nx } and Y = {y 1 , y 2 , ..., y Ny }, define two contours, and · is the Euclidean norm. To compare the MHD between different clusters, we adopt the relative MHD defined by Alonso Asensio et al. (2020): for the distance r, which also in turn defines the contours, X r , Y r , via the angle averaged one-dimensional radial profile. Situations can arise where a level defined by r results in more than one contour, i.e., multiple disconnected regions (see Figure 4). In this case we follow the suggestion of Alonso Asensio et al. (2020) and choose only the largest contour when evaluating equation 5. Thus, this method may not be robust in cases where the spatial distributions exhibit multiple local sub-structures (disconnected from the main body), which are expected in dynamically active clusters, such as Bullet or Coma cluster. For example, the previous study using the MHD method (Alonso Asensio et al. 2020) excluded samples from their analysis, if there are more than one central galaxy in the cluster.

Modified WOC method and Characteristics of different methods
What does it mean for one distribution to trace another distribution? In some cases we may want to quantify both the spatial co-existence and the corresponding signal strength together. The WOC method, however, focuses solely on the spatial mismatch between two distributions, and is not sensitive to the relative strength difference within distributions, such as the shape of the radial profile. This is because we set the contours using the reference distribution, and find the corresponding contours with the same areas from the comparison distribution.
We developed another WOC system (WOC M ), where the contours from the comparison distribution are not found using the same area, but rather by the same mass fraction. For example, under the WOC M method, if the contour of the first level for the dark matter distribution encloses 10% of the total mass of dark matter, it will be compared with a contour which encloses 10% of the total mass of the comparison component. In this case, the overlap area a i (%) in equation 1 will be calculated as the overlapped region area divided by the total area occupied by the contours. Here, the signal strength is already taken into account by setting the contours, thus, among the three weight factors in equation 1, we have only one weight factor, which is the area ratio, w A,i .
We generated mock images that followed NFW profiles, while varying four parameters; the concentration rate, offset, ellipticity, and position angle (PA). We then checked how the different spatial distribution comparison methods reacted to the various situations (see Figure 5 and Table 2). Fixing the four parameters for the reference distribution (Map 1 in Table 2   The original WOC method (WOC A ) and the MHD method (ζ M HD ) showed similar results, which were insensitive to the concentration rate, but were sensitive to the other three parameters changes. The results of the modified WOC method (WOC M ) dropped as the concentration rate in Map 2 differed from Map 1. The WOC M method is highly sensitive to the concentration rate (radial profile) change, so that it could not detect the offset, ellipticity, and PA difference if the comparing distribution had a different concentration rate (the red lines in the middle row in Figure  5). Thus, the results of the modified WOC method could exhibit a degeneracy between the situations where cluster components have different profiles or where they are spatial mismatched, in the low WOC M case. However, the WOC M method has an advantage over WOC A method in that it compares contours with the same mass ratio (more physical meaning assigned to the pair of contours) and has a simplified weighting system. Depending on their scientific interest, ones could choose an appropriate method.

Measurement of Overlap Area before binning and weighting
The WOC method measures by choosing a certain binning in the radial distance from the center of the image, measuring the overlap area percentage and giving weight. Without binning and weighting, we re-arranged the pixels of each distribution in the order of pixel value, and we then measured the overlap area while increasing the selected brightest area pixels. For example, we checked the location of the brightest 100 pixels in Map 1 and Map 2, and counted how many pixels resided in the same position. This measurement was repeated with increasing brightest area (see Figure 6). The maximum area was 1024 × 1024 pixels. The test was conducted for various image convolution factors (upper panel), various galaxy cluster components (middle panel), and various galaxy clusters with different dynamical stages (lower panel).
Convolution factors: We vary the convolution scale in increments of 2 pixels from 0 upto 18 pixels. Depending on the convolution factor the behavior of the overlap area can change significantly, although at larger scales the curves become similar, as we might expect, since convolutions are local transformations. Because Map 2 (in this case, the BCG+ICL of cluster 1) has pixels with a pixel value of zero, the overlap area for the non-convolved (the lime green line in the upper panel) and the lightly convolved images saturate at a certain point.
Cluster components: For various components of an example galaxy cluster (cluster 1), we measured the overlap area with the dark matter distribution by increasing the brightest area (middle panel). The ICL only component (magenta dotted line) had zero overlap until the 200 brightest pixels. Apart from the ICL only component, all stellar particles (BCG + other member galaxies + ICL), all member galaxies, and the BCG+ICL components showed a high overlap area percentage until around 4000 brightest pixels. As the brightest area grows (lowering the pixel value), new islands may emerge for each component and if they are non-overlapping the overall percentage of overlap area will decrease, however this may recover again, resulting in the small dips in we see in Figure 6. Between the 5 × 10 3 and 5 × 10 5 interval, the BCG+ICL component (red dash-dotted line) showed a higher overlap than the stellar particles (blue solid line) or the member galaxies (yellow dashed line). Due to the fixed size of the images, and the fact that the faintest pixels all have a value of zero, the overlap area rise up again after the brightest area of 10 5 pixels.
Dynamical stages: The overlap areas of the BCG+ICL and dark matter for the galaxy clusters with various dynamical stages are shown in the lower panel in Figure 6. Overall, the more relaxed clusters with earlier formation times (clusters 1, 2, and 3) had higher overlap than the dynamically younger clusters with later formation times (clusters 4, 5, and 6). For the case of cluster 6, the peaks of dark matter and BCG+ICL did not match, which caused a zero overlap until the 500 brightest pixels. Moreover, apart from cluster 6 (with a large offset between the dark matter peak and the BCG), the more relaxed clusters (clusters 1, 2, and 3) exhibit a flatter trend than the dynamically younger clusters (clusters 4 and 5) between the area interval corresponding to the 100 kpc and 200 kpc contours (see the vertical gray dashed lines in the lower panel of Figure 6). The ratios of the overlap area calculated for small and large contours (overlap area (%) at 200 kpc contour/ overlap area (%) at 100 kpc contour) are 0.89, 0.88, 0.88, 0.77, 0.68, and 1.25 for cluster 1, 2, 3, 4, 5, and 6, respectively. This may be expected from the fast collapse and slow growth scenario (More et al. 2015;Fujita et al. 2018), where both relaxed and unrelaxed system have well matched distribution between dark matter and BCG+ICL in the inner radius, while only relaxed systems have had enough time to virialize in the outer regions. Thus, this ratio or the steepness of the slop in the lower panel of Figure 6 could be used as an indicator of the dynamical stage. Even cluster 6 could be an example use case for this indicator, showing that the ratio of overlap area fractions for small and large contours greater than one could indicate that the cluster is undergoing extreme merging events.  A smoothing factor of 17 is used for the decenter and binning tests. The labels are in the order of formation time (z m/2 ), where cluster 1 formed early (dynamically relaxed) and cluster 6 formed late (dynamically young). Note that, for spatially matched distributions, the WOC would be close to 1, and the MHD would be close to 0. The range of the MHD result is indefinite, i.e., the ζ-value is unbounded for the case of non-matching distributions. We thus take the logarithm of ζMHD, to allow for a qualitative visual comparison between the tends of the MHD and WOC methods. The MHD method results using the largest contour are plotted as solid lines, while using all contours are plotted as dashed faint lines.

Robustness Test
An ideal estimator should reflect the desired property that is being investigated (in this case, show a clear trend between WOC and the dynamical stage of galaxy clusters) and at the same time, be robust under data processing parameters such as smoothing scale, centroid, bin size, etc. In this section, we test the robustness of the WOC method while varying the convolution factor for the images, the decenter distance, and number of levels in the measurement of the WOC using GRT simulations galaxy clusters, and compare the results with those of the MHD method.
In the MHD method, the procedure choosing the largest contour for each level and excluding all sub-structures (see Section 3.2) could produce larger error, thus we compute both cases, i.e., using the largest contour and all contours.
Note that, a direct comparison between the results of WOC and MHD is difficult because the range of the MHD result is indefinite (see the y-axis of right panel in Figure 7), whereas the WOC result is normalized between 0 and 1 (see the y-axis of left panel in Figure 7). Thus, it is only appropriate to focus on the relative difference between samples within each method individually and look for systematic trends. We then may safely concluded the relative robustness (or biases) of each method under the different observational circumstances mentioned above.
Smoothing: We convolved the images using the Gaussian2DKernel astropy package, where the full-width halfmaximum size of the gaussian kernel is indicated as the smoothing factor. As the smoothing scale applied to the input images increases, the resultant WOC also increases because the convolution fills the gaps between pixels and results in higher overlap fractions. And then, it becomes flat at around a smoothing factor of 16 ∼ 18 (see the upper left panel in Figure 7). On the upper right panel in Figure 7, the same test was conducted with the MHD method using the largest contour for each level (solid line) and all contours (dashed faint line). To express the MHD result (which is a function of radius) as a single number, we calculated the relative MHD (ζ M HD ) by taking the mean of ζ values for the galaxy cluster. To compare the trend of the MHD result with the WOC result, we took the logarithmic scale of ζ M HD . As the smoothing factor increased, the relative MHD value using the largest contour changed slightly, except for cluster 6. For cluster 6, which undergoes merging, the relative MHD value fluctuated severely (see the blue solid line in the upper right panel in Figure 7). This symptom was alleviated as we take account all contours (dashed faint line), however the overall results become less robust.
Decentering: The first step in calculating the WOC or the MHD is to determine the center of the galaxy cluster and compute a radial profile. Depending on the choice of the cluster center, the contours could be altered, which could cause differences in the results for the WOC or MHD. We tested the effect of decentering on the WOC/ MHD result (see middle panel of Figure 7). Both the WOC and MHD result seemed to be relatively robust compared to other parameter changes. Taking all contours for the MHD calculations, the difference between cluster 6 and other clusters becomes milder, however the overall scatter becomes larger.
Binning: The choice of binning for measuring the overlap area at each level (for example n = 3 and r = 100, 200, 300 kpc) is somewhat arbitrary, and could affect the final result. Considering the observational studies of galaxy clusters, we set the outermost radius as 300 kpc, and divided it into n equally spaced bins ranging from n = 3 ∼ 12 (see the lower panel in Figure 7). Both the WOC method and MHD method showed consistent results with the binning change, with the exception of the MHD result for cluster 6. As the number of levels for measuring the overlap area increased, the ζ M HD of the dark matter and BCG+ICL component in cluster 6 rapidly increased (contours became mismatched).
In general, both methods seem to be robust against our tests (except for cluster 6 case in the MHD method), although both methods appear more sensitive to the smoothing size rather than to decentering or the number of radial bins. For all three robustness tests between the galaxy clusters, the relaxed clusters (cluster 1, 2, and 3) had an overall higher WOC than the dynamically younger clusters (cluster 4, 5, and 6). However, the MHD results did not show significant trend with the dynamical stage of sample clusters. Only cluster 6 showed an undesirable fluctuation (trend) as the smoothing factor (binning) changed. Modifying the MHD method, taking account of all contours (substructure), alleviates the difference between cluster 6 and the other clusters, although this ends up producing less robust results. However, the WOC result for cluster 6 did not exhibit these variations, and appeared to be consistent with the other clusters.

RESULTS AND DISCUSSION
Luminous tracer for dark matter: We calculated the WOC of the dark matter and various components of each galaxy cluster to investigate which component had a spatial distribution most similar to the dark matter in the clusters. We prepared maps containing (1) all stellar particles, which includes BCG, other member galaxies, and ICL, (2) cluster galaxies, which includes BCG and other member galaxies, (3) BCG and ICL. Among the various components (stellar particles, galaxies, BCG+ICL), the BCG+ICL component traced the dark matter best, except for the case of cluster 6 (see the upper panel in Figure 8). The three data points with the same symbols come from the three projection angles of each galaxy cluster, which are connected by solid bars. A possible explanation for the fact that the BCG+ICL component (without other member galaxies) traces the dark matter better rather than all stellar particles, could be that the shape of gravitational potential generated by dark matter is smoother than the sharp local minimum of potential produced by the stellar particles from other member galaxies.
Evolutionary stage estimator for galaxy clusters: The more relaxed galaxy clusters, which have had enough time to virialize, could have more aligned distributions between the dark matter and the other components of the cluster. If the BCG+ICL component of more relaxed galaxy clusters has a distribution more similar to the dark matter in them, compared to dynamically younger galaxy clusters, the similarity in spatial distribution between dark matter and BCG+ICL could be used as a dynamical stage indicator. Accordingly, we investigated the relation between the WOC(DM, BCG+ICL) result and one of the relaxation level proxy, formation time (z m/2 ). Between the sample galaxy clusters, the relaxed clusters showed a stronger similarity in spatial distribution between the dark matter and BCG+ICL than the dynamically young clusters (see the middle panel in Figure 8). The WOC -z m/2 relation shows a considerable trend, although cluster 4 seems to be off from the trend. Among the clusters with high WOC (clusters 1, 2, 3 and 4), cluster 4 exhibits a slightly larger variation with viewing angle.
Care should be taken when comparing between the results of the MHD method and the WOC method, since their ranges and scaling are different. In Figure 5, taking the linear scale of ζ M HD seems like the most appropriate choice when comparing with the WOC results. Nevertheless, we took the logarithmic scale of ζ M HD to allow for an easier qualitative visual comparison. The MHD result shows similar trend with formation time but with a larger variation inferred from the viewing angle (see the solid line in lower panel in Figure 8). Specially, cluster 6 showed significantly different ζ M HD (DM,BCG+ICL) results for different projection angles. Modifying the MHD method taking account of all contours in the calculation (dashed faint line) reduces the sensitivity to the viewing angle, however the resulting trend with dynamical stage is broken. The modified MHD improves in the merging cluster case by taking account of sub-structures, while in the relaxed cluster case, the resulting ζ M HD becomes larger compared to the one-to-one contour comparison.
The WOC method is a non-parametric comparison of two distributions, which is easy to use and does not require any fitting that may cause numerical error. It quantifies the similarity of two spatial distributions robustly, under choices of a smoothing factor, center, and binning ( Figure 7). Moreover, the WOC method works well for disconnected contour cases when there are multiple significant sub-structures, such as in the Bullet cluster or Coma cluster. Also it does not require the computation of individual contours, and thus it suffers less bias even when working with masked maps. Furthermore, WOC(DM, BCG+ICL) could be used as the dynamical stage estimator of galaxy cluster, using the WOC(DM, BCG+ICL)-z m/2 relation. Note that, the MHD method could consider the sub-structure in a way similar to the WOC method, if we did not choose the largest contour for each level but rather used all contours for the level (see the dashed faint line in the right panel in Figure 7 and lower panel in Figure 8), however the MHD method does not show a significant trend in the ζ M HD (DM,BCG+ICL)-z m/2 relation.
Application on observational data: In the WOC calculation, we require precise density maps to decide the levels and give weights using the density level at the contour. However, a computation of a density requires a robust measurement of the zero level, which is typically challenging for observations of the ICL. Therefore, when applying the WOC method to real data, we should compute a conservative detection limit of the image, and determine the binning range adequately. If the outermost bin (the largest contour) is securely within the area defined by the detection limit, the dimmer pixels outside the largest contour do not affect the WOC result. Even within the secure region, the ICL density level at each contour has an error from the background fluctuations. In the WOC weight system, all density level values are normalized (see the equation 2), thus the absolute accuracy of each density level is not so important, rather the relative density levels within the secure area are considered. The WOC M (introduced in Section 4.1) could be more problematic for observational data because the measurement of the total mass of dark matter and total luminosity of BGC+ICL requires the pixel values outside the largest contour. In either case, we should determine the detection limit boundary conservatively, to minimize this problem.
Possible applications: Aside from 2-dimensional spatial distributions, the WOC method could also be applicable to any two-parameter spaces, such as color-magnitude diagrams, phase-space diagrams, etc. One could also use our methodology to compare the characteristics of species in different galaxy clusters. Moreover, the WOC method could be simply extended to 3-dimensional spatial distributions by measuring 'overlapping volume'.  The three data points with the same symbols come from the three projection angles of each galaxy cluster, which are connected with lines. The MHD method results using only the largest contour at each level are plotted as solid lines, while the results using all contours are plotted as dashed faint lines. We shift the MHD results with all contours +0.03 in the x-axis for readability. See more details in Section 5.
In future work we will apply the method to galaxy clusters simulated in hydrodynamical simulation, such as Illustris TNG (Nelson et al. 2019) and Horizon Run5 (Lee et al. 2021), and study the relation between WOC(DM, BCG+ICL) and the cluster's dynamical stage (represented by virial ratio, sub-structure ratio, center of mass-BCG offset, formation time, etc.).
Creating mock observation maps that represent different ICL formation origins (such as BCG major merger, tidal stripping or in-situ formation) and comparing them with observed ICL maps could constrain the origin of ICL. Dynamically young clusters, whose ICL and dark matter are not well aligned, may possess a blue colored ICL, indicating freshly stripped stars. Studying the ICL color and the WOC(DM, BCG+ICL) together can provide additional information, allowing us to test such hypotheses.
Furthermore, we could utilize the method to constrain dark matter models such as the self-interacting dark matter model (SIDM) or the cold dark matter model (CDM), which predict different tidal interaction histories (Sirks et al. 2021). Because our methodology is trivially extendable to N -dimensions, we will search for use cases outside of astronomy that also require a spatial comparison of multiple distributions.

SUMMARY
We developed a new methodology to quantify the similarity of two or more 2-dimensional spatial distributions, which we call the Weighted Overlap Coefficient (WOC) method. In this methodology, we calculate the fraction of overlapping area between two distributions at various density threshold levels, while weighting using the signal strength, and normalizing the similarity measure between 0 and 1. We compared this with one of the previous methods used in the literature, the Modified Hausdorff Distance (MHD) method. Our methodology is robust when choosing smoothing factor, center, and binning (Figure 7), and performs well even with the existence of multiple sub-structures, even those far away (∼ 1 Mpc) from the center (Figure 4). Furthermore, the fixed range of the WOC result allows us to interpret the result easily and to compare between different cluster samples in a fair way.
Applying the WOC method to six galaxy clusters simulated in GRT simulations, we found that the spatial distribution of the BCG+ICL component was a better match with dark matter than other cluster components, such as all stellar particles or all member galaxies (upper panel in Figure 8). Most interestingly, galaxy clusters that were more relaxed (with earlier formation time, z m/2 ) showed higher WOC(DM, BCG+ICL) results, which would enable us to use it as cluster's dynamical stage indicator (middle panel in Figure 8). Conversely, we found that the MHD results show weaker trend with the dynamical state of the cluster (lower panel in Figure 8).
Our methodology may have broad application across astronomy and the physical sciences. The software is available on GitHub 1 under a GNU General Public License and version 0.0.1 used in this analysis is archived in Zenodo (Sabiu & Yoo 2022).