© ESO, 2017

1. Introduction

In the era of precision cosmology, the cluster mass function occupies a central role owing to its dependence on the geometrical properties of the universe and the growth of structures (e.g. Vikhlinin et al. 2009; Mantz et al. 2010b; Allen et al. 2011; Böhringer et al. 2014; Planck Collaboration XXIV 2016). The main ingredient to derive cosmological constraints from galaxy clusters is a robust determination of their mass, which can only be achieved individually on small samples. The practical analysis of the cluster mass function thus relies on indirect estimates via mass-observable scaling relations predicted by the hierarchical model of structure formation (Kaiser 1986). Considerable efforts have been made to calibrate such relations (e.g. Böhringer et al. 2012; or Giodini et al. 2013, for a review), since departures from the theoretical predictions provide insight on the non-gravitational physical processes at play in galaxy clusters (e.g. Voit 2005, for a review).

Estimating the mass of a cluster is a complicated task since its main component, the dark matter halo, is not directly observable. Instead, one can infer its mass via its gravitational lensing distortion on the shape of distant background galaxies, or via its effects on the other components of a cluster, such as the dynamics of the cluster galaxies or the X-ray emission of the intra cluster medium. Each of these approaches requires various hypothesis and suffers from different systematics. For instance, non-thermal pressure support, for example, due to turbulence, cosmic rays, or feedback from active galactic nuclei, can lead to a substantial negative bias of 10 − 30% in the hydrostatic estimates (e.g. Lau et al. 2009; Rasia et al. 2012; Nelson et al. 2014). Weak gravitational lensing does not assume dynamical equilibrium but is subject to projection effects, in particular from the triaxial shape of the dark matter halo, which can lead to overestimated masses by up to ~ 40% (e.g. Corless et al. 2009; Becker & Kravtsov 2011; Feroz & Hobson 2012).

The comparison of different mass estimators has been the subject of many studies, in particular between X-ray hydrostatic and lensing measurements. The topic regained interest with the tension found between the cosmological constraints derived by the Planck Collaboration from cluster counts (Böhringer et al. 2014; Planck Collaboration XXIV 2016) and the cosmic microwave background (Planck Collaboration XIII 2016). This discrepancy could possibly be resolved by a ~ 40% hydrostatic mass bias (the Planck masses were derived from a Sunyaev-Zeldovich mass proxy calibrated with hydrostatic masses). Despite a significant amount of work focussing on this hydrostatic bias, no agreement has yet been reached regarding its level. Among the most recent studies, we can mention Donahue et al. (2014), von der Linden et al. (2014), Sereno et al. (2015), Hoekstra et al. (2015), Penna-Lima et al. (2017), and Sereno et al. (2017) who found weak-lensing masses larger than the hydrostatic estimates by ~ 20 − 30%. On the other hand, Gruen et al. (2014), Israel et al. (2014), Applegate et al. (2016), or Smith et al. (2016) obtained results consistent with a vanishing bias. Detailed comparisons between lensing and X-ray masses have also highlighted a radial (e.g. Mahdavi et al. 2008; Donahue et al. 2014), mass (e.g. Hoekstra et al. 2015), and redshift dependence (e.g. Sereno et al. 2015; Smith et al. 2016) of the hydrostatic bias.

This variety of results suggests that, in addition to likely selection effects between the different samples (e.g. Sereno et al. 2015), the systematics of the X-ray and lensing techniques are not yet fully controlled. This seems to be particularly true for lensing estimates, as shown with the discrepant results found by different teams for common cluster subsamples (e.g. Okabe & Smith 2016). Therefore, it is necessary to employ other estimators to further constrain a possible hydrostatic or lensing mass bias. An interesting alternative is to use the dynamics of cluster galaxies. For instance, Rines et al. (2016) compared the scaling relation between Planck masses and galaxy velocity dispersions, ruling out the large bias of ~ 40% required by the Planck cosmological results. Maughan et al. (2016) applied the caustic method of Diaferio & Geller (1997) to infer the mass of 16 massive clusters, finding also no evidence for an hydrostatic bias since their X-ray masses are actually ~ 20% larger than the caustic estimates.

Our goal is to explore this route, alternative to the widely used lensing technique, to constrain the hydrostatic mass bias. To do so, we present here the dynamical analysis of ten galaxy clusters among the most luminous in X-rays at intermediate redshifts. To assess possible systematics in the dynamical mass estimates, we used three approaches: the virial theorem, the caustic method, and the Jeans equation. The latter is used here for the first time to derive statistical constraints on the hydrostatic bias; our results based on the caustic masses are also among the very first of their kind. Thanks to a detailed analysis of each cluster, made possible by the combination of wide field spectroscopy and photometry, we put a particular attention to the treatment of substructures.

The paper is organised as follows. In Sect. 2, we introduce the cluster sample and describe the data sets used for this study. We outline the three dynamical mass estimators in Sect. 3 and we describe the substructure analysis in Sect. 4. The comparison of the mass estimators is presented in Sect. 5, before concluding in Sect. 6. We briefly discuss the clusters individually in the Appendix, which also contains some intermediate results of the dynamical and substructure analyses. Our results are scaled to a flat, ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7 and a Hubble constant H₀ = 70 km s^-1 Mpc^-1.

2. Data: description and reduction

2.1. The DXL sample

Thirteen medium distant X-ray galaxy clusters with luminosities ${\mathit{L}}_{\mathrm{X}}^{\mathrm{bol}}\mathrm{=}\mathrm{0.5}\mathrm{-}\mathrm{4}\mathrm{\times }{\mathrm{10}}^{\mathrm{45}}\mathrm{erg}\hspace{0.17em}{\mathrm{s}}^{-1}$ were selected from the ROSAT-ESO Flux Limited X-ray survey (REFLEX, Böhringer et al. 2001, 2004) to form a statistically complete sample (DXL; see e.g. Zhang et al. 2004, for more details). It contains the clusters that are the most X-ray luminous in the redshift interval z = 0.26 − 0.31 and it covers a mass range M₅₀₀ = 0.48 − 1.1 × 10¹⁵ M_⊙ (Zhang et al. 2006). Its volume completeness can be estimated with the well known selection function of the REFLEX survey (Böhringer et al. 2004). The most luminous clusters of the original DXL sample (ROSAT luminosities L_X> 10⁴⁵ erg s^-1 in the (0.1–2.4 keV) band) were selected for a wide-field photometric and spectroscopic follow-ups, allowing for a comprehensive analysis of their properties; the DXL cluster RXC J2011.3-5727 was also included despite its luminosity of L_X = 6.7 × 10⁴⁴ erg s^-1. Three clusters were added to the spectro-photometric campaign to cover a larger redshift range: RXC J1206.2-0848 and RXC J1347.5-114 at redshifts z ~ 0.45, and RXC J0225.9-4154 at redshift z ~ 0.22. They were also selected according to their high X-ray luminosities, the former two being the most luminous clusters at high redshifts and the latter among the four most luminous in the redshift range z ~ 0.2 − 0.22.

The X-ray properties of the DXL clusters are presented by Zhang et al. (2004−2006) and Finoguenov et al. (2005), along with results on the calibration of scaling relations. Braglia et al. (2007, 2009) and Pierini et al. (2008) studied the galaxy content of RXC J0014.3-3022, RXC J0232.2-4420, and RXC J2308.3-0211, focusing on the star formation activity as function of environment and the properties of the diffuse stellar emission around the brightest cluster galaxies. Ziparo et al. (2012) and Foëx et al. (2017) conducted a detailed analysis of the structure and dynamical state of RXC J0225.9-4154, RXC J0528.9-3927, RXC J1131.9-1955, and RXC J2308.3-0211.

2.2. Optical spectroscopy

Multi-object spectroscopy observations were carried out with the VIMOS instrument mounted at the Nasmyth focus B of VLT-UT3 Melipal at Paranal Observatory (ESO), Chile. The VIMOS provides an array of four identical CCDs separated by a 2′ gap, each with a field of view (FOV) of 7 × 8 arcmin² and a 0.205′′ pixel resolution. Each cluster was observed at three different pointings, covering an extended region along the major axis of the cluster shape (as observed in X-rays) and overlapping over its core. Given the size of VIMOS’ total FOV, the observations span a roughly rectangular area of 9 × 5 Mpc² at z = 0.3 (see e.g. Fig. 1 in Braglia et al. 2009). The selection of targets was done only on the basis of their I-band luminosity, in order to avoid any colour bias for the comparative analysis of passive and star-forming galaxies.

The spectra were obtained with the low resolution LR-Blue grism, which provides a spectral coverage from 3700 to 6700 Å, has a spectral resolution of about 200 for 1′′ width slits, and does not suffer from fringing. Moreover, it allows up to four slits in the direction of dispersion, thus significantly increasing the number of targets per mask. For a galaxy at z ~ 0.3, this grism covers important spectral features such as the [ OII ] , [ OIII ] , H_β, H_δ emission lines, the CaII_{H + K} absorption lines and the 4000 Å break, providing reliable redshift estimates up to z ~ 0.8. The reduction of VIMOS spectra was performed in a standard way with the VIPGI software (Scodeggio et al. 2005).

To estimate spectroscopic redshifts (hereafter z_spec), we first ran the EZ tool (Garilli et al. 2010) in blind mode, restricted to z ∈ [ 0 − 2 ]. In the second step, we reviewed the spectra by eye, and used VIPGI for the manual detection and fit of spectral features in case of probable misidentification. Since EZ does not provide redshift errors, we relied on repeated observations of the same object to estimate a typical uncertainty. We found an average value δ_cz ~ 300 km s^-1 with variations of ~ 50 km s^-1 from cluster to cluster; velocity dispersions were corrected accordingly, following the prescription of Danese et al. (1980).

Some of the spectroscopic data have already been used in previous works, for example for RXC J0014 (Braglia et al. 2009), RXC J1131 (Ziparo et al. 2012), or RXC J1206 (Biviano et al. 2013). However, we re-reduced them for consistency with the other clusters. We completed our data set with redshifts available in the literature. For RXC J1347, we also reduced additional VIMOS masks found in the ESO archive, from the programmes 186.A-0798 (PI: P. Rosati) and 090.A-0958 (PI: A. Von der Linden).

2.3. Optical imaging

The photometric follow-up of the clusters was conducted with the Wide Field Imager (WFI, Baade et al. 1999) mounted on the Cassegrin focus of the ESO/MPG 2.2-m telescope at La Silla, Chile. The WFI is a mosaic camera composed of 4 × 2 CCD chips, each made of 2048 × 4096 pixels with an angular resolution of 0.238′′/pixel. The total FOV is 34′ × 33′, which fully encompasses the region observed with VIMOS. The data reduction was performed with the THELI pipeline (Schirmer 2013), which performs the basic pre-processing steps (bias subtraction, flat-fielding, background modelling and sky subtraction), and uses third party softwares for the astrometry (Scamp, Bertin 2006) and the co-addition of mosaic observations (SWarp, Bertin 2010). The photometry was made with SExtractor. Stars, galaxies, and false detections were sorted according to their position in the magnitude-central flux diagram, their size with respect to that of the PSF, and their stellarity index (CLASS_STAR parameter). Luminosities were estimated from the MAG_BEST parameter, whereas colours were computed with MAG_APER, measured in a fixed aperture of 3′′.

The clusters were observed in at least three of the four (B,V,R,I) pass bands (only two for RXC J2011.3-5725), allowing us to compute photometric redshifts (hereafter z_phot). Our approach, which is described in details in Foëx et al. (2017), relies on the simple technique of the “k-nearest neighbour” fitting (Altman 1992), performed in colour space. For each cluster, we divided the sample of z_spec into training and testing sets, and we assessed the accuracy of z_phot with the fraction of catastrophic errors η = | z_phot − z_spec | / (1 + z_spec) > 0.15, and the redshift accuracy σ_z = 1.48 × med [ | z_phot − z_spec | / (1 + z_spec) ] (Ilbert et al. 2006). For most clusters, we obtained η and σ_z within 0.05–0.1 (see top panel of Fig. 2). The variability of these values between clusters is mostly due to the number of available z_spec. A larger number implies a better sampling of the colour space, hence a better accuracy.

2.4. Selection of cluster members

To select the cluster members from the spectroscopic data, we opted for an iterative 3σ clipping scheme, combined with an iterative radial binning in the projected-phase space (hereafter PPS). This method extends the approach originally proposed by Yahil & Vidal (1977), by accounting for radial variations in velocity dispersion; more details of our implementation are given in Foëx et al. (2017). We used the spectroscopic members to compute the clusters’ redshift with the robust biweight estimator of Beers et al. (1990). The number of cluster members and the clusters’ redshifts, along with their 1σ confidence interval obtained from bootstrapping, are given in Table 1. We present in Fig. 1 the location of galaxies in PPS for the cluster RXC J1206 (see Fig. C.3 for the other clusters).

For the galaxies with only a z_phot estimate, we proceeded as follows. First, we discarded the galaxies with | z_phot − μ_c | > max [ 0.1,3σ_c ], where μ_c represents the mean photometric redshift of the spectroscopically confirmed members and σ_c its scatter. To improve the selection, we then removed galaxies having σ_z,spec> 0.1, where σ_z,spec is the dispersion in z_spec of the kNN ten nearest neighbours. These selection criteria lead to a typical completeness of ~ 70 − 80% and purity of ~ 40 − 60% (see bottom panel of Fig. 2), values that were estimated by dividing the sample of z_spec into training and testing sets. For each testing set, the number of spectroscopic members surviving the z_phot selection criteria was used to defined the purity, as the proportion with respect to the total number of galaxies in the final sample, and the completeness, as the proportion with respect to the number of spectroscopic members present in the initial sample.

Finally, we merged for each cluster the spectroscopic and photometric catalogues, giving priority to the spectroscopic classification when possible. From these combined catalogues, we fitted the clusters’ red sequence with a 2σ clipping method in the magnitude-colour diagram. The locus and scatter, σ_RS, of the red sequence were used to divide the catalogues into two broad populations: the red-sequence galaxies, having a colour residual within 3σ_RS, and the blue members. Additionally, the combined catalogues were cut to a limiting magnitude m ≤ m^∗ + 3 in order to reduce the contamination by faint background galaxies.

Assuming that the spectroscopic samples do not contain interlopers, the purity of the combined catalogues are typically increased by ~ 10% as compared to the photometric catalogues. We can also note in Fig. 2 that the completeness and purity of the red-sequence galaxies are ~ 10% better than the values obtained for the full population. If we assume that the contamination by field galaxies is constant over the field of view, we can use the fraction of spectroscopic members in the combined catalogues to estimate the completeness of the spectroscopic survey. It decreases smoothly from ~ 0.5 ± 0.1 at small radii to ~ 0.3 ± 0.1 (~ 0.1 ± 0.1) at R₂₀₀ (1.5 R₂₀₀) for galaxies with magnitudes m_BCG<m<m ∗ + 3. Finally, we can mention the fraction of red-sequence galaxies within R₂₀₀. With average values of ~ 0.4 ± 0.1 for the photometric catalogues, 0.5 ± 0.1 for the combined catalogues, and ~ 0.6 ± 0.1 for the spectroscopic members, we can see that the spectroscopic data are only slightly biased towards the population of early-type galaxies.

3. Dynamical masses

3.1. Jeans equation

Assuming stationarity, the absence of streaming motions, and spherical symmetry, the first velocity moment of the collisionless Boltzmann equation gives the Jeans equation: $\frac{\mathrm{d}\mathit{p}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{d}\mathit{r}}\mathrm{+}\mathrm{2}\frac{\mathit{\beta }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{r}}\mathit{p}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{\nu }\mathrm{\left(}\mathit{r}\mathrm{\right)}\frac{\mathit{GM}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{r}}^{\mathrm{2}}}\mathit{,}$(1)where $\mathit{p}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{\nu }\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ is the radial dynamical pressure, ν(r) the space density of the tracer used to observe the system (i.e. galaxies), σ_r(r) the radial velocity dispersion of the tracer, and M(r) the total mass profile. The tracer’s velocity anisotropy profile, β(r), reads $\mathit{\beta }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\frac{{\mathit{\sigma }}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{·}$(2)It varies from β → − ∞ for circular orbits to β = 1 for radial orbits; β = 0 for an isotropic velocity field.

The Jeans equation is a linear first-order differential equation, which under the boundary condition p(r) → 0 when r → ∞, has the general solution (van der Marel 1994): $\mathit{p}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{r}}^{\mathrm{\infty }}\exp \left[\mathrm{2}{\mathrm{\int }}_{\mathit{r}}^{\mathit{s}}\mathit{\beta }\mathrm{\left(}\mathit{t}\mathrm{\right)}\frac{\mathrm{d}\mathit{t}}{\mathit{t}}\right]\mathit{\nu }\mathrm{\left(}\mathit{s}\mathrm{\right)}\frac{\mathit{GM}\mathrm{\left(}\mathit{s}\mathrm{\right)}}{{\mathit{s}}^{\mathrm{2}}}\mathrm{d}\mathit{s}\mathit{.}$(3)The tracer’s density profile, ν(r), can be inferred from the observed surface density Σ(R) via Abel de-projection, and a third equation links the projected dynamical pressure $\mathit{P}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{\Sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}{\mathit{\sigma }}_{\mathrm{P}}^{\mathrm{2}}\mathrm{\left(}\mathit{R}\mathrm{\right)}$ (where σ_P(R) is the observed line-of-sight velocity dispersion at projected radius R) to the dynamical pressure (Binney & Mamon 1982). However, since the number of unknowns exceeds the number of equations, one needs to make further assumptions to break the well-known degeneracy between the mass and anisotropy profiles (e.g. Merritt 1987; Merrifield & Kent 1990; van der Marel et al. 2000).

Here we have chosen to use a forward fitting approach, following closely the prescription given by Mamon et al. (2013). The method consists in using a maximum likelihood estimator to find the best-fit parameters of analytical M(r), ν(r), and β(r) profiles. In our case, however, we chose to derive the density profile ν(r) prior to the dynamical analysis since our spectroscopic catalogues suffer from a radial-dependent completeness. We used the combined catalogues defined in 2.4 to fit the observed surface density profiles Σ(R) by the sum of a background contribution plus the projection of a spherical NFW (Navarro et al. 1997) or King distribution; the model giving the smallest χ² was chosen (an example is given in Fig. 3; see Fig. C.1 for the full sample). With the characteristic radius r_c of the tracer’s density in hand, the best-fit parameters for M(r) and β(r) are obtained by minimising $\mathrm{-}\ln \mathrm{ℒ}\mathrm{=}\mathrm{-}\begin{array}{c}\sum \end{array}\left[\ln \mathit{g}\mathrm{\left(}{\mathit{R}}_{\mathit{i}}\mathit{,}{\mathit{\upsilon }}_{\mathit{z,i}}\mathrm{\right|}{\theta }\mathrm{)}\mathrm{-}\ln \mathrm{\Sigma }\mathrm{(}{\mathit{R}}_{\mathit{i}}\mathrm{\left|}{\mathit{r}}_{\mathrm{c}}\mathrm{\right)}\right]\mathit{,}$(4)where θ is the vector containing the free parameters of the mass and anisotropy models. Assuming a Gaussian 3D velocity distribution, the galaxy density in PPS is given by (Mamon et al. 2013): $\mathit{g}\mathrm{\left(}\mathit{R,}{\mathit{\upsilon }}_{\mathit{z}}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}{\mathrm{\int }}_{\mathit{R}}^{\mathrm{\infty }}\frac{\mathit{\nu }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathit{z}}\mathrm{\left(}\mathit{R,r}\mathrm{\right)}}\exp \left[\mathrm{-}\frac{{\mathit{\upsilon }}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{z}}^{\mathrm{2}}\mathrm{\left(}\mathit{R,r}\mathrm{\right)}}\right]\frac{\mathit{r}\mathrm{d}\mathit{r}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{2}}}}\mathit{,}$(5)with ${\mathit{\sigma }}_{\mathit{z}}^{\mathrm{2}}\mathrm{\left(}\mathit{r,R}\mathrm{\right)}\mathrm{=}\left[\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{\left(}\mathit{r}\mathrm{\right)}{\left(\frac{\mathit{R}}{\mathit{r}}\right)}^{\mathrm{2}}\right]{\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}$(6)The mass profile enters the above equation via the general solution of the dynamical pressure, ${\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{p}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{/}\mathit{\nu }\mathrm{\left(}\mathit{r}\mathrm{\right)}$.

For simplicity, we restricted our analysis to the combination of a NFW mass profile with three different anisotropy models: a constant anisotropy (hereafter “C” model), the Mamon & Łokas (2005) model (hereafter “ML”), ${\mathit{\beta }}_{\mathrm{ML}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{r}}{\mathit{r}\mathrm{+}{\mathit{r}}_{\mathit{\beta }}}\mathit{,}$(7)and the simplified Tiret et al. (2007) model (hereafter “T”), ${\mathit{\beta }}_{\mathrm{T}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\beta }}_{\mathrm{\infty }}\frac{\mathit{r}}{\mathit{r}\mathrm{+}{\mathit{r}}_{\mathrm{s}}}\mathit{,}$(8)where r_s is the scale radius of the NFW mass profile. Since we do not assume that “light traces mass”, the NFW scale radius can be different than the characteristic radius r_c of the galaxy spatial distribution.

We performed the minimisation of − lnℒ with the Powell algorithm (e.g. Press et al. 1986), and kept the combination (M₂₀₀,c₂₀₀,β) that gives the smallest value. The parameter space was restricted to c₂₀₀ ∈ [ 2,20 ], β_C ∈ [ 0,1 ] (C model), r_β ∈ [ 0,5 ] Mpc (ML model), and β_∞ ∈ [ 0,1 ] (T model). Only galaxies within the radius R₂₀₀ estimated from the virial theorem (see below) were used to compute the likelihood. The 1σ confidence intervals on the best-fit parameters were estimated by bootstrapping the catalogue of galaxies. We accounted for the uncertainties in r_c by sampling its distribution, assumed to be Gaussian, for each bootstrap realisation. In Fig. 4, we compare the projected velocity dispersion profile derived from the Jeans best-fit model to the observed profile for RXC J1206 (see Fig. C.4 for the other clusters).

3.2. Virial theorem

Our second dynamical mass estimator is based on the scalar virial theorem (e.g. Limber & Mathews 1960; Heisler et al. 1985; Merritt 1988): ${\mathit{M}}_{\mathrm{V}}\mathrm{=}\frac{\mathrm{3}\mathit{\pi }}{\mathit{G}}{\mathit{\sigma }}_{\mathrm{P}}^{\mathrm{2}}{\mathit{R}}_{\mathrm{PH}}\mathit{.}$(9)The projected (line-of-sight) velocity dispersion, σ_P, and the projected harmonic mean radius, R_PH, are given by: $\begin{array}{ccc}& & {\mathit{\sigma }}_{\mathrm{P}}^{\mathrm{2}}\mathrm{=}\frac{{\sum }_{\mathit{i}}{\mathit{\upsilon }}_{\mathit{i}}^{\mathrm{2}}}{\mathit{N}\mathrm{-}\mathrm{1}}\mathit{,}\\ & & {\mathit{R}}_{\mathrm{PH}}\mathrm{=}\frac{\mathit{N}\mathrm{(}\mathit{N}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathrm{2}{\sum }_{\mathit{i}\mathit{>}\mathit{j}}{\mathit{R}}_{\mathit{ij}}^{-1}}\mathit{,}\end{array}$with N the number of galaxies, υ_i the rest-frame line-of-sight velocity, and R_ij the angular-diameter distance between galaxy pairs. In practice, we used the robust biweight estimator of Beers et al. (1990) to evaluate σ_P. We used bootstrap realisations of the galaxy catalogue to estimate 1σ confidence intervals on the mass and velocity dispersion.

Similarly to the Jeans analysis, the virial theorem is based on the hypothesis of dynamical equilibrium and sphericity. Interestingly, its application does not require knowledge of the velocity field anisotropy since one always has $\mathrm{3}{\mathit{\sigma }}_{\mathrm{P}}^{\mathrm{2}}\mathrm{=}{\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathit{\theta }}^{\mathrm{2}}$ for a spherical system. However, the virial theorem relies on the additional assumptions that the galaxies have the same spatial and velocity distribution as dark matter particles, and that all galaxies have the same mass. The latter approximation can be justified by the lack of observational evidence of a strong luminosity/mass segregation in the galaxy population (e.g. Adami et al. 1998; Biviano et al. 2002). However, dark matter haloes are well represented by a cuspy NFW profile, whereas cluster members have typically a cored King-like spatial distribution. Numerical simulations have also shown that a velocity bias exists between galaxies and dark matter particles (e.g. Berlind et al. 2003; Biviano et al. 2006; Munari et al. 2013). Despite these various approximations, the virial theorem has proven to be a good estimator (e.g. Biviano et al. 2006), thus it has been widely used owing to its simple and direct implementation.

The theorem, as expressed above, makes the implicit assumption that the system is isolated and observed in its entirety. However, galaxy clusters are embedded in dense environments, continuously accreting matter from their surroundings. Therefore, the virialised region of a cluster cannot be considered as an isolated system. Consequently, an additional correction must be applied to the virial theorem, the so-called surface pressure term (hereafter SPT; see e.g. Binney & Tremaine 1987; Carlberg et al. 1996). Neglecting the dynamical pressure from the radial infall of matter leads to virial masses typically overestimated by ~ 15 − 20% (Carlberg et al. 1997; Girardi et al. 1998; Biviano et al. 2006). When the spectroscopic survey does not fully cover the virialised region, a larger correction should be used, since the velocity dispersion is in general decreasing with radius. Under the hypothesis that the galaxies and dark matter share the same density profile, ρ(r), the corrected virial mass can be expressed as (Girardi et al. 1998): ${\mathit{M}}_{\mathrm{CV}}\mathrm{=}{\mathit{M}}_{\mathrm{V}}\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\mathrm{1}\mathrm{-}\mathrm{4}\mathit{\pi }{\mathit{b}}^{\mathrm{3}}\frac{\mathit{\rho }\mathrm{\left(}\mathit{b}\mathrm{\right)}}{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{b}}\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}}{\left[\frac{{\mathit{\sigma }}_{\mathit{r}}\mathrm{\left(}\mathit{b}\mathrm{\right)}}{\mathit{\sigma }\mathrm{(}\mathit{<}\mathit{b}\mathrm{)}}\right]}^{\mathrm{2}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{,}$(11)where the projected radius b is the maximal extent of the spectroscopic observation (or the aperture radius within which one wishes to estimate the mass), and σ( <b) is the integrated velocity dispersion within b.

For an isotropic velocity field of decreasing dispersion with radius, one has ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{(}\mathit{<}\mathit{r}\mathrm{)}\mathrm{=}\mathrm{3}{\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{(}\mathit{<}\mathit{r}\mathrm{)}\mathrm{\ge }\mathrm{3}{\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$. In that case, the term involving the velocity dispersion is at most equal to 1 / 3 (1 for radial orbits, 0 for circular ones), which is the value that we used to estimate M_CV. For the term involving the density ρ(r), we used the characteristic radius r_c obtained previously from the fit of the galaxy surface density profile by either a NFW or King model. For the clusters analysed in this work, we found SPT corrections in the range 20 − 30% at R₂₀₀.

In Foëx et al. (2017) we outlined a procedure to estimate the virial mass (or M₂₀₀) without prior knowledge of the virial radius (or R₂₀₀). We followed this approach to derive the radii R₂₀₀ used to select the galaxies from which the likelihood of the Jeans analysis was computed. Otherwise, we ran the virial theorem within a given aperture, for example R₂₀₀ from the Jeans analysis or the X-ray R₅₀₀, in order to make meaningful comparisons between the different mass estimators.

3.3. Caustics

In contrast with the two previous mass estimators, the caustic approach (Diaferio & Geller 1997; Diaferio 1999) does not rely on the hypothesis of dynamical equilibrium. It is based on the simple consideration that particles gravitationally bound to a system cannot have a velocity larger than the escape velocity ${\mathit{\upsilon }}_{\mathrm{esc}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{\phi }\mathrm{\left(}\mathit{r}\mathrm{\right)}$. For a realistic mass distribution, the gravitational potential φ(r) is an increasing function of radius. Therefore, the maximum allowed velocity decreases with radius, producing a trumpet-like pattern in PPS that is delimited by the upper and lower caustics. Accounting for projection effects, that is taking the average line-of-sight component of the escape velocity, one can obtain a simple relation that links the caustic amplitude, , to the mass profile (Diaferio & Geller 1997): $\mathit{GM}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}{\mathrm{𝒜}}^{\mathrm{2}}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathrm{ℱ}}_{\mathit{\beta }}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{d}\mathit{s}\mathit{.}$(12)The projection along the line of sight is captured within the filling factor ℱ_β(r): ${\mathrm{ℱ}}_{\mathit{\beta }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{\pi G}\frac{\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{r}}^{\mathrm{2}}}{\mathit{\phi }\mathrm{\left(}\mathit{r}\mathrm{\right)}}\left(\frac{\mathrm{3}\mathrm{-}\mathrm{2}\mathit{\beta }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{\left(}\mathit{r}\mathrm{\right)}}\right)\mathit{,}$(13)where ρ(r) and β(r) are the mass density and anisotropy profiles, respectively.

Similarly to the Jeans analysis, a degeneracy between the mass and anisotropy profiles prevents estimating a cluster mass from spectroscopic observations alone. Therefore, it is customary to assume a constant value for the filling factor. Serra et al. (2011) showed that setting F_β(r) = 0.7 provides a good approximation, allowing for the recovery of true masses within 10% at radii larger than ~ 0.6 R₂₀₀. As shown in Fig. 5 for a NFW mass profile with M₂₀₀ = 10¹⁵ M_⊙ and c₂₀₀ = 4, this approximation is roughly equivalent to assuming a constant anisotropy β ≈ 0.5, or a β_ML(r) profile with r_β ≈ 0.2 R₂₀₀. At small radii, the filling factor takes smaller values, hence caustic masses obtained with this method are typically overestimated. In principle one could use the β(r) profile derived from the Jeans analysis to obtain a more robust mass profile (e.g. Biviano et al. 2013). However, we chose to set a constant ℱ_β(r) = 0.7, to avoid introducing too much correlation between the two mass estimators, and to facilitate direct comparisons with previous works.

Our implementation of the caustic method follows closely the prescription given by Diaferio (1999). The first step consists in estimating the galaxy density in PPS, after an appropriate rescaling of the coordinates (we chose a scaling factor q = σ_υ/σ_r = 25). The smooth density map, f_q(r,υ), was obtained with the DEDICA algorithm (Pisani 1996), which uses an iterative and adaptive kernel estimator. In the second step, we determined the upper and lower caustics as the loci of pairs (r,υ) where f_q(r,υ) = κ. To determine the adequate iso-density level κ that sets the caustic amplitude, we can assume that the virial condition $\mathrm{⟨}{\mathit{\upsilon }}_{\mathrm{esc}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{R}}\mathrm{=}\mathrm{4}\mathrm{⟨}{\mathit{\upsilon }}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{R}}$ holds within the central region, defined here by R = R₂₀₀ from the virial theorem. Since this condition remains valid with line-of-sight velocities under the hypothesis of isotropic orbits, the density level κ is obtained as the root of the equation $\mathrm{⟨}{\mathit{\upsilon }}_{\mathrm{esc}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{\kappa ,R}}\mathrm{-}\mathrm{4}{\mathit{\sigma }}_{\mathrm{P}}^{\mathrm{2}}\mathrm{(}\mathit{<}\mathit{R}\mathrm{)}\mathrm{=}\mathrm{0}$, with $\mathrm{⟨}{\mathit{\upsilon }}_{\mathrm{esc}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{\kappa ,R}}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{R}}{\mathrm{𝒜}}_{\mathit{\kappa }}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{\varphi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\mathit{/}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{R}}\mathit{\varphi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}$ and $\mathit{\varphi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{}^{\mathrm{\int }}\mathit{f}_{\mathit{q}}\mathrm{\left(}\mathit{r,\upsilon }\mathrm{\right)}\mathrm{d}\mathit{\upsilon }$.

At each radius, the caustic amplitude was set to , where are the lower and upper caustics. This approach, as compared to taking their average value, reduces the contamination by interlopers, and limits the impact of massive high-velocity substructures. The last step of the algorithm consists in computing the logarithmic derivative , which should be ≲ 1 / 4 for a typical cluster. Following the prescription by Serra et al. (2011), we replaced with a new value yielding whenever the original derivative was superior than two (e.g. due to interlopers or substructures). Finally, the error on the caustic amplitude was estimated as , where max [ f_q(r,υ) ] is the maximum density along the υ-axis at fixed r (Diaferio 1999). Serra et al. (2011) found that this recipe gives a 50% confidence level on the true amplitude (see their Fig. 16). Therefore, we multiplied it by 1.4 to get a 1σ error before propagating it to the mass profile. Figure 6 presents the caustic mass profile of RXC J1206 (see Fig. C.5 for the other clusters).

4. Substructures

The strongest prerequisite for obtaining unbiased masses with the Jeans analysis and the virial theorem is to have a tracer that reached dynamical equilibrium. High-velocity substructures associated with mergers disturb the cluster’s velocity field and increase its overall velocity dispersion, thus leading to overestimated masses. The caustic amplitude in PPS depends also, though indirectly, on the galaxy velocity dispersion ${\mathit{\sigma }}_{\mathrm{P}}^{\mathrm{2}}\mathrm{(}\mathit{<}\mathit{R}\mathrm{)}$. Thus we can expect the caustic estimator to fail in some cases, for instance when a major merger is taking place (Diaferio 1999). Therefore, we need to quantify the degree of relaxation of the clusters, in order to assess the robustness of their mass estimates. To do so, we introduce in the following two substructure indicators: a photometric value that is based on the galaxy surface density excess with respect to its best-fit model, and a dynamical indicator that quantifies local deviations of the velocity distribution with respect to the overall dynamics. We also outline our procedure for the identification and removal of individual substructures.

4.1. Galaxy surface density

We constructed galaxy surface density maps, Σ(x,y), from the photometric plus spectroscopic catalogues of cluster members. After smoothing with a Gaussian kernel of width 100 kpc, the maps were fitted with a two-dimensional elliptical King profile: $\begin{array}{ccc}& & {\mathrm{\Sigma }}_{\mathrm{K}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{0}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{c}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{+}\mathit{b,}\\ & & {\mathit{r}}^{\mathrm{2}}\mathrm{=}\mathrm{(}\mathit{x}\cos \mathit{\phi }\mathrm{+}\mathit{y}\sin \mathit{\phi }{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{y}\cos \mathit{\phi }\mathrm{-}\mathit{x}\sin \mathit{\phi }{\mathrm{)}}^{\mathrm{2}}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{e}{\mathrm{)}}^{\mathrm{2}}\mathit{.}\end{array}$The ellipticity and position angle (e,φ) describe the shape of the cluster. The centre was fixed at the highest density peak, corresponding to the origin of the coordinates. r_c is the core radius, Σ₀ the central surface density, and b ≥ 0 a constant background contribution. Defining the residual δ_i,j = Σ(x_i,y_j) − Σ_K(x_i,y_j) between the measured surface density and the best-fit model, we computed the photometric substructure indicator $\mathrm{\Delta }\mathrm{=}\frac{{\sum }_{\mathit{i,j}}\max \mathrm{\left[}\mathrm{0}\mathit{,}{\mathit{\delta }}_{\mathit{i,j}}\mathrm{\right]}}{{\sum }_{\mathit{i,j}}\mathrm{\Sigma }\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathit{,}{\mathit{y}}_{\mathit{j}}\mathrm{\right)}}\mathrm{·}$(15)The sum was limited to pixels within R₂₀₀ and we only considered positive excess in order to focus on deviations produced by the largest substructures. Uncertainties on Δ were obtained by propagating the Poisson noise of the surface density maps. Assuming that the number of galaxies in a substructure scales linearly with its mass, as does the richness of a cluster with its total mass, then Δ provides a crude estimate of the mass fraction contained in substructures.

Weißmann et al. (2013) analysed the X-ray morphology of 80 galaxy clusters with the classical power ratio P3 /P0_max and centre shifts ω_c. They also introduced the maximum power ratio P3 /P0_max, corresponding to the peak of the ratio profile within 0.3 − 1 R₅₀₀. Seven of our clusters are in their sample, therefore we can compare our photometric substructure indicator with their X-ray morphology estimators. The results are presented in Fig. B.2. The Pearson coefficients ρ(Δ,P3 /P0_max) = 0.89 and ρ(Δ,ω_c) = 0.83 indicate a good correlation between the X-ray and optical morphologies.

4.2. Dressler-Shectman test

Our second indicator is based on the Dressler & Shectman test (hereafter DS test, Dressler & Shectman 1988), which combines spatial and line-of-sight velocity information. For each of the N galaxies in the spectroscopic sample, the local velocity, ⟨ υ ⟩ _loc, and projected velocity dispersion, σ_loc, are estimated with the n_NN nearest neighbours (in projected distance). We used ${\mathit{n}}_{\mathrm{NN}}\mathrm{=}\sqrt{\mathit{N}}$, since it provides a better sensitivity to significant substructures while being less affected by Poisson noise (e.g. Silverman 1986). The dynamical deviation of the ith galaxy is quantified as:${\mathit{\delta }}_{\mathit{i}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{NN}}\mathrm{+}\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{P}}^{\mathrm{2}}}\mathrm{\left[}\mathrm{\right(}\mathrm{⟨}\mathit{\upsilon }{\mathrm{⟩}}_{\mathrm{loc}\mathit{,i}}\mathrm{-}\mathrm{⟨}\mathit{\upsilon }\mathrm{⟩}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{\sigma }}_{\mathrm{loc}\mathit{,i}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{P}}{\mathrm{)}}^{\mathrm{2}}\mathrm{]}\mathit{,}$(16)where ⟨ υ ⟩ and σ_P are the global values, obtained with the N galaxies. As pointed out by Pinkney et al. (1996), gradients in the velocity dispersion profile can produce false positive detections of substructures. Therefore, we used a radial-dependent σ_P (best-fit power law of the observed profile) as the “global” value against which σ_loc is compared.

For each cluster, we ran the DS test on 10⁴ random realisations of the galaxy distribution, where positions were kept fixed but velocities shuffled, in order to erase any correlation between them. We used the resulting probability distribution of the deviation δ to specify the threshold δ_min such as P(δ>δ_min) = 0.05. Above this value, the local dynamics of a galaxy differs significantly from that of the cluster, thus it is likely to be part of a substructure. Finally, we defined our dynamical substructure indicator, f_DS, to be the fraction of cluster members satisfying this selection criterion, that is f_DS = N(δ_i>δ_min) /N. As for the photometric indicator, we limited the analysis within R₂₀₀; uncertainties on f_DS were obtained assuming Poisson statistics.

The values of the two substructure indicators are given in Table 3 (see also Fig. B.1). A Pearson coefficient of ρ(Δ,f_DS) = 0.69 indicates a moderate correlation. Since both indicators are based on relative galaxy counts, we used their average value as a third indicator, hereafter Δf. All the clusters show a rather high degree of substructure with Δ > 0.15. RXC J0014 presents the largest value for both indicators, while RXC J1347 appears to be the most regular cluster in the sample.

4.3. Identification and removal of substructures

To identify substructures, we first made use of the DS test results to locate groups of galaxies having δ>δ_min. Following Foëx et al. (2017), we ran the DS test using separately the deviations in the local mean velocity and those in the local velocity dispersion. For each group, we computed its rest-frame velocity, velocity dispersion, average position, and projected shape (using the moment approach, e.g. Carter & Metcalfe 1980). These values were used as initial guess for the 3D version of the Kaye’s mixture model algorithm (KMM, Ashman et al. 1994), which is a typical iterative expectation-maximisation algorithm for the modelling of a mixture of Gaussian distributions. The KMM algorithm was developed initially to separate different components in velocity space. However, it is straightforward to include spatial information by using multivariate distributions (more details are given in Foëx et al. 2017). Starting from the initial guess describing the parameters of each component, the KMM algorithm estimates the probability that a given galaxy belongs to a given component. It then partitions the galaxies, and re-estimates the parameters of each component before the next iteration. Once the algorithm has converged, it provides the list of galaxies associated with the cluster main body and with each additional substructure.

We also used the galaxy surface density maps to identify substructures not detected with the DS test. Here we focused on the spatial distribution of red-sequence galaxies, in particular the bright ones with m<m^∗ + 1, since the z_phot accuracy, completeness, and purity are higher for this galaxy population (see Fig. 2). We located the largest over-density peaks, and defined the substructures as elliptical regions around these peaks.

Finally, we combined the two previous steps to select the galaxies that most likely trace the relaxed component of a cluster. We kept the main KMM partition and discarded the galaxies located within the additional substructures defined from the galaxy surface density maps. An example is provided in Fig. 7 (see Figs. C.6–C.8 for the other clusters). It shows the isopleths of the galaxy surface density (combined catalogues), the spatial distribution of the spectroscopically confirmed cluster members, highlighting those associated with a KMM partition, and the ellipses englobing the extra substructure candidates.

To summarise, we excluded the galaxies which are part of substructures identified due to their specific dynamics or as over-densities in the galaxy surface density. The remaining spectroscopic members were then used to estimate the dynamical mass of the relaxed component of the clusters (see Sect. 5.5).

5. Comparison of the different mass estimators

The three dynamical estimators make use of the same data set, however, they rely on different hypothesis and simplifications. Therefore, we start by looking for differences in their results. We constrain the scaling relation between mass and velocity dispersions, then we examine how the dynamical estimators compare with the X-ray hydrostatic masses. Finally, we investigate how the dynamical measurements are affected by using only the red-sequence galaxies, or by excluding galaxies located in substructures. The results are summarised in Table 6.

5.1. Dynamical masses

As mentioned previously, we used the radius R₂₀₀ derived from the Jeans analysis as the aperture within which we applied the virial theorem, giving the mass M_V(R₂₀₀). We did not fit the caustic mass profile with an NFW model, or directly estimate a non-parametric M₂₀₀. Instead, we simply considered the caustic mass at this radius, labelled M_c(R₂₀₀) hereafter. While this approach introduces a correlation between the different estimators, it has the advantage of comparing masses within the same physical radius. For each pair of estimators, we estimated the arithmetic mean ratio, its standard deviation, and the geometric mean ratio, which is equivalent to the best-fit intercept of a linear regression of slope 1 in log-log space. The latter has been used for instance by Smith et al. (2016) to infer the hydrostatic mass bias with respect to lensing masses. We also estimated the logarithmic orthogonal scatter, σ_⊥, of the points around this best-fit linear regression. Errors on each parameter were obtained from bootstrap realisations of the cluster catalogue. Given the limited number of objects and the small mass range, we did not try to fit a slope for the regression in log-log space. The individual cluster masses are given in Table 2 and the relation between the different estimators are summarised in Table 6.

The top-panel of Fig. 8 presents the comparison between the Jeans and caustic estimators. The average mass ratio is compatible with 1, indicating the absence of systematics between the two estimators. This suggests that having used a constant ℱ_β = 0.7 was, on average, a valid approximation.

The Jeans and virial estimators (bottom panel of Fig. 8) have a mass ratio larger than 1 at the 3σ level, the latter giving masses larger by ~ 15%. One possibility to explain this bias is the contamination by interlopers. Working with N-body numerical simulations, Biviano et al. (2006) have shown that the projected harmonic radius, R_PH, tends to be overestimated because of to this contamination. They found that the virial estimator gives masses that are typically overestimated by ~ 10%, which is close to the value found here. Another possibility could be an underestimation of the SPT. As mentioned previously, we made the hypothesis of isotropy to correct the M_V masses. The results of the Jeans analysis suggest that the clusters’ velocity field has some degree of anisotropy, so we can assume that the SPTs were indeed underestimated, leading to overestimated M_V masses.

5.2. Scaling relation M − σ_P

For a virialised cluster, the gravitational potential energy scales with its kinetic energy, thus one has $\mathit{M}\mathit{/}\mathit{R}\mathrm{\propto }{\mathit{\sigma }}_{\mathit{\upsilon }}^{\mathrm{2}}$, where σ_υ is the 3d velocity dispersion of the dark matter particles. Using the definition of a cluster mass, ${\mathit{M}}_{\mathrm{\Delta }}\mathrm{=}\mathrm{4}\mathit{/}\mathrm{3}\mathit{\pi }{\mathit{R}}_{\mathrm{\Delta }}^{\mathrm{3}}\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{c}}$, one obtain the simple scaling relation $\mathit{h}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{M}}_{\mathrm{\Delta }}\mathrm{\propto }{\mathit{\sigma }}_{\mathit{\upsilon }}^{\mathrm{3}}$, where the reduced Hubble “constant” h(z) = H(z) /H₀ comes from the definition of the critical density ρ_c(z) ∝ h²(z). Since this scaling law has been validated at high accuracy with numerical simulations (e.g. Evrard et al. 2008), it offers an interesting consistency check for our dynamical measurements. Figure 9 (top panel) presents the results obtained for the Jeans M₂₀₀ and the velocity dispersion estimated within the corresponding R₂₀₀. Here we used the BCES(X|Y) estimator of Akritas & Bershady (1996) to fit the regression log [ h(z)M₂₀₀ ] = Alog σ_P + B, treating the velocity dispersion as the response variable. With a slope A = 2.99 ± 0.41, the agreement with the theoretical prediction is excellent. It is also interesting to see that the best-fit normalisation matches very well that of Biviano et al. (2006) after a rescaling to a redshift z = 0.3 and setting ${\mathit{\sigma }}_{\mathit{\upsilon }}\mathrm{=}\sqrt{\mathrm{3}}{\mathit{\sigma }}_{\mathrm{P}}$ (their normalisation was obtained assuming a constant slope of three). The scatter in log [ h(z)M₂₀₀ ] at fixed velocity dispersion is σ_log M = 0.07 ± 0.02. The calibration of the scaling with caustic masses gives similar results, see bottom panel of Fig. 9. Despite a shallower slope A = 2.63 ± 0.35, the best-fit regression agrees, within its errors, with the relation of Biviano et al. (2006) over the mass range probed by our sample. Here the scatter in mass is σ_log M = 0.06 ± 0.01.

We also calibrated the σ_P(M₂₀₀) scaling relation with the Jeans masses. We obtained a slope A = 0.303 ± 0.041 and a normalisation 10^B = 1042 ± 18 km s^-1 for masses in units of 10¹⁵h^-1 M_⊙. While the former agrees well with the results of Evrard et al. (2008), our normalisation is slightly smaller than their value 10^B = 1082.9 ± 4.0 km s^-1. Since the scaling relations of Biviano et al. (2006) or Evrard et al. (2008) were obtained with dark matter particles, our results tend to exclude a strong dynamical segregation between galaxies and dark matter.

5.3. X-rays vs. dynamics

We now turn to the comparison between our dynamical masses and the hydrostatic estimates derived from XMM-Newton observations. The latter were taken from Zhang et al. (2006), from Foëx et al. (2012) for RXC J1206 and RXC J1347, and from Mantz et al. (2010a) for RXC J2011. The comparison is done within the same aperture, the X-ray NFW radius R₅₀₀. We interpolated the Jeans best-fit NFW model to R₅₀₀, while we ran the virial theorem using only the galaxies within this aperture (a new SPT was estimated accordingly); the caustic mass profile provides directly M_c(R₅₀₀). The results are presented in Fig. 10 for the Jeans and caustic estimators.

For the Jeans, virial theorem, and caustic estimators, we obtained a mass ratio ⟨ M_J/M_HE ⟩ = 1.22 ± 0.18, ⟨ M_V/M_HE ⟩ = 1.51 ± 0.26, and ⟨ M_c/M_HE ⟩ = 1.32 ± 0.18, respectively. The dynamical masses are on average larger than the hydrostatic values, however, the results do not provide a significant evidence for a bias between the hydrostatic and dynamical estimators (< 2σ). The clusters present a rather large scatter around the mean ratio, σ ~ 50% for the Jeans an caustic estimators, σ ~ 75% for the virial theorem. We can also note that the geometric means of the mass ratios are slightly smaller than their arithmetic values, hence the intercepts of the regression lines in log-log space are closer to a vanishing fractional bias (≲1σ).

The standard deviation of the mass ratios and the residual scatters σ_⊥, are ~ 2 − 3 times larger than the values obtained when comparing the dynamical estimators to each other. This is not surprising, since the latter involve similar assumptions and make use of the same data set, whereas the hydrostatic estimator is based on different physical mechanisms and hypothesis. In particular, one can expect the dynamical estimators to be more affected by substructures in the galaxy distribution. As a matter of fact, we do obtain strong correlations between the dynamical-to-hydrostatic mass ratio and the value of the substructure indicators (see Fig. 11). For the Jeans estimator, the Pearson coefficients are ρ(M_J/M_HE,Δ) = 0.84 and ρ(M_J/M_HE,f_DS) = 0.90. Defining the average value Δf = (f_DS + Δ) / 2, we obtain a correlation ρ(M_J/M_HE,Δf) = 0.95. Using the caustic masses (bottom panel of Fig. 11) or those from the virial theorem leads to very similar correlation coefficients (see Table 4). These results suggest that the dynamical state of a cluster, as traced by its substructure content, is a potential source of bias for dynamical mass measurements.

Finally, we can note that the large scatter of the mass ratio is mainly driven by two outliers, RXC J0014 and RXC J1347, whose dynamical masses are respectively ~ 2 times larger and smaller than their hydrostatic counterpart. While this is not a surprise for RXC J0014 given its level of substructures, the mass ratio of RXC J1347 is more intriguing, since its dynamical and morphological properties are rather simple (more details given in the appendix). Therefore, we can suppose that the X-ray analysis is failing at estimating the correct mass of this particular cluster. If we exclude it from the sample, we end with average mass ratios larger by ~ 10%: ⟨ M_J/M_HE ⟩ = 1.32 ± 0.17, ⟨ M_c/M_HE ⟩ = 1.41 ± 0.17, and ⟨ M_V/M_HE ⟩ = 1.65 ± 0.25; the rms standard deviations are reduced by ~ 10%. In that case, we obtain a fractional mass bias that is significant at the ≳ 2σ level (the bias has a similar significance when using the geometrical means).

5.4. Red-sequence galaxies

When making a mask for MOS observations, it is tempting to focus on red-sequence galaxies, since they provide a higher success rate for selecting cluster members. From the dynamical point of view, one can also argue that a high-velocity galaxy belonging to the red sequence has a smaller risk of being an interloper than a blue galaxy observed with the same rest-frame velocity. In other words, one can expect a smaller contamination by interlopers for this galaxy population, hence more accurate dynamical mass estimates (e.g. Biviano et al. 2006; Saro et al. 2013). It has also been argued that early-type galaxies are a better tracer of the relaxed component of a cluster, whereas the late-type population has not yet fully reached equilibrium because of galaxies on radial orbits falling into the cluster for the first time (e.g. Biviano et al. 1992; Colless & Dunn 1996; Adami et al. 1998; Biviano & Katgert 2004; Barsanti et al. 2016; however see e.g. Rines et al. 2005, 2013; Girardi et al. 2015, for opposite claims). On the other hand, massive elliptical galaxies are more likely to be affected by dynamical friction (e.g. Merritt 1985). Since this mechanism is not taken into account during the Jeans analysis or with the virial theorem, it may be a possible source of bias when estimating masses with the red-sequence galaxies only.

All these considerations motivated us to investigate the impact of limiting the analysis to these galaxies. We first compared the velocity dispersions of the two populations (see Fig. 12). If late-type galaxies are gravitationally bound to the cluster but not yet virialised, their velocity dispersion should be ≲$\sqrt{\mathrm{2}}$ larger than the value obtained for the red-sequence galaxies. With a mean ratio ⟨ σ_blue/σ_RS ⟩ = 1.13 ± 0.05 (standard deviation σ = 0.15 ± 0.04), we confirm that blue galaxies have a larger velocity dispersion. It is interesting to note that RXC J0014 is again a clear outlier (top-right point in Fig. 12). Its red-sequence galaxies have a velocity dispersion significantly larger than its blue population, which is a somewhat counter-intuitive result. A possible explanation is the presence of high-velocity substructures with a galaxy content dominated by red galaxies. In such a configuration, the bulk velocity of the merging sub-haloes outweighs the increased velocity dispersion of the blue population due to infalling galaxies.

Regarding the dynamical masses obtained with the red-sequence galaxies only, M^RS, we find a global offset with respect to the values derived previously when using both the red and blue populations (note that the new masses were estimated within the same radii R₂₀₀, given in the fourth column of Table 2). For the Jeans estimator, we find an average mass ratio $\mathrm{⟨}{\mathit{M}}_{\mathrm{J}}^{\mathrm{RS}}\mathrm{\left(}{\mathit{R}}_{\mathrm{200}}\mathrm{\right)}\mathit{/}{\mathit{M}}_{\mathrm{J}\mathit{,}\mathrm{200}}\mathrm{⟩}\mathrm{=}\mathrm{0.92}\mathrm{\pm }\mathrm{0.02}$ and small standard deviation σ = 0.08 ± 0.02 (see top panel of Fig. 13). We would like to point out that we re-estimated the characteristic radius of the galaxy density profile, to account for the specific spatial distribution of the red galaxies. For the virial theorem, the ratio is smaller: $\mathrm{⟨}{\mathit{M}}_{\mathrm{V}}^{\mathrm{RS}}\mathrm{\left(}{\mathit{R}}_{\mathrm{200}}\mathrm{\right)}\mathit{/}{\mathit{M}}_{\mathrm{V}}\mathrm{\left(}{\mathit{R}}_{\mathrm{200}}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathrm{0.86}\mathrm{\pm }\mathrm{0.04}$, with a similar scatter. Interestingly, using the red-sequence galaxies leads to a better agreement between the two estimators, $\mathrm{⟨}{\mathit{M}}_{\mathrm{V}}^{\mathrm{RS}}\mathit{/}{\mathit{M}}_{\mathrm{J}}^{\mathrm{RS}}\mathrm{⟩}\mathrm{=}\mathrm{1.06}\mathrm{\pm }\mathrm{0.05}$, whereas we had a ratio greater than 1 at the 3σ level. This result seems to confirm the comments made previously: a smaller contamination by interlopers leading to less underestimated harmonic radii and to a more adequate estimate of the SPT under the assumption of isotropic orbits.

With a mass ratio $\mathrm{⟨}{\mathit{M}}_{\mathrm{c}}^{\mathrm{RS}}\mathrm{\left(}{\mathit{R}}_{\mathrm{200}}\mathrm{\right)}\mathit{/}{\mathit{M}}_{\mathrm{c}}\mathrm{\left(}{\mathit{R}}_{\mathrm{200}}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathrm{0.85}\mathrm{\pm }\mathrm{0.03}$, the caustic estimator presents the largest dependence on the galaxy population (see bottom panel of Fig. 13). This result might appear surprising since the caustic estimator is based on the escape velocity, which does not depend a priori on the galaxy population. Moreover, we used the same assumption regarding the filling factor ℱ_β = 0.7. While this approximation is justified for anisotropies β ~ 0.5, we saw in Fig. 5 that one should use a smaller value for an isotropic velocity field. Consequently, we should obtain even smaller caustic masses. This dependence on the galaxy population comes most likely from the combination of different effects. First of all, we remind here that the caustic amplitude is set by the average velocity dispersion within a certain radius. As shown above, the red-sequence galaxies have a smaller velocity dispersion; therefore, the density level in PPS satisfying the virial condition is smaller, leading to smaller masses. A second possibility could be the mislocating of the caustics due to a lack of galaxies observed at the escape velocity, which is better traced by the blue galaxies on their first orbit, prior to virialisation. Conversely, the presence of high-velocity interlopers in the blue population could also have biased high the caustic masses. We can mention here the work led by Gifford et al. (2013), who analysed the impact of a colour selection on caustic mass estimates with semi-analytical simulations. They found a mass bias ranging from ~− 5% for purely red-sequence galaxies to ~+ 10% when decreasing the fraction of red members to the average value of our sample ~ 0.6. These values are consistent with the 15% decrease in mass observed here.

To summarise, the red-sequence galaxies are characterised by smaller velocity dispersions. Using only these galaxies leads to dynamical masses that are ~ 10 − 15% smaller. As a consequence, the average mass ratios between the dynamical and hydrostatic estimators are slightly smaller than the values obtained previously, but their scatter are essentially the same (see fourth part of Table 6). A precise characterisation of the dynamical properties of the two broad populations of red and blue galaxies is beyond the scope of this paper. However, the present results seem to confirm the existence of a dynamical segregation that can potentially affect mass estimates.

5.5. Impact of substructures

The correlation between the dynamical-to-hydrostatic mass ratio and the substructure indicators suggests that the dynamical estimates might be significantly biased by the presence of substructures. To further investigate this possibility, we now compare the masses obtained before and after excluding the galaxies that are associated with substructures (see Sect. 4.3). The results are presented in Fig. 14 and the new masses are listed in Table 5; the average ratios are given in the fifth part of Table 6. For the three estimators, we obtain ⟨ M^MB/M ⟩ < 1 at the ≳ 3σ level, where M^MB stands for the mass estimated after removing substructures. Their impact appears to be larger for the virial estimator, whose masses are ~ 20% smaller.

We find that the mass ratios exhibit moderate anti-correlations with the substructure indicators: $\mathit{\rho }\mathrm{(}{\mathit{M}}_{\mathrm{J}}^{\mathrm{MB}}\mathit{/}{\mathit{M}}_{\mathrm{J}}\mathit{,}\mathrm{\Delta }\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{0.56}$ and $\mathit{\rho }\mathrm{(}{\mathit{M}}_{\mathrm{J}}^{\mathrm{MB}}\mathit{/}{\mathit{M}}_{\mathrm{J}}\mathit{,}{\mathit{f}}_{\mathrm{DS}}\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{0.61}$ for the Jeans estimator, $\mathit{\rho }\mathrm{(}{\mathit{M}}_{\mathrm{c}}^{\mathrm{MB}}\mathit{/}{\mathit{M}}_{\mathrm{c}}\mathit{,}\mathrm{\Delta }\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{0.24}$ and $\mathit{\rho }\mathrm{(}{\mathit{M}}_{\mathrm{c}}^{\mathrm{MB}}\mathit{/}{\mathit{M}}_{\mathrm{c}}\mathit{,}{\mathit{f}}_{\mathrm{DS}}\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{0.68}$ for the caustic masses, and $\mathit{\rho }\mathrm{(}{\mathit{M}}_{\mathrm{V}}^{\mathrm{MB}}\mathit{/}{\mathit{M}}_{\mathrm{V}}\mathit{,}\mathrm{\Delta }\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{0.66}$ and $\mathit{\rho }\mathrm{(}{\mathit{M}}_{\mathrm{V}}^{\mathrm{MB}}\mathit{/}{\mathit{M}}_{\mathrm{V}}\mathit{,}{\mathit{f}}_{\mathrm{DS}}\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{0.42}$ for the virial theorem (the corresponding p-values are given in Table 4). The virial theorem appears to be more affected by the spatial distribution of galaxies, the harmonic radius depending on the inverse distance of galaxy pairs, whereas the Jeans and caustic estimators show a larger dependence on substructures in velocity space. These results indicate that a larger fraction of galaxies in substructures implies a larger impact on dynamical mass estimates. Interestingly, we also find a mass dependence of the relative substructure content. Using the Jeans masses as reference, we obtain correlation coefficients ρ(M,Δ) = 0.77, ρ(M,f_DS) = 0.71, and ρ(M,Δf) = 0.79. These values show that the more massive a cluster, the higher its level of substructure. This agrees well with the hierarchical model of cluster formation, according to which the largest objects are formed later by accretion and merging of smaller-scale systems, leaving them less time to reach a relaxed and homogenous sate.

To summarise, excluding galaxies that belong to substructures leads to dynamical masses ~ 15% smaller. The impact of substructures increases with a cluster mass, as a consequence of two correlations: the larger the mass, the higher the level of substructure, and the higher this level, the larger the effect on the dynamical mass estimators. Therefore, a proper characterisation of substructures is mostly required when dealing with very massive clusters.

Having dealt with substructures, we can also check how the new dynamical masses compare with the hydrostatic values. The results are presented in Fig. 15 for the Jeans and caustic masses (see also the last part of Table 6). The main difference with the results obtained in Sect. 5.3 concerns the average mass ratios, which now agree within 1σ with a vanishing fractional bias. We can also note that the standard deviations are decreased by roughly a factor of two, suggesting that our approach for treating substructures reduces efficiently their impact on dynamical mass estimates.

We have shown in Sect. 5.3 that, prior to the substructure analysis, excluding RXC J1347 from the sample leads to a mass bias that is significant at the ≳ 2σ level. Here, we find $\mathrm{⟨}{\mathit{M}}_{\mathrm{J}}^{\mathrm{MB}}\mathit{/}{\mathit{M}}_{\mathrm{HE}}\mathrm{⟩}\mathrm{=}\mathrm{1.06}\mathrm{\pm }\mathrm{0.06}$ (σ = 0.18), $\mathrm{⟨}{\mathit{M}}_{\mathrm{c}}^{\mathrm{MB}}\mathit{/}{\mathit{M}}_{\mathrm{HE}}\mathrm{⟩}\mathrm{=}\mathrm{1.07}\mathrm{\pm }\mathrm{0.09}$ (σ = 0.25), and $\mathrm{⟨}{\mathit{M}}_{\mathrm{V}}^{\mathrm{MB}}\mathit{/}{\mathit{M}}_{\mathrm{HE}}\mathrm{⟩}\mathrm{=}\mathrm{1.17}\mathrm{\pm }\mathrm{0.11}$ (σ = 0.30) for the Jeans, caustic, and virial estimator, respectively. Thus, the standard deviations are decreased roughly by a factor of two while the average ratios increase by ~+ 0.07. However, they stay compatible with no bias at the ~ 1σ level (1.5 for the virial theorem estimator). In other words, excluding RXC J1347 does not change our main conclusion: after the substructure analysis, we do not observe any significant bias between the hydrostatic and dynamical mass estimators.

We found in Sect. 5.4 that the blue population tends to have a larger velocity dispersion than the red-sequence galaxies. Girardi et al. (2015) suggested that the dynamical state of a cluster could be responsible for this segregation. This is clearly the case for RXC J0014, but with the opposite result, that is a larger velocity dispersion for the red-sequence galaxies. Now that we have isolated the galaxies belonging to the clusters’ main body, we can check whether the two galaxy populations still present the same difference regarding their velocity dispersion. With an average ratio ⟨ σ_blue/σ_RS ⟩ = 1.27 ± 0.05, we find that the segregation is even more pronounced; the corresponding standard deviation is σ = 0.16 ± 0.03, which is very similar to the value quoted in Sect. 5.4. Since we do not find a significant correlation between the ratio of velocity dispersions and either substructure indicators (| ρ | < 0.4), our results indicate that the two populations have a different velocity distribution regardless of the dynamical state of their host. However, it is possible that our findings are specific to the cluster population analysed here (massive objects at moderate redshifts). A comparison with clusters at lower redshifts would provide more insights on a possible evolutionary trend of the dynamical segregation between early- and late-type galaxies (see e.g. Barsanti et al. 2016).

6. Conclusions

In this paper we presented the dynamical analysis of ten X-ray luminous galaxy clusters at redshifts z = 0.2 − 0.45 based on wide-field VIMOS spectroscopic data. Additional WFI photometric observations were used to compute photometric redshifts, in order to constrain the clusters’ morphology, to fit their galaxy surface density profiles, and to locate substructures. We derived dynamical masses with three different approaches: the classical virial theorem estimator, a fully parametric forward fitting procedure to solve the Jeans equation, and the caustic method.

We introduced two substructure indicators to quantify the clusters’ degree of relaxation. The first indicator is based on the galaxy surface density excess with respect to the best-fit 2D elliptical King model. The second indicator makes use of the Dressler and Shectman test to estimate the fraction of galaxies associated with substructures having a specific dynamics. We also described a methodology to properly identify the galaxies that belong to the main body of the clusters. It combines the visual identification of substructures in the galaxy surface density map with a minimisation-expectation algorithm to partition the galaxies into the different components found with the DS test.

After comparing the results of the three estimators, we constrained the well-known scaling relation between mass and velocity dispersion. Then we confronted our dynamical mass estimates with the hydrostatic values derived from the X-ray data. Finally, we investigated the effect of a colour-based selection of cluster members and the impact of substructures on the dynamical measurements. Our results can be summarised as follows.

• The Jeans and caustic estimators give very similar results, ⟨ M_c/M_J ⟩ = 1.04 ± 0.06, with a standard deviation σ = 0.20 ± 0.05. Masses from the virial theorem are ~ 15% larger, probably because of underestimated SPTs and/or overestimated harmonic radii.

• Our calibration of the scaling relation between mass and velocity dispersion agrees perfectly with the theoretical prediction and the results of numerical simulations for dark matter particles. In particular, the relation calibrated with the Jeans masses is indistinguishable from that of Biviano et al. (2006); the caustic masses lead to a slightly smaller logarithmic slope, but the agreement is nonetheless very good over the mass range probed by our sample. Such results indicate the absence of a strong dynamical segregation between galaxies and dark matter.

• On average, the dynamical masses are larger than the hydrostatic values: ⟨ M_J/M_HE ⟩ = 1.22 ± 0.18, ⟨ M_c/M_HE ⟩ = 1.32 ± 0.18, and ⟨ M_V/M_HE ⟩ = 1.51 ± 0.26 for the Jeans, caustic, and virial theorem estimators, respectively.

• The large scatters of the ⟨ M_dyn/M_HE ⟩ ratios (~ 50% for the Jeans and caustic estimators, ~ 75% for the virial theorem) are mainly driven by two clusters, RXC J0014 and RXC J1347. The former clearly has an overestimated dynamical mass due to substructures (highest value of our indicators), whereas the latter most likely has an overestimated hydrostatic measurement. Excluding RXC J1347 from the sample leads to a positive fractional mass bias that is significant at the ≳ 2σ level.

• We found strong correlations between the ⟨ M_dyn/M_HE ⟩ ratios and the substructure indicators, in particular with f_DS (substructures in velocity space). This highlights the limitations of the dynamical mass estimators for clusters observed in a highly disturbed state.

• With an average ratio ⟨ σ_blue/σ_RS ⟩ = 1.13 ± 0.05 (scatter of ~ 15%), we confirmed previous claims that the blue star-forming galaxies are characterised by a larger velocity dispersion than the passive red-sequence galaxies. The trend is stronger when excluding the galaxies belonging to substructures: ⟨ σ_blue/σ_RS ⟩ = 1.27 ± 0.05. Moreover, the absence of a correlation between this ratio and the substructure indicators suggests an intrinsic dynamical segregation between the two galaxy populations regardless of the dynamical state of their cluster host.

• The dynamical masses obtained with the red-sequence galaxies are 10 − 15% smaller than the value derived with the whole population. The impact of the colour selection is the largest for the caustic estimator, most likely due to a lack of red galaxies observed at the escape velocity.

• The removal of substructures leads to dynamical masses that are ~ 15% percent smaller; the largest impact was found for RXC J0014, whose mass is reduced by a factor of ~ 1.5. We found that more massive clusters have a larger fraction of substructures and that the more substructures, the larger the impact on dynamical mass estimates. These correlations imply that a proper treatment of substructures is mainly required for the most massive galaxy clusters.

• The dynamical masses perfectly agree with the hydrostatic values after the removal of substructures. Namely, we obtained average ratios ⟨ M_J/M_HE ⟩ = 0.98 ± 0.09, ⟨ M_c/M_HE ⟩ = 1.00 ± 0.10, and ⟨ M_V/M_HE ⟩ = 1.08 ± 0.13 for the Jeans, caustic, and virial theorem estimators, respectively. The standard deviations around the means are roughly divided by two, indicating that our procedure properly corrects the impact of substructure on dynamical masses. Excluding RXC J1347 from the sample no longer gives a significant fractional mass bias between the hydrostatic and dynamical estimators.

The takeaway message of this work is that substructures can significantly affect dynamical mass measurements of massive clusters, leading to a spurious mass bias when comparing to the hydrostatic estimates. However, a careful analysis based on large spectroscopic data sets allows to efficiently reduce their impact. Therefore, dynamical mass estimators are an interesting alternative to lensing measurements for the mass calibration of clusters’ scaling relations. Our work was focused on the high-mass end of the cluster population, with results that seem to rule out the hydrostatic mass bias commonly found with lensing studies. Conducting a similar analysis on a larger sample including lower-mass systems would allow for a better characterisation of the possible systematics of the X-ray, dynamical, and lensing mass estimators.

Acknowledgments

The authors thank Mischa Schirmer for precious advice regarding the reduction of WFI images with THELI. We would like to acknowledge support from the Deutsche Forschungsgemeinschaft through the Transregio Programme TR33 and through the Munich Excellence Cluster “Structure and Evolution of the Universe”, as well as from Deutsches Zentrum für Luft- und Raumfahrt through grant No. 50OR1601.

References

Appendix A: Notes on individual clusters

RXC J0014.3-3023

This cluster presents the highest level of substructure in our sample. Its exceptional nature was recently confirmed by Jauzac et al. (2016) with the latest data from the HST-Frontier Field campaign. Our dynamical results are very similar to those obtained by Owers et al. (2011) with the same spectroscopic data set. In particular, we find that the cluster core is made of two high-velocity substructures, whose spatial overlap produces the main peak in the galaxy surface density map. As for the mass estimate of the whole system, we can mention the weak-lensing results of Medezinski et al. (2016), who reported M₂₀₀ = 2.06 ± 0.42 × 10¹⁵ M_⊙ when using a single mass component. Their multi-halo model leads to a smaller total mass of M = 1.76 ± 0.23 × 10¹⁵ M_⊙. Using their R₂₀₀ to rescale our results, we find that the Jeans analysis gives a mass M(r< 2.35 Mpc) = 1.81 ± 0.18 × 10¹⁵ M_⊙ after the removal of substructures, which agrees well with their findings. The caustic analysis returns a smaller mass M = 1.48 ± 0.50 × 10¹⁵ M_⊙ that agrees nonetheless with the other estimates within their 1σ uncertainties. This is because the removal of the two central high-velocity substructures leaves the PPS empty below ~ 0.5 Mpc (see Fig. C.3). As a consequence, the caustic amplitude is largely underestimated within this region, thus giving a smaller mass when integrated to larger radii. Prior removing the substructures, the caustic mass within the same aperture is M = 2.36 ± 1.10 × 10¹⁵ M_⊙, which also agrees with the lensing estimates. Our work presents the first dynamical mass estimates for this cluster.

RXC J0225.9-4154

A detailed analysis of this cluster was presented in Foëx et al. (2017). On the NE side of the cluster core lies an infalling galaxy clump. On the SW part, two additional clumps are found with a small rest-frame velocity, indicating a future merger close to the plane of the sky. The largest of these two clumps has a mass similar to that of the cluster core, hence the dynamical analysis most likely underestimate the total mass of the system by a factor ~ 2. Several other smaller galaxy clumps are found along the NE-SW axis, extending over at least 8 Mpc. The X-ray emission of the cluster core shows a secondary peak, probably due to a recent merger.

RXC J0516.6-5430

This cluster possesses a large amount of substructures in its galaxy surface density map (second largest Δ). It has a bimodal core containing three BCGs and several other galaxy clumps are distributed along the N-S axis (see Fig. C.6). These remarks suggest that RXC J0516 is observed at an early stage of its formation history. This is also confirmed by the presence of high-velocity galaxies in the substructures north and south from the core, and by the complex X-ray morphology that is also elongated on the N-S axis (see Fig. B.3 of Weißmann et al. 2013). McInnes et al. (2009) made a weak-lensing analysis of this cluster, as part of the Southern Cosmology Survey. They obtained a mass M₂₀₀ = 0.56 ± 0.25 × 10¹⁴ M_⊙. In the corresponding radius, R₂₀₀ ≈ 1.5 Mpc, the Jeans analysis gives a mass M(r< 1.5 Mpc) = 1.06 ± 0.13 × 10¹⁵ M_⊙ after the substructure analysis; the caustic mass at the same radius is M = 1.05 ± 0.22 × 10¹⁵ M_⊙. Both estimators return a mass roughly twice larger, even after excluding the high-velocity substructures. The same authors provide a mass estimate based on the cluster optical luminosity, M(L₂₀₀) = 1.7 × 10¹⁵ M_⊙, which is in good agreement with our M₂₀₀ derived from the Jeans analysis. The consistency of the dynamical estimates, their relatively good agreement with the hydrostatic value, and the high richness and optical luminosity indicate that their lensing mass is most likely underestimated. Our work presents the first dynamical mass estimates for this cluster.

RXC J0528.9-3927

A detailed analysis of this cluster was presented in Foëx et al. (2017). It is located in a rather poor environment, and has a simple dynamical structure (no substructure from the DS test, smallest f_DS value). A second BCG 200 kpc north from the central galaxy and a couple of small substructures in the galaxy surface density map suggests an accretion history along the N-S axis. The X-ray emission peak is slightly offset along the same N-S axis from the centroid estimated at larger scale, which also indicates that the cluster has not yet reached a fully relaxed state.

RXC J0658.5-5556

The Bullet cluster is one of the most studied galaxy clusters owing to the peculiar configuration of a merger taking place nearly in the plane of the sky (e.g. Markevitch et al. 2004). The first dynamical analysis of this object was performed by Barrena et al. (2002). They applied the virial theorem estimator (accounting for the SPT) on 71 galaxies located within an aperture of radius ~ 1.5 Mpc. They extrapolated the resulting mass to the radius ${\mathit{R}}_{\mathrm{200}}\mathrm{=}\sqrt{\mathrm{3}}{\mathit{\sigma }}_{\mathit{\upsilon }}\mathit{/}\mathrm{\left(}\mathrm{10}{\mathit{H}}_{\mathit{z}}\mathrm{\right)}\mathrm{\approx }\mathrm{2.65}$ Mpc to obtain M₂₀₀ = 1.24 × 10¹⁵ M_⊙. After the removal of substructures, our Jeans analysis gives a mass M(r< 2.65 Mpc) = 1.91 ± 0.25 × 10¹⁵ M_⊙, while the caustic estimator returns M = 1.26 ± 0.37 × 10¹⁵ M_⊙. The latter agrees well with their estimate, whereas the Jeans mass is somewhat larger. Recently Melchior et al. (2015) conducted a weak-lensing analysis of the Bullet cluster, finding a mass M₂₀₀ = 1.3 ± 0.6 × 10¹⁴ M_⊙. At the corresponding radius R₂₀₀ ≈ 2 Mpc, we find M = 1.49 ± 0.18 × 10¹⁵ M_⊙ and M = 1.11 ± 0.26 × 10¹⁵ M_⊙ for the Jeans and caustic estimator, respectively. The Jeans analysis gives again a larger mass, however, both values agree perfectly with the weak-lensing estimate. Interestingly, the hydrostatic mass agrees with the dynamical mass too, even though the cluster is known for being dynamically young (the merger takes place on the plane of the sky, hence it does not affect the dynamical measurements).

Regarding the structure of the cluster, it is worth making a few comments. On the one hand, the Bullet did not appear in the DS test. We added a KMM partition anyway and found a rest-frame velocity υ ≈ 550 km s^-1, which is slightly smaller than the value quoted by Barrena et al. (2002). On the other hand, the DS test found a group of ~ 10 high-velocity galaxies (υ ≈ − 1200 km s^-1) between the Bullet and the main clump (blue circles in Fig. C.7); this high-velocity substructure is clearly visible in the PPS diagram (Fig. C.3) and is responsible for the bimodal velocity distribution (Fig. C.2). The cluster itself is at the centre of a rich larger-scale environment with a clear filamentary-like extension towards south and an additional large galaxy clump northwest from the Bullet (see Fig. C.7). The dynamics from the southern galaxies suggest that the filament is located behind the cluster (the KMM partition associated with the blue triangles has a rest-frame velocity υ ≈ − 1000 km s^-1). Finally, we can note that Melchior et al. (2015) found a similar NW-SE orientation in the galaxy distribution at even larger scales (see their Fig. 8).

RXC J1131.9-1955

A detailed analysis of this cluster was presented in Ziparo et al. (2012). Our analysis confirms that the cluster possesses a complex structure that is characteristic of a dynamically young object. Several galaxy overdensities are found within its virial radius. Two of them present a specific dynamics that can be associated to two filamentary structures located north (positive rest-frame velocity) and south (negative rest-frame velocity) from the core. The latter filament can be traced further south with two additional KMM partitions (see Fig. C.7). Ziparo et al. (2012) used the scaling relation of Biviano et al. (2006) to estimate the mass of the cluster’s main body from its velocity dispersion. They obtained M₂₀₀ = 1.1 ± 0.2 × 10¹⁵ M_⊙ and a corresponding R₂₀₀ ≈ 1.9 Mpc. Within the same aperture, the Jeans analysis and the caustic estimator give, after the substructures removal, a mass M(r< 1.9 Mpc) = 1.26 ± 0.15 × 10¹⁵ M_⊙ and M = 1.16 ± 0.64 × 10¹⁵ M_⊙, respectively. Both values agree with the prediction of the M(σ_υ) scaling law.

RXC J1206.2-0848

Our analysis of RXC J1206 is based on the same data set that was used for the detailed dynamical analyses of Biviano et al. (2013), Lemze et al. (2013), and Girardi et al. (2015). Therefore, this cluster provides a consistency check for our results. In particular, Biviano et al. (2013) reported a mass M(r< 1.96 Mpc) = 1.37 ± 0.18 × 10¹⁵ M_⊙ from the Jeans analysis, and M(r< 2.08 Mpc) = 1.63 ± 0.58 × 10¹⁵ M_⊙ from the caustic method. Once rescaled to the same radii, our implementation of the Jeans approach gives M = 1.51 ± 0.14 × 10¹⁵ M_⊙, while the caustic mass profile leads to M = 1.70 ± 0.64 × 10¹⁵ M_⊙, that is, a perfect agreement for both estimators. Different lensing analyses of the cluster have found very similar results with M₂₀₀ ≈ 1.5 × 10¹⁵ M_⊙ (see e.g. the latest work by Sereno et al. 2017).

RXC J1347.5-1144

This cluster, which is the brightest in the REFLEX sample (Böhringer et al. 2004), has generated a great deal of discussion regarding its mass estimate using either X-ray, SZ, lensing, or dynamical data sets (see e.g. Cohen & Kneib 2002; Gitti et al. 2007; Bradač et al. 2008). Most of the debate originates from the first dynamical analysis by Cohen & Kneib (2002), who reported a velocity dispersion σ_P = 910 ± 130 km s^-1, whereas lensing analyses consistently find an equivalent σ_P ~ 1300 − 1500 km s^-1 (e.g. Fischer & Tyson 1997; Bradač et al. 2008; Lu et al. 2010; Foëx et al. 2012; Hoekstra et al. 2015). The merger of a substructure located close to the core has been invoked to explain the large discrepancies between the different mass estimates. Such a substructure can boost the strong-lensing signal, leading to an overestimated mass when extrapolated to large radii, while shocks associated with the merger can significantly enhance the gas temperature and X-ray luminosity, leading to a biased high hydrostatic mass. On the other hand, if the merger is taking place on the plane of the sky, a dynamical analysis would fail at estimating the total mass of the system. The lack of substructure found by the DS test and the NE-SW elongation of the galaxy distribution (see Fig. C.7) support the merger configuration described by Kreisch et al. (2016).

Our dynamical results tend to agree with the original work of Cohen & Kneib (2002) since they obtained a velocity dispersion σ_P = 820 ± 110 km s^-1 after a 3σ clipping selection of the cluster members. Using this value and accounting for the SPT correction at a radius of 1 Mpc, their mass estimate becomes M ≈ 4 × 10¹⁴ M_⊙, which is perfectly consistent with our results (interpolated to the same radius) M = 3.9 ± 0.5 × 10¹⁴ M_⊙ and M = 5.2 ± 1.6 × 10¹⁴ M_⊙ for the Jeans analysis and caustic estimators, respectively. More recently, Lu et al. (2010) reported a dynamical mass M₂₀₀ = 1.16 ± 0.3 × 10¹⁵ M_⊙ (at a corresponding R₂₀₀ ≈ 1.85 Mpc) under the assumption of a singular isothermal density profile. They also obtained a caustic mass M ~ 1 × 10¹⁵ M_⊙ within the same aperture. Our analysis gives M = 0.75 ± 0.1 × 10¹⁵ M_⊙ and M = 0.91 ± 0.3 × 10¹⁵ M_⊙ for the Jeans and caustic estimators, respectively. The former is somewhat smaller, as expected given the velocity dispersion σ_P = 1163 ± 97 km s^-1 reported by these authors. They also derived a weak-lensing mass M = 1.47 ± 0.45 × 10¹⁵ M_⊙, which is again larger than the dynamical measurements. Several galaxy clumps found in the vicinity of the cluster, in particular in the SW region, could contribute to the shear signal, thus explaining the larger weak-lensing masses.

RXC J2011.3-5725

This cluster is embedded in a rich star field. It was observed only in the B and V bands, making the photometric redshifts much less accurate, and thus preventing a precise characterisation of the cluster morphology. The DS test finds a small high-velocity substructure right in the cluster centre, so we removed the corresponding KMM partition before re-estimating the dynamical masses. Since the results of the DS test could be the consequence of low statistics, we also verified that our conclusions in Sect. 5.5 are essentially unchanged if we do not remove this KMM partition prior to the mass estimation. Our work presents the first dynamical results for this cluster.

RXC J2308.3-0211

A detailed analysis of this cluster was presented in Foëx et al. (2017). It is part of a supercluster, with several galaxy over-densities distributed along two main axes, N-S and E-W. The cluster core shows signs of a recent merger with significant X-ray residuals and possibly two high-velocity substructures following orbits close to the line of sight. Newman et al. (2013) obtained a combined strong+weak lensing mass M₂₀₀ = 1.32 ± 0.12 × 10¹⁵ M_⊙ for a radius R₂₀₀ ≈ 2.05 Mpc. After removing the substructures, the Jeans analysis gives a mass M(r< 2.05 Mpc) = 1.22 ± 0.13 × 10¹⁵ M_⊙ that agrees perfectly with their result. The caustic estimator gives M = 1.06 ± 0.23 × 10¹⁵ M_⊙. As for RXC J0014, this lower value is the consequence of the PPS being empty in the central region after the removal of the two high-velocity substructures. RXC J0014 and RXC J2308 are two good examples for which our approach to treat substructures is not optimal for the caustic estimator. We note that prior to the substructure analysis, we obtain a caustic mass M(r< 2.05 Mpc) = 1.73 ± 0.41 × 10¹⁵ M_⊙ that agrees with the lensing value.

Appendix B: Substructure indicators

Figure B.1 shows the photometric and dynamical substructure indicators defined in Sects. 4.1 and 4.2. In Fig. B.2 we present the comparison between the photometric indicator and the X-ray morphology estimators; the latter were taken from Weißmann et al. (2013).

Appendix C: Intermediate results

We provide in the following some intermediate results of the photometric and dynamical analyses. In particular, Figs. C.6–C.8 present the spatial distribution of the spectroscopically confirmed members, the galaxy surface density contours, and the different KMM partitions associated with substructures.

Appendix D: Redshifts

Table D.1 contains the list of the photometric and spectroscopic cluster members used for this work. The full list is available at the CDS.

All Tables

All Figures