© ESO 2019

1. Introduction

Strong gravitational lensing is a powerful technique to probe the total projected mass distribution of the inner regions of galaxy clusters and, therefore, the dark matter (DM) distribution once the baryonic mass components are independently mapped. By characterizing the substructure of cluster cores on different scales and the mass density profile of the innermost regions, one can test the ΛCDM structure formation scenario and indirectly constrain DM physical properties by comparing reconstructed mass maps with the latest cosmological simulations (e.g., Diemand & Moore 2011). In recent years, dedicated Hubble Space Telescope (HST) observational campaigns, such as the Cluster Lensing And Supernova survey with Hubble (CLASH, Postman et al. 2012), the Hubble Frontier Field campaign (HFF, Lotz et al. 2017), and the Reionization Lensing Cluster Survey (RELICS, Coe et al. 2019) have provided multi-band data with unprecedented quality on a sizable set of massive galaxy clusters. In parallel, extensive spectroscopic campaigns with highly-multiplexing multi-slit instruments (e.g., CLASH-VLT, Rosati et al., in prep.; Balestra et al. 2016), and especially with the Multi Unit Spectroscopic Explorer (MUSE) integral field spectrograph on the Very Large Telescope (VLT; Bacon et al. 2012), have secured redshifts for hundreds of multiple lensed images in cluster cores, as well identified large samples of cluster galaxies. The combination of these new imaging and spectroscopic data sets have enabled the development of new high-precision parametric strong lensing models based on 50−100 bona fide multiple images per cluster and highly complete samples of cluster members (e.g., Richard et al. 2014; Grillo et al. 2015; Jauzac et al. 2015; Limousin et al. 2016; Kawamata et al. 2016; Lagattuta et al. 2017; Caminha et al. 2016, 2017a,b, 2019; Bonamigo et al. 2018). These new lens studies typically reach a root mean square difference between the observed and predicted positions of multiple images on the lens plane of Δ_rms ≈ 0.5″. The latter is commonly used as a simple figure of merit for model accuracy. Nonetheless, the relatively large number of parameters used to describe the cluster mass distribution still suffer from strong internal degeneracies. The generally complex cluster mass distribution in parametric lens models is generally separated into a cluster-scale diffuse component, made of one or more large halos, and a clumpy distribution traced by cluster galaxies (Natarajan & Kneib 1997). The latter describes the sub-halo population of DM halos (De Lucia et al. 2004; Giocoli et al. 2010) and hence the inner substructure of the cluster mass distribution (Springel et al. 2001). In the effort to limit the number of free parameters, this sub-halo component is generally modeled adopting two scaling relations, which link the internal velocity dispersion and size of each halo with the luminosity of member galaxies. A fixed mass-to-light scaling is thereby assumed for all sub-halos. Since strong lensing models constrain the total projected mass within each family of multiple image positions, a certain amount of degeneracy is always present between the velocity dispersion, the size (or equivalently the profile), and the shape of each halo component of the mass distribution. Given the cross-talks among different halo components in strong lensing models, such a degeneracy can lead to systematics in the reconstruction of the substructure at different scales, for example with a transfer of mass between the clumpy and diffuse mass components. As we demonstrate in this work, such a degeneracy can be broken and significantly reduced, on the scale of the sub-halos, by using an independent measurement of the internal stellar velocity dispersion of cluster members.

A combination of lensing and kinematic measurements has long been exploited in the study of the mass density profile of field early-type galaxies. This combination has proved to be particularly effective since the two diagnostics complement each other (e.g., Treu & Koopmans 2004; Czoske et al. 2008; Barnabè et al. 2009), breaking the mass-anisotropy and mass-sheet (e.g., Falco et al. 1985; Meylan et al. 2006) degeneracies of the dynamical and lensing analyses, respectively. Joint strong lensing and stellar-dynamical studies have been used to determine the average logarithmic density slope of the total mass inside the Einstein radius (Treu & Koopmans 2004; Koopmans et al. 2009) and to decompose the total mass distribution into luminous and dark components of the lens galaxies belonging to the Lenses Structure and Dynamics (LSD) and Sloan Lens ACS (SLACS) surveys. These two surveys have measured a remarkably homogeneous total (luminous and dark) mass density profile that is consistent with an isothermal one (i.e., ρ ∝ 1/r²) out to a few hundreds of kiloparsecs (see Gavazzi et al. 2007; Bolton et al. 2008). A one-component isothermal model is fully characterized by the value of an effective velocity dispersion, approximated within ≲3% by the value of the galaxy central stellar velocity dispersion (i.e., the velocity dispersion of the stars projected within a disk of radius R_e/8; see, e.g., Treu et al. 2006; Bolton et al. 2008). This result is theoretically supported by the Jeans equation for realistic stellar density distributions (e.g., Jaffe 1983; Hernquist 1990) embedded in a globally isothermal distribution (see Kochanek 1993), as observed in samples of nearby and luminous early-type galaxies (e.g., Kochanek 1994; Grillo et al. 2008).

A combination of the spatially resolved kinematics of the brightest cluster galaxy (BCG) and strong lensing modeling in the inner cores of massive clusters has also been used to constrain the central slope of their mass density profiles using simplified dynamical models (Sand et al. 2004; Newman et al. 2013). A first attempt to include cluster galaxy kinematics in cluster strong lensing models was made by Verdugo et al. (2007). A more extensive study was carried out by Monna et al. (2015) in Abell 383 at z = 0.187 with early CLASH data. They integrated the velocity dispersion measurements of 21 member galaxies, obtained with the Hectospec fiber spectrograph at the MMT, into a lens model based on the identification of nine multiply imaged systems (six of which spectroscopically confirmed). Thus, they show how typical degeneracies among lens model parameters describing the galaxy mass-to-light scaling relation can be significantly reduced. In particular, a meaningful scaling relation between galaxy halo truncation radii and velocity dispersions can be derived (see also Monna et al. 2017). Additional constraints on halo sizes can be obtained by modeling the surface brightness distribution of strong lensing features around individual member galaxies (see also Suyu & Halkola 2010; Eichner et al. 2013 who did not use galaxy kinematics).

In our work, we extend the methodology introduced by Monna et al. (2015), by using a much improved HST and spectroscopic data set on three CLASH/HFF clusters, specifically an extended set of velocity dispersion measurements obtained with the MUSE spectrograph at the VLT. A robust characterization of cluster sub-halo populations, particularly the distribution of their sizes and masses, as well as their abundance, offer a critical test of the predictions of cosmological simulations (e.g., Limousin et al. 2009; Grillo et al. 2015; Munari et al. 2016; Natarajan et al. 2009, 2017). Such a comparison can shed light on baryonic processes shaping cluster substructure, dynamical processes leading to halo stripping in different environments, and indirectly on the nature of dark matter (Despali & Vegetti 2017; Chua et al. 2017; Nipoti et al. 2018; Niemiec et al. 2019).

This work is organized as follows. In Sect. 2, we describe our imaging and spectroscopic data sets. In Sect. 3, we detail how internal velocity dispersions of member galaxies are measured, including spectral simulations to assess their robustness. Strong lensing models for the three clusters under study are described in Sect. 4, while the specific methodology to incorporate galaxy kinematics information into our lens models is discussed in Sect. 5. Results are discussed in Sect. 6 where we also present the velocity dispersion functions for the three clusters. In Sect. 7, we summarize the main conclusions of our study.

Throughout this article, we adopt a flat ΛCDM cosmology with Ω_m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹. With these parameters, 1″ corresponds to a physical scale of 5.68, 5.34 and 4.92 kpc at z = 0.439, 0.396, 0.348, respectively the redshift of the three clusters of this study: MACS 1206, MACS 0416 and AS1063 (see below). All magnitudes refer to the AB system.

2. Photometric and spectroscopic data

This section summarizes the photometric and spectroscopic data sets for the three clusters used in this work, namely MACS J1206.2−0847, MACS J0416.1−2403, and Abell S1063 (a.k.a. RXJ 2248.7−4431), hereafter MACS 1206, MACS 0416 and AS1063, at redshifts 0.439, 0.396 and 0.348, respectively (see Table 1). These clusters were observed with HST in 16 broad band filters, from UV to near-IR, as part of the CLASH program. HST imaging of MACS 0416 and AS1063 was significantly augmented with the HFF program by adding deep exposures in seven filters (F435, F606W, F814W, F105W, F125W, F140W, F160W). The three clusters were also part of an extensive spectroscopic campaign with the CLASH-VLT Large program (P.I. P. Rosati), using the VIMOS high-multiplexing spectrograph, which provided over 4000 redshifts in each of the three clusters, over an area of ∼25 × 25 arcmin². These data sets yielded approximately 600 spectroscopic members for MACS 1206 (Biviano et al. 2013; Girardi et al. 2015), 900 members for MACS 0416 (Balestra et al. 2016) and over 1200 members for AS1063 (Mercurio et al., in prep.). Spectroscopic information in the cores of the three clusters has been significantly enhanced with the MUSE integral field spectrograph at the VLT, which is at the basis of the kinematic measurements presented in this work. These data, which are described in more detail below, have enabled new high-precision strong lensing models based on large samples of multiply lensed sources (see Caminha et al. 2017a,b, 2016 and Sect. 4). MUSE has a field of view of 1 arcmin², a spatial sampling of 0.2″, a spectral resolution of ∼2.4 Å over the spectral range 4750−9350 Å, with a spectral sampling of 1.25 Å pix⁻¹.

MACS 1206. MUSE data were obtained between 2015 and 2016¹, the redshift measurements of the member galaxies and a large sample of multiple images were presented in Caminha et al. (2017b; hereafter C17b). Three MUSE pointings cover a total area of 2.63 arcmin², mapping the SE–NW elongation of the cluster. The exposure time is 8.5 h in the central ∼0.5 arcmin² and 4 h in the remaining area.

MACS 0416. MUSE archival observations used in this work were presented in Caminha et al. (2017a; hereafter C17a), along with the redshift catalog, and consist of two pointings², one of 2 h in the NE region, with a seeing of 0.5″, and a SW pointing of 11 h, with a seeing of 1″.

AS1063. MUSE data consist of two pointings, which were presented in Karman et al. (2015, 2017). The SW pointing³ has an exposure of 3.1 h and seeing ∼1.1″, the NE pointing⁴ has an exposure of 4.8 h and seeing of 0.9″.

3. Cluster members: spectral extraction and internal kinematics

In this section, we describe the methodology adopted to extract the spectra of the cluster members from the MUSE data-cubes and to measure their internal stellar velocity dispersions. Catalogs of cluster members for MACS 1206, MACS 0416, AS1063 were presented in C17b, C17a and Caminha et al. (2016; hereafter C16), respectively, with the main objective of identifying the sub-halos to be included in the lens models (see below). Cluster members were defined as galaxies in the redshift intervals 0.425−0.453, 0.382−0.410, and 0.335−0.362 for the three clusters respectively, lying within a rest-frame velocity of approximately ±3000 km s⁻¹ around the median redshift of the cluster. With 114−145 spectroscopic members per cluster in the HST FoV, the extensive CLASH multi-band information was used to obtain highly (∼95%) complete and pure samples of photometric members down to m_F160W = 24, following the method described in Grillo et al. (2015).

The spectra of cluster members are extracted from the MUSE data-cubes within apertures of R_ap = 0.8″ radius. The latter is found to be a good compromise between the signal-to-noise of the spectra and the contamination from other sources in the field (either interlopers along the line-of-sight or close members). All the extractions are visually inspected to assess possible contamination from nearby sources and the apertures are reduced in specific cases down to 0.6″. The spectra in which the contamination from bright nearby members is too strong are discarded as the velocity dispersion of the fainter galaxy is likely biased.

We measure the stellar line-of-sight velocity distribution (hereafter LOSVD) of cluster members using the public software pPXF (penalized pixel-fitting, Cappellari & Emsellem 2004), with the latest improvements included in the 02/2018 python version (Cappellari 2017). With pPXF, we determine the best-fit LOSVD parameters (first and higher moments) by performing a cross-correlation of the observed spectrum with a set of spectral templates. The best-fit is obtained by minimizing a χ² between the template and the observed spectrum, including an additional penalty function whose weight is regulated by a parameter λ. This bias parameter is used to suppress higher velocity moments (h3, h4), when they become unreliable in spectra with low velocity dispersion and low S/N (see Cappellari & Emsellem 2004). For each spectrum, we measure a mean signal-to-noise, ⟨S/N⟩, over the selected wavelength range. We then adopt a relation between λ and ⟨S/N⟩, as suggested in Cappellari et al. (2011), which we tested with extensive spectral simulations described in Appendix A.

To perform the pPXF spectral fits, we use a subset of 105 stellar templates of different spectral types, drawn from the National Optical Astronomy Observatory library (Valdes et al. 2004). To match the typical underlying stellar populations of early-type galaxies in the cluster cores, most of the templates are of G, K, M spectral classes. In addition, we include 10 A-stars to extend spectral fits to a non-negligible fraction of E+A galaxies in our cluster sample and a few O and B stars. The stellar templates cover a wavelength range from 3465 Å to 9469 Å, with a sampling of 0.4 Å pix⁻¹ and have an intrinsic resolution of 1.35 Å full width at half-maximum.

In Appendix A, we describe in detail a set of simulations which reproduce our galaxy spectra with varying input LOSVD parameters, redshift and ⟨S/N⟩, to quantify the accuracy and precision of pPXF in measuring the velocity dispersions σ of member galaxies, thus optimizing pPXF input parameters. Simulated spectra are constructed from model spectra spanning the spectral type diversity of our cluster galaxy populations and by adding noise drawn from the variance map of the reduced MUSE data-cubes.

With these simulations, we check the reliability of the statistical error provided by pPXF on our data, as well as the presence of systematic errors as a function of ⟨S/N⟩ and input velocity dispersion (σ_in). The latter can become important especially when measuring velocity dispersions of low mass galaxies, as we approach the MUSE instrument resolution. We find that statistical errors are generally underestimated by ∼20% for ⟨S/N⟩ > 15, up to ∼25% for ⟨S/N⟩∼10, whereas a positive bias of a few km s⁻¹ becomes evident for σ ≲ 100 km s⁻¹ at high ⟨S/N⟩. In addition, we find that measurements become increasingly uncertain at ⟨S/N⟩ < 10. We therefore include in our galaxy kinematic sample only galaxies with ⟨S/N⟩ > 10 and σ >  80 km s⁻¹. In all cases, we use empirical formulas (see Appendix A) to correct the measured velocity dispersions and their uncertainties in different ⟨S/N⟩ and σ regimes.

Simulations were also used to choose the optimal wavelength range for pPXF fits, for a given galaxy redshift, by excluding regions of low S/N due to the MUSE sensitivity curve, particularly on the red side of the spectrum strongly affected by sky lines residuals. The resulting selected rest-frame wavelength ranges are 3600−4900 Å, 3600−5200 Å and 3600−5300 Å for MACS 1206, MACS 0416 and AS1063, respectively. The LOSVD input velocity value (zero moment) is taken from our redshift catalogs. The measured parameters are the velocity shift (V, typically within 50 km s⁻¹), the velocity dispersion σ, and higher moments (h3, h4), and their one-standard deviation errors. All velocity dispersions measured with pPXF are labeled in the following with the subscript “ap” to emphasize that these are line-of-sight quantities within an aperture, as opposed to σ₀, and σ_LT, which refer to parameters inferred with the lens models (see next session). Figure 1 shows examples of pPXF spectral fitting for three cluster galaxies in MACS 1206, one in the high ⟨S/N⟩ regime, one for the median ⟨S/N⟩ of our sample (∼30), and one at the limiting ⟨S/N⟩∼10.

Fig. 1. Results from line-of-sight velocity dispersion fitting of spectra of three cluster members in MACS 1206 as obtained with pPXF. Galaxy spectra are shown in black; red curves are the pPXF best-fit models, while the green points correspond to the data−model residuals. The blue shaded regions along the wavelength axis were excluded in the fitting procedure due to the presence of noisy sky subtractions around emission lines in the spectra. Corresponding residuals in these regions are marked in blue. The first from the top is a high signal-to-noise E+A galaxy. The second is a passive galaxy spectrum with a ⟨S/N⟩ = 28.1, corresponding approximately to the mean S/N of our galaxy sample. The bottom spectrum has a ⟨S/N⟩ = 10.7, close to our lower limit for reliable velocity dispersion measurements. Coordinates, measured velocity dispersion (σap), mean signal-to-noise and F160W magnitudes are indicated in each panel. Cutouts are HST RGB images, 4″ across, showing in green the apertures of 0.8″ radius used for the spectral extraction.

In Table 1, we quote the number of uncontaminated spectra extracted in each cluster (), for which we can reliably measure velocity dispersions, together with the total number of spectro-photometric members and relevant cluster parameters. In Fig. 2, the data points correspond to measured velocity dispersions of member galaxies as a function of their F160W magnitude, defining the Faber–Jackson relation in the three clusters.

Fig. 2. Data points are measured stellar velocity dispersions of cluster members, color-coded according to their spectral ⟨S/N⟩. Green lines are 300 scaling relations randomly drawn from the posterior distributions of the σap–mag scaling relation parameters, α, , Δσap, obtained from fitting the data points (see Fig. 3). Optimized parameters are also quoted. Red and blue areas are obtained by projecting the 3D σ–L scaling relations from our previous lens models with no kinematic prior (see text). Some velocity dispersion measurements, which deviate significantly from the scaling relations, are marked with squares and circles and discussed in the text.

4. Strong lensing models

Accurate strong lensing models were developed for MACS 1206, MACS 0416 and AS1063 in C17b, C17a and C16, respectively. These models were further refined in Bonamigo et al. (2018; B18 hereafter), who included the mass distribution of the hot gas component in each cluster, as derived from the Chandra X-ray data, dominating the smooth baryonic cluster component. As customary in cluster strong lensing modeling, none of these models included any kinematic information on cluster galaxies. We describe here our methodology which combines the B18 lens models with internal stellar kinematics derived from a large number of velocity dispersions measured with MUSE. For each cluster, we use the same catalogs of multiple images and cluster galaxies as in B18 and the Caminha et al. models, and likewise we employ the public software LensTool (Kneib et al. 1996; Jullo et al. 2007; Jullo & Kneib 2009).

A parametric lens model for the total mass distribution of each cluster is optimized searching for the set of parameters ξ, which minimize the χ² defined on the lens plane as (Jullo et al. 2007):

(1)

where identifies the number of multiple images associated to the same source j (usually called a family), N_fam is the number of families, x^obs are the observed positions of the multiple images on the lens plane, x^pred are their predicted positions, given the set of model parameters ξ, and Δx_i, j represent the uncertainties on the observed positions. Following B18, we perform a first optimization of our lensing models assuming a positional error of 0.5″ for images identified in HST and 1″ for those only found in the MUSE data. We multiply these errors by a constant factor ensuring that the best-fit χ² is close to the number of degree of freedom of the models. These updated errors are then used to sample the posterior distributions of the free parameters.

Fig. 3. Posterior probability distributions for σap–mag scaling relation parameters, obtained as described in Appendix B, from velocity dispersion measurements of 58 cluster members in MACS 1206. The 16th, 50th and 84th percentiles of the marginalized distributions for the slope (α), normalization () and scatter around the scaling relation (Δσap) are quoted and shown as vertical dashed lines.

Following the B18 modeling, the total mass distribution (or equivalently the gravitational potential ϕ) of each cluster is described as the sum of three contributions: (1) an elliptical large-scale smooth halo, which is further decomposed in a DM component and a smooth gas mass component, both described as dual pseudo-isothermal elliptical density (dPIE) profiles; the latter is obtained from deep Chandra observations, as multiple dPIE fits, as described in B18; (2) a clumpy component representing the cluster member galaxies (DM+baryons), modeled as spherical dPIE halos; (3) a shear+foreground-structure term to take into account the presence of massive structures in the outer cluster regions and line-of-sight mass distributions. Specific details on these multiple components are given below in the description of each lens model. The total cluster gravitational potential has therefore the form:

(2)

where the first sum runs over the N_h cluster-scale mass distributions, while the second on the N_g cluster member galaxies.

The general functional form for the spherical dPIE, including the relations between 3D and projected mass densities, as well as the expressions to derive aperture projected line-of-sight velocity dispersions, are given in Appendix C. Following our previous models and the general LensTool methodology, the dPIE parameters for the sub-halo population follow a scaling relation for the central velocity dispersion and the truncation radius:

(3)

(4)

where L_i is the luminosity of the ith cluster member and represents the corresponding truncation radius. A similar scaling relation for the core radius, r_core, is used in LensTool, however it is not relevant here since a vanishing core radius is adopted. is the LensTool fiducial velocity dispersion of each member, which is related to the central velocity dispersion of the dPIE profile by (see Appendix C). It is common practice to fix the slopes α and β_cut, so that the model optimization is performed over only two free parameters, that is the normalizations of the velocity dispersion and the truncation radius corresponding to the reference luminosity L₀. We measure the luminosities L_i and L₀ using the HST F160W Kron magnitudes, which are a good proxy of the stellar mass of the cluster members (see Grillo et al. 2015) and include members down to m_F160W = 24. This leads to a minimum of 193 sub-halos to be included in the lens model (see Table 1), fixed at the galaxy positions.

The total mass of a circular dPIE profile is given by (see Appendix C, Elíasdóttir et al. 2007; Limousin et al. 2005) the relation: .

Assuming a fixed scaling between the cluster members luminosity L_i (in the same band considered by Eqs. (3) and (4)) and its total mass M_tot, i, that is , one can obtain the following relation between the slopes of Eqs. (3) and (4):

(5)

The main goal of this work is to use prior information for the scaling relations derived directly from measured velocity dispersions of cluster galaxies. The latter are light-weighted projected values of the 3D-velocity dispersions within the extracted spectroscopic apertures, hereafter σ_ap. In order to compare σ_LT − L scaling relations with measured quantities, we need to compute the aperture-averaged line-of-sight velocity dispersion from the adopted dPIE mass models. The projection coefficients needed to transform σ_LT into σ_ap are obtained with a numerical integration depending on r_cut, r_core and R_ap. In Appendix C we recall all the equations to perform this projection procedure.

In Fig. 2, we show the best-fit scaling relations obtained from the previous models for each cluster (red and blue curves), which did not include any prior from internal kinematics of cluster members. These curves are computed by projecting the LensTool scaling relations, using the posterior distribution of the parameters r_cut and at different magnitudes, yielding the model σ_ap as a function of m_F160W. In this process, we use an aperture of 0.8″ so that the projected scaling relations can be directly compared with our kinematic measurements.

As customary in previous lens models to date, all the slopes for scaling relations α,  β_cut,  γ (see Eq. (5)) are fixed, here instead we take advantage of the measured stellar velocity dispersions to directly fit the normalization and slope α of the σ − L scaling relation. We then derive β_cut by adopting γ = 0.2, which is consistent with the canonical fundamental plane (Faber et al. 1987; Bender et al. 1992), and that we verify to be appropriate using photometric and morphological data for AS1063 (see below and Mercurio et al., in prep.).

In order to fit the σ–mag relation, a Bayesian approach is used, as described in Appendix B. Thus, we derive the value of α and , as shown in Fig. 3, which also includes the intrinsic scatter Δσ_ap of the measured velocity dispersions around the model (Eq. (3)). Similar behaviors for the parameters derived for the other two clusters can be found in Appendix B. The best-fit scaling relations obtained with this method are shown in Fig. 2, including the corresponding uncertainties (green curves), which are derived by sampling 300 times the posterior distributions. These best-fit parameters for each scaling relation (see Table 2) are then used as dynamical priors in our lens models, as described below. It is interesting to note that the slope that we obtain for the Faber–Jackson relation, L-σ^1/α, is very similar for the three clusters (α = 0.27−0.28) and consistent with several spectro-photometric studies of cluster early-type galaxy populations in the literature (e.g., Kormendy & Bender 2013; Focardi & Malavasi 2012).

We note that some galaxies deviate significantly from the best-fit scaling relations (as marked by boxes and circles in Fig. 2). This is however expected, as the Faber–Jackson relation is one of the projections of the fundamental plane relation among half-light radius R_e, mean surface brightness μ_e within R_e, and velocity dispersion: at a given luminosity, more compact galaxies tend to have higher velocity dispersion. For example, the four cluster galaxies at m_F160W ∼ 19 in AS1063, with a σ significantly higher than the best-fit relation (see boxes in Fig. 2), have μ_e in the 16th highest percentile, however we verified that they still lie on the fundamental plane defined by Jorgensen et al. (1996) (Mercurio et al., in prep.).

On the other hand, the few galaxies lying well below the σ–mag relation (see circled data points in Fig. 2 and the spectrum in the upper panel of Fig. 1) show sign of emission lines and young stellar populations in their spectra, for which lower velocity dispersions are expected when compared to early-type galaxies with similar luminosities.

5. Combining lensing models with kinematics measurements

In Fig. 2, a comparison between the sub-halo scaling relations obtained from the lens models and those directly constrained from kinematic measurements (green curves) shows significant discrepancies in the case of MACS 0416 and AS1063, while they are consistent for MACS 1206. In the case of AS1063, the normalization of the kinematic scaling relations are found to be approximately 100 and 150 km s⁻¹ above the values inferred from the B18 and C16 lens models, respectively. In MACS 0416, the discrepancy is significant (∼50 km s⁻¹) for the C17a model, while it is negligible for the B18 model within the errors, albeit with a slope which deviates from the observed one (α = 0.27 against 0.35 assumed in B18). This shows how inherent degeneracies of sub-halo population parameters in strong lensing models can lead to inferred velocity dispersion normalizations which are inconsistent with kinematic measurements of cluster galaxies. Nevertheless, these lens models can reproduce the positions of the multiple images with high precision, with a root-mean square value between the observed and model-predicted images on the image plane of Δ_rms ≃ 0.45−0.6″ (see Table 3). It is in fact well known that parametric cluster lens models are, in general, affected by some degeneracy between the mass distribution of the macro-halo(s) and that of the sub-halos, even when a large number of constraints are available, as in our case (Meneghetti et al. 2017). Despite that, the projected total mass value within a given cluster-centric radius remains robust. In addition, a significant degeneracy exists between the central velocity dispersions and the cut-off radii of the sub-halos.

In the following, we describe in detail our new lens models for each cluster. We start from the same parameterization and input constraints as in the B18 models, we then proceed to add critical constraints on the sub-halo scaling relation parameters, with priors from our kinematic measurements. The optimization of model parameters is obtained from MCMC chains of approximately 10⁵ samples, excluding the burn-in phase. Model components and parameters are summarized in Table 3. As customary in the literature, no scatter in the scaling relations is assumed in our lens models, as it is currently not possible to include it in LensTool (Bergamini et al., in prep.).

5.1. MACS 1206

Following C17b and B18 models, the cluster smooth mass distribution includes three dark matter halos which can reproduce the apparent elongated asymmetry in the distribution of the cluster galaxies and intra-cluster light. The three halos are described as dPIEs profiles whose values of sky positions, ellipticities, position angles, core radii, and velocity dispersions are left free to vary with flat priors, while their truncation radii are fixed to a large value. Three dPIE gas clumps are used to model the X-ray surface brightness distribution as in B18. Thus, the cluster-scale mass components have a total of 18 free parameters. The presence of an external shear term introduces two extra free parameters in the model. Finally, the clumpy component includes 258 halos describing the cluster members and the BCG, centered on the peaks of their light emission. All these galaxies are described as circular dPIEs profiles, whose values of central velocity dispersion and truncation radius scale with their F160W magnitude, according to Eqs. (3) and (4), where the BCG magnitude is used as reference luminosity (L₀). The normalization is free to vary between 1″ and 50″ (i.e. 5.67−283.9 kpc at the cluster redshift) with a flat prior. A Gaussian prior derived from our kinematic measurements is then introduced for the σ − L scaling relation. To determine this prior, we deproject the best-fit (median) normalization of the σ–mag relation ( km s⁻¹, see Fig. 3) to obtain a LensTool fiducial reference velocity dispersion, . This deprojection is achieved with an iterative procedure by computing the projection coefficient, c_p(r_core, r_cut, R_ap) (see Appendix C), using the best-fit scaling relation parameters from a first model without kinematic prior, as . In the case of MACS 1206, we thus obtain a Gaussian prior on with a mean of 264 km s⁻¹ and a standard deviation of 18 km s⁻¹. The standard deviation also includes the uncertainty in the computation of c_p, due to the allowed range of r_cut. Finally, the slope α is fixed to the value obtained from the σ–mag fit (α = 0.28), while a β_cut = 0.64 is obtained from Eq. (5), with γ = 0.2. The lens model is thus optimized with a total number of 22 free parameters, using the observed positions of 82 spectroscopically confirmed multiple images associated to 27 families (C17b and B18).

5.2. MACS 0416

Referring to Table 3, the cluster-scale mass components of this merging cluster (Balestra et al. 2016) include two massive dark matter halos, whose positions are left free to vary in areas of 30″ × 30″ and 15″ × 15″ around the northern and southern BCG, respectively (BCG,N, BCG,S); plus a third circular halo in the NW region to improve the model accuracy (Δ_rms) in that area. Four dPIEs with fixed parameters are used to describe the complex X-ray emitting gas distribution in MACS 0416, as done in B18. The cluster-scale component has therefore 16 free parameters. No shear is present, however an additional circular dPIE halo with free velocity dispersion and truncation radius is added to account for the presence of a foreground galaxy (Caminha et al. 2017a; Jullo et al. 2007).

The sub-halo mass component, including the BCG, comprises 193 halos in this case. The parameter is left free to vary over a wide range (1−20″, or 5.3−106.8 kpc). As described above for MACS 1206, the deprojection of the best-fit σ_ap–mag relation provides a Gaussian prior for the normalization of km s⁻¹, while we use the measured slope α = 0.27 (see Fig. B.1/top), from which we derive β_cut = 0.66 (from Eq. (5), with γ = 0.2). With a total number of 20 free parameters, 102 spectroscopically confirmed multiple images (associated to 37 families) are then used to constrain the lens model (as in C17a and B18).

5.3. AS1063

The dark-matter macro-halo includes in this case an elliptical dPIE, close to the BCG position, and a circular dPIE whose position is left free to vary in an area of 150″ × 120″ centered in the north-east region of the cluster (Table 3). Three dPIEs are used to describe the gas component (B18), and no shear term is included. There are therefore nine free parameters for the cluster-scale component.

The clumpy sub-halo component is constituted by 222 cluster members, including the BCG. From our kinematic measurements and best-fit scaling relation we obtain a deprojected value of km s⁻¹. As a result, we adopt the median value of 310 km s⁻¹ as a Gaussian prior for the normalization of the σ_ap − L scaling relation, while the standard deviation is reduced to 15 km s⁻¹ to force the lens model to converge to a χ² minimum solution compatible with our galaxy kinematics. The slope is fixed to the measured value α = 0.27 (see Fig. B.1, bottom), which implies β_cut = 0.66. The lens model, which has 11 free parameters in total, is thus optimized using the positions of 55 spectroscopically confirmed multiple images, divided in 20 families, as implemented in B18.

6. Results and discussion

The inclusion of the kinematic prior in the lens models has the main consequence of significantly reducing the intrinsic degeneracies in the scaling relation parameters. Specifically, since the projected mass of the sub-halos depends on the central velocity dispersion, σ₀, and truncation radius r_cut, for a vanishing core radius, there is clearly an inherent degeneracy between the normalization parameters of the σ − L and r_cut − L scaling relation. This can be easily appreciated from the expression of the dPIE projected mass within a radius R, for r_core = 0 (Eq. (C.5) in Appendix C, or Jullo et al. 2007):

(6)

Therefore, several combinations of σ₀ and r_cut will yield similar aperture masses, which are constrained by the multiple image positions. This degeneracy, in combination with those among the parameters describing cluster-scale mass components or shear terms, can lead the lens model to predict a distribution of sub-halo masses not consistent with internal kinematics of member galaxies. This can be noticed in Fig. 2 where we display the case of MACS 0416 and AS1063. The implementation of the kinematic prior on the normalization and slopes of the scaling relation leads to models with a global precision which is similar to that of the previous lens models, as apparent from a comparison of Δ_rms values (see Table 3). The new values compared to the B18 models differ by 0.01″ for MACS 1206, 0.02″ for MACS 0416, and 0.04″ for AS1063, with the same number of degrees of freedom (88 for MACS 1206, 110 for MACS 0416 and 59 for AS1063). The slight increase of Δ_rms in the new models is expected since the freedom of model parameters is generally reduced with the Gaussian priors.

In Fig. 4, we show the posterior probability distributions for the projected values of the normalization of the σ − L and r_cut − L relations, and , respectively, for the new (in green) and previous models. The kinematic priors lead in general to significant smaller uncertainties for these parameters, up to factor of 10 for the truncation radius of AS1063, with respect to the B18 results.

Fig. 4. Panels a–c: posterior probability distributions for the normalization of the σap − L scaling relation, that is aperture projected velocity dispersion , and truncation radius (), obtained with different lens models for clusters under study. The green contours and distributions refer to our best models which include the kinematic prior based on measured velocity dispersions of member galaxies (Fig. 2). For visualization clarity, in panel c we omit the C16 results because they are outside the chosen range of velocity dispersion, as visible in the bottom panel of Fig. 2. The reference magnitudes for each cluster are in Table 2. Panel d: posterior distribution of the reference central velocity () and reference truncation radius () in each cluster, normalized to an absolute mag , close to L* of the early-type galaxy population. The different colors correspond to the models with kinematic priors and contours refer to 1, 2, 3-σ confidence levels. Once normalized to the same luminosity, the three scaling relations for σ and rcut are consistent.

While for MACS 1206 the new solution is consistent with the kinematic data, for the other two clusters the kinematic prior moves the best-fit solution to a different region of the parameter space. In the case of MACS 0416, the B18 model had already found a different solution with a higher σ-normalization with respect to the C17a, by extending the MCMC parameter search and including the hot gas mass component. Our new lens model is consistent within 1σ with B18, further reducing parameters’ uncertainties especially on r_cut. In the case of AS1063, the kinematic prior moves the (for the BCG) to a much larger value, which is 70 km s⁻¹ higher, and provides a much tighter constraint for , moving the solution away from previous values of ∼190 kpc (a dPIE encircles 90% of the total projected mass within 5r_cut).

It is relevant to note that, despite the higher normalization of the sub-halo scaling relations for MACS 0416 and AS1063, the cluster total projected mass profile does not vary appreciably when compared to previous models. This is expected since the multiple image positions provide information about the cluster total mass projected within circles with radii equal to the average distance of the multiple images of each family. In Fig. 5, one can appreciate that the differences in the cumulative projected mass between the new and previous models are well within 10% over the radial region with multiple image constraints. Therefore, although the fraction of total mass in sub-halos differs for MACS 0416 and AS1063, the total encircled mass remains a robust measurement. This however shows that the inclusion of the kinematic prior is important to obtain a reliable estimate of the relative contribution of the diffuse and clumpy (sub-halo) mass components from strong lensing models. In general, an underestimate of the sub-halo mass component will lead to a higher contribution of the cluster-scale halos, in particular the central one.

Fig. 5. Top: projected cumulative mass profiles, as a function of the projected distance from BCGs, corresponding to previous and current lens models for each cluster. The colored regions encompass the 16th and the 84th percentiles; colored solid lines are the median values. The dashed green lines correspond to the mass component associated to cluster members (i.e. sub-halos) from the new lens models with kinematic prior (16th, 50th, 84th percentiles). The multiple image projected distances from the cluster centers are marked with vertical black lines. Bottom: relative variation of the cumulative projected total masses with respect to our reference (green) model.

6.1. Sub-halo scaling relations

The consistency of the sub-halo scaling relations, based on kinematic measurements in the three clusters, suggests a rescaling of the normalization of the σ − L and r_cut − L relations to the same absolute luminosity for a meaningful comparison. To this aim, we choose as reference an absolute magnitude of , which is close to the value of L^*, estimated from the F160W luminosity function of MACS 0416 at z = 0.4 (Connor et al. 2017; Mercurio et al., in prep.). We then rescale the normalization of the scaling relations obtained from the kinematic lens models, anchored on the BCG magnitude (), to the new reference using the distance modulus μ of each cluster, as follows:

(7)

where the central reference velocity dispersion is in Eq. (3), and

(8)

The comparison of the posterior distributions of σ₀, and r_cut, obtained from our lens models with kinematic prior, once renormalized to the same L^* luminosity (i.e. ), is shown in the bottom right panel (d) in Fig. 4. The consistency of these distributions among the three clusters is quite remarkable, particularly for r_cut, and suggests an empirical relation for the truncation radius, which is generally poorly constrained by lens models without kinematic priors. By combining Eqs. (3) and (4), one obtains:

(9)

where the range for each parameter represents the 16th−84th percentiles of the combined posterior distributions of , and , and the slope is obtained by the distribution of β_cut/α. A similar relation was derived by Monna et al. (2015) for Abell 383 (z = 0.18) incorporating the velocity dispersion measurements of 21 member galaxies. However, in that study the slopes of the σ ∝ L^α and r_cut ∝ L^β_cut relations are fixed, based on literature studies, assuming M_tot/L = const. In our case, the slope α is directly measured from the MUSE spectra of a cluster galaxy sample larger by more than a factor of 2, while β_cut is derived assuming galaxies lying on the fundamental plane.

The slope β_cut/α = 2.43 that we find for the r_cut − σ₀ relation is likewise the result of our kinematic measurements in three different clusters, whereas the slope of 1.25 quoted in Monna et al. (2015) is derived from the assumed values of α and β_cut. However, the normalization of the relation in Eq. (9) turns out to be consistent with the Monna et al. (2015) result. Other attempts to constrain halo sizes of cluster galaxies have involved modeling of single strong lensing systems in clusters, with no kinematics information from spectroscopic measurements. For example, Suyu & Halkola (2010) found a kpc for a galaxy with a lens-based velocity dispersion km s⁻¹ in a group, in agreement with our results. Eichner et al. (2013) found instead a higher normalization by modeling the surface brightness distribution of the “snake arc” in MACS 1206. However, we should note that the lens model used in that work was based on a first limited sample of multiple images, when MUSE spectroscopy was not available. Limousin et al. (2007) estimated a scaling relation between halo sizes and σ₀ using galaxy-galaxy weak lensing in 5 clusters at z ∼ 0.2, which is broadly consistent with our results.

The consistent constraints we find on the central velocity dispersions and truncation radii of the cluster sub-halos from the independent analysis of three clusters suggest that we can combine Eqs. (C.6) and (9) to obtain an empirical relation between σ₀ and the total mass of member galaxies. By propagating the uncertainties derived from the posterior distributions of σ₀ and r_cut, we obtain the following M_tot − σ₀ relation:

(10)

where the parameter range refers to the 16th and 84th percentiles.

6.2. Sub-halo velocity dispersion function

In light of our results, we present a central velocity dispersion function of cluster members, which can be derived from the sub-halo component of the lens models. This was first presented, in the form of circular velocity function (, see Appendix C) in Grillo et al. (2015) for MACS 0416 and later extended by B18 to the same sample of three clusters of the present study, including the hot-gas components. The new sub-halo velocity functions presented here are particularly relevant since the new lens models incorporate a prior on the σ − L scaling relation based on measured velocity dispersions of a large number of cluster members.

Each data point in Fig. 6 is obtained by extracting the central velocity dispersion of each member from 500 random realizations of the model posterior distributions, thereby computing the median number of members in each velocity bin 35 km s⁻¹ wide. We adopt the same bin-size as in B18 to allow a direct comparison with previous results. The vertical error bars are obtained by combining the error associated to the 16th and 84th percentiles of these realizations with a Poissonian error generally appropriate for low number counts⁵. In order to compare the velocity functions of the three clusters, we include only member galaxies with cluster-centric distances R <  0.16 × R_200c (see Table 1), which is the maximum aperture within which we have a highly complete sample of member galaxies in all three clusters.

Fig. 6. Sub-halo central velocity dispersion function derived from our lens models for MACS 1206, MACS 0416 and AS1063 with kinematic prior based on measured velocity dispersions of cluster members. Previous models with no kinematic prior for AS1063 and MACS 0416 are shown for comparison. Data points correspond to median galaxy counts in 35 km s−1 bins, drawn from the posterior distributions of sub-halo model parameters. Errors also include Poissonian number statistics in each bin. Only galaxies within a cluster-centric radius of 0.16 × R200c are included (see Table 1). The top axis gives the total mass related to the central velocity dispersion by Eq. (10) (see text).

In Fig. 6, we compare the new velocity functions with the previous determinations by Grillo et al. (2015; for MACS 0416, dashed blue line) and B18 (for AS1063, dashed green line). Not surprisingly, these previous determinations were biased low due to the lower normalization of the scaling relations in previous models, whereas the new lens models with kinematic priors produce velocity functions which are quite consistent with each other.

7. Conclusions

In this work, we have improved our previous strong lensing models of three Hubble Frontier Fields/CLASH clusters, based on a large number of spectroscopically confirmed multiple images, using robust measurements of the internal stellar velocity dispersion of large samples of clusters galaxies, thanks to MUSE integral-field spectroscopy in cluster cores. Such measurements constrain independently the normalization and slope of the Faber–Jackson σ − L relation so that the sub-halo component of the cluster mass distribution in the lens models is now bound to the kinematic measurements, thus allowing a significant improvement on the reconstruction of the cluster internal sub-structure. We summarize below the main results of our study.

Using spectra of cluster early-type galaxies with a mean S/N of 25, we measure robust velocity dispersions of 37, 49 and 58 members in AS1063, MACS 0416, MACS 1206, whose accuracy is tested with extensive spectral simulations. By obtaining kinematics measurements over 4−5 mag in each cluster, down to m_F160W ≃ 21.5 which corresponds to ∼2.5 mag below L^*, we sample well the Faber–Jackson relation. A maximum likelihood modeling of this relation yields well consistent normalizations and slopes for the three clusters.

The new lens models incorporating kinematics information of cluster galaxies reproduce the positions of multiple images with similar precision to previous models. While the total projected cluster mass profile remains essentially unchanged, the mass of sub-halos are now robustly constrained, thus reducing degeneracies with other mass components and parameters of the lens model.

Specifically, the inherent degeneracy of lens models between the central velocity dispersion (σ₀) and truncation radius (r_cut) of sub-halos is strongly reduced, thus providing robust estimates of sub-halo masses and sizes. Once normalized to the same absolute luminosity, the three σ₀ − r_cut scaling relations, independently derived for each cluster, are found in very good agreement. This is particularly interesting in consideration of the different dynamical states of the three systems. As a result, we obtain a statistical determination of sub-halo truncation radii with masses ranging from that of the BCG down to ∼10¹⁰ M_⊙. Our findings extend the initial results derived by Monna et al. (2015) in A383 by using now a high-quality data set which includes precision lens models based on large number bona-fide multiple images and high S/N kinematic measurements of cluster members.

Interestingly, with such a robust determination of the scaling relations of the sub-halo populations, we infer fully consistent velocity dispersion functions for the three clusters, unlike in previous lens models. In addition, our new constraints on the sub-halo masses provide an empirical M_tot − σ₀ relation which can be used to translate the velocity functions into sub-halo mass functions.

We plan to extend this methodology to a larger sample of CLASH clusters with strong lensing models based on MUSE observations (Caminha et al. 2019) to better explore the role of systematics in parametric lens models when constraining different mass components of galaxy clusters. Further constraints on the mass profile of cluster galaxies, including the dark matter fraction, will require the use of galaxy-scale strong lensing systems (e.g., Suyu & Halkola 2010; Grillo et al. 2014; Monna et al. 2015), in combination with the internal kinematics of the lenses to extend in dense environments the extensive work carried out in field early-type galaxies. A complementary analysis based on the statistics of galaxy-galaxy strong lensing events will also be presented in an upcoming paper (Meneghetti et al., in prep.). We defer to future papers a detailed comparison of the sub-halo mass function and halo-size distributions derived from kinematic and lensing data with results from different cosmologicalsimulations.

Acknowledgments

AA was supported by a grant from VILLUM FONDEN (project number 16599). This project is partially funded by PRIM-MIUR 2017WSCC32 “Zooming into dark matter and proto-galaxies with massive lensing clusters” and the Danish council for independent research under the project “Fundamentals of Dark Matter Structures”, DFF–6108-00470. MM acknowledges support from the Italian Space Agency (ASI) through contract “Euclid – Phase D” and from the grant MIUR PRIN 2015 “Cosmology and Fundamental Physics: illuminating the Dark Universe with Euclid”. CG acknowledges support by VILLUM FONDEN Young Investigator Programme through grant no. 10123.

References

Appendix A: pPXF spectral simulations

A.1. Simulated sample

We describe here our simulations of galaxy spectra to assess the robustness of the velocity dispersion measurements with the pPXF software (Cappellari 2017) which are critical in our lens model analysis. A sample of 10 000 simulated spectra of early-type galaxies were created according to the following steps:

• (1)

  From the MACS 1206 MUSE datacube we extracted all the spectra of spectroscopic members with m_F160W <  24 and ⟨S/N⟩ > 20, using an aperture of R = 0.8″.

• (2)

  Using pPXF best-fit to the spectra, we determined the stellar templates and relative weights, whose linear combination reproduce our galaxy spectra, which include younger stellar populations in some cases with significant Balmer lines. From this set of template, we created rest-frame synthetic model spectra that is not yet convoluted with a LOSVD and without noise, with a dispersion of 0.4 Å per pixel and a spectral resolution of 1.35 Å FWHM.

• (3)

  Model spectra are convoluted with a LOSVD up to the 4th moment, with velocity dispersion ranging from 0 to 250 km s⁻¹.

• (4)

  We assigned to each model spectrum a redshift so to reproduce the observed redshift distribution of cluster galaxies. The spectral resolution is then degraded to the lower MUSE instrumental resolution of 2.6 Å FWHM approximately constant over the whole considered wavelength range.

• (6)

  The spectra wavelength-range is matched to the MUSE range and re-sampled to MUSE 1.25 Å pix⁻¹ scale. Gaussian noise, drawn from the MUSE variance spectra, is added to the spectra in such a way to reproduce a ⟨S/N⟩ from 5 to 100. In this way, the high variance due to sky subtraction around sky lines is included in the simulations.

The final set of simulations is divided into six subsets that differ for their input σ and ⟨S/N⟩ (see Table A.1). The first (V_in), third (), and fourth () moments of the LOSVD are uniformly distributed within the intervals (−50.0, 50.0) km s⁻¹, (−0.1, 0.1), and (−0.1, 0.1), respectively.

First, we investigated the possible presence of systematics, especially at low ⟨S/N⟩ and σ ≲ 100 km s⁻¹, where the MUSE instrumental resolution may affect our measurements. Figure A.1 shows that an increasing bias in velocity dispersion measurements at low velocities is detectable in the high ⟨S/N⟩ regime (≫60). Using a polynomial function to fit the median values of the σ − σ_in distributions in seven σ bins, we found a correction for measured velocity dispersions, σ, given by:

(A.1)

Fig. A.1. Top: difference between pPXF measured velocity dispersions (σ) and input velocity dispersions in simulated spectra (σin), as a function of σ. Only the galaxies with ⟨S/N⟩ > 60 are shown here. The data points are color shaded according to the spectra ⟨S/N⟩. Black triangles are the medians of the distributions in 7 different σ bins. The black polynomial line underlines the small increasing bias at low σ. Bottom: black solid (blue dashed) line shows the relative statistical error estimated from simulations (provided by pPXF), as a function of ⟨S/N⟩. The red curve corresponds to the ratio between the two; the statistical relative errors appear to be underestimated by ∼20% in our spectra.

The increase in the σ − σ_in scatter at low ⟨S/N⟩ makes it difficult to detect such a small bias at lower ⟨S/N⟩. The scatter in σ − σ_in increases significantly below ⟨S/N⟩ = 10, with relative errors exceeding 10% and measurements become increasingly sensitive to the presence of sky lines. For this reason, we chose a limiting ⟨S/N⟩ of 10 in our analysis.

Another goal of the simulations was to quantify realistic errors on velocity dispersion measured by pPXF. In our case, we found that relative errors, δσ/σ, are systematically underestimated by ∼20% for ⟨S/N⟩≫15 and of ∼25% for ⟨S/N⟩∼10. In all our measurements, velocity dispersions are corrected for the small bias in the high ⟨S/N⟩ regime (Eq. (A.1)) and statistical errors are corrected using a polynomial fit, as shown in Fig. A.1/bottom.

Appendix B: σ–L scaling relation from kinematic measurements

Given a set of N cluster members of magnitude with measured velocity dispersions, , the σ–mag scaling relation, corresponding to the σ–L^α relation in Eq. (3), can be written as:

(B.1)

where are the model predicted velocity dispersions for a cluster member with F160W magnitude .

To estimate the model parameters of the scaling relation and their uncertainties, we sample the posterior distribution of and α, including the intrinsic scatter Δσ_ap of the measured velocities around the backbone of the scaling relation. Using Bayes’ theorem, the posterior probability function can be written as:

(B.2)

In particular, the posterior is the product of a likelihood function (B.3) and a prior (B.4):

(B.3)

(B.4)

The boundaries , , α_min, α_max, (Δσ_ap)_min and (Δσ_ap)_max were chosen to limit the parameter space around the measured velocity dispersions.

To sample the log-posterior in the 3D parameters space (, α, Δσ_ap), we use the Affine-Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler developed by Goodman and Weare (Goodman & Weare 2010), and in particular its python implementation⁶ (Foreman-Mackey et al. 2013). The parameter space is explored with 100 walkers, with 5000 steps each, which are initialized in a narrow Gaussian sphere around the maximum-likelihood point. To ensure that the final distributions are independent from the initial walker positions, we remove 80 steps in the burn-in phase based on the auto-correlation time computed for each parameter.

The posterior probability distributions of model parameters (Eq. (B.2)) obtained with this procedure are shown in Fig. 3 for MACS 1206 and Fig. B.1 for MACS 0416 and AS1063.

Fig. B.1. Posterior probability distributions for scaling relation parameters in Eq. (B.1) obtained from measured velocity dispersions of cluster members in MACS 0416 and AS1063. The 16th, 50th and 84th percentiles of the marginalized distributions for the slope (α), normalization () and scatter around the scaling relation (Δσap) are quoted and shown as dashed lines in each histogram.

Appendix C: dPIE mass distribution

The dual pseudo-isothermal elliptical mass distribution (dPIE), described in Elíasdóttir et al. (2007) and Limousin et al. (2005), is defined by eight parameters; three for the sky position and redshift (z), two for ellipticity (e) and position angle, and three for core radius (r_core), truncation radius (r_cut) and central velocity dispersion (σ₀). All the halos in our lens models are described as dPIEs profiles, cluster-scale halos are generally included in the limit r_cut → ∞; in this case, the dPIE reduces to a pseudo isothermal elliptical mass distribution (PIEMD), defined by Kassiola & Kovner (1993). Sub-halos associated to cluster members are instead parametrized as circular dPIEs of vanishing r_core, with r_cut and σ₀ following the scaling relations in Eqs. (3) and (4) (see Sect. 4).

The equations below, mostly derived in Elíasdóttir et al. (2007) and Limousin et al. (2005), refer to circular (e = 0) dPIEs. For completeness, we give the expressions for projected masses and surface-mass densities and show how the dPIE σ₀ parameter can be derived from the aperture average line-of-sight (projected) velocity dispersion σ_ap, which is associated to our kinematic measurements.

C.1. dPIE density and mass

The mass density profile for the spherical dPIE, as implemented in LensTool, is defined through the 3D-density (Limousin et al. 2005):

(C.1)

where r is the distance from the center of the mass distribution, while r_core and r_cut are the core and truncation radii respectively (with r_cut >  r_core). Equation (C.1) shows a smooth separation between an isothermal behaviour, ρ ∝ r⁻², where r_core <  r <  r_cut and a steeper decrease, ρ ∝ r⁻⁴, for r >  r_cut.

The relation between the central density ρ₀ and the 1D-central velocity dispersion, σ₀, is (Limousin et al. 2005):

(C.2)

In the Singular Isothermal Sphere (SIS) limit, r_core → 0 and r_cut → ∞, Eq. (C.1) reduces to the familiar expression ). In LensTool, dPIEs are implemented trough a fiducial velocity dispersion σ_LT related to the central velocity dispersion σ₀ by: .

Projecting the 3D-density profile on a plane perpendicular to the line-of-sight, we obtain the projected dPIE surface-mass-density as a function of the projected distance from the center R (Elíasdóttir et al. 2007):

(C.3)

Integrating Eq. (C.1), we can also write the total 3D mass m(r) enclosed by a sphere of radius r:

(C.4)

Similarly, the total projected mass within an aperture of projected radius R, M(R), is (from Eq. (C.3)):

(C.5)

The dPIE profile, contrary to the SIS, has a finite total mass M_tot whose value can be obtained from Eq. (C.4) (or Eq. (C.5)) in the limit of r → ∞ (R → ∞):

(C.6)

In the limit r_core/r_cut → 0, a sphere of radius r = r_cut encloses half of the total 3D dPIE mass. While 60 (90)% of the total projected mass is included within a projected radius of R = r_cut (R = 5r_cut) (see Fig. C.1).

Fig. C.1. Fraction of projected mass over total mass as a function of aperture radii (in units of rcut) for dPIE profiles with different rcore/rcut values. For small rcore/rcut, ∼60% of the total mass is contained within rcut, while the 90% is contained within 5rcut.

Similarly, Elíasdóttir et al. (2007) give an expression for the half-mass radius, R_Mtot/2, that is enclosing half of the total projected mass, from Eq. (C.5):

(C.7)

Finally, from Eq. (C.4) we derive an expression for the circular velocity :

(C.8)

which has a maximum at . In the SIS limit, Eq. (C.8) reduces to . To compare the circular velocity in Grillo et al. (2015) and B18 with our sigma functions in Fig. 6, we adopt , as one can easily verify that for the full range of r_core and r_cut values of our member galaxies, is smaller than by less than the 3%.

C.2. dPIE velocity dispersion

Assuming an isotropic mass distribution, the observed line-of-sight velocity dispersion at a projected distance R from the center is given by (Elíasdóttir et al. 2007; Agnello et al. 2014; Binney & Mamon 1982):

(C.9)

where ν(r) is the 3D-luminosity density of the galaxy which is related to the surface brightness profile, I(R), by:

(C.10)

Similarly, the average line-of-sight velocity dispersion inside an aperture radius R, σ_ap, for an isotropic system, is given by:

(C.11)

where L(R) is the total projected luminosity within R:

(C.12)

σ_ap is the quantity that we can directly compare to our kinematic measurements (see Sects. 4 and 5) and it can also be rewritten in terms of the measurable surface brightness of the elliptical galaxy as (Agnello et al. 2014):

(C.13)

The last line recasts the aperture-averaged σ_ap in terms of the central velocity dispersion σ₀ explicitly, with the projection coefficient c_p(R) encompassing all other terms except the prefactor . The relation between σ_ap and the LensTool fiducial velocity dispersion, σ_LT, can be written as:

(C.14)

To compute c_p(R), a surface brightness profile I(R) for cluster members is needed. For our purpose, we tested two choices for I(R) using: (1) a surface brightness profile scaling with the projected mass density Σ(R) (Eq. (C.3)) derived from the cluster lens model by B18, and (2) a Sérsic profile extracted from the HST data of MACS 0416 in the F814 band:

(C.15)

The Sérsic index n and the effective radius R_e are determined by fitting the light profile of cluster members using the public software GALFIT (Peng et al. 2002, 2010), adopting b_n = 2n − 1/3 (Ciotti & Bertin 1999).

Figure C.2 shows that for our spectral aperture radius, R = 0.8″, the differences between the σ_ap obtained assuming a dPIE or a Sérsic surface brightness profile are within ∼5%. For this reason, in our work we always assume a dPIE surface brightness profile which makes the analysis convenient and independent from Sérsic light profile fitting.

Fig. C.2. Ratio between aperture average line-of-sight velocity dispersions, computed assuming a dPIE or a Sersic surface brightness profile of 95 cluster members in MACS 0416, as a function of F160W magnitudes. The aperture radius is R = 0.8″. Red dashed lines correspond to the 16th, 50th and 84th percentiles of the entire distribution.

Assuming I(R) = Σ(R), the projection coefficient c_p(R) can be numerically computed as

(C.16)

In Fig. C.3, we show the projection coefficients as a function of aperture radius R, for different values of r_core and r_cut. The dashed black line indicates the asynthotic value of , which corresponds to the limit r_core → 0 and r_cut → ∞, where the dPIE reduces to a SIS for which .

Fig. C.3. Projection coefficient cp as a function of aperture radius R in arcseconds. Curves associated to the same rcore are grouped with the same color, with rcut = 5, 10, 20, 100 arcsec, from bottom to top. The thick dashed red and green lines correspond to typical values of rcore and rcut derived from our lensing+dynamics modeling of clusters. The dashed black horizontal line corresponds to .

All Tables

All Figures

Fig. 1. Results from line-of-sight velocity dispersion fitting of spectra of three cluster members in MACS 1206 as obtained with pPXF. Galaxy spectra are shown in black; red curves are the pPXF best-fit models, while the green points correspond to the data−model residuals. The blue shaded regions along the wavelength axis were excluded in the fitting procedure due to the presence of noisy sky subtractions around emission lines in the spectra. Corresponding residuals in these regions are marked in blue. The first from the top is a high signal-to-noise E+A galaxy. The second is a passive galaxy spectrum with a ⟨S/N⟩ = 28.1, corresponding approximately to the mean S/N of our galaxy sample. The bottom spectrum has a ⟨S/N⟩ = 10.7, close to our lower limit for reliable velocity dispersion measurements. Coordinates, measured velocity dispersion (σap), mean signal-to-noise and F160W magnitudes are indicated in each panel. Cutouts are HST RGB images, 4″ across, showing in green the apertures of 0.8″ radius used for the spectral extraction.

In the text

Fig. 2. Data points are measured stellar velocity dispersions of cluster members, color-coded according to their spectral ⟨S/N⟩. Green lines are 300 scaling relations randomly drawn from the posterior distributions of the σap–mag scaling relation parameters, α, , Δσap, obtained from fitting the data points (see Fig. 3). Optimized parameters are also quoted. Red and blue areas are obtained by projecting the 3D σ–L scaling relations from our previous lens models with no kinematic prior (see text). Some velocity dispersion measurements, which deviate significantly from the scaling relations, are marked with squares and circles and discussed in the text.

In the text

	Fig. 3. Posterior probability distributions for σap–mag scaling relation parameters, obtained as described in Appendix B, from velocity dispersion measurements of 58 cluster members in MACS 1206. The 16th, 50th and 84th percentiles of the marginalized distributions for the slope (α), normalization () and scatter around the scaling relation (Δσap) are quoted and shown as vertical dashed lines.
In the text

Fig. 4. Panels a–c: posterior probability distributions for the normalization of the σap − L scaling relation, that is aperture projected velocity dispersion , and truncation radius (), obtained with different lens models for clusters under study. The green contours and distributions refer to our best models which include the kinematic prior based on measured velocity dispersions of member galaxies (Fig. 2). For visualization clarity, in panel c we omit the C16 results because they are outside the chosen range of velocity dispersion, as visible in the bottom panel of Fig. 2. The reference magnitudes for each cluster are in Table 2. Panel d: posterior distribution of the reference central velocity () and reference truncation radius () in each cluster, normalized to an absolute mag , close to L* of the early-type galaxy population. The different colors correspond to the models with kinematic priors and contours refer to 1, 2, 3-σ confidence levels. Once normalized to the same luminosity, the three scaling relations for σ and rcut are consistent.

In the text

Fig. 5. Top: projected cumulative mass profiles, as a function of the projected distance from BCGs, corresponding to previous and current lens models for each cluster. The colored regions encompass the 16th and the 84th percentiles; colored solid lines are the median values. The dashed green lines correspond to the mass component associated to cluster members (i.e. sub-halos) from the new lens models with kinematic prior (16th, 50th, 84th percentiles). The multiple image projected distances from the cluster centers are marked with vertical black lines. Bottom: relative variation of the cumulative projected total masses with respect to our reference (green) model.

In the text

Fig. 6. Sub-halo central velocity dispersion function derived from our lens models for MACS 1206, MACS 0416 and AS1063 with kinematic prior based on measured velocity dispersions of cluster members. Previous models with no kinematic prior for AS1063 and MACS 0416 are shown for comparison. Data points correspond to median galaxy counts in 35 km s−1 bins, drawn from the posterior distributions of sub-halo model parameters. Errors also include Poissonian number statistics in each bin. Only galaxies within a cluster-centric radius of 0.16 × R200c are included (see Table 1). The top axis gives the total mass related to the central velocity dispersion by Eq. (10) (see text).

In the text

Fig. A.1. Top: difference between pPXF measured velocity dispersions (σ) and input velocity dispersions in simulated spectra (σin), as a function of σ. Only the galaxies with ⟨S/N⟩ > 60 are shown here. The data points are color shaded according to the spectra ⟨S/N⟩. Black triangles are the medians of the distributions in 7 different σ bins. The black polynomial line underlines the small increasing bias at low σ. Bottom: black solid (blue dashed) line shows the relative statistical error estimated from simulations (provided by pPXF), as a function of ⟨S/N⟩. The red curve corresponds to the ratio between the two; the statistical relative errors appear to be underestimated by ∼20% in our spectra.

In the text

	Fig. B.1. Posterior probability distributions for scaling relation parameters in Eq. (B.1) obtained from measured velocity dispersions of cluster members in MACS 0416 and AS1063. The 16th, 50th and 84th percentiles of the marginalized distributions for the slope (α), normalization () and scatter around the scaling relation (Δσap) are quoted and shown as dashed lines in each histogram.
In the text

	Fig. C.1. Fraction of projected mass over total mass as a function of aperture radii (in units of rcut) for dPIE profiles with different rcore/rcut values. For small rcore/rcut, ∼60% of the total mass is contained within rcut, while the 90% is contained within 5rcut.
In the text

	Fig. C.2. Ratio between aperture average line-of-sight velocity dispersions, computed assuming a dPIE or a Sersic surface brightness profile of 95 cluster members in MACS 0416, as a function of F160W magnitudes. The aperture radius is R = 0.8″. Red dashed lines correspond to the 16th, 50th and 84th percentiles of the entire distribution.
In the text

	Fig. C.3. Projection coefficient cp as a function of aperture radius R in arcseconds. Curves associated to the same rcore are grouped with the same color, with rcut = 5, 10, 20, 100 arcsec, from bottom to top. The thick dashed red and green lines correspond to typical values of rcore and rcut derived from our lensing+dynamics modeling of clusters. The dashed black horizontal line corresponds to .
In the text