# Escape Velocity Mass of Abell S1063

Alexander Rodriguez University of Michigan, Ann Arbor
Department of Astronomy  Christopher J. Miller University of Michigan, Ann
Arbor
Department of Astronomy  University of Michigan, Ann Arbor
Department of Physics  Vitali Halenka Independent Researcher  Anthony Kremin
Lawrence Berkeley National Laboratory, Berkeley, CA 94720

###### Abstract

We measure the radius-velocity phase-space edge profile for Abell S1063 using
galaxy redshifts from Karman et al. (2015) and Mercurio et al. (2021).
Combined with a cosmological model and after accounting for interlopers and
sampling effects, we infer the escape velocity profile. Using the Poisson
equation, we then directly constrain the gravitational potential profile and
find excellent agreement between three different density models. For the NFW
profile, we find log10(M${}_{200},{\rm crit}$)= $15.40^{+0.06}_{-0.12}$M⊙,
consistent to within $1\sigma$ of six recently published lensing masses. We
argue that this consistency is due to the fact that the escape technique
shares no common systematics with lensing other than radial binning. These
masses are 2-4$\sigma$ lower than estimates using X-ray data, in addition to
the Gómez et al. (2012) velocity dispersion estimate. We measure the 1D
velocity dispersion within r200 to be $\sigma_{v}=1477^{+87}_{-99}$ km/s,
which combined with our escape velocity mass, brings the dispersion for AS1063
in-line with hydrodynamic cosmological simulations for the first time.

Galaxy clusters, Dark energy, Cosmology, Gravitational lensing

††journal: ApJ

##  1 Introduction

Galaxy clusters represent the most massive gravitationally bound structures
having formed to date. Outside of their cores, these systems are unique in
that their evolution is determined only through the combined effects of
gravity and the expanding spacetime from the accelerated universal expansion
on $\sim\text{Mpc}$ scales.

In the current $\Lambda$CDM framework, the cluster potential is dominated by
dark matter (hereafter DM). A galaxy in the cluster feels an inward force from
not only the DM, gas, and galaxies, but also an outward force due to the
accelerated universal expansion. In a non-accelerating universe, dynamical
tracers of the cluster potential can escape a cluster if it can somehow reach
a speed above the escape velocity, given by

$\displaystyle v_{\text{esc}}{\color[rgb]{0,0,0}{(r)}}$
$\displaystyle=\sqrt{-2\Psi(r)},$ (1)

where $\Psi(r)$ is the potential from matter alone.

In an accelerating universe, Nandra et al. (2012) showed that there is a
radius ($r_{\text{{eq}}}$) at which the inward force on a test particle by a
massive object is perfectly balanced by the effective outward force of the
accelerating space time. In the accelerating case, Behroozi et al. (2013)
showed that a form of Equation 1 holds as well, where the matter-only
potential $\Psi(r)$ is replaced by the effective potential, $\Phi(r)$, which
includes contributions from both gravity and the accelerated expansion. The
key difference being that the escape speed relative to the clusters becomes
zero at $r_{\text{{eq}}}$.

Arguably, the most important factor in determining the evolution of a cluster
is its mass. Diaferio & Geller (1997) showed that it is possible to relate the
mass of a cluster to its escape velocity through construction of a
radius/velocity phase-space diagram. Within the phase-space, particles, sub-
halos, and galaxies are confined to regions within which they are bound,
confined by the phase-space edge, a “caustic”, or sharp cut-off in the phase-
space with a velocity at a corresponding radius given by Equation 1.

Simulations have been used to show that the escape velocity can be measured
from a clearly defined edge in the radius/velocity phase-space diagram. Except
for a very small number of escaping tracers, only those with the maximum
possible radial or tangential 1D speed will contribute to this edge (Behroozi
et al., 2013). The power of utilizing the observed $v_{\text{esc}}(r)$ is in
its direct connection to the total potential, enabling cluster mass
estimations, tests of gravity on the largest scales in the weak field limit,
and placing constraints on the $\Lambda$CDM cosmological parameters (Miller et
al., 2016; Stark et al., 2016, 2017). Notably, the cosmic acceleration results
in a smaller escape velocity for a given mass than one would expect in a
universe without acceleration (Diaferio & Geller, 1997; Behroozi et al., 2013;
Miller et al., 2016; Stark et al., 2017).

Prior work associated with utilizing observed phase-spaces for measuring
cluster masses has required cosmological simulations to calibrate “caustic”
surfaces (Diaferio & Geller, 1997; Diaferio, 1999; Diaferio et al., 2005;
Serra et al., 2011; Geller et al., 2013; Gifford et al., 2013; Pizzardo et
al., 2023, 2022). These caustics are then related to cluster masses or mass
profiles. However, Gifford et al. (2017) and Halenka et al. (2022) employed
direct measurements of the observed maximum phase-space peculiar velocity
edge. Gifford et al. (2017) defined the edge as the 90th percentile of the
peculiar velocity distribution to infer average cluster masses. Halenka et al.
(2022) measured the absolute maximum velocity in projected data to infer the
3D escape profile. We will use Abell S1063 as an observational test of how
well the observed phase-space maximum edge profile relates to the cluster
mass.

Abell S1063 (hereafter AS1063)–also referred to as RXC J2248.7-4431— is a
massive, luminous, cluster, located at a redshift of $z=0.3475$ (Abell et al.,
1989; Böhringer et al., 2004), and is part of the Cluster Lensing And
Supernova survey with Hubble (HST CLASH) program (Postman et al., 2012; Lotz
et al., 2017). This cluster has had its mass measured many times using
independent techniques, making it an ideal choice for our analysis.

Gómez et al. (2012) performed an in-depth investigation of properties of
AS1063 by measuring its X-ray temperature and velocity dispersion (given in
Table 1). These were obtained through Gemini observations, and found a
corresponding mass of $\log{M_{200c}}\sim 15.6\log(M_{\odot})$ 111$M_{200c}$
refers to the mass within which the mean density of the cluster is
$200\times\rho_{c}(z)$, where $\rho_{c}(z)$ is the critical density of the
universe at redshift $z$., extremely high relative to most known clusters.

Gruen et al. (2013) subsequently performed an analysis using weak lensing, and
found a mass nearly $\sim 2\times$ lower (again see Table 1 for further
details), attributing the discrepancy to bulk motion of different galaxy
populations along the line-of-sight, evidence for a dynamically non-relaxed
system.

Umetsu et al. (2014, 2016) repeated the analysis of Gruen et al. (2013) with
the same CLASH-VLT data, but with independent shear calibration measurements
and photometry, resulting in findings which more closely matched Gruen et al.
(2013) than Gómez et al. (2012).

Sartoris, B. et al. (2020) then combined the CLASH data with Caminha et al.
(2016), where 16 additional background sources for the cluster had been
discovered. These data were used for a shear + magnification lensing analysis,
consistent with the previous lensing estimates from Gruen et al. (2013);
Umetsu et al. (2014, 2016).

Due to its extreme X-ray luminosity, Comis et al. (2011); Gómez et al. (2012)
also studied the cluster using its X-ray properties, using data from Chandra.
The resulting X-ray masses were similar to the dynamical mass using
$\sigma_{v}$ from Gómez et al. (2012), a systematic offset between the
dynamical and X-ray masses compared to lensing estimates. Donahue et al.
(2014) repeated the X-ray analysis XMM-Newton and Chandra telescopes, and
again found similar results for the X-ray mass.

Subsequently, the prominent SZ signal in the CMB sparked the Planck team (Ade
et al., 2011) to use AS1063 for SZ analyses. Williamson et al. (2011); Penna-
Lima et al. (2017); Capasso et al. (2019) measured SZ masses which are
significantly lower than the X-ray estimates and more consistent with the
lensing results.

Most recently, Sartoris, B. et al. (2020) performed another dynamical analysis
of AS1063, this time making use of MAMPOSSt (Mamon et al., 2013). This
technique involves inversion of the Jean’s equation of dynamical equilibrium
(Binney & Tremaine, 2008) to perform a joint maximum likelihood fit to the
velocity dispersion profile of the BCG and to the velocity distribution of
cluster member galaxies. The resulting dynamical mass is higher than the
lensing and SZ estimates, although in better agreement than the dynamical
estimate from Gómez et al. (2012) using the velocity dispersion.

Mercurio et al. (2021) expanded upon these recent findings by combining CLASH-
VLT observations with Karman et al. (2015), and found strong statistical
evidence of merging subclumps in the system. These subclumps were in addition
to a bi-modal distribution of peculiar redshifts–further evidence that the
cluster is not dynamically relaxed.

Thanks to the CLASH-VLT spectroscopic program, the radius/velocity phase-space
of AS1063 has been mapped with thousands of galaxies. Alongside the numerous
prior independent mass measurements, these data enable us to conduct the first
escape edge analysis of a cluster. Notably, we will not be modeling the
internal structure of the phase-space data. In fact, only a few 10s of
galaxies in the phase-space will ultimately contribute to the edge measurement
and used in the statistical inference. We will use simulations and analytical
phase-space representations of AS1063-like systems to quantify the accuracy
and precision of our edge measurement.

The paper is structured as follows: in §2, we present a theoretical overview
of the methodology of the escape edge framework. In §3, we perform a detailed
measurement of the phase-space edge profile and its error. Mass inference
based on the escape profile is presented in §4. In §5, we present a new
estimate of the velocity dispersion where we model the many contributions to
its uncertainty. In §6, we discuss the astrophysical consequences of our
findings, and implications.

Throughout this paper, we adopt
$H_{0}=72\,\text{km}\,\text{s}^{-1}\,\text{Mpc}^{-1}$, $\Omega_{M}=0.27$,
$\Omega_{\Lambda}=0.73$. We call $r_{\Delta}$ the radius that encloses an
average density $\Delta$ times the critical or mean density at the cluster
redshift. In all plots and results, we use the critical threshold unless
explicitly stated otherwise.

Table 1: Summary of the various mass estimation methods from previous work.
The Intracluster Medium (ICM) masses require a hydrostatic equilibrium
correction, which is incorporated in some of the ICM measurements below.
Masses are measured with respect to $\rho_{\text{crit}}$. In Figure 8, we show
these mass estimates compared to our mass constraint from the escape velocity,
in addition to the theoretical expectation from simulations.

AS1063 lensing log mass estimates | $M_{200c}(M_{\odot})$
---|---
Shear1 (Gruen et al., 2013) | $15.37^{+0.11}_{-0.10}$
Shear (Melchior et al., 2015) | $15.24^{+0.10}_{-0.10}$
Shear (Klein et al., 2019) | $15.22^{+0.12}_{+0.10}$
\+ magnification1 (Umetsu et al., 2014) | $15.32^{+0.13}_{-0.18}$
\+ magnification (Sartoris, B. et al., 2020) | $15.34^{+0.10}_{-0.15}$
\+ mag + strong1 (Umetsu et al., 2016) | $15.28^{+0.13}_{-0.19}$
AS1063 ICM HSE log mass estimates | $M_{200c}(M_{\odot})$
X-ray2,4 (Chandra) (Donahue et al., 2014) | $15.60^{+0.04}_{-0.05}$
X-ray2,4 (XMM) (Donahue et al., 2014) | $15.52^{+0.02}_{-0.02}$
X-ray3 (Chandra) (Sartoris, B. et al., 2020) | $15.57^{+0.07}_{-0.08}$
SZ3,5 (SPT) (Williamson et al., 2011) | $15.34^{+0.08}_{-0.10}$
SZ2,6 (Planck) (Penna-Lima et al., 2017) | $15.25^{+0.03}_{-0.01}$
AS1063 dynamical log mass estimates | $M_{200c}(M_{\odot})$
$\sigma_{v}$7 (Gómez et al., 2012) | $15.60^{+0.15}_{-0.11}$
Jeans8 (Sartoris, B. et al., 2020) | $15.46^{+0.04}_{-0.05}$

  * 1

We convert their $M_{200,c}$ in a $\Omega_{M}=0.27$ universe to a $M_{200,c}$
in a $\Omega_{M}=0.3$ universe.

  * 2

This Work accounts for a hydrostatic bias correction.

  * 3

This Work does not account for a hydrostatic bias correction.

  * 4

We convert their $M_{2500,c}$ to $M_{200,c}$ using the (Duffy et al., 2008a)
mass-concentration relation.

  * 5

From a Mass versus S/N scaling relation calibrated to simulations. We convert
their $M_{200m}$ to w.r.t. to critical as described therein. We use the
average of statistical and systematic errors.

  * 6

Calibrated using X-ray scaling relations. We convert their $M_{500,c}$ to
$M_{200,c}$ using the (Duffy et al., 2008a) mass-concentration relation.

  * 7

They measure $\sigma_{v}$ to $\sim 750kpc$ and estimate what it should be at
$r_{200}$ using N-body simulations. We use this “de-biased core” $\sigma_{v}$
in the Evrard et al. (2008) scaling relation to represent a dynamical mass.

  * 8

This is the total mass: dark matter + gas.

##  2 Theoretical Background

###  2.1 Escape dynamics

Nandra et al. (2012) studied the dynamical effects of an accelerating
spacetime on a galaxy cluster in the weak-field limit of general relativity.
In this paradigm, we can relate the acceleration to the corresponding
effective potential, $\Phi(r)$, with $a_{\text{eff}}=\nabla\Phi(r)$, where
$a_{\text{eff}}$ is the effective acceleration on the particle

$\displaystyle a_{\text{eff}}$
$\displaystyle=\nabla\Phi(r)=\nabla\Psi(r)+q(z)H^{2}(z)r\hat{r}$ (2)

Here, $\Psi(r)$ is the potential due to dark matter, ICM, and galaxies, and
$H(z)=H_{0}\sqrt{\Omega_{\Lambda}+\Omega_{M}(1+z)^{3}}$, with
$q(z)=\frac{1}{2}\Omega_{M}(z)-\Omega_{\Lambda}(z)$ (assuming a flat universe,
or $\Omega_{k}=0$) (Sandage, 1970). The second term in Equation 2 can be
thought of as the deviation from the classical Newtonian regime, quantifying
the curvature produced by the acceleration of the expansion of spacetime. A
dynamic tracer inside a galaxy cluster will experience two forces: an inward
gravitational force due to the internal mass-energy and an additional outward
effective force due to the acceleration of the expansion of spacetime.

Unlike the classical Newtonian regime, in Equation 2, there is a radius
relative to a cluster’s gravitational center where the acceleration due to
gravity balances the acceleration due to the expansion of universe (Behroozi
et al., 2013). This location $r_{\text{eq}}$ is given by
$r_{\text{eq}}=(\frac{GM}{-q(z)H^{2}(z)})^{1/3}$, occurring where
$\nabla\Phi=0$. It is a surface that depends on cosmology and the mass
internal to it. At the redshift of AS1063, we find $r_{\text{eq}}\sim
10\,r_{200}$ for a AS1063-size cluster.

Potentials are only physically meaningful in a relative sense (i.e., escape is
to some point as opposed to infinity). As shown in Miller et al. (2016), we
need to enforce the constraint that a tracer has a relative escape speed of
zero at $r_{\text{{eq}}}$:

$\displaystyle v_{\text{esc}}(r)$
$\displaystyle=\left[-2(\Psi(r)-\Psi(r_{\text{eq}}))-q(z)H^{2}(z)(r^{2}-r^{2}_{\text{eq}})\right]^{1/2},$
(3)

where the Poisson equation relates the gravitational potential and the density
out to $r_{\text{{eq}}}$

$\Psi(r)=-{\rm 4\pi G}\Big{[}\frac{1}{r}\int_{0}^{r}\rho(r^{\prime})r^{\prime
2}dr^{\prime}+\int_{r}^{r_{\text{{eq}}}}\rho(r^{\prime})r^{\prime}dr^{\prime}\Big{]}.$
(4)

Behroozi et al. (2013) used $v_{\text{esc}}$ and $r_{\text{eq}}$ to measure
the bounded fraction of dark matter particles in simulations. Miller et al.
(2016) used $r_{\text{eq}}$ to show that the dark matter and sub-halo escape
speed profiles accurately and precisely match expectations from the cluster
density profiles. They explored three different analytical density profiles
commonly used in the literature: the classic “NFW” (Navarro et al., 1997), the
“Einasto” density profile (Einasto, 1969), and the “Dehnen” profile (Dehnen,
1993), for which the two latter profiles both fall off more quickly at large
radii compared to the NFW.

Miller et al. (2016) found from simulations that while the Einasto and Dehnen
density profiles can accurately predict the escape profile in the phase-space
to within $3-5\%$, the NFW profile will tend to overestimate the potential at
all radii by $\sim 10\%-15\%$, given the NFW is shallower at larger radii than
the Einasto and Dehnen profiles. In this work, we intend to explore these
simulation predictions for AS1063 specifically.

###  2.2 From theory to observable

Prior work has focused on using the “caustic technique” to infer the mass
profile (Diaferio & Geller, 1997). This technique employs the phase-space data
and identifies a number of possible escape edges in the 2D radius/velocity
histogram of a cluster phase-space. From these, a velocity “caustic” is chosen
such that the average of the square of a phase-space isodensity contour is
less than or equal to 4 times the square of the velocity dispersion over some
radial range. This calibration to the dispersion stems from the virial energy
equilibrium of clusters without an accelerating spacetime (see Gifford &
Miller (2013) for a detailed explanation). In this context, the escape edge is
calibrated by the projected velocity dispersion. In order to use it to infer
the radial escape edge, Diaferio & Geller (1997) and Diaferio (1999) suggest a
geometric model such that:

$\displaystyle\langle v_{\text{esc, los}}^{2}(r_{\perp})\rangle$
$\displaystyle=\frac{1-\beta_{r_{\perp}}}{3-2\beta_{r_{\perp}}}\langle
v_{\text{esc}}^{2}(r_{\perp})\rangle$ (5)
$\displaystyle=g^{-1}(\beta_{r_{\perp}})\langle
v_{\text{esc}}^{2}(r_{\perp})\rangle,$

where $\langle\cdot\rangle$ refers to an ensemble average of the galaxies
within some radial bin and the measurements are made along lines-of-sight
(los) with increasing projected radius $r_{\perp}$. Here, $\beta_{r_{\perp}}$
is the velocity anisotropy measured in projection where

$\beta(r)=1-\sigma_{t}^{2}/\sigma_{r}^{2},$ (6)

and $r$ and $t$ refer to radial and tangential components of the velocities.

For the anisotropy-based suppression, typical values of $\beta$ for clusters
range from $\beta=0$ to $\beta=0.5$ (Iannuzzi & Dolag, 2012; Munari et al.,
2013; Stark et al., 2019; Capasso et al., 2019). With this expected level of
$\beta$, the observed $v_{\text{esc, los}}$ should be suppressed compared to
true $v_{\text{esc}}$ by 50% - 57%, regardless of the number of galaxies
observed.

Equation 5 has some interesting limiting properties. First, the suppression
quickly approaches a constant for increasingly negative values of $\beta$ such
that $g^{-1}(\beta\rightarrow-\infty)\sim 0.7$. In other words, all projected
phase-space escape edges must be suppressed by at least 30% compared to the 3D
edge. Second, as the anisotropy becomes radial ($\beta\rightarrow 1$), the
suppression grows to $g^{-1}(\beta)\rightarrow 0$, such that no edge would be
detectable in projected data when the tracers are dominated by motions
corresponding to radial infall. Third, for isotropic systems ($\beta=0$), the
suppression is $g^{-1}(\beta)=1/\sqrt{3}$ (Halenka et al., 2022).

The number of tracers does not appear in equation 5. Yet Serra et al. (2011)
noticed a dependence on the number of tracers in simulation phase-spaces.
Motivated by this issue and prior work by Gifford et al. (2017),  Halenka et
al. (2022) took another approach. They recognized that there is always the
possibility that a tracer is observed escaping along the line-of-sight. As one
fills the phase-space with additional tracers, more such galaxies should
exist. They proposed that the suppression in observed escape edges should
simply be due to low statistical sampling of the phase-space. Halenka et al.
(2022) utilized a purely analytical approach to sample cluster phase-spaces
via action-based modeling (Vasiliev, 2019). They show that in an idealistic
scenario where hundreds of thousands of tracers are projected in velocity and
radius, the true escape profile can be accurately observed in projected data.
They then show quantitatively that the suppression of the escape edge is
dominated by statistical sampling. Anisotropy, cosmology, and cluster mass
play very minor roles (e.g., a maximum of a few percent).

Halenka et al. (2022) provide a model of the suppression of the projected
escape profile :

$\displaystyle\langle v_{\text{esc, los}}\rangle(r_{\perp})$
$\displaystyle=\frac{\langle v_{\text{esc}}\rangle(r_{\perp})}{Z_{v}(N)}$ (7)
$\displaystyle Z_{v}(N)$
$\displaystyle=1+\left(\frac{N_{0}}{N}\right)^{\lambda},$ (8)

where $N$ is the number of tracers in the projected phase-space between
$0.3\,r_{200}\leq r\leq r_{200}$. $N_{0}$ and $\lambda$ are different
depending on the data used to calibrate the suppression. For instance, the
Action-based Galaxy Modeling Architecture (AGAMA) (Vasiliev, 2019) can be used
to generate analytic mock phase spaces. Similarly, semi-analytic galaxies from
massive halos in the Millennium N-body can also be employed. Both datasets
produce phase-spaces which include interlopers. Galaxies with large 3D
distances can sometimes be projected onto smaller radii on the sky plane. The
Millennium galaxies also include realistic sub-structure as well as a
cosmological background of galaxies. While the suppression values inferred
from the two datasets are statistically the same, there are differences in
their means we might be able to constrain observationally.

At the same time, two mechanisms have been proposed to explain the suppression
of escape edge profiles: Equation 5 and Equation 8. These have fundamentally
different physical explanations and we aim to provide the first observed value
for the suppression and address the validity of the models which can explain
it.

##  3 AS1063 Phase-Space Edge Analysis

###  3.1 Inference of the Phase-Space Edge

####  3.1.1 Assignment of Projected Radii and Velocities

Figure 1:  Distribution of the 1305 redshifts for AS1063 on the sky (black
points). The green point indicates the position of the BCG
([R.A.=$22^{\text{h}}48^{\text{m}}44.0^{\text{s}}$,
Dec.$=-44^{\text{d}}31^{\text{m}}51.0^{\text{s}}$]) (Caminha et al., 2016).
These are shown relative to $r_{200}=2.63\,\text{Mpc}$ (the blue, dashed,
line), obtained from Sartoris, B. et al. (2020).

We use data for AS1063 from CLASH-VLT observations, as part of the Cluster
Lensing And Supernova survey with Hubble program (Postman et al., 2012; Lotz
et al., 2017). The galaxy spectra we use are from Karman et al. (2015);
Mercurio et al. (2021), and the full dataset contains 3850 redshifts.
Following Mercurio et al. (2021), we apply a redshift cut between
$0.30905<z<0.3716$ as a starting point for the phase-space analysis. We do not
apply any additional radial cuts, although all galaxies after this initial
redshift cut are within $r_{\perp}<5\,\text{Mpc}$. The above cuts yield 1305
galaxies remaining in the sample.

To obtain the AS1063 phase-space, we start by assigning line-of-sight peculiar
velocities to each galaxy, via

$\displaystyle v_{\text{los}}$ $\displaystyle=c\frac{z_{g}-z_{c}}{1+z_{c}},$
(9)

where $z_{g}$ is the galaxy redshift, $z_{c}$ is the mean cluster redshift,
and $c$ is the speed of light. The projected radius for each galaxy is
calculated for our chosen cosmology and the galaxy redshifts:

$\displaystyle r_{\perp}$
$\displaystyle=r_{\theta}\left(\frac{1}{1+z_{c}}\frac{c}{H_{0}}\int_{0}^{z_{g}}\frac{dz^{\prime}}{E(z^{\prime})}\right),$
(10)

where $r_{\theta}$ and $r_{\perp}$ are the angular and projected physical
separation between the galaxy and the center of the cluster, and
$E(z)=\left[\Omega_{\Lambda}+\Omega_{M}(1+z)^{3}\right]^{1/2}$ for a flat
$\Lambda\text{CDM}$ universe. We adopt redshift and spatial centers of
$z_{c}=0.3452$ and [R.A.=$22^{\text{h}}48^{\text{m}}44^{\text{s}}$ (about
$13"$ offset from the BCG), Dec.$=-44^{\text{d}}31^{\text{m}}44^{\text{s}}$]
respectively, measured from galaxies solely within the virial radius of
$r_{200}=2.4$ Mpc, given the Sartoris, B. et al. (2020) lensing estimate.

With these projected radii and velocities, we cull all galaxies with
velocities $>|4500|$ km/s as these are readily identified as non-cluster
members. Additionally, we restrict projected radii to only be only within
$0.2<r_{\perp}/R_{200}<1$, where the spectroscopic completeness of the sample
is approximately uniform (Mercurio et al., 2021), yielding 645 galaxies.

####  3.1.2 Interloper Analysis

Fadda et al. (1996a) showed that in phase-space data there can be
”interlopers” that exist outside the main body of the cluster. These are
identifiable based on velocity gaps in the data. Like many others, we use a
shifting-gapper technique (Fadda et al., 1996a; Girardi et al., 1996; Adami et
al., 1997; Wing & Blanton, 2013) to identify these interlopers in order to
infer the phase-space edge profile. We then choose a velocity gap size and a
radial bin size to conduct the shifting–gapper technique on the remaining
galaxies and identify interlopers to remove.

In deciding the value for the gap and the binning scheme, an examination of
the literature shows conflicting choices in the gap size, with 500km/s or
1000km/s being the most common choices (Fadda et al., 1996b; Wing & Blanton,
2013; Sifón et al., 2013; Gifford & Miller, 2013; Crawford et al., 2014;
Barsanti et al., 2016; Sifón et al., 2016). In terms of binning, the typical
choices is between 10 and 30 galaxies per bin. We note that Sartoris, B. et
al. (2020) used an 800km/s gap and a minimum of 15 galaxies per bin to define
interlopers in the first step of their membership measurement for AS1063.

Figure 2:  Fractional differences between edges of Millennium halos with
variable bin/GAP parameters, and fixed shifting-gapper parameters (the edge
itself is divided into 5 total bins, while the shifting-gapper parameters are
15 galaxies per bin and a velocity gap of 800 km/s, identical to Halenka et
al. (2022)). Halos are enforced to have at least 400 galaxies within
$0.3<r_{\perp}/R_{200}<1$, and the fractional differences above are averaged
over the 5/100 halos in the sample which satisfy this criterion. We measure a
statistically negligible difference ($\sim 0.5\%$) between the variable and
fixed bin/GAP edge definitions, indicating that the broadening of the
interloper range shown as opposed to fixed parameters does not induce a source
of bias on the edge relative to Halenka et al. (2022)

Figure 3:  Left panel: Predicted radial suppression of a AS1063-type cluster
using the AGAMA analytic model. A smaller bin-width corresponds to lower
sampling per bin and hence more suppression, as shown by the fact that $Z_{v}$
is higher for 10 bins than 5 bins. Right panel: Example analytic phase-space
of a AS1063-type cluster using AGAMA, where the edge is again measured using 5
(red regions) and 10 bins (blue regions). It is clear that a lower bin width
yields more suppression, however the same 3D escape profile is still recovered
(to within uncertainty) so long as the correct bin calibration is used. Figure
4:  The phase-space edge for Abell S1063 inferred from the galaxy data after
radial and velocity cuts from $0.2\,r_{200}<r_{\perp}<r_{200}$ (the range
where the suppression function is used, and also spectroscopic completeness is
approximately uniform). Here, the edge incorporates the scatter (red bands)
from interloper variation.

Different bin sizes and velocity gaps will often identify different
interlopers, and there is no obvious “best” a priori parameter choice. Old et
al. (2014) compared two values for the gap size in clusters from an N-body
simulation, and found that the mass inferred using the velocity dispersion
using the smaller gap size had a smaller mass than the cluster with the larger
gap. The clusters with the smaller gap size were also more discrepant from the
true $M_{200}$ compared to the large gap size. Crawford et al. (2014)
performed a systematic study by varying the number of galaxies per bin as well
as the gap size. They found a weak trend of an increasing mass with higher
galaxies per bin, and a stronger trend of increasing mass with higher velocity
gap size.

Similar to Crawford et al. (2014), we vary the bin size and velocity gap to
identify interlopers and compare the measurement to a fiducial one. While they
studied the dispersion, we study the edge suppression ($Z_{v}$). For a
baseline constraint on accuracy and precision, we use clusters from N-body
simulations where the true edge and suppression is known. We use the halo
sample from Halenka et al. (2022) which was based on semi-analytic galaxies in
the Millennium simulation. We use their galaxy data having the dimmer apparent
magnitude limit for semi-analytic galaxies to maximize the phase-space
sampling. We then choose the five halos with the highest sampling and similar
to AS1063 with $N_{\text{{gal}}}>500$ within $0.2<r_{\perp}/R_{200}<1$. This
exercise could be repeated for halos with lower sampling, but that is beyond
our scope.

Using projected phase-spaces, we then measure the escape edges over a multiple
lines-of-sight for a fixed set of shifting-gapper parameters. We use the same
parameters of 800 km/s and 15 galaxies/bin from Halenka et al. (2022). Using
the underlying 3D edges, we recover the radially averaged suppression from
Halenka et al. (2022) which we label as $Z_{v}^{\rm fixed}$. We then measure
the edge over a wide range of shifting-gapper parameters and calculate new
suppression values, one for each of the $(i,j)$ pairs of bin and gap sizes
ranging from 10-50 galaxies/bin and a velocity gap ranging from 100-1500 km/s.
The accuracy of these edges can then be assessed by calculating the mean and
standard deviation of the absolute value of the fractional difference,
$|f(\Delta Z_{v})|$, between $Z_{v}^{\rm{fid}}$ and $Z_{v}^{i,j}$.

These fractional differences are shown in Figure 2. There is a clear high
variance region below 20 galaxies/bin. Above a velocity gap of $\sim 1000$
km/s and $\sim 35$ galaxies/bin, we see convergence to $|f(\Delta Z_{v})|\sim
0.02$, i.e. that the edges become biased relative to the fixed case. Given a
redshift uncertainty of 150 km/s (Mercurio et al., 2021), we identify a
velocity gap floor of twice this uncertainty or 300 km/s. The upper bound of
the velocity gap is chosen to keep $|f(\Delta Z_{v})|$ low, where 800 km/s is
the maximum gap where $|f(\Delta Z_{v})|<1\%$. This corresponding parameter
space where the edge measurement is robust to velocity gap and bin variation
is shown by the blue box in Figure 2, where we measure $|f(\Delta
Z_{v})|=0.005\pm{0.003}$.

For a system with observed data like AS1063, we have used simulation halos to
show that, for a wide range of gap-size/bin choices, the systematics induced
into the inferred 3D edge from the projected data are at the sub-percent
level. We also use the simulations to constrain the measured edge scatter,
which we carry through the rest of the analyses.

###  3.2 Binning and Smoothing the Edge Profile

In order to calculate the phase-space maxima, we need to choose a bin size.
Typically, a nominal physical bin size is chosen, e.g., 0.2 Mpc (Stark et al.,
2019). However, prior characterizations of cluster escape profiles have used
data with significantly lower sampling. In this section, we explore the effect
of bin size to make an informed as opposed to arbitrary choice for the bin
size.

We explore the bin size using the Action-based Galaxy Modeling Architecture
(AGAMA) framework (Vasiliev, 2019). Distinct from N-body simulations such as
Millennium, This technique relies on forward modeling a cluster using
Action/Angle variables (Binney & Tremaine, 2008). Unlike a Jean’s inversion
scheme (e.g. Sartoris, B. et al. (2020)), AGAMA uses a set of predefined
characteristics of a cluster (e.g. the mass, density profile, velocity
anisotropy, cosmology, etc.) in order to see how these affect the phase-space
characteristics. In our analysis of the analytic halos, we use the same mass
(assuming a Dehnen density profile) and sampling ($\sim 800$ tracers within
$r_{200}$) as AS1063 (although note that Halenka et al. (2022) found the
suppression is not affected by the mass of the system). AGAMA is preferable to
Millennium in this context given the systematic control of parameters it
allows us.

To incorporate the effects of cosmology, we cull all tracers outside the
cosmological escape profile (Equation 3) for the pre-specified cosmology in
§1. As mentioned in §1, AGAMA however lacks a cosmological background as well
as realistic substructure. We trade these effects however for systematic
control of phase-space characteristics. In the projected analytic phase-
spaces, we follow the same projection procedure from Gifford & Miller (2013);
Halenka et al. (2022) and treat the viewing angle as a random variable along
the z-axis. We simulate 100 different realizations of phase-spaces for each
mock AGAMA cluster. The clusters are placed at the same redshift as AS1063,
for which we sample 100 different viewing angles.

In the left panel of Figure 3, we show how the phase-space edge profile
differs for different bin choices, corresponding to 5 bins (blue shaded
region) and 10 bins (red shaded region). Here, the suppression is measured for
the phase-space edge between $0.2\,r_{200}<r_{\perp}<r_{200}$ where 100
different line-of-sight draws are performed. It is apparent that the higher
bin count yields more suppression at all radii, which can be attributed to the
fact that the average sampling per bin is reduced for a smaller bin width.

In the right panel of Figure 3, we show the effect on the edges from the
number of bins. This is for a single line-of-sight draw on a single projected
AGAMA phase-space. The edge from 10 bins is more suppressed than the edge from
5 bins. However, so long as this increased suppression is accounted for using
the analogous $Z_{v}$ (e.g., from the left panel of Figure 3), the 3D edge
profile will still be inferred. The final accuracy and precision on the edge
from binning is contained within the $Z_{v}$ function. In addition to the
scatter from gapsize/bin-size in Section 3.1.2, we incorporate the suppression
uncertainty in the mass inference.

For AS1063, we will use 5 bins for the rest of this paper and the appropriate
$Z_{v}$ from Figure 3. Furthermore, we include the radial shape dependence of
the suppression function, which is higher in the core compared to $r_{200}$.
Specifically, we have used the AGAMA framework to calibrate $Z_{v}(r)$ for 5
bins and for a sample of galaxies selected within the range
$0.2\,r_{200}<r_{\perp}<r_{200}$ where the spectroscopic completeness function
is approximately uniform from (Mercurio et al., 2021).

While potentials are smooth on the scales probed here, the measured edges can
be noisy. There are various techniques used to control this effect (Serra et
al., 2011; Gifford et al., 2013). In our work, we employ smoothing functions,
which of course relates closely to binning. Our smoothing analysis requires us
to choose the number of bins and also a form for the parameterized function to
smooth the raw data. We then measure the statistical bias and variance. In
this context, the square of the bias of the model quantifies the degree to
which the model is overfit, while the variance of the model quantifies the
degree to which the model is underfit. If these quantities are equal, or their
sum is minimized, the fit will be optimized, and can be quantified according
to the statistical risk, or mean square error (MSE), of the model, given by
(Keener, 2010)

risk
$\displaystyle=(\theta-\mathbb{E}[\delta])^{2}+\mathbb{E}[(\delta-\mathbb{E}[\delta])^{2}],$
(11)

Here, $\mathbb{E}[x]$ is the expected value of the variable $x$, $\theta$ is
the escape edge, and $\delta$ is our estimator we compare the model to, or an
NFW fit of the escape profile, obtained through MCMC (see §4).

To measure the MSE, we need to choose a function to make the estimator
$\delta$. We explore a range of different smoothing functions, with polynomial
degrees from one through five (linear, quadratic, etc.). We find quadratic and
cubic smoothing functions minimize the MSE by $\sim 70\%$ relative to the
unsmoothed edge and higher order polynomials. Hence, we choose to use these
fits to avoid over and under-fitting of the edge profile222It is worth noting
that the virial mass is statistically unaffected by the choice of smoothing.
Rather, smoothing has the strongest effect on the edge profile near the core,
which primarily controls the measured concentration..

##  4 Mass Estimation

###  4.1 From the Phase-Space Edge to $M_{200}$

The phase-space edge for AS1063 is shown in Figure 4, where the errorbars
combine interloper and redshift error. To model the effective potential
profile from this edge, we use the NFW potential (Miller et al., 2016)

$\displaystyle\Psi(r)$ $\displaystyle=-\frac{4\pi
G\rho_{s}r_{s}^{2}\log(1+r/r_{s})}{r/r_{s}},$ (12)

where $r_{s}$ is the scale radius and $\rho_{s}$ is the normalization. However
in modeling the effective potential profile, we need to account for the phase-
space edge suppression due to sampling, and also over and under-fitting. We
use the $Z_{v}(r)$ in Figure 3 to account for suppression, and a quadratic
smoothing function to account for over and under-fitting.

Given the 3D phase-space edge (including its error from suppression), we use
equations 3 and 8 to infer the mass profile using the Poisson equation. We use
Bayes’ theorem and Goodman & Weare’s affine invariant Markov chain Monte Carlo
(MCMC) Ensemble sampler to model the posteriors of the NFW density profile
(Foreman-Mackey et al., 2013), corresponding. In all cases we use a Gaussian
likelihood, where the error comes from the phase-space edge uncertainty. To
find $M_{200}$, we interpolate the cumulative mass density to identify the
radius which is 200$\times$ the critical (and also mean) density and report
the corresponding mass within that radius as $M_{200}$.

We consider two fitting options to the NFW: a single parameter model using two
different mass-concentration relations and a two parameter model fitting to
the central density and scale radius. The priors are summarized in Table 3.
Here, we present only the NFW profile, while the other profiles are presented
in the Appendix.

In the top left panel of Figure 5, we show the 2D posterior for the NFW
profile’s parameters. We use the GTC package (Bocquet & Carter, 2016) to
generate the histograms. The covariance of $r_{s}$/$\rho_{s}$ suggests that
only a slim range of these parameters yield physical mass constraints, both in
cases when the edge profile is smoothed and unsmoothed. The top right panel
converts our $r_{s}$/$\rho_{s}$ posterior to $M_{200}/c_{200}$, measured with
respect to the mean density. We find that when smoothing the edge profile, our
$M_{200}/c_{200}$ is consistent with Gruen et al. (2013), shown by the blue
contours. If no smoothing functions are used, the mass profile is relatively
unchanged, but the concentration increases by $\sim 0.5\sigma$. Further, the
covariance in $M_{200}/c_{200}$ yields a different shape when the escape
velocity is used as opposed to lensing, where we find the mass varies hardly
at all with the concentration. This is expected from the fact that the
potential profile is a second derivative of the density profile and therefore
much flatter so that it is less well constrained by the data. This analysis
suggests that future joint analyses using lensing and the escape profiles
could place tight constraints on the mass-concentration relation.

In the bottom panel of Figure 5, we show the escape profile inferred from the
smoothed two-parameter (M200, c200) model, with respect to the critical
density. The agreement with Gruen et al. (2013) is highlighted by the
agreement in our escape profile and their inferred escape profile, given their
$M_{200}/c_{200}$ posterior. While more recent lensing estimates are available
(see Table 1), we use Gruen et al. (2013) as a baseline given their published
mass/concentration posterior contours.

The results (16.5,50,83.5 percentiles of the 1D marginalized posteriors) for
all of our the NFW fits (in addition to the Einasto and Dehnen fits from the
Appendix) are provided in Table 2.

Figure 5:  Left panel: 2D posterior estimates (68% and 95%) of the NFW profile parameters. The covariance of the distribution suggests that only a narrow range of $r_{s}$ and $\rho_{s}$ values yield physical mass constraints. Right panel: Conversion of the $r_{s}$ and $\rho_{s}$ posterior to the concentration and mass, measured with respect to the mean density of the universe. We show the posteriors for both the smoothed and unsmoothed edges, where it is apparent smoothing primarily controls the concentration. The blue contours are taken from the weak-lensing analysis of Gruen et al. (2013). The covariance in mass/concentration has a different shape than that from weak lensing, where our mass posterior is nearly independent of concentration. Right panel: The best-fit two parameter (M200, c200) NFW escape profile relative to the Gruen et al. (2013) escape profile, inferred from their mass. These profiles are shown relative to all the galaxy data (including interlopers). Table 2: AS1063 Potential Profile and Mass Constraints. Central values are medians, while errors are 16.5% and 83.5% percentiles of the marginalized 1D posteriors.  NFW 1 parameter |  |
---|---|---
Edge | log10 M200 |  |  |  |  |  | Notes
unsmoothed | 15.39${}^{+0.03}_{-0.04}$ |  |  |  |  |  | Duffy et al. (2008b)
smoothed | 15.40${}^{+0.03}_{-0.03}$ |  |  |  |  |  | -
unsmoothed | 15.59${}^{+0.04}_{-0.04}$ |  |  |  |  |  |  $200\times\rho_{{\rm mean}}$)
smoothed | 15.60${}^{+0.04}_{-0.04}$ |  |  |  |  |  | -
NFW 2 parameter |  |
Edge | log10 M200 | c200 | r200 | log10 $\rho_{s}$ | log${}_{\rm 10}r_{s}$ |  | Notes
unsmoothed | 15.42${}^{+0.06}_{{\color[rgb]{0,0,0}{-0.09}}}$ | 3.53${}^{+3.81}_{-2.43}$ | 2.50${}^{+0.10}_{-0.16}$ | 14.89${}^{+0.73}_{-0.99}$ | -0.14${}^{+0.48}_{-0.31}$ |  | fit to $\rho_{s}$, rs
smoothed | 15.40${}^{+0.06}_{-0.12}$ | 1.91${}^{+2.8}_{-1.31}$ | 2.46${}^{+0.12}_{-0.2}$ | 14.34${}^{+0.86}_{-0.88}$ | -0.12${}^{+0.46}_{-0.39}$ | - |
unsmoothed | 15.54${}^{{\color[rgb]{0,0,0}{+0.04}}}_{-0.04}$ | 4.94${}^{+5.07}_{-3.28}$ | 3.52${}^{+0.10}_{-0.10}$ |  |  |  | $200\times\rho_{{\rm mean}}$
smoothed | 15.55${}^{{\color[rgb]{0,0,0}{+0.04}}}_{{\color[rgb]{0,0,0}{-0.04}}}$ | 2.77${}^{+3.75}_{-1.81}$ | 3.57${}^{+0.10}_{-0.10}$ |  |  |  | -
Einasto 3 parameter
Edge | log10 M200 | c${}_{200}(r_{-2})$ | r200 | log10 $\rho_{0}$ | log10h | n | Notes
unsmoothed | 15.40${}^{+0.08}_{-0.16}$ | 3.79${}^{+5.15}_{-2.67}$ | 2.46${}^{+0.15}_{-0.29}$ | 19.25${}^{+1.62}_{-1.85}$ | -6.42 ${}^{+2.85}_{-2.23}$ | 5.82${}^{+1.55}_{-1.99}$ | fit to $\rho_{0}$, h, n
smoothed | 15.34${}^{+0.11}_{-0.23}$ | 1.95${}^{+2.86}_{-1.55}$ | 2.35${}^{+0.21}_{-0.38}$ | 18.81${}^{+1.58}_{-1.72}$ | -6.50${}^{+2.82}_{-2.37}$ | 6.11${}^{+1.55}_{{\color[rgb]{0,0,0}{-2.03}}}$ | -
unsmoothed | $15.54^{+0.08}_{-0.16}$ | $3.79^{+5.15}_{-2.67}$ | $2.46^{+0.15}_{-0.29}$ |  |  |  | $200\times\rho_{{\rm mean}}$
smoothed | 15.49${}^{+0.07}_{{\color[rgb]{0,0,0}{-0.12}}}$ | $2.81^{+3.83}_{-2.17}$ | $3.38^{+0.17}_{-0.30}$ |  |  |  | -
Dehnen 3 parameter
Edge | log10 M200 |  | log10 r200 | log10 $M_{{\rm tot}}$ | rs | $\gamma$ | Notes
unsmoothed | 15.43${}^{+0.05}_{{\color[rgb]{0,0,0}{-0.10}}}$ |  | 2.52${}^{{\color[rgb]{0,0,0}{+0.10}}}_{-0.10}$ | $15.88^{+0.25}_{-0.18}$ | 17.27 ${}^{+8.63}_{-9.68}$ | 2.09${}^{+0.12}_{-0.19}$ | fit to Mtot, rs, $\gamma$
smoothed | 15.39${}^{+0.08}_{{\color[rgb]{0,0,0}{-0.17}}}$ |  | 2.45${}^{{\color[rgb]{0,0,0}{+0.08}}}_{-0.17}$ | $16.03^{+0.28}_{-0.20}$ | 14.76 ${}^{+9.99}_{-9.22}$ | 1.74${}^{+0.19}_{-0.36}$ | -
unsmoothed | 15.51${}^{{\color[rgb]{0,0,0}{+0.04}}}_{-0.05}$ |  | 3.44${}^{{\color[rgb]{0,0,0}{+0.11}}}_{{\color[rgb]{0,0,0}{-0.13}}}$ |  |  |  | $200\times\rho_{{\rm mean}}$
smoothed | $15.51^{+0.04}_{-0.08}$ |  | $3.54^{+0.13}_{-0.18}$ |  |  |  | -

###  4.2 Dependence on the Density Profile

Miller et al. (2016) noticed that when constraining against the NFW density
profile of dark matter halos, the inferred escape velocity was higher than
what was actually observed in simulations. They attribute this to the fact
that the NFW density profile is too shallow in the outskirts for realistic DM
halos (see also Diemand et al. (2005)). In our results, we are not
constraining to the density, but instead to the escape edge. By definition,
the NFW, Einasto, and Dehnen potential profiles will be nearly identical to
each other within the fitting errors, due to the fact that they are being fit
to the same edge.

Using the best-fit parameters, we plot the fractional difference in the
density profiles in Figure 6. We see percent level agreement out to $r_{200}$.
Beyond that radius, we also see the Einasto and Dehnen densities become much
smaller compared to the NFW, expected from Miller et al. (2016). The key
points of this figure are that: (a) the dynamically inferred density profiles
are nearly identical within $r_{200}$ for the three functional forms we
examine and (b) that the integral to $r_{\text{{eq}}}$ of the density in
Equation 4 will sum to a larger value in the NFW compared to the integrals for
the Einasto and Dehnen forms. This is reflected in the estimated $M_{200}$
values in Table 2, although the differences are not statistically significant.
This analysis provides a promising new path forward in determining the true
relative steepness of the densities in the outer regions of galaxy clusters,
which is vital in studies of the splashback radius (Contigiani et al., 2019).

Figure 6:  The fraction of the Einasto or Dehnen density with respect to the
NFW functional form for our best parameters in Table 2. While there is close
agreement inside $r_{200}$, the difference in the outer radii reflects the
fact that both the Einasto and Dehnen are steeper. This leads to a smaller
expected mass from the escape velocity profile.

##  5 Velocity Dispersion-Based Mass

One of the main techniques used to infer dynamical masses is through the
velocity dispersion. Gómez et al. (2012) estimated the velocity dispersion of
AS1063 to be $\sigma_{v}=1840^{+230}_{-150}$ km/s within $\sim 730$ kpc. Due
to the fact that they observed only a small fraction of the virial radius,
they suggested a correction to this measurement out to $r_{200}$ based on an
isothermal King model. Their corrected dispersion was reported as
$\sigma_{v}=1660^{+230}_{-150}$ km/s, where they used the Evrard et al. (2008)
mass vs. $\sigma$ relation for dark matter to estimate a mass significantly
higher than lensing and SZ estimates, which we report in Table 1.

More recently, Mercurio et al. (2021) estimated $\sigma_{v}$ using the same
data we use in this work, and found $\sigma_{v}=1380^{+26}_{-32}$ km/s for the
whole observational sample. However, this estimate cannot be easily translated
into a virial mass given that the spectroscopic completeness varies
significantly over the entire radial range, therefore making this dispersion
difficult to interpret.

A more sophisticated dynamical approach was taken by Sartoris, B. et al.
(2020), who performed a reconstruction of the mass profile through a Jean’s
likelihood regression analysis for the Jean’s equation of dynamical
equilibrium using the MAMPOSSt framework (Mamon et al., 2013). From their mass
profile, their virial mass is smaller than the Gómez et al. (2012) dispersion
estimate, but both are systematically still larger than weak lensing and SZ
measurements in Table 1. Moreover, the Jean’s likelihood analysis implicitly
assumes dynamical equilibrium (Mamon et al., 2013), while AS1063 has been
shown to have a notable presence of substructure (Gómez et al., 2012; Mercurio
et al., 2021). This motivates a further exploration of the dispersion-inferred
mass for AS1063.

###  5.1 The Dispersion of AS1063

We start by calculating the dispersion of AS1063, in this case following the
same approach we used to estimate the phase-space edge profile by
incorporating interloper stochasticity. Previous dispersion estimates (Gómez
et al., 2012; Mercurio et al., 2021) of the system used a fixed bin/GAP pair
to identify interlopers. The procedure is analogous to §4, except in this case
we include data within $0.2\,r_{200}$ for the dispersion estimate, although no
interlopers are identified within $0.2\,r_{200}$.

We provide an estimate for the dispersion using standard bi-weighted estimator
(Beers et al., 1990), of $\sigma_{v}=1477^{+31}_{-30}$ km/s. The uncertainty
on the dispersion is quantified using a bootstrapping technique (Beers et al.,
1990). In our case, we use a Monte Carlo approach. First, we draw a new subset
of galaxy velocities from the original data using normal distributions defined
by the mean velocity of each galaxy and a fixed uncertainty of
$150\text{km/s}$ (Mercurio et al., 2021). We then bootstrap resample the data
to create a new distribution of the velocities which accounts for the redshift
measurement uncertainties.

We note that the bi-weight estimator requires a tuning parameter, c which
defines the weights on the outliers. The value of the tuning parameter
corresponds to the multiple of the median absolute deviation of the data. We
examined the dispersion as a function of c and find that it depends strongly
on c below c = 25. A high tuning constant is a better choice when few outliers
are expected. Given that we employ the shifting-gapper technique to remove
interlopers prior to calculating the dispersion, we argue that the tuning
parameter should be defined with a value where it does not drive the actual
value of the measured dispersion. Therefore, we choose c = 25. Values as low
as c = 15 or as high as c = 30 change the measured dispersion by less than
$0.3\%$ from its measurement for c = 25. However, the change in the dispersion
is tens times that level (towards higher dispersion) going from c = 15 to c =
7. When using the default tuning constant used in Astropy (Astropy
Collaboration et al., 2013, 2018, 2022), of c = 9 , not accounting for proper
weighting of outliers, we recover the same dispersion for the whole sample
($\sigma_{v}=1380^{+26}_{-32}$ km/s) as Mercurio et al. (2021).

###  5.2 From $\sigma_{{\rm v,los}}$ to $\sigma_{{\rm v}}$

Figure 7:  Left panel: Predicted virial scaling relation in AGAMA as a
function of velocity anisotropy for $\beta=0.5$ (orange line), $\beta=0$ (blue
line), and $\beta=-0.5$ (green line), where the small errorbars come from
variation in phase-spaces for the $N=10^{5}$ particles. We find $\sim 5\%$
variation between the scalings between $\beta=0.5$ and $\beta=-0.5$. Right
panel: Ratio between the line-of-sight velocity dispersion and 1D velocity
dispersion for $N=10^{3}$ particles. We find that when sampling along the
line-of-sight, anisotropy and cluster mass have little effect, with
$\sigma_{v,{\rm los}}$ recovering $\sigma_{v}$ in all cases (at higher N, all
$\sigma_{v,{\rm los}}/\sigma_{v}$ converge to exactly 1). At this choice of
$N$, we expect variation in this ratio of $4\%$, corresponding to an
uncertainty in the line-of-sight dispersion of $\sim 60$ km/s for AS1063.
Figure 8:  Comparison of previous mass estimates in the literature from Table
1 for lensing (cyan points), SZ (black points), X-ray (green points), and
dynamical masses (lime green points). We note these do not each have measured
dispersions, but are plotted with values close to our $\sigma_{v}$ for
convenience. These are shown compared to the mass constraint of this work
using the escape velocity, shown by the red shaded region. We also compare
these to theoretical predictions from simulations from Anbajagane et al.
(2022). We lastly show our velocity dispersion, measured from (1) The
uncertainty determined through the $150$ km/s redshift error (Mercurio et al.,
2021) determined through bootstrap simulations (2) The uncertainty in using a
line-of-sight dispersion given our sampling of $N=828$ tracers to measure the
dispersion, using the projection model of this work (3) The uncertainty due to
$\beta$, assuming $0<\beta<0.5$, yielding an additional contribution based on
the model from this work.

Next, we perform dynamical modeling of cluster phase-spaces to investigate
their more general properties. Our approach for dynamical modeling relies on
the Action-based Galaxy Modeling Architecture (AGAMA) framework (Vasiliev,
2019) (see §3.2). AGAMA can be contrasted with the Jean’s inversion performed
in Sartoris, B. et al. (2020), as AGAMA uses a predefined set of
characteristics such as the system’s mass, the density profile, the velocity
anisotropy, etc. to infer the phase-space properties, rather than trying to
infer these properties from the phase-space directly. In our case, we vary
these parameters to see how they affect predictions of the velocity
dispersion.

For a given M200, we generate mock phase-space tracers at $z=0$. We use a
variable number of tracers, from a few hundred to $N=10^{5}$. We specify a
Dehnen density profile which is similar in shape to AS1063. We then define a
log-linear grid of masses between $10^{14}$ and $10^{15.5}$ M⊙, and also a
coarse grid of constant velocity anisotropies between $\beta=-0.5$ and
$\beta=0.5$ (Lemze et al., 2012; Stark et al., 2019). For the generated mock
phase-space, the 1D velocity dispersion
${\color[rgb]{0,0,0}{(}}\sigma_{\text{3D}}/\sqrt{3}$) is inferred for the
particles within $r_{200}$.

We first compare the AGAMA 1D dispersions to the simulation prediction for
dark matter in halos from the Bahamas + Mascis Simulation (McCarthy et al.,
2016; Barnes et al., 2017) and The Three Hundred Project (Cui et al., 2018),
using the publicly available fitting functions from Anbajagane et al. (2022).
Three major differences between AGAMA and the N-body simulations are the non-
linear growth patterns, realistic sub-structures, and a cosmological
background, which all exist in the simulation halos. We trade this for the
ability to systematically control the velocity anisotropy of the mock clusters
in AGAMA. This would be difficult in the simulations, where there is a finite
(and often small) number of available halos which have galaxies with the right
mixtures of orbits and dynamical populations.

In the left panel of Figure 7, we show the predicted $\sigma_{v}$ vs.
$M_{200}$ scalings of the mock AGAMA phase-spaces when using $N=10^{5}$
tracers, which is similar to the number of dark matter particles in the
simulations. Overall, we note three interesting findings: (1) the simulation
particles have higher dispersions at fixed mass compared to the AGAMA tracers;
(2) the slopes between the simulations and the AGAMA tracers differ; and (3)
we find an unexpected dependence for the 1D dispersion by the anisotropy.
Explaining these results is beyond the scope of this work. For now, we simply
see how $\beta$ affects the 1D velocity dispersion at a given $M_{200}$.

In the right-hand panel of Figure 7 we show how well we can recover the 1D
dispersion ($\sigma_{v}$) using only projected velocities for cluster with the
same number of tracers (828 between $0.0\leq r_{\perp}/r_{200}\leq 1$) as our
dataset for AS1063. We find that our line-of-sight dispersion is unbiased as
both a function of mass and anisotropy.

###  5.3 $\sigma_{v}$ and spectroscopic completeness

As noted earlier, Mercurio et al. (2021) report a non-uniform sampling
completeness function corresponding roughly to: $100\%$ from $0\,r_{200}$ to
$0.1\,r_{200}$, $80\%$ from $0.1\,r_{200}$ to $0.5\,r_{200}$, and $70\%$ from
$0.5\,r_{200}$ to $r_{200}$. A benefit of AGAMA is that we can create mock
phase-spaces using this sample spectroscopic completeness function in order to
access the magnitude of its affect on the dispersion.

The phase-spaces we construct are analogous to AS1063, where we use $N=828$
particles within the virial radius, and estimate the line-of-sight dispersion
using 100 line-of-sight draws from the specified completeness function. We
find the statistical uncertainty from the phase-space sampling (see §5.4)
dominates over the effect of non-uniform completeness, where we find a
difference of $4\pm 34$ km/s compared to a full completeness function for a
$\beta=0.5$ cluster, where differences in $\beta$ are negligible compared to
the statistical uncertainty. Due to the fact that the statistical uncertainty
dominates over the completeness uncertainty, we neglect the effect of non-
uniform sampling on the dispersion analysis.

###  5.4 $\sigma_{v}$ and its uncertainties

We made an initial estimate of the dispersion in §3.1, where we found
$\sigma_{v}=1477^{+31}_{-30}$ km/s for galaxies within the virial radius. As
noted before, we choose a tuning constant (Beers et al., 1990) of $c=25$,
corresponding to $25\times$ the median absolute deviation from the median
residual in the data. This high tuning constant results in a stable
measurement of the dispersion (see §5.1). The errors are calculated as before
using bootstrap resampling and accounting for galaxy redshift errors of 150
km/s. As shown in the last section, we expect $\sigma_{v,{\rm los}}$ to be an
unbiased measure of the 1D velocity dispersion, $\sigma_{v}$.

We consider two additional effects which introduce uncertainties into our
determination of $\sigma_{v}$: sampling, and radial dependence. We calibrate a
model for how sampling affects the projection. This has been previously
studied in simulations by other authors (Gifford & Miller (2013); Serra et al.
(2011); Caldwell et al. (2016). We find that the statistical scatter as a
function of phase-space sampling follows the same simple power-law obtained by
Caldwell et al. (2016) who used simple Gaussians as the 1D representatives of
the phase-space data. We measure a finite sampling uncertainty of
${}^{+43}_{-43}$ km/s from given our sample size of $N=828$ within r200. The
size of this error is reflected in the right panel of Figure 7.

Since $\sigma_{v}$ is only measured within $r_{200}$, the uncertainty in the
inferred virial radius presents an additional source of error. When
constructing the mock line-of-sight dispersion profile in AGAMA, we find no
considerable differences in $\sigma_{v}(<r)$ in the radial profile between
$0.2\,r_{200}$ and $r_{200}$. However, in the data we find differences of
${}^{+21}_{-34}$ km/s between the inner and outer 1$\sigma$ ranges on r200. We
include these as a systematic uncertainty.

Our final dispersion uncertainty on the inferred 1D $\sigma_{v}$ from the
observed $\sigma_{v,{\rm los}}$ includes (1) the redshift measurement
uncertainty; (2) the systematic error from the uncertainty in $r_{200}$; (3)
the finite sampling uncertainty; and (4) the uncertainty in theory from
anisotropy from Figure 7 (left). The final dispersion we obtain is then
$\sigma_{v}=1477^{+87}_{-99}$ km/s.

In Figure 8, we plot the escape mass using the analytic suppression model,
inferred by the NFW density profile (vertical red bands) against the inferred
1D velocity dispersion (horizontal orange bands). If we used the Millennium
suppression model, we infer more suppression corresponding to lower masses
by$\leavevmode\nobreak\ \sim 0.1$ dex. The colored squares are mass values
from the literature (see Table 1). The literature data are placed at the
position of the central value of $\sigma_{v}$ and offset for clarity. The
diagonal bands are expectations for galaxies with stellar masses greater than
$5\times 10^{9}M_{\odot}$ from two hydrodynamic simulations (Anbajagane et
al., 2022).

When compared to masses inferred from weak (and strong) lensing, the SZ
effect, and the X-ray temperature, we find that $\sigma_{v}$ is consistent
with the expectations from the cosmological hydrodynamic simulations so long
as the escape velocity mass of this work is used (note again the masses are
not directly associated with a $\sigma_{v}$, but are simply plotted for
convenience in a band similar to our own), with most lensing estimates being
within $1\sigma$. The Jean’s estimate is also in good agreement with the
theoretical prediction, while only one of the SZ masses appears to be. The
X-ray masses are slightly above the theoretical expectations, and adds to the
evidence that the system is not in hydrostatic equilibrium.

##  6 Discussion

Our goal in this work is to infer an accurate and precise mass for AS1063
using the escape velocity. As a baseline, we focus our comparisons on non-
dynamical lensing inferred masses. There are other studies comparing masses
(or mass profiles) inferred from the “caustic” technique to weak lensing
(Diaferio et al., 2005; Geller et al., 2013; Rines et al., 2013; Hoekstra et
al., 2015; Herbonnet et al., 2020). However, this is the first attempt at
measuring the mass of a galaxy cluster using the escape edge as a tool and as
derived in Nandra et al. (2012), Behroozi et al. (2013), Miller et al. (2016),
and Halenka et al. (2022). This technique is fundamentally different than the
“caustic” technique (Diaferio & Geller, 1997; Gifford et al., 2013) for three
reasons: (1) cosmology is a required component in our analysis; (2) the
suppression of the edge of the data compared to dark matter simulations is
assumed to be caused from statistical sampling and not from the velocity
anisotropy (Halenka et al., 2022); (3) the velocity dispersion is not used to
calibrate the edge.

In our analyses, we have been careful to study the uncertainties on the escape
mass, which include not only measurement error on the galaxy redshifts, but
also edge error from the interloper removal process and uncertainty in the
radial suppression. We have paid similar attention to the error in the
inferred 1D velocity dispersion.

In terms of precision, we have shown that this new technique provides similar
or better mass constraints compared to current weak lensing data.
Surprisingly, the covariance between M200 and concentration has both a
different shape and also a lesser effect on the mass than in weak lensing
measurements (see Figure 5). We suggest that this lack of sensitivity in the
constraint of the shape of the potential is from the fact that the escape
profile is the square root of the potential profile, which is naturally
flatter than the steep density profile. This opens up a new possibility to
place more precise constraints on the mass-concentration relation through
joint shear/v${}_{\text{esc}}$ posterior constraints on mass and shape. This
idea might also enable better theoretical predictions on the M200-c relation,
since the potential and density profiles of dark matter particles in
simulations are independent measures of the underlying shape of the Poisson-
pair profile.

In order to infer the mass from the edge, we use a statistical suppression
model based on the work of Halenka et al. (2022). The value of this
suppression can be measured using an analytic model for the phase-space or
N-body simulations. In terms of the radially averaged value, the predictions
are similar ($\frac{1}{Z_{v}}=0.78$ and $\frac{1}{Z_{v}}=0.84$ for the
analytic and N-body models respectively). This is an additional systematic
that needs further evaluation. Our reported masses in Table 2 and all of the
figures use the AGAMA calibration. Had we instead used the Millennium-based
suppression, our final mass for AS1063 would be lower by about 0.1dex and
close aligned to the weak-lensing and SZ data in Figure 8.

Prior work meanwhile (see §2.2) had explained the edge suppression through
velocity anisotropy via $g(\beta)$-2 in Equation 5. Thus, we have that
$g(\beta)$${}^{-2}\equiv$ $\frac{1}{Z_{v}}=0.78$${}^{+0.05}_{-0.04}$. However,
a nominal value of $\beta=0.25$ would suppress the observed edge by 0.55,
which is much more than what we observe. Pure isotropy would suppress it a
little less to 0.58. A fully rotating cluster (i.e. with no radial motion of
its members) would suppress the edge to 70% of its 3D value, which is closer
to what we observer, but is unphysical. The difference is exacerbated with the
Millennium-based suppression value. From this we conclude that the explanation
for suppression of the observed edge due to velocity anisotropy is
inconsistent with both the data and the theory.

For Abell S1063, we have shown that it is possible to its dynamical mass
estimate for a non-relaxed system in-line with weak lensing and SZ mass
estimates, as well as with theoretical predictions from the Three-Hundred
Simulation and Bahamas + Mascis hydrodynamic cosmological simulations
(McCarthy et al., 2016; Barnes et al., 2017; Cui et al., 2018)

On a final note, even though the weak lensing and escape masses are consistent
with each other, they may still be individually biased. Consider the weak
lensing mass estimate of Umetsu et al. (2016). Besides the statistical
uncertainty on the shear measurement, the total cluster mass uncertainties
incorporate systematic errors due to the mass-sheet uncertainty, uncorrelated
large-scale structure, halo asphericity, the presence of correlated halos,
model-dependent systematics in the strong lensing, as well as systematics from
their radial binning scheme. While an exhaustive list, there is the
possibility that unknown systematics remain, such as internal cluster
substructure. Gruen et al. (2013) identified a second density peak within 1
Mpc of the center and regressed against a mass model which contained two NFW
halos. However, their combined two halo mass is statistically consistent with
their single NFW model fit.

For the escape mass, the mass uncertainty incorporates statistical
uncertainties on the galaxy redshifts, systematic uncertainties from the
shifting-gapper parameters, intrinsic scatter in the suppression function, and
the effect of realistic phase-space sub-structure on the suppression function.
Like lensing, it is possible that there remain unknown systematics, such as
internal substructure. However, Behroozi et al. (2013) used simulations to
show that the escape velocity has negligible correlation with the presence of
substructure. Therefore, other than radial binning, there are no clear
systematics shared in common between weak lensing and the escape techniques.
To fully quantify the agreement between an escape-based dynamical mass and a
shear-based non-dynamical mass, we need to take the next steps to increase the
sample size and we need clusters with quality data for both the lensing and
the spectroscopy.

##  7 Acknowledgments

This material is based upon work supported by the National Science Foundation
under grant No. 1812739.

## References

  * Abell et al. (1989) Abell, G. O., Corwin, Harold G., J., & Olowin, R. P. 1989, ApJS, 70, 1, doi: 10.1086/191333
  * Adami et al. (1997) Adami, C., Biviano, A., & Mazure, A. 1997, Segregations in clusters of galaxies, arXiv, doi: 10.48550/ARXIV.ASTRO-PH/9709268
  * Ade et al. (2011) Ade, P. A. R., Aghanim, N., Arnaud, M., et al. 2011, Astronomy & Astrophysics, 536, A1, doi: 10.1051/0004-6361/201116464
  * Anbajagane et al. (2022) Anbajagane, D., Aung, H., Evrard, A. E., et al. 2022, MNRAS, 510, 2980, doi: 10.1093/mnras/stab3587
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, apj, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Barnes et al. (2017) Barnes, D. J., Kay, S. T., Henson, M. A., et al. 2017, MNRAS, 465, 213, doi: 10.1093/mnras/stw2722
  * Barsanti et al. (2016) Barsanti, S., Girardi, M., Biviano, A., et al. 2016, A&A, 595, A73, doi: 10.1051/0004-6361/201629012
  * Beers et al. (1990) Beers, T. C., Flynn, K., & Gebhardt, K. 1990, AJ, 100, 32, doi: 10.1086/115487
  * Behroozi et al. (2013) Behroozi, P. S., Loeb, A., & Wechsler, R. H. 2013, Journal of Cosmology and Astroparticle Physics, 2013, 019, doi: 10.1088/1475-7516/2013/06/019
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Bocquet & Carter (2016) Bocquet, S., & Carter, F. W. 2016, The Journal of Open Source Software, 1, 46, doi: 10.21105/joss.00046
  * Böhringer et al. (2004) Böhringer, H., Schuecker, P., Guzzo, L., et al. 2004, A&A, 425, 367, doi: 10.1051/0004-6361:20034484
  * Caldwell et al. (2016) Caldwell, C. E., McCarthy, I. G., Baldry, I. K., et al. 2016, Monthly Notices of the Royal Astronomical Society, 462, 4117, doi: 10.1093/mnras/stw1892
  * Caminha et al. (2016) Caminha, G. B., Grillo, C., Rosati, P., et al. 2016, Astronomy & Astrophysics, 587, A80, doi: 10.1051/0004-6361/201527670
  * Capasso et al. (2019) Capasso, R., Saro, A., Mohr, J. J., et al. 2019, MNRAS, 482, 1043, doi: 10.1093/mnras/sty2645
  * Comis et al. (2011) Comis, B., de Petris, M., Conte, A., Lamagna, L., & de Gregori, S. 2011, MNRAS, 418, 1089, doi: 10.1111/j.1365-2966.2011.19562.x
  * Contigiani et al. (2019) Contigiani, O., Hoekstra, H., & Bahé, Y. M. 2019, MNRAS, 485, 408, doi: 10.1093/mnras/stz404
  * Crawford et al. (2014) Crawford, S. M., Wirth, G. D., & Bershady, M. A. 2014, ApJ, 786, 30, doi: 10.1088/0004-637X/786/1/30
  * Cui et al. (2018) Cui, W., Knebe, A., Yepes, G., et al. 2018, MNRAS, 480, 2898, doi: 10.1093/mnras/sty2111
  * Dehnen (1993) Dehnen, W. 1993, MNRAS, 265, 250, doi: 10.1093/mnras/265.1.250
  * Diaferio (1999) Diaferio, A. 1999, MNRAS, 309, 610, doi: 10.1046/j.1365-8711.1999.02864.x
  * Diaferio & Geller (1997) Diaferio, A., & Geller, M. J. 1997, The Astrophysical Journal, 481, 633–643, doi: 10.1086/304075
  * Diaferio et al. (2005) Diaferio, A., Geller, M. J., & Rines, K. J. 2005, ApJ, 628, L97, doi: 10.1086/432880
  * Diemand et al. (2005) Diemand, J., Moore, B., & Stadel, J. 2005, Nature, 433, 389, doi: 10.1038/nature03270
  * Donahue et al. (2014) Donahue, M., Voit, G. M., Mahdavi, A., et al. 2014, ApJ, 794, 136, doi: 10.1088/0004-637X/794/2/136
  * Duffy et al. (2008a) Duffy, A. R., Schaye, J., Kay, S. T., & Vecchia, C. D. 2008a, Monthly Notices of the Royal Astronomical Society: Letters, 390, L64, doi: 10.1111/j.1745-3933.2008.00537.x
  * Duffy et al. (2008b) —. 2008b, Monthly Notices of the Royal Astronomical Society: Letters, 390, L64, doi: 10.1111/j.1745-3933.2008.00537.x
  * Einasto (1969) Einasto, J. 1969, Astronomische Nachrichten, 291, 97, doi: 10.1002/asna.19682910303
  * Evrard et al. (2008) Evrard, A. E., Bialek, J., Busha, M., et al. 2008, The Astrophysical Journal, 672, 122–137, doi: 10.1086/521616
  * Fadda et al. (1996a) Fadda, D., Girardi, M., Giuricin, G., Mardirossian, F., & Mezzetti, M. 1996a, ApJ, 473, 670, doi: 10.1086/178180
  * Fadda et al. (1996b) —. 1996b, ApJ, 473, 670, doi: 10.1086/178180
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Geller et al. (2013) Geller, M. J., Diaferio, A., Rines, K. J., & Serra, A. L. 2013, ApJ, 764, 58, doi: 10.1088/0004-637X/764/1/58
  * Gifford et al. (2017) Gifford, D., Kern, N., & Miller, C. J. 2017, The Astrophysical Journal, 834, 204, doi: 10.3847/1538-4357/834/2/204
  * Gifford et al. (2013) Gifford, D., Miller, C., & Kern, N. 2013, The Astrophysical Journal, 773, 116, doi: 10.1088/0004-637x/773/2/116
  * Gifford & Miller (2013) Gifford, D., & Miller, C. J. 2013, ApJ, 768, L32, doi: 10.1088/2041-8205/768/2/L32
  * Girardi et al. (1996) Girardi, M., Fadda, D., Giuricin, G., et al. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 94, Mapping, Measuring, and Modelling the Universe, ed. P. Coles, V. Martinez, & M.-J. Pons-Borderia, 221
  * Gómez et al. (2012) Gómez, P. L., Valkonen, L. E., Romer, A. K., et al. 2012, AJ, 144, 79, doi: 10.1088/0004-6256/144/3/79
  * Gruen et al. (2013) Gruen, D., Brimioulle, F., Seitz, S., et al. 2013, Monthly Notices of the Royal Astronomical Society, 432, 1455, doi: 10.1093/mnras/stt566
  * Halenka et al. (2022) Halenka, V., Miller, C. J., & Vansickle, P. 2022, ApJ, 926, 126, doi: 10.3847/1538-4357/ac4786
  * Herbonnet et al. (2020) Herbonnet, R., Sifón, C., Hoekstra, H., et al. 2020, Monthly Notices of the Royal Astronomical Society, 497, 4684, doi: 10.1093/mnras/staa2303
  * Hoekstra et al. (2015) Hoekstra, H., Herbonnet, R., Muzzin, A., et al. 2015, MNRAS, 449, 685, doi: 10.1093/mnras/stv275
  * Iannuzzi & Dolag (2012) Iannuzzi, F., & Dolag, K. 2012, MNRAS, 427, 1024, doi: 10.1111/j.1365-2966.2012.22017.x
  * Karman et al. (2015) Karman, W., Caputi, K. I., Grillo, C., et al. 2015, Astronomy & Astrophysics, 574, A11, doi: 10.1051/0004-6361/201424962
  * Keener (2010) Keener, R. 2010, Theoretical Statistics: Topics for a Core Course, Springer Texts in Statistics (Springer New York).  https://books.google.co.in/books?id=aVJmcega44cC
  * Klein et al. (2019) Klein, M., Israel, H., Nagarajan, A., et al. 2019, MNRAS, 488, 1704, doi: 10.1093/mnras/stz1491
  * Lemze et al. (2012) Lemze, D., Wagner, R., Rephaeli, Y., et al. 2012, ApJ, 752, 141, doi: 10.1088/0004-637X/752/2/141
  * Lotz et al. (2017) Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97, doi: 10.3847/1538-4357/837/1/97
  * Mamon et al. (2013) Mamon, G. A., Biviano, A., & Boué , G. 2013, Monthly Notices of the Royal Astronomical Society, 429, 3079, doi: 10.1093/mnras/sts565
  * McCarthy et al. (2016) McCarthy, I. G., Schaye, J., Bird, S., & Brun, A. M. C. L. 2016, Monthly Notices of the Royal Astronomical Society, 465, 2936, doi: 10.1093/mnras/stw2792
  * Melchior et al. (2015) Melchior, P., Suchyta, E., Huff, E., et al. 2015, MNRAS, 449, 2219, doi: 10.1093/mnras/stv398
  * Mercurio et al. (2021) Mercurio, A., Rosati, P., Biviano, A., et al. 2021, arXiv e-prints, arXiv:2109.03305.  https://arxiv.org/abs/2109.03305
  * Miller et al. (2016) Miller, C. J., Stark, A., Gifford, D., & Kern, N. 2016, The Astrophysical Journal, 822, 41, doi: 10.3847/0004-637x/822/1/41
  * Munari et al. (2013) Munari, E., Biviano, A., Borgani, S., Murante, G., & Fabjan, D. 2013, MNRAS, 430, 2638, doi: 10.1093/mnras/stt049
  * Nandra et al. (2012) Nandra, R., Lasenby, A. N., & Hobson, M. P. 2012, Monthly Notices of the Royal Astronomical Society, 422, 2931–2944, doi: 10.1111/j.1365-2966.2012.20618.x
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Old et al. (2014) Old, L., Skibba, R. A., Pearce, F. R., et al. 2014, MNRAS, 441, 1513, doi: 10.1093/mnras/stu545
  * Penna-Lima et al. (2017) Penna-Lima, M., Bartlett, J. G., Rozo, E., et al. 2017, A&A, 604, A89, doi: 10.1051/0004-6361/201629971
  * Pizzardo et al. (2022) Pizzardo, M., Diaferio, A., & Rines, K. 2022, arXiv e-prints, arXiv:2207.04982, doi: 10.48550/arXiv.2207.04982
  * Pizzardo et al. (2023) Pizzardo, M., Geller, M. J., Kenyon, S. J., Damjanov, I., & Diaferio, A. 2023, arXiv e-prints, arXiv:2307.07398, doi: 10.48550/arXiv.2307.07398
  * Postman et al. (2012) Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25, doi: 10.1088/0067-0049/199/2/25
  * Retana-Montenegro, E. et al. (2012) Retana-Montenegro, E., Van Hese, E., Gentile, G., Baes, M., & Frutos-Alfaro, F. 2012, A&A, 540, A70, doi: 10.1051/0004-6361/201118543
  * Rines et al. (2013) Rines, K., Geller, M. J., Diaferio, A., & Kurtz, M. J. 2013, ApJ, 767, 15, doi: 10.1088/0004-637X/767/1/15
  * Sandage (1970) Sandage, A. R. 1970, Physics Today, 23, 34
  * Sartoris, B. et al. (2020) Sartoris, B., Biviano, A., Rosati, P., et al. 2020, A&A, 637, A34, doi: 10.1051/0004-6361/202037521
  * Serra et al. (2011) Serra, A. L., Diaferio, A., Murante, G., & Borgani, S. 2011, MNRAS, 412, 800, doi: 10.1111/j.1365-2966.2010.17946.x
  * Sifón et al. (2013) Sifón, C., Menanteau, F., Hasselfield, M., et al. 2013, ApJ, 772, 25, doi: 10.1088/0004-637X/772/1/25
  * Sifón et al. (2016) Sifón, C., Battaglia, N., Hasselfield, M., et al. 2016, MNRAS, 461, 248, doi: 10.1093/mnras/stw1284
  * Stark et al. (2016) Stark, A., Miller, C. J., & Gifford, D. 2016, ApJ, 830, 109, doi: 10.3847/0004-637X/830/2/109
  * Stark et al. (2019) Stark, A., Miller, C. J., & Halenka, V. 2019, The Astrophysical Journal, 874, 33, doi: 10.3847/1538-4357/ab06fa
  * Stark et al. (2017) Stark, A., Miller, C. J., & Huterer, D. 2017, Phys. Rev. D, 96, 023543, doi: 10.1103/PhysRevD.96.023543
  * Umetsu et al. (2016) Umetsu, K., Zitrin, A., Gruen, D., et al. 2016, The Astrophysical Journal, 821, 116, doi: 10.3847/0004-637x/821/2/116
  * Umetsu et al. (2014) Umetsu, K., Medezinski, E., Nonino, M., et al. 2014, ApJ, 795, 163, doi: 10.1088/0004-637X/795/2/163
  * Vasiliev (2019) Vasiliev, E. 2019, MNRAS, 482, 1525, doi: 10.1093/mnras/sty2672
  * Williamson et al. (2011) Williamson, R., Benson, B. A., High, F. W., et al. 2011, ApJ, 738, 139, doi: 10.1088/0004-637X/738/2/139
  * Wing & Blanton (2013) Wing, J. D., & Blanton, E. L. 2013, ApJ, 767, 102, doi: 10.1088/0004-637X/767/2/102

##  Appendix A Alternate Density Profiles

In this Appendix, we present the fits for the additional density profiles,
where §4 only presented the fits to the NFW. As outlined in §4, we use the
$Z_{v}(r)$ from Figure 3 to account for radial suppression, and a degree-two
polynomial smoothing function for the phase-space edge to account for over and
under-fitting.

Here, we consider the Einasto (Einasto, 1969; Retana-Montenegro, E. et al.,
2012) and Dehnen (Dehnen, 1993) proifles. The NFW is most distinct from the
other two profiles in that it falls off less quickly in the outskirts (Navarro
et al., 1997). The potential for the Einasto can be obtained by solving the
Poisson equation, and is given by (Retana-Montenegro, E. et al., 2012):

$\displaystyle\Psi(r)$
$\displaystyle=\frac{-GM_{\text{tot}}}{r}\left[1-\frac{\Gamma(3n,(\frac{r}{h})^{1/n})}{\Gamma(3n)}+\frac{r}{h}\frac{\Gamma(2n,(\frac{r}{h}))^{1/n}}{\Gamma(3n)}\right]$
(A1)

where $\Gamma(\alpha,x)$ is the incomplete gamma function, given by
$\Gamma(\alpha,x)=\int_{x}^{\infty}t^{\alpha-1}e^{-t}dt$. The total mass,
$M_{\text{tot}}$, meanwhile is given by
$M_{\text{tot}}=4\pi\rho_{0}h^{3}n\Gamma(3n)$, and here $\rho_{0}$ is a
normalization factor, $h$ is the scale radius, and $n$ is the Einasto index.
Meanwhile, the potential for the Dehnen is (Dehnen, 1993):

$\displaystyle\Psi(r)$
$\displaystyle=\frac{GM_{\text{tot}}}{r_{s}}\frac{-1}{2-\gamma}\left[1-(\frac{r}{r+r_{s}})^{2-\gamma}\right],n\neq
2$ (A2) $\displaystyle\Psi(r)$
$\displaystyle=\frac{GM_{\text{tot}}}{r_{s}}\log(\frac{r}{r+r_{s}}),n=2$ (A3)

where the total mass $M_{\text{tot}}$ is again a normalization factor, $r_{s}$
is the scale radius, and $\gamma$ is the Dehnen index.

The top panel of Figure A1 shows the density profile posteriors for the
Einasto (left) and Dehnen (right) profiles. In each case, priors for the free
parameters are given in Table 3. We use the GTC package (Bocquet & Carter,
2016) to generate the histograms. The posteriors are highly degenerate, yet
yield relatively tightly constrained masses within $\sim 0.1$ dex. They are
similar to the top left panel of Figure 5, in that only a slim range of the
density profile’s parameters yields a physical escape profile, just like the
NFW. However in the case of the Einasto and Dehnen profiles, there is more
degeneracy than the case of the NFW. The bottom panel of Figure A1 is
analogous to the bottom panel of Figure 5, except using the additional density
profiles. Like the NFW, there is clear agreement between the dynamical fits
and the prediction from Gruen et al. (2013), and suggests that there is no
clear distinction between any of the profiles, highlighted as well in Table 2.

Table 3: The priors on the function parameters used in the potential profiles (Equations 12, A1, and A2 for the NFW, Einasto, and Dehnen profiles respectively). The parameters are drawn from a uniform distribution within the specified ranges. NFW: |
---|---
$r_{s}$ | $0.1<r_{s}/\hat{R}_{200}<10$
$\rho_{s}$ | $12<\log_{10}(\rho_{s})<16$
Einasto: |
$\rho_{0}$ |  $16<$log${}_{10}(\rho_{0})<22$
$h$ |  $-10<$log10(h)$<-1$
$n$ | $1<n<10$
Dehnen: |
$\gamma$ | $0.1<\gamma<5$
$r_{s}$ | $0.1<r_{s}/\hat{R}_{200}<10$
$M_{\text{tot}}$ | $14<\log(M_{\text{tot}})<18$

Figure A1:  Top panel: 2D posterior estimates (68% and 95%) of the Einasto
(left) and Dehnen (right) profile parameters. The covariance of the
distribution suggests that only a narrow range of the parameters values yield
physical mass constraints. However, the ranges are both wider than in the case
of the NFW, seen in the top panel of Figure 5. Bottom row: The escape profiles
estimated from drawing from the posterior distributions. The Gruen et al.
(2013) weak lensing estimate (cyan lines) is in-line with our
mass/concentration predictions for the two density profiles, as for the NFW.