© G. A. Mamon et al. 2019

1. Introduction

The orbits of galaxies in galaxy clusters are a useful tool to understand the evolution of clusters. Galaxies detaching themselves from their initial Hubble expansion should enter clusters on fairly radial orbits. In the inner regions of clusters, most galaxies have arrived at early times, and the two-body relaxation time is often thought to be shorter than the age of the Universe deep inside the cluster where the crossing times are very short. Hence, the galaxies in the inner regions should have forgotten their initial trajectories and the inner population should have isotropic velocities. But clusters grow by mergers, from very minor to major. In the limit of very minor cluster mergers, clusters are relatively isolated systems accreting individual galaxies, whose orbits should be fairly radial on their first infall. In the opposite limit of major cluster mergers, galaxies suffer violent relaxation that isotropizes their orbits. Moreover, the angular momentum of the secondary cluster will be transferred into individual galaxies, which may lead to an excess of more circular orbits.

The measure of the elongations of galaxy orbits is therefore a fundamental tool to understand the formation of clusters, and how galaxy orbits vary with cluster mass, elongation, and large-scale environment. One may go further and understand how orbital elongations depend on galaxy stellar mass or luminosity, as well as galaxy specific star formation rate or color or even morphological type, as in the present study.

Galaxy orbits in clusters are best studied through mass-orbit modeling. Galaxies can be considered tracers of the gravitational potential. The prime method to extract orbital shapes from observational data is through the use of the Jeans equation of local dynamical equilibrium (JE), which states that the divergence of the dynamical pressure tensor is the opposite of the product of the tracer density times the potential gradient. The pressure tensor is the tracer density times the tensor of mean squared velocities, which may be anisotropic (i.e., a velocity ellipsoid with unequal eigenvalues). In spherical symmetry, the stationary JE is

$$\begin{array}{c}\hfill \frac{\mathrm{d}\left(\nu \langle v_r^2\rangle \right)}{\mathrm{d}r}+2\frac{\beta }{r}\nu \langle v_r^2\rangle =-\nu \frac{GM}{r^2},\end{array}$$(1)

where ν(r) and M(r) are the radial profiles of respectively the tracer density and the total mass, $\langle v_r^2\rangle (r)$ is the mean squared radial velocity profile, and

$$\begin{array}{c}\hfill \beta =1-\frac{\langle v_{\theta }^2\rangle +\langle v_{\varphi }^2\rangle }{2\langle v_r^2\rangle },\end{array}$$(2)

represents the (spherical, scalar) anisotropy of the mean squared velocities, called the velocity anisotropy parameter (or simply anisotropy), and is expressed in terms of the root-mean-squared (rms) velocities (β usually varies with radius). Radial, isotropic, and circular orbits respectively have β = 1, β = 0, and β → −∞. For equilibrium systems, there are no net meridional streaming motions, hence $\langle v_{\theta }^2\rangle ={\sigma }_{\theta }^2$. For galaxy clusters, one neglects rotation leading to $\langle v_{\varphi }^2\rangle ={\sigma }_{\varphi }^2$. By symmetry, σ_ϕ = σ_θ. On the other hand, the JE contains radial streaming motions, for example from infall.

In addition to such a Jeans analysis, another class of analysis uses the collisionless Boltzmann (Vlasov) equation (CBE), which states that the six-dimensional (6D) fluid is incompressible in the absence of galaxy collisions. In fact, the JE is a direct consequence (first velocity moment) of the CBE. One can attempt to find a well-behaved and realistic 6D distribution function (DF), expressed in terms of energy and angular momentum, whose moments match the data. In particular, the distribution of tracers in projected phase space (projected radii and relative line-of-sight velocities, PPS) can be expressed as a triple integral of the DF (Dejonghe & Merritt 1992).

There are, however, many hurdles in extracting the orbital shapes from either Jeans or DF analyses. 1) Clusters are observed in projection, with only two positional (sky) coordinates and one (line-of-sight, LOS) velocity. 2) The lack of information on the depth coordinate causes observers to mix clusters along the LOS. 3) Clusters tend to be prolate systems (Skielboe et al. 2012) (although not far from spherical symmetry), as well as prolate in phase space (Wojtak 2013). 4) The spherical JE contains a single equation linking the unknown mass profile (M(r) = (r²/G) dΦ/dr) and anisotropy (linked to the orbital shapes), which is an issue called the mass-anisotropy degeneracy (MAD, Binney & Mamon 1982). The DF analysis suffers from an analogous degeneracy between potential and DF. 5) Streaming motions complicate the analysis. 6) The JE includes partial time derivatives, which are very difficult to estimate (see Falco et al. 2013).

There are many variants of Jeans and DF analyses (see Chap. 5 of Courteau et al. 2014 for a partial review with a focus on galaxy masses). A first class of methods attempts to fit models to the data. The simplest method is to fit models of the LOS velocity dispersion profile to the observed one. While the LOS velocity dispersion profile is a double integral of the tracer density and total mass profiles, it can be simplified to single integral formulations for simple anisotropy profiles (Mamon & Łokas 2005). In a non-rotating system, the velocity dispersion anisotropy affects the shapes of the distribution of LOS velocities (Merritt 1987). The MAD can thus be (partially) lifted by folding in the LOS velocity kurtosis profile (Łokas 2002; Łokas & Mamon 2003; Richardson & Fairbairn 2013; Read & Steger 2017). Another class of methods inverts the data to recover models. For example, assuming a mass profile, which may be deduced from other types of observations (e.g., X-rays, strong and weak gravitational lensing), one can derive a non-parametric anisotropy profile (Binney & Mamon 1982; Tonry 1983; Bicknell et al. 1989; Solanes & Salvador-Solé 1990; Dejonghe & Merritt 1992), a process called anisotropy inversion¹.

A disadvantage of Jeans analysis methods is that they usually require the radial binning of the data². One way around this is to specify the form of the LOS velocity distribution function. Starting with Walker et al. (2009), it has become popular to assume a Gaussian LOS velocity distribution function (in studies of dwarf spheroidal galaxies, but this has not yet been done for clusters). But since the LOS velocity distribution function depends on the velocity anisotropy (Merritt 1987), it is not desirable to measure the anisotropy in this manner. In DF methods, one specifies a form for the DF written as f = f(E, L), where E is energy and L is angular momentum. In particular, models with constant anisotropy β have f ∝ L^−2β. One can then compute not only moments in radial bins, but also at specific positions of the tracers in PPS. Early studies used fairly arbitrary choices for the DF. Wojtak et al. (2008) found that a separable form for f(E, L) provided a good match to the halos in cosmological N-body simulations. This method, adapted to observational data by Wojtak et al. (2009), is powerful, but slow as it involves computing triple integrals (Dejonghe & Merritt 1992). Nevertheless, it has been successfully applied to clusters (Wojtak & Łokas 2010; Wojtak et al. 2011) and galaxies (Wojtak & Mamon 2013). A promising method is to express the DF in terms of action-angle variables (Vasiliev 2019).

In the Modeling Anisotropy and Mass Profiles of Observed Spherical Systems (MAMPOSSt) method (Mamon et al. 2013, hereafter MBB13), the DF is no longer expressed in terms of E and L, but in terms of the three-dimensional (3D) velocity distribution function, the simplest form being a Gaussian distribution. This greatly accelerates the method as it only involves single integrals to predict the observed distribution of tracers in PPS. It is a hybrid model, because while MAMPOSSt does not involve radial binning and assumes a (velocity) distribution function, it uses parametric forms for the total mass and velocity anisotropy profiles and solves the JE for $\langle v_r^2(r)\rangle$ to compute the likelihood of the distribution of tracers in PPS. Using mock clusters from cosmological simulations, MBB13 found that MAMPOSSt lifts the MAD, with slightly comparable accuracy on the mass normalization and scale as the dispersion-kurtosis method of Łokas & Mamon (2003) (according to the tests by Sanchis et al. 2004) and the DF method of Wojtak et al. (2009)³. In both comparisons, MAMPOSSt did much better on the velocity anisotropy, reaching twice the accuracy on $log({\langle v_r^2\rangle }^{1/2}/{\sigma }_{\theta })$.

There have been many attempts to measure the anisotropy of galaxy orbits in clusters. In a pioneering study, Merritt (1987) attempted anisotropy inversion on the Coma cluster with 300 tracers, but was not able to determine the orbital anisotropy of galaxies, given his uncertainty on the mass profile. Łokas & Mamon (2003) considered both LOS dispersion and kurtosis profiles of the Coma cluster and determined a slightly tangential (assumed constant) anisotropy (with large uncertainty).

Another way to lift the MAD is to adopt the mass profile from other methods. Applying anisotropy inversion (from Bicknell et al. 1989) to the mass profile of Abell 1689 (z = 0.18) derived from weak lensing, Natarajan & Kneib (1996) found that the velocities are isotropic in the core, and become radial (β = 0.5) at 800 kpc, which corresponds to roughly r₂₀₀/3 given the published virial masses (Sereno et al. 2013; Lemze et al. 2008)⁴. Benatov et al. (2006) applied the same anisotropy inversion to five clusters from the Cluster And Infall Region Nearby Survey (CAIRNS) (z = 0.03 to 0.3) whose mass profiles were obtained from X-rays and/or lensing. They derived anisotropy profiles that were radial at the very center, and showed a diversity of profiles in their bodies, with two clusters showing slightly tangential or isotropic velocities at r₂₀₀, one mildly radial (β = 0.3) and two fully radial (β ≃ 0.95 at r₂₀₀). Analyzing a mere 64 dwarf galaxies in Coma, Adami et al. (2009) adopted the mass profile derived by Geller et al. (1999) with the caustic method (Diaferio 1999), then fit a constant β to the LOS dispersion profile, to obtain β = 0.4 ± 0.2 or 0.7 ± 0.1 depending on their fit of the number density profile. So for the Coma cluster, the orbital anisotropy may be a function of galaxy mass (since the sample analyzed by Łokas & Mamon 2003 was composed of more luminous galaxies than that analyzed by Adami et al. 2009).

Different galaxy types are often thought to have different anisotropies. For example, early-type galaxies prefer dense regions of clusters, and are expected to have fallen into the cluster at early times and relaxed to isotropic velocities. In contrast, spiral galaxies are thought to be falling in clusters on fairly radial orbits, and perhaps bouncing out on similarly radial orbits. It is thus important to separate the two populations, which are expected to have very different kinematics.

Biviano & Katgert (2004) used anisotropy inversion (following the method of Solanes & Salvador-Solé 1990, hereafter SS90) on a joint analysis of 59 ESO Nearby Abell Cluster Survey (ENACS) clusters with over 20 member velocities per cluster, and were the first to measure the radial variations of the orbits of early versus late spiral morphological types. They assumed that early-type galaxies had isotropic velocities as inferred from the roughly Gaussian LOS velocity distribution that Katgert et al. (2004) previously found for non-central early types, following the predictions of Merritt (1987), which enabled Katgert et al. (2004) to first determine the mass profile. Biviano & Katgert (2004) found that, for early-type spirals (Sa, Sb), β rises to 0.7 at half the cluster virial radius, r₂₀₀, then falls to 0.35 at r₂₀₀. In late-type spirals, they found that β = 0 (isotropic velocities) out to 0.6 r₂₀₀ and then rises to β = 0.3 at r₂₀₀. They did not consider lenticular (i.e., S0) galaxies, nor has anybody else until now. The results of Biviano & Katgert (2004) were confirmed by Munari et al. (2014), who studied Abell 2142 (z = 0.09), first determining the mass profile from a combination of X-ray, lensing and dynamical studies, and using anisotropy inversion to deduce that red galaxies have isotropic orbits at all radii, while the orbits of blue galaxies are isotropic in the inner regions and more radial outside. Biviano & Poggianti (2009) analyzed two stacks of cluster galaxies finding that the orbits of non-emission-line galaxies and emission-line galaxies are similar in the z ∼ 0.56 stack, while non-emission-line galaxies move on more isotropic orbits in the z ∼ 0.07 stack. But the statistical evidence for this evolution is very weak.

Other studies of Biviano and collaborators point to somewhat radial outer orbits for passive galaxies. Analyzing a z = 0.4 cluster with 600 member velocities, Biviano et al. (2013) concluded (using MAMPOSSt and performing SS90 anisotropy inversion from the MAMPOSSt mass profile) that both star forming and passive galaxies have isotropic orbits inside and radial orbits outside. Working on stacked z ∼ 1 clusters from the Gemini Cluster Astrophysics Spectroscopic Survey (GCLASS), Biviano et al. (2016) determined the mass profile with MAMPOSSt and performed SS90 anisotropy inversion and found that both passive and star forming galaxies show radial outer anisotropy (β = 0.4, with large error bars). Annunziatella et al. (2016) performed SS90 anisotropy inversion on a parametric Navarro-Frenk-White (NFW, Navarro et al. 1996) model fit to lensing data for Abell 209 (z = 0.21), and found that passive galaxies display radial outer anisotropy, while their inner anisotropy depends on their stellar mass (slightly radial at high mass and tangential at low mass). Capasso et al. (2019) analyzed non-emission-line galaxies in 110 SZ-selected clusters at 0.26 <  z <  1.32, using both MAMPOSSt and anisotropy inversion. They concluded that passive galaxies have isotropic inner orbits and more radial outer orbits (their outer anisotropy varies with increasing redshift in an oscillatory manner).

Other anisotropy inversions of clusters led to different conclusions. Hwang & Lee (2008) studied ten nearby (z <  0.12) clusters, determining the mass profiles determined from X-ray analyses, and performed anisotropy inversion (using the technique of Bicknell et al. 1989). They deduced that galaxies typically have radial orbits in the core, sometimes (e.g., for Abell 1795) dropping to very tangential β ≈ −3 at r₂₀₀/2. Aguerri et al. (2017) analyzed Abell 85 (z = 0.06) using a parametric NFW model for the mass profiles obtained by the caustic method (Diaferio & Geller 1997) and from X-ray data, and applying their own anisotropy inversion equations (which turn out to be equivalent to those of SS90, as shown in Appendix A). They found isotropic outer orbits (β = 0.0 ± 0.3) for blue dwarf galaxies, but very radial outer orbits for red galaxies (β = 0.7 ± 0.2). This is the first study to point to red (or passive or elliptical) galaxies having more radial outer orbits than blue (or star forming or spiral) galaxies. Could the hierarchy of outer radial anisotropy versus the morphological type or specific star formation class depend on the cluster?

One can alternatively blame the high sensitivity of anisotropy inversion to the required extrapolation of both the data and the model tracer density and mass profiles both outwards to r → ∞ and inwards to r = 0. Moreover, all anisotropy inversion algorithms involve differentiating the observational data – they require the knowledge of $\text{d}\left[\mathrm{\Sigma }(R){\sigma }_{\text{los}}^2(R)\right]/\text{d}R$, where Σ is the surface density⁵. This is where the DF methods have an advantage. Wojtak & Łokas (2010) studied nearby clusters with their state-of-the-art DF method, which they analyzed individually and then jointly assuming common anisotropy profiles. In eight out of ten clusters, they found that the distribution of galaxies in PPS implied isotropic inner orbits and radial outer orbits (β = 0.7 at 1 to 1.5 r₂₀₀). Using the same technique, Wojtak & Mamon (2013) studied the kinematics of satellite galaxies around galaxies themselves (i.e., in small groups), and deduced that satellites around red galaxies lied on orbits with radial outer anisotropy, but they did not separate the satellites according to color or morphological type.

In this article, we study the dependence of the velocity anisotropy profiles of galaxies in clusters on their morphological type, distinguishing between elliptical, spiral, and (for the first time) lenticular galaxies. We thus seek to settle the debate on the different orbital shapes of elliptical and spiral galaxies in clusters, but also wish to better understand S0 galaxies through their orbital shapes. Should the S0 galaxies resemble spiral galaxies in their orbital shapes, because they both have disks, or should they resemble more elliptical galaxies, because they have large bulges and old stellar populations? Answering this may help understand whether S0s originate from spiral galaxies that saw their disks fade, grew their bulges by mergers, or possibly even originate from ellipticals that accreted disks.

We use the WIde field Nearby Galaxy clusters Survey (WINGS), which contains X-ray-selected (median luminosity $L_X^{0.1-2.4\text{keV}}={10}^{43.75}\text{erg}{\text{s}}^{-1}$) clusters at redshifts 0.04 <  z <  0.07 and that are located on the sky at least 20° from the Galactic Plane (Fasano et al. 2006). The WINGS spectroscopic dataset (Cava et al. 2009; Moretti et al. 2014) has been complemented with redshifts from literature data collected through the SDSS-DR7⁶ and NED⁷ databases. Thanks to the stacking procedure, the size of the sample we analyze here is roughly an order of magnitude larger than those used in previous studies of individual clusters (which were typically limited to a few hundred members). In the first article (Cava et al. 2017, hereafter Paper I) of the present series of articles on WINGS clusters, the sample of WINGS clusters was split between regular and irregular clusters. Here, we focus on the regular clusters, and consider the sub-populations of elliptical, S0, and spiral-irregular galaxies as different tracers of the same gravitational potential.

We adopted MAMPOSSt as our primary tool to simultaneously extract the mass profile and the anisotropy profiles of the different classes of galaxies. This will be the first application of MAMPOSSt to a large sample of stacked nearby clusters. Stacking clusters reduces the intrinsic effects of triaxiality that create unavoidable biases in the derived mass profile parameters of individual halos (MBB13).

The outline of this article is the following. We present the MAMPOSSt mass-orbit modeling algorithm in Sect. 2. In Sect. 3, we explain the data sample and the stacking method, while in Sect. 4 we explain the practical implementation of MAMPOSSt, in particular the radial profiles adopted for number density, surface number density, mass, and anisotropy. In Sect. 5, we determine the mass and anisotropy profiles of the stacked samples. We discuss our results in Sect. 6 and provide our conclusions in Sect. 7. We assume a Λ Cold Dark Matter (ΛCDM) cosmological model with Ω_m, 0 = 0.3, Ω_Λ, 0 = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. The MAMPOSSt algorithm

In its standard implementation, MAMPOSSt (MBB13) performs a maximum likelihood fit of the distribution of galaxies in PPS, using parameterized models for the radial profiles of total mass and velocity anisotropy, as well as for the radial number profile and its corresponding surface number density and projected number profiles. MAMPOSSt assumes spherical symmetry, negligible streaming motions, and a form for the 3D velocity distribution of the tracers (taken to be a Gaussian distribution in its current implementation).

In Paper I, we determined the number density profiles of the three morphological types by fits of NFW plus constant background models on the photometric data. Given the known distribution of projected radii R, the MAMPOSSt likelihood of the distribution of galaxies in PPS (projected radii and LOS velocities v_z) is (MBB13)

$$-\ln \mathcal{L}=-{\displaystyle \sum \ln p(v_z|R,r_{\nu },\mathit{\eta }),}$$(3)

$$p(v_z|R,\mathit{\eta })=\frac{\sqrt{2/\pi }}{\Sigma (R)}{\displaystyle {\int }_R^{\infty }\frac{\nu }{{\sigma }_z}}\exp \left[-\frac{v_z^2}{2{\sigma }_z^2}\right]\frac{r}{\sqrt{r^2-R^2}}\text{d}r,$$(4)

$${\sigma }_z^2(r,R)=\left[1-\beta (r)\frac{R^2}{r^2}\right]\langle v_r^2(r)\rangle ,$$(5)

where η is the vector of parameters describing the radial profiles of mass and velocity anisotropy, while ν and Σ are the model number density and corresponding surface density profiles, respectively. The mean squared radial velocity, $\langle v_r^2\rangle$, in the right-hand side of Eq. (5) is previously determined for a given η, by solving the spherical stationary JE of Eq. (1), which can be inverted by solving for dlnK_β/dlnr = 2 β(r), yielding (van der Marel 1994; Mamon & Łokas 2005)

$$\begin{array}{c}\hfill \nu (r)\langle v_r^2(r)\rangle =\frac{1}{K_{\beta }(r)}{\displaystyle {\int }_r^{\infty }{K}_{\beta }(s)\nu (s)\frac{GM(s)}{s^2}\mathrm{d}s,}\end{array}$$(6)

where $K_{\beta }(r)/K_{\beta }(s)=exp[2{\int }_r^s\beta (t)\text{d}t/t]$ depends on the anisotropy model, and its values are given in Appendix A of MBB13 for simple anisotropy models.

LOS velocity uncertainties ϵ_v are accounted for by MAMPOSSt by replacing σ_z in Eq. (4) by $\sqrt{{\sigma }_z^2+{ϵ}_v^2}$. In our dataset, the typical velocity uncertainties are ϵ_v ≈ 53 km s⁻¹, which turn out to be negligible relative to the cluster velocity dispersions. Hence these typical velocity uncertainties have little effect on the MAMPOSSt results. MBB13 tested that splitting the PPS into separate determinations of the number density profile from the distribution of projected radii on one hand, and the mass and anisotropy profiles from the distribution of LOS velocities at given projected radius on the other hand, leads to virtually the same parameters as the standard joint fit of PPS.

Finally, MAMPOSSt can jointly analyze the positions in PPS of several independent tracers, such as elliptical (E), lenticular (S0), and spiral-irregular (S) galaxies. Since the three populations of E, S0, and S galaxies move in the same gravitational potential but have different spatial and velocity distributions (see Fig. 1), making a joint analysis of the three populations by allowing a different β(r) for each of them results in a more stringent constraint on the remaining parameters for M(r). The joint likelihood is the product of those from each of the tracers: ln ℒ = ln ℒ_E + ln ℒ_S0 + ln ℒ_S. In their comparison of mass modeling methods on mock clusters from a semi-analytical model, Old et al. (2015) found MAMPOSSt to perform slightly better on the measure of the virial mass when jointly analyzing red vs. blue tracers instead of grouping them together. Biviano & Poggianti (2009) adopted this approach in their analysis of two sets of clusters.

3. Data preparation

3.1. Sample and interloper removal

In Paper I, we applied the substructure test of Dressler & Shectman (1988) to our initial sample of 73 WINGS clusters containing over 10 000 galaxies. This test led to 15 irregular clusters (at 95 percent confidence), leaving 58 regular clusters. Irregular clusters violate the condition of smooth gravitational potential upon which mass-orbit modeling techniques are based. Foëx et al. (2017) found that discarding clusters identified as irregular with the Dressler & Shectman (1988) statistic leads to a much better match between MAMPOSSt and X-ray based masses. We thus discard the irregular clusters and only consider the 58 regular clusters.

The median apparent magnitude of galaxies in our spectroscopic sample is V = 17.7 (16.8–18.9 quartiles), which translates to an absolute magnitude M_V = −19.2 at the median cluster redshift z = 0.054, that is a luminosity satisfying log(L/L_⊙) = 9.6. Assuming galaxy mass-to-light ratios of M/L_V = 3 (for ellipticals, from Fig. 6 of Auger et al. 2010, who adopted a Chabrier 2003 initial mass function) and 2 (for spirals), our median stellar mass is log(M_stars/M_⊙) = 10.0 for spirals and 10.4 for ellipticals. The galaxy morphologies were determined with the MORPHOT automatic tool (Fasano et al. 2012), and following Paper I, we assigned galaxies to classes E (ellipticals, T_M ≤ −4), S0 (lenticulars, −4 <  T_M ≤ 0), and S (spirals, T_M >  0).

Following Paper I, we assume that the clusters are centered on their Brightest Cluster Galaxy (BCG), defined within 0.5 r₂₀₀. We ran the Clean algorithm (MBB13) to remove interlopers in LOS velocity $v_{\mathrm{LOS}}=c_{\mathrm{light}}(z-\overline{z})/(1+\overline{z})$, where $\overline{z}$ is the median cluster redshift (not that of the BCG) and c_light is the speed of light.

Clean begins by searching for significant gaps in the distribution of LOS velocities using the gapper technique of Wainer & Thissen (1976) with a gapper parameter (not a concentration) C = 4 (Girardi et al. 1996). Clean then iterates the membership defined by the criterion $|v-\mathrm{median}(v)|<2.7{\sigma }_{\mathrm{LOS}}^{\mathrm{NFW}}(R)$, where the factor 2.7 was found to optimally recover the LOS velocity dispersion profile of pure NFW models (Mamon et al. 2010). The term ${\sigma }_{\mathrm{LOS}}^{\mathrm{NFW}}(R)$ requires the knowledge of the scale radius r₋₂ and the mass within, M(r₋₂), or equivalently the virial radius and the concentration c₂₀₀ = r₂₀₀/r₋₂. Clean estimates the virial radius from the aperture velocity dispersion, assuming an NFW model and the Mamon & Łokas (2005) velocity anisotropy profile that goes from isotropic in the inner regions to somewhat radial in the outer regions (Mamon et al. 2010) with a transition radius equal to the NFW scale radius. On first pass, the aperture velocity dispersion is measured by the robust median absolute deviation technique (e.g., Beers et al. 1990) and the concentration of the NFW model is taken as c = 4, typical of rich clusters. On subsequent passes, the aperture velocity dispersion is measured using the biweight estimator (e.g., Beers et al. 1990) and the concentration is taken from the concentration-mass relation that Macciò et al. (2008) fit to the halos of dissipationless cosmological N-body simulations. We then restricted our cluster sample to the 54 clusters with at least N_m = 30 members within R₂₀₀ (median velocity dispersion 763 km s⁻¹).

3.2. Stacking

We stacked the clusters into a pseudo-cluster by rescaling the projected radii R_i, j of galaxies in cluster j to the mass-weighted average cluster with $R_{i,j}^{\prime }=R_{i,j}\langle r_{200}\rangle /r_{200,j}$ and the LOS velocities v_i, j to $v_{i,j}^{\prime }=v_{i,j}\langle v_{200}\rangle /v_{200,j}$ (see footnote 4). For each cluster, we estimated r₂₀₀ in three different manners:

1. A velocity dispersion based estimator, sigv, obtained from the Clean algorithm (MBB13).

2. A richness based estimator called Num (Mamon et al., in prep., see Old et al. 2014), which performs a linear fit between log richness and log virial radius. Num performed the best among over 20 algorithms in recovering the value of M₂₀₀, hence that of r₂₀₀ (Old et al. 2014, 2015).

3. An X-ray temperature-based estimator derived from the mass-temperature relation of Arnaud et al. (2005), which we call tempX. But this is limited to the 40 regular clusters with observed X-ray temperatures.

We then applied the Clean procedure one last time to each of these three stacks to remove remaining undetected interlopers. We finally discarded galaxies with projected radii beyond r₁₀₀ ≃ 1.35 r₂₀₀ of our stacked clusters, where we adopted the Clean values of r₂₀₀ found in Paper I. Indeed, r₁₀₀ is the theoretical virial radius, which is thought to be the maximum physical radius where dynamical equilibrium is achieved (i.e., no net radial velocities). Moreover, the Jeans equation is not valid beyond r = 2 r₁₀₀ (Falco et al. 2013), so the limiting projected radius must satisfy R <  2 r₁₀₀, hence our conservative choice of r₁₀₀.

This leaves us with up to nearly 5000 galaxies for our three stacks. We also excluded the very inner region since it is dominated by the internal dynamics of the BCG, rather than by the overall cluster (see, e.g., Biviano & Salucci 2006). Using a minimum projected radius of 50 kpc (roughly 0.03 r₂₀₀) leaves us now with a total of up to 4682 galaxies (for the sigv stack), as displayed in Table 1.

While, according to the mass challenge of Old et al. (2015), Num recovers M₂₀₀ (hence r₂₀₀) with much less scatter than Clean (sigv), it has the drawback that the recovered log mass varies as roughly 0.5 log M_true, thus leading to positive (negative) bias for clusters of low (respectively, high) mass. It thus seems preferable to avoid this non-unity recovered vs. true log mass slope of Num, which may bias our results. We thus adopt sigv as our main stacking method, but will compare in Sect. 6 the sigv results with those from the Num and tempX stacking methods.

3.3. Brief overview of stacked data

Figure 1 displays the projected phase-space diagram, highlighting the distributions of galaxies of different morphological types. Figure 2 shows the LOS velocity dispersion profiles for the three morphological types of galaxies. A major part of the differences in these velocity dispersion profiles arises from the different scale radii of the three number density profiles. The solid curves in Fig. 2 show the predicted LOS velocity dispersion profiles for an NFW model, with virial radius and tracer radii adopted from Table 2 of Paper I, also assuming isotropic orbits and a mass concentration of c₂₀₀ = 4. The curves match well the observed velocity dispersions for elliptical and S0 galaxies, suggesting that these early types do not depart much from velocity isotropy. The sharp reader may notice that at projected radii R <  330 kpc, the three S0 galaxy LOS velocity dispersions are all above the isotropic prediction, while the three elliptical LOS velocity dispersions are all below the isotropic prediction. This suggests that the inner anisotropies of S0s may be somewhat radial while those of ellipticals may be somewhat tangential.

On the other hand, the LOS velocity dispersions of spiral galaxies are significantly higher than expected from the isotropic model. This may be the signature of anisotropic velocities of the spiral population. Alternatively, this may signify that the spiral population is not in dynamical equilibrium. We discuss this possibility in Sect. 6.2.6.

4. Practical implementation of MAMPOSSt

We now present the practical implementation of the version of MAMPOSSt that we used here⁸. We present in turn our adopted radial profiles for number density, mass, and velocity anisotropy.

4.1. Tracer number density profiles

MAMPOSSt allows for the density profiles of the observed tracers to have up to two parameters (scale and shape; no normalization is required for tracers with negligible mass). In this work, we assume one-parameter NFW number density profiles, following the fits to the photometric data performed in Paper I. We express the NFW number density profile as

$$\nu (r)=\frac{N(r_{\nu })}{4\pi r_{\nu }^3}\stackrel{\sim }{\nu }\left(\frac{r}{r_{\nu }}\right)$$(7)

$$\stackrel{\sim }{\nu }(x)=\frac{x^{-1}{(x+1)}^{-2}}{\ln 2-1/2},$$(8)

where r_ν is the tracer scale radius, defined as the radius of logarithmic slope −2, which matches the usual scale radius for the NFW model (but not for the other models presented below). The NFW surface density profile is expressed as (see Bartelmann 1996 and Łokas & Mamon 2001 for similar formulations)

$$\begin{array}{cc}\hfill \mathrm{\Sigma }(R)& =\frac{N_{\mathrm{p}}(r_{\nu })}{\pi r_{\nu }^2}\stackrel{\sim }{\mathrm{\Sigma }}(X)\hfill \\ \hfill & =\frac{1/(x^2-1)-C^{-1}(1/X)/{|X^2-1|}^{3/2}}{2ln2-1},\hfill \end{array}$$(9)

where

$$\begin{array}{c}\hfill C^{-1}(x)=\{\begin{array}{cc}{cos}^{-1}x\hfill & (x<1),\hfill \\ {cosh}^{-1}x\hfill & (x>1),\hfill \end{array}\end{array}$$(10)

and $\stackrel{\sim }{\mathrm{\Sigma }}(1)=(1/3)/(2ln2-1)$.

We adopted the values of the decimal logarithm of the NFW scale radii of the stacked clusters, obtained from fits of a (projected) NFW cluster model plus constant background to the photometric data of the stacked clusters. These values (and their uncertainties) are provided in Table 1 for all three stacks, for galaxy morphological types E (ellipticals), S0 (lenticulars), and S (spirals). We adapted MAMPOSSt to take priors on these log scale radii assuming a Gaussian distribution centered on the value and with a dispersion equal to the uncertainty on log scale radius, and truncated at 3σ.

4.2. Mass profile

In MAMPOSSt, the total mass profile is generally specified by a dark component, possibly massive tracer components, and a possible central black hole. In the present work, the tracer components are generally massless and the black hole mass is assumed to be negligible (fixed to zero), hence the “dark” component generally refers to the total mass, and we call it the “cluster” component, which refers to the total mass unless we include a massive tracer for the BCG. The cluster mass model can have up to three parameters (normalization, scale, and possible shape).

We express the mass profile in terms of the mass at the scale radius times a dimensionless function of radius over scale radius:

$$\begin{array}{c}\hfill M(r)=M(r_{\rho })\stackrel{\sim }{M}\left(\frac{r}{r_{\rho }}\right),\end{array}$$(11)

where r_ρ is the mass density radius, defined as that where the logarithmic slope of the mass density profile is −2. In virial units, one would then have $M(r)/M_{\mathrm{vir}}=\stackrel{\sim }{M}(cr/r_{\mathrm{vir}})/\stackrel{\sim }{M}(c)$, where c  =  r_vir/r_ρ is the concentration parameter. We consider the following dimensionless mass profiles (noting that $\stackrel{\sim }{M}(1)=1$ by definition):

NFW. The usual NFW model, whose density profile is given in Eqs. (7) and (8), for which the inner and outer logarithmic slopes are respectively −1 and −3:

$$\begin{array}{c}\hfill \stackrel{\sim }{M}(x)=\frac{ln(x+1)-x/(x+1)}{ln2-1/2}·\end{array}$$(12)

cNFW. The cored NFW model (generalized NFW with zero inner slope),

$$\begin{array}{c}\hfill \stackrel{\sim }{M}(x)=\frac{ln(2x+1)(3x+1)/{(2x+1)}^2}{ln3-8/9}·\end{array}$$(13)

In the cored-NFW model, r_ρ is equal to twice the scale radius.

gNFW. The generalized NFW model, whose density profile is ρ(r)∝r^γ (r + r_s)^{−3 − γ}, for which the inner and outer logarithmic slopes are respectively γ (e.g., −1 for NFW) and −3 again:

$$\begin{array}{cc}\hfill \stackrel{\sim }{M}(x)& =\frac{I_{-(\gamma +2)x}[\gamma +3,-(\gamma +2)]}{I_{-(\gamma +2)}[\gamma +3,-(\gamma +2)]},\hfill \\ \hfill & =x^{\gamma +3}\frac{{}^2{F}_1[\gamma +3,\gamma +3,\gamma +4,-(\gamma +2)x]}{{}^2{F}_1[\gamma +3,\gamma +3,\gamma +4,-(\gamma +2)]},\hfill \end{array}$$(14)

where I_x(a, b) is the regularized incomplete beta function, while ²F₁(a, b, c, x) is the ordinary (Gaussian) hypergeometric function. In the generalized NFW model, r_ρ = (γ + 2) r_s.

Hernquist. The Hernquist (1990) model with density profile ρ(r)∝r⁻¹ (r + r_s)⁻³, for which the inner and outer logarithmic slopes are respectively −1 and −4:

$$\begin{array}{c}\hfill \stackrel{\sim }{M}(x)=9{\left(\frac{x}{x+2}\right)}^2.\end{array}$$(15)

In the Hernquist model, r_ρ = r_s/2.

Einasto. The Einasto (1965) model was introduced for stellar distributions in the Milky Way, but was found by Navarro et al. (2004) and confirmed by many to fit the density profiles of halos even better than the NFW model:

$$\begin{array}{c}\hfill \stackrel{\sim }{M}(x)=\frac{\gamma (3n,2n{x}^{1/n})}{\gamma (3n,2n)}.\end{array}$$(16)

4.3. Velocity anisotropy profile

MAMPOSSt also allows for a wide variety of velocity anisotropy profiles for each of the tracer components, based on up to three parameters (inner anisotropy, outer anisotropy, and transition radius). We consider the following models for the anisotropy profile, β(r):

T. Tiret et al. (2007) profile,

$$\begin{array}{c}\hfill {\beta }_{\mathrm{T}}(r)={\beta }_0+({\beta }_{\infty }-{\beta }_0)\frac{r}{r+r_{\beta }},\end{array}$$(17)

T₀. The same as the Tiret profile, but with β₀ = 0.

gOM. Generalized Osipkov-Merritt model (Osipkov 1979; Merritt 1985),

$$\begin{array}{c}\hfill {\beta }_{\mathrm{gOM}}(r)={\beta }_0+({\beta }_{\infty }-{\beta }_0)\frac{r^2}{r^2+r_{\beta }^2},\end{array}$$(18)

where the usual Osipkov-Merritt anisotropy is recovered for β₀ = 0 and β_∞ = 1. Note that the usual constant anisotropy model can be retrieved as a singular case of the above T and gOM models, when assuming β₀ = β_∞.

We generally assume that the anisotropy scale radius, r_β, matches the radius of slope −2, r_ν, of the tracer in consideration, which we call the Tied Anisotropy Number Density (TAND) assumption. This is indeed the case for halos of dark matter particles (Mamon et al. 2010).

4.4. Assumptions

Our modeling assumes spherical symmetry for the visible and mass components, neglecting any possible rotation or other non-radial streaming motions. With our choice of regular clusters, we are in a better position to assume that the galaxies are non-interacting tracers of the gravitational potential.

We fit our analytical models to the total mass profile, thus neglecting the contributions of the galaxy and gas components to the cluster. The three galaxy populations are considered as massless independent tracers of the same potential while performing a joint likelihood analysis with MAMPOSSt. This is the first time that such a joint analysis is performed considering galaxy morphological types, and in particular the first study of the orbits of S0 galaxies in clusters.

4.5. Maximum physical radius in integrals

We integrated the inversion of the JE (Eq. (6)) out to 120 Mpc and the LOS integral of Eq. (4) out to 40 Mpc. The former is integrated three times further than the second, to ensure that the radial velocity dispersion is obtained with sufficient accuracy for the LOS integral.

4.6. Free parameters

The following parameters can be set free in MAMPOSSt:

• logarithm of the mass normalization (log M₂₀₀ or equivalently log r₂₀₀);

• logarithm of the mass scale radius (log r_ρ or equivalently of the concentration c₂₀₀ = r₂₀₀/r_ρ);

• inner slope of the mass density profile (from −1.99 to 0);

• logarithms of the tracer scale radii (for each of the three galaxy types;

• inner (r = 0) and outer (r → ∞) symmetrized velocity anisotropies (for each of the three galaxy types)

  $$\begin{array}{c}\hfill {\beta }_{\mathrm{sym}}=2\frac{\langle v_r^2\rangle -{\sigma }_{\theta }^2}{\langle v_r^2\rangle +{\sigma }_{\theta }^2}=\frac{\beta }{1-\beta /2},\end{array}$$(19)

  which can be as low as −2 for circular orbits and as high as +2 for radial orbits, and where β_sym → β for |β|≪1 (we allow −1.8 to 1.8);

• anisotropy transition radius r_β (see Eqs. (17) and (18)) for each of the three galaxy types (unless we assumed TAND).

This amounts to a maximum of 15 free parameters. We also allow ourselves an extra mass component, treated as a massive tracer, potentially adding two extra free parameters (but we then forced an NFW cluster mass model, thus subtracting the free inner slope, for a net single extra parameter).

Given our lack of knowledge on the parameters, we adopt flat priors for all parameters, except for the log scale radii of the E, S0, and S tracer density profiles, determined (externally) from the photometric data (see Paper I), for which we adopt Gaussian priors (using our mean values of log r_ν and their uncertainties).

In our basic set of 30 MAMPOSSt runs, we assumed that the cluster mass concentration follows the ΛCDM relation found by Dutton & Macciò (2014) for massive halos in dissipationless cosmological N-body simulations in a Planck cosmology:

$$\begin{array}{c}\hfill log{c}_{200}=2.13-0.10log\left(\frac{M_{200}}{M_{\odot }}\right).\end{array}$$(20)

We assumed a Gaussian prior on this concentration-mass (c − M) relation with 0.1 dex uncertainty. We also performed additional MAMPOSSt runs with free cluster mass concentration, with different minimum and maximum allowed projected radii, different minimum number of galaxies per cluster in building the stacks, as well as for the different stacks.

4.7. Marginal distributions from the MCMC analysis

We determined the marginal distributions of the k free parameters using the Markov chain Monte Carlo (MCMC) technique. We used the public COSMOMC code (Lewis & Bridle 2002) in the generic mode. COSMOMC uses the Metropolis-Hastings algorithm to probe the distribution of posteriors (likelihoods times priors) in the k-dimensional parameter space. A chain of points in this space is initialized by selecting a point in k-dimensional parameter space from a k-dimensional Gaussian distribution centered on the k-dimensional hyper-cube, with standard deviations σ_k = [Max(θ_k)−Min(θ_k)]/5. We then advanced each chain, by moving (“walking”) from position θ_old to θ_new using a k-dimensional Gaussian proposal density function, with standard deviations equal σ_k/2, that is one-tenth of the allowed range of parameters. Because this proposal density function is symmetric between two consecutive steps, the Metropolis-Hastings algorithm advances the position θ_old in k-dimensional parameter space to with probability

$$\begin{array}{c}\hfill P=\mathrm{Min}[\frac{p\left({\mathit{\theta }}_{\mathrm{new}}\right)}{p\left({\mathit{\theta }}_{\mathrm{old}}\right)},1],\end{array}$$(21)

where the ps are the posteriors. In other words, the point is kept if the posterior is greater than for the previous point. If the posterior of the new point is lower, a random number is drawn and the new point is kept if the ratio of posteriors (between 0 and 1) is lower than the random number. Otherwise, the new point is discarded and the walker remains on the previous point, whose weight is increased by unity.

We ran six chains in parallel on an eight-core computer, each one run for 10 000 k steps in parallel. While the proposal density function is initially circular in k dimensions, COSMOMC periodically updates it to the (elliptical-shaped) parameter covariances of the previous elements of the chain. We then discarded the first 2000 k steps, where the posterior distribution is dependent on the starting points of the chains (the burn-in phase). We computed the radial profiles of mass density and of the velocity anisotropy of the three galaxy types from the marginal distributions of the free parameters (after discarding the burn-in phase).

5. Results

5.1. Model comparison

5.1.1. Preamble

Most studies employ a single set of priors and show their results. Some consider a few extra choices for their priors. Here, we have choices to make on the inner slope of the cluster mass profile, on an additional mass profile for a possible BCG, and on fixed or free inner and anisotropy profiles for all three components. This led us to consider a large number of prior sets.

Table 2 displays the 30 MAMPOSSt runs on our WINGS clusters, stacked according to their velocity dispersions (sigv), and sorted in decreasing maximum likelihood estimate (MLE), that is increasing − ln ℒ_MLE (see Col. 3). Our values of ℒ_MLE are really posteriors, but are close to likelihoods since all our priors are flat (within a wide range), except for Gaussian priors on the tracer and cluster mass log scale radii (roughly 0.1 dex – see Table 1 – and exactly 0.1 dex, respectively).

We mainly considered mass priors using the NFW or gNFW models. But, we later added two priors with the Einasto mass model, with either free index or fixed to n = 6, as found for ΛCDM halos (Navarro et al. 2004). We did not run more Einasto models, given that the n = 6 Einasto is very similar to the NFW model, while free index Einasto models resemble the gNFW models of a given inner slope. Indeed, we found that the mass profile of the n = 6 Einasto model fits that of the NFW model to 8.5 percent relative precision in the range of 0.135 to 13.5 r₋₂ (0.03 to 3 r₂₀₀), equally spaced in log r, for typical cluster concentrations (the index n = 4.4 provides the best fit – 4.8 percent relative precision – to the NFW mass profile in this radial range). Similarly, the mass profiles of Einasto models with free indices (up to n = 25) fit those of given gNFW models to better than 8.5 percent rms relative precision in the same range of radii (6.1 percent for γ ≥ −1.9).

5.1.2. Bayesian evidence methods

Using different priors can lead to different results, so MAMPOSSt results should be analyzed carefully. The runs leading to the highest likelihood ℒ_MLE naturally tend to have the largest number of free parameters (Table 2). But one can ask whether the addition of extra free parameters improves the likelihood significantly, or whether one is over-fitting the data instead. For this, we use both the Akaike Information Criterion (AIC, Akaike 1973):

$$\begin{array}{c}\hfill \mathrm{AIC}=-2ln{\mathcal{L}}_{\mathrm{MLE}}+2N_{\mathrm{pars}}\end{array}$$(22)

and the Bayes Information Criterion (BIC, Schwarz 1978):

$$\begin{array}{c}\hfill \mathrm{BIC}=-2ln{\mathcal{L}}_{\mathrm{MLE}}+ln{N}_{\mathrm{data}}{N}_{\mathrm{pars}}.\end{array}$$(23)

Given our data sample with N_data = 4682 (for the sigv stack), each extra parameter must decrease − ln ℒ_MLE by 4.23 to lead to a better (i.e., lower) BIC value, while with AIC it must decrease by only 1.

According to Kass & Rafferty (1995), when a model has a BIC lower than another model by 2, 6, or 10 units, one can conclude that there is respectively positive, strong, or very strong evidence in favor of the former one. Since BIC penalizes extra parameters much more than AIC, BIC seems preferable to AIC when our model is built with a small number of parameters. Equations (22) and (23) lead to ΔAIC = ΔBIC − (lnN_data − 2) ΔN_pars. Therefore, the condition for strong AIC evidence for the simpler model compared to another one with ΔN_pars extra free parameters, and given our data sample sizes (see Table 1), is ΔAIC >  6 − 6.45 ΔN_pars, that is ΔAIC >  −0.45 compared to a more complex model with a single extra parameter. Conversely, a more complex model is strongly favored if ΔBIC <  −6, leading to ΔAIC <  −6 − 6.45 ΔN_pars, that is ΔAIC <  −12.45 compared to a simpler model with one less free parameter.

5.1.3. Preferred models

The models with the highest likelihoods from the MCMC analysis (lowest −ℒ_MLE) are complex models with more free parameters. Among our 30 models, the models ranked first and fourth in likelihood (models 1 and 4) all have a cluster mass profile that is steeper (i.e., gNFW with γ <  −1) than NFW at inner radii. Note that the free Einasto model with free velocity anisotropy with TAND (not shown in Table 2) has an even lower − ln ℒ_MLE than model 1 (though only slightly). In comparison, our models 2 and 3, which have an additional mass component for a central BCG, also fit well as they rank second and third.

The seven-parameter model 16 (NFW cluster with isotropic ellipticals and T₀ anisotropy with TAND for S0s and spirals) has the lowest AIC. The best BIC is reached for the six-parameter model 23, with n = 6 Einasto mass density, isotropic orbits for ellipticals and lenticulars, and T₀ anisotropy with TAND for spirals⁹. This model is also the second best model for AIC.

5.1.4. Rejection of mass models

There is strong BIC evidence that Einasto model 23 is preferable to all other mass models, except the equivalent NFW model 24, model 25 (with its cored NFW mass), and (marginally) model 16 (which is identical to model 24, but with free outer anisotropy for S0 galaxies).

In particular, there is strong evidence (i.e., ΔBIC >  7) favoring the n = 6 Einasto or NFW cluster mass density profiles of models 23 and 24 compared to the Hernquist (1990) mass profile (model 26, with identical velocity anisotropy). This conclusion is unchanged if we relax the concentration-mass relation of Eq. (20).

There is also strong BIC evidence against replacing the n = 6 Einasto and NFW cluster mass models with a free Einasto index model or a gNFW (free inner slope) model, as both models (17 and 21) have higher BICs by over 7 for the same velocity anisotropy priors. The evidence against a gNFW cluster mass model is even stronger (ΔBIC >  10) in moderately complex anisotropy models 11 vs. 13, but the evidence against gNFW is only weak for more complex anisotropy models, for example 1 vs. 9 and 4 vs. 10. And while the best AIC is reached for NFW cluster model 16, the second best model for AIC is gNFW model 4, which has ΔAIC nearly −2 relative to the analogous NFW model 10. Hence, the case against gNFW is less clear with AIC than with BIC.

Similarly, there is very strong BIC evidence against the need for a specific BCG component in comparison to a single NFW cluster mass model, with ΔBIC = 13 between models 3 and 9 as well as between models 2 and 12. The posteriors and BIC values are even worse (i.e., higher) when assuming that the BCG mass follows the stellar mass, with an n = 4 Prugniel & Simien (1997) model (an excellent approximation to the deprojection of the de Vaucouleurs 1948 surface density model): as seen by comparing models 8 and 3. There is also strong evidence in favor of the gNFW cluster mass model compared to the NFW cluster plus NFW BCG (which has an extra free parameter), with ΔBIC ≃ 9 between models 3 and 1 and ΔBIC >  7 between models 2 and 7. These conclusions are unaltered when using AIC in place of BIC.

5.1.5. Rejection of velocity anisotropy models

Moreover, according to its returned BIC values, MAMPOSSt shows no need for complex velocity anisotropy models. Indeed, there is strong BIC evidence that model 24 with NFW mass profile and T anisotropic outer orbits for the spirals (and isotropic orbits for E and S0 galaxies) is preferable to a) isotropic orbits for the spirals as well as the E and S0 galaxies (model 28); b) anisotropic outer orbits only for ellipticals (model 30) or lenticulars (model 29); c) anisotropic outer orbits for all three morphological types (model 13); d) anisotropic inner orbits only for spirals (model 22); and e) the generalized gOM anisotropy model with inner isotropic orbits for the spirals (model 27). There is also strong evidence against the need for allowing a free anisotropy radius for the spirals (model 20) instead of TAND (and very strong evidence against free anisotropy radius comparing the more complex anisotropy models 9 and 12, both with NFW cluster mass).

Similarly, compared to the lowest BIC model with gNFW mass profile (model 21, with isotropic orbits for ellipticals and S0s and T₀ anisotropy for spirals), there is strong evidence against all variations on the velocity anisotropy, that is a) allowing for outer anisotropy for E and S0 galaxies (in addition to spirals, model 11), b) allowing for fully anisotropic models for all morphological types (model 1), and c) freeing the transition radius of the spiral anisotropy profile (model 18). For more complex anisotropy models, there is no preference for the T model compared to gOM (e.g., comparing models 1 and 6 for gNFW, or models 3 and 5 for an extra NFW BCG). Finally, there is strong evidence that the anisotropy radius is close to the TAND value (models 1 vs. 7 for gNFW mass, 3 vs. 2 for an extra NFW BCG, and 8 vs. 14 for an n = 4 Sérsic BCG).

AIC is more forgiving than BIC for extra parameters. Its preferred anisotropy is with the relatively simple model 16, with isotropic elliptical orbits and T₀ anisotropy for S0s and spirals (both with TAND). However, the second best non-Einasto model for AIC is model 4, where only the inner orbits of S0s and spirals are fixed to isotropic. The third best model for AIC is model 13, which is intermediate in its complexity, with inner isotropy and free outer anisotropy for all three morphological types. It also fails to distinguish between T and gOM anisotropy (both of which involve the same number of free parameters) and also prefers TAND.

5.1.6. Summary of model comparison

In summary, while AIC slightly prefers the NFW mass model over gNFW and the n = 6 Einasto mass model over the free one, BIC strongly rejects the more complex gNFW and free Einasto mass models (except for complex velocity anisotropy models). Both AIC and BIC point to isotropic elliptical orbits and radial outer spiral orbits, but AIC prefers radial orbits for the lenticulars, while BIC finds moderate evidence against anisotropic outer velocities for S0s. Both AIC and BIC prefer the anisotropy radii to be set by TAND, and fail to distinguish between T and gOM anisotropy.

5.2. Best-fitting parameters

We now focus on just a few of the MAMPOSSt runs. We consider the highest likelihood model (model 1), the model with the strongest AIC evidence (model 16), the NFW model with strongest BIC evidence (model 24), the model with gNFW mass with strongest BIC evidence (model 21), and the 2-component mass model with the strongest AIC and BIC evidences (model 3).

Figure 3 shows the MAMPOSSt MCMC posterior marginal distributions and covariances for model 24. All parameters are well fit within their allowed range (except that the outer anisotropy of spirals can reach the physical limit of pure radial orbits). The lower left panel of Fig. 3 shows that the virial radius is anti-correlated with the outer anisotropy of the spiral population (recall that E and S0 galaxies are assumed here isotropic). The mass-anisotropy degeneracy is more acute on the anisotropy (0.5 on β_sym amounting to 0.11 dex on ${\langle v_r^2\rangle }^{1/2}/{\sigma }_{\theta }$) than on the virial mass (3 × 0.006 = 0.017 dex). This precise measurement of the cluster mass with little influence from the anisotropy (of the spirals) is a consequence of the mass inversion of NFW models being most insensitive to the anisotropy near the virial radius (Mamon & Boué 2010, their Fig. 3) where cluster mass is measured¹⁰.

Figure C.1 is similar to Fig. 3, but for model 11, with gNFW mass and T₀ anisotropy for all morphological types. This figure indicates that the inner density slope is weakly correlated with the mass normalization (left panel of third row) and concentration of the mass profile (second panel of third row), as well as with the outer anisotropies for all three morphological types (third panel of the three bottom rows), which are all three correlated.

Table 3 shows the MLE values and the uncertainties from the marginal distributions derived from the MAMPOSSt MCMC for models 1, 21, and 24. For model 24 with an NFW mass model, the virial radius is very well measured leading to an MLE value of r₂₀₀ = 1690 ± 20 kpc, thus with an uncertainty of only 0.005 dex. This value of r₂₀₀ is consistent with the value 1749 ± 64 kpc given in Paper I using the less accurate (see Old et al. 2015) Clean method (which assumes the NFW mass model), as it should be. The gNFW model 21 leads to r₂₀₀ = 1675 ± 23 kpc, which is still consistent with the Clean value. On the other hand, the gNFW model 1 leads to r₂₀₀ = 1507 ± 59 kpc, which is significantly smaller than the Clean value. For models 21 and 1, the inner slope is consistent with the −1 value for NFW.

When it is a free parameter, the mass concentration of model 24 is c₂₀₀ = 3.8 ± 0.4 Given the mass at r₂₀₀ = 1698 ± 24 kpc of 10^14.8 M_⊙, our concentrations are consistent with the values found for relaxed ΛCDM halos (c₂₀₀ = 4.4 according to Dutton & Macciò 2014, see Eq. (20); or 4.0 according to Child et al. 2018) and for relaxed clusters of galaxies using weak lensing (c₂₀₀ ≃ 4.2 Okabe & Smith 2016, which is the median of the 13 measurements within a mass range of 0.3 dex, see Fig. 12 of Child et al. 2018). We further discuss the c − M relation in Sect. 6.1.3.

Moreover our concentration (set free) for model 24 leads to log(r_ρ/kpc) = 2.65 ± 0.05 in comparison with log(r_E/kpc) = 2.67 ± 0.08, log(r_S0/kpc) = 2.82 ± 0.08, and log(r_S/kpc) = 3.19 ± 0.08. Thus, the elliptical galaxies appear to follow the mass, while the distribution of S0s is very slightly (one-third) but quite significantly more extended. In contrast, the distribution of spirals is nearly four times more extended than that of the mass or of the ellipticals. We return to this issue in Sect. 5.4.

5.3. Goodness of fit

Figure 4 shows the LOS velocity dispersion profiles for the elliptical, lenticular, and spiral galaxies predicted from model 1 (gNFW mass with T anisotropy for the three morphological types, all with TAND). The MAMPOSSt model predictions reproduce the data very well.

5.4. Radial profiles

We now show the radial profiles of mass density, mass over number density, and velocity anisotropy. These profiles were computed in radial bins of width 0.2 dex. Extracting the free parameters from 1001 random chain elements (after the burn-in phase) among the 6 (chains) × (10 000–2000) × (# free parameters), that is typically half a million or more chain elements, we computed the set of three radial profiles at each radial bin.

Figure 5 displays the mass density profiles for models 21 (gNFW) and 3 (NFW + NFW for BCG). In the top panel, a gNFW model was assumed by MAMPOSSt, and the density profile prefers to be steeper than NFW, but not significantly (γ = −1.51 ± 0.42 according to Table 3). Only 85% of all chain elements past burn-in produce γ <  −1. Figure 5 and the constraints on the inner slope from Table 3 for model 21 both suggest that the cNFW model (blue) is ruled out. However, as seen in Table 2, model 25, which is the same as model 24, replacing NFW by cNFW, leads to Min (− ln ℒ) that is only 1.6 higher than for model 24. Considering the cNFW model to be a physical one, it has the same number of free parameters as the NFW model and its BIC is only 3.2 higher than that of model 24. So, one cannot reject the cored NFW model for clusters with our WINGS data.

The bottom panel of Fig. 5 shows the mass density profile for model 3 (NFW cluster + NFW BCG). MAMPOSSt was not able to constrain well the BCG mass profile given the minimum allowed projected radius of 50 kpc. MAMPOSSt prefers a BCG with a tiny scale radius, which is not physical.

One may also wonder which morphological type has a number density profile closest to the mass density profile. Figure 6 displays the ratios of mass density over number density for the three morphological types for models 24 and 21. We normalize the number density profiles by eliminating N(r_ν) between Eq. (7) and the average number of galaxies of a given morphological type per cluster, N_tot, between the minimum and maximum allowed projected radii, R_min and R_max, which we model as

$$\begin{array}{cc}\hfill N_{\mathrm{tot}}& =N_{\mathrm{p}}(R_{\max })-N_{\mathrm{p}}(R_{\min })\hfill \end{array}$$(24)

$$\begin{array}{cc}& =N(r_{\nu })[{\stackrel{\sim }{N}}_{\mathrm{p}}(R_{\max })-{\stackrel{\sim }{N}}_{\mathrm{p}}(R_{\min })],\hfill \end{array}$$(25)

where

$$\begin{array}{c}\hfill {\stackrel{\sim }{N}}_{\mathrm{p}}(X)=\frac{1}{ln2-1/2}\frac{C^{-1}(1/X)}{\sqrt{|X^2-1|}}+ln\left(\frac{X}{2}\right)\end{array}$$(26)

for the NFW model, with N_p(1)=(1 − ln2)/(ln2 − 1/2) and where C⁻¹ is given in Eq. (10).

It should be recalled that the NFW model was assumed for the number density profiles of the three morphological types and that these number density profiles were obtained from fits to the photometric data, and thus do not suffer from any spectroscopic incompleteness. In both panels of Fig. 6, the elliptical galaxies trace the mass almost perfectly; S0 galaxies are nearly as good mass tracers as ellipticals (but somewhat more extended), while spirals are much more extended. In runs with free concentration, the ellipticals trace the mass even better. Since model 21 (see the bottom panel of Fig. 6) is based on a gNFW mass profile, there is less agreement between the NFW number density profile of the ellipticals and the gNFW mass profile, but ellipticals remain the best tracers of mass among the three morphological types. Indeed, the mass over elliptical number density ratio is nearly consistent with being constant (horizontal line in Fig. 6), although there may be a need for extra mass in the BCG.

Figure 7 displays the anisotropy profiles for models 24, 9, 1, and 7. Model 24 (our best non-Einasto model in terms of BIC), which assumes isotropic orbits for the E and S0 galaxies, and also inner isotropy for the spirals, shows that the spiral galaxies clearly have radial orbits in the outer regions of clusters. The other three models, with fully free priors on inner and outer velocity anisotropy, confirm that spiral galaxies have increasingly radial orbits at large distances. Early-type galaxies show moderately radial outer orbits for these three models, but all are consistent with isotropy. The similarity in the anisotropy profiles between models 1 and 7 (which only differ in that the latter has a free anisotropy radius) confirms that this radius is close to the scale radius of the tracer density as in the TAND approximation. The inner and outer anisotropies are displayed in Table 4.

5.5. Outer versus inner velocity anisotropies

Table 4 illustrates the details of the anisotropy for four models. If we free the inner and outer anisotropies (but tie the anisotropy radii to the scale radii, as in model 1), we find that the anisotropy at r₂₀₀ of the E and S0 galaxies are also typically somewhat radial (see also Table 3 and Figs. 7 and 8). However, the uncertainties on outer anisotropy are much larger (almost double in β_sym for E vs. S) for these early-type galaxies compared to spirals (Figs. 7 and 8), which explains why BIC evidence prefers having radial outer orbits for the spiral population only. Hence, there is only marginal evidence that the S0 population has radial outer orbits, while the moderately radial orbits of ellipticals is not statistically significant (Table 4). The outer anisotropies of the elliptical and S0 galaxies are less radial when the anisotropy radii are set free (model 7) compared to the analogous TAND model 1.

The inner anisotropies of the three morphological types are always consistent with isotropy (Table 4), where the uncertainties for the spiral population are much smaller for the TAND assumption. But a close inspection of Table 4 indicates that elliptical galaxies have slightly tangential inner values of anisotropy, as expected from our overview of the LOS velocity dispersion profile (Fig. 2). However, this tangential anisotropy is not statistically significant.

Figure 8 provides a clearer way to view the anisotropy profiles, by plotting the value at r₂₀₀ as a function of the value at 0.03 r₂₀₀. We restricted these plots to models with free inner and outer anisotropies for all morphological types. When the anisotropy radius is forced to the scale radius, as favored by the Bayesian (both AIC and BIC) evidence (top panel of Fig. 8), the 95 percent confidence contours for β(r₂₀₀) for spirals are above zero for all values of the inner anisotropy (at 0.03r₂₀₀), which is almost the case for S0 galaxies, but not the case for ellipticals. Moreover, the 95th percent confidence level is always in the direction of increasingly radial anisotropy for the spirals, which is not the case for the S0s and ellipticals. Also, the inner anisotropy of the ellipticals is somewhat tangential (though not significantly), while those of the S0s and spirals appear to be even more isotropic.

On the other hand, by freeing the anisotropy radii (bottom panel of Fig. 8), the outer anisotropies become independent of the inner values, for all three morphological types. The free vs. fixed anisotropy radii have a stronger effect on the contours of outer vs. inner anisotropy than does the mass model. Nevertheless, only spiral galaxies show clearly radial anisotropy at r₂₀₀ (Table 4). The lack of correlation between inner and outer anisotropies is probably due to the wide range of anisotropy radii allowed by the data. Indeed, while the log anisotropy radii (in units of kpc) are allowed to span between 1 and 4, the uncertainty on the best-fit anisotropy radii for the non-TAND runs are typically as high as 1 dex for all three morphological types. Nevertheless, as for the TAND case, spirals are the sole morphological type for which the orbits systematically become more radial from the inner regions to near the virial radius (last column of Table 4, and as seen by the contours for spirals lying above the oblique line in Fig. 8).

6. Discussion

This work represents the largest analysis of velocity anisotropy in cluster galaxies and the first to distinguish the orbits of ellipticals, spirals, and lenticulars using a Bayesian model to predict the distribution of these three morphological types in PPS. We have constructed a stacked cluster, which helps us avoid departures from spherical symmetry, although it introduces artificial phase mixing.

Our conclusions depend on our choice of priors. We have presented 30 choices of priors (and tried many more). We can restrict our conclusions to the simpler set of priors that lead to the highest Bayesian evidence measures (within ΔBIC = 6 of the lowest BIC), or we can analyze the detailed radial profiles expected from the models that reach the highest likelihoods (which are really posteriors), although their BIC Bayesian evidence is so high that they can be strongly rejected relative to the lowest BIC model.

6.1. Mass density profiles of WINGS clusters

6.1.1. General trends

Our highest BIC Bayesian evidence is reached for models 23 and 24, where the mass profile is n = 6 Einasto or NFW with isotropic velocities for ellipticals and S0s, while for the spirals they are isotropic at the center and fairly radial in the outer regions of clusters (Table 2). There is strong Bayesian evidence against the Hernquist model (26), whose outer slope is steeper (−4) then that of NFW (−3). There is only positive (but not strong) evidence against a cored NFW model (relative to model 24). The case against the gNFW and free index Einasto models is less clear with AIC evidence: the lowest AIC value is reached for an NFW model, but the second lowest among non-Einasto models is for a gNFW model.

On the other hand, our best fitting models prefer a mass profile with a free inner slope of order −1.6 ± 0.4, which is marginally consistent with the −1 slope of the NFW model (Table 3 and Fig. C.1), and rejected by BIC evidence. High likelihoods are also attained by summing an NFW model for the cluster with another smaller NFW model for a central BCG (Table 2). But surprisingly, the BCG would require a very high concentration NFW model that is essentially a −3 power law in the innermost regions of the cluster that we analyzed (Fig. 5). The BCG dominates the cluster within the inner 20 kpc, thus well inside the minimum projected radius for which we are confident of our cluster centers (assumed to be at the BCG location) before we stack them. We thus simply do not have enough tracers to constrain the mass density profile within this radius.

6.1.2. Robustness

We ran several models for the two other cluster stacks (Num and tempX). We found that for model 21 (which uses gNFW for the cluster mass), the constraints on the inner slope range from $\gamma =-1.4_{-0.3}^{+0.5}$ for sigv to $-1.7_{-0.2}^{+0.3}$ for tempX and $-1.8_{-0.1}^{+0.4}$ for Num. The steeper inner mass density slopes for the tempX and Num stacks lead to different Bayesian evidence for or against gNFW. Indeed, comparing model 24 (NFW) to 21, we find that AIC prefers gNFW for Num (ΔAIC = AIC(gNFW) − AIC(NFW) = −5.9) and tempX (ΔAIC = −2.5), whereas it slightly favors NFW for sigv (ΔAIC = 1.2). However, BIC does not favor gNFW despite the much steeper gNFW inner slopes. While there is strong BIC evidence against gNFW with sigv (ΔBIC >  7), there is still positive evidence against gNFW with tempX (ΔBIC = 3.7) and Num (ΔBIC = 0.5).

We also test the robustness of the radial profiles of mass over number density for the three morphological types to the choice of stack. Figure 9 shows that for model 24, the (NFW) mass density profile is almost exactly proportional to the elliptical number density profile for all three stacks, and close to proportional to the S0 number density profile. However, the spirals poorly trace the mass profile as they are more extended, thus ρ/ν is more concentrated.

For model 21 (which is the same as model 24, but with gNFW mass instead of NFW), ellipticals and S0s show a U-shaped mass-over-number profile. A close look indicates that the elliptical number density profile traces the mass density profile slightly better (within r₂₀₀) than do the S0 galaxies, except for the tempX stack where the two types trace the mass with similar accuracy. Again, the spiral galaxies poorly trace the mass profile.

6.1.3. Comparison with other work

In comparison, by combining weak lensing at large radii, strong lensing at intermediate radii, and stellar kinematics at low radii to study seven regular clusters, Newman et al. (2013a) deduced that the total mass density profile is close to a gNFW with inner slope −1.2 ± 0.1, while the dark matter follows a gNFW with a shallow slope of −0.5 ± 0.2 (Newman et al. 2013b). Their total mass profile is consistent with ours (i.e., NFW for lowest BIC model 24 as well as $\gamma =-1.4_{-0.3}^{+0.5}$ for model 21).

Figure 10 shows the constraints on the concentration-mass relation obtained with NFW model 24. Interestingly, with free concentration (i.e., flat prior 0 <  log c₂₀₀ <  1), the magenta contours are well matched to the ΛCDM relation of Dutton & Macciò (2014) given in Eq. (20), especially for stacks tempX and sigv. It is not surprising that folding in this relation as a prior, we recover similar contours that are closer to this relation. The concentrations obtained with MAMPOSSt are well matched to those obtained by weak lensing, except for one (fairly old) weak lensing study that strongly underestimates the concentration. Biviano et al. (2017) performed a MAMPOSSt analysis of 49 WINGS clusters, individually, and found a cluster-to-cluster scatter in concentration that was greater than the uncertainties returned from MAMPOSSt (≃0.3 dex) and from the effect of a range of cluster masses (≃0.3 dex dispersion) combined with the −0.1 slope of the concentration-mass relation (leading to a dispersion of 0.03 dex). Their best fit at the median cluster mass yields c₂₀₀ = 3.34, close to the value of our Num stack (Fig. 10). Our other two stacks lie within the confidence band of the concentration-mass relation of Biviano et al. (2017).

With our assumptions of NFW mass and number density profiles, we find that our lowest NFW BIC model 24 indicates that the elliptical population follows the mass, while the lenticulars are slightly less concentrated and the spirals are much less concentrated (top panel of Fig. 6). These results can be compared to the mass traced by the red (van der Marel et al. 2000 for CNOC clusters), early spectral type (Biviano & Girardi 2003 for a stack of 43 2dFGRS clusters), and non-BCG E and S0 (Katgert et al. 2004 for ENACS clusters) galaxies. The much weaker concentration of the spiral population agrees with the much weaker concentration of the blue galaxies relative to the red ones observed by Collister & Lahav (2005).

Allowing for a gNFW mass model, the scale-radius (radius of slope −2) (and thus the concentration of the mass) remains consistent with the corresponding values of the elliptical population, but not of the lenticulars or spirals (bottom panel of Fig. 6). There is a discrepancy between mass and elliptical number at very small radii (bottom panel of Fig. 6) because of the steeper rise with decreasing radius of the gNFW mass profile compared to the NFW number profile of the ellipticals. Since the BCG contributes a large stellar mass at the cluster center, we expect the E stellar mass density profile to follow the total mass density profile even better than does the E number density profile.

Finally, it is surprising that the distribution of ellipticals, which may be assimilated to the dwarf spheroidals orbiting the Milky Way, follows the dark matter in contrast with the subhalos in the Aquarius dark matter-only simulations (Springel et al. 2008). This discrepancy might be attributed to the missing dissipative gas in Aquarius, and the uncertain link between the subhalo and satellite radial distributions given the uncertain radially-dependent link between the minimum subhalo and galaxy stellar masses. One could also ascribe the discrepancy to the NFW assumption for the E, S0, and S radial distributions, but these are consistent with the data (Cava et al. 2017).

6.2. Velocity anisotropy profiles of WINGS clusters

6.2.1. General trends

For the single component mass models, all anisotropy models that differ from that of T₀ anisotropy for the spirals and isotropy for the E and S0 galaxies are strongly rejected by BIC Bayesian evidence, with the sole exception of model 16 (the best one using AIC evidence) where S0 galaxies and spirals have outer anisotropy (Table 2). But the Bayesian evidence of model 24 (with outer anisotropy only for the spirals) against model 16 is ΔBIC = 5.5, that is “positive” and almost “strong” (ΔBIC >  6). Hence, we are confident that the spiral population has radial orbits at r₂₀₀, which are significantly more radial than at 0.03 r₂₀₀. (e.g., Table 4). On the other hand, while early-type galaxies appear to prefer radial orbits at the virial radius, the trend is not statistically significant. The BIC Bayesian evidence suggests that there is no need for radial outer anisotropy of the E, and marginally so for S0 populations, while AIC evidence prefers to also have radial outer orbits for the S0s.

6.2.2. Robustness

We now test the robustness of our results on anisotropy. Figure 11 shows the effect on the inner and outer velocity anisotropies of moving from free mass concentration (i.e., thin contours) to the ΛCDM concentration (Eq. (20), thick contours). Since the concentrations obtained when they have wide (i.e., free) priors end up in the realm of the ΛCDM concentration-mass relation (Fig. 10), there is virtually no difference between the anisotropies obtained with free or ΛCDM concentrations.

Figure 12 indicates that our results are fully robust to our choice of minimum and maximum projected radii. The top panel of Fig. 12 also highlights the effect of changing the definition of the individual cluster centers before the stacking. Indeed, the BCGs can be displaced from the cluster center (e.g., Skibba et al. 2011), although centering on the BCG or the X-rays leads to cuspier cluster profiles than using the barycenter (Beers & Tonry 1986). The figure shows that centering clusters on X-rays instead of BCGs hardly affects the velocity anisotropy of the spirals, but leads to more tangential (respectively radial) inner anisotropy for the ellipticals (resp. S0s).

Figure 13 compares the outer vs. inner velocity anisotropies when we change the minimum number of member galaxies in clusters that we stack. The orbits of spirals are virtually unaffected by the minimum number of cluster members, whereas the ellipticals and S0s both allow somewhat more radial inner orbits with 81 minimum members per cluster.

Figure 14 compares the outer vs. inner anisotropies from stacks computed using three different methods to estimate the r₂₀₀ radii of the individual clusters (see Paper I for details). The outer orbits of S0s are radial for sigv, quasi-radial for tempX, and isotropic for Num. For ellipticals, the outer orbits are only slightly radial for sigv, but isotropic for the other two stacks.

In the sigv and Num stacks, the ellipticals show signs of tangential inner anisotropy, while they do not in the tempX stack, which is consistent with isotropic velocities for the ellipticals at all radii. The Num stack shows isotropic outer velocities for the S0s, while the tempX and especially our standard sigv stacks indicate radial outer orbits. On the other hand, the radial outer orbits of the spirals are robust to the stacking method (but the strongest radial anisotropy is seen in the sigv stack).

6.2.3. Comparison with previous studies

Figure 15 compares our constraints on inner and outer anisotropy with those from the literature. We first compare to previous studies that did not separate galaxies into different classes. The two anisotropy measures of Abell 2218 by Natarajan & Kneib (1996) are consistent with our anisotropies. However, the inner anisotropies of the 10 clusters measured by Hwang & Lee (2008) (only two are shown in Fig. 15: Abell 85 for comparison with another study and Abell 1689 which has the richest data sample) are much more radial than we (and most others) found. Aguerri et al. (2017) and Benatov et al. (2006) found much lower inner radial anisotropy for Abell 85 and Abell 2199, respectively, than Hwang & Lee (2008). Still, Benatov et al. (2006) and Wojtak & Łokas (2010) both find a very wide range of anisotropies at r₂₀₀, some of which are perfectly radial.

Biviano et al. (2013) found very radial anisotropy at r₂₀₀ for the star forming galaxies in MACS J1206, which is consistent with our analysis given the important uncertainty of their anisotropy. The inner and outer anisotropies of the passive galaxies agree with those of our early-type galaxies. The results of Munari et al. (2014) on red and blue galaxies in Abell 2142 are marginally consistent with our respective results on E+S0 and S galaxies, with their anisotropy of the red galaxies at r₂₀₀ appearing too tangential. There is tension between the very radial (β >  0.55) anisotropy of the passive galaxies at r₂₀₀ found by Annunziatella et al. (2016) for Abell 209 (for both low and high stellar mass) and our much lower values of β for the ellipticals. On the other hand, the inner and outer isotropy of passive galaxies found by Capasso et al. (2019) is consistent with the orbits we find for E and S0 galaxies. The anisotropies found in Abell 85 by Aguerri et al. (2017) for blue galaxies are consistent with ours, except that the outer isotropy that they found for their blue dwarf galaxies is in some tension with our results for spirals. On the other hand, they found significantly more radial orbits for red galaxies in Abell 85 compared to our findings for early type galaxies.

The discrepancies in orbital anisotropies between studies of individual clusters and our stacked analysis may be caused by a possible diversity of clusters, either natural or modulation with mass or redshift (the alternative is that the methods and/or priors were different). It is therefore interesting to compare the anisotropies from the analyses of stacked or clusters or from joint analyses of clusters.

Wojtak & Łokas (2010) also performed a joint analysis of their 31 relaxed clusters that shows a small range of radial anisotropies. They found significantly more radial orbits at the virial radius (for which they used the rather large value of 7 times the mass scale radius) than our anisotropies at r₂₀₀ averaged over the three morphological types (see Table 4 and Fig. 15).

Our results for spirals agree with the analysis of ENACS clusters by Biviano & Katgert (2004) for early-type spirals (Sa and Sb), while they are marginally inconsistent with their results for late-type spirals, for which they found tangential inner anisotropies. Note that for these Sa and Sb galaxies, Biviano & Katgert (2004) found β(r) to increase and then decrease. Our results are consistent with those for z ∼ 1 GCLASS clusters of Biviano et al. (2016), who found that both passive (i.e., E and S0) and star forming (S) galaxies show isotropic orbits inside and radial orbits outside.

The uncertainties on anisotropy are always greater at inner radii than at the virial radius (e.g., Table 4), which suggests that there may be a greater range of anisotropies deep inside clusters rather than near their virial) radius. This appears to contradict the wide range of outer anisotropies found by Benatov et al. (2006), Wojtak & Łokas (2010), as well as from Wojtak & Mamon (2013) for the satellites of galaxies that may or may not be brightest group galaxies. On the other hand, it is consistent with the work of Annunziatella et al. (2016) who find that the inner anisotropy of passive galaxies (i.e., ellipticals and S0s) depends on their stellar mass, while their outer anisotropy does not.

6.2.4. Infall

These constraints on inner and outer anisotropy help us to understand the mechanisms and timescales for the transformation of morphological types for the quenching of star formation. This is done by comparing the orbital anisotropy of star forming vs. passive populations.

The simplest view is that spiral galaxies fall onto clusters on nearly radial orbits, and are fairly rapidly transformed into S0 and E galaxies as they orbit through the cluster. Such morphological transformation may occur through processes of galaxy merging (in the cluster envelope, Mamon 1992; or in infalling groups), by galaxy harassment from numerous minor flybys (Moore et al. 1996), by starving the galaxy of its supply of infalling gas either by tidal stripping (Larson et al. 1980) or by ram pressure stripping (Gunn & Gott 1972). In this picture of fairly rapid morphological transformation, assuming a monolithic evolution of clusters (as in the onion-ring model of Gott 1975), galaxies that enter the cluster later (as the spirals) will lie at larger radii. This picture is confirmed in cosmological N-body simulations (e.g., Fig. 11 of Haines et al. 2015). It is also confirmed by the much higher scale radius (lower number concentration) observed for spirals relative to the S0s and ellipticals (Paper I and Table 3 of the current work). The rapid transformation of spirals into early-type morphologies might imply that infalling spirals may not have time to exchange energy and acquire angular momentum from the other cluster galaxies, hence they should not isotropize. Indeed, we find that the outer anisotropy of spirals is greater than that of ellipticals (Fig. 14).

The signs of some radial outer anisotropy for the lenticular and possibly even the elliptical galaxies may indicate that, at the virial radius, the early-type galaxies are a mixture of an isotropized virialized population with other early-type galaxies that are infalling for the first time (mostly the central galaxies and quenched satellites in galaxy groups). This is consistent with the narrower range of outer anisotropies of spirals relative to that of ellipticals and S0s (see Figs. 8 and 14). As one moves to smaller physical radii, galaxies first entered the cluster at earlier times, and thus have had more time to isotropize (see Sect. 6.2.5 below), which would explain the positive gradients in β(r). However, the kinematical evidence for this natural scenario is thin: only spiral galaxies show statistical evidence of increasingly radial anisotropy profiles as one moves from 0.03 r₂₀₀ to r₂₀₀ (see Table 4).

6.2.5. Isotropization

At small radii, the great majority of early-type galaxies is expected to have entered the cluster long enough ago to have been morphologically transformed from their spiral progenitors (again, as in the onion model of cluster growth of Gott 1975). This raises the question: Should early-type galaxies isotropize or retain the radial orbits of their spiral progenitors?

The natural way for them to isotropize is by two-body relaxation with other galaxies. The typical timescale for two-body relaxation roughly scales as N/(8lnN) times the orbital time (Eq. (4.9) of Binney & Tremaine 1987), which for NFW models is never less than e (i.e., 2.718) times the crossing time (at any radius). However, this formula assumes that the system is self gravitating, whereas galaxies in clusters account for a small portion of the cluster mass. In other words, the dominant dark matter in cluster leads to much greater galaxy velocities than expected from their number and mean mass.

In Appendix D, we perform a more precise quasi-analytical measurement of the two-body relaxation time of galaxies infalling into clusters. For the relatively low median galaxy mass of our sample (10^10.0 M_⊙, Sect. 3.1), the mean number of galaxies per cluster in our sample (87), and for reasonable choices of the pericenter radius and the ratio of apocenter to pericenter, we find (see Fig. D.1) that at the very least 30 orbits are required to isotropize. In considering a growing NFW cluster, Tollet et al. (2017) found that (see Fig. B1 therein) it takes 3 Gyr for a galaxy to move from pericenter to its second apocenter (assuming hereafter that r_apo/r_vir(t_apo) = 3.5), and thus over 4 Gyr, between the last two pericentric passages. These orbital times were shorter at early times, but the number of orbits is expected to be less than a dozen. This suggests that two-body relaxation suffered by a single galaxy is insufficient to explain the isotropy of the early-type galaxies in the inner regions (or further out in some of the stacks).

Galaxies also lose energy by their encounters with other galaxies and especially with dark matter particles (since dark matter dominates the cluster mass distribution). The dynamical friction (df) time is on the order of [M(r)/m]/ln(1 + M(r)/m) times the orbital time (Jiang et al. 2008), and the ratios of cluster mass M(r) to galaxy subhalo mass m are too high for df to be effective (especially since tidal stripping of subhalos by the cluster gravitational field leads to much lower subhalo masses). However, df affects the groups that fall into clusters, hence the galaxies within these groups will lose their radial velocities by df on their host groups. Moreover, the infalling groups will be distorted by the cluster tidal field, and this tidal heating will lead to tidal braking, meaning the transfer of orbital energy into internal energy. However, simulations indicate that galaxies bounce out of clusters to 1 to 2.5 virial radii (Mamon et al. 2004; Gill et al. 2005), suggesting that not all galaxies lose their orbital energy by the df and tidal braking of their host groups.

Since clusters do not evolve monolithically, but grow by mergers, galaxies may see their orbits perturbed by the rapidly varying gravitational potential. This violent relaxation (Lynden-Bell 1967) occurring during major cluster mergers should transfer angular momentum from the other cluster into the galaxies, leading to isotropization. According to cosmological N-body simulations (Fig. 3 of Fakhouri et al. 2010), cluster-mass halos typically undergo 0.8 major mergers since z = 1 (7 Gyr for the cosmology of the simulation studied). Thus, since z = 1, 1 − e^−0.8≃ one-half of clusters undergo a major merger, hence half of the galaxies in our stacked cluster would have gone through rapid isotropization. However, this fraction is an overestimation in our case, because we selected our cluster sample to be composed of regular clusters, thus avoiding clusters that have gone through recent major mergers – although they may have suffered major mergers in the fairly recent past. We further discuss irregular clusters in Sect. 6.3. Finally, the inner isotropy of galaxy orbits may be the consequence of the artificial phase mixing that is inherent in our stacked cluster, although many studies of individual clusters also find isotropic inner orbits (Fig. 15).

6.2.6. Inner isotropy of spiral galaxies: a selection effect due to single orbits?

It is harder to explain why the inner orbits of spirals are isotropic. The time for spirals to morphologically transform into lenticulars or ellipticals should be at least as large as the quenching time for star formation, which is expected to be slow for massive spirals within clusters. This was first shown by Mahajan et al. (2011), who compared the distributions of galaxies in PPS with predictions from cosmological simulations, and concluded that star formation in infalling galaxies is only quenched around the time when these galaxies cross the virial radius of the cluster on their first passage out of the cluster. According to Fig. B1 of Tollet et al. (2017), this occurs ∼3 Gyr after pericenter or ∼4 Gyr after cluster entry, while Wetzel et al. (2013, in their Fig. 8) found 4.5 Gyr since cluster entry for our median stellar and halo masses. Therefore, the mean radial velocities of the spiral population should be near zero, and the MAMPOSSt analysis is based on a valid Jeans equation. Also, the inner anisotropy of the spiral population should be roughly as radial as the outer anisotropy, in contrast to the isotropic inner velocities that we find for the spirals, for all three stacks.

However, many of our galaxies are not so massive and may be quenched at pericenter. The simplest explanation for the isotropic inner orbits of spirals would then be that spiral morphologies are destroyed at or before their first pericentric passage in the cluster. But there may not be sufficient time for the cluster to alter the morphologies of infalling spirals. Indeed, if ram pressure stripping is at the origin of morphological transformation of spirals into S0s (with depleted disks), the timescale for such morphological transformation should be at least the gas consumption time, which is typically 2 Gyr (Bigiel et al. 2011). If, instead intermediate mergers are the cause of transforming spirals into S0s (by bloated bulges), the timescale for violent relaxation during the merger should be on the order of a few internal galaxy crossing times, roughly 1 Gyr. This is comparable to the time of ∼1.3 Gyr from entry through the cluster virial radius to the first pericenter (see Fig. B1 of Tollet et al. 2017). It is difficult to imagine that spirals begin their morphological transformation as soon as their first entry into the cluster virial sphere.

On the other hand, spiral morphologies may be transformed before their return into the cluster on their second passage. Given the typical orbital times of 4 Gyr today (see Sect. 6.2.5), and that orbital times scale as the age of the Universe, galaxies that reach their second pericenter today have had an orbit lasting ∼3 Gyr. Therefore, present-day spirals should have time to complete their morphological transformation in a single orbit.

If spirals orbit only once around the clusters with their original morphology, their range of log pericenters will be much narrower than if they orbit many times. Those that fell in the cluster long ago at very small pericenters (given their small apocenters, as was the virial radius of the cluster’s most massive progenitor) should have different morphologies now¹¹. When the range of (log) pericenters is wide, the velocity anisotropy at a given radius r is dominated by the orbits with pericenters much smaller than r, which are near radial at r. But in the limit of a unique pericenter, the velocity anisotropy of spirals would be fully tangential (i.e., circular) at r = r_peri, rapidly increasing with radius to the radial values caused by infall (if the apocenters were all equal, one would return to circular at r = r_apo). This rapid transition in the velocity anisotropy profile is seen in recent hydrodynamical simulations (Lotz et al. 2019). But, as illustrated in Fig. 16, if the spirals only orbit once through the cluster, the radial orbits at r contribute less in comparison to the quasi-circular obits for r ≳ r_peri, leading to a more isotropic velocity distribution at r.

Therefore, while early-type galaxies may owe their inner isotropic orbits to isotropization from violent relaxation, as well as df and tidal braking (all over a few orbits), the quasi-isotropic orbits in late-type galaxies may simply be a selection effect leading to a narrow range for their pericenters, thus missing the contribution of radial orbits at the radius of study.

At high redshift, the orbital times are shorter, while the morphological transformation times should be roughly the same. Hence high redshift spirals may survive several orbits in their clusters, and we would then predict that the inner velocity anisotropy of spirals in high-redshift clusters will be somewhat radial. This is indeed found in the mass-orbit analysis of Biviano & Poggianti (2009). However, one could argue that at high redshift, regular quasi-spherical clusters do not yet exist and galaxy motions are set by the more filamentary geometry of proto-clusters.

6.2.7. Somewhat tangential inner orbits of ellipticals: tidal selection effects from BCGs?

There are signs of preferentially tangential inner orbits for the elliptical population (Figs. 7, 8, and 15). Admittedly, the evidence is weak (see Table 4) and is only seen for the sigv stack. This is in agreement with the evidence from a much smaller galaxy system, the Milky Way. Indeed, using proper motions from the second release of the Gaia astrometric mission, the 3D orbits of the dwarf spheroidals of the Milky Way have been reconstructed. The left panel of Fig. D.3 from Gaia Collaboration (2018) indicates that the apocenters of all nine dwarfs, except Leo II, are less than three times their pericenters. This implies somewhat tangential velocities, as demonstrated in Appendix E, where we found ⟨r_apo/r_peri⟩ = 3.79 ± 0.01 for an isotropic Hernquist model (the differences with the NFW model are negligible within the radial extent of the dwarf spheroidals).

The somewhat tangential orbits of ellipticals may be caused by those with small pericenters being tidally stripped by the BCG to the point that their masses fall below the limit of our sample (see Annunziatella et al. 2016), which we call tidal selection. For example, returning to the Milky Way, globular clusters, which are much more compact than dwarf spheroidals, have much more elongated 3D orbits (Gaia Collaboration 2018), which presumably are also caused by less tidal selection than for dwarf spheroidals.

6.2.8. Orbits of S0 galaxies: a clue to their formation?

We finally discuss whether the velocity anisotropy profiles of S0 galaxies provide clues to their formation. Our analysis points to the S0 population having intermediate outer orbits relative to the ellipticals and spirals for the sigv stack. This is not only seen in our complex models (Figs. 7, 8, and 15 and Table 4), but also from the simpler models, where Bayesian evidence favors radial outer orbits for the spirals, but marginally disfavors S0 outer anisotropy, while it strongly disfavors outer anisotropy of the ellipticals. A closer look at Figs. 8 and 15 reveals that, for the sigv stack, the positions of S0s in outer vs. inner velocity anisotropy space lie closer to those of the spirals, in particular because of signs of tangential inner anisotropy of the ellipticals.

But the position of S0s relative to the E and S galaxies depends on the stack. As seen in Fig. 14, which shows model 1, in the Num stack S0s display isotropic outer orbits, the ellipticals prefer slightly tangential outer orbits (but which remain consistent with isotropy), while the spirals show mildly radial outer orbits, and the S0s appear to lie somewhat closer to the ellipticals than to the spirals. On the other hand, in the tempX stack, the S0s appear to lie closer to the spirals, in particular for their outer orbits.

In fact, the 68% contours of the S0s fully encompass those of the spirals for the sigv and tempX stacks, but this is not seen in the Num stack, nor between the other types in any stack. We quantified the correspondence of the velocity anisotropies between morphological types of a given stack by computing the Pearson correlation coefficient C_1, 2 between two types as

$$\begin{array}{c}\hfill C_{1,2}=\frac{\langle (f_1-\langle f_1\rangle )(f_2-\langle f_2\rangle )\rangle }{\sigma (f_1)\sigma (f_2)},\end{array}$$(27)

where f_i is the fraction of MCMC points (after burn-in) for type i that lie in a given cell of β_sym(r₂₀₀) vs. β_sym(0.3 r₂₀₀) of width Δβ_sym = 0.03 and Eq. (27) is estimated over the entire set of cells. The Pearson coefficients between pairs of morphological types, displayed in Table 5, for model 1, imply that S0 galaxies have orbits closer to the spirals for the sigv and tempX stacks, while their orbits are closer to ellipticals for the Num stack. The orbits of spirals and ellipticals are the furthest apart for all three stacks.

The difference in inner anisotropy between S0s and ellipticals in the sigv stack may be the consequence of the possible shorter timescale for transformation from S0 to E than for tidal stripping, in which case the BCGs would have less time to tidally strip the S0s than the ellipticals, and therefore ellipticals on radial orbits are more subject to the selection effect against radial inner orbits (Sect. 6.2.7) than are S0s.

S0s appear to form relatively late in clusters, the fraction of S0s is 0.1 at z ∼ 1 (Smith et al. 2005; Postman et al. 2005), that is 8 Gyr ago, while our sample has ∼40% of S0s (see Table 1). Roughly speaking, one then expects that the great majority of the S0s were formed over 2 Gyr ago, which is the rough time scale for S→S0 evolution. One could then argue that orbital isotropization should take longer than the S→S0 evolution, otherwise S0 orbits would have mostly isotropized and resemble those of ellipticals more than those of spirals. Finally, our results on S0s should be taken with caution, because the S0 morphological class is notoriously difficult to cleanly classify for ranges of inclinations and apparent magnitudes (e.g., van den Bergh 2009) and may thus be contaminated by both spirals and ellipticals.

6.3. Perspectives

A larger dataset is needed to obtain better constraints on the orbital anisotropy of galaxies of different morphological types. This dataset should be more complete in projected radius and in stellar mass. This ought to lead to smaller differences in the orbital anisotropies inferred using different stacking methods.

It would also be interesting to compare our results for regular clusters with those for irregular clusters. Indeed, since irregular clusters are merging galaxy systems, the orbits of galaxies are affected by the changes in the gravitational potential. This should cause violent relaxation, leading to orbit isotropization, which should occur on the timescale of the cluster-cluster merger. This is on the order of the orbital timescale for galaxies that have recently entered their cluster. According to Paper I, the scale radii of the ellipticals and S0s in irregular clusters are double those in (comparable mass) regular clusters, while the scale radii of the spiral population is somewhat smaller than in regular clusters. Thus, while in regular clusters spirals had scale radii 4 (3) times greater than those of ellipticals (respectively S0s), the scale radii of the three morphological types are much more similar in irregular clusters. Perhaps the violent relaxation in irregular clusters perturbs the orbits more efficiently as compared to how the different histories of the three morphological types differentiate the orbits.

7. Conclusions

We ran the MAMPOSSt mass-orbit modeling algorithm on three stacks of 54 regular, nearby (z   ≈   0.05) clusters from the WINGS sample, composed of up to 4682 galaxies located between 0.03 r₂₀₀ and r₂₀₀, split between ellipticals, lenticulars, and spirals (including irregulars). MAMPOSSt is a Bayesian method that jointly fits the distribution of galaxies of these three morphological types in projected phase space, fitting for the shape of the total mass profile on one hand and of the three velocity anisotropy profiles on the other. We ran MAMPOSSt with 30 different sets of priors. Our results for the sigv stack are as follows:

• There is no compelling evidence for a mass density profile steeper than NFW or n = 6 Einasto at 0.03 r₂₀₀ (in fact a cored-NFW profile is only weakly rejected), even though our highest likelihoods are reached with total density profiles that are steeper than NFW (with an inner slope of roughly −1.5 ± 0.5) or with an NFW profile for the BCG in addition to an NFW profile for the remaining part of the cluster.

• An outer slope as steep as −4 (Hernquist model) is ruled out.

• The concentration of the mass distribution, when set free, is consistent with those in massive halos within dissipationless cosmological N-body simulations as well as with those measured in similar-mass clusters using weak gravitational lensing.

• The number density profile of elliptical galaxies traces the total mass density profile very well, while that of S0s only does so marginally and the spiral profile clearly does not.

• The velocity anisotropy of spirals rises from near isotropic in the inner regions to mildly radial (β ≃ 0.45 ± 0.08) at r₂₀₀.

• The velocity anisotropy of the lenticulars also rises from near isotropic in the inner regions to somewhat less radial (β = 0.31 ± 0.17) at r₂₀₀ than for the spirals.

• The velocity anisotropy of the ellipticals is consistent with isotropic velocities everywhere, even though the highest likelihoods are reached for slightly tangential inner orbits and mildly radial anisotropy (β = 0.19 ± 0.25) at r₂₀₀.

• BIC Bayesian evidence (which prefers simpler models) favors isotropy everywhere for the ellipticals and S0s, but does not strongly reject having outer radial anisotropy for the S0s. This is actually the preferred model using AIC evidence.

• Bayesian evidence (both BIC and AIC) suggests that the anisotropy radius (transitioning from the lowest to highest values) is not different from the scale radius of the considered morphological type.

• For simple priors, Bayesian evidence favors mild increases to the velocity anisotropy (T model) compared to the sharp increase of the generalized Osipkov-Merritt model. For complex priors, the two models lead to similar likelihoods.

Some of these conclusions are marginally different for the other two stacks:

• There is marginal evidence for a steeper inner mass density than NFW with the Num and tempX stacks.

• The outer anisotropy of spirals is less pronounced.

• The outer orbits of E and S0 galaxies are consistent with being isotropic for Num, and ellipticals also show isotropic outer orbits for tempX, while they are moderately radial with sigv (although not favored by Bayesian evidence).

• There is no weak evidence of tangential inner orbits for ellipticals.

• S0 orbits resemble more those of spirals for the sigv and tempX stacks and more those of ellipticals for the Num stack.

The velocity anisotropies of the three morphological types provide important clues to their transformations as they orbit clusters.

The very large radial extent of spiral galaxies suggests that they are infalling. Such infall should lead to fairly radial outer orbits for spirals (as seen in the sigv stack, but less so in the other two). Near r₂₀₀, E and S0 galaxies should be a mix of the virialized (i.e., isotropized) population and the infalling members, hence with less radial orbits than the spirals. The inner isotropy of the early-type galaxies cannot be produced by two-body relaxation, which is too slow. One possibility is that inner isotropy of the E and S0 galaxies is the consequence of violent relaxation occurring during major mergers of clusters, which appear to occur at a sufficient rate. Alternatively, galaxies may lose their orbital energy by a combination of dynamical friction and tidal braking suffered by the host groups that they may live in.

The inner isotropy of spirals cannot be explained in this manner, because spirals should be transformed into S0s over an orbital time (as confirmed by their much wider spatial distribution). If spiral galaxies only pass once through pericenter, there is a selection against radial orbits at a given small distance to the cluster center, explaining their quasi-isotropic inner orbits. Finally, although only marginally significant in the sigv stack and not in the others, we conjecture that the possible tangential anisotropy of the ellipticals may be caused by tidal selection where those on small pericenters are tidally stripped and fall below the sample mass threshold.

---

Acknowledgments

GAM thanks Radek Wojtak and Stefano Ettori for useful comments and Clifford M. Hurvich for statistical advice. We also thank the anonymous referee for her/his constructive and enlightening comments. AB is grateful to the IAP for its hospitality. The work of AC was supported by the STARFORM Sinergia Project funded by the Swiss National Science Foundation. GAM’s version of the MAMPOSSt code (the version used in this work) is publicly available at https://gitlab.com/gmamon/MAMPOSSt.

References

Appendix A: Equivalence of the Solanes & Salvador-Solé (1990) and Aguerri et al. (2017) anisotropy inversions

In their Sect. 2.6, Aguerri et al. (2017, hereafter, A+17) present a formula (but do not derive it) for anisotropy inversion that appears similar to that of SS90, partially using the notations of SS90, but without referring to these authors where they present their equations. In this appendix, we confirm that the two formulae are indeed equivalent and provide a more robust version of it.

Following the notations of mass inversion study of Mamon & Boué (2010), we express the radial dynamical pressure as $p=\nu \langle v_r^2\rangle$, and introduce a projected pressure of $P=\mathrm{\Sigma }{\sigma }_{\mathrm{los}}^2$ (which does not have the dimension of pressure, but that of pressure times size). We denote $q=\nu v_c^2$, which has the dimension of pressure. The H function introduced by both SS90 and A+17 is equal to P/2, while Ψ of A+17 (Ψ₁ of SS90) satisfies Ψ(r) = −q/r.

The kinematic projection equation (Binney & Mamon 1982)

$$\begin{array}{c}\hfill P(R)=2{\displaystyle {\int }_R^{\infty }\frac{r}{\sqrt{r^2-R^2}}[1-\beta (r){\left(\frac{R}{r}\right)}^2]p(r)\mathrm{d}r}\end{array}$$(A.1)

can be inverted for general anisotropy (Mamon & Boué 2010; Wolf et al. 2010). For isotropy (β = 0) this simply involves the standard Abel inversion applied to determine p from P (instead of ν from Σ) to yield

$$\begin{array}{c}\hfill p(r)\equiv p_{\mathrm{iso}}(r)=-\frac{1}{\pi }{\displaystyle {\int }_r^{\infty }\frac{P^{\prime }(R)}{\sqrt{R^2-r^2}}\mathrm{d}R=-\frac{1}{\pi r}\frac{\mathrm{d}I}{\mathrm{d}r},}\end{array}$$(A.2)

where (using the notation of SS90)

$$\begin{array}{c}\hfill I(r)={\displaystyle {\int }_r^{\infty }\frac{R}{\sqrt{R^2-r^2}}P(R)\mathrm{d}R}\end{array}$$(A.3)

(see Eq. (22) of Mamon & Boué 2010, who denoted I by J, while A+17 used K in place of I).

Both SS90 and A+17 provide two equations, one for β p and one for (3 − 2β)p. They can then deduce β(r) from the ratio of the two equations (thus eliminating p). The first of the two final equations of SS90 (their Eq. (23)) and A+17 (their Eq. (10)) are clearly the same apart from notation. In our notation, this is

$$\begin{array}{c}\hfill [3-2\beta (r)]p(r)={\displaystyle {\int }_r^{\infty }\frac{q(s)}{s}\mathrm{d}s+2p_{\mathrm{iso}}(r).}\end{array}$$(A.4)

The second of the final equations of SS90 and A+17 are also clearly the same, and, transformed to our notation, are

$$\begin{array}{cc}\hfill \beta (r)p(r)=& -\frac{1}{r^3}{\displaystyle {\int }_0^r{s}^{2}q(s)\mathrm{d}s-p_{\mathrm{iso}}(r)-\frac{3}{\pi r^2}I(r)}\hfill \\ \hfill & +\frac{3}{\pi r^3}{\displaystyle {\int }_0^{r}I(R)\mathrm{d}R.}\hfill \end{array}$$(A.5)

However, the last term of Eq. (A.5) is an integral of an integral, which can be transformed to

$$\begin{array}{cc}\hfill {\displaystyle {\int }_0^{r}I(s)\mathrm{d}s=}& {\displaystyle {\int }_0^r\mathrm{d}s{\int }_s^{\infty }\frac{R}{\sqrt{R^2-s^2}}P(R)\mathrm{d}R}\hfill \\ \hfill =& \frac{\pi }{2}{\displaystyle {\int }_0^{r}RP(R)\mathrm{d}R+{\int }_r^{\infty }{sin}^{-1}\left(\frac{r}{R}\right)RP(R)\mathrm{d}R,}\hfill \end{array}$$

where the last integral is obtained after reversing the order of integration in the double integral. It thus seems preferable to rewrite Eq. (A.5) as

$$\begin{array}{cc}\hfill {(\beta p)}_{\mathrm{new}}(r)=& -\frac{1}{r^3}{\displaystyle {\int }_0^r{s}^{2}q(s)}\mathrm{d}s-p_{\mathrm{iso}}(r)\hfill \\ \hfill & +\frac{3}{2r^3}{\displaystyle {\int }_0^{r}RP(R)\mathrm{d}R}\hfill \\ \hfill & +\frac{3}{\pi r^3}{\displaystyle {\int }_r^{\infty }[{sin}^{-1}\left(\frac{r}{R}\right)-\frac{r}{\sqrt{R^2-r^2}}]RP(R)\mathrm{d}R.}\hfill \end{array}$$(A.6)

Appendix B: Equivalence of the Bicknell et al. (1989) and Dejonghe & Merritt (1992) anisotropy inversions

Bicknell et al. (1989, hereafter B+89) and Dejonghe & Merritt (1992, hereafter DM92) both first computed the radial dynamical pressure to then solve the Jeans equation for the anisotropy profile. Using the same notations as in Appendix A, we can express the pressure of B+89 (their Eqs. (3.5)–(3.8)) as

$$\begin{array}{cc}\hfill p_{\mathrm{B}+89}(r)=& \frac{1}{3}{\displaystyle {\int }_r^{\infty }\frac{q(s)}{s}\mathrm{d}s-\frac{2}{3r^3}{\int }_0^r{s}^{2}q(s)\mathrm{d}s}\hfill \\ \hfill & +\frac{1}{r^3}{\displaystyle {\int }_0^{r}RP(R)\mathrm{d}R}\hfill \\ \hfill & -\frac{2}{\pi r^3}{\displaystyle {\int }_r^{\infty }[\frac{r}{\sqrt{R^2-r^2}}-{sin}^{-1}\left(\frac{r}{R}\right)]RP(R)\mathrm{d}R,}\hfill \end{array}$$(B.1)

where we assumed that their $A=(3/2){\int }_0^{\infty }RP(R)\mathrm{d}R/{\int }_0^{\infty }q{r}^2\mathrm{d}r$ is unity, as expected from virial equilibrium (see the discussion by DM92)¹². With our notation, the radial pressure of DM92 (their Eq. (57a)) is

$$\begin{array}{cc}\hfill p_{\mathrm{DM}92}(r)& =\frac{1}{3}{\displaystyle {\int }_r^{\infty }\frac{q(s)}{s}\mathrm{d}s+\frac{2}{3r^3}{\int }_r^{\infty }{s}^{2}q(s)\mathrm{d}s}\hfill \\ \hfill & -\frac{2}{\pi r^3}{\displaystyle {\int }_r^{\infty }[\frac{r}{\sqrt{R^2-r^2}}+{cos}^{-1}\left(\frac{r}{R}\right)]RP(R)\mathrm{d}R}\hfill \end{array}$$(B.2)

where we assumed again the relation A = 1 (which translates to G(∞) = 0 in DM92’s notation). Replacing sin⁻¹ by π/2 − cos⁻¹ x, one can re-express the pressure of B+89 of Eq. (B.1) as

$$\begin{array}{c}{p}_{\text{B+89}}(r)=\frac{1}{3}{\displaystyle {\int }_r^{\infty }\frac{q(s)}{s}}\text{d}s-\frac{2}{3r^3}{\displaystyle {\int }_0^{\infty }{s}^2}q(s)\text{d}s\\ +\frac{2}{3r^3}{\displaystyle {\int }_r^{\infty }{s}^2}q(s)\text{d}s+\frac{1}{r^3}{\displaystyle {\int }_0^{\infty }R}P(R)\text{d}R\\ -\frac{2}{\pi r^3}{\displaystyle {\int }_r^{\infty }\left[\frac{r}{\sqrt{R^2-r^2}}+{\cos }^{-1}\left(\frac{r}{R}\right)\right]}RP(R)\text{d}R\\ =p_{\text{DM92}}(r)\\ +\frac{1}{6\pi r^3}\left[{\displaystyle {\int }_0^{\infty }6}\pi RP(R)\text{d}R-{\displaystyle {\int }_0^{\infty }4}\pi r^{2}q(r)\text{d}r\right].\end{array}$$(B.3)

Now, the term in brackets in Eq. (B.3) is zero from the virial theorem, as the first integral represents the kinetic energy, while the second integral represents the absolute value of the potential energy (see Eq. (51) of DM92). This can be checked by inserting P(R) from the equation of kinematic projection (A.1) into the first integral, yielding (after inversion of the order of the integrals)

$${\displaystyle {\int }_0^{\infty }RP(R)\mathrm{d}R=\frac{2}{3}{\int }_0^{\infty }[3-2\beta (r)]p(r)r^2\mathrm{d}r}$$

and by inserting q(r) from the Jeans equation, which in our simplified notation is

$$\begin{array}{c}\hfill \frac{\mathrm{d}p}{\mathrm{d}r}+2\beta \frac{p(r)}{r}=-\frac{q(r)}{r},\end{array}$$(B.4)

yielding (again after inversion of the order of integration)

$${\displaystyle {\int }_0^{\infty }{r}^{2}q(r)\mathrm{d}r={\int }_0^{\infty }[3-2\beta (r)]p(r)r^2\mathrm{d}r.}$$

Therefore, the B+89 and DM92 anisotropy inversions are equivalent, as expected. The anisotropy is obtained by solving the Jeans Eq. (B.4) for the anisotropy, yielding

$$\begin{array}{c}\hfill \beta (r)=-\frac{r{p}^{\prime }(r)+q(r)}{2p(r)}=-\frac{1}{2}[\frac{q(r)}{p(r)}+\frac{\mathrm{d}lnp}{\mathrm{d}lnr}].\end{array}$$(B.5)

Thus both B+89 and DM92 algorithms involve differentiation of the data, as seen from the derivative p′(r) in Eq. (B.5) and the data terms involving P(R) in Eqs. (B.1) and (B.2). The DM92 algorithm seems preferable, as it is simpler and avoids inward extrapolation of the data to R = 0¹³.

The derivative of the radial pressure appearing in Eq. (B.5) can be written as an integral involving the derivative of an analytical fit to the observed projected pressure, as follows. Differentiating Eq. (B.2) leads to

$$\begin{array}{c}\hfill p\prime (r)=-\frac{q(r)}{r}-\frac{2}{\pi }\frac{\mathrm{d}(J/r^3)}{\mathrm{d}r}=-\frac{q(r)}{r}+\frac{6}{\pi r^4}J(r)-\frac{2}{\pi r^3}\frac{\mathrm{d}J}{\mathrm{d}r},\end{array}$$(B.6)

where

$$\begin{array}{c}\hfill J(r)={\displaystyle {\int }_r^{\infty }[{cos}^{-1}\left(\frac{r}{R}\right)+\frac{r}{\sqrt{R^2-r^2}}]RP(R)\mathrm{d}R.}\end{array}$$(B.7)

With the substitution R = r/cosθ, we can write

$$\begin{array}{c}\hfill J(r)=r^2{\displaystyle {\int }_0^{\pi /2}(\theta +cot\theta )P\left(\frac{r}{cos\theta }\right)\frac{sin\theta }{{cos}^3\theta }\mathrm{d}\theta ,}\end{array}$$(B.8)

which leads to

$$\begin{array}{c}\hfill \frac{\text{d}J}{\text{d}r}=2r{\displaystyle {\int }_0^{\pi /2}\left(\theta +\cot \theta \right)}P\left(\frac{r}{\cos \theta }\right)\frac{\sin \theta }{{\cos }^3\theta }\text{d}\theta \\ \hfill +r^2{\displaystyle {\int }_0^{\pi /2}\left(\theta +\cot \theta \right)}{P}^{\prime }\left(\frac{r}{\cos \theta }\right)\frac{\sin \theta }{{\cos }^4\theta }\text{d}\theta \end{array}$$(B.9)

$$\begin{array}{c}=\frac{2}{r}{\displaystyle {\int }_r^{\infty }\left[{\cos }^{-1}\left(\frac{r}{R}\right)+\frac{r}{\sqrt{R^2-r^2}}\right]}RP(R)\text{d}R\\ +\frac{1}{r}{\displaystyle {\int }_r^{\infty }\left[{\cos }^{-1}\left(\frac{r}{R}\right)+\frac{r}{\sqrt{R^2-r^2}}\right]}{R}^2{P}^{\prime }(R)\text{d}R.\end{array}$$(B.10)

Inserting the expression of dJ/dr of Eq. (B.10) into the second equality of Eq. (B.6) yields

$$\begin{array}{cc}\hfill rp\prime (r)+q(r)& =\frac{2}{\pi r^3}{\displaystyle {\int }_r^{\infty }[{cos}^{-1}\left(\frac{r}{R}\right)+\frac{r}{\sqrt{R^2-r^2}}]}\hfill \\ \hfill & \times R[P(R)-R{P}^{\prime }(R)]\mathrm{d}R.\hfill \end{array}$$(B.11)

Inserting equation (B.11) into the first equality of equation (B.5) then yields the anisotropy

$$\beta (r)=\frac{A}{B},$$

where

$$\begin{array}{ccc}A& =\hfill & {\displaystyle {\int }_r^{\infty }[{cos}^{-1}\left(\frac{r}{R}\right)+\frac{r}{\sqrt{R^2-r^2}}]R[P(R)-R{P}^{\prime }(R)]\mathrm{d}R,}\\ B& =& 2{\displaystyle {\int }_r^{\infty }[{cos}^{-1}\left(\frac{r}{R}\right)+\frac{r}{\sqrt{R^2-r^2}}]RP(R)\mathrm{d}R}\hfill \\ & & -\pi \frac{r^2}{3}{\displaystyle {\int }_r^{\infty }q(s)(\frac{r}{s}+2\frac{s^2}{r^2})\mathrm{d}s.}\hfill \end{array}$$

Appendix C: Additional figure

Appendix D: Two-body relaxation time for galaxies falling into clusters

D.1. Transverse velocity formalism

We consider equal-mass galaxies of total (i.e., subhalo) mass m infalling into a pre-existing spherical NFW cluster of scale radius a. For simplicity, we will (incorrectly) assume that the cluster is stationary. We also assume that the galaxy number density profile is proportional to the cluster mass density profile, and that the galaxy masses are constant, unaffected by tidal stripping.

The change of velocity of a galaxy in a single encounter can be written as (Eq. (4–2) of Binney & Tremaine 1987)

$$\begin{array}{c}\hfill |\delta v_{\perp }|=2\frac{\mathit{Gm}}{p{v}_{\mathrm{rel}}},\end{array}$$(D.1)

where v_rel is the relative velocity of the two galaxies, p is the impact parameter, and G is Newton’s gravitational constant. During a time dt, the galaxy will suffer a number of encounters with other galaxies with impact parameter between p and p + dp of

$$\begin{array}{c}\hfill {\mathrm{d}}^2\mathcal{N}=(\nu v_{\mathrm{rel}}\mathrm{d}t)\times (2\pi p\mathrm{d}p).\end{array}$$(D.2)

Equations (D.1) and (D.2) lead to a change in squared transverse velocity over an orbit of

$$\begin{array}{c}\hfill \mathrm{\Delta }{v}_{\perp }^2={\displaystyle {\int }_{\mathrm{orbit}}\mathrm{d}t{\int }_{p_{\min }}^{p_{\max }}\frac{{\mathrm{d}}^2\mathcal{N}}{\mathrm{d}p\mathrm{d}t}\delta v_{\perp }^2\mathrm{d}p=8\pi G^2{m}^2{\int }_{\mathrm{orbit}}\frac{\nu (t)}{v_{\mathrm{rel}}(t)}ln\mathrm{\Lambda }(t)\mathrm{d}t,}\end{array}$$(D.3)

where Λ = p_max/p_min. One generally assumes $p_{\min }=Gm/v_{\mathrm{rel}}^2$ and p_max = r, yielding $\mathrm{\Lambda }=r{v}_{\mathrm{rel}}^2/(Gm)$.

We now need to write the equation of motion of infalling galaxies to change the outer integration variable from time to radius.

D.2. Equation of motion of orbit of known pericenter and apocenter

We begin by expressing the conservation of energy E and angular momentum L:

$$\begin{array}{cc}& E=\frac{1}{2}{v}^2+\mathrm{\Phi }(r)=\frac{1}{2}{v}_{\mathrm{peri}}^2+\mathrm{\Phi }(r_{\mathrm{peri}})=\frac{1}{2}{v}_{\mathrm{apo}}^2+\mathrm{\Phi }(r_{\mathrm{apo}}),\hfill \end{array}$$(D.4)

$$\begin{array}{cc}& L^2=r^2{v}_{\mathrm{t}}^2=r_{\mathrm{peri}}^2{v}_{\mathrm{peri}}^2=r_{\mathrm{apo}}^2{v}_{\mathrm{apo}}^2,\hfill \end{array}$$(D.5)

where “peri” and “apo” respectively denote the pericenter and apocenter of the orbit, while $v_{\mathrm{t}}^2$ is the squared tangential velocity at any time. Applying $v^2={\dot{r}}^2+v_{\mathrm{t}}^2$ to the first two equalities of Eqs. (D.4) and (D.5), we deduce the equation of motion

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}r}{\mathrm{d}t}\right)}^2-[1-{\left(\frac{r_{\mathrm{peri}}}{r}\right)}^2]v_{\mathrm{peri}}^2-2[\mathrm{\Phi }(r_{\mathrm{peri}})-\mathrm{\Phi }(r)]=0.\end{array}$$(D.6)

Eliminating v_apo in the third equalities of Eqs. (D.4) and (D.5) yields

$$\begin{array}{c}\hfill v_{\mathrm{peri}}^2=\frac{2}{1-1/{\lambda }^2}[\mathrm{\Phi }(\lambda r_{\mathrm{peri}})-\mathrm{\Phi }(r_{\mathrm{peri}})],\end{array}$$(D.7)

where λ = r_apo/r_peri. Combining Eqs. (D.6) and (D.7) leads to an equation of motion expressed in terms of the pericenter and the ratio of apo- to pericenter:

$$\begin{array}{cc}\hfill {\left(\frac{\mathrm{d}r}{\mathrm{d}t}\right)}^2& -\frac{2}{1-1/{\lambda }^2}[1-{\left(\frac{r_{\mathrm{peri}}}{r}\right)}^2]\mathrm{\Phi }(\lambda r_{\mathrm{peri}})\hfill \\ \hfill & -2\{1-\frac{1}{1-1/{\lambda }^2}[1-{\left(\frac{r_{\mathrm{peri}}}{r}\right)}^2]\}\mathrm{\Phi }(r_{\mathrm{peri}})+2\mathrm{\Phi }(r)=0.\hfill \end{array}$$(D.8)

D.3. Change of transverse velocity over an orbit

Expressing dt = ϵ dr/|dr/dt|, where ϵ = −1 and 1 for galaxies falling in and bouncing out, respectively, we can express the squared change in transverse velocity over an orbit as

$$\begin{array}{c}\hfill \mathrm{\Delta }{v}_{\perp }^2=16\pi G^2{m}^2{\displaystyle {\int }_{r_{\mathrm{peri}}}^{r_{\mathrm{apo}}}\frac{\nu (r)}{v_{\mathrm{rel}}(r)}ln\left(\frac{r{v}_{\mathrm{rel}}^2(r)}{Gm}\right)\frac{\mathrm{d}r}{\left|\mathrm{d}r/\mathrm{d}t\right|},}\end{array}$$(D.9)

where |dr/dt| is deduced from Eq. (D.8). Equation (D.9) can be rewritten in dimensionless form by dividing the squared velocities by the circular velocity $\sqrt{GM(a)/a)}$ at the scale radius a and using Eq. (7), which leads to

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{v}_{\perp }^2}{GM(a)/a}=4N(a){\left(\frac{m}{M(a)}\right)}^2{\displaystyle {\int }_{r_{\mathrm{peri}}/a}^{r_{\mathrm{apo}}/a}\frac{\stackrel{\sim }{\nu }(x)}{y_{\mathrm{rel}}(x)}ln\left(\frac{x{y}_{\mathrm{rel}}^2}{m/M(a)}\right)\frac{\mathrm{d}x}{|y_r|},}\end{array}$$(D.10)

where N(a) is the number of galaxies in the sphere of radius a, x = r/a, $y_{\mathrm{rel}}=v_{\mathrm{rel}}/\sqrt{GM(a)/a}$, and $y_r=\mathrm{d}r/\mathrm{d}t/\sqrt{GM(a)/a}$, noting that the NFW gravitational potential can be written as

$$\begin{array}{cc}\hfill \mathrm{\Phi }(r)& =-\frac{GM(a)}{a}\stackrel{\sim }{\mathrm{\Phi }}\left(\frac{r}{a}\right)\hfill \end{array}$$(D.11)

$$\begin{array}{cc}\hfill \stackrel{\sim }{\mathrm{\Phi }}(x)& =\frac{1}{ln2-1/2}\frac{ln(1+x)}{x}·\hfill \end{array}$$(D.12)

Since the typical velocities of galaxies in clusters are of the order of the circular velocity at the scale radius, Eq. (D.10) directly provides the inverse of the typical number of orbits for a galaxy to isotropize by two-body relaxation.

Infalling galaxies encounter outgoing galaxies as well as virialized galaxies and other infalling galaxies. Thus, the relative velocities will, on average, be greater than the velocities of the test galaxies, that is |v_rel|>|v|. Hence, assuming v_rel = v, meaning that the galaxies that the test galaxies encounter are static, produces an upper limit to the amount of two-body relaxation in pumping angular momentum into the infalling galaxies, thus building up $\mathrm{\Delta }{v}_{\perp }^2$. We thus write

$$\begin{array}{c}\hfill y_{\mathrm{rel}}=\sqrt{y_r^2+\frac{2}{1-1/{\lambda }^2}{\left(\frac{x_{\mathrm{peri}}}{x}\right)}^2[\stackrel{\sim }{\mathrm{\Phi }}(x_{\mathrm{peri}})-\stackrel{\sim }{\mathrm{\Phi }}(\lambda x_{\mathrm{peri}})]},\end{array}$$(D.13)

where we used x_peri = r_peri/a, the second equality in Eq. (D.5) and Eq. (D.11).

According to our MAMPOSSt fits of the distribution of E, S0, and S galaxies in PPS, we have for model 24, r₂₀₀ = 1.7 Mpc, c₂₀₀ = 3.9, leading to M₂₀₀ = 10^14.7 M_⊙, and a = 437 kpc. Solving Eq. (25) for N(a) = N(r_ν) yields N(a) = 13.8. Given the mean stellar mass of 10¹⁰ M_⊙, we deduced the subhalo mass m by solving the stellar mass as a function of halo mass from the abundance matching formula of Behroozi et al. (2013). For our mean stellar mass of 10¹⁰ M_⊙ (Sect. 3.1), this yields a subhalo mass of m = 4.6 × 10¹¹ M_⊙ at z = 0 and 5.4 × 10¹¹ M_⊙ at z = 1. We take the higher subhalo mass to be conservative. This in turns gives a mass ratio of m/M(a) = 0.004.

Figure D.1 provides the log number of orbits to isotropize, $[GM(a)/a]/\mathrm{\Delta }{v}_{\perp }^2$, whose inverse is given in Eq. (D.10), given our estimates of m/M(a) and N(a) for our cluster sample and its associated galaxies. The figure indicates that two-body relaxation is slow for the bulk of our galaxies, with 10^1.5 ≃ 30 orbits required for deeply penetrating orbits with λ = 5 (the typical value of apo- to pericenter ratio found in ΛCDM halos, Ghigna et al. 1998), and even more for orbits with greater pericenters. In comparison, a naïve application of the N/(8 ln N) formula (Eq. (4)–(9) of Binney & Tremaine 1987), would yield only 5.7 orbits at r = a for N = M(a)/m = 1/0.004 = 250.

More massive galaxies do not relax faster with Eq. (D.10), because they encounter too few galaxies of comparable mass to isotropize. Recall that this calculation is conservative as many of the encounters of infalling galaxies involve galaxies moving in the other directions with relative velocities that are double the velocity of the test galaxy.

Appendix E: Apocenter-to-pericenter ratio for isotropic Hernquist models

Given a system of particles orbiting in a fixed gravitational potential, one can determine the pericentric and apocentric radii by expressing the conservation of energy and angular momentum. The first equalities of Eqs. (D.4) and (D.5) imply that the pericenter and apocenter are the roots of

$$\begin{array}{c}\hfill \frac{1}{2}\frac{L^2}{r^2}+\mathrm{\Phi }(r)-E=0.\end{array}$$(E.1)

We built an isotropic Hernquist (1990) model following the method of Kazantzidis et al. (2004), where we first draw random radii, compute the gravitational potential, and then we draw velocities from f(v|r) = v²f(v²/2 + Φ(r), where f ≡ f(E) is the 6D distribution function of the isotropic Hernquist model, given by Hernquist (1990). We checked that the velocity anisotropy profile was near zero at all radii (median value of β = −0.01). Once we drew 100 000 6D coordinates, we solved Eq. (E.1), where the two roots correspond to the pericenter and apocenter.

The extraction of the roots of Eq. (E.1) for each of the 10⁵ particles was performed in vectorial fashion: for a list of 6001 geometrically spaced radii r_i between 0.001 and 1000 (in units of the Hernquist scale radius), we estimate the left-hand-side (LHS) of Eq. (E.1) with r = r_i. We first set r_peri = 0 and r_apo = 10⁶a (where a is the scale radius of the Hernquist model) for all the particles. Noting that the LHS of Eq. (E.1) must be less than or equal to 0, since it represents $-\frac{1}{2}{(\mathrm{d}r/\mathrm{d}t)}^2$, we then vectorially adjusted r_peri and r_apo with the conditions

where |LHS| < |oldLHS| & LHS < oldLHS
    update rperi
where |LHS| < |oldLHS| & LHS > oldLHS
    update rapo
save LHS to oldLHS

It took 15 (1.5) seconds to process 100 000 particles with 0.001 (respectively 0.01) dex precision in this manner with a script language (SM, also known as SuperMongo) on a single processor. We found a median r_apo/r_peri of 3.79 with an uncertainty of 0.02 (from ten trials).

All Tables

All Figures