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

G.A. Mamon, A. Cava, A. Biviano, A. Moretti, B. Poggianti, D. Bettoni

► To cite this version:

G.A. Mamon, A. Cava, A. Biviano, A. Moretti, B. Poggianti, et al. Structural and dynamical modeling of WINGS clusters. II. The orbital anisotropies of elliptical, spiral and lenticular galaxies. Astronomy and Astrophysics - A&A, 2019, 631, pp.A131. 10.1051/0004-6361/201935081. hal-02008897

HAL Id: hal-02008897

https://hal.science/hal-02008897

Submitted on 1 Jun 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Structural and dynamical modeling of WINGS clusters
## II. The orbital anisotropies of elliptical, spiral, and lenticular galaxies

G. A. Mamon<fnref n="1" />, A. Cava<fnref n="2" />, A. Biviano<fnref n="3" />,<fnref n="1" />,<fnref n="4" />, A. Moretti<fnref n="5" />, B. Poggianti<fnref n="5" />, and D. Bettoni<fnref n="5" />

Received 17 January 2019 / Accepted 2 August 2019

## ABSTRACT

The orbital shapes of galaxies of different classes are a probe of their formation and evolution. The Bayesian MAMPOSSt mass-orbit modeling algorithm is used to jointly fit the distribution of elliptical, spiral-irregular, and lenticular galaxies in projected phase space, on three pseudo-clusters (built by stacking the clusters after re-normalizing their positions and velocities) of 54 regular clusters from the Wide-field Nearby Galaxy-clusters Survey (WINGS), with at least 30 member velocities. Our pseudo-clusters (i.e., stacks) contain nearly 5000 galaxies with available velocities and morphological types. Thirty runs of MAMPOSSt with different priors are presented. The highest MAMPOSSt likelihoods are obtained for generalized Navarro-Frenk-White (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 best trace the mass distribution while S0s are close. Spiral galaxies show increasingly elongated orbits at increasing radii, as do S0s on two stacks, and ellipticals on one stack. The inner orbits of all three types in the three stacks are consistent with isotropy. Spiral galaxies should transform rapidly into early-types given their much larger extent in clusters. Elongated outer orbits are expected for the spirals, a consequence of their recent radial infall into the cluster. The less elongated 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.

**Key words.** galaxies: kinematics and dynamics – dark matter – galaxies: clusters: general

## 1. Introduction

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

The measure of the elongations of galaxy orbits is therefore a fundamental tool to understand the formation of clusters, and how galaxy orbits vary with cluster mass, elongation, and large-scale environment. One may go further and understand

how orbital elongations depend on galaxy stellar mass or luminosity, as well as galaxy specific star formation rate or color or even morphological type, as in the present study.

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

$$ \frac{d(v\langle v_r^2 \rangle)}{dr} + 2\frac{\beta}{r}v\langle v_r^2 \rangle = -v\frac{GM}{r^2}, \quad (1) $$

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

$$ \beta = 1 - \frac{\langle v_{\theta}^2 \rangle + \langle v_{\phi}^2 \rangle}{2 \langle v_r^2 \rangle}, \quad (2) $$

<footnote n="1">¹ Institut d'Astrophysique de Paris (UMR 7095: CNRS & Sorbonne Université), 98 bis Bd Arago, 75014 Paris, France
e-mail: gam@iap.fr</footnote>

<footnote n="2">² Department of Astronomy, University of Geneva, 51 Ch. des Maillettes, 1290 Versoix, Switzerland</footnote>

<footnote n="3">³ INAF-Osservatorio Astronomico di Trieste, Via Tiepolo 11, 34143 Trieste, Italy</footnote>

<footnote n="4">⁴ IFPU – Institute for Fundamental Physics of the Universe, via Beirut 2, 34014 Trieste, Italy</footnote>

<footnote n="5">⁵ INAF-Osservatorio Astronomico di Padova, Vicolo Osservatorio 5, 35122 Padova, Italy</footnote>
represents the (spherical, scalar) anisotropy of the mean squared
velocities, called the velocity anisotropy parameter (or simply
anisotropy), and is expressed in terms of the root-mean-squared
(rms) velocities (β usually varies with radius). Radial, isotropic,
and circular orbits respectively have $\beta = 1, \beta = 0$, and $\beta \to -\infty$.
For equilibrium systems, there are no net meridional streaming
motions, hence $\langle v_{\theta}^2 \rangle = \sigma_{\theta}^2$. For galaxy clusters, one neglects
rotation leading to $\langle v_{\phi}^2 \rangle = \sigma_{\phi}^2$. By symmetry, $\sigma_{\phi} = \sigma_{\theta}$. On the
other hand, the JE contains radial streaming motions, for example
from infall.

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

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

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

A disadvantage of Jeans analysis methods is that they usually
require the radial binning of the data<fnref n="2" />. One way around this is to
specify the form of the LOS velocity distribution function. Start-
ing with Walker et al. (2009), it has become popular to assume a
Gaussian LOS velocity distribution function (in studies of dwarf
spheroidal galaxies, but this has not yet been done for clusters).
But since the LOS velocity distribution function depends on the
velocity anisotropy (Merritt 1987), it is not desirable to mea-
sure 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 $\beta$ have $f \propto L^{-2\beta}$. One can then compute not only
moments in radial bins, but also at specific positions of the trac-
ers in PPS. Early studies used fairly arbitrary choices for the
DF. Wojtak et al. (2008) found that a separable form for $f(E, L)$
provided a good match to the halos in cosmological N-body sim-
ulations. This method, adapted to observational data by Wojtak
et al. (2009), is powerful, but slow as it involves computing triple
integrals (Dejonghe & Merritt 1992). Nevertheless, it has been
successfully applied to clusters (Wojtak & Łokas 2010; Wojtak
et al. 2011) and galaxies (Wojtak & Mamon 2013). A promising
method is to express the DF in terms of action-angle variables
(Vasiliev 2019).

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

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

Another way to lift the MAD is to adopt the mass profile from other methods. Applying anisotropy inversion (from Bicknell et al. 1989) to the mass profile of Abell 1689 ($z = 0.18$) derived from weak lensing, Natarajan & Kneib (1996) found that the velocities are isotropic in the core, and become radial ($\beta =$

<footnote n="1">¹ Conversely, one can assume an anisotropy profile and deduce a non-parametric mass profile (Mamon & Boué 2010; Wolf et al. 2010), although this mass inversion is less relevant for the study of orbits of the present work.</footnote>

<footnote n="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.</footnote>

<footnote n="3">³ MAMPOSSt was the second 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.</footnote>
0.5) at 800 kpc, which corresponds to roughly $r_{200}/3$ given the published virial masses (Sereno et al. 2013; Lemze et al. 2008)$^4$. Benatov et al. (2006) applied the same anisotropy inversion to five clusters from the Cluster And Infall Region Nearby Survey (CAIRNS) ($z = 0.03$ to $0.3$) whose mass profiles were obtained from X-rays and/or lensing. They derived anisotropy profiles that were radial at the very center, and showed a diversity of profiles in their bodies, with two clusters showing slightly tangential or isotropic velocities at $r_{200}$, one mildly radial ($\beta = 0.3$) and two fully radial ($\beta \approx 0.95$ at $r_{200}$). Analyzing a mere 64 dwarf galaxies in Coma, Adami et al. (2009) adopted the mass profile derived by Geller et al. (1999) with the caustic method (Diaferrio 1999), then fit a constant $\beta$ to the LOS dispersion profile, to obtain $\beta = 0.4 \pm 0.2$ or $0.7 \pm 0.1$ depending on their fit of the number density profile. So for the Coma cluster, the orbital anisotropy may be a function of galaxy mass (since the sample analyzed by Łokas & Mamon 2003 was composed of more luminous galaxies than that analyzed by Adami et al. 2009).

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

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

Other studies of Biviano and collaborators point to somewhat radial outer orbits for passive galaxies. Analyzing a $z = 0.4$ cluster with 600 member velocities, Biviano et al. (2013) concluded (using MAMPOSSt and performing SS90 anisotropy inversion

from the MAMPOSSt mass profile) that both star forming and passive galaxies have isotropic orbits inside and radial orbits outside. Working on stacked $z \sim 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 ($\beta = 0.4$, with large error bars). Annunziatella et al. (2016) performed SS90 anisotropy inversion on a parametric Navarro-Frenk-White (NFW, Navarro et al. 1996) model fit to lensing data for Abell 209 ($z = 0.21$), and found that passive galaxies display radial outer anisotropy, while their inner anisotropy depends on their stellar mass (slightly radial at high mass and tangential at low mass). Capasso et al. (2019) analyzed non-emission-line galaxies in 110 SZ-selected clusters at $0.26 < z < 1.32$, using both MAMPOSSt and anisotropy inversion. They concluded that passive galaxies have isotropic inner orbits and more radial outer orbits (their outer anisotropy varies with increasing redshift in an oscillatory manner).

Other anisotropy inversions of clusters led to different conclusions. Hwang & Lee (2008) studied ten nearby ($z < 0.12$) clusters, determining the mass profiles determined from X-ray analyses, and performed anisotropy inversion (using the technique of Bicknell et al. 1989). They deduced that galaxies typically have radial orbits in the core, sometimes (e.g., for Abell 1795) dropping to very tangential $\beta \approx -3$ at $r_{200}/2$. Aguerri et al. (2017) analyzed Abell 85 ($z = 0.06$) using a parametric NFW model for the mass profiles obtained by the causatic 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 ($\beta = 0.0 \pm 0.3$) for blue dwarf galaxies, but very radial outer orbits for red galaxies ($\beta = 0.7 \pm 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 \to \infty$ and inwards to $r = 0$. Moreover, all anisotropy inversion algorithms involve differentiating the observational data – they require the knowledge of $d[\Sigma(R)\sigma_{\text{los}}^2(R)]/dR$, where $\Sigma$ 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 eight out of ten clusters, they found that the distribution of galaxies in PPS implied isotropic inner orbits and radial outer orbits ($\beta = 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

<footnote>$^4$ The radius $r_\Delta$ is where the mean density of a system is $\Delta$ times the critical density of the Universe at the cluster redshift. We will call $r_{200}$ the virial radius, but will also refer to the theoretical virial radius, which is close to $r_{100} \approx 1.35 r_{200}$ with the analogous definition. We also define $M_\Delta = M(r_\Delta)$ and $v_\Delta = \sqrt{G M_\Delta / r_\Delta}$.</footnote>

<footnote>$^5$ The algorithms of Bicknell et al. (1989) and Dejonghe & Merritt (1992) also involve differentiating the dynamical pressure, which is the sum of integrals of the model and of the data (see Appendix B, which shows that the two algorithms are equivalent for the specific case of virial equilibrium).</footnote>
type, distinguishing between elliptical, spiral, and (for the first time) lenticular galaxies. We thus seek to settle the debate on the different orbital shapes of elliptical and spiral galaxies in clusters, but also wish to better understand S0 galaxies through their orbital shapes. Should the S0 galaxies resemble spiral galaxies in their orbital shapes, because they both have disks, or should they resemble more elliptical galaxies, because they have large bulges and old stellar populations? Answering this may help understand whether S0s originate from spiral galaxies that saw their disks fade, grew their bulges by mergers, or possibly even originate from ellipticals that accreted disks.

We use the W1de field Nearby Galaxy clusters Survey (WINGS), which contains X-ray-selected (median luminosity $L_X^{0.1-2.4 \text{ keV}} = 10^{43.75} \text{ erg s}^{-1}$) clusters at redshifts $0.04 < z < 0.07$ and that are located on the sky at least $20^{\circ}$ from the Galactic Plane (Fasano et al. 2006). The WINGS spectroscopic dataset (Cava et al. 2009; Moretti et al. 2014) has been complemented with redshifts from literature data collected through the SDSS-DR7⁶ and NED<fnref n="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 simultaneously extract the mass profile and the anisotropy profiles of the different classes of galaxies. This will be the first application of MAMPOSSt to a large sample of stacked nearby clusters. Stacking clusters reduces the intrinsic effects of triaxiality that create unavoidable biases in the derived mass profile parameters of individual halos (MBB13).

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

## 2. The MAMPOSSt algorithm

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

In Paper I, we determined the number density profiles of the three morphological types by fits of NFW plus constant back-

ground models on the photometric data. Given the known distribution of projected radii $R$, the MAMPOSSt likelihood of the distribution of galaxies in PPS (projected radii and LOS velocities $v_z$) is (MBB13)

$$ - \ln \mathcal{L} = - \sum \ln p(v_z|R, r_v, \eta), \quad (3) $$

$$ p(v_z|R, \eta) = \frac{\sqrt{2/\pi}}{\Sigma(R)} \int_R^\infty \frac{v}{\sigma_z} \exp\left[-\frac{v_z^2}{2\sigma_z^2}\right] \frac{r}{\sqrt{r^2 - R^2}} dr, \quad (4) $$

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

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

$$ v(r) \langle v_r^2(r) \rangle = \frac{1}{K_\beta(r)} \int_r^\infty K_\beta(s) v(s) \frac{G M(s)}{s^2} ds, \quad (6) $$

where $K_\beta(r)/K_\beta(s) = \exp[2 \int_r^s \beta(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 $\epsilon_v$ are accounted for by MAMPOSSt by replacing $\sigma_z$ in Eq. (4) by $\sqrt{\sigma_z^2 + \epsilon_v^2}$. In our dataset, the typical velocity uncertainties are $\epsilon_v \approx 53 \text{ km s}^{-1}$, which turn out to be negligible relative to the cluster velocity dispersions. Hence these typical velocity uncertainties have little effect on the MAMPOSSt results. MBB13 tested that splitting the PPS into separate determinations of the number density profile from the distribution of projected radii on one hand, and the mass and anisotropy profiles from the distribution of LOS velocities at given projected radius on the other hand, leads to virtually the same parameters as the standard joint fit of PPS.

Finally, MAMPOSSt can jointly analyze the positions in PPS of several independent tracers, such as elliptical (E), lenticular (S0), and spiral-irregular (S) galaxies. Since the three populations of E, S0, and S galaxies move in the same gravitational potential but have different spatial and velocity distributions (see Fig. 1), making a joint analysis of the three populations by allowing a different $\beta(r)$ for each of them results in a more stringent constraint on the remaining parameters for $M(r)$. The joint likelihood is the product of those from each of the tracers: $\ln \mathcal{L} = \ln \mathcal{L}_E + \ln \mathcal{L}_{S0} + \ln \mathcal{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.

## 3. Data preparation

### 3.1. Sample and interloper removal

In Paper I, we applied the substructure test of Dressler & Shectman (1988) to our initial sample of 73 WINGS clusters containing over 10 000 galaxies. This test led to 15 irregular

<footnote n="6">⁶ http://www.sdss.org/</footnote>

<footnote n="7">⁷ http://ned.ipac.caltech.edu/</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Projected phase space diagram of sigv-stacked pseudo-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} \approx r_{100}$. The typical velocity errors are 53 km s$^{-1}$. The gray shaded region denotes the inner projected radii that are not considered in the MAMPOSSt analysis. The curves indicate the $\pm 2.7 \sigma_{LOS}(R)$ conditions obtained from the Clean algorithm (Sect. 3.1).</figcaption></figure>

clusters (at 95 percent confidence), leaving 58 regular clusters. Irregular clusters violate the condition of smooth gravitational potential upon which mass-orbit modeling techniques are based. Foëx et al. (2017) found that discarding clusters identified as irregular with the Dressler & Shectman (1988) statistic leads to a much better match between MAMPOSSt and X-ray based masses. We thus discard the irregular clusters and only consider the 58 regular clusters.

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

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

Clean begins by searching for significant gaps in the distribution of LOS velocities using the gapper technique of Wainer & Thissen (1976) with a gapper parameter (not a concentration) $C = 4$ (Girardi et al. 1996). Clean then iterates the membership defined by the criterion $|v - \text{median}(v)| < 2.7 \sigma_{\text{NFW}}(R)$, where the factor 2.7 was found to optimally recover the LOS velocity dispersion profile of pure NFW models (Mamon et al. 2010). The term $\sigma_{\text{LOS}}^{\text{NFW}}(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 (Mamon et al. 2010) with a transition radius equal to the NFW scale radius. On first pass, the aperture velocity dispersion is measured by the robust median absolute deviation technique (e.g., Beers et al. 1990) and the concentration of the NFW model is taken as $c = 4$, typical of rich clusters. On subsequent passes, the aperture velocity dispersion is measured using the biweight estimator (e.g., Beers et al. 1990) and the concentration is taken from the concentration-mass relation that Macciò et al. (2008) fit to the halos of dissipationless cosmological N-body simulations. We then restricted our cluster sample to the 54 clusters with at least $N_m = 30$ members within $R_{200}$ (median velocity dispersion 763 km s$^{-1}$).

## 3.2. Stacking

We stacked the clusters into a pseudo-cluster by rescaling the projected radii $R_{i,j}$ of galaxies in cluster $j$ to the mass-weighted average cluster with $R'_{i,j} = R_{i,j} \langle r_{200} \rangle / r_{200,j}$ and the LOS velocities $v_{i,j}$ to $v'_{i,j} = v_{i,j} \langle v_{200} \rangle / v_{200,j}$ (see footnote 4). For each cluster, we estimated $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, 2015).

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

We then applied the Clean procedure one last time to each of these three stacks to remove remaining undetected interlopers. We finally discarded galaxies with projected radii beyond $r_{100} \approx 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 three stacks. We also excluded the very inner region since it is dominated by the internal dynamics of the BCG, rather than by the overall cluster (see, e.g., Biviano & Salucci 2006). Using a minimum projected radius of 50 kpc (roughly $0.03 r_{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_{\text{true}}$, thus leading to positive (negative) bias for clusters of low (respectively, high) mass. It thus seems preferable to avoid this non-unity recovered vs. true log mass slope of Num, which may bias our results. We thus adopt sigv as our main stacking method, but will compare in Sect. 6 the sigv results with those from the Num and tempX stacking methods.
<table><caption>**Table 1. Stacked clusters.**</caption><thead><tr><th rowspan="2">Method<br>(1)</th><th rowspan="2">N<sub>min</sub><br>(2)</th><th rowspan="2">r<sub>200</sub><br>(kpc)<br>(3)</th><th rowspan="2">log M<sub>200</sub><br>(M<sub>☉</sub>)<br>(4)</th><th rowspan="2">log r<sub>E</sub><br>(kpc)<br>(5)</th><th rowspan="2">log r<sub>S0</sub><br>(kpc)<br>(6)</th><th rowspan="2">log r<sub>S</sub><br>(kpc)<br>(7)</th><th colspan="3">Number (all)</th><th colspan="3">Number (R<sub>min</sub> ↔ R<sub>max</sub>)</th></tr><tr><th>Total<br>(8)</th><th>E<br>(9)</th><th>S0<br>(10)</th><th>S<br>(11)</th><th>Total<br>(12)</th><th>E<br>(13)</th><th>S0<br>(14)</th><th>S<br>(15)</th></tr></thead><tbody><tr><td>sigv</td><td>30</td><td>1749</td><td>14.79</td><td>2.690 ± 0.097</td><td>2.813 ± 0.107</td><td>3.314 ± 0.105</td><td>4816</td><td>1701</td><td>1871</td><td>1244</td><td>4682</td><td>1607</td><td>1845</td><td>1230</td></tr><tr><td>sigv</td><td>81</td><td>1846</td><td>14.86</td><td>2.708 ± 0.128</td><td>2.813 ± 0.087</td><td>3.322 ± 0.122</td><td>3361</td><td>1164</td><td>1334</td><td>863</td><td>3289</td><td>1109</td><td>1321</td><td>859</td></tr><tr><td>Num</td><td>30</td><td>1629</td><td>14.69</td><td>2.690 ± 0.106</td><td>2.799 ± 0.103</td><td>3.270 ± 0.103</td><td>4730</td><td>1676</td><td>1838</td><td>1216</td><td>4593</td><td>1583</td><td>1807</td><td>1203</td></tr><tr><td>tempX</td><td>30</td><td>1710</td><td>14.76</td><td>2.602 ± 0.152</td><td>2.699 ± 0.078</td><td>3.369 ± 0.126</td><td>3607</td><td>1264</td><td>1441</td><td>902</td><td>3523</td><td>1198</td><td>1426</td><td>899</td></tr></tbody></table>

**Notes.** Columns are as follows. (1): stacking method; (2): minimum number of galaxies per cluster; (3): $r_{200}$ obtained from Clean; (4): log $M_{200}$ corresponding to $r_{200}$; (5–7): 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.

### 3.3. Brief overview of stacked data

Figure 1 displays the projected phase-space diagram, highlighting the distributions of galaxies of different morphological types. Figure 2 shows the LOS velocity dispersion profiles for the three morphological types of galaxies. A major part of the differences in these velocity dispersion profiles arises from the different scale radii of the three number density profiles. The solid curves in Fig. 2 show the predicted LOS velocity dispersion profiles for an NFW model, with virial radius and tracer radii adopted from Table 2 of Paper I, also assuming isotropic orbits and a mass concentration of $c_{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 three S0 galaxy LOS velocity dispersions are all above the isotropic prediction, while the three elliptical LOS velocity dispersions are all below the isotropic prediction. This suggests that the inner anisotropies of S0s may be somewhat radial while those of ellipticals may be somewhat tangential.

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

## 4. Practical implementation of MAMPOSSt

We now present the practical implementation of the version of MAMPOSSt that we used here<fnref n="8" />. We present in turn our adopted radial profiles for number density, mass, and velocity anisotropy.

### 4.1. Tracer number density profiles

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

$$v(r) = \frac{N(r_v)}{4\pi r_v^3} \left(\frac{r}{r_v}\right) \quad (7)$$

$$\bar{v}(x) = \frac{x^{-1} (x + 1)^{-2}}{\ln 2 - 1/2}, \quad (8)$$

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Line-of-sight velocity dispersion profiles for elliptical (red), lenticular (green), and spiral (blue) galaxies for sigv-stacked pseudo-cluster. Radial bins of ≈150 galaxies were used. Error bars are $\sigma_{LOS}/\sqrt{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 second-degree polynomial fits in log-log space.</figcaption></figure>

where $r_v$ 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)

$$\Sigma(R) = \frac{N_p(r_v)}{\pi r_v^2} \bar{\Sigma}(X) = \frac{1/(x^2 - 1) - C^{-1}(1/X)/|X^2 - 1|^{3/2}}{2 \ln 2 - 1}, \quad (9)$$

where

$$C^{-1}(x) = \begin{cases} \cos^{-1} x & (x < 1), \\ \cosh^{-1} x & (x > 1), \end{cases} \quad (10)$$

and $\bar{\Sigma}(1) = (1/3)/(2 \ln 2 - 1)$.

<footnote n="8">⁸ https://gitlab.com/gmamon/MAMPOSSt</footnote>
We adopted the values of the decimal logarithm of the NFW scale radii of the stacked clusters, obtained from fits of a (projected) NFW cluster model plus constant background to the photometric data of the stacked clusters. These values (and their uncertainties) are provided in Table 1 for all three stacks, for galaxy morphological types E (ellipticals), S0 (lenticulars), and S (spirals). We adapted MAMPOSSSt 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\sigma$.

## 4.2. Mass profile

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

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

$$M(r) = M(r_{\rho}) \tilde{M} \left( \frac{r}{r_{\rho}} \right), \qquad (11)$$

where $r_{\rho}$ 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} = \tilde{M}(cr/r_{vir})/\tilde{M}(c)$, where $c = r_{vir}/r_{\rho}$ is the concentration parameter. We consider the following dimensionless mass profiles (noting that $\tilde{M}(1) = 1$ by definition):

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

$$\tilde{M}(x) = \frac{\ln(x+1) - x/(x+1)}{\ln 2 - 1/2}. \qquad (12)$$

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

$$\tilde{M}(x) = \frac{\ln(2x+1)(3x+1)/(2x+1)^2}{\ln 3 - 8/9}. \qquad (13)$$

In the cored-NFW model, $r_{\rho}$ is equal to twice the scale radius.

**gNFW.** The generalized NFW model, whose density profile is $\rho(r) \propto r^{\gamma} (r+r_{s})^{-3-\gamma}$, for which the inner and outer logarithmic slopes are respectively $\gamma$ (e.g., -1 for NFW) and -3 again:

$$\begin{aligned} \tilde{M}(x) &= \frac{I_{-(\gamma+2)x} [\gamma + 3, -(\gamma + 2)]}{I_{-(\gamma+2)} [\gamma + 3, -(\gamma + 2)]}, \\ &= x^{\gamma+3} \frac{2 F_1 [\gamma + 3, \gamma + 3, \gamma + 4, -(\gamma + 2)x]}{2 F_1 [\gamma + 3, \gamma + 3, \gamma + 4, -(\gamma + 2)]}, \end{aligned} \qquad (14)$$

where $I_x(a, b)$ is the regularized incomplete beta function, while $2F_1(a, b, c, x)$ is the ordinary (Gaussian) hypergeometric function. In the generalized NFW model, $r_{\rho} = (\gamma + 2) r_s$.

**Hernquist.** The Hernquist (1990) model with density profile $\rho(r) \propto r^{-1} (r+r_s)^{-3}$, for which the inner and outer logarithmic slopes are respectively -1 and -4:

$$\tilde{M}(x) = 9 \left( \frac{x}{x+2} \right)^2 . \qquad (15)$$

In the Hernquist model, $r_{\rho} = r_s/2$.

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

$$\tilde{M}(x) = \frac{\gamma(3n, 2n x^{1/n})}{\gamma(3n, 2n)}. \qquad (16)$$

## 4.3. Velocity anisotropy profile

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

T. Tiret et al. (2007) profile,

$$\beta_T(r) = \beta_0 + (\beta_\infty - \beta_0) \frac{r}{r + r_\beta}, \qquad (17)$$

$T_0$. The same as the Tiret profile, but with $\beta_0 = 0$.

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

$$\beta_{gOM}(r) = \beta_0 + (\beta_{\infty} - \beta_0) \frac{r^2}{r^2 + r_{\beta}^2}, \qquad (18)$$

where the usual Osipkov-Merritt anisotropy is recovered for $\beta_0 = 0$ and $\beta_{\infty} = 1$. Note that the usual constant anisotropy model can be retrieved as a singular case of the above T and gOM models, when assuming $\beta_0 = \beta_{\infty}$.

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

## 4.4. Assumptions

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

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

## 4.5. Maximum physical radius in integrals

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

The following parameters can be set free in MAMPOSSSt:

- logarithm of the mass normalization (log $M_{200}$ or equivalently log $r_{200}$);

- logarithm of the mass scale radius (log $r_ρ$ 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 three galaxy types;

- inner ($r = 0$) and outer ($r \to \infty$) symmetrized velocity anisotropies (for each of the three galaxy types)

$$ \beta_{\text{sym}} = 2 \frac{\langle v_r^2 \rangle - \sigma_{\theta}^2}{\langle v_r^2 \rangle + \sigma_{\theta}^2} = \frac{\beta}{1 - \beta/2}, \quad (19) $$

which can be as low as -2 for circular orbits and as high as +2 for radial orbits, and where $\beta_{\text{sym}} \to \beta$ for $|\beta| \ll 1$ (we allow -1.8 to 1.8);

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

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

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

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

$$ \log c_{200} = 2.13 - 0.10 \log \left( \frac{M_{200}}{M_{\odot}} \right). \quad (20) $$

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

## 4.7. Marginal distributions from the MCMC analysis

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

k-dimensional parameter space to with probability

$$ P = \text{Min} \left[ \frac{p(\theta_{\text{new}})}{p(\theta_{\text{old}})}, 1 \right], \quad (21) $$

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

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

# 5. Results

## 5.1. Model comparison

### 5.1.1. Preamble

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

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

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

### 5.1.2. Bayesian evidence methods

Using different priors can lead to different results, so MAMPOSSSt results should be analyzed carefully. The runs leading to the highest likelihood $L_{MLE}$ naturally tend to have the largest
<table><caption>**Table 2.** MAMPOSStr runs.</caption><thead><tr><th rowspan="2">Model<br>(1)</th><th colspan="2">Mass model</th><th>Anis.</th><th colspan="3">Inner anisotropy</th><th colspan="3">Outer anisotropy</th><th rowspan="2">TAND<br>(11)</th><th rowspan="2">R<sup>-1</sup><br>(12)</th><th rowspan="2">- ln L<sub>MLE</sub><br>(13)</th><th rowspan="2">#<br>(14)</th><th rowspan="2">AIC<br>(15)</th><th rowspan="2">BIC<br>(16)</th></tr><tr><th>Cluster<br>(2)</th><th>BCG<br>(3)</th><th>model<br>(4)</th><th>E<br>(5)</th><th>S0<br>(6)</th><th>S<br>(7)</th><th>E<br>(8)</th><th>S0<br>(9)</th><th>S<br>(10)</th></tr></thead><tbody><tr><td>1</td><td>gNFW</td><td>-</td><td>T</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.004</td><td><strong>33526.28</strong></td><td>12</td><td>67076.62</td><td>67154.32</td></tr><tr><td>2</td><td>NFW</td><td>NFW</td><td>T</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>N</td><td><strong>0.065</strong></td><td>33526.47</td><td>16</td><td>67085.05</td><td>67188.62</td></tr><tr><td>3</td><td>NFW</td><td>NFW</td><td>T</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.012</td><td>33526.50</td><td>13</td><td>67079.08</td><td>67163.24</td></tr><tr><td>4</td><td>gNFW</td><td>-</td><td>T</td><td>F</td><td>0</td><td>0</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.005</td><td>33526.50</td><td>10</td><td><strong>67073.05</strong></td><td>67137.80</td></tr><tr><td>5</td><td>NFW</td><td>NFW</td><td>gOM</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.007</td><td>33526.55</td><td>13</td><td>67079.18</td><td>67163.34</td></tr><tr><td>6</td><td>gNFW</td><td>-</td><td>gOM</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.011</td><td>33526.79</td><td>12</td><td>67077.65</td><td>67155.34</td></tr><tr><td>7</td><td>gNFW</td><td>-</td><td>T</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>N</td><td><strong>0.040</strong></td><td>33526.91</td><td>15</td><td>67083.92</td><td>67181.01</td></tr><tr><td>8</td><td>NFW</td><td>PS4</td><td>T</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.005</td><td>33528.33</td><td>13</td><td>67082.74</td><td>67166.90</td></tr><tr><td>9</td><td>NFW</td><td>-</td><td>T</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.002</td><td><strong>33528.36</strong></td><td>11</td><td><strong>67078.77</strong></td><td>67150.00</td></tr><tr><td>10</td><td>NFW</td><td>-</td><td>T</td><td>F</td><td>0</td><td>0</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.001</td><td>33528.41</td><td>9</td><td>67074.86</td><td>67133.14</td></tr><tr><td>11</td><td>gNFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.005</td><td>33528.41</td><td>9</td><td>67074.86</td><td>67133.14</td></tr><tr><td>12</td><td>NFW</td><td>-</td><td>T</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>N</td><td><strong>0.031</strong></td><td>33528.50</td><td>14</td><td>67085.09</td><td>67175.72</td></tr><tr><td>13</td><td>NFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.002</td><td>33528.54</td><td>8</td><td>67073.11</td><td>67124.92</td></tr><tr><td>14</td><td>NFW</td><td>PS4</td><td>T</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>N</td><td><strong>0.043</strong></td><td>33528.55</td><td>16</td><td>67089.21</td><td>67192.77</td></tr><tr><td>15</td><td>NFW</td><td>-</td><td>gOM</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>Y</td><td>0.003</td><td>33528.92</td><td>11</td><td>67079.90</td><td>67151.12</td></tr><tr><td>16</td><td>NFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>F</td><td>Y</td><td>0.003</td><td>33529.20</td><td>7</td><td><strong>67072.42</strong></td><td>67117.76</td></tr><tr><td>17</td><td>E (free n)</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>F</td><td>Y</td><td>0.002</td><td>33529.74</td><td>6</td><td>67073.50</td><td>67118.84</td></tr><tr><td>18</td><td>gNFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>N</td><td>0.015</td><td>33529.90</td><td>8</td><td>67075.83</td><td>67127.64</td></tr><tr><td>19</td><td>gNFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>F</td><td>0</td><td>0</td><td>F</td><td>Y</td><td>0.003</td><td>33530.03</td><td>8</td><td>67076.09</td><td>67127.90</td></tr><tr><td>20</td><td>NFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>N</td><td>0.007</td><td>33530.23</td><td>7</td><td>67074.48</td><td>67119.82</td></tr><tr><td>21</td><td>gNFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>Y</td><td>0.002</td><td>33530.27</td><td>7</td><td>67074.56</td><td><strong>67119.90</strong></td></tr><tr><td>22</td><td>NFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>F</td><td>0</td><td>0</td><td>F</td><td>Y</td><td>0.002</td><td>33530.35</td><td>7</td><td>67074.72</td><td>67120.06</td></tr><tr><td>23</td><td>E (n = 6)</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>Y</td><td>0.002</td><td>33530.50</td><td>6</td><td>67073.02</td><td><strong>67111.88</strong></td></tr><tr><td>24</td><td>NFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>Y</td><td>0.001</td><td>33530.68</td><td>6</td><td>67073.38</td><td><strong>67112.24</strong></td></tr><tr><td>25</td><td>cNFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>Y</td><td>0.003</td><td>33532.30</td><td>6</td><td>67076.62</td><td>67115.48</td></tr><tr><td>26</td><td>H</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>Y</td><td>0.001</td><td>33534.44</td><td>6</td><td>67080.90</td><td>67119.76</td></tr><tr><td>27</td><td>NFW</td><td>-</td><td>gOM</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>Y</td><td>0.002</td><td>33537.09</td><td>6</td><td>67086.20</td><td>67125.06</td></tr><tr><td>28</td><td>NFW</td><td>-</td><td>iso</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>0.002</td><td>33538.52</td><td>5</td><td>67087.05</td><td>67119.44</td></tr><tr><td>29</td><td>NFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>0</td><td>F</td><td>0</td><td>Y</td><td>0.001</td><td>33539.46</td><td>6</td><td>67090.94</td><td>67129.80</td></tr><tr><td>30</td><td>NFW</td><td>-</td><td>T</td><td>0</td><td>0</td><td>0</td><td>F</td><td>0</td><td>0</td><td>Y</td><td>0.002</td><td>33539.56</td><td>6</td><td>67091.14</td><td>67130.00</td></tr></tbody></table>

**Notes.** All runs are for the *sigv* stack with the concentration-mass relation of Eq. (20). The runs are listed in decreasing maximum likelihood (really posterior), i.e., increasing $-\ln L_{MLE}$. The columns are as follows. (1): run number; (2): cluster mass model (i.e., “E” for Einasto, “H” for Hernquist); (3): brightest cluster galaxy mass model (“PS4” for $n = 4$ Prugniel & Simien 1997, which is a good approximation of the deprojection of a Sersic 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 (“Y” = yes, “N” = no); (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.

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):

$$AIC = -2 \ln L_{MLE} + 2 N_{pars} \quad (22)$$

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

$$BIC = -2 \ln L_{MLE} + \ln N_{data} N_{pars}. \quad (23)$$

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 (i.e., lower) BIC value, while with AIC it must decrease by only 1.

According to Kass & Rafferty (1995), when a model has a BIC lower than another model by 2, 6, or 10 units, one can conclude that there is respectively positive, strong, or very strong

evidence in favor of the former one. Since BIC penalizes extra parameters much more than AIC, BIC seems preferable to AIC when our model is built with a small number of parameters. Equations (22) and (23) lead to $\Delta AIC = \Delta BIC - (\ln N_{data} - 2)\Delta N_{pars}$. Therefore, the condition for strong AIC evidence for the simpler model compared to another one with $\Delta N_{pars}$ extra free parameters, and given our data sample sizes (see Table 1), is $\Delta AIC > 6 - 6.45 \Delta N_{pars}$, that is $\Delta AIC > -0.45$ compared to a more complex model with a single extra parameter. Conversely, a more complex model is strongly favored if $\Delta BIC < -6$, leading to $\Delta AIC < -6 - 6.45 \Delta N_{pars}$, that is $\Delta AIC < -12.45$ compared to a simpler model with one less free parameter.

### 5.1.3. Preferred models

The models with the highest likelihoods from the MCMC analysis (lowest $-\ln 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 $\gamma < -1$) than NFW at inner radii. Note that the free Einasto model with free velocity anisotropy with TAND (not shown in Table 2) has an even lower $-\ln L_{MLE}$ than model 1 (though only slightly). In comparison, our models 2 and 3, which have an additional mass component for a central BCG, also fit well as they rank second and third.

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

### 5.1.4. Rejection of mass models

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

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

There is also strong BIC evidence against replacing the $n = 6$ Einasto and NFW cluster mass models with a free Einasto index model or a gNFW (free inner slope) model, as both models (17 and 21) have higher BICs by over 7 for the same velocity anisotropy priors. The evidence against a gNFW cluster mass model is even stronger ($\Delta$BIC > 10) in moderately complex anisotropy models 11 vs. 13, but the evidence against gNFW is only weak for more complex anisotropy models, for example 1 vs. 9 and 4 vs. 10. And while the best AIC is reached for NFW cluster model 16, the second best model for AIC is gNFW model 4, which has $\Delta$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 $\Delta$BIC = 13 between models 3 and 9 as well as between models 2 and 12. The posteriors and BIC values are even worse (i.e., higher) when assuming that the BCG mass follows the stellar mass, with an $n = 4$ Prugniel & Simien (1997) model (an excellent approximation to the deprojection of the de Vaucouleurs 1948 surface density model): as seen by comparing models 8 and 3. There is also strong evidence in favor of the gNFW cluster mass model compared to the NFW cluster plus NFW BCG (which has an extra free parameter), with $\Delta$BIC $\approx 9$ between models 3 and 1 and $\Delta$BIC > 7 between models 2 and 7. These conclusions are unaltered when using AIC in place of BIC.

### 5.1.5. Rejection of velocity anisotropy models

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

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

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

### 5.1.6. Summary of model comparison

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

## 5.2. Best-fitting parameters

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

<footnote n="9">⁹ We found an even better BIC with the five-parameter priors of NFW or $n=6$ Einasto mass and velocity anisotropy model $\beta(r) = (1/2) r/(r+r_\beta)$ proposed by Mamon & Łokas (2005), which is a special case of the T model, with $\beta \to 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).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** MAMPOSSt marginal distributions (diagonal panels) and covariances for model 24 (non-Einasto model with the strongest BIC evidence) in sigv-stacked pseudo-cluster. 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 the priors are Gaussian distributions with means in the middle and extending to $\pm 3\sigma$ on the panel edges.</figcaption></figure>

mass model with the strongest AIC and BIC evidences (model 3).

Figure 3 shows the MAMPOSSt MCMC posterior marginal distributions and covariances for model 24. All parameters are well fit within their allowed range (except that the outer anisotropy of spirals can reach the physical limit of pure radial orbits). The lower left panel of Fig. 3 shows that the virial radius is anti-correlated with the outer anisotropy of the spiral population (recall that E and S0 galaxies are assumed here isotropic). The mass-anisotropy degeneracy is more acute on the anisotropy (0.5 on $\beta_{sym}$ amounting to 0.11 dex on $(v_r^2)^{1/2}/\sigma_\theta$) than on the virial mass ($3 \times 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 is similar to Fig. 3, but for model 11, with gNFW mass and $T_0$ anisotropy for all morphological types. This figure indicates that the inner density slope is weakly correlated with the mass normalization (left panel of third row) and concentration of the mass profile (second panel of third row), as well as with the outer anisotropies for all three morphological types (third panel of the three bottom rows), which are all three correlated.

<table><caption>**Table 3.** Best-fitting parameters and uncertainties from MAMPOSSt for models 1 (most likely), 21 (strongest BIC evidence among models with gNFW mass), and 24 (strongest BIC evidence non-Einasto model).</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Range</th><th>Model 24</th><th>Model 21</th><th>Model 1</th></tr><tr><td></td><td></td><td></td><td>NFW</td><td>gNFW</td><td>gNFW</td></tr><tr><td></td><td></td><td></td><td>T<sub>0</sub> (S)</td><td>T<sub>0</sub> (S)</td><td>T</td></tr></thead><tbody><tr><td>log r<sub>200</sub></td><td>kpc</td><td>3.0–3.6</td><td>3.23 ± 0.01</td><td>3.23 ± 0.01</td><td>3.19 ± 0.02</td></tr><tr><td>γ</td><td>—</td><td>−2–0</td><td>—</td><td>−1.39<sup>+0.50</sup><sub>−0.33</sub></td><td>−1.64<sup>+0.74</sup><sub>−0.27</sub></td></tr><tr><td>log r<sub>E</sub></td><td>kpc</td><td>2.69 ± 0.10</td><td>2.66 ± 0.08</td><td>2.67 ± 0.08</td><td>2.66 ± 0.08</td></tr><tr><td>log r<sub>S0</sub></td><td>kpc</td><td>2.81 ± 0.11</td><td>2.82 ± 0.08</td><td>2.82 ± 0.08</td><td>2.79 ± 0.08</td></tr><tr><td>log r<sub>S</sub></td><td>kpc</td><td>3.31 ± 0.10</td><td>3.20 ± 0.08</td><td>3.19 ± 0.08</td><td>3.27 ± 0.08</td></tr><tr><td>β<sup>E</sup><sub>sym,0</sub></td><td>—</td><td>−1.8–1.8</td><td>—</td><td>—</td><td>−0.66 ± 0.44</td></tr><tr><td>β<sup>S0</sup><sub>sym,0</sub></td><td>—</td><td>−1.8–1.8</td><td>—</td><td>—</td><td>−0.18 ± 0.39</td></tr><tr><td>β<sup>S</sup><sub>sym,0</sub></td><td>—</td><td>−1.8–1.8</td><td>—</td><td>—</td><td>−0.08 ± 0.17</td></tr><tr><td>β<sup>E</sup><sub>sym,∞</sub></td><td>—</td><td>−1.8–1.8</td><td>—</td><td>—</td><td>1.19 ± 0.55</td></tr><tr><td>β<sup>S0</sup><sub>sym,∞</sub></td><td>—</td><td>−1.8–1.8</td><td>—</td><td>—</td><td>0.95 ± 0.41</td></tr><tr><td>β<sup>S</sup><sub>sym,∞</sub></td><td>—</td><td>−1.8–1.8</td><td>1.24 ± 0.31</td><td>1.11 ± 0.32</td><td>1.61 ± 0.33</td></tr></tbody></table>

**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 $\sigma$ equal to the uncertainty and cut at $\pm 3\sigma$. The quoted values for the three 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 both E and S0 galaxies.

Table 3 shows the MLE values and the uncertainties from the marginal distributions derived from the MAMPOSSt MCMC for models 1, 21, and 24. For model 24 with an NFW mass model, the virial radius is very well measured leading to an MLE value of $r_{200} = 1690 \pm 20$ kpc, thus 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 \pm 23$ kpc, which is still consistent with the Clean value. On the other hand, the gNFW model 1 leads to $r_{200} = 1507 \pm 59$ kpc, which is significantly smaller than the Clean value. For models 21 and 1, the inner slope is consistent with the $-1$ value for NFW.

When it is a free parameter, the mass concentration of model 24 is $c_{200} = 3.8 \pm 0.4$ Given the mass at $r_{200} = 1698 \pm 24$ kpc of $10^{14.8} M_\odot$, 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} \approx 4.2$ Okabe & Smith 2016, which is the median of the 13 measurements within a mass range of 0.3 dex, see Fig. 12 of Child et al. 2018). We further discuss the $c-M$ relation in Sect. 6.1.3.

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

## 5.3. Goodness of fit

Figure 4 shows the LOS velocity dispersion profiles for the elliptical, lenticular, and spiral galaxies predicted from model 1

<footnote n="10">10 This “sweet spot” at the virial radius where mass is most independent of velocity anisotropy is a blessing from nature allowing us to use the cluster mass function derived from spectroscopic surveys as a cosmological tool.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Best-fit, line-of-sight velocity dispersion profiles for model 1 (gNFW with free T outer anisotropy with TAND) for sigv-stacked pseudo-cluster. The symbols are the data (150 galaxies per radial bin) as in Fig. 2, while the curves and shaded regions respectively show the median predictions of model 1 and their MCMC uncertainties. The vertical gray-shaded region represents $r_{200}$ and its MCMC uncertainty.</figcaption></figure>

(gNFW mass with T anisotropy for the three morphological types, all with TAND). The MAMPOSSt model predictions reproduce the data very well.

## 5.4. Radial profiles

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

Figure 5 displays the mass density profiles for models 21 (gNFW) and 3 (NFW + NFW for BCG). In the top panel, a gNFW model was assumed by MAMPOSSt, and the density profile prefers to be steeper than NFW, but not significantly ($\gamma = -1.51 \pm 0.42$ according to Table 3). Only 85% of all chain elements past burn-in produce $\gamma < -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$) that is only 1.6 higher than for model 24. Considering the cNFW model to be a physical one, it has the same number of free parameters as the NFW model and its BIC is only 3.2 higher than that of model 24. So, one cannot reject the cored NFW model for clusters with our WINGS data.

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

One may also wonder which morphological type has a number density profile closest to the mass density profile. Figure 6 displays the ratios of mass density over number density for the three morphological types for models 24 and 21. We normalize the number density profiles by eliminating $N(r_v)$ between Eq. (7) and the average number of galaxies of a given morphological type per cluster, $N_{tot}$, between the minimum and maximum allowed projected radii, $R_{min}$ and $R_{max}$, which we

<figure><img src="image_6.png" /><figcaption>Fig. 5. Radial mass density profiles for model 21 (gNFW with isotropic orbits for the ellipticals and S0s, and $T_0$ velocity anisotropy for spirals, top) and model 3 (NFW cluster + NFW BCG, with T velocity anisotropy, bottom) for sigv-stacked pseudo-cluster. 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 gray) 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, in order 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.</figcaption></figure>

model as

$$N_{\text{tot}} = N_{\text{p}}(R_{\text{max}}) - N_{\text{p}}(R_{\text{min}}) \quad (24)$$

$$= N(r_v) [\tilde{N}_{\text{p}}(R_{\text{max}}) - \tilde{N}_{\text{p}}(R_{\text{min}})], \quad (25)$$

where

$$\tilde{N}_{\text{p}}(X) = \frac{1}{\ln 2 - 1/2} \frac{C^{-1}(1/X)}{\sqrt{|X^2 - 1|}} + \ln\left(\frac{X}{2}\right) \quad (26)$$

for the NFW model, with $N_p(1) = (1 - \ln 2)/(\ln 2 - 1/2)$ and where $C^{-1}$ is given in Eq. (10).

It should be recalled that the NFW model was assumed for the number density profiles of the three morphological types and that these number density profiles were obtained from fits to the photometric data, and thus do not suffer from any spectroscopic incompleteness. In both panels of Fig. 6, the elliptical galaxies trace the mass almost perfectly; S0 galaxies are nearly as good mass tracers as ellipticals (but somewhat more extended), while spirals are much more extended. In runs with free concentration, the ellipticals trace the mass even better. Since model 21 (see
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Radial profiles of mass over number density ratios from MAM-POSSSt for model 24 (NFW mass profile, *top*), and model 21 (gNFW mass profile, *bottom*), for E (red), S0 (green), and S (blue) galaxies, for sigv-stacked pseudo-cluster. The normalization is explained in Eqs. (24)–(26). The horizontal lines are shown to highlight how well mass follows number.</figcaption></figure>

the bottom panel of Fig. 6) is based on a gNFW mass profile, there is less agreement between the NFW number density profile of the ellipticals and the gNFW mass profile, but ellipticals remain the best tracers of mass among the three morphological types. Indeed, the mass over elliptical number density ratio is nearly consistent with being constant (horizontal line in Fig. 6), although there may be a need for extra mass in the BCG.

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

## 5.5. Outer versus inner velocity anisotropies

Table 4 illustrates the details of the anisotropy for four models. If we free the inner and outer anisotropies (but tie the anisotropy radii to the scale radii, as in model 1), we find that the anisotropy at $r_{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 $\beta_{sym}$ for

<table><caption>**Table 4.** Velocity anisotropies from MAMPOSSSt at 0.03 $r_{200}$ and at $r_{200}$.</caption><thead><tr><th rowspan="2">Model<br>(1)</th><th rowspan="2">Mass<br>(2)</th><th rowspan="2">β(r)<br>(3)</th><th rowspan="2">TAND<br>(4)</th><th rowspan="2">Type<br>(5)</th><th colspan="3">β<sub>sym</sub></th><th rowspan="2">&#9650;<br>(8)</th></tr><tr><th>(0.03 r<sub>200</sub>)<br>(6)</th><th>(r<sub>200</sub>)<br>(7)</th></tr></thead><tbody><tr><td rowspan="3">1</td><td rowspan="3">gNFW</td><td rowspan="3">T</td><td rowspan="3">Y</td><td>E</td><td>-0.33<sup>+0.49</sup><sub>-0.29</sub></td><td>0.22 ± 0.24</td><td>0.81</td></tr><tr><td>S0</td><td>0.06<sup>+0.41</sup><sub>-0.30</sub></td><td>0.36 ± 0.20</td><td>0.73</td></tr><tr><td>S</td><td>0.01 ± 0.17</td><td>0.48 ± 0.14</td><td>1.00</td></tr><tr><td rowspan="3">3</td><td rowspan="3">2NFW</td><td rowspan="3">T</td><td rowspan="3">Y</td><td>E</td><td>-0.12<sup>+0.39</sup><sub>-0.31</sub></td><td>0.14 ± 0.26</td><td>0.70</td></tr><tr><td>S0</td><td>0.21<sup>+0.40</sup><sub>-0.31</sub></td><td>0.30 ± 0.20</td><td>0.58</td></tr><tr><td>S</td><td>-0.02 ± 0.17</td><td>0.47 ± 0.14</td><td>0.99</td></tr><tr><td rowspan="3">7</td><td rowspan="3">gNFW</td><td rowspan="3">T</td><td rowspan="3">N</td><td>E</td><td>-0.11<sup>+0.39</sup><sub>-0.46</sub></td><td>0.14 ± 0.21</td><td>0.71</td></tr><tr><td>S0</td><td>0.16<sup>+0.37</sup><sub>-0.49</sub></td><td>0.31 ± 0.19</td><td>0.63</td></tr><tr><td>S</td><td>-0.15<sup>+0.36</sup><sub>-0.53</sub></td><td>0.49 ± 0.15</td><td>0.97</td></tr><tr><td>24</td><td>NFW</td><td>T<sub>0</sub></td><td>Y</td><td>S</td><td>0.02 ± 0.01</td><td>0.50 ± 0.12</td><td>1.00</td></tr></tbody></table>

**Notes.** Columns are: (1): model number (see Table 2); (2): mass model (2NFW stands for NFW + NFW (BCG)); (3): anisotropy model; (4): anisotropy radius set equal to the tracer scale radius (TAND); (5): morphological type; (6): $\beta_{sym}$ (Eq. (19)) at $r = 0.03 r_{200}$; (7): $\beta_{sym}$ at $r = r_{200}$; (8): fraction of anisotropy profiles that increase with radius. The uncertainties are $(p_{84} - p_{16})/2$, where $p_i$ are the ith percentiles.

E vs. S) for these early-type galaxies compared to spirals (Figs. 7 and 8), which explains why BIC evidence prefers having radial outer orbits for the spiral population only. Hence, there is only marginal evidence that the S0 population has radial outer orbits, while the moderately radial orbits of ellipticals is not statistically significant (Table 4). The outer anisotropies of the elliptical and S0 galaxies are less radial when the anisotropy radii are set free (model 7) compared to the analogous TAND model 1.

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

Figure 8 provides a clearer way to view the anisotropy profiles, by plotting the value at $r_{200}$ as a function of the value at 0.03 $r_{200}$. We restricted these plots to models with free inner and outer anisotropies for all morphological types. When the anisotropy radius is forced to the scale radius, as favored by the Bayesian (both AIC and BIC) evidence (top panel of Fig. 8), the 95 percent confidence contours for $\beta(r_{200})$ for spirals are above zero for all values of the inner anisotropy (at 0.03 $r_{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 three morphological types. The free vs. fixed anisotropy radii have a stronger effect on the contours of outer vs. inner anisotropy than does the mass model. Nevertheless, only spiral galaxies show clearly radial anisotropy at $r_{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
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Velocity anisotropy (Eq. (19)) profiles of E, S0, and S galaxies from MAMPOSSt for model 24 (NFW with T-TAND anisotropy for spirals only, *upper left*), model 9 (NFW with T-TAND anisotropy, *upper right*), model 1 (gNFW with T-TAND anisotropy, *lower left*), and model 7 (gNFW with T anisotropy and free anisotropy radius, *lower right*) for sigv-stacked pseudo-cluster. The hashed regions indicate the 68% confidence zone, while the curves display the 5th and 95th percentiles. The thick vertical gray 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}$.</figcaption></figure>

to span between 1 and 4, the uncertainty on the best-fit anisotropy radii for the non-TAND runs are typically as high as 1 dex for all three morphological types. Nevertheless, as for the TAND case, spirals are the sole morphological type for which the orbits systematically become more radial from the inner regions to near the virial radius (last column of Table 4, and as seen by the contours for spirals lying above the oblique line in Fig. 8).

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 $\Delta$BIC = 6 of the lowest BIC), or we can analyze the detailed radial profiles expected from the models that reach the highest likelihoods (which are really posteriors), although their BIC Bayesian evidence is so high that they can be strongly rejected relative to the lowest BIC model.

# 6. Discussion

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

## 6.1. Mass density profiles of WINGS clusters

### 6.1.1. General trends

Our highest BIC Bayesian evidence is reached for models 23 and 24, where the mass profile is $n = 6$ Einasto or NFW with isotropic velocities for ellipticals and S0s, while for the spirals they are isotropic at the center and fairly radial in the outer regions of clusters (Table 2). There is strong Bayesian evi-
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Outer ($r = r_{200}$) vs. inner ($r = 0.03 r_{200}$) velocity anisotropy (Eq. (19)) from MAMPOSS for three morphological types and for model 1 (gNFW-T anisotropy with TAND), model 7 (same as model 1, but without TAND), 9 (NFW-T-TAND), and model 12 (same as model 9, but without TAND) for sigv-stacked pseudo-cluster. The contours are 68% and 95% confidence (with pixel resolution $\Delta\beta_{\text{sym}} = 0.03$ and smoothed with a Gaussian kernel of $\sigma = 2$ pixels).</figcaption></figure>

dence against the Hernquist model (26), whose outer slope is steeper (–4) then that of NFW (–3). There is only positive (but not strong) evidence against a cored NFW model (relative to model 24). The case against the gNFW and free index Einasto models is less clear with AIC evidence: the lowest AIC value is reached for an NFW model, but the second lowest among non-Einasto models is for a gNFW model.

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

stack them. We thus simply do not have enough tracers to constrain the mass density profile within this radius.

### 6.1.2. Robustness

We ran several models for the two other cluster stacks (Num and tempX). We found that for model 21 (which uses gNFW for the cluster mass), the constraints on the inner slope range from $\gamma = -1.4_{-0.3}^{+0.5}$ for sigv to $-1.7_{-0.2}^{+0.3}$ for tempX and $-1.8_{-0.1}^{+0.4}$ for Num. The steeper inner mass density slopes for the tempX and Num stacks lead to different Bayesian evidence for or against gNFW. Indeed, comparing model 24 (NFW) to 21, we find that AIC prefers gNFW for Num ($\Delta$AIC = AIC(gNFW) – AIC(NFW) = –5.9) and tempX ($\Delta$AIC = –2.5), whereas it slightly favors NFW for sigv ($\Delta$AIC = 1.2). However, BIC does not favor gNFW despite the much steeper gNFW inner slopes. While there is strong BIC evidence against gNFW with sigv ($\Delta$BIC > 7), there is still positive evidence against gNFW with tempX ($\Delta$BIC = 3.7) and Num ($\Delta$BIC = 0.5).
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Same as Fig. 6 but for different cluster stacking methods, again for models 24 and 21. Only median trends are shown, for clarity.</figcaption></figure>

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

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

### 6.1.3. Comparison with other work

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

Figure 10 shows the constraints on the concentration-mass relation obtained with NFW model 24. Interestingly, with free concentration (i.e., flat prior $0 < \log c_{200} < 1$), the magenta contours are well matched to the $\Lambda$CDM relation of Dutton & Macciò (2014) given in Eq. (20), especially for stacks tempX and sigv. It is not surprising that folding in this relation as a prior, we recover similar contours that are closer to this relation. The concentrations obtained with MAMPOSSt are well matched to those obtained by weak lensing, except for one (fairly old

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Cluster mass-concentration vs. mass for all three stacks with model 24. The magenta contours indicate the run with free concentration (flat prior $0 < \log c_{200} < 1$), while the yellow contours display the run with the $\Lambda$CDM relation of Eq. (20), highlighted in the shaded regions for 1 and $2\sigma$ constraints. The contours are 68% constraints. The points are the weak-lensing analyses by Johnston et al. (2007) (triangles), Mandelbaum et al. (2008) (downwards triangles), Okabe et al. (2010) (curly squares), Oguri et al. (2012) (cross), Okabe et al. (2013) (open square), Sereno & Covone (2013) (open circle), Umetsu et al. (2014) (curly diamond), Umetsu et al. (2016) (diamond), Okabe & Smith (2016) (filled square), Umetsu & Diemer (2017) (filled diamond), and Cibirka et al. (2017) (filled circle), all corrected to be $(1+z)^{0.38}c_{200}$ following Child et al. (2018). The error bar at the bottom is from Mandelbaum et al. (2008).</figcaption></figure>

weak lensing study that strongly underestimates the concentration. Biviano et al. (2017) performed a MAMPOSSt analysis of 49 WINGS clusters, individually, and found a cluster-to-cluster scatter in concentration that was greater than the uncertainties returned from MAMPOSSt ($\approx 0.3$ dex) and from the effect of a range of cluster masses ($\approx 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 (Katgert et al. 2004 for ENACS clusters) galaxies. The much weaker concentration of the spiral population agrees with the much weaker concentration of the blue galaxies relative to the red ones observed by Collister & Lahav (2005).

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

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

## 6.2. Velocity anisotropy profiles of WINGS clusters

### 6.2.1. General trends

For the single component mass models, all anisotropy models that differ from that of $T_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 model 16 (the best one using AIC evidence) where S0 galaxies and spirals have outer anisotropy (Table 2). But the Bayesian evidence of model 24 (with outer anisotropy only for the spirals) against model 16 is $\Delta$BIC = 5.5, that is “positive” and almost “strong” ($\Delta$BIC > 6). Hence, we are confident 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.

### 6.2.2. Robustness

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

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

Figure 13 compares the outer vs. inner velocity anisotropies when we change the minimum number of member galaxies in

<figure><img src="image_12.png" /><figcaption>Fig. 11. Same as Fig. 8, but comparing priors on cluster concentration for models 1 (TAND, top) and 7 (free anisotropy radius, bottom). Only 68% contours are shown.</figcaption></figure>

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><img src="image_13.png" /><figcaption>**Fig. 12.** Top: same as Fig. 11 (but for model 1) for different minimum projected radii 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_{\text{min}} = 50$ kpc. Bottom: same as top panel, for different maximum projected radii for galaxy selection: $0.67r_{200} \approx r_{100}/2$ (thin) and $1.35r_{200} \approx r_{100}$ (thick, our standard case).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Same as Fig. 12, but for different minimum numbers of members in individual clusters used for the stacks: 30 (thick) and 81 (thin).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Same as Fig. 12, but for the all three stacks of clusters according to different estimates of 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 mass-temperature relation from X-ray measurements (tempX, thin).</figcaption></figure>

### 6.2.3. Comparison with previous studies

Figure 15 compares our constraints on inner and outer anisotropy with those from the literature. We first compare to previous studies that did not separate galaxies into different classes. The two anisotropy measures of Abell 2218 by Natarajan & Kneib (1996) are consistent with our anisotropies. However, the inner anisotropies of the 10 clusters measured by Hwang & Lee (2008) (only two are shown in Fig. 15: Abell 85 for comparison with another study and Abell 1689 which has the richest data sample) are much more radial than we (and most others) found. Aguerri

et al. (2017) and Benatov et al. (2006) found much lower inner radial anisotropy for Abell 85 and Abell 2199, respectively, than Hwang & Lee (2008). Still, Benatov et al. (2006) and Wojtak & Łokas (2010) both find a very wide range of anisotropies at $r_{200}$, some of which are 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
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Comparison with previous measurements (shown in symbols) of 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 sigv-stacked pseudo-cluster. Top panels: comparison with the full literature over the entire possible range of anisotropies (for Hwang & Lee (2008), only Abell 85 (black) and Abell 2199 (red and blue) are shown). Bottom panels: zooms, restricted to previous works that differentiated between galaxy types (with the same symbol 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 early-type or passive, blue for late-type or star forming, and black for all galaxies (mixed types). 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 ($\beta(0.03 r_{200}) = 0.2$) and for the later type spirals ($\beta(0.03 r_{200}) = -1.7$). The two for Annunziatella et al. refer to low ($\beta(0.03 r_{200}) \approx -0.6$) and high ($\beta(0.03 r_{200}) \approx 0.1$) mass galaxies. Finally, the two blue symbols for Aguerri et al. are for all ($\beta(r_{200}) \approx 0.5$) and dwarf ($\beta(r_{200}) \approx 0$) galaxies.</figcaption></figure>

Munari et al. (2014) on red and blue galaxies in Abell 2142 are marginally consistent with our respective results on E+S0 and S galaxies, with their anisotropy of the red galaxies at $r_{200}$ appearing too tangential. There is tension between the very radial ($\beta > 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) and our much lower values of $\beta$ for the ellipticals. On the other hand, the inner and outer isotropy of passive

galaxies found by Capasso et al. (2019) is consistent with the orbits we find for E and S0 galaxies. The anisotropies found in Abell 85 by Aguerri et al. (2017) for blue galaxies are consistent with ours, except that the outer isotropy that they found for their blue dwarf galaxies is in some tension with our results for spirals. On the other hand, they found significantly more radial orbits for red galaxies in Abell 85 compared to our findings for early type galaxies.
The discrepancies in orbital anisotropies between studies of individual clusters and our stacked analysis may be caused by a possible diversity of clusters, either natural or modulation with mass or redshift (the alternative is that the methods and/or pri-ors were different). It is therefore interesting to compare the anisotropies from the analyses of stacked or clusters or from joint analyses of clusters.

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

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

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

### 6.2.4. Infall

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

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

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

### 6.2.5. Isotropization

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

The natural way for them to isotropize is by two-body relaxation with other galaxies. The typical timescale for two-body relaxation roughly scales as $N/(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_{\odot}$, Sect. 3.1), the mean number of galaxies per cluster in our sample (87), and for reasonable choices of the pericenter radius and the ratio of apocenter to pericenter, we find (see Fig. D.1) that at the very least 30 orbits are required to isotropize. In considering a growing NFW cluster, Tollet et al. (2017) found that (see Fig. B1 therein) it takes 3 Gyr for a galaxy to move from pericenter to its second apocenter (assuming hereafter that $r_{\text{apo}}/r_{\text{vir}}(t_{\text{apo}}) = 3.5$), and thus over 4 Gyr, between the last two pericentric passages. These orbital times were shorter at early times, but the number of orbits is expected to be less than a dozen. This suggests that two-body relaxation suffered by a single galaxy is insufficient to explain the isotropy of the early-type galaxies in the inner regions (or further out in some of the stacks).

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

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

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

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

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

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Illustration of wide (salmon) and narrow (blue) ranges of pericenters, leading to different mixes of velocity anisotropies at $r = r_{min}$ (light green). At pericenter and apocenter, orbits are necessarily circular, hence $\beta_{sym} = -2$. 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}$.</figcaption></figure>

transformation as soon as their first entry into the cluster virial sphere.

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

If spirals orbit only once around the clusters with their original morphology, their range of log pericenters will be much narrower than if they orbit many times. Those that fell in the cluster long ago at very small pericenters (given their small apocenters, as was the virial radius of the cluster's most massive progenitor) should have different morphologies now<fnref n="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 fully tangential (i.e., circular) at $r = r_{peri}$, rapidly increasing with radius to the radial values caused by infall (if the apocenters were all equal, one would return to circular at $r = r_{apo}$). This rapid transition in the velocity anisotropy profile is seen in recent hydrodynamical simulations (Lotz et al. 2019). But, as illustrated in Fig. 16, if the spirals only orbit once through the cluster, the radial orbits at $r$ contribute less in comparison to the quasi-circular orbits for $r \ge r_{peri}$, leading to a more isotropic velocity distribution at $r$.

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

At high redshift, the orbital times are shorter, while the morphological transformation times should be roughly the same. Hence high redshift spirals may survive several orbits in their clusters, and we would then predict that the inner velocity

<footnote n="11">¹¹ If spirals lose their morphology in a single orbit, their number density profile is no longer an NFW model as assumed, but a similar model truncated at small radii. This does not affect our analysis since the integrations in MAMPOSS are outwards.</footnote>
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 protoclusters.

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

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

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

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

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

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

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

<table><caption>**Table 5.** Pearson correlation coefficients of inner-outer anisotropies between morphological types (for model 1).</caption><thead><tr><th>Stack</th><th>E vs. S0</th><th>E vs. S</th><th>S0 vs. S</th></tr></thead><tbody><tr><td>sigv</td><td>0.51</td><td>0.23</td><td>0.69</td></tr><tr><td>Num</td><td>0.72</td><td>0.17</td><td>0.50</td></tr><tr><td>tempX</td><td>0.44</td><td>0.18</td><td>0.82</td></tr></tbody></table>

correlation coefficient $C_{1,2}$ between two types as

$$C_{1,2} = \frac{\langle\langle (f_1 - \langle f_1 \rangle)(f_2 - \langle f_2 \rangle)\rangle\rangle}{\sigma(f_1)\sigma(f_2)}, \quad (27)$$

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

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

S0s appear to form relatively late in clusters, the fraction of S0s is 0.1 at $z \sim 1$ (Smith et al. 2005; Postman et al. 2005), that is 8 Gyr ago, while our sample has $\sim 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 \to S0$ evolution. One could then argue that orbital isotropization should take longer than the $S \to S0$ evolution, otherwise S0 orbits would have mostly isotropized and resemble those of ellipticals more than those of spirals. Finally, our results on S0s should be taken with caution, because the S0 morphological class is notoriously difficult to cleanly classify for ranges of inclinations and apparent magnitudes (e.g., van den Bergh 2009) and may thus be contaminated by both spirals and ellipticals.

### 6.3. Perspectives

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

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

## 7. Conclusions

We ran the MAMPOSSt mass-orbit modeling algorithm on three stacks of 54 regular, nearby ($z \approx 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 three morphological types in projected phase space, fitting for the shape of the total mass profile on one hand and of the three velocity anisotropy profiles on the other. We ran MAMPOSSt with 30 different sets of priors. Our results for the sigv stack are as follows:

- There is no compelling evidence for a mass density profile steeper than NFW or $n = 6$ Einasto at 0.03 $r_{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 (with an inner slope of roughly $-1.5 \pm 0.5$) or with an NFW profile for the BCG in addition to an NFW profile for the remaining part of the cluster.

- An outer slope as steep as $-4$ (Hernquist model) is ruled out.

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

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

- The velocity anisotropy of spirals rises from near isotropic in the inner regions to mildly radial ($\beta \approx 0.45 \pm 0.08$) at $r_{200}$.

- The velocity anisotropy of the lenticulars also rises from near isotropic in the inner regions to somewhat less radial ($\beta = 0.31 \pm 0.17$) at $r_{200}$ than for the spirals.

- The velocity anisotropy of the ellipticals is consistent with isotropic velocities everywhere, even though the highest likelihoods are reached for slightly tangential inner orbits and mildly radial anisotropy ($\beta = 0.19 \pm 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. This is actually the preferred model using AIC evidence.

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

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

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

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

- The outer anisotropy of spirals is less pronounced.

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

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

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

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

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

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

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

## References

Adami, C., Le Brun, V., Biviano, A., et al. 2009, A&A, 507, 1225

Aguerri, J. A. L., Agulli, I., Diaferio, A., & Dalla Vecchia, C. 2017, MNRAS, 468, 364

Akaike, H. 1973, International Symposium on Information Theory (Budapest: Akadémiae Kiadó), 267

Annunziatella, M., Mercurio, A., Biviano, A., et al. 2016, A&A, 585, A160

Arnaud, M., Pointecouteau, E., & Pratt, G. W. 2005, A&A, 441, 893

Ascasibar, Y., & Gottlöber, S. 2008, MNRAS, 386, 2022

Auger, M. W., Treu, T., Bolton, A. S., et al. 2010, ApJ, 724, 511

Bartelmann, M. 1996, A&A, 313, 697

Beers, T. C., & Tonry, J. L. 1986, ApJ, 300, 557

Beers, T. C., Flynn, K., & Gebhardt, K. 1990, AJ, 100, 32

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Benatov, L., Rines, K., Natarajan, P., Kravtsov, A., & Nagai, D. 2006, MNRAS, 370, 427

Bicknell, G. V., Bruce, T. E. G., Carter, D., & Killeen, N. E. B. 1989, ApJ, 336, 639

Bigiel, F., Leroy, A. K., Walter, F., et al. 2011, ApJ, 730, L13

Binney, J., & Mamon, G. A. 1982, MNRAS, 200, 361

Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton, NJ: Princeton University Press)

Biviano, A., & Girardi, M. 2003, ApJ, 585, 205

Biviano, A., & Katgert, P. 2004, A&A, 424, 779

Biviano, A., & Poggianti, B. M. 2009, A&A, 501, 419

Biviano, A., & Salucci, P. 2006, A&A, 452, 75

Biviano, A., Rosati, P., Balestra, I., et al. 2013, A&A, 558, A1

Biviano, A., van der Burg, R. F. J., Muzzin, A., et al. 2016, A&A, 594, A51
Biviano, A., Moretti, A., Paccagnella, A., et al. 2017, A&A, 607, A81

Capasso, R., Saro, A., Mohr, J. J., et al. 2019, MNRAS, 482, 1043

Cava, A., Bettoni, D., Poggianti, B. M., et al. 2009, A&A, 495, 707

Cava, A., Biviano, A., Mamon, G. A., et al. 2017, A&A, 606, A108 (Paper I)

Chabrier, G. 2003, PASP, 115, 763

Child, H. L., Habib, S., Heitmann, K., et al. 2018, ApJ, 859, 55

Cibirka, N., Cypriano, E. S., Brimiouille, F., et al. 2017, MNRAS, 468, 1092

Collister, A. A., & Lahav, O. 2005, MNRAS, 361, 415

Courteau, S., Cappellari, M., de Jong, R. S., et al. 2014, Rev. Mod. Phys., 86, 47

de Vaucouleurs, G. 1948, Ann. Astrophys., 11, 247

Dejonghe, H., & Merritt, D. 1992, ApJ, 391, 531

Diaferio, A. 1999, MNRAS, 309, 610

Diaferio, A., & Geller, M. J. 1997, ApJ, 481, 633

Diakogiannis, F. I., Lewis, G. F., Ibata, R. A., et al. 2017, MNRAS, 470, 2034

Dressler, A., & Shectman, S. A. 1988, AJ, 95, 985

Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359

Einasto, J. 1965, Trudy Inst. Astroz. Alma-Ata, 51, 87

Fakhouri, O., Ma, C.-P., & Boylan-Kolchin, M. 2010, MNRAS, 406, 2267

Falco, M., Mamon, G. A., Wojtak, R., Hansen, S. H., & Gottlöber, S. 2013, MNRAS, 436, 2639

Fasano, G., Marmo, C., Varela, J., et al. 2006, A&A, 445, 805

Fasano, G., Vanzella, E., Dressler, A., et al. 2012, MNRAS, 420, 926

Foëx, G., Böhringer, H., & Chon, G. 2017, A&A, 606, A122

Gaia Collaboration (Helmi, A., et al.) 2018, A&A, 616, A12

Geller, M. J., Diaferio, A., & Kurtz, M. J. 1999, ApJ, 517, L23

Ghigna, S., Moore, B., Governato, F., et al. 1998, MNRAS, 300, 146

Gill, S. P. D., Knebe, A., & Gibson, B. K. 2005, MNRAS, 356, 1327

Girardi, L., Bressan, A., Chiosi, C., Bertelli, G., & Nasi, E. 1996, A&AS, 117, 113

Gott, J. R. 1975, ApJ, 201, 296

Gunn, J. E., & Gott, III., J. R. 1972, ApJ, 176, 1

Haines, C. P., Pereira, M. J., Smith, G. P., et al. 2015, ApJ, 806, 101

Hernquist, L. 1990, ApJ, 356, 359

Hwang, H. S., & Lee, M. G. 2008, ApJ, 676, 218

Jiang, C. Y., Jing, Y. P., Faltenbacher, A., Lin, W. P., & Li, C. 2008, ApJ, 675, 1095

Johnston, D. E., Sheldon, E. S., Wechsler, R. H., et al. 2007, ArXiv e-prints [arXiv:0709.1159]

Kass, R. E., & Rafferty, A. E. 1995, J. Am. Stat. Assoc., 90, 773

Katgert, P., Biviano, A., & Mazure, A. 2004, ApJ, 600, 657

Kazantzidis, S., Magorrian, J., & Moore, B. 2004, ApJ, 601, 37

Larson, R. B., Tinsley, B. M., & Caldwell, C. N. 1980, ApJ, 237, 692

Lemze, D., Barkana, R., Broadhurst, T. J., & Rephaeli, Y. 2008, MNRAS, 386, 1092

Lemze, D., Wagner, R., Rephaeli, Y., et al. 2012, ApJ, 752, 141

Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511

Łokas, E. L. 2002, MNRAS, 333, 697

Łokas, E. L., & Mamon, G. A. 2001, MNRAS, 321, 155

Łokas, E. L., & Mamon, G. A. 2003, MNRAS, 343, 401

Lotz, M., Remus, R.-S., Dolag, K., Biviano, A., & Burkert, A. 2019, MNRAS, 488, 5370

Lynden-Bell, D. 1967, MNRAS, 136, 101

Macciò, A. V., Dutton, A. A., & van den Bosch, F. C. 2008, MNRAS, 391, 1940

Mahajan, S., Mamon, G. A., & Raychaudhury, S. 2011, MNRAS, 416, 2882

Mamon, G. A. 1992, ApJ, 401, L3

Mamon, G. A., & Boué, G. 2010, MNRAS, 401, 2433

Mamon, G. A., & Łokas, E. L. 2005, MNRAS, 363, 705

Mamon, G. A., Sanchis, T., Salvador-Solé, E., & Solanes, J. M. 2004, A&A, 414, 445

Mamon, G. A., Biviano, A., & Murante, G. 2010, A&A, 520, A30

Mamon, G. A., Biviano, A., & Boué, G. 2013, MNRAS, 429, 3079

Mandelbaum, R., Seljak, U., & Hirata, C. M. 2008, J. Cosmol. Astropart. Phys., 8, 006

Merritt, D. 1985, ApJ, 289, 18

Merritt, D. 1987, ApJ, 313, 121

Moore, B., Katz, N., Lake, G., Dressler, A., & Oemler, A. 1996, Nature, 379, 613

Moretti, A., Poggianti, B. M., Fasano, G., et al. 2014, A&A, 564, A138

Munari, E., Biviano, A., & Mamon, G. A. 2014, A&A, 566, A68

Natarajan, P., & Kneib, J.-P. 1996, MNRAS, 283, 1031

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039

Newman, A. B., Treu, T., Ellis, R. S., et al. 2013a, ApJ, 765, 24

Newman, A. B., Treu, T., Ellis, R. S., & Sand, D. J. 2013b, ApJ, 765, 25

Oguri, M., Bayliss, M. B., Dahle, H., et al. 2012, MNRAS, 420, 3213

Okabe, N., & Smith, G. P. 2016, MNRAS, 461, 3794

Okabe, N., Takada, M., Umetsu, K., Futamase, T., & Smith, G. P. 2010, PASJ, 62, 811

Okabe, N., Smith, G. P., Umetsu, K., Takada, M., & Futamase, T. 2013, ApJ, 769, L35

Old, L., Skibba, R. A., Pearce, F. R., et al. 2014, MNRAS, 441, 1513

Old, L., Wojtak, R., Mamon, G. A., et al. 2015, MNRAS, 449, 1897

Osipkov, L. P. 1979, Sov. Astron. Lett., 5, 42

Postman, M., Franx, M., Cross, N. J. G., et al. 2005, ApJ, 623, 721

Prugniel, P., & Simien, F. 1997, A&A, 321, 111

Read, J. I., & Steger, P. 2017, MNRAS, 471, 4541

Richardson, T., & Fairbairn, M. 2013, MNRAS, 432, 3361

Sanchis, T., Łokas, E. L., & Mamon, G. A. 2004, MNRAS, 347, 1198

Schwarz, G. 1978, Ann. Stat., 6, 461

Sereno, M., & Covone, G. 2013, MNRAS, 434, 878

Sereno, M., Ettori, S., Umetsu, K., & Baldi, A. 2013, MNRAS, 428, 2241

Sersic, J. L. 1968, 1968, Atlas de galaxias australes (Observatorio Astronomico: Cordoba, Argentina)

Skibba, R. A., van den Bosch, F. C., Yang, X., et al. 2011, MNRAS, 410, 417

Skielboe, A., Wojtak, R., Pedersen, K., Rozo, E., & Rykoff, E. S. 2012, ApJ, 758, L16

Smith, G. P., Treu, T., Ellis, R. S., Moran, S. M., & Dressler, A. 2005, ApJ, 620, 78

Solanes, J. M., & Salvador-Solé, E. 1990, A&A, 234, 93

Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685

Tiret, O., Combes, F., Angus, G. W., Famaey, B., & Zhao, H. S. 2007, A&A, 476, L1

Tollet, É., Cattaneo, A., Mamon, G. A., Moutard, T., & van den Bosch, F. C. 2017, MNRAS, 471, 4170

Tonry, J. L. 1983, ApJ, 266, 58

Umetsu, K., & Diemer, B. 2017, ApJ, 836, 231

Umetsu, K., Medezinski, E., Nonino, M., et al. 2014, ApJ, 795, 163

Umetsu, K., Zitrin, A., Gruen, D., et al. 2016, ApJ, 821, 116

van den Bergh, S. 2009, ApJ, 694, L120

van der Marel, R. P. 1994, MNRAS, 270, 271

van der Marel, R. P., Magorrian, J., Carlberg, R. G., Yee, H. K. C., & Ellingson, E. 2000, AJ, 119, 2038

Vasiliev, E. 2019, MNRAS, 482, 1525

Wainer, H., & Thissen, D. 1976, Psychometrica, 41, 9

Walker, M. G., Mateo, M., Olszewski, E. W., et al. 2009, ApJ, 704, 1274

Wetzel, A. R., Tinker, J. L., Conroy, C., & van den Bosch, F. C. 2013, MNRAS, 432, 336

Wojtak, R. 2013, A&A, 559, A89

Wojtak, R., & Łokas, E. L. 2010, MNRAS, 408, 2442

Wojtak, R., & Mamon, G. A. 2013, MNRAS, 428, 2407

Wojtak, R., Łokas, E. L., Mamon, G. A., et al. 2008, MNRAS, 388, 815

Wojtak, R., Łokas, E. L., Mamon, G. A., & Gottlöber, S. 2009, MNRAS, 399, 812

Wojtak, R., Hansen, S. H., & Hjorth, J. 2011, Nature, 477, 567

Wolf, J., Martinez, G. D., Bullock, J. S., et al. 2010, MNRAS, 406, 1220
## Appendix A: Equivalence of the Solanes & Salvador-Solé (1990) and Aguerri et al. (2017) anisotropy inversions

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

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

The kinematic projection equation (Binney & Mamon 1982)

$$P(R) = 2 \int_{R}^{\infty} \frac{r}{\sqrt{r^2 - R^2}} \left[ 1 - \beta(r) \left( \frac{R}{r} \right)^2 \right] p(r) dr \quad (A.1)$$

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

$$p(r) \equiv p_{\text{iso}}(r) = -\frac{1}{\pi} \int_{r}^{\infty} \frac{P'(R)}{\sqrt{R^2 - r^2}} dR = -\frac{1}{\pi r} \frac{dI}{dr}, \quad (A.2)$$

where (using the notation of SS90)

$$I(r) = \int_{r}^{\infty} \frac{R}{\sqrt{R^2 - r^2}} P(R) dR \quad (A.3)$$

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

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

$$[3 - 2\beta(r)] p(r) = \int_{r}^{\infty} \frac{q(s)}{s} ds + 2 p_{\text{iso}}(r). \quad (A.4)$$

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

$$\beta(r) p(r) = -\frac{1}{r^3} \int_0^r s^2 q(s) ds - p_{\text{iso}}(r) - \frac{3}{\pi r^2} I(r) \\ + \frac{3}{\pi r^3} \int_0^r I(R) dR. \quad (A.5)$$

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

$$\int_{0}^{r} I(s) ds = \int_{0}^{r} ds \int_{s}^{\infty} \frac{R}{\sqrt{R^2 - s^2}} P(R) dR \\ = \frac{\pi}{2} \int_{0}^{r} R P(R) dR + \int_{r}^{\infty} \sin^{-1}\left(\frac{r}{R}\right) R P(R) dR,$$

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

$$(\beta p)_{\text{new}}(r) = -\frac{1}{r^3} \int_0^r s^2 q(s) ds - p_{\text{iso}}(r) \\ + \frac{3}{2r^3} \int_0^r R P(R) dR \\ + \frac{3}{\pi r^3} \int_r^{\infty} \left[ \sin^{-1}\left(\frac{r}{R}\right) - \frac{r}{\sqrt{R^2 - r^2}} \right] R P(R) dR. \quad (A.6)$$

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

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

$$p_{\text{B+89}}(r) = \frac{1}{3} \int_{r}^{\infty} \frac{q(s)}{s} ds - \frac{2}{3r^3} \int_{0}^{r} s^2 q(s) ds \\ + \frac{1}{r^3} \int_{0}^{r} R P(R) dR \\ - \frac{2}{\pi r^3} \int_{r}^{\infty} \left[ \frac{r}{\sqrt{R^2 - r^2}} - \sin^{-1}\left(\frac{r}{R}\right) \right] R P(R) dR, \quad (B.1)$$

where we assumed that their $A = (3/2) \int_0^\infty R P(R) dR / \int_0^\infty qr^2 dr$ is unity, as expected from virial equilibrium (see the discussion by DM92)$^{12}$. With our notation, the radial pressure of DM92 (their Eq. (57a)) is

$$p_{\text{DM92}}(r) = \frac{1}{3} \int_{r}^{\infty} \frac{q(s)}{s} ds + \frac{2}{3r^3} \int_{r}^{\infty} s^2 q(s) ds \\ - \frac{2}{\pi r^3} \int_{r}^{\infty} \left[ \frac{r}{\sqrt{R^2 - r^2}} + \cos^{-1}\left(\frac{r}{R}\right) \right] R P(R) dR \quad (B.2)$$

where we assumed again the relation $A = 1$ (which translates to $G(\infty) = 0$ in DM92's notation). Replacing $\sin^{-1} x$ by $\pi/2 - \cos^{-1} x$, one can re-express the pressure of B+89 of Eq. (B.1) as

$$\begin{aligned} p_{\text{B+89}}(r) &= \frac{1}{3} \int_{r}^{\infty} \frac{q(s)}{s} ds - \frac{2}{3r^3} \int_{0}^{\infty} s^2 q(s) ds \\ &\quad + \frac{2}{3r^3} \int_{r}^{\infty} s^2 q(s) ds + \frac{1}{r^3} \int_{0}^{\infty} R P(R) dR \\ &\quad - \frac{2}{\pi r^3} \int_{r}^{\infty} \left[ \frac{r}{\sqrt{R^2 - r^2}} + \cos^{-1}\left(\frac{r}{R}\right) \right] R P(R) dR \\ &= p_{\text{DM92}}(r) \\ &\quad + \frac{1}{6\pi r^3} \left[ \int_{0}^{\infty} 6\pi R P(R) dR - \int_{0}^{\infty} 4\pi r^2 q(r) dr \right]. \end{aligned} \quad (B.3)$$

Now, the term in brackets in Eq. (B.3) is zero from the virial theorem, as the first integral represents the kinetic energy, while the second integral represents the absolute value of the potential energy (see Eq. (51) of DM92). This can be checked by

<footnote>$^{12}$ The elegant anisotropy inversion of DM92 suffers from several typos: in their Eqs. (43) and (44), $\psi$ should be replaced by $\Phi = -\psi$, and in their Eq. (45) $D_z$ should be $D_x$.</footnote>
inserting $P(R)$ from the equation of kinematic projection (A.1) into the first integral, yielding (after inversion of the order of the integrals)

$$ \int_{0}^{\infty} R P(R) \, dR = \frac{2}{3} \int_{0}^{\infty} [3 - 2\beta(r)] p(r) r^2 \, dr $$

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

$$ \frac{dp}{dr} + 2\beta \frac{p(r)}{r} = -\frac{q(r)}{r}, \qquad (B.4) $$

yielding (again after inversion of the order of integration)

$$ \int_{0}^{\infty} r^2 q(r) \, dr = \int_{0}^{\infty} [3 - 2\beta(r)] p(r) r^2 \, dr. $$

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

$$ \beta(r) = - \frac{r p'(r) + q(r)}{2 p(r)} = - \frac{1}{2} \left[ \frac{q(r)}{p(r)} + \frac{d \ln p}{d \ln r} \right]. \qquad (B.5) $$

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

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

$$ p'(r) = - \frac{q(r)}{r} - \frac{2}{\pi} \frac{d(J/r^3)}{dr} = - \frac{q(r)}{r} + \frac{6}{\pi r^4} J(r) - \frac{2}{\pi r^3} \frac{dJ}{dr}, \quad (B.6) $$

where

$$ J(r) = \int_{r}^{\infty} \left[ \cos^{-1} \left( \frac{r}{R} \right) + \frac{r}{\sqrt{R^2 - r^2}} \right] R P(R) \, dR. \qquad (B.7) $$

With the substitution $R = r/\cos\theta$, we can write

$$ J(r) = r^2 \int_{0}^{\pi/2} (\theta + \cot \theta) P\left(\frac{r}{\cos \theta}\right) \frac{\sin \theta}{\cos^3 \theta} d\theta, \qquad (B.8) $$

which leads to

$$ \frac{dJ}{dr} = 2r \int_{0}^{\pi/2} (\theta + \cot \theta) P\left(\frac{r}{\cos \theta}\right) \frac{\sin \theta}{\cos^3 \theta} d\theta + r^2 \int_{0}^{\pi/2} (\theta + \cot \theta) P'\left(\frac{r}{\cos \theta}\right) \frac{\sin \theta}{\cos^4 \theta} d\theta \qquad (B.9) $$

$$ = \frac{2}{r} \int_{r}^{\infty} \left[ \cos^{-1}\left(\frac{r}{R}\right) + \frac{r}{\sqrt{R^2 - r^2}} \right] R P(R) \, dR + \frac{1}{r} \int_{r}^{\infty} \left[ \cos^{-1}\left(\frac{r}{R}\right) + \frac{r}{\sqrt{R^2 - r^2}} \right] R^2 P'(R) \, dR. \qquad (B.10) $$

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

$$ r p'(r) + q(r) = \frac{2}{\pi r^3} \int_{r}^{\infty} \left[ \cos^{-1}\left(\frac{r}{R}\right) + \frac{r}{\sqrt{R^2 - r^2}} \right] \times R [P(R) - R P'(R)] \, dR. \qquad (B.11) $$

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

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

where

$$ A = \int_{r}^{\infty} \left[ \cos^{-1} \left( \frac{r}{R} \right) + \frac{r}{\sqrt{R^2 - r^2}} \right] R [P(R) - R P'(R)] \, dR, $$

$$ B = 2 \int_{r}^{\infty} \left[ \cos^{-1} \left( \frac{r}{R} \right) + \frac{r}{\sqrt{R^2 - r^2}} \right] R P(R) \, dR - \pi \frac{r^2}{3} \int_{r}^{\infty} q(s) \left( \frac{r}{s} + 2 \frac{s^2}{r^2} \right) \, ds. $$

<footnote n="13">13 Dejonghe & Merritt (1992) discuss at length the possibility of an extra $C/r^3$ term to the radial pressure, where the constant C is an integral over physical radius from zero to infinity, meaning that inward integration is required as well as the outward integrals of Eq. (B.2). On one hand, they argue that $C = 0$ to ensure that 1) both the radial and tangential velocity variances remain positive at large radii, and 2) the radial velocity dispersion does not reach unphysically large values in cases where the tracer density falls much faster than $1/r^3$. On the other hand, they argue that some choices for the gravitational potential (or equivalently the mass profile) may lead to situations where the virial theorem is inconsistent with the second moment equations, in which case the term $C/r^3$ should be incorporated into the solution of the radial pressure.</footnote>
Appendix C: Additional figure

<figure><img src="image_18.png" /><figcaption>Fig. C.1. Same as Fig. 3, but for model 11.</figcaption></figure>
# Appendix D: Two-body relaxation time for galaxies falling into clusters

## D.1. Transverse velocity formalism

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

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

$$|\delta v_{\perp}| = 2 \frac{Gm}{p v_{\text{rel}}}, \qquad (\text{D.1})$$

where $v_{\text{rel}}$ is the relative velocity of the two galaxies, $p$ is the impact parameter, and $G$ is Newton’s gravitational constant. During a time $dt$, the galaxy will suffer a number of encounters with other galaxies with impact parameter between $p$ and $p + dp$ of $d^2N = (v v_{\text{rel}} dt) \times (2\pi p dp)$. (D.2)

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

$$\Delta v_{\perp}^2 = \int_{\text{orbit}} dt \int_{p_{\text{min}}}^{p_{\text{max}}} \frac{d^2 N}{dp dt} \delta v_{\perp}^2 dp = 8\pi G^2 m^2 \int_{\text{orbit}} \frac{v(t)}{v_{\text{rel}}(t)} \ln \Lambda(t) dt, \qquad (\text{D.3})$$

where $\Lambda = p_{\text{max}}/p_{\text{min}}$. One generally assumes $p_{\text{min}} = Gm/v_{\text{rel}}^2$ and $p_{\text{max}} = r$, yielding $\Lambda = rv_{\text{rel}}^2/(Gm)$.

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

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

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

$$E = \frac{1}{2} v^2 + \Phi(r) = \frac{1}{2} v_{\text{peri}}^2 + \Phi(r_{\text{peri}}) = \frac{1}{2} v_{\text{apo}}^2 + \Phi(r_{\text{apo}}), \qquad (\text{D.4})$$

$$L^2 = r^2 v_{\perp}^2 = r_{\text{peri}}^2 v_{\text{peri}}^2 = r_{\text{apo}}^2 v_{\text{apo}}^2, \qquad (\text{D.5})$$

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

$$\left(\frac{dr}{dt}\right)^2 - \left[1 - \left(\frac{r_{\text{peri}}}{r}\right)^2\right] v_{\text{peri}}^2 - 2 \left[\Phi(r_{\text{peri}}) - \Phi(r)\right] = 0. \qquad (\text{D.6})$$

Eliminating $v_{\text{apo}}$ in the third equalities of Eqs. (D.4) and (D.5) yields

$$v_{\text{peri}}^2 = \frac{2}{1 - 1/\lambda^2} \left[ \Phi(\lambda r_{\text{peri}}) - \Phi(r_{\text{peri}}) \right], \qquad (\text{D.7})$$

where $\lambda = r_{\text{apo}}/r_{\text{peri}}$. Combining Eqs. (D.6) and (D.7) leads to an equation of motion expressed in terms of the pericenter and the ratio of apo- to pericenter:

$$\left(\frac{dr}{dt}\right)^2 - \frac{2}{1 - 1/\lambda^2} \left[1 - \left(\frac{r_{\text{peri}}}{r}\right)^2\right] \Phi(\lambda r_{\text{peri}}) - 2 \left\{1 - \frac{1}{1 - 1/\lambda^2} \left[1 - \left(\frac{r_{\text{peri}}}{r}\right)^2\right]\right\} \Phi(r_{\text{peri}}) + 2 \Phi(r) = 0. \qquad (\text{D.8})$$

## D.3. Change of transverse velocity over an orbit

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

$$\Delta v_{\perp}^2 = 16\pi G^2 m^2 \int_{r_{\text{peri}}}^{r_{\text{apo}}} \frac{v(r)}{v_{\text{rel}}(r)} \ln \left( \frac{r v_{\text{rel}}^2(r)}{G m} \right) \frac{dr}{|dr/dt|}, \qquad (\text{D.9})$$

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

$$\frac{\Delta v_{\perp}^2}{GM(a)/a} = 4N(a) \left(\frac{m}{M(a)}\right)^2 \int_{r_{\text{peri}}/a}^{r_{\text{apo}}/a} \frac{\tilde{v}(x)}{y_{\text{rel}}(x)} \ln \left( \frac{x y_{\text{rel}}^2}{m/M(a)} \right) \frac{dx}{|y_r|}, \qquad (\text{D.10})$$

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

$$\Phi(r) = -\frac{GM(a)}{a} \tilde{\Phi}\left(\frac{r}{a}\right) \qquad (\text{D.11})$$

$$\tilde{\Phi}(x) = \frac{1}{\ln 2 - 1/2} \frac{\ln(1 + x)}{x}. \qquad (\text{D.12})$$

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

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

$$y_{\text{rel}} = \sqrt{y_r^2 + \frac{2}{1 - 1/\lambda^2} \left(\frac{x_{\text{peri}}}{x}\right)^2 \left[\tilde{\Phi}(x_{\text{peri}}) - \tilde{\Phi}(\lambda x_{\text{peri}})\right]}, \qquad (\text{D.13})$$

where we used $x_{\text{peri}} = r_{\text{peri}}/a$, the second equality in Eq. (D.5) and Eq. (D.11).

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

Figure D.1 provides the log number of orbits to isotropize, $[GM(a)/a] / \Delta v_{\perp}^2$, whose inverse is given in Eq. (D.10), given our estimates of $m/M(a)$ and $N(a)$ for our cluster sample and its associated galaxies. The figure indicates that two-body relaxation is slow for the bulk of our galaxies, with $10^{1.5} \approx 30$ orbits required for deeply penetrating orbits with $\lambda = 5$ (the typical value of apo- to pericenter ratio found in $\Lambda$CDM halos,
<figure><img src="image_19.png" /><figcaption>Fig. D.1. Contours of inverse of $log_{10}[(\Delta v_{\perp}^2/[GM(a)/a])]$. This represents the decimal logarithm of the number of orbits for an infalling galaxy to isotropize by two-body relaxation, using Eq. (D.10) with Eqs. (8), (D.12), and (D.13), and assuming $m/M(a) = 0.004$ and $N(a) = 14$.</figcaption></figure>

Ghigna et al. 1998), and even more for orbits with greater pericenters. In comparison, a naïve application of the $N/(8 \ln N)$ formula (Eq. (4)–(9) of Binney & Tremaine 1987), would yield only 5.7 orbits at $r = a$ for $N = M(a)/m = 1/0.004 = 250$.

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

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

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

$$ \frac{1}{2} \frac{L^2}{r^2} + \Phi(r) - E = 0. \qquad (E.1) $$

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

The extraction of the roots of Eq. (E.1) for each of the $10^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 Eq. (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 Eq. (E.1) must be less than or equal to 0, since it represents $-\frac{1}{2}(dr/dt)^2$, we then vectorially adjusted $r_{peri}$ and $r_{apo}$ with the conditions

where $|LHS| < |oldLHS|$ & LHS < oldLHS
update $r_{peri}$

where $|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 (respectively 0.01) dex precision in this manner with a script language (SM, also known as SuperMongo) on a single processor. We found a median $r_{apo}/r_{peri}$ of 3.79 with an uncertainty of 0.02 (from ten trials).