Structural and dynamical modeling of WINGS clusters. II. The orbital anisotropies of elliptical, spiral and lenticular galaxies

The Bayesian MAMPOSSt mass/orbit modeling algorithm is used to jointly fit the distribution of elliptical, spiral (and irregular), and lenticular galaxies in projected phase space, on 3 stacked clusters (with normalized positions and velocities) of 54 regular clusters from the WINGS survey, with at least 30 member velocities. Our stacked clusters contain ~5000 galaxies with available velocities and morphological types. 30 runs of MAMPOSSt with different priors are presented. The highest MAMPOSSt likelihoods are obtained for generalized NFW models with steeper inner slope, free-index Einasto models, and double NFW models for the cluster and the brightest cluster galaxy. However, there is no strong Bayesian evidence for a steeper profile than the NFW model. The mass concentration matches the predictions from cosmological simulations. Ellipticals usually trace best the mass distribution, while S0s are close. Spiral galaxies show increasingly radial orbits at increasing radius, as do S0s on two stacks, and ellipticals on one stack. The inner orbits of all three types in the 3 stacks are consistent with isotropy. Spiral galaxies should transform rapidly into early-types given their much larger extent in clusters. Outer radial orbits are expected for the spirals, a consequence of their recent radial infall into the cluster. The less radial orbits we find for early-types could be related to the longer time spent by these galaxies in the cluster. We demonstrate that two-body relaxation is too slow to explain the inner isotropy of the early types, which suggests that inner isotropy is the consequence of violent relaxation during major cluster mergers or dynamical friction and tidal braking acting on subclusters. We propose that the inner isotropy of the short-lived spirals is a selection effect of spirals passing only once through pericenter before being transformed into early-type morphologies.


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-Send offprint requests to: G. A. Mamon, e-mail: gam@iap.fr 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 where ν(r) and M(r) are the radial profiles of respectively the tracer density and the total mass, v 2 r (r) is the mean squared radial velocity profile, and 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). For equilibrium systems, there are no net meridional streaming motions, hence v 2 θ = σ 2 θ . For galaxy clusters, one neglects rotation leading to v 2 φ = σ 2 φ . By symmetry, σ φ = σ θ . On the other hand, the JE contains radial streaming motions, e.g. from infall. Radial, isotropic and circular orbits respectively have β = 1, β = 0, and β → −∞.
In additional to such a Jeans analysis, another class of analysis uses the collisionless Boltzmann (Vlasov) equation (CBE), which states that the six-dimensional fluid is incompressible in the absence of galaxy collisions (see, e.g., chap. 5 of Courteau et al. 2014 for a quick review of these stellar dynamics). In fact, the JE is a direct consequence (1st 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 or in action angle space, 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 to extract orbital shapes from either Jeans or DF analyses. 1) Clusters are observed in projection, with only 2 positional (sky) coordinates and one (lineof-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 (radial) equation linking the unknown mass profile (M(r) = (r 2 /G) dΦ/dr) and anisotropy (linked to the orbital shapes), an issue called the mass / anisotropy degeneracy (MAD, Binney & Mamon 1982), while 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, focused 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. 1 A disadvantage of Jeans analysis methods is that they usually require the radial binning of the data. 2 One way around this is to specify the form of the LOS velocity distribution function, and starting with Walker et al. (2009), it is 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) assumed a separable form for f (E, L) and found that it matched well 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 . A promising method is to express the DF in terms of action-angle variables (Vasiliev 2019).
In a hybrid method called MAMPOSSt (Mamon, Biviano, & Boué 2013, hereafter MBB13), the DF is no longer expressed in terms of E and L, but in terms of the three-dimensional velocity distribution function, the simplest form being a Gaussian. 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 v 2 r (r) 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 of Sanchis, Łokas, & Mamon 2004) and the DF method of Wojtak et al. (2009). 3 In both comparisons, MAMPOSSt did much better on the velocity anisotropy, reaching double the accuracy on log v 2 r 1/2 /σ θ .
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 settle whether the orbits were circular, radial or isotropic, given his uncertainty on the mass profile. Łokas & Mamon (2003) considered both LOS dispersion and kurtosis profiles of the same 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 200 /3 given 2 One can instead bin the model, i.e. the radial profile of velocity variance (Diakogiannis et al. 2017) or mass density (Read & Steger 2017), but the choice of bins and regularization are complex issues. 3 MAMPOSSt was the 2nd most efficient among 23 algorithms in a recent cluster mass challenge (Old et al. 2015) to determine the mass normalizations of 1000 realistic mock clusters The Num richness-based method (Mamon et al., in prep.) that was the most efficient does not compute mass and anisotropy profiles. the published virial masses Lemze et al. 2008). 4 Benatov et al. (2006) applied the same anisotropy inversion to 5 clusters from the CAIRNS survey (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 2 clusters showing slightly tangential or isotropic velocities at r 200 , 1 mildly radial (β = 0.3) and 2 fully radial (β 0.95 at r 200 ). Analyzing a mere 64 dwarf galaxies in Coma, Adami et al. (2009) adopted the mass profile derived by Geller, Diaferio, & Kurtz (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 was composed of more luminous galaxies than that analyzed by Adami et al.).
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.  used anisotropy inversion (following the method of Solanes & Salvador-Solé 1990, hereafter SS90) on a joint analysis of 59 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, Biviano, & Mazure 2004 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 found that, for early type spirals (Sa, Sb), β rises to 0.7 at half the cluster 'virial' radius, r 200 , then falls to 0.35 at r 200 . In late-type spirals, they found that β = 0 (isotropic velocities) out to 0.6 r 200 and then rises to β = 0.3 at r 200 . They did not consider lenticular (S0) galaxies, nor has anybody else until now. The results of Biviano & Katgert were confirmed by Munari, Biviano, & Mamon (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-emissionline 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 4 The radius r ∆ is where the mean density of a system is ∆ times the critical density of the Universe at the cluster redshift. We will call r 200 the 'virial' radius (in quotes), but will also refer to the theoretical virial radius, which is close to r 100 1.35 r 200 with the analogous definition. We also define M ∆ = M(r ∆ ) and v ∆ = √ G M ∆ /r ∆ . 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 nonemission-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 indicate different conclusions. Hwang & Lee (2008) studied Abell 1795 (z = 0.06) with several mass profiles determined from X-ray analyses, and performed anisotropy inversion (using the technique of Bicknell et al. 1989) to deduce that early-and late-type galaxies had similar radial profiles of anisotropy starting radial in the core, dropping to very tangential β ≈ −3 at r 200 /2 and roughly flat beyond; but their study suffered from having only 160 galaxies. 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 d Σ(R) σ 2 los (R) /dR, where Σ is the surface density. 5 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 8 clusters out of 10, 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 200 ). 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 understand better 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? 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 0.1−2.4 keV X = 10 43.75 erg 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 6 and NED 7 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 extract simultaneously 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, the stacking method, while we explain in Sect. 4 the practical implementation of MAM-POSSt, 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 ΛCDM cosmological model with Ω m,0 = 0.3, Ω Λ,0 = 0.7, and H 0 = 70 km s −1 Mpc −1 .

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 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 6 http://www.sdss.org/ 7 http://ned.ipac.caltech.edu/ distribution of galaxies in PPS (projected radii and LOS veloci- 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, v 2 r , in the righthand-side of equation (5) is previously determined for a given η, by solving the spherical stationary JE of equation (1), which can be inverted by solving for d ln K β /d ln r = 2 β(r), yielding (van der Marel 1994; Mamon & Łokas 2005) where K β (r)/K β (s) = exp 2 s r β(t)dt/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 MAM-POSSt by replacing σ z in equation (4) by σ 2 z + 2 v . In our dataset, the typical velocity uncertainties are v ≈ 53 km s −1 , which turn out to be negligible relative to the cluster velocity dispersions, hence have little effect on the MAMPOSSt results.
MBB13 have tested that splitting the PPS into the 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 (S) 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, i.e. ln L = ln L E + ln L S0 + ln L 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.

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 per cent confidence), leaving 58 regular clusters. Irregular clusters violate the condition of smooth gravitational potential upon which mass-orbit modeling techniques are based. Foëx, : scale radii of the spatial distributions of elliptical (5), lenticular (6) and spiral (7) galaxies, determined from fits of the photometric data; (8-11): total number of galaxies in sample; (12-15): number of galaxies restricted to the projected radii analyzed by MAMPOSSt. Fig. 1. Projected phase space diagram of the sigv stacked cluster. Each symbol is a galaxy, with shapes and colors provided in the legend. The maximum projected radius corresponds to our maximum allowed value of 1.35 r 200 r 100 . The typical velocity errors are 53 km s −1 . The grey shaded region denotes the inner projected radii that are not considered in the MAMPOSSt analysis. The curves indicate the ±2.7 σ LOS (R) conditions obtained from the Clean algorithm (Sect. 3.1). (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, i.e. a luminosity satisfying log(L/L ) = 9.6. Assuming mass-to-light ratios of M/L V = 3 (ellipticals, from fig. 6 of Auger et al. 2010, who adopted a Chabrier 2003 initial mass function) and 2 (spirals), our median stellar mass is log(M stars /M ) = 10.0 for spirals and 10.4 for ellipticals.
Following Paper I, we assume that the clusters are centered on their Brightest Cluster Galaxy (BCG), defined within 0.5 r 200 . We run the Clean algorithm (MBB13) to remove interlopers in LOS velocity v LOS = c light (z − z)/(1 + z), where 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 gapper parameter (not a concentration) C = 4 (Girardi et al. 1996). Clean then iterates on the membership defined by the criterion |v − median(v)| < 2.7 σ NFW LOS (R), where the factor 2.7 was found to optimally recover the LOS velocity dispersion profile of pure NFW models (Mamon, Biviano, & Murante 2010). The term σ NFW LOS (R) requires the knowledge of the scale radius r −2 and the mass within, M(r −2 ), or equivalently the virial radius and the concentration c 200 = r 200 /r −2 . 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 ) 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, Flynn, & Gebhardt 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.) and the concentration is taken from the concentration-mass relation that Macciò, Dutton, & van den Bosch (2008) fit to the haloes of dissipationless cosmological N-body simulations.
We then restrict our cluster sample to the 54 clusters with at least N m = 30 members within R 200 (median velocity dispersion 763 km s −1 ).

Stacking
We stack 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 = R i, j r 200 /r 200, j and the LOS velocities v i, j to v i, j = v i, j v 200 /v 200, j (see footnote 4). For each cluster, we estimate r 200 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 200 , hence that of r 200 (Old et al. 2014(Old et al. , 2015. 3. An X-ray temperature based estimator derived from the mass-temperature relation of Arnaud, Pointecouteau, & Pratt (2005), which we call tempX. But this is limited to the 40 regular clusters with observed X-ray temperatures.
We then apply the Clean procedure one last time to each of these 3 stacks to remove yet undetected interlopers. We finally discard galaxies with projected radii beyond r 100 1.35 r 200 of our stacked clusters, where we adopted the Clean values of r 200 found in Paper I. Indeed, r 100 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 100 (Falco et al. 2013), so the limiting projected radius must satisfy R < 2 r 100 , hence our conservative choice of r 100 .
This leaves us with up to nearly 5000 galaxies for our 3 stacks. We also exclude 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 200 ), 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 200 (hence r 200 ) 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 (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. 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 3 morphological types of galaxies. A major part of the differences in these velocity dispersion profiles arises from the different scale radii of the 3 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 200 = 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 3 S0 galaxy LOS velocity dispersions are all above the isotropic prediction, while the 3 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.

Quick look at the stacked data
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 an anisotropic velocities of the spiral population. Alternatively, this may signify that the spiral population is not in dynamical equilibrium. We will discuss this possibility in Sect. 6.2.6. Line-of-sight velocity dispersion profiles for the elliptical (red), lenticular (green), and spiral (blue) galaxies in the sigv stacked cluster. Radial bins of 150 galaxies were used. Error bars are σ LOS / √ 2(N − 1) (expected for Gaussian LOS velocity distributions). The solid curves are predictions assuming the tracer number density profiles given in Paper I, isotropic velocities, for an NFW total mass profile of concentration c 200 = 4, while the dashed curves are 2nd-degree polynomial fits in log-log space.

Practical implementation of MAMPOSSt
We now present the practical implementation of the version of MAMPOSSt that we used here, 8 starting with our adopted radial profiles for number density, mass and anisotropy.

Tracer number density profiles
MAMPOSSt allows for the density profiles of the observed tracers to have up to 2 parameters (scale and shape -no normalization is required for tracers with negligible mass). In this work, we assume 1-parameter NFW number density profiles, following the fits to the photometric data performed in Paper I. We express the NFW number density profile as 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) where and Σ(1) = (1/3)/(2 ln 2 − 1). We adopt 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 σ.

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 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 3 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: 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 vir = M(cr/r vir )/ M(c), where c = r vir /r ρ is the concentration parameter. We consider the following dimensionless mass profiles (noting that M(1) = 1 by definition): NFW: the usual NFW model, whose density profile is given in equations (7) and (8), for which the inner and outer logarithmic slopes are respectively -1 and -3: cNFW: the cored NFW model (generalized NFW with zero inner slope), 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: where I x (a, b) is the regularized incomplete beta function, while 2 F 1 (a, b, c, x) is the ordinary (Gaussian) hypergeometric function. In the generalized NFW model, r ρ = (γ + 2) r s .
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 even better the density profiles of haloes than the NFW model:

Velocity anisotropy profile
MAMPOSSt also allows for a wide variety of velocity anisotropy profiles for each of the tracer components, based on up to 3 parameters (inner anisotropy, outer anisotropy and transition radius). We consider the following models for the anisotropy profile, β(r): T 0 : the same as the Tiret profile, but with β 0 = 0.
gOM: the generalized Osipkov-Merritt model (Osipkov 1979;Merritt 1985), where the usual Osipkov-Merritt anisotropy is recovered for β 0 = 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 β 0 = β ∞ . 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 ).

Assumptions
Our modeling assumes spherical symmetry for the visible and mass components, neglecting any possible rotation or other nonradial streaming motions. With our choice of regular clusters, we are in a better position to assume that the galaxies are noninteracting 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.

Maximum physical radius in integrals
We integrate the inversion of the JE (eq. [6]) out to 120 Mpc and the LOS integral of equation (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.

Free parameters
The following parameters can be set free in MAMPOSSt: logarithm of the mass normalization (log M 200 or equivalently log r 200 ); logarithm of the mass scale radius (or equivalently of the concentration c 200 = r 200 /r ρ ); inner slope of the mass density profile (from −1.99 to 0); logarithms of the tracer scale radii (for each of the 3 galaxy types; inner (r = 0) and outer (r → ∞) symmetrized velocity anisotropies (for each of the 3 galaxy types) 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 3 galaxy types (unless we assume 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 2 extra free parameters (but we then force 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 (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: We assumed a Gaussian prior on this relation with 0.1 dex uncertainty. We also performed extra MAMPOSSt runs with free cluster mass concentration, with different minimum and maximum allowed projected radii, with different minimum number of galaxies per cluster in building the stacks, and for the different stacks.

Marginal distributions (MCMC)
We determine the marginal distributions of the k free parameters using the Markov Chain Monte Carlo (MCMC) technique. We use 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 centered on the k-dimensional hyper-cube, with standard deviations σ k = [Max(θ k ) − Min(θ k )]/5. We then advance each chain, by moving from position θ old to θ new using a k-dimensional Gaussian proposal density function, with standard deviations equal σ k /2, i.e. 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 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 (obviously 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 run 6 chains in parallel on an 8-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 discard the first 2000 k steps, where the posterior distribution is dependent on the starting points of the chains (the burn-in phase).
We compute the radial profiles of mass density and of the velocity anisotropy of the 3 galaxy types from the marginal distributions of the free parameters (after discarding the burn-in phase).

Preamble
Most studies employ a single set of priors and show their results. Some will consider a few extra choices for their priors. Here, we have choices to make on the inner slope of the cluster mass profile, an additional mass profile for a possible BCG, fixed or free inner and anisotropy profiles for all three components. This led us to consider a large number of sets of priors. 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), i.e. increasing − ln L MLE (column 3). Our values of L MLE are really posteriors, but are close to likelihoods since all of 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 2 priors with the Einasto mass model, with either free index or fixed to n = 6, as found for ΛCDM haloes (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 given inner slope. Indeed, we found that the mass profile of the n=6 Einasto fits the NFW one to 8.5 per cent relative precision in the range 0.135 to 13.5 r −2 (0.03 to 3 r 200 ), equally spaced in log r, for typical cluster concentrations (the index n = 4.4 provides the best fit -4.8 per cent relative precision -to the NFW mass profile in this radial range). Similarly, the  (1997), which is a good approximation to the deprojection of a Sérsic (1968) model); (4): velocity anisotropy model ('iso' for isotropic, 'T' for generalized Tiret et al. and 'gOM' for generalized Osipkov-Merritt); (5-7): inner velocity anisotropy for E, S0 and S galaxies ('F' = free, 0 = isotropic); (8-10): outer velocity anisotropy (β) for E, S0 and S galaxies ('F' = free, 0 = isotropic); (11): velocity anisotropy radius tied to tracer scale radius? (12): MCMC convergence criterion (R −1 < 0.02 is considered as properly converged, worse convergence runs are shown in red italics); (13): minus log (maximum likelihood estimate, which really is a maximum posterior); (14): number of free parameters; (15): corrected Akaike Information Criterion; (16): Bayes Information Criterion. The best values for − ln L MLE , AIC and BIC are highlighted in bold, while blue italics and green italics respectively highlight the best NFW and gNFW models that do not have an extra BCG component.
mass profiles of Einasto models with free indices (up to n = 25) fit those of given gNFW models to better than 8.5 per cent rms relative precision in the same range of radii (6.1 per cent for γ ≥ −1.9).

Bayesian evidence methods
Using different priors can lead to different results, so one has to be careful in analyzing MAMPOSSt results. The runs leading to the highest likelihood L 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): and the Bayes Information Criterion (BIC, Schwarz 1978): Given our data sample with N data = 4682 (for the sigv stack), each extra parameter must decrease − ln L MLE by 4.23 to lead to a better (lower) BIC value, while with AIC it must decrease by only 1. According to Kass & Rafferty (1995), when a model has BIC lower than another model by 6 (10, 2) units, one can conclude that there is strong (very strong, positive) 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 − (ln N data − 2) ∆N pars . Therefore, the condition for strong AIC evidence for the simpler model compared to a another one with ∆N pars extra free parameters, given our data sample sizes (Table 1), is ∆AIC > 6 − 6.45 ∆N pars , i.e. ∆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 , i.e. ∆AIC < −12.45 compared to a simpler model with one less free parameter.

Preferred models
The models with the highest likelihoods from the MCMC analysis (lowest −L 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 (slightly) lower − ln L MLE than model 1. 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 7-parameter model 16 (NFW cluster with isotropic ellipticals and T 0 anisotropy with TAND for S0s and spirals) has the lowest AIC. The best BIC is reached for the 6-parameter model 23, with n=6 Einasto mass density, isotropic orbits for ellipticals and lenticulars, and T 0 anisotropy with TAND for spirals. 9 This model is also the 2th best model for AIC.

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 (cored NFW mass), and (marginally) model 16 (identical to model 24, but with free outer anisotropy for S0 galaxies).
In particular, there is strong evidence (∆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 equation (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, i.e. very strong) in moderately complex anisotropy models 11 vs. 13, but the evidence against gNFW is only weak for more complex anisotropy models, e.g. 1 vs. 9 and 4 vs. 10. And while the best AIC is reached for NFW cluster model 16, the 2nd 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 (higher) when assuming that the BCG mass follows the stellar mass, with an n=4 Prugniel & Simien (1997) 9 We found an even better BICs with the 5-parameter priors of NFW or n=6 Einasto mass and velocity anisotropy model β(r) = (1/2) r/(r + r β ) proposed by Mamon & Łokas (2005), which is a special case of the T model, with β → 1/2 at large radii. But this model has no physical basis, for example the anisotropy at r 200 increases with halo mass (Ascasibar & Gottlöber 2008;Lemze et al. 2012). 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.

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), 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 0 anisotropy for spirals), there is strong evidence against all variations on the velocity anisotropy, i.e. 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 0 anisotropy for S0s and spirals (both with TAND). However, the 2nd best non-Einasto model for AIC is model 4, where only the inner orbits of S0s and spirals are fixed to isotropic. The 3rd best model for AIC is model 13, which is intermediate in its complexity, with inner isotropy and free outer anisotropy for all 3 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.

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.

Best-fitting parameters
We now focus on just a few of the MAMPOSSt runs, by considering the highest likelihood model (model 1), 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). Fig. 3. MAMPOSSt marginal distributions (diagonal panels) and covariances for model 24 (the non-Einasto model with the strongest BIC evidence) in the sigv stack. The red stars and arrows show the parameters with the highest likelihoods. The priors are flat within the panels and zero outsize, except for the tracer NFW scale radii r E , r S0 , and r S , for which they are Gaussians with means in the middle and extending to ±3 σ on the panel edges. 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 v 2 r 1/2 /σ θ ) 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. 10 Figure C.1 in Appendix C is similar to Figure 3, but for model 19, which is the gNFW model with the 3rd lowest BIC model (same as model 24, except that it allows for free inner slope and free inner anisotropy for the spirals). This figure indicates that the inner density slope is weakly correlated with the mass normalization (left panel of 3rd row) and concentration of the mass profile (2nd panel of 3rd row), as well as with the outer spiral anisotropy (3rd panel of bottom row), and also correlated with the outer spiral anisotropy (3rd panel of bottom row). Figure C.2 in Appendix C shows the marginal distributions and covariances for model 11, with gNFW mass and T 0 anisotropy for all morphological types. One notices that the outer anisotropies of the 3 types are correlated. Notes: the parameters are uniformly distributed in the given ranges, except for the scale radii of the E, S0 and S galaxies, for which the mean and uncertainty are given and MAMPOSSt assumes Gaussian priors with dispersion σ equal to the uncertainty and cut at ±3 σ. The quoted values for the 3 models are the MLE estimates and (p 84 − p 16 )/2 estimates from the MCMC chains, where p i are the ith percentiles. Models 21 and 24 are isotropic for E and S0 galaxies. Table 3 shows, for models 1, 21, and 24, the MLE values and the uncertainties from the marginal distributions derived from the MAMPOSSt MCMC. For model 24 with an NFW mass model, the 'virial' radius is very well measured leading to a MLE value of r 200 = 1690 ± 20 kpc, i.e. with an uncertainty of only 0.005 dex. This value of r 200 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 200 = 1675 ± 23 kpc, still consistent with the Clean value, while the gNFW model 1 leads to r 200 = 1507 ± 59 kpc, 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 200 = 3.8 ± 0.4 Given the mass at r 200 = 1698 ± 24 kpc of 10 14.8 M , our concentrations are consistent with the values found for relaxed ΛCDM halos (c 200 = 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 200 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.). We will further discuss the concentration-mass 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 4 times more extended than that of the mass or of the ellipticals. We will return to this issue in Sect. 5.4.  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 3 morphological types, all with TAND). The MAMPOSSt model predictions reproduce very well the data.

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), i.e. 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 Fig. 5. Radial mass density profiles for models 21 (gNFW with isotropic orbits for the ellipticals and S0s, and T 0 velocity anisotropy for the spirals, top) and 3 (NFW cluster + NFW BCG, with T velocity anisotropy, bottom) in the sigv stack. In both models, the anisotropy radius is tied to the scale radius of the galaxy distribution. The shaded regions show the MAMPOSSt constraints for the cluster (light and dark grey) and the BCG (light and dark purple), where the light and dark zones respectively delimit 90% and 68% confidence intervals. The curves are the predictions from various analytical models, normalized to have the mass scale radii and the same density at the scale radius, simply to guide the eye. The scale of the bottom panel is different, and the curves to the left of the vertical line, denoting the minimum considered projected radius, are extrapolations. 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 L) only 1.6 higher than for model 24. Considering the cNFW model to be a physical one, it has the same number of 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 3 morphological types for models 24 and 21. We normalize the number density profiles by eliminating N(r ν ) between equation (7) and the average number of galaxies of given morphological type per cluster, N tot , between the minimum and maximum Fig. 6. Radial profiles of mass over number density ratios from MAM-POSSt for models 24 (NFW mass profile, top), and 21 (gNFW mass profile, bottom), for E (red), S0 (green), and S (blue) galaxies, for the sigv stack. The normalization is explained in equations (24)-(26). The horizontal lines are shown to highlight how well mass follows number. allowed projected radii, R min and R max , which we model as where N p (X) = 1 ln 2 − 1/2 for the NFW model, with N p (1) = (1 − ln 2)/(ln 2 − 1/2) and where C −1 is given in equation (10).
Recall that the NFW model was assumed for the number density profiles of the 3 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 almost perfectly the mass, 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 even better the mass. Since model 21 (bottom panel) 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 3 morphological types. Indeed, the mass over elliptical number density ratio is nearly consistent with being constant (horizontal line), 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 inner isotropy also for the spirals, shows that the spiral galaxies clearly have radial orbits in the outer regions of clusters. The other 3 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 3 models, but all 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.  Table 2 Table 4 illustrates the details of the anisotropy for 4 models. If we free the inner and outer anisotropies (but tie the anisotropy radii to the scale radii, model 1), we find that the anisotropy at r 200 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.

Outer versus inner velocity anisotropies
The inner anisotropies of the 3 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 quick look at the LOS velocity dispersion Fig. 7. Velocity anisotropy (eq. [19]) profiles of the E, S0 and S galaxies from MAMPOSSt for models 24 (NFW with T-TAND anisotropy for spirals only, upper left), 9 (NFW with T-TAND anisotropy, upper right), 1 (gNFW with T-TAND anisotropy, lower left), and 7 (gNFW with T anisotropy and free anisotropy radius, lower right), for the sigv stack. The hashed regions indicate the 68% confidence zone, while the curves display the 5th and 95th percentiles. The thick vertical grey line indicates the position of r 100 = 1.35 r 200 , which is close to the theoretical virial radius, where its width shows the uncertainty on log r 200 . 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 200 as a function of the value at 0.03 r 200 . We restrict 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 200 ) for spirals are above zero for all values of the inner anisotropy (at 0.03r 200 ), 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 3 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 200 (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 3 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).

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 3 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 (really posteriors), although their BIC Bayesian evidence is so high that they can be strongly rejected relative to the lowest BIC model.

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 2nd 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 Figs. C.1 and C.2), 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 analyze (Fig. 5). The BCG dominates the cluster within the inner 20 kpc, i.e 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.

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), we find that the constraints on the inner slope range from γ = −1.4 +0.5 −0.3 for sigv to −1.7 +0.3 −0.2 for tempX and −1.8 +0.4 −0.1 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: whereas 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 3 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 3 stacks, and close to proportional to the S0 number density profile, while the spirals trace poorly the mass profile, as they are more extended, hence ρ/ν 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 Ushaped mass-over-number profile. A close look indicates that the elliptical number density profile traces slightly better the mass density profile (within r 200 ) then do the S0 galaxies, except for the tempX stack where the two types trace the mass with similar accuracy. Again, the spiral galaxies trace poorly the mass profile.

Comparison with other work
In comparison, combining weak lensing at large radii, strong lensing at intermediate radii and stellar kinematics at low radii to study 7 regular clusters, Newman et al. (2013b) 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. 2013a). Their total mass profile is consistent with ours (NFW for lowest BIC model 24 as well as γ = −1.4 +0.5 −0.3 for model 21).  Figure 10 shows the constraints on the concentration-mass relation obtained with NFW model 24. Interestingly, with free concentration (flat prior 0 < log c 200 < 1), the contours match well the ΛCDM relation of Dutton & Macciò (2014) given in equation (20), especially for stacks tempX and sigv. It is not surprising that folding in this relation as a prior, we recover similar contours, simply closer to the relation itself. The MAM-POSSt analysis of the cluster kinematics matches well the concentrations obtained by weak lensing, except for one (fairly old) weak lensing study that strongly underestimates the concentration. Biviano et al. (2017) performed 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 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 200 = 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  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), hence 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 that the E stellar mass density profile follows 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 satllite radial distributions given the uncertain radiallydependent link between the minimum subhalo and galaxy stellar masses. One could also blame the NFW assumption for the E, S0 and S radial dirtibutions, but these are consistent with the data (Cava et al. 2017).

General trends
For the single component mass models, all anisotropy models that differ from that of T 0 anisotropy for the spirals and isotropy for the E and S0 galaxies are strongly rejected by BIC Bayesian evidence, with the sole exception of the case (model 16, the best one using AIC evidence) where S0 galaxies have outer anisotropy as do the spirals (Table 2). But the Bayesian evidence of model 24 (with outer anisotropy only for the spirals) against model 16 is ∆BIC = 5.5, i.e. "positive" and almost "strong" (∆BIC > 6). Hence, we have good confidence that the spiral population has radial orbits at r 200 , which are significantly more radial than at 0.03 r 200 . (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. Fig. 11. Same as Figure 8 comparing the priors on the cluster concentration for models 1 (TAND, top) and 7 (free anisotropy radius, bottom). Only 68% contours are shown.

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 (thin contours) to the ΛCDM concentration (eq. [20], thick contours). Since the concentrations obtained when they have wide (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 11 (for model 1) varying the minimum projected radius for galaxy selection: 25 (thin), 50 (medium, our standard case), and 100 kpc (thick). The thick dotted contour is for the case where clusters were stacked after centering on their X-ray positions instead of their BCGs, with R min = 50 kpc. Bottom: Same as top panel, for the maximum projected radius for galaxy selection: 0.67r 200 r 100 /2 (thin) and 1.35r 200 r 100 (thick, our standard case). Figure 12 indicate that our results are fully robust to our choice of minimum and maximum projected radii. Figure 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 (radial) inner anisotropy for the ellipticals (S0s). Fig. 13. Same as Figure 12 comparing the minimum number of members in individual clusters used for the stacks: 30 (thick) and 81 (thin). 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 200 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). 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 6 clusters measured by Hwang & Lee (2008) are much more radial than we (or others) found. Aguerri et al.  Figure 12 for 3 different stacks of clusters according to different estimates of the virial radii of individual clusters: our standard Clean method based on the distribution of galaxies in PPS (sigv, thick), and a richness in PPS method (Num, medium), and the masstemperature relation from X-ray measurements (tempX, thin).

Comparison with previous studies
(2017) and Benatov et al. (2006) found much lower inner radial anisotropy for Abell 85 and Abell 2199, respectively, than Hwang & Lee. Still, Benatov et al. and Wojtak & Łokas (2010) both find a very wide range of anisotropies at r 200 , sometimes perfectly radial. Biviano et al. (2013) found very radial anisotropy at r 200 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 (their anisotropy of the red galaxies at r 200 appears too tangential). There is tension between the very radial (β > 0.55) anisotropy of the passive galaxies at r 200 found by Annunziatella et al. (2016) for Abell 209 (for both low and high stellar mass). 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 us for early types 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 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 200 averaged over the 3 morphological types (Table 4 and Fig. 15).
Our results for spirals agree with the analysis of ENACS clusters by  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 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  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.

Infall
These constraints on inner and outer anisotropy help understand the mechanisms and timescales for the transformation of morphological types for the quenching of star formation (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), galaxy harassment from numerous minor flybys (Moore et al. 1996) or starving the galaxy of its supply of infalling gas either by tidal stripping (Larson, Tinsley, & Caldwell 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 Nbody 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). 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 (Figs. 8 and 14). As one moves to smaller physical radii, galaxies first entered the cluster at earlier times, and thus has 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 Fig. 15. Comparison with previous measurements (symbols) of the velocity anisotropy (eq. [19]) from MAMPOSSt (68% confidence contours) at 0.03 r 200 and r 200 for models 1 (gNFW, T, TAND, left) and 7 (same as 1 but with free anisotropy radii, right), all for the sigv stack. The top panels show the comparison with the full literature over the entire possible range of anisotropies, while the bottom panels show zooms, restricted to previous works that differentiated between galaxy types (with the same symbols meanings as in the top panel). The open and filled symbols respectively correspond to single clusters and stacks of clusters. The symbols are color-coded by the galaxy population: red for earlytype or passive, blue for late-type or star forming, and black for all galaxies. In the legend (top panels), the purple symbols denote studies that separately analyzed both passive (or red) and star forming (or blue) galaxies. The vertical and horizontal lines respectively indicate isotropic inner and outer velocities. The two symbols for Biviano & Katgert are for the Sa and Sb spirals (β(0.03 r 200 ) = 0.2) and for the later type spirals (β(0.03r 200 ) = −1.7). The two for Annunziatella et al. refer to low (β(0.03 r 200 ) ≈ −0.6) and high (β(0.03r 200 ) ≈ 0.1) mass galaxies. Finally, the two blue symbols for Aguerri et al. are for all (β(r 200 ) ≈ 0.5) and dwarf (β(r 200 ) ≈ 0) galaxies. spiral galaxies show statistical evidence of increasingly radial anisotropy profiles as one moves from 0.03 r 200 to r 200 (Table 4).

Isotropization
At small radii, the great majority of early-type galaxies is expected to have entered the cluster sufficiently long ago to have been morphologically transformed from their spiral progenitors (again, as in the onion model of cluster growth of Gott 1975).
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/(8 ln N) 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), 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 (Fig. D.1) that at the very least 30 orbits are required to isotropize. According to fig. B1 of Tollet et al., who considered a growing NFW cluster, 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), hence 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 of order [M(r)/m]/ ln(1+ M(r)/m) times the orbital time (Mamon 1995;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, i.e. the transfer of orbital energy into internal energy. However, simulations indicate that galaxies bounce out of clusters to 1 to 2.5 virial radii 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 (figure 3 of Fakhouri, Ma, & Boylan-Kolchin 2010), cluster-mass halos typically undergo 0.8 major mergers since z = 1 (7 Gyr for the cosmology of the simulation studied). At z = 1, the Hubble constant is 1.75 times greater, hence the orbital time is 1.75 times shorter, i.e. a little over 2 Gyr. Thus, since z = 1, roughly one-third of clusters undergo a major merger, hence one-third of galaxies in our stacked cluster would have gone through rapid isotropization. However, this fraction is an overestimate 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. The inner isotropy of spiral galaxies: a selection effect due to single orbits?
It is more difficult 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, Mamon, & Raychaudhury (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 find (their fig. 8) 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 3 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 of 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 2nd 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 2nd 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. 11 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 full tangential (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. r βsym 0 -2 2 rmin circular radial isotropic Fig. 16. Illustration of a wide (salmon) and narrow (blue) range of pericenters, leading to different mixes of velocity anisotropies at r = r min (light green): in the case of a narrow range of pericenters just below r min , the orbits at r min cannot be radial, leading to less radial velocity anisotropy at r min . At pericenter and apocenter, orbits are necessarily circular, hence β sym = −2.
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. The 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 (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 2nd release of the Gaia astrometric mission, the 3D orbits of the dwarf spheroidals of the Milky Way have been reconstructed. The left panel of figure D.3 of Gaia Collaboration et al. (2018) indicates that the apocenters of all the 9 dwarfs, except Leo II, are less than 3 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 et al. 2018), which presumably is also caused by tidal selection.
6.2.8. The orbits of S0 galaxies: a clue to their formation?
Finally, one may wonder 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 Figures 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 Figure 14 displaying model 1, in the Num stack S0s display isotropic outer orbits, the ellipticals prefer slightly tangential outer orbits (but 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 where f i is the fraction of MCMC points (after burn-in) for type i that lie in a given cell of β sym (r 200 ) vs. β sym (0.3 r 200 ) and equation (27) is estimated over the entire set of cells. The Pearson  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 3 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), i.e. 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 should 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.

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 on 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, which is of 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 (S0s), the scale radii of the 3 morphological types are much more similar in irregular clusters. Perhaps the violent relaxation in irregular clusters perturbs the orbits more efficiently than the different histories of the 3 morphological types differentiate the orbits.

Conclusions
We ran the MAMPOSSt mass/orbit modeling algorithm on 3 stacks of 54 regular, nearby (z ≈ 0.05) clusters from the WINGS sample, composed of up to 4682 galaxies located between 0.03 r 200 and r 200 , split between ellipticals, lenticulars, and spirals (including irregulars). MAMPOSSt is a Bayesian method that jointly fits the distribution of galaxies of these 3 morphological types in projected phase space, fitting for the shape of the total mass profile on one hand and of the 3 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 200 (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 (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 very well the total mass density profile, while that of S0s only marginally does so and the spiral one 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 200 . -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 200 than for the spirals. -The velocity anisotropy of the ellipticals is consistent with isotropic anywhere, even though the highest likelihoods are reached for slightly tangential inner orbits and mildly radial anisotropy (β = 0.19 ± 0.25) at r 200 . -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, which 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 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 3 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 200 , E and S0 galaxies should be a mix of the virialized (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.  ⊥ , whose inverse is given in equation (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, highlighting the need for orbit averaging, as done here.
More massive galaxies do not relax faster with equation (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 1 2 L 2 r 2 + Φ(r) − E = 0 . (E.1) We built an isotropic Hernquist (1990) model following the method of Kazantzidis, Magorrian, & Moore (2004), where we first draw random radii, compute the gravitational potential and then we draw velocities from f (v|r) = v 2 f (v 2 /2 + Φ(r)), where f ≡ f (E) is the 6D distribution function of the isotropic Hernquist model, given by Hernquist (1990). We tested 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 equation (E.1), where the two roots correspond to the pericenter and apocenter.
The extraction of the roots of equation (E.1) for each of the 10 5 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 equation (E.1) with r = r i . We first set r peri = 0 and r apo = 10 6 a (where a is the scale radius of the Hernquist model) for all the particles. Noting that the LHS of equation (E.1) must be less than or equal to 0, since it represents − 1 2 (dr/dt) 2 , we then vectorially adjusted r peri and r apo with the conditions if |LHS| < |oldLHS| & LHS < oldLHS update r peri if |LHS| < |oldLHS| & LHS > oldLHS update r apo save LHS to oldLHS It took 15 (1.5) seconds to process 100 000 particles with 0.001 (0.01) dex precision in this manner with a script language (SM, aka SuperMongo) on a single processor. We found a median r apo /r peri of 3.79 with an uncertainty of 0.02 (from 10 trials).