EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 650 (June 2021) A&A, 650 (2021) A105 Full HTML
Free Access
Issue
A&A
Volume 650, June 2021
Article Number A105
Number of page(s) 17
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/202140564
Published online 15 June 2021

© ESO 2021

1. Introduction

Cosmological halos from Λ cold dark matter (ΛCDM) numerical simulations are known to have an internal mass distribution described by a universal shape, well represented by the Navarro-Frenk-White (NFW) model (Navarro et al. 1996, 1997). Alternative models of this universal profile have been proposed (e.g., Moore et al. 1998; Ricotti et al. 2007; Del Popolo 2010; Navarro et al. 2004; Stadel et al. 2009). The inner slope of the mass distribution might deviate from the NFW model because of several processes such as dynamical friction, central condensation of cooled gas, and active galactic nucleus (AGN) feedback (e.g., Blumenthal et al. 1986; Martizzi et al. 2012; Schaller et al. 2015; Laporte & White 2015; Peirani et al. 2017; He et al. 2020). The external slope of the mass distribution in individual clusters might deviate from the NFW model depending on the mass accretion rate in the cluster outskirts (e.g., Diemer & Kravtsov 2014; More et al. 2015; Vallés-Pérez et al. 2020). The shape of the cluster mass profile also depends on the properties of the dark matter (DM) component. Models that deviate from the classical cold DM (CDM) model generally predict mass density profiles, ρ(r), that are flatter than the NFW model near the halo center (see, e.g., Hu et al. 2000; Bode et al. 2001; Peter et al. 2010; Spergel & Steinhardt 2000; Rocha et al. 2013). An example of these profiles is the cored profile of Burkert (1995).

The two quantities that characterize most models for the mass profile of cosmological halos, M(r), are a scale parameter, the mass at a given overdensity, MΔ1, and a shape parameter, the concentration cΔ ≡ rΔ/r−2, with r−2 the radius where the logarithmic derivative of the mass density profile ρ(r), γ ≡ dlnρ/dlnr = −2. A prediction of CDM cosmological simulations is the existence of a relation between cΔ and MΔ (see, e.g., Navarro et al. 1996; Bullock et al. 2001) that depends on the halo’s redshift, z. More massive halos are less concentrated because they form later when the density of the Universe is lower (White & Rees 1978). The concentration is predicted to decrease with increasing z at a given MΔ, and the concentration–mass relation is predicted to flatten with z (e.g., Navarro et al. 1996; Bullock et al. 2001; Zhao et al. 2003; Neto et al. 2007). The earliest investigations predicted a strong z dependence of the concentration–mass relation, but more recent works indicate a very mild z dependence, with cΔ increasing by only ∼30% from z ∼ 2 to z ∼ 0 (e.g., De Boni et al. 2013; Dutton & Macciò 2014). The normalization and evolution of the concentration–mass relation depends on the cosmological model, in particular on the Hubble and density parameters h and Ωm; the dispersion of the mass fluctuation within spheres of comoving radius equal to 8 h−1 Mpc, σ8; and the dark energy equation of state parameter w (e.g., Klypin et al. 2003; Dolag et al. 2004; Macciò et al. 2008; Carlesi et al. 2012; De Boni et al. 2013; Kwan et al. 2013). The concentration–mass relation also depends on baryonic processes. Hydrodynamical simulations have shown that, at a given mass, gas cooling and star formation tend to increase the halo concentration compared to the case where only gravitational processes are considered (Fedeli 2012; Martizzi et al. 2012; Rasia et al. 2013; Cui et al. 2016). However, feedback from AGN has the opposite effect, so that the concentrations of cosmological halos turn out to be similar in DM-only and full hydrodynamical simulations (King & Mead 2011; Teyssier et al. 2011; Martizzi et al. 2012, 2013; Rasia et al. 2013; Schaller et al. 2015; Shirasaki et al. 2018).

The universal shape of the mass profiles of cosmological halos might be the result of an initial fast assembly phase (e.g., Huss et al. 1999; El-Zant 2008; Lapi & Cavaliere 2011), characterized by the chaotic mixing and violent relaxation processes (Hénon 1964; Lynden-Bell 1967). Taylor & Navarro (2001) suggested that an even more universal quantity is the pseudo-phase-space density profile, Q(r)  ≡  ρ/σ3, where σ(r) is the total velocity dispersion profile of DM particles. Numerical simulations indicate that Q(r), for cosmological halos with a wide range of masses, follows a power-law behavior with a universal slope (e.g., Taylor & Navarro 2001; Dehnen & McLaughlin 2005; Knollmann et al. 2008). An alternative formulation of the pseudo-phase density profile is given by (Dehnen & McLaughlin 2005), where σr is the radial component of σ. The Qr profile also follows a power law with radius. The lack of a particular scale radius suggests that Q(r) and Qr have a purely gravitational nature, and are more fundamental than M(r) in describing the internal dynamics of cosmological halos. The profiles Q(r) and Qr(r) have been identified with a sort of gravitational entropy, K ≡ σ2/ρ2/3 = Q−2/3 (e.g., Lapi & Cavaliere 2009), for its formal analogy with the commonly adopted definition of entropy of the hot intra-cluster medium (ICM), (e.g., Biffi et al. 2017). Very little evolution is predicted for Q(r), with its power-law slope steepening by ≲15% from z ∼ 0 to z ∼ 2 (Lapi & Cavaliere 2009).

A relation between the power-law behavior of Q(r) and the NFW shape of M(r) (Navarro et al. 1996, 1997) can be established through the Jeans equation of dynamical equilibrium (Binney & Tremaine 1987) if the logarithmic slope of ρ(r) has a linear relation with the velocity anisotropy profile

(1)

where σθ and σϕ are the two tangential components of the velocity dispersion, usually assumed to be identical in clusters of galaxies. The existence of such a linear relation was suggested by Hansen & Moore (2006) and interpreted by Hansen (2009) in terms of the relative shapes of the radial and tangential velocity distribution functions of bound particles in a halo. Studying the velocity anisotropy profile of cosmological halos is thus important for an understanding of their internal dynamics.

Constraining M(r), Q(r), and β(r) and their evolution can in principle provide useful information on how and when cosmological halos reach dynamical equilibrium, which physical processes are involved, and the nature of DM. Clusters of galaxies are observational targets of particular interest in this sense because they are expected to be the last halos to achieve virial equilibrium in the CDM scenario (White & Rees 1978), so we can hope to trace their evolution close to their formation epoch by observing them at relatively low z. The transition epoch between the early fast accretion phase and the late slow accretion phase is indeed predicted to occur at z ≲ 1 for the most massive halos (Lapi & Cavaliere 2009).

Another advantage of studying clusters of galaxies is that their M(r) can be determined over a wide range of scales in several ways, that is, via X-ray observations of the hot ICM, via the Sunyaev-Zeldovich effect (Sunyaev & Zeldovich 1970), by the gravitational lensing of galaxies in the cluster background, and by analyzing the distribution of cluster members in projected phase-space. (see Pratt et al. 2019, for a review of these methods). The last method determines M(r) from the central ∼50 kpc to very large radii (e.g., Rines & Diaferio 2006; Biviano et al. 2013), and determines the β(r) of cluster galaxies (e.g., Biviano & Katgert 2004). Knowledge of β(r) allows us to determine σ(r) and σr from the line-of-sight velocity dispersion profile of cluster galaxies, σlos. If ρ(r) is known from any of the methods mentioned above, Q(r) and Qr(r) of the DM components can then be derived, modulo an assumption about the similarity of the galaxies and DM σlos since the DM σlos (and σ and σr) are not observables.

There have been many determinations of the M(r) of clusters of galaxies using cluster galaxies as tracers of the potential, both for individual clusters (e.g., Geller et al. 1999; Rines et al. 2000, 2003; Łokas & Mamon 2003; Łokas et al. 2006; Rines & Diaferio 2006; Wojtak & Łokas 2010; Biviano et al. 2013, 2017a; Guennou et al. 2014; Munari et al. 2014; Balestra et al. 2016; Maughan et al. 2016; Sartoris et al. 2020) and for stacks of several clusters (e.g., Carlberg et al. 1997a; van der Marel et al. 2000; Biviano & Girardi 2003; Katgert et al. 2004). Results from these studies have been used to constrain the concentration–mass relation, in most cases by adopting the NFW model (Groener et al. 2016; Biviano et al. 2017b). In general, these studies confirm that the NFW model provides an acceptable fit to the cluster M(r), albeit with considerable variance from cluster to cluster, and there is a reasonable agreement between the observed concentration–mass relation and that predicted by ΛCDM cosmological numerical simulations. However, all these studies involve clusters at z < 1. There are only a few determinations of the M(r) of clusters of galaxies at z ∼ 1 or above, the one by Biviano et al. (2016, hereafter B+16) based on the phase-space distribution of cluster galaxies in the GCLASS sample (Muzzin et al. 2012), the ones by Babyk et al. (2014) and Amodeo et al. (2016) based on X-ray observations, and the one by Sereno et al. (2015) based on gravitational lensing.

Previous observational determinations of Q(r) and Qr(r) for clusters of galaxies are those of Munari et al. (2014) at z = 0.09, Biviano et al. (2013) at z = 0.44, and B+16 at z ∼ 1. All these studies found results that agree with the theoretical predictions, namely that Q(r) and Qr(r) follow power laws with the predicted slopes.

Observational determinations of β(r) for clusters of galaxies are also in general agreement with predictions from numerical simulations (e.g., Diaferio 1999; Munari et al. 2013; Lotz et al. 2019), with β(r) ≃ 0 near the cluster center and increasing outside. This means cluster galaxies are on isotropic orbits near the center and increasingly radial orbits outside (e.g., Natarajan & Kneib 1996; Mahdavi et al. 1999; Biviano & Katgert 2004; Biviano & Poggianti 2009; Lemze et al. 2009; Biviano et al. 2013; Munari et al. 2014; Annunziatella et al. 2016; Capasso et al. 2019; Mamon et al. 2019; Stark et al. 2019). There is, however, considerable variance in the shape of β(r) from cluster to cluster (Hwang & Lee 2008; Benatov et al. 2006; Aguerri et al. 2017), as also found in cluster-sized halos from cosmological simulations (Mamon et al. 2013, hereafter MBB). At z ∼ 0, determinations of cluster β(r) indicate a dichotomy in the orbits of early-type, red, and quiescent galaxies on one side, and late-type, blue, star-forming galaxies, on the other. Early-type, red, quiescent galaxies move on isotropic orbits also at large distances from the cluster center, while late-type, blue, star-forming galaxies display more radially elongated orbits (e.g., Mamon et al. 2019, and references therein). At higher-z, the orbits of red, quiescent galaxies are more similar to those of blue, star-forming galaxies, that is, more radial outside the center of the cluster (e.g., B+16, and references therein).

In this paper we use the GOGREEN spectroscopic data set (Balogh et al. 2017, 2021) complemented with the GCLASS spectroscopic data set (Muzzin et al. 2012) to probe in detail the internal dynamics of clusters of galaxies at an unprecedented high z. In particular, we determine the M(r), β(r), Q(r), and Qr(r) of a stack of 14 clusters at redshift 0.87 ≤ z ≤ 1.37 by using the MAMPOSSt method (MBB) and the Jeans equation inversion technique (Binney & Mamon 1982; Solanes & Salvador-Solé 1990). Our analysis is similar to that of B+16, which was limited to the GCLASS data set. As in B+16 we consider the subsamples of star-forming and quiescent galaxies, and in addition we consider subsamples of clusters in two bins of cluster z and M200, and in two bins of galaxy stellar mass M⋆. With respect to our previous analysis, here we increase the number of clusters (from 10 to 14) and the total number of cluster members (from 418 to 581).

The structure of this paper is the following. In Sect. 2 we describe our data set, how we identify cluster members (Sect. 2.1), the construction of the ensemble cluster from stacking (Sect. 2.2), and the completeness of the spectroscopic sample (Sect. 2.3). In Sect. 3 we describe the methodology we used to determine M(r) (Sect. 3.1), and provide the results of our analysis for the M(r) of the ensemble cluster and its subsamples. We compare the results with theoretical predictions and previous observational results (Sect. 3.2). We also compare ρ(r) with the number density and M⋆ density profiles of cluster members (Sect. 3.2). In Sect. 4 we describe the methodology by which we determine β(r) (Sect. 4.1) and provide the results of our analysis for the β(r) of the ensemble cluster and its subsamples (Sect. 4.2). In Sect. 5 we determine the Q(r) and Qr(r) of the ensemble cluster and compare them to the theoretical predictions. In Sect. 6 we discuss our results, separately for M(r) (Sect. 6.1), β(r) (Sect. 6.2), and Q(r),Qr(r) (Sect. 6.3). In Sect. 7 we provide our conclusions.

Throughout this paper we adopt the following cosmological parameter values: a Hubble constant H0 = 70 km s−1 Mpc−1, a present-day matter density Ωm = 0.3, and a curvature parameter value Ωk = 0.

2. The data set

The cluster sample studied in this work is drawn from the GOGREEN (Balogh et al. 2017, 2021) and GCLASS (Muzzin et al. 2012) surveys. The GOGREEN survey targeted 27 clusters and groups at 1.0 < z < 1.5. Three of these clusters were discovered by the South Pole Telescope (SPT) survey (Brodwin et al. 2010; Foley et al. 2011; Stalder et al. 2013), and another nine were taken from the Spitzer Adaptation of the Red-sequence Cluster Survey (SpARCS, Muzzin et al. 2009; Wilson et al. 2009; Demarco et al. 2010). Nine groups in the COSMOS and Subaru-XMM Deep Survey (Finoguenov et al. 2007, 2010; George et al. 2011) complete the GOGREEN survey. The GCLASS survey targeted ten clusters at 0.8 < z < 1.3, all taken from the SpARCS; five of them were also targeted in the GOGREEN survey.

The core observations of both GCLASS and GOGREEN were extensive multi-object spectroscopy of cluster galaxies with the GMOS spectrograph at the Gemini telescopes (Hook et al. 2004), with a spectroscopic resolving power R = 440 and with a typical individual cz uncertainty of 278 km s−1, corresponding to a cluster rest-frame uncertainty δz < 154 km s−1. This uncertainty is small compared to the σlos of the GCLASS and GOGREEN clusters, and thus fully adequate for a dynamical analysis, given that δz contributes to the observed σlos in quadrature (Harrison & Noonan 1979; Danese et al. 1980). The combined GCLASS and GOGREEN sample contains 2257 unique objects with good quality z measurements (as defined in Balogh et al. 2021): 1529 from GOGREEN and 728 from GCLASS. To this data set we add 112 unique objects with good z measurements from the literature (references are given in Table 1).

Table 1.

List of cluster properties.

In this paper we only consider the 14 clusters with 20 spectroscopic members or more (membership determination is described in Sect. 2.1). Cluster velocity distributions based on fewer than 20 members can be biased low by the effect of dynamical friction (Old et al. 2013; Saro et al. 2013), and cluster velocity dispersion estimates based on fewer than 20 members have been shown to be statistically unreliable (Girardi et al. 1993). The list of the 14 clusters and their properties are given in Table 1.

2.1. Cluster membership

To identify cluster members we proceeded in three steps. First we identified the main cluster peak in z space (following Beers et al. 1991; Girardi et al. 1993) by selecting those galaxies in the cluster field with c ∣ z − zc ∣ ≤6000 km s−1, where c is the speed of light and zc is the initial cluster redshift estimate, taken from Balogh et al. (2017) for the GOGREEN clusters and B+16 for the GCLASS clusters.

We then applied the Kernel Mixture Model (KMM) algorithm (McLachlan & Basford 1988; Ashman et al. 1994) to the distribution of redshifts located in the main peak. KMM estimates the probability that the z distribution is better represented by k Gaussians rather than a single one. Given the limited number of galaxies with redshifts available for each cluster, we only considered the case k = 2. This part of the procedure is meant to identify cases of merging subclusters close to the line of sight, and to separate these subcluster components from the main cluster. Secondary peaks are identified only in two clusters, SpARCS1051, and SpARCS1616. These secondary peaks correspond to two small groups of four galaxies each, and they are removed from the samples of cluster members. Overall, the KMM analysis suggests that the GOGREEN clusters have no strong contamination from fore- or background groups nor are they undergoing major mergers along the line of sight.

The membership selection of the galaxies that are left in the sample after the main-peak and KMM selections was then refined by using two methods: Clean (MBB) and CLUMPS, a new algorithm that we describe in detail in Appendix A. Both Clean and CLUMPS identify cluster members based on location in projected phase-space R, vrf, where R is the projected radial distance from the cluster center and is the rest-frame velocity, with the mean cluster redshift estimated on the identified cluster members (see below). For the cluster center we use the position of the brightest cluster galaxy (BCG; van der Burg et al. 2020).

These two methods for cluster membership are conceptually very different. Clean is theoretically motivated, its parameters fixed by our understanding of the properties of cluster-sized halos extracted from cosmological numerical simulations. CLUMPS is instead less model-dependent, and is based only on the fact that a cluster of galaxies is a concentration in projected phase-space. Clean ensures stability thanks to the use of theoretically motivated models, but its results depend on how close the chosen models are to reality. CLUMPS, on the other hand, depends on the choice of parameters that are not linked to a specific model, thus chosen with considerable freedom, and this can cause instability in the results. The combined use of Clean and CLUMPS helps in reducing possible systematics in the process of membership selection (Wojtak et al. 2007) as we expect the two methods to compensate each other for their relative weaknesses.

More specifically, the Clean method uses an estimate of the cluster line-of-sight velocity dispersion, σlos, to provide a first estimate of the cluster mass from a scaling relation. It then adopts the NFW profile, a theoretical concentration–mass relation (Macciò et al. 2008), and a velocity anisotropy profile model (Mamon et al. 2010), to predict σlos(R) and to iteratively reject galaxies with ∣vrf∣ > 2.7 σlos at any radius (for more details, see MBB). The adopted profiles are known to be a good description of real cluster profiles on average. Adopting different models would not make a large difference in the membership selection since they would predict similar velocity dispersion profiles (for an example, see Fig. 5), but the adopted profiles might not be a good description of the individual cluster profile. For this reason we supplement Clean with CLUMPS. CLUMPS evaluates the density of galaxies in projected phase-space, by counting galaxies in bins of R and vrf, and applies a convolution of this density distribution with a Gaussian filter in Fourier space. Then for each bin in the radial direction CLUMPS identifies the main peak in velocity space. The minima around this peak define the velocity boundaries outside which interlopers are removed in any given bin (see Appendix A for more details).

We ran each of these two methods twice. On the first run we defined as the average redshift of the galaxies that passed the main-peak and KMM selection. On the second run we re-defined as the average redshift of the galaxies selected as members in the first run.

In Col. 5 of Table 1 we list NM, the weighted number of cluster members, assigning weight = 1.0 to objects selected as members by both membership algorithms (described in the text), weight = 0.5 to objects selected as members by only one of the two membership algorithms, and weight = 0.0 to objects that are not selected as members by either of the two algorithms. In the same Col. 5 we also list in parentheses Nm, the number of galaxies identified as members by either of the two membership algorithms. In Cols. 6 and 7 we provide the cluster mean redshift and, respectively, velocity dispersion σlos, obtained on the sample of NM members with the weighted mean and weighted dispersion estimators, using the galaxy membership weights.

2.2. The ensemble cluster

We did not have a sufficient number of spectroscopically identified member galaxies per cluster (typically ≳200) to allow the determination of the mass and anisotropy profile of each cluster independently. To improve the statistics we built an ensemble cluster by stacking the data of our 14 clusters, a procedure widely used in the literature (e.g., Carlberg et al. 1997b; Biviano & Girardi 2003; Mahdavi & Geller 2004; Katgert et al. 2004; Rines et al. 2013; Cava et al. 2017). We based our stacking procedure on the theoretical expectation that the mass profile of massive halos are almost homologous. Both theory and observations indicate that, on average, the M(r) of clusters have very similar concentrations at different redshifts and masses (e.g., De Boni et al. 2013; Biviano et al. 2016), and therefore depend on a single parameter, r200. We could then build the ensemble cluster by stacking together all the members of our 14 clusters in projected phase-space, by normalizing the clustercentric distances R by r200 and the rest-frame velocities vrf by v200. We therefore needed to estimate r200 for each of our 14 clusters. To check for possible systematic errors, we provided three estimates of r200 of each individual cluster. These three estimates are listed in Col. 8 of Table 1.

The first estimate is a by-product of the Clean procedure. Clean does in fact use σlos to predict the cluster r200 in an iterative way (see Appendix B of MBB). We call this r200 estimate r200, C. The second r200 estimate, called r200, Mt, comes from application of the MAMPOSSt method (described in Sect. 3.1) to the sample of members defined by Clean and CLUMPS, using membership weights. Since the individual cluster data sets are not large enough to constrain several parameters of the mass and anisotropy profiles, we only allowed the value of r200 as a free parameter in the MAMPOSSt procedure. We assumed the light-traces-mass hypothesis, by forcing the galaxy number density profile to have the same shape as the mass density profile, and fixed the latter to the NFW profile with a concentration c200 = 3.5, a typical value for clusters at z ∼ 1 according to cosmological numerical simulations (e.g., Duffy et al. 2008). We also assumed the β(r) model of Tiret et al. (2007), β(r) = β∞ r/(r + r−2), where r−2 is the radius where the logarithmic derivative of the mass density profile equals −2, and we took β∞ = 0.5. This β(r) model was found to provide a good fit to the β(r) of simulated clusters (e.g., Mamon et al. 2010). Like the previous one, our third r200 estimate also comes from an application of the MAMPOSSt method. Here we do not use theoretical estimates of M(r) and β(r); instead, we adopt the best-fit profiles obtained by B+16 on the GCLASS sample of clusters, namely a Burkert profile (Burkert 1995) for M(r) and an Opposite profile (Biviano et al. 2013) for β(r). We call this r200 estimate r200, Mp.

The three r200 estimates are rather similar; on average, r200, C and r200, Mt are equal, r200, C and r200, Mp differ by 6%, and r200, Mt and r200, Mp differ by 6%. The maximum difference found among the three r200 estimates for any of our 14 clusters is 13%. We adopted r200 = r200, Mt to build our ensemble cluster in normalized projected phase-space, R/r200 and vrf/v200. We verified that our conclusions are not sensitive to this choice. The normalized projected phase-space distribution obtained using r200 = r200, Mt is shown in Fig. 1.

thumbnail Fig. 1.

Normalized projected phase-space diagram of the ensemble cluster obtained using r200 = r200, Mt. Large dots have membership weights of 1.0, small dots membership weights of 0.5. Escape-velocity curves are shown based on a NFW M(r) with c200 = 3.5 and with a Tiret β(r) with β∞ = 0.5.

To increase the spectroscopic completeness of the sample (important for the dynamical analysis; see Sect. 2.3) we restrict the dynamical analysis to galaxies in the ensemble clusters with stellar mass M∗ ≥ 109.5 M⊙. This is the M* limit to which the exposure times, and the associated depths, of the GOGREEN photometry were tailored (van der Burg et al. 2020). This M* selection only removes 7% of the cluster members.

Properties of this ensemble cluster (and subsamples extracted from it) are listed in Table 2. They are evaluated as weighted averages of the corresponding properties of the 14 clusters that make the ensemble cluster, using the number of members of each cluster as a weight. We considered the following subsamples of the ensemble sample in the current analysis:

  • low-M200, 7 clusters with M200 below the median M200 of the ensemble sample, M200 = 2.8 × 1014 M⊙;

  • high-M200, 7 clusters with M200 above the median M200 of the ensemble sample;

  • low-z, 7 clusters with z below the median z of the ensemble sample, z = 1.10;

  • high-z, 7 clusters with z above the median z of the ensemble sample.

Table 2.

Properties of the ensemble cluster and its subsamples.

2.3. Spectroscopic completeness

Our analysis is impacted by the different survey strategies, and particularly by which galaxies were targeted spectroscopically. Spectroscopic incompleteness is not likely to affect the observed distribution of cluster galaxy velocities since velocity segregation with galaxy mass has a very small effect in clusters (Biviano et al. 1992, 2002; Adami et al. 1998; Goto 2005). However, since the number density profiles are affected by spectroscopic incompleteness, it is essential to correct for it in the Jeans dynamical analysis. More precisely, the global value of spectroscopic incompleteness does not affect the dynamical analysis because the number density profile only enters the Jeans equation through its logarithmic derivative. What matters is the radial dependence of the spectroscopic incompleteness. If a sample has different levels of completeness at different distances from the cluster center, the inferred number density profile will differ from the intrinsic one.

It is not feasible to evaluate the spectroscopic correction as a function of R/r200 for each individual cluster in a robust way because of poor statistics. However, the spectroscopic completeness (including the dependence on stellar mass and radial distance) is very different for the GOGREEN and GCLASS surveys. Therefore, to evaluate the spectroscopic completeness of our ensemble cluster, we divided our ensemble cluster into three subsamples: one with only GOGREEN spectroscopic data, another with both GOGREEN and GCLASS spectroscopic data, and a third with only GCLASS spectroscopic data (clusters belonging respectively to the GOGREEN, GGGC, and GCLASS samples in Table 1). We determined a spectroscopic completeness correction separately for the three subsamples of the ensemble cluster.

We evaluated the spectroscopic completeness of each of the three subsamples following van der Burg et al. (2020). Specifically, in each of several radial bins, we counted the number of galaxies, Np, with photometric redshifts within a given range of the mean cluster redshift2, |Δzp|≤0.08, and the number of these galaxies with good quality redshift measurements, Ns. Photometric redshifts were evaluated with the code EAZY of Brammer et al. (2008), as described in van der Burg et al. (2020). We considered radial bins in units of r200, and only galaxies with stellar mass M* ≥ 109.5 M⊙ (see Sect. 2.2). In the same radial bins we estimated the fraction fM of the covered area that is masked because of the presence of bright stars, diffraction spikes, and reflective halos, or because of artifacts in any photometric band (see van der Burg et al. 2020, for details). The spectroscopic completeness in each radial bin was then evaluated as C = (1 − fM) Ns/Np.

We plot in Fig. 2 the spectroscopic completeness profiles for the three subsamples of the ensemble cluster, in the radial range that is adopted for the fit of an NFW model to the number density profile (see Sect. 3.1). To avoid over-interpreting completeness oscillations due to Poisson uncertainties, we smoothed these profiles with the LOWESS technique (e.g., Gebhardt et al. 1994) and adopted the smoothing curves to correct the observed number density profiles for incompleteness.

thumbnail Fig. 2.

Points with 1σ error bars: radial spectroscopic completeness profiles for galaxies with M∗ ≥ 109.5 M⊙ in three subsamples of the ensemble cluster: GOGREEN (green dots), GCLASS (blue squares), and GGGC (red crosses). The symbols representing the three subsamples have been slightly displaced along the abscissa for display purposes. Solid, dashed, and dash-dotted curves represent LOWESS smoothing of the data represented by the dots, squares, and crosses, respectively.

3. The mass profile

3.1. Methodology: MAMPOSSt

The MAMPOSSt method (MBB) performs a maximum likelihood fit of the M(r) and β(r) models to the projected phase-space distribution of cluster members. MAMPOSSt assumes a shape for the 3D velocity distribution (a Gaussian in the current implementation of the method). It has been successfully tested on simulated halos from cosmological simulations that included both dynamically relaxed and unrelaxed halos, and it has already been applied to several data sets (e.g., Biviano et al. 2013, 2017b; Guennou et al. 2014; Munari et al. 2014; Verdugo et al. 2016; Pizzuti et al. 2017; Mamon et al. 2019; Sartoris et al. 2020).

The MAMPOSSt method is based on the Jeans equation of dynamical equilibrium (Binney & Tremaine 1987; Pratt et al. 2019), and therefore it should be applied to clusters in dynamical equilibrium. We cannot exclude the possibility that some of the 14 clusters forming the ensemble cluster are not in dynamical equilibrium. However, there is no evidence for ongoing major mergers from the KMM analysis (see Sect. 2.1), and we limited our analysis to R ≤ r200 to exclude unvirialized regions as completely as possible from our analysis3. We also exclude the very inner region, R ≤ 0.05 Mpc, where the gravitational potential of the BCG becomes dominant relative to the gravitational potential of the cluster as a whole (see, e.g., Biviano & Salucci 2006)4.

MAMPOSSt assumes spherical symmetry. While clusters are not spherically symmetric systems (e.g., Chiu et al. 2018), the ensemble cluster is spherically symmetric by construction unless the clusters that compose it have been selected with a preferential orientation along the line of sight (see also Sect. 4 in van der Marel et al. 2000, for a discussion of the spherical assumption in the case of an ensemble cluster).

We chose to run MAMPOSSt in the Split mode (see Sect. 3.4 in MBB). We first determined the best fit to the number density profile, ν(r), of the ensemble cluster members, and then we determined the best fit to the distribution of cluster member velocities as a function of their radial positions. We used the projected NFW model (Bartelmann 1996) to fit the projected number density profile, N(R), with a maximum likelihood technique. The only free parameter is the scale radius rν/r200 where the logarithmic derivative of the number density profile equals −2, or, equivalently, cν = r200/rν. In the maximum likelihood fitting technique, the normalization of the fitting function is not a free parameter, as it is constrained by the requirement that the total (completeness corrected) number of observed galaxies is identical to the integral of the fitting function over the radial range of the fit (Sarazin 1980). In the fitting procedure we used the spectroscopic radial completeness curves shown in Fig. 2 to correct for incompleteness by weighting each galaxy by 1/Cj(Ri), where Cj(Ri) is the value of the completeness curve at the radial position of the galaxy, Ri, and j identifies the subsample of clusters the galaxy belongs to, GOGREEN, GGGC, or GCLASS (see Sect. 2.3 and Fig. 2). In the fit we also weighted for membership, so that the total weight is mi/Cj(Ri), where mi = 1.0 when the galaxy is classified as a cluster member by both the Clean and the CLUMPS algorithms, and mi = 0.5 when it is classified as a cluster member by only one of the two algorithms.

While not needed for MAMPOSSt, we also determined the fit of an NFW model to the stellar mass density profile, by weighting each galaxy by M⋆, i mi/Cj(Ri). We did this to compare the stellar mass concentration to the total mass concentration (see Sect. 3.2).

There are no limitations in the number of parameters that can be used to describe M(r) and β(r) in MAMPOSSt. Given the size of our data set we considered only three free parameters, two to describe M(r) and one to describe β(r). For M(r) the free parameters we considered are r200 and r−2. The ratio c200 = r200/r−2 defines the concentration, and therefore the shape of the mass profile, while r200 is related to M200 and this defines the M(r) normalization. Given that we were working on an ensemble cluster, where radii are in units of r200 and velocities in units of v200, rather than allowing r200 to be a free parameter, we also ran MAMPOSSt by fixing r200 to the weighted average value (listed in Table 2) of the 14 cluster r200 used to stack the clusters.

We considered the following M(r) models (see Eqs. (4)–(6) in Biviano et al. 2013 and Eq. (5) in Sartoris et al. 2020 for the analytic expressions of these models):

  • the Burkert model (Burkert 1995);

  • the Einasto model (Einasto 1965; Mamon et al. 2010), with m = 5 (as in Biviano et al. 2013);

  • the Hernquist model (Hernquist 1990);

  • gNFW, the generalized NFW model, which differs from the NFW model because the inner slope of the mass density profile, γ, is allowed to be different from 1. The characteristics of our data set do not allow us to constrain γ as a free parameter. We chose instead to consider four different values, γ = 0.0, 0.5, 1.0, 1.5, with γ = 1.0 corresponding to the traditional NFW model.

For β(r), the only free parameter of MAMPOSSt depends on the chosen model. We considered the following β(r) models:

  • C: Constant anisotropy with radius, β(r) = βC;

  • O: anisotropy of opposite sign at the center and at large radii (Biviano et al. 2013), β(r) = β∞ (r − r−2)/(r + r−2);

  • OM: from Osipkov (1979) and Merritt (1985), rapidly increasing radial anisotropy with radius, ;

  • T: from Tiret et al. (2007), increasing radial anisotropy with radius, β(r) = β∞ r/(r + r−2);

  • Tl: see Tiret et al. (2007) and Mamon et al. (2019), similar to the T model, but using the scale radius of the galaxy distribution rather than the scale radius of the mass distribution, β(r) = β∞ r/(r + rν).

Model C has been frequently used in previous studies because of its simplicity (e.g., Merritt 1987; van der Marel et al. 2000; Łokas & Mamon 2003), but models T and Tl appear to better represent the β(r) of cosmological halos (e.g., Mamon et al. 2010, 2013; Munari et al. 2013). Model O (respectively OM) provides the best fit to the β(r) of GCLASS clusters (B+16) (respectively EDisCS clusters; Biviano & Poggianti 2009).

In total we have 7 × 5 = 35 MAMPOSSt runs for each combination of the 7 M(r) and 5 β(r) models. We find the best-fit parameters using the NEWUOA software for unconstrained optimization (Powell 2006). Using the Bayes information criterion (BIC, Schwarz 1978) we find that the results obtained by fixing r200 to the weighted average of the 14 cluster r200 values are statistically favored compared to the results obtained by leaving r200 as a free parameter. In the following we therefore present only the results obtained with two free parameters, r−2 in units of r200 or, equivalently, c200, for M(r); and βC, rβ, or β∞, depending on the chosen β(r) model.

3.2. Results

In Fig. 3 we show N(R) of the ensemble sample corrected for spectroscopic incompleteness, and the best-fit (projected) NFW model, for all member galaxies, and for the subsamples of quiescent and star-forming galaxies. Following van der Burg et al. (2020), we identify quiescent and star-forming galaxies based on their U − V and V − J rest-frame colors, with quiescent galaxies satisfying the following criteria:

(2)

thumbnail Fig. 3.

Projected galaxy number density profile of the ensemble cluster (dots with 1σ Poissonian error bars), corrected for spectroscopic incompleteness, and best fit with the NFW model (curve). Black dots and solid curve: all member galaxies; red squares and dash-dotted curve: quiescent galaxies; blue stars and dashed curve: star-forming galaxies. Best-fit concentration values are listed.

We list in Table 3 the best-fit values of cν and their uncertainties for the ensemble cluster and all the subsamples extracted from it.

Table 3.

Concentrations of the ensemble cluster and its subsamples.

In Fig. 4 we show the results of the MAMPOSSt analysis for M(r). In particular, we display the weighted average of all 35 best-fit model M(r), using the MAMPOSSt likelihoods ℒ as weights; the best-fit NFW model; and the model with the highest likelihood, a Burkert M(r). Taken at face value, these results suggest that the preferred solutions have ρ(r) with inner slopes that are less steep than in the NFW model. A similar result was suggested by B+16 based on the GCLASS data. However, the evidence is not statistically significant since all models provide statistically acceptable fits to the data, according to the likelihood-ratio test (Meyer 1975).

thumbnail Fig. 4.

MAMPOSSt results for the M(r) of the ensemble cluster. Red curve and orange shading: ℒ-weighted average of all 35 best-fit model M(r) from the MAMPOSSt analysis, and 1σ error. Blue dashed curve: highest-ℒ among NFW models. Black dot-dashed curve: highest-ℒ model (Burkert).

In general, there is no guarantee that a maximum likelihood solution also provides a good fit to the data if the chosen models are a poor choice. To allow a direct comparison with the data we project the M(r) + β(r) solutions of MAMPOSSt on the velocity dispersion profile and find very good agreement (see Fig. 5). Different models fit almost equally well to the observed velocity dispersion profile, a confirmation that we are unable to distinguish different mass profile models with the current data set.

thumbnail Fig. 5.

Line-of-sight velocity dispersion profile of the ensemble cluster (dots and 1σ error bars, estimated using Eq. (16) in Beers et al. 1990). Also shown are the MAMPOSSt results for M(r) and β(r) projected onto the velocity dispersion profile. Red curve and orange shading: ℒ-weighted average of all 35 best-fit model M(r) from the MAMPOSSt analysis, and 1σ error (reduced χ2 = 0.6). Blue dashed curve: highest-ℒ among NFW models (reduced χ2 = 0.7). Black dot-dashed curve: highest-ℒ model (Burkert, reduced χ2 = 0.4). The last point (empty dot) represents the velocity dispersion of galaxies that were not used in the MAMPOSSt analysis because they are at R ≥ r200.

We evaluate a ℒ-weighted average c200 ≡ r200/r−2 from the 35 MAMPOSSt runs for each (sub)sample (see Table 3). The 1σ errors on c200 come from the addition of three sources of errors, δc = δc, s + δc, ν + δc, j, listed separately in Table 3. The first term, δc, s, is the statistical error on the weighted average. The second term, δc, ν, accounts for our incomplete knowledge of ν(r). It is half the difference between the two values of c200 obtained by running MAMPOSSt with either the upper or the lower limit on cν as input parameter; these limits are obtained in the NFW maximum likelihood fit to N(R). Finally, the third term, δc, j, accounts for sample variance. It is estimated by a jackknife technique (Efron 1982), by taking the rms of the c200 obtained by running MAMPOSSt Nc times on Nc − 1 clusters, eliminating one different cluster at each MAMPOSSt run. Sample variance is not a dominant source of error in our analysis since the amplitude of δc, j is similar to the amplitudes of the other two error terms.

In Fig. 6 we show the ℒ-weighted average c200 for the ensemble cluster. In the top panel of the same figure we also display three other observational values of c200 from the literature (gray diamonds). One of them is the average value estimated by B+16 for the GCLASS sample. Next to this value we also plot our new determination of c200 for the same GCLASS sample (gray dot), including new GOGREEN survey data from the GCLASS clusters. Our new determination is very close to the original one by B+16.

thumbnail Fig. 6.

The concentration of the ensemble cluster compared to previous observational results and theoretical predictions. Top panel: green horizontal line: ℒ-weighted average c200 for the ensemble cluster. Shadings represent the 1σ error bar on the c200, with contributions from the three terms, δc, s, δc, ν, and δc, j (inner dark green shading, intermediate light green shading, and external dark green shading, respectively). The two symbols with 1σ error bars labeled “GCLASS” represent the results of the dynamical analysis of GCLASS clusters, as obtained by B+16 (diamond) and as obtained in this work by including the new GOGREEN survey data for the GCLASS clusters. The symbols and 1σ error bars labelled “Lensing” and “X-ray” represent mean values from the literature for clusters in the same r200 and z range as the clusters in our sample, obtained through the analysis of gravitational lensing and X-ray data, respectively (references are given in the text). Bottom panel: green horizontal line and shadings as in the top panel. Symbols with 1σ error bars are theoretical predictions from Bhattacharya et al. (2013), Dutton & Macciò (2014), Correa et al. (2015), Child et al. (2018), De Boni et al. (2013), Ragagnin et al. (2021) for a halo of the same M200 and z as our ensemble cluster, and rescaled to the cosmological parameters used in this paper.

The other two literature values shown in Fig. 6 are obtained by considering the 20 clusters with gravitational lensing determinations of c200 and the 8 clusters with X-ray data based determinations of c200, within the z and r200 ranges of our sample. We refer to these two literature values as “Lensing” and “X-ray”, respectively. When a cluster has more than one c200 determination listed, we assign the cluster the error-weighted average of its c200 values. We then compute the weighted average of all the literature values and its error, using the inverse of the listed c200 errors as weights. The c200 values on which these two averages are computed come from Amodeo et al. (2016), Babyk et al. (2014), Jee et al. (2006, 2011), Margoniner et al. (2005), and Sereno et al. (2015) and the compilations of Sereno & Covone (2013) and Groener et al. (2016).

In the bottom panel of Fig. 6 we also show the theoretical predictions for c200 of a halo with the same M200 and z as our ensemble cluster, as obtained from numerical simulations. We consider, in particular, the theoretical predictions of De Boni et al. (2013) and Ragagnin et al. (2021), which are based on hydrodynamical simulations, and those of Bhattacharya et al. (2013), Dutton & Macciò (2014), Correa et al. (2015), Child et al. (2018), which are based on gravity-only simulations. Since these simulations adopted different cosmological parameters, we use Eq. (7) and Table 2 in Ragagnin et al. (2021) to convert the theoretical c200 values to those expected in the assumed cosmology in this paper. The value of σ8 does not enter in our dynamical analysis; we assume σ8 = 0.8. We neglect the effect of a change in the baryonic mass density, as this is sub-dominant with respect to changes in Ωm and σ8. After transformation to the same cosmological parameters, the modifications to the published theoretical c200 values range from −2 to +11% and the rms among the different values decreases by 10%, from 0.36 to 0.33. The remaining differences can result from the inclusion or not of baryonic physics in the simulations (De Boni et al. 2013) and/or from the way the parameter c200 has been measured (De Boni et al. 2013; Child et al. 2018). According to Bhattacharya et al. (2013) and Child et al. (2018) the distribution of c200 among cosmological halos has a dispersion of ≃c200/3, so the expected theoretical uncertainty for a c200 value based on a sample of Nc halos is , which defines the error bars on the theoretical values shown in Fig. 6.

The comparison with previous observational determinations show that our new determination of c200 is in good agreement with that of B+16, only slightly lower, and more precise, thanks to our larger data set. This agreement is not surprising, given that the sample of GCLASS clusters investigated by B+16 is included in our sample. Our c200 value is also in good agreement with the Lensing and X-ray literature values, and it lies between the two estimates.

Our new determination of c200 is also in excellent agreement with most theoretical predictions, in particular with the most recent ones, based on either DM-only or on hydrodynamical simulations. Theoretically predicted c200 values of cluster-size halos have a very mild dependence on baryonic processes.

We do not find any significant difference between the c200 values of the low-z and high-z subsamples, nor between the c200 values of the low-M200 and high-M200 subsamples (see Table 3 and Fig. 7). This is in line with theoretical expectations (Ragagnin et al. 2021), according to which c200 changes by < 0.2 in the z and M200 ranges of our four subsamples (see Table 2, Cols. 5 and 7); that is, the predicted changes are less than our observational uncertainties.

thumbnail Fig. 7.

Concentrations (c200) for the ensemble cluster and its subsamples (black dots with 1σ error bars). The green, magenta, red, and blue squares with dot-dashed 1σ error bars represent c⋆, cν, cν, Q, and cν, SF, respectively.

We then compare the c200 values found for the different subsamples with the concentrations obtained for the number density profiles of all, quiescent, and star-forming galaxies, cν, cν, Q, cν, SF, respectively, and with the concentration obtained for the stellar mass density profiles using all galaxies, c⋆ (see Table 3 and Fig. 7). The values of c200 are consistent within 2σ with the values of cν and also with the values of cν, Q, for all (sub)samples. However, c200 > cν, SF at > 3σ level for the ensemble cluster and the low-M200 and low-z subsamples. Cluster galaxies therefore seem, in general, to trace the mass distribution if they are not star-forming.

Previous results for clusters at lower-z have found only a mild radial dependence of the total mass per galaxy in clusters (e.g., Biviano & Girardi 2003; Lin et al. 2004; Annunziatella et al. 2014), which is not in conflict with our result, given the large uncertainties in our concentration values. B+16 found a higher stellar mass concentration compared to the total mass concentration, but we do not confirm their finding.

We do not find any significant difference between cν and c⋆ in any of our subsamples, in agreement with van der Burg et al. (2014) who did not find a significant difference between cν and c⋆ in GCLASS clusters at z ∼ 1. For the ensemble cluster and its low-z and high-M200 subsamples, cν, Q > cν, SF at > 3σ. Quiescent galaxies tend to have a more concentrated spatial distribution than star-forming galaxies, as expected from previous works (e.g., Gisler 1978; Biviano et al. 1997; Dressler et al. 2013; van der Burg et al. 2014). This difference is not found to be statistically significant in all of our subsamples, likely due to the large uncertainties.

4. The velocity anisotropy profiles

4.1. Method: The Jeans inversion technique

With MAMPOSSt it is possible to constrain M(r) and β(r) at the same time, but only within the limits allowed by the restricted set of models chosen. The choice of M(r) models is facilitated by the large number of works that have determined cluster and halo M(r) from observations and cosmological simulations (e.g., Ludlow et al. 2013; Pratt et al. 2019, and references therein). On the other hand, simulations suggest that there is a large variance in the shape of β(r) among different cluster-size halos (see Fig. 1 in Mamon et al. 2013). There is no guarantee that the four β(r) models we adopt in the MAMPOSSt analysis (see Sect. 3.1) are a close fit to the real, average cluster β(r).

We can go beyond the adopted β(r) models by performing the inversion of the Jeans equation by the method of Binney & Mamon (1982) in the implementations of Solanes & Salvador-Solé (1990) and Dejonghe & Merritt (1992). The inversion of the Jeans equation requires knowledge of M(r). We adopt the ℒ-weighted average of all 35 MAMPOSSt best-fit model M(r). The inversion procedure also requires knowledge of two observables: the projected number density profile N(R) and the line-of-sight velocity dispersion profile σlos(R). In practice, some smoothing of the observed profiles is needed to run the Jeans inversion algorithm, and we use the LOWESS smoothing technique (see, e.g., Gebhardt et al. 1994). The number density profile is then deprojected numerically (using Abel’s equation; see Binney & Tremaine 1987). Since the equations to be solved contain integrals up to infinity, we extrapolate the profiles with the functional form of Eq. (10) in Biviano et al. (2013), out to 30 Mpc; we confirmed that this radius is large enough for the results to be stable.

Uncertainties in the β(r) are estimated by performing the Jeans inversion on 100 bootstrap resamplings of the original data sets. The average of the 100 bootstrap resamplings is finally taken as our fiducial β(r) to correct for possible bias in the original estimate (Efron & Tibshirani 1986). The resulting β(r) is finally compared with the ℒ-weighted average of all 35 MAMPOSSt best-fit model β(r) to check for consistency in the procedure.

4.2. Results

We show in Fig. 8 the velocity anisotropy profile for all galaxies in the ensemble cluster (green line and shading). Rather than displaying β(r), we display (σr/σθ)(r) since this represents the relative importance of the radial and tangential components of the velocity dispersion in a more linear way. Figure 8 shows that (σr/σθ)(r) increases from ≃1 near the cluster center to slightly > 1 at r200, which means that the orbits of cluster members are isotropic near the cluster center and become more radially elongated at larger radii. The Jeans inversion result is fully consistent with the MAMPOSSt result (magenta dash-dotted curves). In Fig. 8 we also compare our new result with that found for the GCLASS sample by B+16 (gray shading). We find full consistency between the two profiles, the new one being more precise thanks to the increased statistics, now showing stronger evidence for non-isotropic orbits at large radii.

thumbnail Fig. 8.

Velocity anisotropy profile (σr/σθ)(r) for all galaxies of the ensemble cluster (pale green line) and 1σ confidence level (dark green shading). Dot-dashed magenta lines: 1σ intervals of the ℒ-weighted average of the 35 MAMPOSSt best-fit model (σr/σθ)(r). Gray shading: 1σ intervals of the result of B+16 for the GCLASS sample.

Our new β(r) determination is the highest-z so far obtained for galaxies in clusters. It is very similar to other β(r) found for galaxies in lower-z clusters (e.g., Natarajan & Kneib 1996; Biviano & Katgert 2003; Biviano & Poggianti 2009; Lemze et al. 2012; Biviano et al. 2013; Annunziatella et al. 2016; Capasso et al. 2019; Mamon et al. 2019).

We also determine the β(r) of eight subsamples of the ensemble cluster (see Fig. 9). Splitting the ensemble cluster into its subsamples reduces the statistics and increases the error bars on the profiles, so in general it is impossible to find significant differences among the subsamples. There is no significant difference between the orbits of galaxies in low-M200 and high-M200 clusters (see upper left panel of Fig. 9) nor between the orbits of galaxies in low-z and high-z clusters (see lower left panel of Fig. 9). These results are in agreement with those of Capasso et al. (2019). Theoretical expectations are for more radial orbits in more massive and higher-z cluster-sized halos (see Fig. 10 in Munari et al. 2013); we might be unable to distinguish any such trend because of the limited statistics and the limited z- and M200-range of our clusters.

thumbnail Fig. 9.

Velocity anisotropy profiles (σr/σθ)(r) for several subsamples. Confidence intervals are 1σ. Upper left panel: low-M200 (cyan line, blue shading) and high-M200 (red line, orange shading). Lower left panel: low-z (cyan line, blue shading) and high-z (red line, orange shading). Upper right panel: star-forming galaxies (cyan line, blue shading) and quiescent galaxies (red line, orange shading) in the ensemble sample. Lower right panel: 9.5 ≤ log M⋆/M⊙ < 10.5 galaxies (cyan line, blue shading) and log M⋆/M⊙ ≥ 10.5 galaxies (red line, orange shading) in the ensemble sample.

We find that star-forming galaxies have more radially elongated orbits than quiescent galaxies (see upper right panel of Fig. 9), in agreement with previous results (Biviano & Katgert 2003; Biviano & Poggianti 2009; Biviano et al. 2013; Munari et al. 2014; Mamon et al. 2019), but this difference is only marginally significant. Slightly more significant is the result that low-M⋆ (9.5 ≤ log M⋆/M⊙ < 10.5) galaxies have more radially elongated orbits than high-M⋆ (log M⋆/M⊙ ≥ 10.5) galaxies (see lower right panel of Fig. 9). Given that star-forming galaxies are, on average, less massive than quiescent ones, the two results given above are not independent. Taken at face value, the apparent difference we find between low- and high-M⋆ galaxy orbits is in the opposite sense of that found by Annunziatella et al. (2016) in a z = 0.2 cluster. However, the sample of Annunziatella et al. (2016) only considered quiescent galaxies (defined based on galaxy colors) and had a lower cut in M⋆, so a direct comparison with the current result is difficult and perhaps not very meaningful.

The orbital difference between quiescent and star-forming galaxies has generally been interpreted as being the consequence of a different time of accretion into the cluster of the two populations, with the star-forming, less-massive galaxies being recent infallers (Lotz et al. 2019). Those that we observe today as massive quiescent galaxies might have entered the cluster at an earlier phase of its evolution when collective collisions were able to isotropize galaxy orbits (Lapi & Cavaliere 2011). The lower-velocity dispersion of clusters at an earlier phase might have been more suitable to the occurrence of galaxy-galaxy mergers, which could have played a role in dissipating part of the orbital energy of these galaxies.

5. The pseudo phase-space density profiles

By combining the results of Sect. 3.2 on M(r) and Sect. 4.2 on β(r), in this section we determine the pseudo phase-space density profiles Q(r) and Qr(r). These are obtained as the ratio of either the total mass density or the number density profile to the cube of either the total (for Q(r)) or the radial (for Qr(r)) velocity dispersion profile.

We determine the total mass density profile ρ(r) from differentiation of the ℒ-weighted average of the 35 MAMPOSSt best-fit M(r) (see Sect. 3.2), and obtain its 1σ uncertainty from the ℒ-weighted dispersion of the corresponding 35 ρ(r). We determine the number density profile ν(r) by Abel inversion of the projected number density profile N(R) (see Sect. 4.1). We assume for ν(r) the same fractional uncertainties of N(R). Within the Jeans inversion procedure (see Sect. 4.1) we not only determine β(r), but also the total 3D velocity dispersion profile σ(r) and radial velocity dispersion profile σr(r) (see Eqs. (20) and (21) in Solanes & Salvador-Solé 1990), and their uncertainties via the already described bootstrap technique.

In Fig. 10 we show the ensemble cluster Q(r) and Qr(r) divided by power laws of the radius Qth ∝ r−α and , respectively. We show in green the results for α = 1.84 and αr = 1.92, the values found by Dehnen & McLaughlin (2005) based on the cosmological DM-only simulations of Diemand et al. (2004a,b). According to Lapi & Cavaliere (2009, see their Fig. 6) these values evolve with z, increasing by a factor ≃1.1 at the mean redshift of our sample. The results for these evolved values of α and αr are shown in red in Fig. 10. The power laws have been normalized to the median value of Q(r) and Qr over the radial range 0.05−1 r200. We show the results for ρ(r) in the left panels and those for ν(r) in the right panels.

thumbnail Fig. 10.

The pseudo phases-space density profiles of the ensemble cluster. The solid lines show Q(r) and Qr(r) of the ensemble cluster (top and bottom panels, respectively) derived from the total mass density profile ρ(r) and the number density profile ν(r) (left and right panels, respectively), divided by power laws of the radius Qth ∝ r−1.84 and Qr,th ∝ r−1.92 (green line and shadings) or by power laws Qth ∝ r−2.02 and Qr,th ∝ r−2.11 (red line and orange shading). The shadings indicate 1σ uncertainties. The dashed black line is unity.

We see from Fig. 10 that the observational profiles do not deviate significantly from the theoretical power laws. Given the current observational uncertainties we are unable to claim which of the considered theoretical slopes is more consistent with our data.

6. Discussion

Our dynamical analysis indicates that the velocity distribution of the members of our z ≃ 1.1 ensemble cluster (Fig. 5) can be nicely described by a spherical system of collisionless tracers in dynamical equilibrium on mildly radial anisotropic orbits (Fig. 8) in a gravitational potential that is quite close to theoretical predictions for cluster-size halos in a ΛCDM cosmology (Figs. 4 and 6).

Our results have been obtained for an ensemble cluster built by stacking 14 different clusters. Even if theory predicts that clusters form a nearly homologous set modulated mostly by rΔ and only mildly by cΔ, there is observational evidence that suggests some variance in the dynamical structure of clusters of similar mass and redshift (Donahue et al. 2002; Biviano et al. 2017a; Lopes et al. 2018). Simulations indicate that the scatter in the cluster mass density profiles has a Gaussian distribution, so the technique of stacking halos should not suffer from bias (Reed et al. 2011). Additional support for the use of a stack sample comes from the analysis of simulated halos by the MAMPOSSt technique. MBB have shown that the recovered c200 from a stack of 11 simulated clusters does not suffer from a statistically significant bias. To address the importance of the variance in the cluster mass density profiles, we ran a jackknife analysis. We found that sample variance does not dominate the error budget on the derived c200 (see Sect. 3.2 and Table 3). However, our results for the ensemble cluster must be considered valid only in an average sense.

6.1. The mass profile

The ensemble cluster M(r) is very close to the NFW model (Fig. 4). Although the Burkert (core) model is preferred over the NFW model by the MAMPOSSt analysis, the difference between the two models is not significant, and other models (Einasto, Hernquist, gNFW) are in fact equally acceptable. It therefore looks like NFW is an acceptable description of the mass distribution in clusters of galaxies from z ∼ 0 (e.g., Carlberg et al. 1997c; Katgert et al. 2004) to z ≃ 1.1, although other models could be acceptable as well.

The ensemble cluster c200 is in very good agreement with, albeit slightly smaller than, the corresponding value obtained by B+16 for the GCLASS clusters by a similar analysis. Our new c200 value is also in good agreement with previous estimates for clusters in the same z and r200 range of our 14 clusters, and based on gravitational lensing and X-ray data. It has been argued that c200 estimated from gravitational lensing (Meneghetti et al. 2011, 2014) and X-ray (Rasia et al. 2013) data might be biased high because of selection effects. More concentrated clusters, and clusters with their major axis aligned along the line of sight, are stronger gravitational lenses and have higher X-ray luminosity per given mass. The similarity of the c200 values obtained via three different methods suggests that the current sample of z ≳ 1 clusters for which c200 have been measured is representative of the whole population in the given r200 range.

Our new c200 value is in excellent agreement with several predictions from cosmological simulations for a halo with the same mass and redshift as our ensemble cluster (see Fig. 6, Tables 2 and 3). The agreement is particularly good with the most recent theoretical predictions, based on either DM-only or hydrodynamical simulations. We note that there is very little difference among the c200 values predicted by DM-only (Child et al. 2018) or by hydrodynamical simulations, indipendently of whether the latter account (Ragagnin et al. 2021) or not (De Boni et al. 2013) for feedback by an active galactic nucleus (AGN) in the central cluster region. The addition of baryonic effects to DM-only simulations can, in principle, lead to a substantial change in the halo c200 (Fedeli 2012; Cui et al. 2016). However, adiabatic contraction and stellar and AGN feedback have the opposite effect on c200 (Teyssier et al. 2011), so the difference between a cluster-size halo c200 in DM-only and full hydrodynamical simulations turns out to be quite small King & Mead (2011), Rasia et al. (2013).

Given the current observational uncertainties on c200 (∼20%; see Table 3) we cannot use it to constrain the relative importance of different feedback schemes. Since the effects of AGN feedback are strongest near the cluster center (e.g., Martizzi et al. 2013; Schaller et al. 2015; Peirani et al. 2017; Shirasaki et al. 2018), better constraints can be obtained by the determination of the inner slope of the total mass density profile. As already pointed out, our ensemble cluster M(r) is better represented by a cored than a cuspy model (Fig. 4). Taken at face value, this result would suggest that AGN feedback is required. However, we cannot make the difference between a cored and a cuspy M(r) in a statistically significant way because we lack information on the kinematics in the inner ∼50 kpc that could be provided by observation of the stellar velocity distribution of the BCG (Sand et al. 2004; Newman et al. 2013; Sartoris et al. 2020).

The concentration of the total mass profile of our ensemble cluster c200 is not significantly different from either the concentration of the number density profile cν or from the concentration of the M⋆ density profile c⋆. B+16 claimed c200 < c⋆ for the sample of GCLASS clusters and attributed this to an excess of massive galaxies near the center of z ∼ 1 clusters. Our low-z subsample has z = 0.96, similar to that of the GCLASS sample analyzed by B+16, yet we do not see any difference between c200 and c⋆ in this subsample either. Therefore, the difference with respect to B+16 cannot be attributed to a difference in z between our sample and theirs. However, the significance of the claimed c200 versus c⋆ difference in B+16 is only ≲2σ, and our values of c200 and c⋆ are consistent with the respective values of B+16 within ≤1.3σ. So we conclude that there is no real tension between our results and those of B+16, in terms of relative total and stellar mass concentrations.

The similarity of the number density, stellar mass, and total mass concentrations can be rephrased by saying that the distribution of cluster galaxies appears to follow that of the total mass (or vice versa). Lower-z cluster studies have found only a mild difference in the two distributions (e.g., Biviano & Girardi 2003; Lin et al. 2004; Annunziatella et al. 2014), so our result suggests that galaxies retain a similar distribution to that of the mass from z ∼ 0 to z ≃ 1.1. This similarity is broken when only star-forming galaxies are considered. We find cν, Q > cν, SF, as also seen in lower-z clusters (e.g., Gisler 1978; Biviano et al. 1997; Dressler et al. 2013; van der Burg et al. 2014), an indication that the spatial segregation of quiescent and star-forming galaxies is already in place at z ≳ 1.

6.2. The velocity anisotropy profile

The ensemble cluster β(r) is fully consistent (but with smaller error bars thanks to the better statistics) with the result obtained by B+16 for the GCLASS sample. Orbits are isotropic near the center, and β(r) > 0 (σr/σθ > 1) at r > 0.4 r200. The evidence for non-isotropic, radial orbits, at r ≃ r200 is partially significant.

Our ensemble cluster β(r) is very similar to that of DM particles in cluster-size cosmological halos (see Fig. 8 and Mamon et al. 2010; Munari et al. 2013). The central isotropy (β(r) ≃ 0 or σr/σθ ≃ 1) is predicted to be the result of collective collisions taking place in the initial phases of the cluster formation, with the radial extent of the isotropic region being dependent on the efficiency of these collective collisions (Lapi & Cavaliere 2011). At large radii, where cluster growth is mostly governed by smooth accretion, the radial accretion of material is reflected in β(r) > 0 (Lapi & Cavaliere 2011).

Different subsamples of the ensemble cluster do not show significantly different β(r) (see left panels of Fig. 9). On the other hand, different classes of cluster galaxies do show partially significant different β(r) (see right panels of Fig. 9). The orbits of star-forming and low-M⋆ (9.5 < log M⋆/M⊙ < 10.5) galaxies are more radially elongated than those of quiescent and high-M⋆ (log M⋆/M⊙ ≥ 10.5) galaxies. A similar difference is seen in most (but not all; see Aguerri et al. 2017) lower-z clusters considering galaxies of different morphological (Mamon et al. 2019) and spectral type (Biviano & Katgert 2003; Biviano & Poggianti 2009; Biviano et al. 2013), and different colors (Munari et al. 2014).

A possible interpretation of the orbital difference between quiescent and star-forming galaxies relies on a different time of accretion into the cluster of the two populations, star-forming less-massive galaxies being recent infallers (Lotz et al. 2019). At an earlier time the cluster gravitational potential was rapidly changing, and the cluster velocity dispersion was sufficiently small to make galaxy-galaxy mergers possible. So galaxies residing in clusters for a long time could have had their orbits isotropized by collective collisions and mergers (Lapi & Cavaliere 2011). These galaxies would then evolve into massive quiescent galaxies.

6.3. The pseudo phase-space density profile

The pseudo-phase space density profiles, Q(r) and Qr(r), of the ensemble cluster are fully consistent with the predicted power-law behavior for cluster-sized halos (Taylor & Navarro 2001; Rasia et al. 2004; Dehnen & McLaughlin 2005; Lapi & Cavaliere 2009, see Fig. 10), both when using ρ(r) and when using ν(r) in the definition of Q(r) and Qr(r). Our result confirms the previous results obtained for clusters in the range 0 ≲ z ≲ 0.9 (Biviano et al. 2013; Munari et al. 2014; Aguerri et al. 2017; Capasso et al. 2019) and for the z ∼ 1 GCLASS cluster sample (B+16).

The scale-free behavior of Q ∝ rα and , and even the value of α, were predicted by the Bertschinger (1985) self-similar solution for secondary infall of material onto a point-mass perturber in a uniformly expanding universe. Cluster assembly is thought to proceed in a more chaotic way, however, so other explanations are needed for the power-law behavior of Q(r) and Qr(r). According to Dehnen & McLaughlin (2005), the scale-free behavior of Qr(r) emerges from an initial scale-free phase-space density that keeps its scale-free behavior after violent relaxation because this process is driven by gravity alone. The system evolves to dynamical equilibrium (as described by the Jeans equation) αr = 1.94 + 0.22 β(0), thus for β(0) ≃ 0, which is what we find (see Fig. 8), αr is very close to the value we observe in our ensemble cluster (see Fig. 10).

The analytical solution found by Dehnen & McLaughlin (2005) for Qr(r) is based on the assumption of a linear relation between γ(r) (the logarithmic slope of ρ(r)) and β(r), as indicated by the analysis of simulated halos (Hansen & Moore 2006). In Fig. 11 we show that β(r) and γ(r) of our ensemble cluster are indeed consistent with a linear relation, although slightly below that predicted by Hansen & Moore (2006)5.

thumbnail Fig. 11.

β − γ relation of Hansen & Moore (2006) with its ±1σ confidence interval (dashed red line and orange shading). The green shading indicates the ±1σ confidence interval for our ensemble cluster.

On average, our results for M(r), β(r), Q(r), and Qr(r) indicate that the clusters in our sample have already reached dynamical equilibrium. The clusters in our sample are not representative of the whole cluster population at z ≃ 1.1 since they are selected based on the presence of a red-sequence galaxy population (the SpARCS clusters), and a hot intra-cluster gas (the SPT clusters), and both characteristics are typical of mature galaxy systems. In this respect we conclude that at least those z ≃ 1.1 clusters that contain a hot intra-cluster medium and an evolved galaxy population are also dynamically relaxed.

The highest-z clusters known are at z ≲ 2 (Stanford et al. 2012; Andreon et al. 2014) and the z > 2 regime is generally considered to be the realm of proto-clusters (Overzier 2016). In the hypothesis that, on average, massive clusters start assembling at z ∼ 2, our results indicate that the ∼2 Gyr between z ∼ 2 and z ∼ 1 is sufficient time for clusters to reach dynamical equilibrium after initial assembly. This time is similar to the dynamical time at z = 1.1, tdyn ≃ (GΔρc(z))−1/2 = 1.5 Gyr (Sarazin 1986), for an overdensity Δ = 200, where ρc(z) is the critical density at redshift z. The origin of the universal shape of M(r), Q(r), and Qr could therefore result from a process that ensures fast dynamical relaxation of clusters at their formation, such as violent relaxation and chaotic mixing (Hénon 1964; Lynden-Bell 1967; Lapi & Cavaliere 2011; Beraldo e Silva et al. 2019). The orbital shape of cluster galaxies also seems to be largely imprinted after this violent relaxation process since the β(r) we find for our ensemble cluster is quite similar to those found for lower-z clusters, at least down to z ≃ 0.2 (Annunziatella et al. 2016; Capasso et al. 2019).

Our finding that z ≃ 1.1 clusters are dynamically old appears to fit well in the scenario supported by the analyses of the cluster galaxy population based on GOGREEN survey data. These studies (Old et al. 2020, 2021; van der Burg et al. 2020; Webb et al. 2020) have in fact shown that the GOGREEN clusters are also old in terms of the age of their member galaxies.

7. Conclusions

We investigated the internal dynamics of 14 clusters at 0.9 ≤ z ≤ 1.4 drawn from the GOGREEN spectroscopic data set (Balogh et al. 2017, 2021) complemented with data from the GCLASS survey (Muzzin et al. 2012). We stacked the 14 clusters to build an ensemble cluster containing 581 member galaxies with M⋆ ≥ 109.5 M⊙. We used MAMPOSSt (MBB) to determine the ensemble cluster mass profile M(r), and we inverted the Jeans equation with the method of Solanes & Salvador-Solé (1990) to determine the ensemble cluster velocity anisotropy profile β(r). Using the results of the MAMPOSSt and Jeans inversion analyses, we determined the pseudo phase-space density profiles Q(r) and Qr(r). We also considered four subsamples of the ensemble cluster by separating the 14 clusters that compose it into two groups of two, split by M200 or by z. Our results are the following:

  • Several M(r) models are acceptable and we cannot discriminate between cored and cuspy models, presumably because we lack dynamical tracers in the very inner region (< 0.05 Mpc) that could be provided by observations of the BCG stellar kinematics (as in Sartoris et al. 2020).

  • The concentration c200 of the ensemble cluster mass profile is close to theoretical predictions, and to previous observational determinations from the literature for clusters at similar redshifts, obtained using gravitational lensing and X-ray data.

  • The value of c200 does not depend on the cluster redshift or mass, as expected from theory, which predicts a very small change in c200 over the redshift and mass ranges spanned by our cluster sample.

  • The total mass concentration is not significantly different from the concentration of the spatial distribution of galaxies, or from the stellar mass concentration.

  • The star-forming galaxies have a less concentrated distribution than the quiescent galaxies.

  • The orbits of cluster galaxies are isotropic near the center and more radial outside.

  • High-M⋆ and quiescent galaxies have more isotropic orbits than, respectively, low-M⋆ and star-forming galaxies, but this orbital difference is only marginally significant.

  • The profiles Q(r) and Qr(r), determined either using the total mass or the number density profile, are very close to theoretical power-law predictions.

We conclude that the internal dynamics of clusters at the highest redshift probed so far in detail, do not differ from the internal dynamics of lower-redshift clusters, and confirm theoretical predictions. The fundamental dynamical properties of clusters, such as the shape of their mass density profile, the orbits of their galaxies, and the power-law behavior of the pseudo phase-space density profile have been in place since z ≲ 1.4. Given that very few clusters are known to exist beyond z ∼ 2, clusters must reach their dynamical equilibrium configuration in ∼2 Gyr, presumably by some rapid process of dynamical relaxation (Lynden-Bell 1967; Lapi & Cavaliere 2011).

Our work is based on a limited sample of 14 clusters, selected because of the presence of a red sequence of cluster members, or of a hot intra-cluster medium. The similarity of our c200 value to those determined using different methods suggests that we are not dealing with a biased cluster sample selection. Nevertheless, we caution that our conclusions might not be applicable to the whole population of z ≃ 1.1 clusters, and future investigations are needed to confirm our results with a statistically representative sample.

1

The mass MΔ is the mass contained within a sphere of radius rΔ within which the mean mass overdensity is Δ times the critical density at the cluster’s redshift. In this paper we adopt Δ = 200, so we have , where Hz is the Hubble parameter at the cluster’s redshift z. The circular velocity at the radius r200 is v200 = 10 Hz r200.

2

Changing the Δzp range from 0.08, the value used in van der Burg et al. (2020), to 0.12 does not modify the completeness curves in a significant way.

3

We checked that further restricting the radial range to R ≤ 0.8 r200 (corresponding roughly to an overdensity Δ = 300) does not change our M(r) in a significant way.

4

It is possible to extend MAMPOSSt by adding the BCG stellar kinematics as a constraint (Sartoris et al. 2020), but a precise determination of the velocity dispersion profile of the BCG is required, and it is not available for any of our 14 clusters.

5

It has been argued by Arora & Williams (2020) that the scale-free behavior of Q(r) is a mere coincidence that results from the Einasto shape of M(r) (we do find that the Einasto model is a good representation of our ensemble cluster M(r); see Sect. 3.2) and the linear β − γ relation. However, we do not see why the Einasto model or the linear β − γ relation should be considered more fundamental than the scale-free behavior of Q(r).

Acknowledgments

This work was enabled by observations made from the Gemini North, Subaru and CFHT telescopes, located within the Maunakea Science Reserve and adjacent to the summit of Maunakea. We are grateful for the privilege of observing the Universe from a place that is unique in both its astronomical quality and its cultural significance. We thank the referee for her/his useful comments that helped improving the scientific content of this paper. RD gratefully acknowledges support from the Chilean Centro de Excelencia en Astrofísica y Tecnologías Afines (CATA) BASAL grant AFB-170002. JN received support from Universidad Andrés Bello research grant DI-12-19/R. LJO acknowledges the support of a European Space Agency (ESA) Research Fellowship. GR acknowledges support from the National Science Foundation grants AST-1517815, AST-1716690, and AST-1814159 and NASA HST grant AR-14310. GR also acknowledges the support of an ESO visiting science fellowship. BV acknowledges financial contribution from the grant PRIN MIUR 2017 n.20173ML3WW_001 (PI Cimatti) and from the INAF main-stream funding programme (PI Vulcani). MLB acknowledges support from the National Science and Engineering Research Council (NSERC) Discovery Grant program. GW gratefully acknowledges support from the National Science Foundation through grant AST-1517863, from HST program number GO-15294, and from grant number 80NSSC17K0019 issued through the NASA Astrophysics Data Analysis Program (ADAP). Support for program number GO-15294 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. We thank the International Space Science Institute (ISSI) for providing financial support and a meeting facility that inspired insightful discussions for team “COSWEB: The Cosmic Web and Galaxy Evolution”. The Millennium Simulations databases used in this paper and the web application providing online access to them were constructed as part of the activities of the German Astrophysical Virtual Observatory. For this work we made use of the softwares Numpy (Harris et al. 2020) and Scipy (Virtanen et al. 2020).

References

  1. Adami, C., Biviano, A., & Mazure, A. 1998, A&A, 331, 439 [NASA ADS] [Google Scholar]
  2. Aguerri, J. A. L., Agulli, I., Diaferio, A., & Dalla Vecchia, C. 2017, MNRAS, 468, 364 [NASA ADS] [CrossRef] [Google Scholar]
  3. Amodeo, S., Ettori, S., Capasso, R., & Sereno, M. 2016, A&A, 590, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Andreon, S., Newman, A. B., Trinchieri, G., et al. 2014, A&A, 565, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Annunziatella, M., Biviano, A., Mercurio, A., et al. 2014, A&A, 571, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Annunziatella, M., Mercurio, A., Biviano, A., et al. 2016, A&A, 585, A160 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Arora, A., & Williams, L. L. R. 2020, ApJ, 893, 53 [CrossRef] [Google Scholar]
  8. Ashman, K. M., Bird, C. M., & Zepf, S. E. 1994, AJ, 108, 2348 [NASA ADS] [CrossRef] [Google Scholar]
  9. Babyk, I. V., Del Popolo, A., & Vavilova, I. B. 2014, Astron. Rep., 58, 587 [CrossRef] [Google Scholar]
  10. Balestra, I., Mercurio, A., Sartoris, B., et al. 2016, ApJS, 224, 33 [Google Scholar]
  11. Balogh, M. L., Gilbank, D. G., Muzzin, A., et al. 2017, MNRAS, 470, 4168 [NASA ADS] [CrossRef] [Google Scholar]
  12. Balogh, M. L., van der Burg, R. F. J., Muzzin, A., et al. 2021, MNRAS, 500, 358 [CrossRef] [Google Scholar]
  13. Bartelmann, M. 1996, A&A, 313, 697 [NASA ADS] [Google Scholar]
  14. Beers, T. C., Flynn, K., & Gebhardt, K. 1990, AJ, 100, 32 [NASA ADS] [CrossRef] [Google Scholar]
  15. Beers, T. C., Gebhardt, K., Forman, W., Huchra, J. P., & Jones, C. 1991, AJ, 102, 1581 [NASA ADS] [CrossRef] [Google Scholar]
  16. Benatov, L., Rines, K., Natarajan, P., Kravtsov, A., & Nagai, D. 2006, MNRAS, 370, 427 [NASA ADS] [CrossRef] [Google Scholar]
  17. Beraldo e Silva, L., de Siqueira Pedra, W., & Valluri, M. 2019, ApJ, 872, 20 [CrossRef] [Google Scholar]
  18. Bertschinger, E. 1985, ApJS, 58, 39 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bhattacharya, S., Habib, S., Heitmann, K., & Vikhlinin, A. 2013, ApJ, 766, 32 [NASA ADS] [CrossRef] [Google Scholar]
  20. Biffi, V., Planelles, S., Borgani, S., et al. 2017, MNRAS, 468, 531 [NASA ADS] [CrossRef] [Google Scholar]
  21. Binney, J., & Mamon, G. A. 1982, MNRAS, 200, 361 [NASA ADS] [CrossRef] [Google Scholar]
  22. Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton: Princeton University Press), 747 [Google Scholar]
  23. Biviano, A., & Girardi, M. 2003, ApJ, 585, 205 [NASA ADS] [CrossRef] [Google Scholar]
  24. Biviano, A., & Katgert, P. 2003, Ap&SS, 285, 25 [NASA ADS] [CrossRef] [Google Scholar]
  25. Biviano, A., & Katgert, P. 2004, A&A, 424, 779 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Biviano, A., & Salucci, P. 2006, A&A, 452, 75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Biviano, A., & Poggianti, B. M. 2009, A&A, 501, 419 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Biviano, A., Girardi, M., Giuricin, G., Mardirossian, F., & Mezzetti, M. 1992, ApJ, 396, 35 [CrossRef] [Google Scholar]
  29. Biviano, A., Katgert, P., Mazure, A., et al. 1997, A&A, 321, 84 [NASA ADS] [Google Scholar]
  30. Biviano, A., Katgert, P., Thomas, T., & Adami, C. 2002, A&A, 387, 8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Biviano, A., Rosati, P., Balestra, I., et al. 2013, A&A, 558, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Biviano, A., van der Burg, R. F. J., Muzzin, A., et al. 2016, A&A, 594, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Biviano, A., Popesso, P., Dietrich, J. P., et al. 2017a, A&A, 602, A20 [CrossRef] [EDP Sciences] [Google Scholar]
  34. Biviano, A., Moretti, A., Paccagnella, A., et al. 2017b, A&A, 607, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Blumenthal, G. R., Faber, S. M., Flores, R., & Primack, J. R. 1986, ApJ, 301, 27 [Google Scholar]
  36. Bode, P., Ostriker, J. P., & Turok, N. 2001, ApJ, 556, 93 [Google Scholar]
  37. Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503 [NASA ADS] [CrossRef] [Google Scholar]
  38. Brodwin, M., Ruel, J., Ade, P. A. R., et al. 2010, ApJ, 721, 90 [NASA ADS] [CrossRef] [Google Scholar]
  39. Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559 [Google Scholar]
  40. Burkert, A. 1995, ApJ, 447, L25 [NASA ADS] [CrossRef] [Google Scholar]
  41. Capasso, R., Saro, A., Mohr, J. J., et al. 2019, MNRAS, 482, 1043 [NASA ADS] [CrossRef] [Google Scholar]
  42. Carlberg, R. G., Yee, H. K. C., & Ellingson, E. 1997a, ApJ, 478, 462 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Carlberg, R. G., Yee, H. K. C., Ellingson, E., et al. 1997b, ApJ, 476, L7 [NASA ADS] [CrossRef] [Google Scholar]
  44. Carlberg, R. G., Yee, H. K. C., Ellingson, E., et al. 1997c, ApJ, 485, L13 [NASA ADS] [CrossRef] [Google Scholar]
  45. Carlesi, E., Knebe, A., Yepes, G., et al. 2012, MNRAS, 424, 699 [NASA ADS] [CrossRef] [Google Scholar]
  46. Cava, A., Biviano, A., Mamon, G. A., et al. 2017, A&A, 606, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Child, H. L., Habib, S., Heitmann, K., et al. 2018, ApJ, 859, 55 [NASA ADS] [CrossRef] [Google Scholar]
  48. Chiu, I. N., Umetsu, K., Sereno, M., et al. 2018, ApJ, 860, 126 [NASA ADS] [CrossRef] [Google Scholar]
  49. Correa, C. A., Wyithe, J. S. B., Schaye, J., & Duffy, A. R. 2015, MNRAS, 452, 1217 [NASA ADS] [CrossRef] [Google Scholar]
  50. Cui, W., Power, C., Knebe, A., et al. 2016, MNRAS, 458, 4052 [CrossRef] [Google Scholar]
  51. Danese, L., de Zotti, G., & di Tullio, G. 1980, A&A, 82, 322 [NASA ADS] [Google Scholar]
  52. De Boni, C., Ettori, S., Dolag, K., & Moscardini, L. 2013, MNRAS, 428, 2921 [NASA ADS] [CrossRef] [Google Scholar]
  53. Dehnen, W., & McLaughlin, D. E. 2005, MNRAS, 363, 1057 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  54. Dejonghe, H., & Merritt, D. 1992, ApJ, 391, 531 [NASA ADS] [CrossRef] [Google Scholar]
  55. Del Popolo, A. 2010, MNRAS, 408, 1808 [NASA ADS] [CrossRef] [Google Scholar]
  56. Demarco, R., Wilson, G., Muzzin, A., et al. 2010, ApJ, 711, 1185 [NASA ADS] [CrossRef] [Google Scholar]
  57. Diaferio, A. 1999, MNRAS, 309, 610 [NASA ADS] [CrossRef] [Google Scholar]
  58. Diemand, J., Moore, B., & Stadel, J. 2004a, MNRAS, 353, 624 [NASA ADS] [CrossRef] [Google Scholar]
  59. Diemand, J., Moore, B., & Stadel, J. 2004b, MNRAS, 352, 535 [NASA ADS] [CrossRef] [Google Scholar]
  60. Diemer, B., & Kravtsov, A. V. 2014, ApJ, 789, 1 [NASA ADS] [CrossRef] [Google Scholar]
  61. Dolag, K., Bartelmann, M., Perrotta, F., et al. 2004, A&A, 416, 853 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Donahue, M., Scharf, C. A., Mack, J., et al. 2002, ApJ, 569, 689 [NASA ADS] [CrossRef] [Google Scholar]
  63. Dressler, A., Oemler Augustus, J., Poggianti, B. M., et al. 2013, ApJ, 770, 62 [NASA ADS] [CrossRef] [Google Scholar]
  64. Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008, MNRAS, 390, L64 [NASA ADS] [CrossRef] [Google Scholar]
  65. Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359 [Google Scholar]
  66. Efron, B. 1982, The Jackknife, the Bootstrap and Other Resampling Plans (Philadelphia: Society for Industrial and Applied Mathematics (SIAM)) [Google Scholar]
  67. Efron, B., & Tibshirani, R. 1986, Stat. Sci., 1, 54 [NASA ADS] [CrossRef] [Google Scholar]
  68. Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87 [NASA ADS] [Google Scholar]
  69. El-Zant, A. A. 2008, ApJ, 681, 1058 [NASA ADS] [CrossRef] [Google Scholar]
  70. Fedeli, C. 2012, MNRAS, 424, 1244 [NASA ADS] [CrossRef] [Google Scholar]
  71. Finoguenov, A., Guzzo, L., Hasinger, G., et al. 2007, ApJS, 172, 182 [NASA ADS] [CrossRef] [Google Scholar]
  72. Finoguenov, A., Watson, M. G., Tanaka, M., et al. 2010, MNRAS, 403, 2063 [NASA ADS] [CrossRef] [Google Scholar]
  73. Foley, R. J., Andersson, K., Bazin, G., et al. 2011, ApJ, 731, 86 [NASA ADS] [CrossRef] [Google Scholar]
  74. Gebhardt, K., Pryor, C., Williams, T. B., & Hesser, J. E. 1994, AJ, 107, 2067 [NASA ADS] [CrossRef] [Google Scholar]
  75. Geller, M. J., Diaferio, A., & Kurtz, M. J. 1999, ApJ, 517, L23 [NASA ADS] [CrossRef] [Google Scholar]
  76. George, M. R., Leauthaud, A., Bundy, K., et al. 2011, ApJ, 742, 125 [Google Scholar]
  77. Girardi, M., Biviano, A., Giuricin, G., Mardirossian, F., & Mezzetti, M. 1993, ApJ, 404, 38 [NASA ADS] [CrossRef] [Google Scholar]
  78. Gisler, G. R. 1978, MNRAS, 183, 633 [NASA ADS] [CrossRef] [Google Scholar]
  79. Goto, T. 2005, MNRAS, 359, 1415 [NASA ADS] [CrossRef] [Google Scholar]
  80. Groener, A. M., Goldberg, D. M., & Sereno, M. 2016, MNRAS, 455, 892 [NASA ADS] [CrossRef] [Google Scholar]
  81. Guennou, L., Biviano, A., Adami, C., et al. 2014, A&A, 566, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Hansen, S. H. 2009, ApJ, 694, 1250 [NASA ADS] [CrossRef] [Google Scholar]
  83. Hansen, S. H., & Moore, B. 2006, New Astron., 11, 333 [NASA ADS] [CrossRef] [Google Scholar]
  84. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [CrossRef] [PubMed] [Google Scholar]
  85. Harrison, E. R., & Noonan, T. W. 1979, ApJ, 232, 18 [NASA ADS] [CrossRef] [Google Scholar]
  86. He, Q., Li, H., Li, R., et al. 2020, MNRAS, 496, 4717 [CrossRef] [Google Scholar]
  87. Hénon, M. 1964, Ann. Astrophys., 27, 83 [Google Scholar]
  88. Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663 [NASA ADS] [CrossRef] [Google Scholar]
  89. Hernquist, L. 1990, ApJ, 356, 359 [Google Scholar]
  90. Hook, I. M., Jørgensen, I., Allington-Smith, J. R., et al. 2004, PASP, 116, 425 [NASA ADS] [CrossRef] [Google Scholar]
  91. Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158 [Google Scholar]
  92. Huss, A., Jain, B., & Steinmetz, M. 1999, MNRAS, 308, 1011 [NASA ADS] [CrossRef] [Google Scholar]
  93. Hwang, H. S., & Lee, M. G. 2008, ApJ, 676, 218 [NASA ADS] [CrossRef] [Google Scholar]
  94. Jee, M. J., White, R. L., Ford, H. C., et al. 2006, ApJ, 642, 720 [NASA ADS] [CrossRef] [Google Scholar]
  95. Jee, M. J., Dawson, K. S., Hoekstra, H., et al. 2011, ApJ, 737, 59 [NASA ADS] [CrossRef] [Google Scholar]
  96. Katgert, P., Biviano, A., & Mazure, A. 2004, ApJ, 600, 657 [NASA ADS] [CrossRef] [Google Scholar]
  97. King, L. J., & Mead, J. M. G. 2011, MNRAS, 416, 2539 [NASA ADS] [CrossRef] [Google Scholar]
  98. Klypin, A., Macciò, A. V., Mainini, R., & Bonometto, S. A. 2003, ApJ, 599, 31 [NASA ADS] [CrossRef] [Google Scholar]
  99. Knollmann, S. R., Knebe, A., & Hoffman, Y. 2008, MNRAS, 391, 559 [NASA ADS] [CrossRef] [Google Scholar]
  100. Kwan, J., Bhattacharya, S., Heitmann, K., & Habib, S. 2013, ApJ, 768, 123 [NASA ADS] [CrossRef] [Google Scholar]
  101. Lapi, A., & Cavaliere, A. 2009, ApJ, 692, 174 [NASA ADS] [CrossRef] [Google Scholar]
  102. Lapi, A., & Cavaliere, A. 2011, ApJ, 743, 127 [Google Scholar]
  103. Laporte, C. F. P., & White, S. D. M. 2015, MNRAS, 451, 1177 [NASA ADS] [CrossRef] [Google Scholar]
  104. Lemson, G., & the Virgo Consortium 2006, ArXiv e-prints [arXiv:astro-ph/0608019] [Google Scholar]
  105. Lemze, D., Broadhurst, T., Rephaeli, Y., Barkana, R., & Umetsu, K. 2009, ApJ, 701, 1336 [NASA ADS] [CrossRef] [Google Scholar]
  106. Lemze, D., Wagner, R., Rephaeli, Y., et al. 2012, ApJ, 752, 141 [NASA ADS] [CrossRef] [Google Scholar]
  107. Lin, Y.-T., Mohr, J. J., & Stanford, S. A. 2004, ApJ, 610, 745 [NASA ADS] [CrossRef] [Google Scholar]
  108. Łokas, E. L., & Mamon, G. A. 2003, MNRAS, 343, 401 [NASA ADS] [CrossRef] [Google Scholar]
  109. Łokas, E. L., Wojtak, R., Gottlöber, S., Mamon, G. A., & Prada, F. 2006, MNRAS, 367, 1463 [NASA ADS] [CrossRef] [Google Scholar]
  110. Lopes, P. A. A., Trevisan, M., Laganá, T. F., et al. 2018, MNRAS, 478, 5473 [NASA ADS] [CrossRef] [Google Scholar]
  111. Lotz, M., Remus, R.-S., Dolag, K., Biviano, A., & Burkert, A. 2019, MNRAS, 488, 5370 [CrossRef] [Google Scholar]
  112. Ludlow, A. D., Navarro, J. F., Boylan-Kolchin, M., et al. 2013, MNRAS, 432, 1103 [NASA ADS] [CrossRef] [Google Scholar]
  113. Lynden-Bell, D. 1967, MNRAS, 136, 101 [NASA ADS] [CrossRef] [Google Scholar]
  114. Macciò, A. V., Dutton, A. A., & van den Bosch, F. C. 2008, MNRAS, 391, 1940 [NASA ADS] [CrossRef] [Google Scholar]
  115. Mahdavi, A., & Geller, M. J. 2004, ApJ, 607, 202 [NASA ADS] [CrossRef] [Google Scholar]
  116. Mahdavi, A., Geller, M. J., Böhringer, H., Kurtz, M. J., & Ramella, M. 1999, ApJ, 518, 69 [NASA ADS] [CrossRef] [Google Scholar]
  117. Mamon, G. A., Biviano, A., & Murante, G. 2010, A&A, 520, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  118. Mamon, G. A., Biviano, A., & Boué, G. 2013, MNRAS, 429, 3079 [NASA ADS] [CrossRef] [Google Scholar]
  119. Mamon, G. A., Cava, A., Biviano, A., et al. 2019, A&A, 631, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  120. Margoniner, V. E., Lubin, L. M., Wittman, D. M., & Squires, G. K. 2005, AJ, 129, 20 [NASA ADS] [CrossRef] [Google Scholar]
  121. Martizzi, D., Teyssier, R., Moore, B., & Wentz, T. 2012, MNRAS, 422, 3081 [NASA ADS] [CrossRef] [Google Scholar]
  122. Martizzi, D., Teyssier, R., & Moore, B. 2013, MNRAS, 432, 1947 [NASA ADS] [CrossRef] [Google Scholar]
  123. Maughan, B. J., Giles, P. A., Rines, K. J., et al. 2016, MNRAS, 461, 4182 [NASA ADS] [CrossRef] [Google Scholar]
  124. McLachlan, G. J., & Basford, K. E. 1988, Mixture Models: Inference and Applications to Clustering (New York: Marcel Dekker) [Google Scholar]
  125. Meneghetti, M., Fedeli, C., Zitrin, A., et al. 2011, A&A, 530, A17 [CrossRef] [EDP Sciences] [Google Scholar]
  126. Meneghetti, M., Rasia, E., Vega, J., et al. 2014, ApJ, 797, 34 [Google Scholar]
  127. Merritt, D. 1985, MNRAS, 214, 25P [NASA ADS] [CrossRef] [Google Scholar]
  128. Merritt, D. 1987, ApJ, 313, 121 [NASA ADS] [CrossRef] [Google Scholar]
  129. Meyer, S. L. 1975, Data Analysis for Scientists and Engineers (New York: John Wiley& Sons Inc.) [Google Scholar]
  130. Moore, B., Governato, F., Quinn, T., Stadel, J., & Lake, G. 1998, ApJ, 499, L5 [NASA ADS] [CrossRef] [Google Scholar]
  131. More, S., Diemer, B., & Kravtsov, A. V. 2015, ApJ, 810, 36 [Google Scholar]
  132. Munari, E., Biviano, A., Borgani, S., Murante, G., & Fabjan, D. 2013, MNRAS, 430, 2638 [Google Scholar]
  133. Munari, E., Biviano, A., & Mamon, G. A. 2014, A&A, 566, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  134. Muzzin, A., Wilson, G., Yee, H. K. C., et al. 2009, ApJ, 698, 1934 [NASA ADS] [CrossRef] [Google Scholar]
  135. Muzzin, A., Wilson, G., Yee, H. K. C., et al. 2012, ApJ, 746, 188 [NASA ADS] [CrossRef] [Google Scholar]
  136. Nantais, J. B., van der Burg, R. F. J., Lidman, C., et al. 2016, A&A, 592, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  137. Natarajan, P., & Kneib, J.-P. 1996, MNRAS, 283, 1031 [NASA ADS] [CrossRef] [Google Scholar]
  138. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
  139. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
  140. Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039 [NASA ADS] [CrossRef] [Google Scholar]
  141. Neto, A. F., Gao, L., Bett, P., et al. 2007, MNRAS, 381, 1450 [NASA ADS] [CrossRef] [Google Scholar]
  142. Newman, A. B., Treu, T., Ellis, R. S., & Sand, D. J. 2013, ApJ, 765, 25 [NASA ADS] [CrossRef] [Google Scholar]
  143. Old, L., Gray, M. E., & Pearce, F. R. 2013, MNRAS, 434, 2606 [NASA ADS] [CrossRef] [Google Scholar]
  144. Old, L. J., Balogh, M. L., van der Burg, R. F. J., et al. 2020, MNRAS, 493, 5987 [CrossRef] [Google Scholar]
  145. Old, L. J., Balogh, M. L., van der Burg, R. F. J., et al. 2021, MNRAS, 500, 355 [CrossRef] [Google Scholar]
  146. Osipkov, L. P. 1979, Sov. Astron. Lett., 5, 42 [NASA ADS] [Google Scholar]
  147. Overzier, R. A. 2016, A&ARv, 24, 14 [NASA ADS] [CrossRef] [Google Scholar]
  148. Peirani, S., Dubois, Y., Volonteri, M., et al. 2017, MNRAS, 472, 2153 [CrossRef] [Google Scholar]
  149. Peter, A. H. G., Moody, C. E., & Kamionkowski, M. 2010, Phys. Rev. D, 81, 103501 [NASA ADS] [CrossRef] [Google Scholar]
  150. Pizzuti, L., Sartoris, B., Amendola, L., et al. 2017, JCAP, 7, 023 [CrossRef] [Google Scholar]
  151. Powell, M. J. D. 2006, in Large-Scale Nonlinear Optimization, eds. G. Di Pillo, & M. Roma (USA: Springer), 255 [Google Scholar]
  152. Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25 [NASA ADS] [CrossRef] [Google Scholar]
  153. Ragagnin, A., Saro, A., Singh, P., & Dolag, K. 2021, MNRAS, 500, 5056 [CrossRef] [Google Scholar]
  154. Rasia, E., Tormen, G., & Moscardini, L. 2004, MNRAS, 351, 237 [NASA ADS] [CrossRef] [Google Scholar]
  155. Rasia, E., Borgani, S., Ettori, S., Mazzotta, P., & Meneghetti, M. 2013, ApJ, 776, 39 [NASA ADS] [CrossRef] [Google Scholar]
  156. Reed, D. S., Koushiappas, S. M., & Gao, L. 2011, MNRAS, 415, 3177 [NASA ADS] [CrossRef] [Google Scholar]
  157. Ricotti, M., Pontzen, A., & Viel, M. 2007, ApJ, 663, L53 [NASA ADS] [CrossRef] [Google Scholar]
  158. Rines, K., & Diaferio, A. 2006, AJ, 132, 1275 [NASA ADS] [CrossRef] [Google Scholar]
  159. Rines, K., Geller, M. J., Diaferio, A., Mohr, J. J., & Wegner, G. A. 2000, AJ, 120, 2338 [NASA ADS] [CrossRef] [Google Scholar]
  160. Rines, K., Geller, M. J., Kurtz, M. J., & Diaferio, A. 2003, AJ, 126, 2152 [NASA ADS] [CrossRef] [Google Scholar]
  161. Rines, K., Geller, M. J., Diaferio, A., & Kurtz, M. J. 2013, ApJ, 767, 15 [NASA ADS] [CrossRef] [Google Scholar]
  162. Rocha, M., Peter, A. H. G., Bullock, J. S., et al. 2013, MNRAS, 430, 81 [NASA ADS] [CrossRef] [Google Scholar]
  163. Sand, D. J., Treu, T., Smith, G. P., & Ellis, R. S. 2004, ApJ, 604, 88 [NASA ADS] [CrossRef] [Google Scholar]
  164. Sarazin, C. L. 1980, ApJ, 236, 75 [NASA ADS] [CrossRef] [Google Scholar]
  165. Sarazin, C. L. 1986, Rev. Mod. Phys., 58, 1 [NASA ADS] [CrossRef] [Google Scholar]
  166. Saro, A., Mohr, J. J., Bazin, G., & Dolag, K. 2013, ApJ, 772, 47 [CrossRef] [Google Scholar]
  167. Sartoris, B., Biviano, A., Rosati, P., et al. 2020, A&A, 637, A34 [CrossRef] [EDP Sciences] [Google Scholar]
  168. Schaller, M., Frenk, C. S., Bower, R. G., et al. 2015, MNRAS, 452, 343 [NASA ADS] [CrossRef] [Google Scholar]
  169. Schwarz, G. 1978, Ann. Stat., 6, 461 [Google Scholar]
  170. Sereno, M., & Covone, G. 2013, MNRAS, 434, 878 [NASA ADS] [CrossRef] [Google Scholar]
  171. Sereno, M., Giocoli, C., Ettori, S., & Moscardini, L. 2015, MNRAS, 449, 2024 [NASA ADS] [CrossRef] [Google Scholar]
  172. Shirasaki, M., Lau, E. T., & Nagai, D. 2018, MNRAS, 477, 2804 [CrossRef] [Google Scholar]
  173. Sifón, C., Battaglia, N., Hasselfield, M., et al. 2016, MNRAS, 461, 248 [NASA ADS] [CrossRef] [Google Scholar]
  174. Solanes, J. M., & Salvador-Solé, E. 1990, A&A, 234, 93 [NASA ADS] [Google Scholar]
  175. Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760 [Google Scholar]
  176. Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  177. Stadel, J., Potter, D., Moore, B., et al. 2009, MNRAS, 398, L21 [NASA ADS] [CrossRef] [Google Scholar]
  178. Stalder, B., Ruel, J., Šuhada, R., et al. 2013, ApJ, 763, 93 [NASA ADS] [CrossRef] [Google Scholar]
  179. Stanford, S. A., Brodwin, M., Gonzalez, A. H., et al. 2012, ApJ, 753, 164 [NASA ADS] [CrossRef] [Google Scholar]
  180. Stark, A., Miller, C. J., & Halenka, V. 2019, ApJ, 874, 33 [CrossRef] [Google Scholar]
  181. Sunyaev, R. A., & Zeldovich, Y. B. 1970, Comm. Astrophys. Space Phys., 2, 66 [Google Scholar]
  182. Taylor, J. E., & Navarro, J. F. 2001, ApJ, 563, 483 [NASA ADS] [CrossRef] [Google Scholar]
  183. Teyssier, R., Moore, B., Martizzi, D., Dubois, Y., & Mayer, L. 2011, MNRAS, 414, 195 [CrossRef] [Google Scholar]
  184. Tiret, O., Combes, F., Angus, G. W., Famaey, B., & Zhao, H. S. 2007, A&A, 476, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  185. Vallés-Pérez, D., Planelles, S., & Quilis, V. 2020, MNRAS, 499, 2303 [CrossRef] [Google Scholar]
  186. van der Burg, R. F. J., Muzzin, A., Hoekstra, H., et al. 2014, A&A, 561, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  187. van der Burg, R. F. J., Rudnick, G., Balogh, M. L., et al. 2020, A&A, 638, A112 [CrossRef] [EDP Sciences] [Google Scholar]
  188. van der Marel, R. P., Magorrian, J., Carlberg, R. G., Yee, H. K. C., & Ellingson, E. 2000, AJ, 119, 2038 [NASA ADS] [CrossRef] [Google Scholar]
  189. Verdugo, T., Limousin, M., Motta, V., et al. 2016, A&A, 595, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  190. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261 [Google Scholar]
  191. Webb, K., Balogh, M. L., Leja, J., et al. 2020, MNRAS, 498, 5317 [CrossRef] [Google Scholar]
  192. White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341 [NASA ADS] [CrossRef] [Google Scholar]
  193. Wilson, G., Muzzin, A., Yee, H. K. C., et al. 2009, ApJ, 698, 1943 [NASA ADS] [CrossRef] [Google Scholar]
  194. Wojtak, R., & Łokas, E. L. 2010, MNRAS, 408, 2442 [NASA ADS] [CrossRef] [Google Scholar]
  195. Wojtak, R., Łokas, E. L., Mamon, G. A., et al. 2007, A&A, 466, 437 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  196. Zhao, D. H., Jing, Y. P., Mo, H. J., & Börner, G. 2003, ApJ, 597, L9 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: CLUMPS

In this section we present a new method, called CLUster Membership in Phase Space (CLUMPS), that is calibrated on numerical simulations, and for which particular care is given to poor systems. The main idea of this approach is that the exact position of a galaxy in the phase space plane is not fundamental. Rather, what matters is that many galaxies cluster around the zero-velocity value of this plane.

For this reason the galaxy density field in the phase-space plane of a cluster is computed in a grid, counting the number of galaxies that lie in each bin of the grid. To construct this grid we need to specify the step sizes dR and dV in the radial and velocity directions, respectively. These values are converted into an integer number of bins by rounding up the ratio of the extent to the step in each direction. Specifically, the density field is computed well outside the range of the data: in the radial direction the grid starts at 0.999999 times the minimum radius up to twice the maximum radius, while in the velocity direction from twice the minimum velocity up to twice the maximum, doubling the number of bins in both directions. This density field is then smoothed by convolving it with a two-dimensional Gaussian kernel. To define the kernel we need to specify the widths σR and σV of the Gaussian in the radial and velocity directions, respectively.

Once the density field is smoothed, we consider the smoothed velocity distribution in each radial bin. Since we assume that the velocity center of the cluster is known exactly, the algorithm looks for the maximum closest to zero and the two minima that enclose this peak. The galaxies in each radial bin that have a velocity value within these two minima are considered members. In all radial bins, except for the one closest to the cluster center, instead of requiring that the maximum must be the closest to zero, it is required that the maximum must be the closest to the velocity position of the maximum in the previous radial bin. This is done to avoid jumps in the peak position from bin to bin due to the presence of structures close to the zero-velocity. To find the minima we make use of the SciPy (Virtanen et al. 2020) routine argrelmin, which has a parameter to set how deep a minimum must be to actually be considered as such. This is important for small fluctuations in the density not to be considered minima. We set this value to four.

Here we note that the algorithm performance depends on four free parameters, namely dR, dV, σR, and σV. They are chosen by applying the algorithm with different values of such parameters to a mock catalog, for which the true membership is known. By comparing the latter with the one predicted by the algorithm it is possible to choose the set of parameters that provides the best results. In particular, we retrieved data from the Millennium database (Springel et al. 2005; Lemson & the Virgo Consortium 2006) searching the Henriques et al. (2015) light cones. We first downloaded a catalog of halos satisfying the conditions M200 > 1014 M⊙h−1 and 0.1 < z < 0.5. We then created fields of view (FOVs) that are centered on halos with 1 < RA < 89 deg and −24 < Dec < 24 deg, and considered all the galaxies that satisfy the following condition: 0 < RA < 90 deg and −25 < Dec < 25 deg, 0.1 ≤ zobs ≤ 0.5, zobs being the observed redshift, therefore including peculiar velocities. The center of the FOVs is chosen in a smaller area so that we do not consider FOVs too close to the border, avoiding the risk of a halo that is cut by the area edge. For the same reason, the FOV center in redshift is restricted to 0.11 < z < 0.49. For each FOV we retain all the galaxies that lie within a circle of radius 6 Mpc. To avoid FOVs with too few objects, we retain FOVs that have at least ten galaxies within 1 virial radius, and with |vrf|< 6000 km s−1. This catalog contains 10094 FOVs.

To find the fiducial set of parameters, we run the algorithm on the simulated sample varying the parameters. For each set of parameters and for each FOV we consider all the galaxies with |vrf|< 4000 km s−1 and within 3 Mpc from the cluster center, and compute the completeness C and the purity P, defined as

(A.1)

(A.2)

where Nmem ∩ < Rvir ∩ CLUMPS is the number of members within 1 virial radius that are also identified as members by CLUMPS, Nmem ∩ < Rvir is the number of members within 1 virial radius, and NCLUMPS ∩ < Rvir is the number of galaxies within 1 virial radius that are considered members by CLUMPS. We also measure the velocity dispersion of the candidate members and that of the actual members. Then, for each run, we divide the sample in bins of richness (number of candidate members within 1 virial radius), and for each bin we compute the median completeness, the median purity, and the median value of the velocity dispersion divided by the true velocity dispersion.

We find a good compromise among these quantities when using dR = 0.075 Mpc, dV = 105 km s−1. For σR and σV, instead of using single values, we proceed as follows. For each run the code is actually run 100 times, varying σR from 1 to 6 in ten linearly separated values, and σV from 1 to 5 in ten linearly separated values. The galaxies identified as members in at least half of these runs are considered as final members.

With such a configuration we find that for the systems with (estimated) richness between 30 and 50, the completeness is 0.95, the purity 0.71, and the velocity dispersion ratio 0.97. When considering richer systems (between 50 and 100) the completeness increases to 0.97, the purity decreases to 0.67, and the velocity dispersion ratio is 0.96.

CLUMPS therefore appears to perform reasonably well for systems with a number of members in the range of our 14 clusters, 28 ≤ Nm ≤ 86 (see Table 1). In particular, the fact that the cluster velocity dispersion is recovered with 5% accuracy is particularly important for the dynamical analysis that relies on the velocity distribution of cluster galaxies.

All Tables

Table 1.

List of cluster properties.

In the text
Table 2.

Properties of the ensemble cluster and its subsamples.

In the text
Table 3.

Concentrations of the ensemble cluster and its subsamples.

In the text

All Figures

thumbnail Fig. 1.

Normalized projected phase-space diagram of the ensemble cluster obtained using r200 = r200, Mt. Large dots have membership weights of 1.0, small dots membership weights of 0.5. Escape-velocity curves are shown based on a NFW M(r) with c200 = 3.5 and with a Tiret β(r) with β∞ = 0.5.
In the text
thumbnail Fig. 2.

Points with 1σ error bars: radial spectroscopic completeness profiles for galaxies with M∗ ≥ 109.5 M⊙ in three subsamples of the ensemble cluster: GOGREEN (green dots), GCLASS (blue squares), and GGGC (red crosses). The symbols representing the three subsamples have been slightly displaced along the abscissa for display purposes. Solid, dashed, and dash-dotted curves represent LOWESS smoothing of the data represented by the dots, squares, and crosses, respectively.
In the text
thumbnail Fig. 3.

Projected galaxy number density profile of the ensemble cluster (dots with 1σ Poissonian error bars), corrected for spectroscopic incompleteness, and best fit with the NFW model (curve). Black dots and solid curve: all member galaxies; red squares and dash-dotted curve: quiescent galaxies; blue stars and dashed curve: star-forming galaxies. Best-fit concentration values are listed.
In the text
thumbnail Fig. 4.

MAMPOSSt results for the M(r) of the ensemble cluster. Red curve and orange shading: ℒ-weighted average of all 35 best-fit model M(r) from the MAMPOSSt analysis, and 1σ error. Blue dashed curve: highest-ℒ among NFW models. Black dot-dashed curve: highest-ℒ model (Burkert).
In the text
thumbnail Fig. 5.

Line-of-sight velocity dispersion profile of the ensemble cluster (dots and 1σ error bars, estimated using Eq. (16) in Beers et al. 1990). Also shown are the MAMPOSSt results for M(r) and β(r) projected onto the velocity dispersion profile. Red curve and orange shading: ℒ-weighted average of all 35 best-fit model M(r) from the MAMPOSSt analysis, and 1σ error (reduced χ2 = 0.6). Blue dashed curve: highest-ℒ among NFW models (reduced χ2 = 0.7). Black dot-dashed curve: highest-ℒ model (Burkert, reduced χ2 = 0.4). The last point (empty dot) represents the velocity dispersion of galaxies that were not used in the MAMPOSSt analysis because they are at R ≥ r200.
In the text
thumbnail Fig. 6.

The concentration of the ensemble cluster compared to previous observational results and theoretical predictions. Top panel: green horizontal line: ℒ-weighted average c200 for the ensemble cluster. Shadings represent the 1σ error bar on the c200, with contributions from the three terms, δc, s, δc, ν, and δc, j (inner dark green shading, intermediate light green shading, and external dark green shading, respectively). The two symbols with 1σ error bars labeled “GCLASS” represent the results of the dynamical analysis of GCLASS clusters, as obtained by B+16 (diamond) and as obtained in this work by including the new GOGREEN survey data for the GCLASS clusters. The symbols and 1σ error bars labelled “Lensing” and “X-ray” represent mean values from the literature for clusters in the same r200 and z range as the clusters in our sample, obtained through the analysis of gravitational lensing and X-ray data, respectively (references are given in the text). Bottom panel: green horizontal line and shadings as in the top panel. Symbols with 1σ error bars are theoretical predictions from Bhattacharya et al. (2013), Dutton & Macciò (2014), Correa et al. (2015), Child et al. (2018), De Boni et al. (2013), Ragagnin et al. (2021) for a halo of the same M200 and z as our ensemble cluster, and rescaled to the cosmological parameters used in this paper.
In the text
thumbnail Fig. 7.

Concentrations (c200) for the ensemble cluster and its subsamples (black dots with 1σ error bars). The green, magenta, red, and blue squares with dot-dashed 1σ error bars represent c⋆, cν, cν, Q, and cν, SF, respectively.
In the text
thumbnail Fig. 8.

Velocity anisotropy profile (σr/σθ)(r) for all galaxies of the ensemble cluster (pale green line) and 1σ confidence level (dark green shading). Dot-dashed magenta lines: 1σ intervals of the ℒ-weighted average of the 35 MAMPOSSt best-fit model (σr/σθ)(r). Gray shading: 1σ intervals of the result of B+16 for the GCLASS sample.
In the text
thumbnail Fig. 9.

Velocity anisotropy profiles (σr/σθ)(r) for several subsamples. Confidence intervals are 1σ. Upper left panel: low-M200 (cyan line, blue shading) and high-M200 (red line, orange shading). Lower left panel: low-z (cyan line, blue shading) and high-z (red line, orange shading). Upper right panel: star-forming galaxies (cyan line, blue shading) and quiescent galaxies (red line, orange shading) in the ensemble sample. Lower right panel: 9.5 ≤ log M⋆/M⊙ < 10.5 galaxies (cyan line, blue shading) and log M⋆/M⊙ ≥ 10.5 galaxies (red line, orange shading) in the ensemble sample.
In the text
thumbnail Fig. 10.

The pseudo phases-space density profiles of the ensemble cluster. The solid lines show Q(r) and Qr(r) of the ensemble cluster (top and bottom panels, respectively) derived from the total mass density profile ρ(r) and the number density profile ν(r) (left and right panels, respectively), divided by power laws of the radius Qth ∝ r−1.84 and Qr,th ∝ r−1.92 (green line and shadings) or by power laws Qth ∝ r−2.02 and Qr,th ∝ r−2.11 (red line and orange shading). The shadings indicate 1σ uncertainties. The dashed black line is unity.
In the text
thumbnail Fig. 11.

β − γ relation of Hansen & Moore (2006) with its ±1σ confidence interval (dashed red line and orange shading). The green shading indicates the ±1σ confidence interval for our ensemble cluster.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.