GalWeight: A New and Effective Weighting Technique for Determining Galaxy Cluster and Group Membership

The problem of contamination of kinematic samples of galaxies in clusters by foreground and background galaxies is longstanding. It arises because of the fact that only the projected positions and velocities of galaxies are measured in redshift surveys. Due to the lack of knowledge about the motion perpendicular to the line of sight, it is difficult to judge a priori which of the galaxies found close to a cluster in projected space are actually bound to it and a good tracer of the underlying potential. Excluding fiducial members or including unbound galaxies, or interlopers, may lead to significantly incorrect estimates of the cluster mass.

Several methods have been suggested in the literature to address this problem. All these methods aim at cleaning the galaxy sample by removing nonmembers before attempting a dynamical analysis of the cluster. Some algorithms utilize only the redshift information, such as (i) the 3σ-clipping method (Yahil & Vidal 1977), which iteratively eliminates interlopers with velocities greater than 3σ; (ii) the fixed gapper technique (Beers et al. 1990; Zabludoff et al. 1990) by which any galaxy that is separated by more than a fixed value (e.g., 1σ of the sample or 500–1000 km s⁻¹) from the central body of the velocity distribution is rejected as a nonmember; or (iii) the jackknife technique (Perea et al. 1990), which removes the galaxy whose elimination causes the largest change in the virial mass estimator. These methods are primarily based on statistical rules and some selection criteria. Other algorithms utilize both position and redshift information, such as (i) the shifting gapper technique (Fadda et al. 1996) which applies the fixed gapper technique to a bin shifting along the distance from the cluster center, or (ii) the den Hartog & Katgert (1996) technique that estimates the maximum (escape) velocity as a function of distance from the cluster center calculated either by the virial or projected mass estimator (e.g., Bahcall & Tremaine 1981; Heisler et al. 1985).

In addition to the techniques described above, the spherical infall models (hereafter referred to as SIMs, e.g., Gunn & Gott 1972; Yahil 1985; Regős & Geller 1989; Praton & Schneider 1994) can determine the infall velocity as a function of distance from the cluster center. The SIM in phase-space has the shape of two trumpet horns glued face to face (Kaiser 1987), which enclose the cluster members. However, studies show that clusters are not well fit by SIMs in the projected phase-space diagram, because of the random motion of galaxies in the cluster outer region caused by the presence of substructure or ongoing mergers (van Haarlem & van de Weygaert 1993; Diaferio 1999). A recent investigation (Abdullah et al. 2013) showed that SIMs can be applied to a sliced phase-space by taking into account the distortion of phase-space due to transverse motions of galaxies with respect to the observer and/or rotational motion of galaxies in the infall region in the cluster rest frame. However, that is out of the scope of the current paper.

Another sophisticated method is the caustic technique described by Diaferio (1999) that, based on numerical simulations (Serra & Diaferio 2013), is estimated to be able to identify cluster membership with ∼95% completeness within 3r_v (r_v is the virial radius defined below). The caustic technique depends on applying the two-dimensional adaptive kernel method (hereafter, 2DAKM, e.g., Pisani 1993, 1996) to galaxies in phase-space (R_p, v_z), with the optimal smoothing length ${h}_{\mathrm{opt}}=(6.24/N)\sqrt{({\sigma }_{{R}_{p}}^{2}+{\sigma }_{{v}_{z}}^{2})/2}$, where ${\sigma }_{{R}_{p}}$ and ${\sigma }_{{v}_{z}}$ are the standard deviations of projected radius and line-of-sight velocity, respectively, and N is the number of galaxies. ${\sigma }_{{R}_{p}}$ and ${\sigma }_{{v}_{z}}$ should have the same units and therefore the coordinates (R_p, v_z) should be rescaled such that $q={\sigma }_{v}/{\sigma }_{{R}_{p}}$, where q is a constant that is usually chosen to be 25 or 35 (additional details about the application of this technique may be found in Serra et al. 2011).

One more technique that should be mentioned here is the halo-based group finder (Yang et al. 2005, 2007). Yang et al. (2007) were able to recover true members with ∼95% completeness in the case of poor groups ($\sim {10}^{13}\,{M}_{\odot }$). However, they found that the completeness dropped to ∼65% for rich massive clusters ($\sim {10}^{14.5}\,{M}_{\odot }$). Also, theirs is an iterative method that needs to be repeated many times to obtain reliable members. Moreover, its application depends on some assumptions and empirical relations to identify the group members.

This paper introduces a simple and effective new technique to constrain cluster membership that avoids some issues facing other techniques, e.g., selection criteria, statistical methods, assumption of empirical relations, or need for multiple iterations. The paper is organized as follows. The simulations used in the paper are described in Section 2. In Section 3, the GalWeight technique is introduced and its efficiency at identifying bona fide members is tested on MultiDark N-body simulations. In Section 4, we compare GalWeight with four well-known existing cluster membership techniques (shifting gapper, den Hartog, caustic, SIM). We apply GalWeight to 12 Abell clusters (including the Coma cluster) in Section 5, and present our conclusions in Section 6. Throughout this paper, we adopt ΛCDM with ${{\rm{\Omega }}}_{m}=0.3$, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.7$, and ${H}_{0}\,=100$ h km s⁻¹ Mpc⁻¹, h = 1.

In this section, we describe the simulated data that we use in this work in order to test the efficiency of the GalWeight technique to recover the true membership of a galaxy cluster.

1. MDPL2: The MDPL2⁴ simulation is an N-body simulation of 3840³ particles in a box of comoving length 1 h⁻¹ Gpc, mass resolution of $1.51\times {10}^{9}$ h⁻¹ M_⊙, and gravitational softening length of 5 h⁻¹ kpc (physical) at low redshifts from the suite of MultiDark simulations (see Table 1 in Klypin et al. 2016). It was run using the L-GADGET-2 code, a version of the publicly available cosmological code GADGET-2 (Springel 2005). It assumes a flat ΛCDM cosmology, with cosmological parameters ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.692$, ${{\rm{\Omega }}}_{m}=0.307$, Ω_b = 0.048, n = 0.96, σ₈ = 0.823, and h = 0.678 (Planck Collaboration et al. 2014). MDPL2 provides a good compromise between numerical resolution and volume (Favole et al. 2016). It also provides us with a large number of clusters of different masses extended from $0.7\times {10}^{14}$ to $37.4\times {10}^{14}$ h⁻¹ M_⊙ (the range used to test the efficiency of GalWeight).

2. Bolshoi: The Bolshoi simulation is an N-body simulation of 2048³ particles in a box of comoving length 250 h⁻¹ Mpc, mass resolution of $1.35\times {10}^{8}$ h⁻¹ M_⊙, and gravitational softening length of 1 h⁻¹ kpc (physical) at low redshifts. It was run using the Adaptive Refinement Tree code (Kravtsov et al. 1997). It assumes a flat ΛCDM cosmology, with cosmological parameters (${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.73$, Ω_m = 0.27, Ω_b = 0.047, n = 0.95, σ₈ = 0.82, and h = 0.70). Bolshoi provides us with clusters of higher mass resolution than MDPL2. Thus, we use both simulations to test the efficiency of GalWeight to recover the true membership.

For both simulations, halos are identified using the bound density maximum (BDM) algorithm (Klypin & Holtzman 1997; Riebe et al. 2013) that was extensively tested (e.g., Knebe et al. 2011) and identifies local density maxima, determines a spherical cut off for the halo with overdensity equal to 200 times the critical density of the universe (ρ = 200ρ_c) for MDPL2 and 360 times the background matter density of the universe (ρ = 360ρ_bg), and removes unbound particles from the halo boundary. Among other parameters, BDM provides virial masses and radii. Virial mass is defined as ${M}_{v}=\tfrac{4}{3}\pi 200{\rho }_{c}{r}_{v}^{3}$ for MDPL2 and ${M}_{v}=\tfrac{4}{3}\pi 360{\rho }_{\mathrm{bg}}{r}_{v}^{3}$ for Bolshoi (see Bryan & Norman 1998; Klypin et al. 2016). The halo catalogs are complete for halos with circular velocity ${v}_{c}\geqslant 150$ km s⁻¹ for MDPL2 (Klypin et al. 2016) and ${v}_{c}\geqslant 100$ km s⁻¹ for Bolshoi (e.g., Busha et al. 2011; Klypin et al. 2011; Old et al. 2015).

For both MDPL2 and Bolshoi, the phase-space (line-of-sight velocity v_z versus projected radius R_p) of a distinct halo (cluster) is constructed as follows. We assume the line of sight to be along the z-direction and the projection to be on the x–y plane. We select a distinct halo of coordinates (${x}^{h},{y}^{h},{z}^{h}$) and velocity components (${v}_{x}^{h},{v}_{y}^{h},{v}_{z}^{h}$), and then we calculate the observed line-of-sight velocity of a subhalo, taking the Hubble expansion into account, as ${v}_{z}=({v}_{z}^{g}-{v}_{z}^{h})+{H}_{0}({z}^{g}-{z}^{h})$, where (${x}^{g},{y}^{g},{z}^{g}$) and (${v}_{x}^{g},{v}_{y}^{g},{v}_{z}^{g}$) are the coordinates and velocity components of the subhalo, respectively. Finally, we select all subhalos within a projected radius of⁵ ${R}_{p,\max }\,=10$ h⁻¹ Mpc from the center of distinct halo and within a line-of-sight velocity interval of $| {v}_{z,\max }| =3500$ km ${{\rm{s}}}^{-1}$. These values are chosen to be sufficiently large to exceed both the turnaround radius and the length of the Fingers-of-God (hereafter, FOG) effect, which are typically $\sim 7\mbox{--}8\,{h}^{-1}$ Mpc and ∼6000 km s⁻¹, respectively, for massive clusters. The turnaround radius r_t is the radius at which a galaxy's peculiar velocity (v_pec) is canceled out by the global Hubble expansion. In other words, it is the radius at which the infall velocity vanishes (${v}_{\inf }={v}_{\mathrm{pec}}-H\,r=0$).

In this section, we describe the GalWeight technique in detail and demonstrate its use by applying it interactively to a simulated cluster of mass $9.37\times {10}^{14}$ h⁻¹ M_⊙ selected from the Bolshoi simulation. Figure 1 shows the phase-space distribution of subhalos (galaxies) near the center of the simulated cluster.

Zoom In Zoom Out Reset image size

The GalWeight technique works by assigning a weight to each galaxy i according to its position (${R}_{p,i}$, ${v}_{z,i}$) in the phase-space diagram. This weight is the product of two separate two-dimensional weights, which we refer to as the dynamical and phase-space weights. The dynamical weight (see Section 3.1 parts A.1 and A.2, and Figure 4(a), which is the product of Figures 2(b) and 3(b)) is calculated from the surface number density Σ(R_p), velocity dispersion ${\sigma }_{{v}_{z}}({R}_{p})$, and standard deviation ${\sigma }_{{R}_{p}}({v}_{z})$ profiles of the cluster. The phase-space weight (see Section 3.1 part B and Figure 4(b)) is calculated from the two-dimensional adaptive kernel method that estimates the probability density underlying the data and consequently identification of clumps and substructures in the phase space (Pisani 1996). The total weight is then calculated as the product of the dynamical and phase-space weights (see Section 3.1 part C and Figure 4(c)). The advantage of using the total weight rather than the dynamical weight or the phase-space weight alone is discussed in Section 3.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Galaxy Weighting Functions

A. Dynamical Weighting, ${{ \mathcal W }}_{\mathrm{dy}}({R}_{p},{v}_{z}):$ In calculating the dynamical weighting function, we assume that the weighting we apply should be larger at the cluster center, i.e., at the origin in phase space (Figure 1), and decay along both the R_p and v_z phase space axes. This weighting function is, therefore, a product of two individual weighting functions; one that decays along the R_p-axis and the other along the v_z-axis as described below.

A.1. R_p-axis Weighting Function, ${{ \mathcal W }}_{{R}_{p}}({R}_{p}):$ In order to calculate the projected radius weighting function, ${{ \mathcal W }}_{{R}_{p}}({R}_{p})$, we select two properties that are strongly correlated with projected radius and with the dynamical state of a cluster.

The first property is the surface number density profile Σ(R_p), defined as the number of galaxies per unit area as a function of distance from the cluster center. It has its maximum value at the cluster center and decreases with radial distance, and is also strongly correlated with the mass distribution of the cluster. The significance of introducing Σ(R_p) for calculating ${{ \mathcal W }}_{{R}_{p}}$ is that the velocities of member galaxies in the core of some clusters can be as high as $\approx 3000$ km s⁻¹. It produces the Kaiser or FOG effect (see Kaiser 1987). This FOG distortion is the main reason that many membership techniques fail to correctly identify galaxies in the core with high line-of-sight velocities as members. Thus, Σ(R_p) is essential to recover the members in the cluster core. In other words, ignoring Σ(R_p) means missing some of the cluster members in the core.

The second property is the projected velocity dispersion profile, ${\sigma }_{{v}_{z}}({R}_{p})$. The significance of introducing ${\sigma }_{{v}_{z}}({R}_{p})$ for calculating ${{ \mathcal W }}_{{R}_{p}}$ is that it characterizes the dynamical state of a cluster from its core to its infall region. Specifically, there is random motion of member galaxies in the infall region due to the presence of substructures and recent mergers (e.g., van Haarlem 1992; van Haarlem & van de Weygaert 1993; Diaferio & Geller 1997). This effect of random motion can be taken into account implicitly in ${\sigma }_{{v}_{z}}({R}_{p})$. This is the main reason why the SIM technique fails in the cluster outskirts in the projected phase-space. Thus, ${\sigma }_{{v}_{z}}({R}_{p})$ is essential to recover the members in the cluster infall region. In other words, ignoring ${\sigma }_{{v}_{z}}({R}_{p})$ means missing some of the cluster members in the infall region.

Thus, the weighting ${{ \mathcal W }}_{{R}_{p}}({R}_{p})$ in the projected radius direction can be calculated by introducing the function ${{ \mathcal D }}_{{R}_{p}}({R}_{p})$ that is given by

$$\begin{eqnarray}&&{{ \mathcal D }}_{{R}_{p}}({R}_{p})=\displaystyle \frac{{\rm{\Sigma }}({R}_{p}){\sigma }_{{v}_{z}}({R}_{p})}{{R}_{p}^{\nu }},\end{eqnarray} \tag{ 1 }$$

with the normalization

$$\begin{eqnarray}&&{{ \mathcal N }}_{{R}_{p}}={\int }_{0}^{{R}_{p,\max }}{{ \mathcal D }}_{{R}_{p}}({R}_{p}){{dR}}_{p},\end{eqnarray} \tag{ 2 }$$

where ${R}_{p,\max }$ is the maximum projected radius in phase space. The denominator ${R}_{p}^{\nu }$, where the slope of the power law ν is a free parameter in the range $-1\lesssim \nu \lesssim 1$. The factor ${R}_{p}^{\nu }$ is introduced in Equation (1) to provide flexibility and generalization for the technique. The free parameter ν is selected to adjust the effect of the distortion of FOG in the core and the distortion of the random motion in the outer region. It is defined as

$$\begin{eqnarray}&&\nu =\displaystyle \frac{{\sigma }_{\mathrm{FOG}}(R\leqslant 0.25)}{{\sigma }_{\mathrm{rand}}(0.25\lt R\leqslant 4))}-1,\end{eqnarray} \tag{ 3 }$$

where σ_FOG is the velocity dispersion of the core galaxies and ${\sigma }_{\mathrm{rand}}$ is the velocity dispersion of the galaxies outside the core.

The function ${{ \mathcal D }}_{{R}_{p}}({R}_{p})$, calculated from the data, is contaminated by interlopers that cause scattering, especially at large projected distances (see black points in the left panel of Figure 2). Therefore, in order to apply a smooth weighting function, we fit ${{ \mathcal D }}_{{R}_{p}}({R}_{p})$ with an analytical function. Any analytical function that is a good fit to ${{ \mathcal D }}_{{R}_{p}}({R}_{p})$ could be utilized. In this paper, we choose to use the function

$$\begin{eqnarray}&&{{ \mathcal W }}_{{R}_{p}}({R}_{p})={{ \mathcal A }}_{0}{\left(1+\displaystyle \frac{{R}_{p}^{2}}{{a}^{2}}\right)}^{\gamma }+{{ \mathcal A }}_{\mathrm{bg}},\end{eqnarray} \tag{ 4 }$$

which has four parameters: a is a scale radius ($0\lt a\lesssim 1$), γ is a slope of the power law ($-2\lesssim \gamma \lt 0$), and ${{ \mathcal A }}_{0}$ and ${{ \mathcal A }}_{\mathrm{bg}}$ are the central and background weights in the R_p-direction. These parameters are determined by applying the chi-squared algorithm using the Curve Fitting MatLab Toolbox. Note that the analytical function we selected here has the same form as the generalized King model (King 1972; Adami et al. 1998).

Thus, the weight ${{ \mathcal W }}_{{R}_{p}}({R}_{p,i})$ of each galaxy can be calculated according to its projected radius R_p from the cluster center. The weighting along R_p is shown in Figure 2(a), where the function ${{ \mathcal D }}_{{R}_{p}}({R}_{p})$ is normalized using Equation (2). The data are smoothed and approximated using Equation (4) (shown by the red line). The right panel (b) shows the projected radius weight function in phase space.

A.2. v_z-axis Weighting Function, ${{ \mathcal W }}_{{v}_{z}}({v}_{z}):$ In phase-space, most members are concentrated near the line v_z = 0 and the number of members decreases with increasing absolute line-of-sight velocity. The weighting function along the v_z-axis can, therefore, be approximated by the histogram of the number of galaxies per bin, ${N}_{\mathrm{bin}}({v}_{z})$, or equivalently the standard deviation of projected radius, ${\sigma }_{{R}_{p}}({v}_{z})$, directed along the line-of-sight velocity axis, normalized by the total number of galaxies N_tot in the cluster field. In other words, the weighting in the line-of-sight velocity direction can be calculated by introducing the function ${{ \mathcal D }}_{{v}_{z}}({v}_{z})$ that is given by

$$\begin{eqnarray}&&{{ \mathcal D }}_{{v}_{z}}({v}_{z})={\sigma }_{{R}_{p}}({v}_{z}),\end{eqnarray} \tag{ 5 }$$

with the normalization

$$\begin{eqnarray}&&{{ \mathcal N }}_{{v}_{z}}={\int }_{-{v}_{z,\max }}^{{v}_{z,\max }}{{ \mathcal D }}_{{v}_{z}}({v}_{z}){{dv}}_{z},\end{eqnarray} \tag{ 6 }$$

where ${v}_{z,\max }$ is the maximum line-of-sight velocity of phase space. As above, to obtain a smooth weighting function in v_z, the histogram or equivalently ${{ \mathcal D }}_{{v}_{z}}({v}_{z})$ can be fitted by an analytical function. In this paper, we select an exponential model that is given by

$$\begin{eqnarray}&&{{ \mathcal W }}_{{v}_{z}}({v}_{z})={{ \mathcal B }}_{0}\exp (b\ {v}_{z})+{{ \mathcal B }}_{\mathrm{bg}},\end{eqnarray} \tag{ 7 }$$

where ${{ \mathcal B }}_{0}$ is the central weight, ${{ \mathcal B }}_{\mathrm{bg}}$ is the background weight in v_z, and b is scale parameter ($-0.01\lesssim b\lt 0$). Again, these parameters are determined by applying the chi-squared algorithm using the Curve Fitting MatLab Toolbox. The weighting along v_z is shown in Figure 3(a), where the function ${{ \mathcal D }}_{{v}_{z}}({v}_{z})$ (black points) is normalized using Equation (6). The data are smoothed and approximated by Equation (7) for an exponential model (blue curve). The right panel (b) shows the resulting exponential-model weight as a function of location in line-of-sight velocity/projected radius phase space. As shown in (a) and (b), the applied weight is greatest at v_z = 0 and decreases as the absolute line-of-sight velocity increases.

We can now construct a two-dimensional dynamical weight ${{ \mathcal W }}_{\mathrm{dy}}({R}_{p},{v}_{z})$ by multiplying ${{ \mathcal W }}_{{R}_{p}}({R}_{p})$ and ${{ \mathcal W }}_{{v}_{z}}({v}_{z})$ together:

$$\begin{eqnarray}&&{{ \mathcal W }}_{\mathrm{dy}}({R}_{p},{v}_{z})={{ \mathcal W }}_{{R}_{p}}({R}_{p}){{ \mathcal W }}_{{v}_{z}}({v}_{z}),\end{eqnarray} \tag{ 8 }$$

${{ \mathcal W }}_{\mathrm{dy}}({R}_{p},{v}_{z})$ is shown in the left panel of Figure 4, and is the product of the weights shown in Figures 2(b) and 3(b). The weight is maximum at the origin and decreases along both R_p and v_z.

To sum up, the dynamical weight is calculated from three properties (surface number density Σ(R_p) and velocity dispersion ${\sigma }_{{v}_{z}}({R}_{p})$ along R_p, and standard deviation of projected radius ${\sigma }_{{R}_{p}}({v}_{z})$ along v_z), which are correlated strongly with the dynamics of the cluster. This weight takes into account the effects of the FOG in the cluster core and the random motion of galaxies in the infall region.

B. Phase-Space Weighting, ${{ \mathcal W }}_{\mathrm{ph}}({R}_{p},{v}_{z}):$ This weighting is the coarse-grained phase-space density, which can be simply calculated by the two-dimensional adaptive kernel method (2DAKM, e.g., Silverman 1986; Pisani 1996). The kernel density estimator is the estimated probability density function of a random variable. For N galaxies with coordinates (x, y) = (R_p, v_z), the density estimator is given by

$$\begin{eqnarray}&&f(x,y)=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}\displaystyle \frac{1}{{h}_{x,i}{h}_{y,i}}K\left(\displaystyle \frac{x-{X}_{i}}{{h}_{x,i}}\right)K\left(\displaystyle \frac{y-{Y}_{i}}{{h}_{y,i}}\right),\end{eqnarray} \tag{ 9 }$$

where the kernel K(t) is given by Gaussian distribution

$$\begin{eqnarray}&&K(t)=\displaystyle \frac{1}{\sqrt{(2\pi )}}\exp \left(-\displaystyle \frac{1}{2}{t}^{2}\right)\end{eqnarray} \tag{ 10 }$$

and ${h}_{j,i}={\lambda }_{i}{h}_{j}$ is the local bandwidth, ${h}_{j}={\sigma }_{j}{N}^{-1/6}$ is the fixed bandwidth for two-dimensional space and ${\sigma }_{j}$ is the standard deviation for j = {x, y}. The term ${\lambda }_{i}={[\gamma /{f}_{0}({x}_{i},{y}_{i})]}^{0.5}$ and $\mathrm{log}\gamma ={\sum }_{i}\mathrm{log}{f}_{0}({x}_{i},{y}_{i})/N$, where ${f}_{0}({x}_{i},{y}_{i})$ is given by Equation (9) for ${\lambda }_{i}=1$ (see also, Shimazaki & Shinomoto 2010).

Consequently, applying 2DAKM for the distribution of galaxies in the phase-space demonstrates high weights for positions of high-density distribution of galaxies. Therefore, the main purpose of introducing the phase-space weight is to take into account the effect of the presence of any clump or substructure in the field that cannot be counted by the dynamical weight. Also, the phase-space weight is introduced to reduce the excessive increase of dynamical weight near the center (see Section 3.3). The phase-space weight ${{ \mathcal W }}_{\mathrm{ph}}({R}_{p},{v}_{z})$ is shown in Figure 4(b) that gives more weights for galaxies in clumps and substructures, and from the distribution of galaxies in the cluster field this weighting function is maximum around the cluster center.

C. Total Weighting, ${{ \mathcal W }}_{\mathrm{tot}}({R}_{p},{v}_{z}):$ The total weighting function is calculated as

$$\begin{eqnarray}&&{{ \mathcal W }}_{\mathrm{tot}}({R}_{p},{v}_{z})={{ \mathcal W }}_{\mathrm{dy}}({R}_{p},{v}_{z}){{ \mathcal W }}_{\mathrm{ph}}({R}_{p},{v}_{z}),\end{eqnarray} \tag{ 11 }$$

and shown in Figure 4(c) for the simulated cluster. It shows the probability distribution function (PDF) of the total weight ${{ \mathcal W }}_{\mathrm{tot}}({R}_{p},{v}_{z})$. The weighting in Figure 4(c) is then applied to individual galaxies. Figure 5(a) shows Figure 1 once again, but now after applying the "total weighting." We still need to separate cluster members from interlopers. We explain how to do that in Section 3.2.

Zoom In Zoom Out Reset image size

3.2. Membership Determination

3.3. Why Do We Use Total Weight Rather Than Dynamical or Phase-space Weights?

One may ask why we depend on the total weight to assign cluster membership rather than using the dynamical weight or phase-space weight alone. We present Figure 6 to help answer this question. It shows the phase space of three Bolshoi simulated clusters (see Section 2). Using simulated clusters provides the advantage that true members are known definitively. Figure 6 shows the optimal contour lines determined by applying, separately, the dynamical weight (green line), the phase-space weight (blue line), and the total weight (black line). The red points show true members within 3r_v.

Zoom In Zoom Out Reset image size

In Figure 6, the dynamical weight ${{ \mathcal W }}_{\mathrm{dy}}({R}_{p},{v}_{z})$ (green; see also Figure 4(a)) is seen to be very smooth and idealized. In other words, ${{ \mathcal W }}_{\mathrm{dy}}({R}_{p},{v}_{z})$ describes well an isolated galaxy cluster in phase-space. It does not take into account the effects of nearby clusters, clumps, and/or substructures. Also, it shows an excessive increase near the cluster center ($\sim 1\,{h}^{-1}$ Mpc) and incorrectly includes interlopers near the center that have very high velocities. This effect is due to introducing Σ(R_p) in ${{ \mathcal W }}_{\mathrm{dy}}({R}_{p},{v}_{z})$, where the surface number density is very high near the cluster center. However, ignoring ${{ \mathcal W }}_{\mathrm{dy}}({R}_{p},{v}_{z})$ leads to missing some cluster members, especially those that close to the center in phase space. Thus, ${{ \mathcal W }}_{\mathrm{dy}}({R}_{p},{v}_{z})$ cannot be used on its own to assign cluster membership, but it is very important for correctly identifying members with high line-of-sight velocities.

Figure 6 demonstrates that, on its own, phase-space weighting ${{ \mathcal W }}_{\mathrm{ph}}({R}_{p},{v}_{z})$ also has some difficulty in recovering true cluster members (blue; see also Figure 4(b)). This is because it does not take into account the FOG effect in the cluster core, where those members that have high velocities do not have high concentration, so they are assigned low weights in phase space and not counted as members. Also, the presence of nearby clusterings and substructures have the effect of widening the "optimal" contour line. Consequently, it is very difficult to separate true members from galaxies (interlopers) located in nearby clumps. This results in the inclusion of some interlopers in the infall region. In summary, using ${{ \mathcal W }}_{\mathrm{ph}}({R}_{p},{v}_{z})$ alone, simultaneously excludes some true members near the cluster center and includes some interlopers in the infall region.

We have shown that both the dynamical weight and phase-space weight have issues in identifying true members when applied alone. However, as the black solid line in Figure 6 shows, the total weight (the product of the dynamical and phase-space weights), is very effective. It can simultaneously identify cluster members moving with high velocities in the core (${R}_{p}\lesssim 1$ Mpc h⁻¹) as well as members moving with random motions in the infall regions (${R}_{p}\sim 3{r}_{v}$).

3.4. Testing the Efficiency of GalWeight on MDPL and Bolshoi Simulations

To further demonstrate and quantify the GalWeight technique when assigning membership, we again utilize the MDPL2 and Bolshoi⁶ simulations from the suite of MultiDark simulations. The efficiency of GalWeight can be quantified by calculating two fractions defined as follows. The first is the completeness f_c, which is the fraction of the number of fiducial members identified by GalWeight as members in the projected phase-space relative to the actual number of 3D members projected in the phase-space. The second is the contamination f_i, which is the fraction of interlopers incorrectly assigned to be members, projected in the phase space (see, e.g., Wojtak et al. 2007; Serra & Diaferio 2013). Ideally, of course, GalWeight would return fractions of f_c = 1 and f_i = 0.

MDPL2 provides us with 1500 simulated clusters with masses ranging from $0.73\times {10}^{14}\,{h}^{-1}\,{M}_{\odot }$ to $37.4\times {10}^{14}\,{h}^{-1}\,{M}_{\odot }$ to which we can apply GalWeight. We calculate the fractions f_c and f_i at three radii—${r}_{v},2{r}_{v}$, and 3r_v. As shown in Table 1, the mean values of f_c and f_i within r_v are 0.993 and 0.112, respectively, for the 1500 clusters overall. Also, the fraction f_c decreases from 0.993 at r_v to 0.981 at 3r_v.

For Bolshoi, we have about 500 clusters with masses greater than $0.70\times {10}^{14}\,{h}^{-1}\,{M}_{\odot }$. In order to increase the cluster sample of Bolshoi to 1500 clusters, we randomly select different lines of sight or ordinations for each distinct halo in additional to the original line of sight along the z-direction (see column 3 in Table 1 for Bolshoi). Then, we apply GalWeight to each cluster. The mass range of the sample is $0.70\times {10}^{14}\,{h}^{-1}\,{M}_{\odot }$ to $10.92\times {10}^{14}\,{h}^{-1}\,{M}_{\odot }$ as shown in Table 1. The mean values of f_c and f_i within r_v are 0.995 and 0.126, respectively, for the 1500 clusters overall. Also, the fraction f_c decreases from 0.995 at r_v to 0.971 at 3r_v.

The main reason that some interlopers are assigned as members (f_i = 0.113 for MDPL2 and f_i = 0.226 for Bolshoi, as maximal value) is because of the triple-value problem (Tonry & Davis 1981). That is, there are some foreground and background interlopers that appear to be part of the cluster body because of the distortion of phase-space. The effect of the triple-value problem is apparent in Figure 7 (discussed below), where most of the interlopers assigned as members are embedded in the cluster body. We defer a discussion of how GalWeight may be adapted to overcome this problem to a future work.

Zoom In Zoom Out Reset image size

In order to demonstrate the ability of GalWeight to assign membership in the case of both poor and massive clusters, we divide the 1500 clusters (for each simulation) into four mass bins as shown in Table 1. The fraction f_c varies from 0.998 (0.996) for the poor clusters of mean mass $1.44\,\times {10}^{14}\,(1.13\times {10}^{14})\,{h}^{-1}\,{M}_{\odot }$ to 0.988 (0.997) for the massive clusters of mean mass $11.34\times {10}^{14}\,(9.68\times {10}^{14})\,{h}^{-1}\,{M}_{\odot }$ at r_v for MDPL2 (Bolshoi). We conclude that GalWeight can be applied effectively to a range of cluster masses with high efficiency.

Figure 7 shows examples of GalWeight being applied to 12 simulated Bolshoi clusters (solid black lines), where red and gray points show fiducial members and interlopers, respectively, within 3r_v. The 12 clusters shown in Figure 7 are ranked by virial mass, with the most massive cluster ($10.92\times {10}^{14}$ h⁻¹ M_⊙) shown in the upper left corner and the least massive one ($1.06\times {10}^{14}$ h⁻¹ M_⊙) shown in the lower right corner. The figure demonstrates that GalWeight can effectively recover cluster membership for rich massive galaxy clusters as well as small or poor groups of galaxies with the same efficiency.

In summary, applying GalWeight to the suite of MDPL2 and Bolshoi simulations demonstrates that GalWeight can successfully recover cluster membership with high efficiency. It also further demonstrates that it can simultaneously identify members in both the virial and infall regions by taking into account the FOG effect and the random motion of galaxies in the infall region. Furthermore, it can be applied to both rich galaxy clusters and poor groups of galaxies with the same efficiency (see Table 1).

In this section, we perform a general comparison between GalWeight and four other well-known techniques (shifting gapper, caustic, den Hartog, and SIM) without doing any quantitative comparison. We defer testing the efficiency of different membership techniques to recover the 3D true members of clusters and the influence of the determining their dynamical masses to a future work (see, e.g., Wojtak et al. 2007).

We begin by showing how each technique fares when it is applied in turn to two simulated clusters with masses of $10.92\times {10}^{14}$ h⁻¹ M_⊙ and $4.24\times {10}^{14}$ h⁻¹ M_⊙ from the Bolshoi simulation, shown in Figure 8. Making the assumption that the cluster is spherical, fiducial members are assumed to lie within three virial radii, 3r_v, and are shown as 2D members in the phase-space (red points) in each panel of Figure 8. We select this radius (3r_v) in order to examine the ability of each technique to recover true members not only within the virial radius but also in the infall region, i.e., the region of a cluster that extends from the viral radius r_v to the turnaround radius r_t, where ${r}_{t}\sim \,2\mbox{--}4\,{r}_{v}$. Shown in each panel by the solid black line is the optimal choice of demarcation contour separating members and field galaxies identified by our GalWeight technique. For reasons of space, we do not describe each of the four techniques (shifting gapper, den Hartog, caustic, and SIM) in detail here. However, we summarize them below and refer the reader to the references for more information.

Zoom In Zoom Out Reset image size

The shifting gapper technique (Fadda et al. 1996) works by first placing galaxies into bins according to their projected radial distance from the cluster center. The user has the freedom to choose the number of galaxies per bin that they believe is best-suited to each application of the technique. Commonly chosen values are N_bin = 10 or 15. For each bin, the galaxies are sorted according to their velocities, then any galaxy separated by more than a fixed value (e.g., 1σ of the sample or 500–1000 km s⁻¹) from the previous one is considered an interloper and removed. Fadda et al. (1996) used a gap of 1000 km s⁻¹ and a bin of 0.4 h⁻¹ Mpc or larger in order to have at least 15 galaxies. The open blue circles in panels (a), (e) and (b), (f) of Figure 8 represent the members identified by this technique, where the number of galaxies utilized per bin was N_bin = 10 and N_bin = 15, respectively. The gray points symbolize interlopers. Clearly, membership identification depends heavily upon the choice of N_bin, as there are many differences between the galaxies identified as members in panels (a), (e) and (b), (f). Additionally, in both cases, some true members of the two clusters are missed, especially at small projected radius. Furthermore, the shifting gapper technique depends on the choice of the velocity gap used to remove interlopers in each bin. A choice of a high-velocity gap results in the identification of a large fraction of interlopers as cluster members, while the choice of a low-velocity gap results in missing true cluster members (Aguerri et al. 2007).

The application of the caustic technique (e.g., Alpaslan et al. 2012; Serra & Diaferio 2013) is shown in panels (c) and (g) of Figure 8 for two rescale parameters, q = 25 (cyan lines) and q = 35 (pink lines). Although this technique is quite successful when applied to the cluster outskirts, it misses some of the true members located within the core, which are the most important galaxies affecting the dynamics of the clusters. They are missed because the caustic technique does not take into account the effects of the FOG distortion. Also, the caustic technique cannot be applied to small galaxy groups. Furthermore, applying the caustic technique is rather subjective and relies upon how the caustics can be inferred from the data (Reisenegger et al. 2000; Pearson et al. 2014). Nonetheless, it is still a powerful technique for estimating cluster masses.

The application of the den Hartog technique (den Hartog & Katgert 1996) is also shown by the dotted black lines in Figure 8, panels (c) and (g). This technique estimates the escape velocity as a function of distance from the cluster center by calculating the virial mass profile (see Section 5), ${v}_{\mathrm{esc}}(R)=\sqrt{\tfrac{2{{GM}}_{\mathrm{vir}}(R)}{R}}$, where G is the gravitational constant. The figure demonstrates that this technique is very biased toward including many far interlopers. In addition, its application relies on assumptions of hydrostatic equilibrium and spherical symmetry.

Panels (d) and (h) in Figure 8 show the application of two spherical infall models (SIMs). The Yahil (1985) and Regős models (Regős & Geller 1989) are shown by dark green and light green lines respectively. Note that one needs to determine the mass density profile and the background mass density in order to apply the SIM technique and determine the infall velocity profile (e.g., van Haarlem & van de Weygaert 1993). We determine the mass density profile for the simulated cluster from the NFW model (Navarro et al. 1996, 1997, Equations (15) and (16)), knowing its concentration c, virial radius r_v, and scale radius ${r}_{s}={r}_{v}/C$. Also, the background mass density is given by ${\rho }_{\mathrm{bg}}={{\rm{\Omega }}}_{m}\,{\rho }_{c}$.

As shown in Figures 8(d) and (h), SIMs have difficulty identifying true members in the infall region in projected phase space. This is due to the fact that the effect of random motion of galaxies in the infall regions (van Haarlem & van de Weygaert 1993; Diaferio 1999) causes some members in the cluster outskirts to be missed. A recent investigation by our own team (Abdullah et al. 2013) has shown that SIMs can successfully be applied to sliced phase space by taking into account some kinds of distortions such as the transverse motion of galaxies with respect to the observer and/or rotational motions of galaxies inside the cluster. However, this is out of the scope of the current paper.

In this section, we apply GalWeight to a sample of 12 Abell galaxy clusters, with galaxy coordinates and redshifts taken from SDSS-DR12⁷ (hereafter, SDSS-DR12 Alam et al. 2015). In order to demonstrate the technique for both massive and poor clusters, we selected clusters with Abell richness parameters ranging from 0 to 3 (Abell et al. 1989). We deliberately selected some clusters that were almost isolated and others that had clumps or groups of galaxies nearby in order to demonstrate how the technique performs under these different scenarios. We apply the GalWeight technique only to this pilot sample of 12 clusters in this paper, deferring application to the entire SDSS-DR13 sample of ∼800 clusters to a later paper.

The data sample is collected as follows. The NASA/IPAC Extragalactic Database (NED)⁸ provides us with a first approximation of the angular coordinates and redshift of the center of our cluster sample (${\alpha }_{c}$, ${\delta }_{c}$, z_c). We then download the coordinates and redshifts (right ascension α, declination δ, and spectroscopic redshift z) for objects classified as galaxies near the center of each cluster from SDSS-DR12 (Alam et al. 2015). The next step is to apply the binary tree algorithm (e.g., Serra et al. 2011) to accurately determine the cluster center (${\alpha }_{c}$, ${\delta }_{c},{z}_{c}$) and create a line-of-sight velocity (v_z) versus projected radius (R_p) phase-space diagram. R_p is the projected radius from the cluster center and v_z is the line-of-sight velocity of a galaxy in the cluster frame, calculated as ${v}_{z}=c(z-{z}_{c})/(1+{z}_{c})$, where z is the observed spectroscopic redshift of the galaxy and z_c is the cluster redshift. The term $(1+{z}_{c})$ is a correction due to the global Hubble expansion (Danese et al. 1980) and c is the speed of light.

We then apply GalWeight to the 12 Abell clusters as described in detail in Section 3 in order to get the optimal contour line. The final step is to determine the virial radius, r_v, at which ρ = 200ρ_c and the turnaround radius, r_t, at which ρ = 5.55ρ_c (e.g., Nagamine & Loeb 2003; Busha et al. 2005; Dünner et al. 2006) from all galaxies located inside the optimal contour line of a cluster

In order to calculate these two radii, we should first determine the cluster mass profile. The cluster mass can be estimated from the virial mass estimator and NFW mass profile (Navarro et al. 1996, 1997) as follows.

The viral mass estimator is given by

$$\begin{eqnarray}&&M(\lt r)=\displaystyle \frac{3\pi N{\displaystyle \sum }_{i}{v}_{z,i}{(\lt r)}^{2}}{2G{\displaystyle \sum }_{i\ne j}\tfrac{1}{{R}_{{ij}}}},\end{eqnarray} \tag{ 12 }$$

where ${v}_{z,i}$ is the galaxy line-of-sight velocity and R_ij is the projected distance between two galaxies (e.g., Limber & Mathews 1960; Binney & Tremaine 1987; Rines et al. 2003).

If a system extends beyond the virial radius, Equation (12) will overestimate the mass due to external pressure from matter outside the virialized region (The & White 1986; Carlberg et al. 1997; Girardi et al. 1998). The corrected virial mass can be determined using the following expression:

$$\begin{eqnarray}&&{M}_{v}(\lt r)=M(\lt r)[1-S(r)],\end{eqnarray} \tag{ 13 }$$

where S(r) is a term introduced to correct for surface pressure. For an NFW density profile and for isotropic orbits (i.e., the projected, ${\sigma }_{v}$, and angular, ${\sigma }_{\theta }$, velocity dispersion components of a galaxy in the cluster frame are the same, or equivalently the anisotropy parameter $\beta =1-\tfrac{{\sigma }_{\theta }^{2}}{{\sigma }_{r}^{2}}=0$), S(r) can be calculated by

$$\begin{eqnarray}&&S(r)={\left(\displaystyle \frac{x}{1+x}\right)}^{2}{\left[\mathrm{ln}(1+x)-\displaystyle \frac{x}{1+x}\right]}^{-1}{\left[\displaystyle \frac{{\sigma }_{v}(r)}{\sigma (\lt r)}\right]}^{2},\end{eqnarray} \tag{ 14 }$$

where $x=r/{r}_{s}$, r_s is the scale radius, $\sigma (\lt r)$ is the integrated three-dimensional velocity dispersion within r, and ${\sigma }_{v}(r)$ is a projected velocity dispersion (e.g., Koranyi & Geller 2000; Abdullah et al. 2011).

The mass density within a sphere of radius r introduced by NFW is given by

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{\rho }_{s}}{x{(1+x)}^{2}},\end{eqnarray} \tag{ 15 }$$

and its corresponding mass is given by

$$\begin{eqnarray}&&M(\lt r)=\displaystyle \frac{{M}_{s}}{\mathrm{ln}(2)-(1/2)}\left[\mathrm{ln}(1+x)-\displaystyle \frac{x}{1+x}\right],\end{eqnarray} \tag{ 16 }$$

where ${M}_{s}=4\pi {\rho }_{s}{r}_{s}^{3}[\mathrm{ln}(2)-(1/2)]$ is the mass within r_s, ρ_s = δ_sρ_c is the characteristic density within r_s and ${\delta }_{s}=({{\rm{\Delta }}}_{v}/3){c}^{3}{\left[\mathrm{ln}(1+c)-\tfrac{c}{1+c}\right]}^{-1}$, and the concentration $c={r}_{v}/{r}_{s}$ (e.g., Navarro et al. 1997; Rines et al. 2003; Mamon et al. 2013).

The projected number of galaxies within a cylinder of radius R is given by integrating the NFW profile (Equation (15)) along the line of sight (e.g., Bartelmann 1996; Zenteno et al. 2016)

$$\begin{eqnarray}&&N(\lt R)=\displaystyle \frac{{N}_{s}}{\mathrm{ln}(2)-(1/2)}g(x),\end{eqnarray} \tag{ 17 }$$

where N_s is the number of galaxies within r_s that has the same formula as M_s, and g(x) is given by (e.g., Golse & Kneib 2002; Mamon et al. 2010)

$$\begin{eqnarray}g(x)=\left\{\begin{array}{ll}\mathrm{ln}(x/2)+\tfrac{{\cosh }^{-1}(1/x)}{\sqrt{1-{x}^{2}}} & {\rm{if}}\ x\ \lt \ 1\\ 1-\mathrm{ln}(2) & {\rm{if}}\ x=1\\ \mathrm{ln}(x/2)+\tfrac{{\cos }^{-1}(1/x)}{\sqrt{{x}^{2}-1}} & {\rm{if}}\ x\ \gt \ 1.\end{array}\right.\end{eqnarray} \tag{ 18 }$$

Thus, we can fit r_s for each cluster to get S(r) from Equation (14) and calculate the corrected mass profile M_v(r) from Equation (13). Also, the NFW mass profile is calculated from Equation (16). Then, r_v, at which Δ = 200ρ_c, can be calculated from the viral or NFW mass profiles. While r_t, at which Δ = 5.55ρ_c, can be determined from NFW mass profile only. We cannot determine r_t from the virial mass profile because the assumption of hydrostatic equilibrium is invalid.

Finally, after we calculate r_v and r_t (from NFW mass profile) the cluster membership can be defined as all galaxies enclosed by the optimal contour line and within r_t, as shown by the red points in Figure 9. It is worth noting once again that GalWeight is effective at taking into account the effects of the FOG distortion in the innermost regions and the random motion of galaxies in the cluster infall region. Moreover, GalWeight is not affected by the presence of substructures or nearby clusters or groups as demonstrated, for example, for A2063 and A2065. Furthermore, GalWeight can be applied both to rich clusters such as A2065 and A1656 and to poor galaxy groups such as A1459 and A2026.

Zoom In Zoom Out Reset image size

In order to compare our results with the literature, we calculate the radii and their corresponding masses at three overdensities, Δ₅₀₀ = 500ρ_c, Δ₂₀₀ = 200ρ_c, and Δ₁₀₀ = 100ρ_c as shown in Table 2. The sample is displayed in order of decreasing NFW M₂₀₀ mass. A complete list of NFW parameters is also provided in Table 2.

In Table 3, we list ratios of radii and masses for each of the 12 Abell clusters using our GalWeight-determined method (assuming an NFW profile) divided by previously published values, $({r}_{\mathrm{NFW}}/{r}_{\mathrm{ref}})$ and $({M}_{\mathrm{NFW}}/{M}_{\mathrm{ref}})$ respectively, at overdensities of Δ = 500, 200, and 100ρ_c. Column 8 of Table 3 also lists the ratio of GalWeight-determined masses relative to those estimated from the caustic technique (Rines et al. 2016), (${M}_{\mathrm{NFW}}/{M}_{\mathrm{caus}}{)}_{200}$, at Δ = 200ρ_c. Table 3 clearly shows that the radii and masses estimated for a cluster are strongly dependent on the technique used to assign membership and remove interlopers (see Wojtak et al. 2007). The ratio $({r}_{\mathrm{NFW}}/{r}_{\mathrm{ref}})$ ranges between 0.63 and 1.55, while the ratio $({M}_{\mathrm{NFW}}/{M}_{\mathrm{ref}})$ ranges between 0.58 and 2.18.

The cluster masses from the literature tabulated in Table 3 have been calculated in various ways. Below, we explicitly compare our values to those obtained from applying the shifting gapper, SIM, and caustic methods.

First, comparing to the shifting gapper technique (see (${}^{\bullet }$) in Table 3, Girardi et al. 2002; Sifón et al. 2015), we find that the ratio $({M}_{\mathrm{NFW}}/{M}_{\mathrm{ref}})$ is larger than unity in some cases (A2065, A1185) and smaller than unity in others (A2029, A2142). This is because members assigned by this technique, and consequently the mass calculated, depend on the selection criteria of number of galaxies and velocity gap per bin. As discussed before, the choice of a high-velocity gap includes more members and consequently larger mass and vice versa.

Second, comparing to the SIM method (see (^⋆) in Table 3, Abdullah et al. 2011), we note that the mass ratio $({M}_{\mathrm{NFW}}/{M}_{\mathrm{ref}})$ is less than unity for the three clusters A0117, A1436, and A1459. This is because SIM includes more galaxy members inside the virial region even though they are very far from the cluster body. This is due to the assumption of conservation of mass that influences the validity of SIM in the innermost region (see Figure 6 in Abdullah et al. 2011).

Third, comparing to the caustic technique (see (⁺) in Table 3, Rines et al. 2003; Rines & Diaferio 2006; Rines et al. 2016), we specifically calculate the ratio ${({M}_{\mathrm{NFW}}/{M}_{\mathrm{caus}})}_{200}$ as listed in Table 3, column 8. It demonstrates that this ratio is larger than unity for seven clusters, and the highest ratio is for A2065, for which the estimated mass from NFW is four times that expected from the caustic technique. As described above, the main reason for this discrepancy is that the caustic technique does not take into consideration the effect of FOG. Thus, it misses more members inside the virial region and consequently expects lower masses.

We compare again GalWeight with the four well-known techniques (shifting gapper, caustic, den Hartog, and SIM) for the Coma cluster as shown in Figure 10. The figure (see also Figure 8) demonstrates that the GalWeight performs very favorably against established methods, taking into account as it does the effects of the FOG distortion at small projected radius well as the random motion of galaxies in the infall region. In order to apply SIM to the Coma cluster, the spatial number density profile is calculated from the NFW model (Navarro et al. 1996, 1997). Also, we assume that the background number density ${\rho }_{\mathrm{bg}}=0.0106$ h³ ${{\rm{Mpc}}}^{-3}$, which is calculated using the parameters of the Schechter luminosity function (${\phi }^{* }=0.0149$ h³ ${{\rm{Mpc}}}^{-3}$, ${M}^{* }-5\mathrm{log}h=-20.44$, and α = −1.05 for r magnitude, Blanton et al. 2003).

Zoom In Zoom Out Reset image size

Because of the presence of interlopers, estimates of cluster mass tend to be biased too high and estimates of cluster concentration tend to be biased too low. Our work suggests that applying GalWeight rather than another technique to determine cluster membership before applying a dynamical mass estimator (virial theorem, NFW model etc.) likely results in a more accurate estimate of the true cluster mass and concentration. In a future work, we will compare the efficiency of different membership techniques to assign membership and their influence on estimating cluster mass using different mass estimators.

In this paper, we introduced the Galaxy Weighting Function Technique (GalWeight), a powerful new technique for identifying cluster members. specifically designed to simultaneously maximize the number of bona fide cluster members, while minimizing the number of contaminating interlopers.

GalWeight takes into account the causes of different distortions in the phase-space diagram and is independent of statistical or selection criteria. It can recover membership in both the virial and infall regions with high accuracy and is minimally affected by substructure and/or nearby clusters.

We first demonstrated GalWeight's use by applying it interactively to a simulated cluster of mass $9.37\times {10}^{14}$ h⁻¹ M_⊙ selected from Bolshoi simulation. Next, we tested the efficiency of the technique on ∼3000 clusters selected from the MDPL2 and Bolshoi simulations with masses ranging from $0.70\times {10}^{14}\,{h}^{-1}\,{M}_{\odot }$ to $37.4\times {10}^{14}\,{h}^{-1}\,{M}_{\odot }$. The completeness and interloper fractions for MDPL2 are ${f}_{c}=0.993,0.992$ and 0.981 and ${f}_{i}=0.096,0.098$ and 0.118, while for Bolshoi ${f}_{c}=0.995,0.981$ and 0.971 and ${f}_{i}=0.126,0.217$ and 0.226 within r_v, 2r_v, and 3r_v, respectively. We then compared its performance to four well-known existing cluster membership techniques (shifting gapper, den Hartog, caustic, and SIM). Finally, we applied GalWeight to a sample of 12 Abell clusters of varying richnesses taken from SDSS-DR12. By assuming an NFW model and applying the virial mass estimator, we determined the radius and corresponding mass at overdensities of Δ₅₀₀, Δ₂₀₀, and Δ₁₀₀. The virial mass (at Δ₂₀₀) of the sample ranged from $0.82\times {10}^{14}$ h⁻¹ M_⊙ to $12.97\times {10}^{14}$ h⁻¹ M_⊙, demonstrating that GalWeight is effective for poor and massive clusters. In the future, we plan to apply GalWeight to a larger SDSS sample of galaxy clusters at low and high redshifts.

We believe that GalWeight has the potential for astrophysical applications far beyond the identification of cluster members, e.g., identifying stellar members of nearby dwarf galaxies, or separating star-forming and quiescent galaxies. We also plan to investigate these applications in a future work.

We thank Brian Siana for useful discussions. We also thank Gary Mamon for his useful comments on the paper. Finally, we thank the reviewer for suggestions that improved this paper. G.W. acknowledges financial support for this work from NSF grant AST-1517863 and from NASA through programs GO-13306, GO-13677, GO-13747, and GO-13845/14327 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555, and grant number 80NSSC17K0019 issued through the Astrophysics Data Analysis Program (ADAP).