{'fraction_non_alphanumeric': 0.03900687640722936, 'fraction_numerical': 0.020279316010466744, 'mean_word_length': 4.582420382165605, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
MAMPOSSt: Modelling Anisotropy and Mass Profiles of Observed Spherical Systems. I. Gaussian 3D velocities

Mass modelling of spherical systems through internal motions is hampered by the mass/velocity anisotropy (VA) degeneracy inherent in the Jeans equation, as well as the lack of techniques that are both fast and adaptable to realistic systems. A new fast method, called MAMPOSSt, which performs a maximum likelihood fit of the distribution of observed tracers in projected phase space, is developed and thoroughly tested. MAMPOSSt assumes a shape for the gravitational potential, but instead of postulating a shape for the distribution function in terms of energy and angular momentum, or supposing Gaussian line-of-sight velocity distributions, MAMPOSSt assumes a VA profile and a shape for the 3D velocity distribution, here Gaussian. MAMPOSSt requires no binning, differentiation, nor extrapolation of the observables. Tests on cluster-mass haloes from LambdaCDM cosmological simulations show that, with 500 tracers, MAMPOSSt is able to jointly recover the virial radius, tracer scale radius, dark matter scale radius and outer or constant VA with small bias (<10% on scale radii and<2% on the two other quantities) and inefficiencies of 10%, 27%, 48% and 20%, respectively. MAMPOSSt does not perform better when some parameters are frozen, and even worse when the virial radius is set to its true value, which appears to be the consequence of halo triaxiality. The accuracy of MAMPOSSt depends weakly on the adopted interloper removal scheme, including an efficient iterative Bayesian scheme that we introduce here, which can directly obtain the virial radius with as good precision as MAMPOSSt. Our tests show that MAMPOSSt with Gaussian 3D velocities is very competitive with, and up to 1000x faster than other methods. Hence, MAMPOSSt is a very powerful and rapid tool for the mass and anisotropy modeling of systems such as clusters and groups of galaxies, elliptical and dwarf spheroidal galaxies.

dust components, as well as the luminosity in different wavebands. These studies can be performed as a function of system mass and other attributes such as galaxy colour (e.g., Wojtak & Mamon 2013 and references therein).
Mass profiles can be derived from internal motions, or alternatively from X-ray or lensing observations. This paper focuses on mass profiles from internal kinematics. In this class of mass modeling, one has to deal with a degeneracy between the unknown radial profiles of total mass and of the velocity anisotropy (hereafter 'anisotropy') β(r) = 1 − σ 2 θ (r) + σ 2 φ (r) 2 σ 2 r (r) (1) (note that, in spherical symmetry, one must have σ φ = σ θ ). While radial outer orbits are expected for structures in an expanding universe (e.g. Ascasibar & Gottlöber 2008 for dark matter particles and Ludlow et al. 2009 for subhaloes), the dissipative nature of the gas dynamics is expected to produce tangential orbits in the inner regions of systems formed from gas-rich mergers or collapse. Therefore, lifting the Mass -Anisotropy Degeneracy can provide useful constraints on the formation of the structure under study.
A common method to extract the mass profile is to assume that the line-of-sight (hereafter, LOS) velocity distribution, at given projected radius, is Gaussian (Strigari et al. 2008;Battaglia et al. 2008;Wolf et al. 2010). These methods perform adequately on the mass profile, but provide weak constraints on the anisotropy (Walker et al. 2009). Merritt (1987) pointed out that anisotropic models have non-Gaussian LOS velocity distributions. Therefore, the observed kurtosis of the distribution of LOS velocities serves as a powerful constraint to the anisotropy (Gerhard 1993;van der Marel & Franx 1993;Zabludoff, Franx, & Geller 1993). In fact, if one assumes that the anisotropy is constant throughout the system, the fourth order Jeans equation can be used to express the LOS velocity kurtosis as an integral of the tracer density, anisotropy and total mass profiles (Łokas 2002). Moreover, Richardson & Fairbairn (2012) were able to generalize the expression for the LOS velocity kurtosis for radially varying anisotropy in the framework of separable augmented density and 4th order anisotropy equal to the standard anisotropy (the latter appears to be an excellent approximation for ΛCDM haloes, see Fig. 10 of Wojtak et al. 2008). One can then perform a joint fit of the observed LOS velocity dispersion and kurtosis profiles. This was found to (partially) lift the mass-anisotropy degeneracy when applied to dwarf spheroidal galaxies (Łokas 2002) and the Coma cluster, (Łokas & Mamon 2003): the joint constraint of LOS velocity dispersion and kurtosis profiles allows the estimation of both the mass profile (i.e., normalization and concentration) and the anisotropy of the cluster, contrary to the case when the LOS kurtosis profile is ignored.
An interesting route is to perform non-parametric inversions of the data assuming either the mass profile to obtain the anisotropy profile (anisotropy inversion, pioneered by Binney & Mamon 1982) or the anisotropy profile to obtain the mass profile (mass inversion, independently developed by Boué 2010 andWolf et al. 2010). These inversion methods are powerful in that they are nonparametric, but they suffer from their requiring the user to bin the data, smooth it, and extrapolate it beyond the range of data.
Hence, one would like to go one step further and constrain the full information contained in the observed projected phase space (projected radii and LOS velocities, hereafter, PPS) of LOS velocities as a function of projected radii. In other words, rather than using the 0th, 2nd and possibly 4th moments of the LOS velocity distribution, we wish to use the full set of even moments.
The traditional way to analyze the distribution of particles in PPS is to assume a form for the six-dimensional distribution function (DF) in terms of energy (E) and angular momentum (J) and fit the triple integral of equation (5) below, using this DF for f , to the distribution of particles in PPS. The worry is that we have no good a priori knowledge of the shape of the DF, f (E, J). One clever idea is to throw orbits in a gravitational potential, since each orbit is a Dirac delta function in energy and angular momentum. One then seeks a linear combination of these orbits, with positive coefficients, to match the data. This orbit model (Schwarzschild 1979;Richstone & Tremaine 1984;Syer & Tremaine 1996) is very powerful (and can handle non-spherical gravitational potentials), but too slow to obtain meaningful errors on the parameters. A similar, and in principle faster, technique is to assume that the DF is the linear combination (again with positive coefficients) of elementary DFs (Dejonghe 1989;Merritt & Saha 1993;Gerhard et al. 1998), but only one such study has been made (Kronawitter et al. 2000), and it is not clear whether the elementary DFs, although numerous, constitute a basis set.
An important step forward has been performed by Wojtak et al. (2008), who analyzed the haloes in ΛCDM cosmological simulations, to show that the DF can be approximated to be separable in energy and angular momentum, with a simple analytical approximation for the angular momentum term. In a sequel, Wojtak et al. (2009) have shown that it is feasible to fit the distribution of particles in PPS with equation (5), using the approximation of the DF found by Wojtak et al. (2008).
However, it is not yet clear whether self-gravitating quasispherical astrophysical systems have the DF of ΛCDM haloes: In particular, if the dynamical evolution of these systems is influenced by the dissipation of their gaseous component, the DF may not be separable in terms of energy and angular momentum. Dissipation is not expected to affect much the internal kinematics of large systems such as galaxy clusters 1 , but is expected to be increasingly important in smaller systems such as galaxy groups, and especially galaxies themselves. For this reason, it is useful to consider a massmodeling method that is independent of the dependence of the DF on energy and angular momentum.
In this work, we present an alternative method, in which we fit the distribution of particles in PPS making assumptions on the radial profiles of mass and anisotropy as well as the radial variations of the distribution of space-velocities. We call this method Modelling of Anisotropy and Mass Profiles of Observed Spherical Systems, or MAMPOSSt for short. 2 The MAMPOSSt method is described in Sect. 2.1, its Gaussian approximation is described in Sect. 2.2. Tests on haloes derived from a cosmological N body simulation are presented in Sect. 3. A discussion follows in Sect. 4.

General method
The observed tracer population of a spherical system has a DF f (r, v) = ν(r) fv(v|r) , where ν(r) is the tracer number density profile. MAMPOSSt fits the distribution of objects in PPS (projected radius R and LOS velocity vz), assuming parametrized forms for (i) the gravitational potential (or equivalently a total mass density profile, through the Poisson equation), (ii) the anisotropy profile (eq. [1]) (iii) the distribution of 3D velocities, fv(v|r).
Consider a point P at distance r from the centre, O, of the spherical system, with projected radius R r and consider the spherical coordinates where the unit vectors er and e θ are in the plane containing OP and the LOS, while e φ is perpendicular to this plane. Consider also the cylindrical coordinate system (vz, v ⊥ , v φ ), where ez is the axis along the LOS and e ⊥ is the axis perpendicular to the LOS, but in the plane containing O and P and the LOS. The Jacobian of the transformation from the spherical coordinate system to the new one is unity, hence one can write The distribution of LOS velocities at P is then obtained by integrating velocities over the two perpendicular axes (dropping {r, R} from fv for clarity): Note that dynamical systems have maximum velocities set by the escape velocity, −2 Φ(r) (where Φ(r) is the gravitational potential), on one hand, and by the maximum allowed (observable) absolute LOS velocity on the other hand. In what follows, we will neglect both limits, unless explicitly mentioned otherwise. The surface density of observed objects (the tracer) in PPS is then obtained by integrating along the LOS is the tracer surface density at projected radius R. Equation (5) is equivalent to equation (2) of Dejonghe & Merritt (1992). If the tracer number density profile ν(r), appearing in equation (4), is not known and if the incompleteness of the data is independent of the projected radius, then one can estimate ν(r) by Abel inversion of Σ(R) of equation (6): But this is not necessary, as we shall see below.
In MAMPOSSt, rather than replace the velocities by energy and angular momentum and numerically solve the triple integral of equation (5) Wojtak et al. 2009), we analytically derive h(vz|R, r) from equation (3) for known 3D velocity distributions. With the analytical form of h(vz|R, r), equation (4) provides the surface density distribution of tracers in PPS through a single integral. Note, however, that another single integral is required because the expression for h(vz|R, r) will involve σr(r) (see eqs. [25] and [26], below, for the Gaussian case), which is obtained by solving the spherical Jeans equation where β is the anisotropy of (eq. [1]) for our given choices of total mass and anisotropy profiles. We thus need to insert the solution (van der Marel 1994; in the expression for h(vz|R, r) (eq. [3]) to derive g(R, vz), via equation (4), where β(t) is given, while ν(r) is obtained with equation (7). In equation (9), M (s) = (s 2 /G) dΦ/ds is the radial profile of the total mass (this is the only instance where the gravitational potential enters MAMPOSSt). For a given choice of parameters, the single integral of equation (4) must be evaluated for every data point (R, vz), whereas the other integral (eq. [9]) for σr(r) need only be evaluated once, on an adequate grid of r.
Note that for projected radii extending from Rmin to Rmax and absolute LOS velocities extending from 0 to a maximum velocity, which for projected radius R is theoretically equal to vesc(R) = −2Φ(R), and in practice is possibly specified by a cut of obvious velocity interlopers, vcut(R), one can write where we used equation (4) for g(R, vz), assumed that h(vz|R, r) is normalised, reversed the order of the integrals in r and vz, and where Np(R) is the predicted number of objects within projected radius R, while ∆Np = Np(Rmax) − Np(Rmin). Equation (10) then implies that the probability density of observing an object at position (R, vz) of PPS is where equation (11) arises from equation (4), while equation (12) is obtained by writing r = R cosh u. Here, Np(R) is the number of tracers in a cylinder of projected radius R, the terms ν and Np are given by where N (r) is the cumulative tracer number density profile, while rν is the characteristic radius of the tracer. One easily verifies that q(R, vz) dR dvz = 1. The values of Rmin and Rmax appearing in ∆Np (eq. [10]) can be hard limits, or alternatively the respective minimum and maximum projected radii of the observed tracers if no hard limits are specified. We fit the parameters (mass scale or concentration and possibly normalization, anisotropy level or radius, as well as the tracer scale rν -if not previously known) that enter the determination of g(R, vz) to the observed surface density, using maximum likelihood estimation (MLE), i.e. by minimizing for the N -parameter vector θ, where n is the number of data points, with q given by equation (11). Writing θ = {rν , η}, where η is the vector of the N − 1 parameters other than rν, one has where and Combining the last equality of equation (18) with equation (5), integrating over LOS velocities, reversing the order of the two outer integrals of the resulting quadruple integral, and using equation (6) yields p(vz|R) dvz = 1. So, if the scale of the tracer distribution is already known, then, according to equations (15) and (16), maximizing the likelihood amounts to minimizing Now, if rν is not known, then one may be tempted to solve for it by minimizing − ln L0 = − ln p0(R|rν), and then proceed with equation (19) to minimize for the N − 1 remaining parameters, η. However, since − ln L = − ln L ′ − ln L0 (from eqs. [15] and [16]), the most likely solution for θ that minimizes − ln L will not in general be that which minimizes at the same time − ln L ′ and − ln L0. Moreover, if one seeks to obtain the distributions of parameters η and rν consistent with the MLE solution (for example with Markov Chain Monte-Carlo techniques), the joint analysis of equations (11) and (15) is required. On the other hand, if rν is known from other data, while the current dataset is known to have a completeness, C(R), that is a function of projected radius, then one could indeed minimize ln L ′ of equation (19). The proper solution is then to minimize − ln L weighting the data points by the inverse completeness, i.e. minimizing For computational efficiency, we perform the following tasks: (i) For each run of parameters, we first compute log σr(rj) from equation (9) on a logarithmic grid of rj, and compute cubic-spline coefficients at these radii. Then, when we compute the LOS integral of equation (4) for each (Ri, vz,i), we evaluate σr(r) with cubic spline interpolation (in log-log space, using the cubic spline coefficients determined at the start).
(ii) For simple anisotropy models, the exponential term in equation (9) is given by equations (A2) and (A3).
(iii) We terminate the LOS integration in equation (11) at roughly 15 virial radii, 3 rv, instead of infinity, as the Hubble flow pushes the velocities of the material beyond this distance to values over 3 σv above the mean of the system (see Mamon, Biviano, & Murante 2010, hereafter MBM10). The LOS integration varies only very slightly with the number of virial radii, so as long as the virial radius is correct to a factor of two, this choice of integration limit is not an issue.
We now need to choose a model for the shape of the 3D velocity distribution. While MAMPOSSt, can, in principle, be run with any model, the simplest one is the (possibly anisotropic) Gaussian distribution, which we describe in Sect. 2.2 below.

Gaussian 3D velocity distributions
The simplest assumption for the 3D velocity distribution is that it is Gaussian: where the velocity dispersions σi are functions of r. This Gaussian distribution assumes no streaming motions: e.g. no rotation, and no mean radial streaming, which is adequate for Rmax < rv in high-mass haloes (i.e. groups and clusters) and Rmax < 4 rv in galaxy-mass haloes (Cuesta et al. 2008). Inserting equation (21) into equation (3) and integrating over v φ leads to Calling θ the angle between the line-of-sight (direction z) and the radial vector r, one has with which the integral over v ⊥ in equation (22) yields a Gaussian distribution of LOS velocities at point P: of squared dispersion The integral of h(vz|R, r) along the LOS is obtained from equations (4) and (25): According to equations (18) and (27), the probability of measuring a velocity vz at given projected radius R is We remind the reader that β is a chosen function of r, ν is a function of r given by equation (7), while σr is a function of r given by equation (9). For isotropic systems (β = 0), equation (27) leads to a Gaussian distribution of LOS velocities. However, for anisotropic velocity tensors, the distribution of LOS velocities will generally not be Gaussian (as Merritt 1987 found when starting from distribution functions instead of Gaussian 3D velocities). Hence, the Gaussian nature of h(vz|R, r) is not equivalent to the popular assumption that g(R, vz) is Gaussian on vz: even if h(vz|R, r) is a Gaussian at point P, its integral along the LOS is not Gaussian, unless β = 0 and σr is constant.
If one of the parameters to determine with MAMPOSSt is the normalization of the mass profile, one should not be tempted in expressing the radii in terms of the virial radius rv, the velocities in terms of the virial velocity vv, the tracer densities in terms of what we wish (as they appear in both the numerator and denominator of eq. [28]). Doing so, equation (28) becomes where the quantities with tildes are in virial units. Equation (29) indicates that when one varies the rv (and the virial velocity in proportion as vv = 2/∆ H0 rv), the highest probabilities are reached for the highest normalizations: vz becomes very small, while σz is unaffected to first order. This unphysical result is the consequence of using a parameter (the virial radius) as part of the data variable. On the other hand, using equation (28), one sees that the highest probabilities p(vz|R) are reached at intermediate values of the normalization.

Simulated haloes
To test MAMPOSSt, we use cluster-mass haloes extracted by Borgani et al. (2004) from their large cosmological hydrodynamical simulation performed using the parallel Tree+SPH GADGET-2 code of Springel et al. (2005). The simulation assumes a cosmological model with Ω0 = 0.3, ΩΛ = 0.7, Ω b = 0.039, h = 0.7, and σ8 = 0.8. The box size is L = 192 h −1 Mpc. The simulation used 480 3 DM particles and (initially) as many gas particles, for a DM particle mass of 4.62 × 10 9 h −1 M⊙. The softening length was set to 22.5 h −1 comoving kpc until z = 2 and fixed afterwards (i.e., 7.5 h −1 kpc). The simulation code includes explicit energy and entropy conservation, radiative cooling, a uniform time-dependent UV background (Haardt & Madau 1996), the self-regulated hybrid multi-phase model for star formation (Springel & Hernquist 2003), and a phenomenological model for galactic winds powered by Type-II supernovae.
DM haloes were identified by Borgani et al. (2004) at redshift z = 0 with a standard Friends-of-friends (FoF) analysis applied to the DM particle set, with linking length 0.15 times the mean inter-particle distance. After the FoF identification, the centre of the halo was set to the position of its most bound particle. A spherical overdensity criterion was then applied to determine, for each halo, our proxy for the virial radius, r200, where the mean density is 200 times the critical density of the Universe.
To save computing time, we worked on a random subsample of roughly two million particles among the 480 3 . We have extracted 11 cluster-mass haloes from these simulations, among which, ten are about logarithmically spaced in virial radius, r200, while the 11 th halo is the most massive in the entire simulation. Their properties are listed in Table 1. We made no effort to omit irregular haloes, but among the list of 12 irregular haloes out of 105 extracted by MBM10 from the same simulation, 2 are in our sample (haloes 17283 and 434). We list the characteristic radii rs, rH, rB of three models fitted by MLE to the mass density profiles of the particle data (from 0.03 to 1 r200), namely: where rs ≡ r−2 is the radius of slope −2 in the mass density profile, related to the concentration c ≡ r200/r−2; (ii) the Hernquist density profile (Hernquist 1990) where rH = 2 r−2, (iii) the Burkert density profile (Burkert 1995) where rB ≃ 0.657 r−2.
Denoting the scales rs, rH and rB by the generic rρ, the mass profiles of these models (required for eq. [9]) can be written where The NFW model has long been known to fit well the density profiles of ΛCDM haloes (Navarro et al. 1996), and while Navarro et al. (2004) found that Einasto models fit them even better, MBM10 found that the NFW model describes the outer LOS velocity dispersion profile of the DM component of their stacked cluster-mass halo in Borgani et al.'s hydrodynamical cosmological simulation even (slightly) better than the Einasto model. The Hernquist model differs from the NFW one because it has a steeper logarithmic slope at large radii, γ ≡ d ln ρ/d ln r = −4 rather than −3. The Burkert model, on the other hand, has the same asymptotic γ = −3 as the NFW model, but a core at the centre, γ = 0, rather than a cusp (γ = −1 in both the NFW and Hernquist models).
In Table 1, we also list the values of the parameters characterizing different velocity-anisotropy models, namely: where we assume spherical symmetry and therefore σ θ = σ φ ; (ii) the model ('ML' model hereafter) of ; characterized by the anisotropy radius r β ; (iii) a generalization of the ML model, which is also a simplified version of the model of Tiret et al. (2007), isotropic at r = 0 and with anisotropy radius identical to r−2 (hereafter called 'T' model): characterized by the anisotropy value at large radii, β∞. In our T model, the anisotropy radius is set to the scale radius of the mass density profile. Note also that in the following we provide the values of A∞ ≡ (σr/σ θ )∞ = (1 − β∞) −1/2 , rather than β∞. Note that the ML and the T models used here are identical for β∞ = 0.5 and r β = r−2. With these values, the ML and T models provide a good fit to the average anisotropy profile of a set of cluster-mass cosmological haloes (MBM10). In Fig. 1, we show the individual halo velocity anisotropy profiles and the ML anisotropy model with r β = r−2 (or, equivalently, the T anisotropy model with β∞ = 0.5). There is a huge scatter in the β(r) of the individual haloes, as already observed by, e.g., Wojtak et al. (2008), especially at r > 0.3 r200, while 8 out of 11 haloes have β(0) = 0 ± 0.15.

Observing cones and interloper removal
To test MAMPOSSt, we select 500 DM particles around each halo, out to a maximum projected distance Rmax from the halo centre, for which we consider three values: r500 ≃ 0.66 r200, r200, and r100 ≃ 1.35 r200. We analyze three orthogonal projections for each halo -these are in fact cones with an observer at D = 90 h −1 Mpc away, but the opening angle being very small has no noticeable effect on our results. The particles in these cones are used by MAM-POSSt as tracers of the halo gravitational potential.
However, these 500-particle samples include interlopers, i.e. DM particles that are located in projection at R Rmax, but are effectively outside Rmax in real (3D) space, i.e. with r > Rmax. It is impossible to remove all these interlopers in the observed redshift space, where only 3 of the 6 phase-space coordinates of the tracers are known (e.g., MBM10). Moreover, since the LOS velocity distribution of interlopers in mock cones around ΛCDM haloes is the sum of a Gaussian component and a uniform one (see MBM10 for a quantified view), and since the Gaussian one resembles that of the particles in the virial sphere, it is important to remove the flat LOS velocity component, at least at high absolute LOS velocity, where it dominates. It is possible to remove these high |vz| objects with suitable interloper removal algorithms.
To see how MAMPOSSt depends on the choice of the interloper removal algorithm, we here consider three different algorithms.
The first one is a new, iterative algorithm, that we name "Clean", which is fully described in Appendix B. Clean first looks for gaps in the LOS velocities, then estimates the virial radius r200 from the aperture velocity dispersion, assuming an NFW model with ML anisotropy with an anisotropy radius r β = r−2 and a concentration depending on the estimate of r200 via the relation of Macciò, Dutton, & van den Bosch (2008), then only considers the galaxies within 2.7 σz(R) from the median LOS velocity, and finally iterates. Our assumed anisotropy profile fits reasonably well the anisotropy profiles of DM haloes (MBM10), as is clear for our 11 haloes (see Fig. 1). The factor 2.7 was found by MBM10 to best preserve the local LOS velocity dispersion for the assumed density and anisotropy models.
We also consider two other interloper removal algorithms, namely: (i) the method (hereafter, dHK) of den Hartog & Katgert (1996), a widely used procedure that works reasonably well on clustermass haloes from cosmological simulations (Biviano et al. 2006;Wojtak et al. 2007), despite its crude underlying physics; (ii) the method (hereafter, KBM) of Katgert, Biviano, & Mazure (2004, see their Appendix A), in which a galaxy is flagged as an interloper under the condition vz/σz > 1.85 (R/r200) −0.3 , with r200 derived from σz using eq. (8) of Carlberg, Yee, & Ellingson (1997). This method was invented as a poor-man's proxy for the dHK method when the observational sampling of the halo projected phase-space is poor.

The general 4-parameter case
There is no a priori limitation on the number of free parameters that can be used in MAMPOSSt to characterise the mass and velocity anisotropy profiles. With samples of 500 tracers (assumed massless throughout these tests) it is appropriate to consider ∼ 4 free parameters, two for the mass distribution, one for the velocity anisotropy distribution, and one for the spatial distribution of the tracers. All these models are characterised by the two free parameters, the 'virial' radius r200 and a characteristic scale-radius (rs, rH, and rB for the NFW, Hernquist, and Burkert models, respectively). Herafter, we generically use rρ to refer to this characteristic scaleradius of the mass density profile.
We use the NFW model, in projection (Bartelmann 1996;Łokas & Mamon 2001), to fit the projected number density profile of the tracer. Note that the normalization of this profile does not enter the MAMPOSSt equations, so the only free parameter is r−2. Herafter we call this parameter rν , to avoid confusion with the characteristic radius of the NFW mass density profile. We only consider one model for the number density profile of the tracer, because this is a direct observable, unlike the mass density profile. While one should not be too restrictive in the model choice for the mass density profile, the observer is generally able to choose the best-suited model for the tracer number density profile by direct examination of the data before running MAMPOSSt. We choose the NFW model because it provides a reasonable description of the number density profiles of the DM particles in our simulated haloes.
For the velocity anisotropy profile, we consider the three models described above, Cst, ML, and T, each characterised by a single anisotropy parameter, A, r β , and A∞, respectively. In equation (38), we use r−2 = rρ.
To search for the best-fit solution, we run the MAMPOSSt algorithm in combination with the NEWUOA 4 minimization routine of Powell (2006). For estimating error bars on the best fit parameters, as well as confidence contours on pairs of parameters, we fit our model parameters using the Markov Chain Monte Carlo (MCMC) technique (e.g., Lewis & Bridle 2002). In MCMC, the kdimensional parameter space is populated with proposals, for each of which the likelihood is computed. The new proposal is accepted if the ratio of new to previous likelihood is either greater than unity or else greater than a uniform [0, 1] random number. The proposal is found by assuming a k-dimensional Gaussian probability distribution around the previous proposal. We adopt the publicly available CosmoMC code by A. Lewis. 5 We run 6 chains in parallel using Message Parsing Interface (MPI), and the covariance matrix is used to update the parameters of the Gaussian proposal density to ensure faster convergence. Fig. 2 illustrates the MAMPOSSt analysis via MCMC for the general case with four free parameters, using the NFW model for the mass density profile, and the constant (free parameter) anisotropy model. In particular, it shows that the different parameters are not correlated, except for a positive correlation between rρ and A.
Our FORTRAN code takes roughly 1 second to find the MAMPOSSt 4-parameter solution for a 500 particle sample run in scalar on a decent desktop or laptop computer, and 4 minutes to produce confidence limits for this solution with the CosmoMC (Lewis & Bridle 2002) MCMC code, with 6 chains of 40 000 elements run in parallel (MPI) on a PC equipped with a 4-core 8-thread Intel Core-I7 2600 processor.
The results for the different interloper rejection methods, mass density and velocity anisotropy models, and for the different maximum projected radii used in the selection of the 500 particles are listed in Table 2 and displayed in Fig. 3. We list and show the biweight measures (see, e.g., Beers, Flynn, & Gebhardt 1990) of Figure 2. Illustration of MAMPOSSt analysis for the general case, with β independent of radius, for a 500 particle sample from axis x of halo 25174 (grey broken line in Fig. 1), using MCMC (with 6 chains of 40 000 elements). The contours are 1, 2, and 3 σ. The red arrows and stars indicate the maximum likelihood solution, while the green arrows and crosses show the true solution (Table 1). The priors for the MCMC were uniform within the boxes of each panel and zero beyond the boxes. mean and dispersion of log(o/t) where o is the recovered value of the parameter and t its true value, because, according to our tests, they perform better than standard statistical estimators of location and scale when the parent distributions are not pure Gaussians. We call 'bias' and 'inefficiency' the mean and dispersion of log(o/t). If the dispersion in true values of a given parameter is small, one can spuriously obtain low values of the log(o/t) dispersion when the MAMPOSSt and true values show no correlation. We therefore also list the Spearman rank correlation coefficient between o and t, marking in boldface those correlations that are significant at the 99% confidence level. We list the results for all interloper rejection methods and all anisotropy models only for the NFW mass density model and for the R r200 radial selection. For simplicity, we only show a limited set of results for the other mass density models and for the other radial selections.
Remarkably, as seen in Table 2, the results for the four parameters are almost independent of the interloper removal algorithm, the Clean and KBM algorithm performing slightly better than dHK. The results for r200, rν, and rρ also depend very little on the cho-sen anisotropy model. On average, the values of the r200 parameter are recovered with almost no bias (from −1 to +4%) and with only ∼ 10% inefficiency. The rν parameter estimates are always slightly positively biased (4-7%), and are recovered with ∼ 25% efficiency. Also, the rρ parameter estimates generally display a slight positive bias, except for the KBM interloper removal method, and overall the bias ranges from −2 to +15%, while the efficiency ranges from ∼ 50 to ∼ 90%.
As far as the anisotropy parameter is concerned, the ML model behaves very differently from the Cst and T models, in that it is virtually impossible to constrain the anisotropy radius of the former, r β , while it is possible to obtain quite good constraints on the anisotropy parameters of the other two models, A and A∞. More precisely, the r β estimates are always negatively biased (by ∼ 60%) and are affected by a huge dispersion (almost one order of magnitude). On the other hand, A and A∞ are almost unbiased (the bias ranges from −10 to +5%) and they are affected by dispersions of, typically, ∼ 20%, if we consider the Clean and KBM interloper removal algorithms. So, apparently, it is much easier to constrain  the 'normalisation' of a given anisotropy profile, than to constrain the characteristic radius at which the anisotropy changes, particularly so if this change is mild, as in the ML model. Note, however, that the difficulty of MAMPOSSt in constraining the anisotropy parameter of the ML model does not mean that the ML model is a poor representation of reality, and in fact Fig. 1 suggests the opposite. Moreover, constraints obtained on the r200, rν, and rρ parameters are equally good with the Cst and T anisotropy models, as with the ML one. As seen in Table 2, correlations between recovered and observed values of the parameters r200, rν, and rρ are almost always signficant. This is also true for the A∞ anisotropy parameter, but not for A and r β . In Fig. 4, we show the correlations existing between the true and recovered values of the different parameters, using the T anisotropy model, for the 11 haloes along the 3 different orthogonal projections. Projections effects render the determination of the mass and anisotropy profile of a single 500-particle halo very uncertain. However, Fig. 4 shows that ∼ 500 tracers are sufficient to rank haloes for each of the different parameters considered here, i.e. by mass (r200), scale radius of the tracer distribution (rν) and of the total mass distribution (rρ), and outer velocity anisotropy A∞.
The importance of projection effects is also very clear from Fig. 5, where we display the ratio of the recovered to true values of the parameters for each halo along the three different projection axes. This figure also shows that there is no trend of under-or overpredicting the parameter values with halo mass.
All the above considerations apply for the NFW mass profile. Our tests with the Hernquist and Burkert mass profiles give similar results, as can be seen in Table 2, where for simplicity, we only list the results for the Clean interloper-removal algorithm and for the Tiret anisotropy model. The results are very similar for the NFW and Hernquist models. Results are similar also for the Burkert model, except for the scale rν of the number density profile, for which the bias and inefficiency are both higher than those obtained using the NFW and Hernquist mass models. Since the model we use for the number density profile has not changed (a projected NFW), this result suggests that it is difficult to accomodate a tracer with a central cuspy spatial distribution in a potential with a central core.
All the results described so far were obtained for a selection of 500 particles within Rmax = r200. Changing the value of Rmax is not without effects on the results. The inefficiency on r200 decreases when Rmax is gradually increased from r500 to r100. The inefficiency on anisotropy is largest for the smallest Rmax, and statistically similar for the two larger values. Increasing the aperture to the virial radius or above increases the number of tracers near the virial radius where r200 is estimated. Moreover, increasingly larger apertures will capture better the asymptotic value of the anisotropy profile A(r). In contrast, rρ is less efficiently determined when Rmax is increased from r500 to r100, while rν has its worst inefficiency for Rmax = r100, with statistically equivalent values for the two smaller apertures. This might be due to the increasing fraction of unidentified interlopers, and/or to the presence of (sub)structures at larger radii.
To assess the sensitivity of the MAMPOSSt technique to the number of tracers, besides the 500-particle selection, we have also considered samples of 100 particle tracers, randomly extracted from the same projections of the same 11 cosmological haloes. Results of the MAMPOSSt analysis are listed at the bottom of Table 2 and displayed in Fig. 3. Also in this case, for simplicity, we only . The mean (dots) and dispersion (error bars) of log(o/t) are respectively illustrated as filled circles and error bars, for the 33 samples of 500 ('Clean', dHK and KBM) and 100 ('Clean' (N =100)) particles. Results for the anisotropy models Cst, T, and ML are shown left to right in green, blue, and red, respectively. display a limited set of results. When compared to the results for the 500-particle samples, there is no significant change in the average values of the bias with which the different parameters are recovered, except for rρ, where the bias becomes negative, while it was mostly positive for the 500-particle samples. The rρ parameter value underestimation is not very severe, however, 25%. The efficiencies with which the different parameters are estimated are significantly affected by the reduction in number of tracers. The dispersion increases from ∼ 10 to 15% for r200, from ∼ 25 to 60% for rν, from ∼ 60 to 100% for rρ, and from ∼ 20 to 33% for A and A∞. There is no significant change in the dispersion for r β , but this was already extremely large for the 500-particle samples.

Cases with constraints on parameters
What is the effect of reducing the number of free parameters on the performance of the MAMPOSSt algorithm? To assess this point we consider several cases that reproduce what observers do in practice when faced with the problem of determining the internal dynamics of cosmological haloes. In all cases, we consider 500 particles selected within r200 in each halo. We only apply the Clean interloper removal algorithm, and we only consider the NFW mass density model, for simplicity.
Specifically, we consider the following Cases: A) General [Gen]: r200, rν, rρ, and the anisotropy parameter (one among the following: A, r β , A∞, depending on the anisotropy model considered) are all free MAMPOSSt parameters. This is the case considered so far. B) General with rν fitted outside MAMPOSSt [Split]: the free parameters are the same as in the Gen case, but rν is fitted outside MAMPOSSt, via MLE. We thus split the minimization of the parameters into two parts. C) Known virial mass or radius [KVir]: r200 is fixed and assumed to be exactly known, rρ and the anisotropy parameter are free parameters in MAMPOSSt, rν is an external free parameter, as in the Split case. D) Estimated virial mass or radius [EVir]: similar to KVir, except that r200 is not the true value, but the value estimated from the LOS aperture velocity dispersion (after interloper removal, see Appendix B). E) ΛCDM: rρ is estimated from r200 using the theoretical relation between these two quantities provided by Macciò et al. (2008); r200 and the anisotropy parameter are free parameters in MAM-POSSt, rν is an external free parameter, as in the Split case.
F) Mass follows Light [MfL]: r200 and the anisotropy parameter are free parameters in MAMPOSSt, rν is an external free parameter, as in the Split case, and rρ is assumed to be identical to rν. G) Tied Light and Mass [TLM]: r200 and the anisotropy parameter are free parameters in MAMPOSSt, while rρ and rν are assumed to be an identical free parameter.
H) Isotropic [β-iso]: A = 1 is assumed, r200, rρ are free parameters in MAMPOSSt, rν is an external free parameter, as in the Split case.   The Cases outlined above correspond to different observational situations. Case Gen is the most general situation in which the observer ignores all the dynamical characteristics of the system, and it is the case analysed so far (see Table 2). We repeat it here for the sake of comparison with the other cases.
Case Split is closer than case Gen to the typical observational situation, since the number density profile of the tracers of the potential (stars, galaxies) is generally determined directly (e.g. by counting them in concentric radial annuli, if the system is assumed to be spherical). The advantage of Split relative to Gen is that any radial incompleteness that might affect the determination of rν can be easily corrected for, if known, before running MAM-POSSt (see, e.g. Biviano & Poggianti 2009). However, as outlined in Sect 2.1, an interesting alternative is to include this incompleteness as in equation (20). Here, in all other cases (except TLM), rν has been determined outside MAMPOSSt. Anyway, the results (see below) obtained for the Gen and Split cases are very similar, hence it makes little difference if rν is fitted within or outside MAM-POSSt, when the analysed sample (as in our case) does not suffer from radial-dependent incompleteness.
In the KVir and EVir cases, r200 is not a free parameter, but is fixed externally. In the KVir case, a perfect knowledge of r200 is assumed, and in fact we take the true r200 of the cosmological haloes. This case corresponds to the situation in which r200 estimates are available from other data than the projected-phase space distribution of galaxies, e.g. from weak-lensing or X-ray observations for clusters, although in the real world also these mass estimates will be affected by uncertainties, so KVir is a rather idealised case. Case EVir corresponds to the situation in which r200 estimates are directly obtained from the Clean interloper removal scheme, as described in Appendix B.
In the ΛCDM and MfL cases, rρ is not a free parameter. In the ΛCDM case, we determine rρ from r200 using the relation of Macciò et al. (2008), which is based on the analysis of haloes in ΛCDM numerical simulations. Case MfL corresponds to the situation in which the observer has good reasons (or a priori theoretical prejudice) to assume that 'Mass follows Light', i.e. that the tracer is spatially distributed like the mass. Therefore the observer first determines rν from the distribution of the tracer and then makes the assumption rρ = rν. In case TLM, the observer is unable to constrain the tracer scale radius from its spatial distribution. This may happen when dealing with an incomplete sample with unknown incompleteness. Hence, rν and rρ are both determined from the dynamical analysis, assuming they are identical.
In the β-iso, β-MBM, and β-HM cases, the anisotropy is no longer a free parameter. The β-iso case corresponds to the situation where, for lack of better knowledge, the velocity anisotropy profile is assumed to be isotropic, A = 1 (β = 0) (like in, e.g. Biviano & Girardi 2003). Since the velocity anisotropy profiles of cluster-mass haloes have been shown to be well represented by an ML profile with r β = rρ (MBM10, see also Fig. 1), it also makes sense to fix this anisotropy model in the fitting. This is the β-MBM case. Finally, in the β-HM case, we adopt an anisotropy profile that varies linearly with the logarithmic slope of the number density profile: with a = −0.15 and b = −0.19, as was determined on a vari- Figure 6. Same as Fig. 3 for the different cases of parameter constraints (see text and Table 3). As in Fig. 3, results for the anisotropy models Cst, T, and ML are shown, left to right, in green, blue, and red, respectively, and here the results for the β-HM anisotropy model are shown in black. Open symbols refer to those parameters that are not free parameters in the MAMPOSSt analysis for the Case considered. ety of simulations by . Note that Sparre & Hansen (2012) have argued that cosmological haloes are better described by this relation than by the attractor of Hansen, Juncher, & Sparre (2010), although Lemze et al. (2012) find that the validity of this relation is limited to the central regions of haloes. By fixing the anisotropy to determine the parameters of the mass profile, these cases are similar in spirit to mass inversion (Mamon & Boué 2010;Wolf et al. 2010), except that the former requires a parametric form for the mass profile, while the latter suffers from binning and extrapolation.
So, while the anisotropy parameter is fixed in the β-iso case, it changes according to rρ (a free parameter here) in the β-MBM case, and according to the logarithmic slope of the tracer number density profile (which is related to the free parameter rν) in the β-HM case.
In Table 4, we list the results of our analysis for the different observational cases for the 500-particle samples. These results are graphically displayed in Fig. 6. For the sake of comparison, we list and plot again here the results for the Gen case, already displayed in Table 2 and Fig. 3.
The results for r200 are essentially independent of the case considered. While the r200 values for the EVir case are directly obtained from the Clean interloper removal scheme (App. B), they are  Table 1. measured with similar bias and inefficiency as for the Gen case. In other words, σz can be used to provide an estimate of r200 with a comparable accuracy to that provided by MAMPOSSt. The advantage of MAMPOSSt is of course that it provides estimates for the other dynamical parameters at the same time.
Fitting rν from the R distribution only, i.e. externally from MAMPOSSt, does not significantly alter the accuracy returned for this parameter. This means that incompleteness in the sample of tracers, if properly accounted for, is not a significant issue for MAMPOSSt. Moreover, this lifts our concern that the Split method does not find the same minimum for − ln L as the Gen (joint) method (see end of Sect. 2.1). Using rρ to predict rν (the TLM case) leads to larger bias and inefficiency on rν, but this occurs because the bias and inefficiency on rρ are worse than those on rν .
Also, the results for rρ depend little on the case considered. The best results are obtained for the ΛCDM case, where rρ is not a free parameter of the MAMPOSSt analysis, but it is estimated from r200 using a theoretical relation. The good results obtained for the ΛCDM case are not surprising, given that r200 is better constrained than rρ in MAMPOSSt, and that the test haloes considered here are extracted from a ΛCDM cosmological simulation, and have therefore similar properties to those used by Macciò et al. (2008) to establish the mass-concentration (and hence the r200 − rρ) relation.
As for the anisotropy parameters, the ML model remains impossible to constrain in all cases. For the Cst and T models, the bias and inefficiency do not depend strongly on the case considered. However, the o vs. t correlation values do depend on the case: a significant correlation is obtained for A∞ for the Gen, Split, ΛCDM and TLM cases, but not for KVir, EVir, or MfL. For the ΛCDM and β-MBM cases, a significant correlation is found even for the poorly constrained r β parameter, but note that in the β-MBM case, r β is not a free parameter of the MAMPOSSt analysis, but is forced to match rρ.   Table 3 It is surprising that MAMPOSSt obtains worse inefficiencies in the KVir case where r200 is assumed perfectly known in comparison to the Gen, Split or EVir cases. The reason for this is probably related to the halo triaxialities. When we observe a halo along a given line-of-sight, we are sampling its velocity distribution only along that line-of-sight. If the components of the halo velocity distribution along different axes are different, this difference will be reflected in our results.

Gen
This point is demonstrated in Fig. 7. There, we show the correlation between the ratio of r200 measured by MAMPOSSt to the true value measured in 3D versus the ratio of velocity dispersions within the true virial sphere, measured along the projection axis over that measured globally. The very high correlation (Spearman rank: r = 0.82) shows that the error on r200 is related to the velocity dispersion ratio above, which directly measures the triaxiality of the halo (without mixing with the effects of interlopers beyond the virial sphere).

Stacked haloes
To extract all the possible information from the available data, it is a common practice to construct stacked samples (e.g. Carlberg et al. 1997;Biviano & Girardi 2003;Katgert et al. 2004;Biviano & Poggianti 2009). In the case of clusters of galaxies, this is done by scaling the galaxy distances from their cluster centres and the galaxy velocities with respect to their cluster mean velocity, by their cluster r200 and v200, respectively, where v200 ≡ (GM200/r200) 1/2 . To do this, prior knowledge of the individual cluster r200 values is needed. Trying to mimic the observational situation as close as possible, we stack our haloes using the values of r200 directly obtained from the Clean interloper removal scheme, as described in Appendix B, i.e. the values used in what we called the EVir case in Sect. 3.4. We thus build three stacked haloes, one for each projection axis, from the eleven 500-particle haloes, each first passed through the Clean interloper removal scheme (app. B).
Since we need to fix the r200 values of the haloes before stacking, r200 cannot be a free parameter of the MAMPOSSt analysis. Moreover, in an effort to mimick the observational procedure, we estimate the individual rν values of the 11 haloes also before stacking, and then we take the biweight average of these values as representative of the rν of the stack. In the real world, this is done  Beers et al. 1990).
because individual haloes often suffer from different spectroscopic incompleteness levels. More precisely, we take the rν estimates obtained in the EVir case (see Sect. 3.4). The remaining two parameters, rρ and the anisotropy parameter are estimated via the MAM-POSSt analysis. Results for the six stacked haloes are listed in Table 5 and displayed in Fig. 8. We use the gapper estimate of dispersion (Wainer & Thissen 1976, see Beers et al. 1990) given our small sample (N = 3). In Fig. 9, we show 1-, 2-, and 3-σ contours in the log rρ − log A∞ plane (for the T anisotropy model). In this figure, we indicate with a cross the expected solution for the stacked haloes. This Figure 8. MAMPOSSt residuals for stacked haloes (labelled 'S500' and 'S100', when built from the 11 haloes, respectively sampled with 500 and 100 particles), and biweight mean of the values obtained for the 11 haloes in the EVir case (see Sect. 3.4, labelled 'M500' and 'M100', for the 500 and 100 particle haloes, respectively). Different symbols identify different projections, x-axis: diamonds, y-axis: squares, z-axis: circles. In the upperleft panel (rρ), the green, blue, and red symbols are the results obtained using, respectively, the Cst, T, and ML anisotropy models.
is the biweight-average (with its error) of the true rρ and A∞ values of the 11 individual haloes from which the stacked haloes are constructed. These values are also listed in the last line of Table 1.
The average values of r−2 and r β are identical, in substantial agreement with what was found by MBM10. Note also that the average value of A∞ corresponds to β∞ = 0.5, meaning that the average ML and T anisotropy models are the same. Fig. 9 shows that the expected solution for the stack sample is always within the 1-σ contour of the MAMPOSSt result for the T model. While not shown, this is also the case for the Cst model, but not for the ML model, where r β remains essentially unconstrained, as was the case for the individual haloes. The estimates of the parameter rρ do not depend on the assumed anisotropy model. They are biased low for the stacks built from the 100-particle samples, even if not significantly so (see Fig. 9). Interestingly, while the bias is stronger for the 100-particle stack than for the 500-particle one, the dispersion of the recovered values appears to be lower. The dispersion is instead higher, as expected, for the anisotropy parameters, except for A∞ where the dispersion is comparable for the 500-particle and the 100-particle stacks. However, with only three results per anisotropy model for each Stack sample, these differences in bias and dispersions do not appear to be statistically significant.
In Table 5, we also list the biweight average values of the 11 individual haloes along each projection, obtained in the EVir case (see Sect. 3.4), and these values are also plotted in Fig. 8 (labelled 'M') together with the best-fit solutions for the stacked haloes (labelled 'S'). Comparing the best-fit values of rρ and the anisotropy parameters obtained from the stacked haloes with those found from the mean of the individual 11 haloes (along the same projection), Figure 9. MAMPOSSt confidence contours in the log rρ − log A∞ plane for the stacked haloes constructed from single-axis projections of 11 haloes from the 500-particle (left panels) and 100-particle (right panels) samples. The x-axis, y-axis, and z-axis stacks are shown in the top, middle, and bottom panels, respectively. The contours are drawn at 1-, 2-and 3-σ. The best-fit solutions are indicated by the filled circles, while the crosses represent the biweight average of the true values of the parameters for the 11 haloes. The cross lengths across the xand y-axes represent the errors on the two parameter averages (see Table 1).
it is difficult to decide whether it is better to run MAMPOSSt on the stacked halo or to adopt means of the individual MAMPOSSt estimates. A close inspection of Table 5 and Fig, 8 indicates that, for N = 500, stacks are more accurate for estimating r β and A∞, while means of individual MAMPOSSt estimates appear more accurate for estimating rρ and A. But the statistics are poor. In the real observational situation, it therefore seems advisable to consider both the average of the results of the individual haloes and the result obtained for the stacked sample.

A large halo
It is interesting to consider whether we can achieve the same accuracy using a single halo, but with a total number of tracer particles similar to that of the stack. We extract ∼ 5000 particles within r200 for the more massive halo in our sample, 5726, remove interlopers with the Clean method, and run MAMPOSSt for the general case, adopting the NFW mass density and different anisotropy models. Table 6 compares the results for the 3 cones of 500 particles (counted before interloper removal) and the 3 with ≃5000 particles. We use again the gapper dispersion, given our small sample (N = 3). As expected, the biases and inefficiencies are generally reduced when selecting all (≃5000) particles instead of randomly selecting 500. For the T anisotropy model, the inefficiencies in log(o/t) are reduced by 0.1 dex for rρ and anisotropy, 0.15 dex for rν, but virtually unchanged for r200. Surprisingly, for Cst and ML anisotropy, r200 is reached with worse inefficiency with the full 5000-particle sample in comparison with the 500-particle one.
However, one should not over-interpret these comparisons. In- deed, our Monte-Carlo tests on 10 000 random samples of 3 objects arising from a Gaussian distribution indicate that while the gapper standard deviation is unbiased even for as few as 3 objects, it has an inefficiency of 0.53 times its true value for N = 3. In other words, an inefficiency of log(o/t) of 0.1 would have roughly 0.05 accuracy. Therefore, any improvement or worsening of inefficiency in log(o/t) of less than a factor of 2 is clearly not statistically significant. Nevertheless, the fact that most inefficiencies are reduced when increasing the sample size from 500 to 5000 is much more significant, and is indeed what is qualitatively expected. However, the expected quantitative improvement of log √ 10 = 0.5 does not appear to be attained for the inefficiency of log(o/t).
The bias and inefficiencies for the 5000-particle sample are typically higher for the single halo than for the 11 × 500 particle stack (with just 10% more particles, although the extraction there was for case EVir instead of Gen). But again, the inefficiencies are within a factor 2, so these changes are individually not statistically significant. But the combination of all of them appears to be statistically significant.
The worse inefficiencies for the large halo in comparison with the stacked halo of similar number of tracers can be due to either our assumption of Gaussian velocity distributions (see Sect. 2.2) or to the triaxiality of the halo velocity ellipsoids. The velocity distribution of the stack halo approaches a Gaussian by the central limit theorem, so our Gaussian assumption is better verified in a stack sample than in individual haloes. As a matter of fact, a number of studies have shown the 3D velocity distributions of individual haloes to deviate from Gaussianity (Wojtak et al. 2005;. But perhaps the dominant effect is that of triaxiality. As discussed at the end of Sect. 3.4, the results obtained by MAM-POSSt for r200 are influenced by the choice of the projection axis, and in the present case, the r200 values obtained for 5726 along the x, y, and z projections when using 5000 particles, are ordered in the same way as those obtained when using 500 particles.

CONCLUSIONS AND DISCUSSION
We have presented the formalism for a new method, called MAM-POSSt, for the determination of the mass and anisotropy profiles of spherical systems, assumed to be in dynamical equilibrium, in particular, as stationary systems with no streaming motions. In MAM-POSSt, the distribution of tracers in projected phase space are fit by the predicted distribution arising from assumed 3D radial profiles of the tracer density, total mass, and velocity anisotropy, as well from an assumed family of shapes for the radial and tangential components of the 3D velocity distribution.
MAMPOSSt has several important advantages over other methods for mass/anisotropy modeling: 1) MAMPOSSt does not involve binning the data, hence its results are independent of the choice of such radial bins, in contrast with methods based upon velocity moments, including anisotropy and mass inversions; 2) MAMPOSSt does not involve interpolations and extrapolations of binned radial profiles of the observables, again in contrast with anisotropy and mass inversions; 3) MAMPOSSt does not require differentiation of the observed projected pressure, Σ(R) σ 2 z (R), once more in contrast with anisotropy and mass inversions; 4) MAMPOSSt extracts accurate constraints on the velocity anisotropy, in contrast with methods that assume Gaussian LOS velocity distributions; 5) MAMPOSSt is very fast, as it involves a single integral (for popular β(r) profiles). Indeed, the calculations for the best fit parameters and their marginal and correlated distributions through MCMC for a 500-tracer system (as displayed in Fig. 2) require roughly 4 minutes of CPU time on a standard desktop personal computer. By contrast, distribution function methods involve triple integrals, and therefore take typically 1000 times longer to run, i.e. a few days for 500-tracer systems. Orbit modeling techniques are even slower and cannot properly probe parameter space.
In this work, we have extensively tested MAMPOSSt, for the case of Gaussian 3D velocity distributions, on a set of 11 clustermass haloes extracted from cosmological simulations. The results of these tests indicate that, for systems with 500 velocities, MAM-POSSt provides essentially unbiased estimates of the relevant mass and velocity anisotropy profile parameters, with inefficiencies of 10% for the virial radius, 20% for the constant or outer velocity anisotropy, 27% for the tracer scale radius, but as high as 48% for the scale radius of the DM. However, MAMPOSSt seems unable to set constraints on the radius of transition of a gently rising anisotropy model such as ML.
We have noted that the results of MAMPOSSt are similar when we inserted the virial radius among the parameters to be jointly fit or when we derived the virial radius using the new interloper removal scheme that we presented in Appendix B. We found that this interloper-removal algorithm performs as well or better than other methods. As shown in Fig. 10, the errors on r200 from MAMPOSSt are highly correlated with those obtained from the Clean procedure. Hence, the new interloper removal scheme produces extremely fast estimates of the virial radius, based on the LOS velocity dispersion, that are very close to those obtained with the full MAMPOSSt procedure.
The correlation highlighted in Fig. 10 suggests that the quality of MAMPOSSt estimates is limited by the accuracy of the interloper removal. However, we also saw (Fig. 7) that the error on log r200 is correlated with the ratio of velocity dispersions within the true virial sphere, measured along the projection axis over that measured globally. This second correlation indicates that our capability of recovering the true values of r200 is limited by the effects of triaxiality, i.e. by the fact that we only observe one component of the 3-dimensional velocity distribution. The effect of triaxiality may be even more important than the uncertainties in the interloper cleaning process.
Our additional tests lead to the following conclusions: 1) Best results are obtained setting all 4 parameters free and using the T anisotropy model.
2) Setting the virial radius to the true value (e.g. measured through other techniques, such as X-rays or lensing) leads to worse results on the other parameters than when the virial radius is an additional free parameter. We argue that this is caused by the triaxiality of the haloes (Fig. 7).
3) MAMPOSSt works surprisingly well with samples of 100 tracers instead of 500, only somewhat better with samples of 5000 tracers, even better with stacked haloes of 11 × 500 tracers, similar to the mean of these individual tracers, but not well as expected from N −1/2 arguments. 4) MAMPOSSt is not very sensitive to the chosen aperture, when this is to within ≈35% of the true value of r200. Table 7 compares the accuracy of MAMPOSSt with that of the dispersion-kurtosis method of Łokas (2002, tested by Sanchis, Łokas, & Mamon 2004 and the DF method of Wojtak et al. (2009). In comparison with the dispersion-kurtosis method, MAMPOSSt does slightly less well on the inefficiencies of the virial radius and worse on that of the DM scale radius (but better on the biases of both radii), but MAMPOSSt performs better for the constant anisotropy (both bias and inefficiency). Moreover, MAMPOSSt has the advantage of being flexible enough to infer the outer anisotropy, while the dispersion-kurtosis method of Łokas (2002) is limited to constant anisotropy (but this limitation to the dispersion-kurtosis method has been recently lifted by Richardson & Fairbairn 2012). Comparing with the much slower DF method, MAMPOSSt does somewhat worse on the inefficiency on the scale radius, with the same bias, yet MAMPOSSt is respectively slightly (not significantly) and much better on the inefficiencies on the virial radius and outer anisotropy, with less bias on the former. Note, however, that we forced here the inner velocities to be isotropic, while these were free parameters in the DF method and that the DF method computes the outer anisotropies at 5 rs instead of the asymptotic value. Nevertheless, we found ( Table 4) that MAMPOSSt always performs as well or better when all parameters are set free, and we would expect the same for the DF method. Moreover, we see no reason why the inefficiency of the anisotropy measured near the virial radius should be worse than that of its asymptotic value.
Despite the very good results of our tests, it must be noted that although the Gaussian MAMPOSSt recovers the 2nd velocity moment (eq. [30]) , it does not recover the fourth moment. Indeed, taking the fourth moment of the velocity distribution of equation (28) leads to Σ(R) v 4 z (R) = Notes: The biases and inefficiencies on virial masses for the dispersion-kurtosis and distribution function methods were divided by 3 to convert to biases and inefficiencies for virial radii. All analyses are for case TLM, and for MAMPOSSt the Clean interloper removal scheme is used. Dispersion-kurtosis log(o/t) for the DM scale radius were taken from the analogous values that Sanchis et al. (2004) obtained for the concentration parameter (basically assuming that in comparison, the errors on the virial radius are small). For the distribution function method, values of log(o/t) for the outer anisotropy were measured at roughly the virial radius, while the MAMPOSSt comparison (last row) is for the T anisotropy model (eq. [38]). The biases and inefficiencies are estimated with the biweight mean and dispersion, except for the dispersion-kurtosis method, where they are taken from the 2nd line of in each of the three cartesian directions. We limited these 3 × 105 = 315 haloes to projected radii within Rmax = Xmax r200 from the