Weak-lensing Mass Bias in Merging Galaxy Clusters

Wonki Lee; Sangjun Cha; M. James Jee; Daisuke Nagai; Lindsay King; John ZuHone; Urmila Chadayammuri; Sharon Felix; Kyle Finner

doi:10.3847/1538-4357/acb76b

1. Introduction

In the hierarchical structure formation scenario, a merger is one of the major channels through which galaxy clusters grow their mass. During the collision between clusters, an order of ∼10⁶⁴ erg (e.g., Ricker & Sarazin 2001) gravitational energy is transferred to the cluster environment, which creates a unique experimental setup for testing various astrophysical processes. This includes energy transport via thermal conduction (e.g., Markevitch et al. 1999, 2003), development of the intracluster magnetic field (e.g., Markevitch et al. 2005; ZuHone et al. 2011), the electron-ion equilibrium timescale (e.g., Markevitch 2006; Russell et al. 2012), acceleration of high-beta plasma (e.g., Brunetti & Lazarian 2007; Kang et al. 2012), and self-interaction cross section of dark matter (e.g., Markevitch et al. 2004; Randall et al. 2008).

However, astrophysical constraints based on merging galaxy clusters generally suffer from uncertainties in the merger configuration (such as viewing angle, collision speed, time since the collision, pericenter distance, mass, and mass ratio). Thus, it is crucial to robustly characterize the merging scenario to utilize the cluster merger as an astrophysical laboratory.

Mass is one of the most critical parameters that affect the merging scenario. However, determining the mass of a dynamically active cluster is very challenging. Mass estimators that use the properties of thermalized plasma from X-ray (e.g., Vikhlinin et al. 2006) or Sunyaev-Zeldovich effect (Sunyaev & Zeldovich 1970, 1972; Carlstrom et al. 2002; Mroczkowski et al. 2019, for a recent review) observations are biased since the condition of hydrostatic equilibrium (HSE) does not hold during the merger (e.g., Nagai et al. 2007; Krause et al. 2012; Nelson et al. 2012, 2014; Shi et al. 2016). In this sense, weak lensing (WL) is a promising tool because its signal is only sensitive to the projected potential, not depending on the dynamical properties of the system. This unique capability has been used to map complex substructures in merging clusters (e.g., Hoekstra et al. 2000; Clowe et al. 2006; Okabe & Umetsu 2008; Jee et al. 2016; Kim et al. 2019; Yoon et al. 2020), which aid numerical simulations by providing constraints on the amount and location of the subcluster halos (e.g., Milosavljević et al. 2007; Springel & Farrar 2007; Lage & Farrar 2014; Molnar & Broadhurst 2015; Wittman et al. 2018; Zhang et al. 2018; Lee et al. 2020; Chadayammuri et al. 2022).

The weakness of WL, however, is the so-called mass-sheet degeneracy, where the signal is invariant under the transformation of the surface mass density: κ → λ κ + 1 − λ (κ: dimensionless mass density, λ: free parameter, which are described in more detail in Section 2.3). Therefore, WL practitioners generally estimate cluster masses by assuming that the halos follow analytic profiles. The most popular choice among many analytic models is the Navarro–Frenk–White (NFW; Navarro et al. 1996) model with the assumption that the mass correlates with the halo concentration (e.g., Duffy et al. 2008). Since both the NFW profile and the mass–concentration (M–c) relation are valid only for the description of the average halo properties (e.g., Jing 2000), systematic errors are inevitable when the method is applied to individual clusters. Previous studies have discussed both scatter and bias from these assumptions, showing that the triaxial shape of the dark halo (e.g., Clowe et al. 2004), projected structures along the line of sight (e.g., Hoekstra et al. 2011), miscentering (e.g., Johnston et al. 2007), and halo concentration (e.g., Schrabback et al. 2018) can cause uncertainties at the 10%–50% level.

The departure from the NFW profile and M–c relation is most pronounced in merging clusters, where two systems severely disrupt each other. Numerical studies have suggested that the halo compactness can sharply increase after the first closest passage (Roediger & ZuHone 2012; Richardson & Corasaniti 2022; Chadayammuri et al. 2022). Observational studies showed that merging clusters can have either higher (e.g., Broadhurst et al. 2008) or lower (e.g., Okabe et al. 2019) concentration values than the M–c relation prediction. Chadayammuri et al. (2022) reported that in the case of A2146, the WL mass could be greater than the pre-merger mass by ∼2 when the merger-driven concentration change is not accounted for in the model fitting. However, no detailed studies have been conducted to quantify the merging cluster WL mass bias for a wide range of masses and mass ratios at various merger phases.

In this study, we use a suite of idealized simulations of cluster mergers to understand the bias of WL mass estimation in post-mergers when the cluster masses are derived by the conventional method based on the NFW fitting with the M–c relation. To mimic an observational procedure, we generate mock WL catalogs by matching the source galaxy properties, such as intrinsic shape dispersion, measurement noise, and source densities, to those from real telescope data.

This paper is organized as follows. The setups for the idealized cluster merger simulations and the mock WL analysis are described in Section 2. We present the time evolution of concentration in Section 3, the bias of WL mass estimates in Section 4, and the application to previous observation studies in Section 5. The discussion is presented in Section 6 before the conclusions in Section 7. We adopt a Λ cold dark matter cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7. Throughout the paper, we assume that the clusters are at z = 0.2. R_Δ describes the radius where the average density becomes Δ times the critical density of the universe. M_Δ is the total mass within R_Δ.

2. Simulations

We use GAMER-2 (Schive et al. 2018), a GPU-accelerated Adaptive MEsh Refinement code, to perform idealized hydrodynamical simulations of merging clusters. We solve Euler equations using the MUSCL-Hancock method (van Leer 1984), the HLLC scheme (Toro et al. 1994), and implement a van Leer-type limiter in the slope computation. We do not implement radiative physics, including cooling and baryonic feedback, as we mainly focus on the dark halo properties on the scale greater than 0.15 r₅₀₀ (e.g., Kravtsov et al. 2006). The simulation refines each subregion of interest by defining a "patch" (also known as a "grid" or "block"). Each patch consists of a fixed number of cells (8³), and the refinement is triggered when the number of dark matter particles exceeds a threshold. This allows us to adaptively refine our merging cluster inside a 15 Mpc-side box with the sparsest resolution of ∼200 kpc in the outskirts to the finest level ∼15 kpc at the cluster center. We performed a resolution test by running simulations with the finest cell size of ∼7 kpc and found that the results, including the halo concentration and the WL mass estimate, are consistent with the ∼15 kpc resolution results.

2.1. Simulation Setups

For the description of the impact of various collision parameters on the halo concentration and the WL mass estimation, we define a reference run. The reference run is a collision between 5 × 10¹⁴ and 1 × 10¹⁴ M_⊙ clusters. We set the baryon fraction to 0.13 at R₂₀₀ based on the cosmic baryon fraction of Hinshaw et al. (2013) and the stellar mass fraction of Giodini et al. (2009). The initial distribution of dark matter particles and intracluster medium (ICM) profiles is generated using the cluster_generator ⁸ package. The position of dark matter particles is randomly generated based on the input mass profile. Then, we calculate the velocity of each particle at its position using the Eddington Formula (Eddington 1916). Similarly, for the ICM, we derive a pressure profile that satisfies HSE. The gas properties of each cell are interpolated based on the derived pressure profile and the assigned gas density profile.

We let the dark halo follow an NFW profile, which we also assumed in our mock WL analysis. The density is defined as

$\begin{eqnarray}&&{\rho }_{\mathrm{DM}}(r)=\displaystyle \frac{200{c}^{3}g(c){\rho }_{c}}{3(r/{r}_{{\rm{s}}}){\left(1+r/{r}_{{\rm{s}}}\right)}^{2}},\end{eqnarray} \tag{ 1 }$

where r_s is the scale radius, c is the concentration (c = r₂₀₀/r_s), and $g(c)={[\mathrm{ln}(1+c)-c/(1+c)]}^{-1}$ (Navarro et al. 1996). We use the M–c relation from Duffy et al. (2008) and adopt the coefficients derived at z = 0. We vary the initial concentrations within c = [2.5, 5.5] so that our study explores potential systematics due to the different coefficients in the M–c relations (e.g., Diemer & Kravtsov 2015).

We model the ICM density with the modified beta profile suggested by Vikhlinin (2006). These profiles are defined as

$\begin{eqnarray}&&{n}_{{\rm{p}}}{n}_{{\rm{e}}}={n}_{0}^{2}\displaystyle \frac{{(r/{r}_{c})}^{-\alpha }}{{(1+{r}^{2}/{r}_{{\rm{c}}}^{2})}^{3\beta -\alpha /2}}\displaystyle \frac{1}{{(1+{r}^{\gamma }/{r}_{{\rm{s}},{\rm{g}}}^{\gamma })}^{\epsilon /\gamma }},\end{eqnarray} \tag{ 2 }$

where n_p and n_e are the number densities of the protons and electrons, respectively, r_c is the core radius, r_s,g is the scale radius of the gaseous halo, and α, β, γ, describe the slope of the gas density. We set r_c (r_s,g) as 0.2 × r₂₅₀₀ (0.67 × r₂₀₀). The slope parameters are set to (α, β, γ, ) = (1.0, 0.67, 3.0, 3.0), which generates a cool-core cluster.

We place the two clusters at an initial separation of 3 Mpc, which is larger than the sum of the two R₂₀₀ values. The cluster profiles are defined up to a radius of 5 Mpc. For the cells located within 5 Mpc from both cluster centers, we add the gas densities, momenta, and pressures to populate the cell properties. The background gas density and temperature are set to ∼10⁻³⁰ g cm⁻³ and ∼10⁷ K, respectively. The initial infall velocity is set as 1200 km s⁻¹, which is a factor of ∼1.3 greater than the circular velocity at the initial separation. The reference run is set as a head-on cluster merger (i.e., zero impact parameter) in the plane of the sky. Hereafter, we refer to the more (less) massive cluster as the main (sub)cluster.

2.1.1. Variation of Collision Parameters and Cluster Masses

We test the dependence of the three collision parameters: impact parameter, velocity, and concentration by varying one parameter at a time from the reference run. The impact parameter ranges from 0 (i.e., head-on collision) to 1 Mpc, which results in the pericenter distance (the separation at the first closest passage) ranging from 0 to ∼0.35 Mpc. We vary the initial velocity from 600–1800 km s⁻¹, which is comparable to the escape velocity. This translates to a range in the collision velocity 3300–3600 km s⁻¹. The pre-merger concentration of each cluster is varied from 2.5–5.5.

Apart from the above, in order to examine the impact of the mass and mass ratios, we vary the mass of the main cluster from 10¹⁴–10¹⁵ M_⊙. For each mass of the main cluster, the subcluster mass is either fixed to 10¹⁴ M_⊙ or modified in such a way that the mass ratio is maintained to be 5:1. To create similar collision parameter setups to the reference run, we modify the impact parameter and initial velocity as follows. The initial velocity is adjusted to have the same virial ratio as the reference run. Also, the impact parameter b is modified so that b/R₂₀₀ remains the same. As in the reference run, the initial concentration is determined from the M–c relation of Duffy et al. (2008).

2.2. Definition of Dark Halo Properties

Since the main goal of the current study is to track the time evolution of halo properties and WL mass bias, it is necessary to robustly define the halo center and profile. Here we follow the Power et al. (2003) approach. First, we divide the dark matter particles into two groups based on their membership at the initial condition setup. Second, we compute the center of mass of each particle group to make an initial guess of the halo center. Third, we recursively refine the center for each group using the dark matter particles residing within a shrinking sphere. At each step, we choose the center of mass defined in the previous step and shrink the radius by 20%. The initial radius of the shrinking circle is 2 Mpc, and the search stops when the value becomes less than 0.2 Mpc. We visually verify that the center determined in this way robustly follows the density peak of the dark matter distribution.

Once the halo center is determined, we obtain R₂₀₀ and M₂₀₀ from the total (DM+ICM) density profile. As before, we compute the density using the dark matter particles and the ICM belonging to each group. For ICM, we define a scalar field, describing each cluster's gas fraction. Then, we generate ICM density profiles only using the gas mass that belongs to the cluster at the initial condition setup. However, when it comes to the measurement of the dark halo scale radius r_s, we do not fit the model (NFW+modified-beta profile) on the total density profile as the ICM drags out from the dark halo center during the major merger. Instead, r_s is calculated by fitting an NFW model (Equation (1)) on the dark matter density profile. The dark matter density profile is log-binned from r = 10 kpc to 3 Mpc, and Poisson noise is assumed on each bin. The concentration parameter is estimated from the ratio R₂₀₀/r_s.

2.2.1. Evolution of M₂₀₀ during Post-merger

Figure 1 shows the time evolution of M₂₀₀ of the main and subclusters in the reference run. The M₂₀₀(t) value at each merger stage t is normalized by the initial (pre-merger) mass M_200,IC. The epoch is described with the time since collision (TSC; the time since the two cores overlap). This first encounter (TSC = 0) happens ∼1.56 Gyr after the initial setup where they are separated by 3 Mpc. Both cluster masses stay close to their initial values at TSC < 0 and gradually increase at TSC > 0. The main cluster mass reaches its maximum at the first apocenter (∼1.2 Gyr) and slowly decreases afterward. On the other hand, the maximum of the subcluster mass occurs sooner (∼0.8 Gyr), and the decrease is faster than its increase.

The difference in the pattern of the mass change between the two clusters is attributed to the contraction of the halos due to the increased gravitational field from the pericenter passage (e.g., Roediger & ZuHone 2012). One can conjecture that this contraction will also affect the time evolution of the halo concentration in a similar way. A quantitative comparison with the halo contraction will be discussed in Section 3. The subcluster mass decrease is much more significant because its shallower potential cannot retain the outflowing mass after the core passage.

2.3. Mock WL Analysis

For our cluster merger simulations (Section 2), we carry out a mock WL analysis to investigate how the M₂₀₀ variation translates to WL mass bias. At each merger phase, we generate a synthetic WL ellipticity catalog based on the projected mass distribution and determine the halo masses by fitting two NFW profiles simultaneously.

We represent a 5 Mpc × 5 Mpc projected mass distribution with a 1000 × 1000 grid. The projected mass density Σ is converted to a convergence κ, which is the dimensionless surface mass density:

$\begin{eqnarray}&&\kappa =\displaystyle \frac{{\rm{\Sigma }}}{{{\rm{\Sigma }}}_{c}},\end{eqnarray} \tag{ 3 }$

where Σ_c is the critical surface mass density. The critical surface mass density is defined as:

$\begin{eqnarray}&&{{\rm{\Sigma }}}_{c}=\displaystyle \frac{{c}^{2}{D}_{s}}{4\pi {{GD}}_{l}{D}_{{ls}}},\end{eqnarray} \tag{ 4 }$

where c is the speed of light, G is the gravitational constant, D_l(D_s) is the angular diameter distance to the cluster (source), and D_ls is the angular diameter distance from the cluster to the source.

We generate mock WL observational data from the converted convergence above (Bartelmann & Schneider 2001). First, from the convergence map, we compute a reduced shear g field by the following:

$\begin{eqnarray}&&g=\displaystyle \frac{\gamma }{1-\kappa }.\end{eqnarray} \tag{ 5 }$

The shear γ is calculated by:

$\begin{eqnarray}&&{\boldsymbol{\gamma }}({\boldsymbol{x}})=\displaystyle \frac{1}{\pi }\int {\boldsymbol{D}}({\boldsymbol{x}}-{{\boldsymbol{x}}}^{{\prime} })\kappa ({{\boldsymbol{x}}}^{{\prime} })d{{\boldsymbol{x}}}^{{\prime} },\end{eqnarray} \tag{ 6 }$

where D is the convolution kernel:

$\begin{eqnarray}&&{\boldsymbol{D}}=-\displaystyle \frac{1}{{({x}_{1}-{{ix}}_{2})}^{2}}.\end{eqnarray} \tag{ 7 }$

Note that the computation of γ through Equation (6) is in principle biased because the convergence κ outside the field (5 Mpc × 5 Mpc) is ignored (in a typical WL study, the extent of the NFW profile is assumed to be infinite). The bias due to this finite-field effect is noticeable only at r ≳ 2.2 Mpc. In this study, we choose to exclude the WL data in this regime. In addition to this finite-field effect, the finite resolution of the grid creates some (∼10%) bias in reduced shear prediction near the cluster peak (r ≲ 250 kpc). Hence, we also avoid using the WL data at such small radii.

The reduced shear g field provides the observed ellipticity variation across the field if the intrinsic shape of the source is a perfect circle. To mimic the real observational WL data, we need to randomize the intrinsic source ellipticity, measurement noise, and galaxy position. Under the reduced shear g , the intrinsic source ellipticity e at the random position (x,y) is modified to the observed ellipticity according to the following rules:

$\begin{eqnarray}&&{\boldsymbol{\epsilon }}=\displaystyle \frac{{\boldsymbol{e}}+{\boldsymbol{g}}}{1+{{\boldsymbol{g}}}^{* }{\boldsymbol{e}}}\,(| {\boldsymbol{g}}| \lt 1)\end{eqnarray} \tag{ 8 }$

and

$\begin{eqnarray}&&{\boldsymbol{\epsilon }}=\displaystyle \frac{1+{\boldsymbol{g}}{{\boldsymbol{e}}}^{* }}{{{\boldsymbol{e}}}^{* }+{{\boldsymbol{g}}}^{* }}\,(| {\boldsymbol{g}}| \gt 1),\end{eqnarray} \tag{ 9 }$

where the asterisk (*) indicates the complex conjugate. When randomizing e , we assume σ_e = 0.25 per component. The final source galaxy ellipticity ${{\boldsymbol{\epsilon }}}^{{\prime} }$ is obtained after we add measurement noise σ_m to each component of . The measurement noise σ_m depends on the signal-to-noise ratio (S/N) of a source. We utilize the HST WL catalog of Kim et al. (2021) to derive the relation between source S/N and ellipticity measurement error. We use the following χ²-minimization to find the best-fit NFW models from the mock WL catalog:

$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{I}\displaystyle \sum _{j=1}^{2}\displaystyle \frac{{({g}_{j}^{i}-{\epsilon }_{j}^{{\prime} i})}^{2}}{{{\sigma }_{s,i}}^{2}+{{\sigma }_{m,i}}^{2}},\end{eqnarray} \tag{ 10 }$

where I is the total number of background source galaxies, ${\epsilon }_{j}^{{\prime} i}$ is the jth component of the noisy ellipticity for the ith source, and g_j ⁱ is the jth component reduced shear expected at the location of the ith source. We set the source redshift to z_s = 1 and the source density to 100 (25) galaxies $\,{\mathrm{arcmin}}^{-2}$ to follow the HST (Subaru) WL catalog. At the redshift of the cluster z = 0.2, the field size corresponds to 25 farcm 86 × 25 farcm 86. If not specified, the WL mass estimate assumes the Duffy et al. (2008) M–c relation in the NFW model fitting.

3. Evolution of Concentration

Figure 2 presents the time evolution of the concentration of the main cluster (red) and the subcluster (blue) in the reference run. As with the M₂₀₀ investigation, we normalize the concentration with its initial (pre-merger) value. The halo concentration significantly evolves with time since the collision. Before the first encounter, similar to the M₂₀₀ evolution, the concentration roughly stays consistent with its initial value. Shortly after the collision, we observe a rapid increase in concentration for both clusters, which is followed by a gradual decrease. This behavior is consistent with the findings of Chadayammuri et al. (2022). At the maximum, the concentration values of the main cluster and the subcluster are factors of ∼1.5 and ∼2.0 greater than their initial values, respectively. The concentration changes faster in the subcluster; for the main (sub)cluster, it takes ∼0.23 Gyr (∼0.13 Gyr) since the first encounter to reach the maximum. Hereafter, we refer to the TSC for the concentration to reach its maximum as t_peak.

In the later phase of the merger, the concentration decreases below its initial value, which indicates that the halo has become more diffuse than in the pre-merger state. The concentration of the main cluster reaches its minimum value (TSC ∼ 1.7 Gyr, ∼60% of its initial value) during the second infall (see the subplot of Figure 2). The main cluster concentration recovers its initial value at around the second closest passage (∼2.5 Gyr), which is followed by small oscillations due to the repetitive passages of the subcluster.

This evolutionary trend of the concentration can be explained by a halo contraction. As mentioned above, the halo contracts because of the increased gravitational field during the first core passage (e.g., Roediger & ZuHone 2012; Chadayammuri et al. 2022), which increases the concentration. Physically, a larger concentration change is expected in the subcluster, as it experiences a more significant relative increase in gravity.

The halo contraction can also explain the difference in t_peak. When we regard the halo contraction as the collective behavior of the modified orbits of dark matter particles, the concentration should peak at half the dynamical time t_dyn, which dark matter particles need to reach the center. We can estimate t_dyn using the mass within the scale radius at the closest passage. Then, t_dyn becomes ∼0.6(0.3) Gyr for the main (sub)cluster. This gives a t_dyn/t_peak ratio of ∼0.4 in both clusters, which implies that the difference in t_peak can be attributed to the different orbital time of the two halos. The actual ratio is smaller than 0.5 because the matter beyond the scale radius contracts and contributes to the enhancement of the gravitational field in practice. This result is consistent with the analysis of Richardson & Corasaniti (2022), which shows that a "universal evolution" is obtained when the time is normalized with the dynamical time.⁹

3.1. Dependence on Collision Parameters

We show the dependence of the concentration variation on collision parameters in Figure 3. In both clusters, a larger impact parameter results in a smaller peak value, which happens at a later epoch (top left). The dependence on the initial velocity is insignificant (top middle) within the range explored here. A cluster with a higher pre-merger concentration leads to faster evolution, a lower amplitude of the concentration variation, and higher amplitude of the concentration variation in the companion cluster (bottom left and bottom right).

The dependence on the collision parameters can be explained with the same halo-contraction-based reasoning mentioned in the discussion of the reference run. Faster evolution with a larger concentration change is expected whenever the revised collision parameter provides a condition for an enhanced gravitational interaction at the closest passage. For example, a larger impact parameter lessens the interaction and thus leads to a smaller and slower concentration change. A higher initial concentration, which gives a larger mass within the fixed radius, imparts a stronger gravitational force to the companion cluster.

As mentioned earlier, variation in the initial velocity shows a negligible impact on the time evolution of concentration. This is because the collision velocities are insensitive (≲10%) to the initial velocity change explored in the current study (600–1800 km s⁻¹); even if we assume a freefall (v_init = 0), the resulting collision velocity is close to that in the v_init = 600 km s⁻¹ case. In summary, at the given mass, the initial velocity has an insignificant impact on the collision velocity and the time evolution of concentration.

The right panels of Figure 3 shows the concentration evolution for different masses and mass ratios. When the subcluster mass is fixed (top-right panel), the amplitude of the main (sub)cluster concentration variation decreases (increases) for increasing the main cluster mass. Again, this trend can be explained by the halo contraction argument because the gravitational impact on the subcluster increases with the main cluster mass while the impact on the main cluster itself decreases. As a result, a larger main cluster mass increases (decreases) the amplitude of the (main) subcluster concentration variation.

The bottom-right panel of Figure 3 shows that the concentration evolution is independent of the cluster masses as long as the mass ratio is fixed. This is not surprising because the relative gravitational impacts also remain unchanged in this case; we verified this argument by repeating the experiment with a 2:1 mass ratio.

4. Impact on Weak-lensing Mass

4.1. Bias of WL Mass Estimation

We present WL mass estimates of simulated cluster mergers at different phases in Figure 4. We discuss two types of mass biases. The left panel presents the bias with respect to the initial (pre-merger) mass.¹⁰ The middle panel presents the bias with respect to the mass at the current (observed) epoch. Since the intrinsic mass defined by the density contrast, such as M₂₀₀, indeed changes during the merger, this bias is critical when we are mainly interested in the cluster properties at the current phase. Hereafter, we compare these mass biases with the statistical uncertainty derived assuming that the source galaxy properties follow the HST observations. Assuming Subaru observations instead gives a consistent mass estimate, albeit with a factor of ∼2 greater mass uncertainty due to its lower source density (see dark and light regions in Figure 4).

The evolutionary trend of the WL mass estimate scaled by the initial mass is similar to that of concentration. Before the collision, the WL mass is a fair representation of the initial mass for both clusters. Then, the WL mass rapidly increases after the first passage, reaches the maximum, and decreases. At its maximum, the WL mass is ∼30% (∼60%) higher than the initial mass of the main (sub)cluster. Similar to the concentration result, the subcluster WL mass peaks earlier (∼0.23 Gyr) than the main cluster (∼0.34 Gyr). This large difference between the initial and WL mass at the observed epoch is consistent with the findings of Chadayammuri et al. (2022). Interestingly, this positive correlation in the evolutionary pattern between the WL mass and halo concentration contradicts the negative correlation implied by the M–c relation.

The middle panel of Figure 4 shows that even when compared with the true mass at the same epoch (M₂₀₀(t)), the WL mass is still biased, and its overall variation resembles the pattern in the left panel in the sense that the bias quickly rises after the first passage, reaches the maximum, and decreases. The largest M_WL/M₂₀₀(t) ratio for the main cluster mass is ∼20%, similar to the M_WL/M_200,IC ratio. We cannot reliably estimate the equivalent ratio for the subcluster because the maximum bias seems to happen when the two halos are too close, and thus the simultaneous two-halo fitting becomes unstable. In late phases, the WL mass bias shifts from overestimation to underestimation. The underestimation appears earlier with a more significant drop for the subcluster, reaching ≲60% at TSC ∼ 1.0 Gyr.

Interestingly, although the WL mass of each cluster is significantly biased during the merger, the mass ratio between the two clusters based on the WL results is relatively stable at the early (TSC ≲ 0.6 Gyr) phases (right panel). The mass ratio derived from WL increases at late (TSC ≳ 0.6 Gyr) phases, as the WL mass estimate of the subcluster is significantly underestimated.

To trace the source of the WL mass bias, we compare the projected mass distributions from the best-fit WL model and the values taken directly from the simulation. The top panel of Figure 5 shows the result for the reference run. The best-fit WL model provides a good description of the true mass distribution during the pre-merger phase (TSC = −0.6 Gyr). When the TSC is 0.25 Gyr, the WL result nicely recovers the central density of both halos, whereas it overestimates the density at large radii. As mentioned earlier, the M_WL/M₂₀₀(t) ratio at this epoch is ∼1.2 (∼1.1) for the main (sub)cluster. This mass overestimation is attributed to the departure from the M–c relation due to the merger-induced halo contraction; the main (sub)cluster concentration increases by ∼40% (∼70%), as shown in Figure 2. The increased concentration enhances the WL signal near the cluster center, which dominates the overall WL S/N. Therefore, when we fit an NFW profile to the signal with an M–c relation, a combination of low-concentration and high-mass results.

The mass bias at TSC = 0.75 Gyr can be explained by a similar logic. In Figure 4, we saw that at this epoch, the WL analysis recovers the main cluster mass within 10%, but significantly (≲60%) underestimates the subcluster mass. The WL mass estimate of the main cluster is close to the truth because the concentration at TSC = 0.75 Gyr is also close to the initial value (Figure 2). On the other hand, the subcluster concentration drops by ∼70%, which is a significant departure from the M–c relation. This suggests that the halo concentration is a critical parameter for accurate mass estimation during the post-merger. One may argue that an NFW model might not adequately describe the mass distribution at a late merger phase because of the asymmetric mass distribution. Nevertheless, the main cluster's M₂₀₀ is reproduced from the mock WL analysis when the halo concentration returns to the initial value (i.e., when the M–c relation holds).

The model-to-true mass ratio in Figure 5 is large near the edges. This is because the mass distribution has contracted after the collision, and the density profile becomes steeper than the NFW model prediction (∝r⁻³). Although these regions present a large ratio, they do not contribute significantly to the total mass. The WL signal from this region is small.

The bottom panel of Figure 5 displays the ratio of the WL model mass to the DM-only mass. This ratio of the WL model mass to the DM-only mass helps us to identify the bias due to the (gas) model mismatch as our WL analysis fits an NFW model to the lensing signal contributed by the total (DM+ICM) mass profile, while the ICM profile does not follow the NFW profile. Comparison between the top and bottom panels illustrates that the difference due to the model mismatch at the cluster core is substantially reduced when the WL mass is compared with the DM-only mass.

4.2. WL Mass Bias Reduction

In Section 4.1, we find that the WL mass bias is reduced at the epoch when the halo concentration is consistent with the M–c prediction. We can, therefore, mitigate the WL mass bias by either not relying on the M–c relation (e.g., High et al. 2012; Finner et al. 2021; Kim et al. 2021) or by obtaining a good prior knowledge of the true concentration from an independent observation (e.g., Biviano et al. 2017).

Figure 6 presents the WL results from our two-parameter (mass and concentration) fitting (top) and one-parameter (mass) fitting with the concentration fixed at the truth (bottom); the true concentration is measured from the simulation data. When we simultaneously fit concentration and mass, the overestimation of the main cluster WL mass at TSC = 0.25 Gyr¹¹ is reduced from ∼20% to ∼4% (top left). The mass ratio map also shows that the simultaneous fitting provides a better description of the true mass distribution (compare the top-middle panel in Figure 6 with the top-middle panel in Figure 5). However, this simultaneous fitting does not always improve the accuracy of the mass estimation. At a later phase of the merger (TSC = 0.75 Gyr), the WL mass bias of the main cluster is greater than the result from the previous M–c relation case. Since the level of the maximum overestimation is similar (∼20%), we conclude that this simultaneous fitting is not a viable solution to resolve the WL mass bias in merging clusters.

On the other hand, we find that using the accurate concentration can mitigate the WL mass bias. As shown in the bottom panel of Figure 6, the WL mass estimate of the main cluster is accurate within ∼3% throughout the merger. For the subcluster, the accuracy is ≲10% at TSC ≲ 0.8 Gyr.

This result suggests that WL can be an unbiased mass estimator in merging clusters with the aid of the true concentration. The halo concentration value can independently be determined with different observational tracers, such as the distribution of cluster galaxies (e.g., Nagai & Kravtsov 2005; Biviano et al. 2017) or strong lensing (e.g., Oguri et al. 2012; Giocoli et al. 2014). In particular, we expect strong lensing to be a suitable solution in recent massive mergers, as the merger-boosted concentration will enhance the strong lensing signals (Fedeli et al. 2008).

4.3. Mass Dependence of WL Mass Bias

Figure 7 shows how the WL mass bias changes as we vary the initial masses. Since the true and WL masses change as the merger progresses, we use the maximum ratio with respect to the initial-true (top) and current-true (bottom) masses to quantify the dependence. Since our snapshot time resolution is finite, the maximum value chosen from the snapshot data is always underestimated. Thus, we refine the maximum ratio by fitting a second-order polynomial to the three values within 20 Myr from the tentative maximum. Note that, sometimes, the epoch where the maximum occurs likely resides within the time interval where our simultaneous fitting is unstable (i.e., the halo separation is too small). For these cases, we mark the value with an asterisk.

With respect to the initial-true mass, the WL mass bias of both clusters increases with the companion cluster mass. For the main cluster, the mass bias decreases with its mass, whereas this self-dependence is insignificant for the subcluster. It is remarkable that for the 1:1 merger of two 10¹⁵ M_☉ clusters, the maximum bias is nearly ∼80%.

When the bias is measured with respect to the current-true mass, the dependence on the main cluster mass is insignificant in both clusters. For instance, with a 10¹⁴ M_☉ subcluster, the 10¹⁴–10¹⁵ M_☉ variation in the main cluster mass affects the bias of the main cluster at the ∼4% level, which is ∼10% of the uncertainty. This insensitivity can be understood from the evolution of the current-true mass M₂₀₀(t). The main (sub)cluster shows an increased time variation of M₂₀₀(t) for decreasing (increasing) main cluster mass and increasing (decreasing) subcluster mass, which is similar to the mass dependence of M_WL that was discussed above. Therefore, the dependence on the main cluster mass is canceled in the evaluation of M_WL/M₂₀₀(t). We note that a similar mass dependence is observed when we simultaneously fit the mass and concentration (as in Section 4.2).

Figure 7 presents significant implications for interpreting various merger observations. For instance, WL mass bias will be small or negligible between low-mass (∼10¹⁴ M_⊙) systems or in a massive cluster with a group-size collider. Even at its maximum, the level of bias is insignificant, compared to the statistical uncertainty of the WL mass (∼50%; Finner et al. 2021). In contrast, the bias is more significant in massive mergers. The difference between the WL mass estimate and the current-true mass can reach up to ∼40% in the collision of two moderately massive (5 × 10¹⁴ M_⊙) clusters. As this bias is significant relative to the statistical uncertainty (≲30%; Kim et al. 2021), the bias correction is required.

Another important takeaway from this investigation is that the maximum bias occurs in an early phase of the merger (0.2–0.4 Gyr after the first encounter), which includes several merger cases reported in the literature. Thus, some existing WL studies may have overestimated the masses, mainly when the cluster selection is based on their prominent merger features.

5. Application to Observations

Our mock WL analysis shows that a WL mass can be significantly overestimated in a massive recent merger. The large discrepancy between the WL estimate and the true pre-merger mass significantly impacts our interpretation of the merger, including the collision scenario reconstruction. In this section, we investigate the impact for the three observational cases: A2034, MACS1752, and ZwCl1856, which have been recently identified as exemplary binary mergers with a near head-on collision.

5.1. Review of Observations and Simulation Setup

A2034 is a dissociative merger where its merger shock is detected with the X-ray feature discontinuity (Owers et al. 2014), whereas MACS1752 and ZWCL1856 are the mergers with double radio relics (Bonafede et al. 2012; de Gasperin et al. 2014). A radio relic is a megaparsec-scale extended diffuse radio emission feature, which often traces the merger shock front, as seen in X-rays (e.g., Feretti et al. 2012). Radio relics come in various shapes due to the diversity in merger configurations and in cluster environments. The double radio relic systems in MACS1752 and ZWCL1856 are remarkably symmetric and highly consistent with the hypothesized collision axis and the locations of the merger shocks, as predicted by numerical simulations (e.g., Ha et al. 2018).

The A2034 WL analysis was performed by Finner (2021). Finner et al. (2021) presented WL studies of MACS1752 and ZWCL1856. Based on the source density of ≳20 galaxies arcmin⁻², these WL studies identified two halos in each system, and their masses were derived by NFW profile fitting with the M–c relation, as is done in our mock WL analysis. Finner (2021) showed that A2034 is a ∼4:1 merger between 3.6 ± 0.7 and 0.9 ± 0.5 × 10¹⁴ M_⊙ halos. According to Finner et al. (2021), MACS1752 (ZWCL1856) is an equal-mass merger with halos of mass 4.7 ± 0.8 and 3.6 ± 0.7 × 10¹⁴ M_⊙ (1.2 ± 0.6 and 1.0 ± 0.6 × 10¹⁴ M_⊙). Readers are referred to Finner (2021) and Finner et al. (2021) for the analysis details.

For simplicity, we describe the three systems with the combinations of 5 × 10¹⁴ M_⊙ and 1 × 10¹⁴ M_⊙ halos. Specifically, we represent A2034 with the head-on collision between 5 × 10¹⁴ and 1 × 10¹⁴ M_⊙ halos (i.e., reference run) whereas the MACS1752 (ZwCl1856) analog is a merger between two 5(1) × 10¹⁴ M_⊙ halos. X-ray observations clearly show that both MACS1752 and ZwCl1856 maintain their gas cores. Thus, we set a nonzero impact parameter of 0.6(0.3) Mpc for MACS1752 (ZwCl1856) so that the gas core survives after the first passage. One can fine-tune the simulation setup so that our mock WL mass closely matches the reported values in Finner et al. (2021) and Finner (2021). However, we do not attempt to do it here since our primary goal is to estimate the level of bias to first order.

5.1.1. Shock-based Time Estimation

Before we discuss the time evolution of the WL mass bias, we introduce the shock-based time τ_sh:

$\begin{eqnarray}&&{\tau }_{\mathrm{sh}}=\displaystyle \frac{{d}_{\mathrm{sh}}}{{V}_{200}},\end{eqnarray} \tag{ 11 }$

where d_sh and V₂₀₀ are the shock-to-shock distance and the virial velocity based on the summation of the two M₂₀₀ values, respectively. We assume that the merger is happening on the plane of the sky. Thus, d_sh and τ_sh are the lower limits. Since we can identify the shock locations in both observation and simulation, τ_sh is a convenient measure of the TSC applicable to both observation and simulation.

For the observations of MACS1752 and ZwCl1856, we adopt the distance between the two radio relics as d_sh. Since the A2034 shock was detected only on one side so far, we use the distance between the shock and the main halo as the time proxy.

For the simulation analogs, we create a gas velocity divergence map and find the locations of the minimum (maximum) convergence (Ryu et al. 2003; Skillman et al. 2008). In the case of the MACS1752 and ZwCl1856 analogs, we adopt τ_sh as the distance between the two local minima. For A2034, we use the distance between the main halo-side shock and the main halo to match the observational procedure.

Figure 8 displays the relation between our shock-based time measure τ_sh and the true time (TSC) from the simulation analogs. It is important to note that the relation is quasi-linear for the time interval shown here. This is because V₂₀₀ is nearly linear with time, although the mass sum is significantly nonlinear. The τ_sh − TSC linearity demonstrates that τ_sh can serve as a valuable measure of the TSC, as its true value is not directly accessible from observations.

We can read off the TSC values for the three clusters with Figure 8. They are estimated to be ∼0.2, ∼0.5, and ∼0.7 Gyr for A2034, MACS1752, and ZWCL1856, respectively. These estimates are comparable to those in the previous studies. In the case of A2034, our estimate ∼0.2 Gyr is slightly lower (earlier merger phase) than the literature values. Owers et al. (2014), Moura et al. (2021), and Monteiro-Oliveira et al. (2018) quoted ∼0.3 Gyr, ∼0.3 Gyr, and ∼0.5 Gyr based on the X-ray shock velocity, the offset between X-ray and DM peaks, and the spectroscopic and WL analysis, respectively. The difference mainly comes from the viewing angle since our plane-of-the-sky-collision assumption provides a lower limit. We will discuss the impact of the viewing angle on the WL mass bias further in Section 6.1. Nevertheless, the key characteristics, such as the X-ray front over the subcluster in A2034 or the X-ray tails of MACS1752, are well reproduced at τ_sh (top panel of Figure 8). This reassures us that τ_sh is a reasonable proxy for the TSC for the first-order interpretation of the observations.

5.2. WL Mass Bias at Observed Epoch

Based on the WL bias versus TSC relation from the simulation analogs, we can determine the WL bias for the three systems. With respect to the current mass, the WL masses are likely to be overestimated by ∼20% and ∼35% for A2034 and MACS1752, respectively, whereas a ∼10% underestimation is expected for ZwCl1856. For A2034, the bias is close to the maximum since the estimated TSC ∼ 0.2 Gyr is near the epoch of the maximum bias (TSC ∼ 0.25 Gyr; see the middle panel of Figure 4). In the case of MACS1752, the WL bias is estimated to be greater than that of A2034 mainly because of the higher system mass. Finally, we expect a mass underestimation for ZWCL1856 due to its late merger phase. As mentioned earlier, the mass bias in MACS1752 is comparable to the statistical uncertainty.

5.3. WL Mass Bias at Pre-merger Epoch and Method for Correction

Regarding following up merger observations with simulations, priority should be given to accurate mass estimation at the pre-merger epoch. We find that the WL mass estimate at the observed epoch is larger than the initial mass at the pre-merger epoch in all three cases. We expect ∼30%, ∼70%, and ∼20% overestimations for the main cluster mass of A2034, MACS1752, and ZwCl1856, respectively.

Our scheme for the derivation of the bias-corrected pre-merger mass is as follows. We utilize the notion that the WL analysis can reproduce the initial mass ratio although each mass is biased (Section 4.1). First, we generate a series of merger simulations with the observed mass ratio while varying the total system mass. This allows us to map out the time evolution of the WL masses for different total masses. Then, we find the simulation set whose predicted WL mass matches the observation. The initial mass of the matching simulation set can be adopted as the bias-corrected pre-merger mass.

Figure 9 illustrates our scheme for the equal-mass merger samples: MACSJ1752, ZWCL1856, and ACT J0102-4915. For MACS1752, we find that the observed WL mass estimate ∼4.1 × 10¹⁴ M_⊙ at τ_sh ∼ 1.6 Gyr (TSC ∼ 0.5 Gyr) is reproduced from the collision of two 3 × 10¹⁴ M_⊙ halos. Similarly, we can explain the observed WL mass of ZWCL1856 (∼1.1 × 10¹⁴ M_⊙) with an equal-mass merger between two 10¹⁴ M_⊙ halos. The ACT J0102-4915 system, nicknamed "El Gordo," is one of the most massive mergers to date (Kim et al. 2021). Figure 9 suggests that the current observed WL masses (∼10¹⁵ and ∼6.5 × 10¹⁴ M_☉) can result from the collision of two ∼5 × 10¹⁴ M_⊙ halos when the merger is observed at TSC ∼ 0.3 Gyr. The expected pre-merger mass of El Gordo is still massive, but it is ∼40% less massive than the WL mass estimate (Kim et al. 2021). The estimated TSC of El Gordo is a factor of ∼2 larger than the result from Zhang et al. (2015), who used a set of idealized collisions between two 10¹⁵ M_⊙ halos to represent the cluster merger.¹²

6. Discussion

6.1. Dependence on Viewing Angle

Throughout the paper, we assumed a plane-of-the-sky merger (i.e., zero viewing angle). This assumption might not be entirely invalid for double radio relic systems, because a projection with a large viewing angle is known to hinder the identifications of radio relics (see Figure 7 of Skillman et al. 2013). However, it is reasonable to assume that the real merger axis is somewhat tilted with respect to the plane of the sky. This section examines the impact of nonzero viewing angles on the WL mass bias. Figure 10 shows that the viewing angle has a negligible impact on the time evolution of WL mass bias. With increasing viewing angle, the M_WL/M₂₀₀(t) ratio marginally increases, yet the difference is within the uncertainty. The WL mass bias is consistent among different viewing angles because of the model fitting. In our mock WL analysis, the mass distribution is simultaneously fitted with the model of two projected spherical halos. This allows the WL analysis to measure the overlapping mass from the two halos, so the viewing angle does not significantly change the final mass estimate.

**Figure 10.** Time evolution of the WL mass bias with a different viewing angle (i). Darker colors represent a larger viewing angle, and the shaded region marks the uncertainty in the plane-of-the-sky collision (i = 0). The WL mass estimate marginally increases with the increasing viewing angle, but the variation is within the uncertainty.
Download figure:
Standard image High-resolution image

Nevertheless, we caution that the viewing angle is still a critical parameter because the constraint on TSC using shock separation can vary with the viewing angle. Nonetheless, we can correct the WL mass bias to first order with the method described in Section 5.3 by deprojecting the observed shock separation with the viewing angle constrained from independent observations (e.g., radio polarization; Ensslin et al. 1998). Note that our discussions on the observational examples in Section 5.3 still hold because these mergers might be happening nearly in the plane of the sky (Wittman et al. 2018).

6.2. WL Mass Estimate in a Late Merger Phase

In Section 4, we showed that the WL analysis underestimates the halo mass as the clusters reach their first apocenter. At the same epoch, we also found that the subcluster quickly disrupts. Therefore, the observations will barely recognize the subcluster center as in the simulations, and so the two halo mass fittings will be unavailable in a late merger phase.

In Figure 11, we perform the WL analysis in a late merger phase of our reference run with a single halo model. The single halo model is fitted with its center fixed on the main cluster. We scale the WL mass estimate, M_WL,Single, with the current mass, M₂₀₀(t), that we measured using all dark matter particles and gas mass density.

We find that the merger can have a long-lasting effect on the WL mass estimate. In the plane-of-the-sky collision, the WL analysis moderately underestimates the current mass, yet the difference is comparable to the mass uncertainty. On the other hand, the WL mass estimate overestimates the mass by ∼20% in the line-of-sight mergers. This overestimation is more significant than the mass uncertainty until TSC ∼4.5 Gyr.

We can explain the mass overestimation in the line-of-sight mergers with the discussions in Section 4. The WL analysis over(under)estimates the current mass when the halo concentration is larger (smaller) than what the M–c relation predicts. Thus, mass overestimation is expected in the line-of-sight collision, as the mass distribution is compact in the center by the projection. This result is consistent with the WL mass bias of fitting triaxial halos, where mass is overestimated when projected along the major axis (e.g., Clowe et al. 2004).

An important takeaway is that an old merger event can affect the current WL mass estimate. Clusters can experience multiple mergers in ∼5 Gyr. Our experiment suggests that a collision within the past ∼5 Gyr can bias the WL mass estimate at the current epoch. Moreover, the WL mass bias was comparable to the mass uncertainty in a single system. However, the bias will become significant compared to the uncertainty of the mass averaged in a cosmological context. The WL mass bias throughout cosmic history is outside the scope of this work, and we plan to address it in future studies.

7. Conclusions

In this study, we have investigated the bias of WL mass estimation in deriving the pre-merger and the current mass of merging clusters. We performed a realistic mock WL analysis with observed source galaxy properties on a suite of cluster merger simulations. The main findings of our paper can be summarized as follows.

1.
We found that the concentration of the dark halo increases after the first encounter, gradually decreases, and recovers after the second closest passage. This evolutionary trend can be explained by gravitational contraction of the halo. We found that the time evolution of concentration varies with different collision parameters, and the trend was consistent with the total mass when the mass ratio is fixed.
2.
In a collision between 5 and 1 × 10¹⁴ M_⊙ clusters, the WL analysis can overestimate the current mass by up to ∼20% at TSC ∼0.2 Gyr. The current mass is underestimated in a late merger phase. Despite the bias in WL mass estimate, we found that the mass ratio between the two halos can be reproduced during the early stage of the merger (TSC ≲ 0.6 Gyr).
3.
We found a positive correlation between the time evolution of the WL mass and that of the halo concentration. Based on this correlation, we attribute the WL mass bias to the halo concentration that departs from the M–c relation during the cluster merger. We can control the WL mass bias at ∼3% level by providing the accurate concentration value using the distribution of cluster galaxies or lensing observations.
4.
We found that the WL mass bias increases with the mass of the companion cluster. The difference between the WL and current mass can reach ∼60% in a collision between 10¹⁵ M_⊙ halos. Based on this mass dependence, we suspect that the bias has a negligible impact on minor or low-mass cluster mergers. In contrast, a proper bias correction is necessary for massive cluster mergers.
5.
We estimated the WL mass bias in the observational examples of A2034, MACS1752, and ZwCl1856. We demonstrated that the spatial separation of the merger shocks could provide a reliable reference for estimating TSC. We found that the WL mass estimate can be over- and underestimated depending on the merger phase and the total mass. We further derive the bias-corrected initial mass by utilizing the fact that WL analysis can reproduce the initial mass ratio. We found that the difference between the WL mass estimate and the pre-merger mass increases with the halo mass. The difference becomes ∼40% in 10¹⁵ M_⊙ mass clusters, which can modify the merger history reconstructed based on the WL mass estimate.

Our results suggest that the previous WL mass estimates of massive merging clusters, using the parametric model fitting, may be biased high. Since the bias depends on the TSC, an accurate reconstruction of merger history is essential for alleviating the WL mass bias of merging clusters. This study was limited to merging clusters inside an isolated box. For a more direct comparison with the observations, future studies need to examine the WL mass bias in a cosmological volume.

We thank the anonymous referee for the useful comments and suggestions. M.J.J. acknowledges support for the current research from the National Research Foundation (NRF) of Korea under the program 2022R1A2C1003130. J.A.Z. acknowledges support from the Chandra X-ray Center, operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060.

Weak-lensing Mass Bias in Merging Galaxy Clusters

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction