11institutetext: Leibniz-Institute for Astrophysics, An der Sternwarte 16,
14482 Potsdam, Germany
11email<EMAIL_ADDRESS>22institutetext: Institute of Astronomy, National Central
University, Taoyuan 32001, Taiwan 33institutetext: Institut für Physik und
Astronomie, Universität Potsdam, Karl-Liebknecht-Straße 24/25, D-14476
Potsdam, Germany 44institutetext: INAF – Arcetri Astrophysical Observatory,
Largo Enrico Fermi 5, I-50125, Firenze, Italy 55institutetext: Department of
Astronomy, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH
44106, USA 66institutetext: Department of Physics, University of Oregon,
Eugene, OR 97403, USA 77institutetext: Department of Physics, University of
Surrey, Guildford GU2 7XH, UK 88institutetext: College of General Studies,
Yuan-Ze University, Taoyuan 32003, Taiwan

# Measuring galaxy cluster mass profiles into low acceleration regions with
galaxy kinematics

Pengfei Li , Humboldt fellow11 Yong Tian 22 Mariana P. Júlio 1133 Marcel S.
Pawlowski 11 Federico Lelli 44 Stacy S. McGaugh 55 James M. Schombert 66
Justin I. Read 77 Po-Chieh Yu 88 Chung-Ming Ko 22

(Received xxx; accepted xxx)

We probe the outer dynamical mass profile of 16 galaxy clusters from the
HIghest X-ray FLUx Galaxy Cluster Sample (HIFLUGCS) using galaxy kinematics.
Numerically solving the spherical Jeans equation, we parameterize the
dynamical mass profile and the galaxy velocity anisotropy profile using two
general functions to ensure that our results are not biased towards any
specific model. The mass-velocity anisotropy degeneracy is broken by using two
‘virial shape parameters’ that depend on the fourth moment of velocity
distribution. The resulting velocity anisotropy estimates consistently show a
nearly isotropic distribution in the inner regions, with an increasing radial
anisotropy towards large radii. We compare our derived dynamical masses with
those calculated from X-ray gas data assuming hydrostatic equilibrium, finding
that massive and rich relaxed clusters generally present consistent mass
measurements, while unrelaxed or low-richness clusters have systematically
larger total mass than hydrostatic mass by an average of 50%. This might help
alleviate current tensions in the measurement of $\sigma_{8}$, but it also
leads to cluster baryon fractions below the cosmic value. Finally, our
approach probes accelerations as low as $10^{-11}$ m s-2, comparable to the
outskirts of individual late-type galaxies. We confirm that galaxy clusters
deviate from the radial acceleration relation defined by galaxies.

###### Key Words.:

cosmology: observations — dark matter — galaxies: clusters: general —
galaxies: clusters: intracluster medium — X-rays: galaxies: clusters

##  1 Introduction

Clusters of galaxies present the most prominent structures in the Universe
given they are the largest self-gravitating objects. Their distribution and
abundance accommodate critical information on the geometry of the Universe and
the growth of the structures (Böhringer et al., 2014; Planck Collaboration et
al., 2016b), which often require the measurements of their total mass. Robust
estimates of cluster mass rely on the strong constraints from spatially
resolved dynamics, which also makes clusters excellent laboratories for
testing various dark matter (DM) models and alternative theories of gravity.

Hot X-ray emitting gas is the most commonly used dynamical tracer to derive
the dynamical mass profiles, given galaxy clusters are generally X-ray
selected and their X-ray surface brightness profiles are usually available.
Recent observations have managed to measure radially varying temperature
profiles (e.g. Ghirardini et al., 2019; Liu et al., 2022), resulting in more
accurate estimates for the total mass. The assumption of hydrostatic
equilibrium is sometimes in question due to the existence of non-thermal
pressure (Lau et al., 2009; Nelson et al., 2014). Perhaps more interestingly,
the cosmological constraints from cluster counts (Planck Collaboration et al.,
2016b) is in tension with that from the cosmic microwave background (CMB,
Planck Collaboration et al., 2016a), which can be resolved by introducing a
40% bias in mass. These discrepancies have motivated studies comparing
hydrostatic mass with that measured with other approaches, usually
gravitational lensing. Gravitational lensing (Postman et al., 2012; Umetsu et
al., 2016) measures the deflection of photons from background sources which
does not assume dynamical equilibrium. As such, it is considered the most
reliable method, but it has only be applied to a subset of galaxy clusters and
subject to projection effects.

Historically, using cluster galaxies as dynamical tracers was the very first
method to measure the dynamical mass of galaxy clusters and uncover their DM
content (Zwicky, 1933). However, this approach is relatively less developed
and applied, since it requires expensive spectroscopic measurements for a
large number of galaxies. On the technical side, it is also challenging to
derive 3D dynamical profiles from projected galaxy distributions and line-of-
sight velocities. There are three mass estimators that have been used in the
literature: virial theorem, caustic approach, solving Jeans equation.
Employing the virial theorem (Limber & Mathews, 1960; Biviano et al., 2006) is
the easiest, but it only provides the estimate of the virial mass. The caustic
approach (Diaferio & Geller, 1997; Diaferio, 1999; Serra et al., 2011)
considers the fact that gravitationally bounded galaxies cannot exceed the
escaping velocity, so one would expect a trumpet-shape distribution of
galaxies in the projected phase space, given the gravitational field decreases
towards large radii. This approach does not rely on dynamical equilibrium. The
Jeans equation (van der Marel, 1994) relates the phase-space distribution of
galaxies to the total gravitational potential. It assumes that galaxies are in
dynamical equilibrium. Deviations from the equilibrium are typically
accompanied by disturbed gas distributions and thereby can be partially
identified from X-ray images (Nagai et al., 2007). Both the caustic approach
and the Jeans modeling have to deal with the degeneracy between the density
profile and the velocity anisotropy. Strategies that have been implemented
include assuming a constant velocity anisotropy (e.g. Biviano et al., 2013),
or parameterizing its profile with a single parameter. The latter generally
assumes perfect isotropy either in the center (Mamon & Łokas, 2005; Tiret et
al., 2007) or at the scale radius of DM halos (Biviano et al., 2013). These
inevitably introduce some biases. Foëx et al. (2017) compared the cluster mass
estimated with these three approaches and found they are generally higher than
hydrostatic mass by 20-50%.

With more spectra of cluster galaxies becoming available and especially the
advent of SDSS-V (Almeida et al., 2023), it is now worth further developing
the approach with galaxy kinematics. In this paper, we numerically solve the
Jeans equation and adopt more general functions to describe the dynamical mass
profile and velocity anisotropy. Our results are thus not biased towards any
specific halo model. We break the degeneracy between the density profile and
velocity anisotropy with fourth order of velocity dispersion, i.e. two virial
shape parameters. We test the approach using 16 HIFLUGCS clusters and compare
with their hydrostatic mass profiles. Section 2 describes our cluster sample;
Section 3 introduces the approach we employed; Section 4 presents the
comparison with hydrostatic mass for six disturbed clusters; Section 5
presents the detailed dynamical analysis for ten X-ray relaxed clusters. We
discuss our results and conclude the paper in Section 6.

##  2 HIFLUGCS clusters

###  2.1 X-ray data

In this paper, we adopted a subsample of HIFLUGCS clusters, which were
observed by both the ROSAT All-Sky Survey and XMM-Newton, leading to high-
quality X-ray data (Zhang et al., 2011). The complete sample of 63 clusters is
X-ray selected and constructed by Reiprich & Böhringer (2002). Chen et al.
(2007) provide fits of their surface brightness profiles with the $\beta$
model (Cavaliere & Fusco-Femiano, 1976),

$S(R)=S_{0}\Big{(}1+R^{2}/r_{c}^{2}\Big{)}^{-3\beta+\frac{1}{2}},$ (1)

where $R$ is the projected radius, $r_{c}$ is the core radius, $S_{0}$ is the
central surface brightness and $\beta$ is the power index. Assuming spherical
symmetry, the above function can be de-projected to derive the electron number
density profile assuming that X-rays are produced by thermal bremsstrahlung,

$n_{e}=n_{c}\Big{(}1+r^{2}/r_{c}^{2}\Big{)}^{-3\beta/2}$ (2)

where $r$ is the 3D radius, and $n_{c}$ is the central electron number
density. Chen et al. (2007) adopted $H_{0}=50$ km s-1 Mpc-1 for distance
estimates, which is significantly lower than current values. In order to
maintain consistency with current estimates, we adopt $H_{0}=70$ km s-1 Mpc-1
and re-scale his results for core radius and electron number density:
$n_{c}=n_{\rm c,original}/1.4$ and $r_{c}=1.4r_{\rm c,original}$. The enclosed
gas mass profile is then

$M_{\rm
gas}(<r)=4\pi\frac{A}{Z}m_{p}n_{c}r_{c}^{3}\int_{0}^{x}(1+x^{2})^{-\frac{3}{2}\beta}x^{2}{\rm
d}x,$ (3)

where $x=r/r_{c}$; $m_{p}$ is the mass of protons; the mean nuclear mass and
charge numbers are $A\simeq 1.4$ and $Z\simeq 1.2$ for the intracluster medium
(ICM) with 0.3 solar abundance (Anders & Grevesse, 1989). Assuming the ICM is
in hydrostatic equilibrium, the dynamical mass profile can be derived through

$M_{\rm hydro}(<r)=\frac{3\beta kT_{h}r}{G\mu
m_{p}}\frac{(r/r_{c})^{2}}{1+(r/r_{c})^{2}},$ (4)

where $T_{h}$ is the temperature of the hotter bulk component in the two-phase
model. The cooler phase is included to account for the possible cooling core.
Chen et al. (2007) pointed out $T_{h}$ provides a better measure for the
gravitational potential and total mass than the single emission-weighted
temperature $T_{m}$. We, therefore, use $T_{h}$ when deriving the hydrostatic
mass in this paper. The mean molecular weight $\mu$ is given by

$\mu=\rho_{\rm gas}/(m_{p}n_{\rm gas})\simeq(2X+0.75Y+0.56Z)^{-1},$ (5)

where $n_{\rm gas}$ is the total number density including electrons, protons,
ionized Helium, and other ionized elements; $X$, $Y$ and $Z$ are the mass
fractions of Hydrogen, Helium, and other elements, respectively. For 0.3 solar
abundance, $X\simeq 0.716$ and $Y\simeq 0.278$. This gives $\mu\simeq 0.6$.

###  2.2 Optical data

The collection of high-quality spectro-photometric data for the HIFLUGCS
clusters is an achievement of decades’ efforts, including the compilations in
Andernach et al. (2005); Zhang et al. (2011). In this paper, we utilize the
compiled sample from Tian et al. (2021), which includes the comprehensive
memberships assembled in the literature and organized from SIMBAD (Wenger et
al., 2000). This sample is cleaned by excluding uncertain or repeated members.
Each cluster in the sample has a single brightest cluster galaxy (BCG), and
its position is defined as the optical center. Tian et al. (2021) required the
offset between the optical center and X-ray weighted center to be smaller than
60 kpc. Therefore, we do not distinguish between optical or X-ray centers in
this paper.

In order to effectively constrain dynamical mass profiles, we require at least
three radial bins for each cluster with at least 25 galaxies in each bin. We,
therefore, exclude clusters with less than 75 member galaxies. We also remove
two galaxy groups, NGC 4636 and A1060, from our cluster sample. Eventually, we
retain 16 clusters in the analysis. According to Zhang et al. (2011), six of
them are disturbed, so their galaxy kinematics may not trace the gravitational
potential. Therefore, we will study disturbed and undisturbed clusters
separately.

##  3 Methodology

###  3.1 Galaxy kinematics

Using galaxy kinematics to determine the dynamics of clusters assumes that
galaxies obey the collisionless Boltzmann equation (Binney & Tremaine, 2008),

$\frac{{\rm d}f}{{\rm d}t}=\frac{\partial f}{\partial
t}+\nabla_{x}f\cdot\@vec{v}-\nabla_{v}f\cdot\nabla_{x}\Phi,$ (6)

where $f(\@vec{x},\@vec{v})$ is the distribution function of galaxies, and
$\Phi$ is the total gravitational potential. The dynamical equilibrium is a
critical assumption. It is also a central question we aim to answer by
comparing the results with those using hydrostatic equilibrium. Since we
assume spherical symmetry for all the considered clusters, the Boltzmann
equation can be simplified as the spherical Jeans equation (Binney & Tremaine,
2008),

$\frac{1}{\nu}\frac{\partial}{\partial
r}(\nu\sigma^{2}_{r})+\frac{2\beta(r)\sigma^{2}_{r}}{r}=-\frac{GM(<r)}{r^{2}},$
(7)

where $\nu=\int f{\rm d^{3}}v$ is the number density of tracer galaxies,
$\sigma_{r}$ is the radial velocity dispersion,
$\beta=1-\sigma^{2}_{t}/\sigma^{2}_{r}$ is the velocity anisotropy with
$\sigma_{t}$ the tangential anisotropy, and $M(<r)$ is the enclosed total
mass. Once $\nu$, $\sigma_{r}$ and $\beta$ are determined, equation 7 provides
a direct measurement for the enclosed total mass. In fact, it is more common
to use the integrated equation (van der Marel, 1994),

$\sigma^{2}_{r}=\frac{1}{\nu(r)g(r)}\int_{r}^{\infty}\frac{GM(<r^{\prime})\nu(r^{\prime})}{r^{\prime
2}}g(r^{\prime}){\rm d}r^{\prime},$ (8)

where $g(r)=\exp{(\int\frac{2\beta(r)}{r}{\rm d}r)}$. In observations, only
projected distributions of galaxies ($\Sigma_{\rm gal}$) and line-of-sight
velocity dispersion ($\sigma_{\rm los}$) are observable. Therefore, it is
important to have an equation relating $\sigma_{\rm los}$ to $\sigma_{\rm r}$,
which has been derived by Binney & Mamon (1982),

$\sigma^{2}_{\rm los}(R)=\frac{2}{\Sigma_{\rm
gal}(R)}\int_{R}^{\infty}\Big{(}1-\beta\frac{R^{2}}{r^{2}}\Big{)}\frac{\nu(r)\sigma^{2}_{r}r}{\sqrt{r^{2}-R^{2}}}{\rm
d}r.$ (9)

###  3.2 GravSphere code

In this paper, we use the GravSphere code by Read & Steger (2017), which was
written for star clusters and dwarf galaxies using individual stars as
dynamical tracers. Since it solves the same equations (eq. 7 and 9), we can
transplant it to galaxy clusters by using individual galaxies as dynamical
tracers and changing the corresponding parameters and priors. GravSphere uses
parameterized functions to describe the unknown variables in equation 7 and 9.
To guarantee that the chosen functions can provide good fits for diverse
objects, they generally have many free parameters. Galaxy number density is
modeled with three Plummer spheres (Plummer, 1911),

$\nu(r)=\sum^{3}_{i=1}\frac{3N_{i}}{4\pi
a^{3}_{i}}\Big{(}1+\frac{r^{2}}{a^{2}_{i}}\Big{)}^{-5/2}$ (10)

where the values of $N_{i}$ are chosen to make sure $\int\nu(r){\rm
d^{3}}\@vec{r}=1$. We set broad boundaries for the parameters:
$10^{-4}<N_{i}<100$ and 50 kpc $<a_{i}<$ 2000 kpc. It is possible to include
more Plummer spheres and hence more free parameters. However, we found that
three Plummer spheres with six parameters are enough for all the considered
systems. The values of $(N_{i},a_{i})$ are determined by fitting the projected
number density to observed values, which is analytically given by

$\Sigma_{\rm
gal}=\sum^{3}_{i=1}\frac{N_{i}a^{2}_{i}}{\pi(a^{2}_{i}+R^{2})^{2}}.$ (11)

An example fit with cluster A0085 is given in the left panel of Figure 1.

The total enclosed mass $M(<r)$ is modeled with the cNFWt profile (Read et
al., 2018), which is based on the Navarro-Frenk-White model (NFW, Navarro et
al., 1996),

$\rho_{\rm
NFW}=\frac{\rho_{s}}{\frac{r}{r_{s}}\Big{(}1+\frac{r}{r_{s}}\Big{)}^{2}}.$
(12)

The cNFWt model modifies the NFW profile at both small and large radii. At
small radii, the density profile is described as the cNFW model (Read et al.,
2016) and its enclosed mass is given by

$M_{\rm cNFW}(<r)=M_{\rm NFW}(<r)f^{n},$ (13)

where $f=\tanh{\Big{(}\frac{r}{r_{c}}\Big{)}}$, and the index $n$ controls how
cuspy the halo is. A NFW-like cusp has $n=0$, while a complete core requires a
minimum value of $n=1$. The corresponding density profile reads as

$\rho_{\rm cNFW}(<r)=\rho_{\rm NFW}f^{n}+\frac{nf^{n-1}(1-f^{2})}{4\pi
r^{2}r_{c}}M_{\rm NFW}(<r).$ (14)

At $r>r_{t}$, the enclosed mass is given by

$M_{\rm cNFWt}(<r)=M_{\rm cNFW}(<r_{t})+4\pi\rho_{\rm
cNFW}(r_{t})\frac{r_{t}^{3}}{3-\delta}\Big{[}\Big{(}\frac{r}{r_{t}}\Big{)}^{3-\delta}-1\Big{]},$
(15)

where $\delta$ is the outer slope. In total, the cNFWt profile has six
parameters, which can accommodate different inner and outer slopes. We set
loose boundaries as: $13<\log(M_{200}/M_{\odot})<20$, $0.1<c_{200}<100$, 10
kpc $<r_{c}<$ 3000 kpc, $-0.5<n<2.0$, 10 kpc $<r_{t}<$ 5000 kpc and
$0<\delta<3$. The total mass $M_{200}$ of galaxy clusters generally ranges
from 1014 to 3$\times$1015 $M_{\odot}$, and the halo concentration $c_{200}$
is typically between 1 and 10 (Umetsu et al., 2016). Therefore, the allowed
ranges are much wider than the expected values. We allow $r_{c}$ and $r_{t}$
to vary from galaxy scale to beyond cluster scale. The core-cusp control
parameter $n$ is typically within (0, 1). The additional outer slope parameter
$\delta$ is introduced to accommodate density profiles less steep than the NFW
model, which can be switched off by setting a large value of $r_{t}$ if
necessary.

Figure 1: Example fits of projected galaxy number density profiles and line-
of-sight velocity dispersion profiles. The cluster shown here is A0085. The
galaxy surface number density profile is fitted with three Plummer spheres,
which is then used as input in the projected Jeans equation. The line-of-sight
velocity dispersion profile is used to determine the total mass profile,
parameterized using the cNFWt profile with six parameters. Dark and light
shadow regions show the 1$\sigma$ and 2$\sigma$ confidence intervals,
respectively.

By default, the output of the GravSphere code is the enclosed dark matter halo
mass, which is calculated by subtracting baryonic mass from the total mass.
Since we will have to model stellar mass and gas mass distributions
separately, we simply set a negligible baryonic mass in the code. As a result,
the output mass profile, as described by the cNFWt model, is the total
enclosed mass. The cNFWt model has sufficient freedom to accommodate various
mass profiles (see an example fit in the right panel of Figure 1). This is a
critical reason why we do not use the NFW model, which is generally believed
to work well in galaxy clusters. Therefore, we treat the cNFWt model simply as
a fitting function rather than a halo model, and set loose boundaries for its
parameters without imposing cosmological priors.

The velocity anisotropy parameter $\beta(r)$ is also modelled by a fitting
function to avoid irregular behavior,

$\beta(r)=\beta_{0}+(\beta_{\infty}-\beta_{0})\frac{1}{1+\big{(}\frac{r_{0}}{r}\big{)}^{n}},$
(16)

where $\beta_{0}$ and $\beta_{\infty}$ are the values at $r=0$ and $r=\infty$,
respectively; $r_{0}$ and $n$ are used to characterize the radial shape. The
boundary for $n$ is set as: $1<n<3$; while $r_{0}$ is allowed to vary within
$0.5R_{\rm half}<r_{0}<2R_{\rm half}$, where $R_{\rm half}$ is the radius
enclosing 50% of the total galaxies. In order to avoid the infinite value for
a full tangential velocity dispersion, Read et al. (2006) defines a
symmetrized anisotropy parameter,

$\tilde{\beta}=\frac{\beta}{2-\beta}.$ (17)

According to this definition, the full tangential and radial velocity
dispersion correspond to $\tilde{\beta}=-1$ and $\tilde{\beta}=1$,
respectively. Since the velocity distribution is expected to be nearly
isotropic in the center of galaxy clusters due to the strong gravitational
field, we set the boundary as $-0.2<\tilde{\beta}_{0}<0.2$ in the center,
while we allow a larger anisotropy at large radii by setting
$-0.2<\tilde{\beta}_{\infty}<1.0$. We will find the set ranges are
sufficiently large for all the considered clusters.

Notoriously, the velocity anisotropy parameter is degenerated with the density
profile, as a larger anisotropy would lead to a larger radial velocity
dispersion, which is positively correlated with the enclosed mass profile. The
degeneracy can be broken by proper motions, which provide two additional
projected equations, similar to eq. 9. However, proper motions are hardly
measurable for cluster galaxies. Merrifield & Kent (1990) found that the
fourth order of the velocity dispersion can be related to the enclosed mass
profiles in two separate and independent ways, i.e.

$\displaystyle v_{s1}$ $\displaystyle=$
$\displaystyle\frac{2}{5}\int_{0}^{\infty}GM\nu(5-2\beta)\sigma^{2}_{r}r{\rm
d}r$ (18) $\displaystyle=$ $\displaystyle\int_{0}^{\infty}\Sigma_{\rm
gal}\langle v^{4}_{\rm los}\rangle R{\rm d}R,$ $\displaystyle v_{s2}$
$\displaystyle=$
$\displaystyle\frac{4}{35}\int_{0}^{\infty}GM\nu(7-6\beta)\sigma^{2}_{r}r^{3}{\rm
d}r$ (19) $\displaystyle=$ $\displaystyle\int_{0}^{\infty}\Sigma_{\rm
gal}\langle v^{4}_{\rm los}\rangle R^{3}{\rm d}R,$

where $\langle v^{4}_{\rm los}\rangle=\int v^{4}_{\rm los}f{\rm
d^{3}}\@vec{v}$. Therefore, the two Virial Shape parameters can help break the
degeneracy.

In total, GravSphere solves the Jeans equation with ten parameters (six in
cNFWt and four in $\beta(r)$). We impose flat priors on these fitting
parameters within the aforementioned hard boundaries. The parameter space is
explored using the Markov Chain Monte Carlo method with the $emcee$ hammer by
Foreman-Mackey et al. (2013). We use 250 walkers and run 20000 steps to ensure
the convergence of the chain. We check that changing the number of steps does
not affect our results.

Figure 2: Dynamical mass profiles of six disturbed clusters of galaxies. Solid
lines are the hydrostatic mass profiles derived from the surface brightness
fits using the $\beta$ function from Chen et al. (2007). Points with errorbar
show the total mass profiles from the line-of-sight velocity dispersion. The
3D radii are chosen according to the binned projected radii to avoid
oversampling. That the dynamical mass profiles determined from different
tracers are inconsistent presumably indicates non-equilibrium conditions
stemming from recent or ongoing mergers

.

###  3.3 Binning data with Binulator

The robust estimate of the dynamical mass is contingent on the accurate
calculation of the velocity dispersion, which typically requires a large
number of tracers in each bin. This is difficult to achieve in many cases. As
such, Collins et al. (2021) introduced a separate routine for data binning,
Binulator, prior to running GravSphere. Instead of directly calculating
velocity dispersion from measured line-of-sight velocities, Binulator fits a
generalized Gaussian function within each bin and estimates the mean, variance
and kurtosis. The generalized function is given by

$p_{i}=\frac{\beta}{2\alpha\Gamma(1/\beta)}\exp{\big{[}-\big{(}|v_{\rm
los,i}-\mu|/\alpha\big{)}^{\beta}\big{]}},$ (20)

where $\alpha$, $\beta$, $\mu$ are fitting parameters, the index $i$ labels
each individual galaxy, and $\Gamma(x)$ is the Gamma function. As before, we
set loose boundaries for these parameters: 50 km/s $<\alpha<$ 2000 km/s,
$1<\beta<10$, $|\mu|<$ 1000 km/s. Based on the fit, the velocity dispersion is
given by $\sigma^{2}_{\rm los}=\alpha^{2}\Gamma(3/\beta)/\Gamma(1/\beta)$.
Since there is no error measurement for line-of-sight velocity in the
available catalogs, we assume a 10% uncertainty on the velocity for all
cluster galaxies. This is a conservative assumption, since galaxy velocities
are measured via spectroscopy with high accuracy. This uncertainty generates
an error probability function, which is convolved with the generalized
Gaussian function. Eventually, the parameters are determined for each bin by
maximizing the likelihood function (with the $emcee$ hammer),

$\displaystyle\mathcal{L}=\prod_{i=1}^{N}p_{i},$ (21)

where $N$ is the total number of galaxies in each bin. Collins et al. (2021)
showed that the Binulator routine can robustly estimate velocity dispersion
with 25 galaxies/bin. We hence set the minimum bin size as 25 galaxies, while
for rich clusters we can also choose to use a larger bin size.

##  4 Disturbed clusters

Galaxy kinematics traces the gravitational potential only if the cluster is
relaxed. This is not always the case, as clusters could be undergoing merging
and thereby the distribution and motions of galaxies are disturbed. Unrelaxed
clusters can have some signatures in their X-ray images. Nagai et al. (2007);
Ventimiglia et al. (2008) examined some simulated clusters and classified
relaxed and unrelaxed clusters based on the morphology of the mock X-ray
images. Vikhlinin et al. (2009) applied their procedure to observed clusters,
and identified unrelaxed clusters as those with secondary maxima, filamentary
X-ray structures, or significant isophotal centroid shifts in their X-ray
images. These clusters may deviate from hydrostatic equilibrium, as they
present a systematic offset from the $M_{\rm tot}-T_{X}$ relation (Mathiesen &
Evrard, 2001; Kravtsov et al., 2006). Kravtsov et al. (2006) found that the
total mass is higher than the expected mass by $17\pm 5\%$ at the same
temperature. Therefore, Vikhlinin et al. (2009) suggested scaling up the total
mass from X-ray temperature by 17%. However, the actual offset of observed
clusters in the $M_{\rm tot}-T_{X}$ is unclear, which has to be calibrated
through weak lensing analysis.

In our sample, six clusters are disturbed according to their X-ray images
(Vikhlinin et al., 2009; Zhang et al., 2011). We derive their hydrostatic mass
from the surface brightness fits using the $\beta$ function (eq. 4), and plot
them in Figure 2. We also plot the mass profiles derived from the line-of-
sight velocity dispersion, which are drawn from the best-fit cNFWt function,
but we only show the discrete points with errors rather than the full curve.
This is intended to avoid over-constraining where there is no actual data. The
points are chosen at the projected radii of the binned data, although the 3D
and projected radii do not exactly correspond to each other. The errors are
estimated from the Markov Chain, so they are formal uncertainties that should
be taken with a grain of salt.

Figure 2 shows that the total mass profiles derived from galaxy kinematics for
five clusters are above the hydrostatic mass curves at all radii. The
differences are larger than 17%, so even if we scale up the hydrostatic mass
as for simulated clusters, the two measurements would not agree. This suggests
that the merging process affects cluster galaxies and X-ray gas with different
significance. A3526 is the only cluster that presents consistent measurements.
Noticeably, its size is significantly smaller than other clusters. This may
imply that X-ray gas and satellite galaxies are equally distorted for small
clusters, though the sample is too small to make a solid conclusion.

Both measurements of the total mass should not be thought reliable given their
X-ray images present the signatures of un-relaxation. We include these
clusters in the paper to illustrate how unrelaxed clusters may look in their
dynamical mass profiles. We will exclude them from our dynamical tests.

Figure 3: Mass profiles of six undisturbed clusters of galaxies. Points with
errorbar are the total mass profiles from velocity dispersion. Blue solid and
dashed lines are the hydrostatic mass and gas mass profiles from the X-COP
project, respectively; while black solid and dashed lines show the
corresponding results from the surface brightness fits using the $\beta$
functions in Chen et al. (2007). Green lines show the stellar mass
distributions. For these undisturbed clusters, the dynamical mass profiles
from galaxy kinematics are generally consistent with hydrostatic mass
profiles, but extend to larger radii in general. Possible reasons for the
inconsistency in some clusters are discussed in the text.

Figure 4: Continuation of Figure 3 for the remaining four clusters.

##  5 Undisturbed clusters

###  5.1 Baryonic mass profiles

The ten clusters identified as undisturbed by Zhang et al. (2011) comprise our
main sample of interest. We plot their gas mass, stellar mass, hydrostatic
mass, and dynamical mass from galaxy kinematics in Figure 3 and 4. The gas
mass profiles are derived by integrating their deprojected surface brightness
fits from Chen et al. (2007). Five clusters have measurements in the XMM
Cluster Outskirts Project (X-COP; Eckert et al., 2019; Ettori et al., 2019;
Ghirardini et al., 2019), so we overplot their results for comparison. The gas
mass distributions from X-COP are supposed to be more extended, given the aim
of the X-COP is to explore the outer regions of clusters and study the growth
of structures. Three clusters in our plot present more extended gas profiles
from the $\beta$ function fits. These are extrapolations rather than true gas
distributions. We extrapolate all the density profiles from the $\beta$
function fits to the outermost radii with the binned galaxy data. We notice
that the gas mass from the $\beta$ function is systematically higher than that
in X-COP for all five clusters. This is expected given the well-known fact
that gas density decreases more quickly than the $\beta$ function predicts at
large radii (Vikhlinin et al., 2006). Therefore, the $\beta$ function should
be used with caution towards larger radii.

Though gas mass dominates the total baryonic mass in galaxy clusters, the
stellar mass contained in galaxies provide significant contributions at small
radii. We distinguish BCG and satellite galaxies. Since galaxy kinematics
cannot probe the inner regions, we do not need to model the stellar mass
distributions of BCGs, but can simply treat them as a point mass. We derive
the stellar mass of each BCG based on their $K_{s}$-band absolute magnitude
using the relation from Cappellari (2013),

$\log M_{\rm BCG}/M_{\odot}\simeq 10.58-0.44\times(M_{K_{s}}+23),$ (22)

where the absolute magnitude is calculated from the observed apparent
magnitude with the Two Micron All Sky Survey (2MASS, Skrutskie et al., 2006)
using the distance measured by Hubble flow. The Cappellar relation corresponds
to a stellar mass-to-light ratio between 1.3 and 1.6 $M_{\odot}/L_{\odot}$ for
our sample, which is higher than the predicted value (1.2
$M_{\odot}/L_{\odot}$) of stellar population models (Schombert et al., 2022).

Table 1: Mass budget of each relaxed cluster. The values of $r_{500}$ and $M_{\rm gas,500}$ are taken from Zhang et al. (2011); $M_{\star,500}$ is calculated using the $M_{\rm gas,500}-M_{\star,500}$ relation from Chiu et al. (2018); $M_{\star,\rm BCG}$ is derived from the $K_{s}$ band magnitude using the approach by Cappellari et al. (2013); $M_{\rm dyn,500}$ is interpolated from the best-fit cNFWt mass profile; $M_{\rm hydro,500}$ is recalculated at the same $r_{500}$ to ensure consistency; baryonic fraction is derived based on the dynamical mass $f_{b}=M_{b}/M_{\rm dyn,500}$. Cluster | $r_{500}$ | $M_{\rm gas,500}$ | $M_{\star,500}$ | $M_{\star,\rm BCG}$ | $M_{\rm dyn,500}$ | $M_{\rm hydro,500}$ |  $M_{\rm dyn,500}$/$M_{\rm hydro,500}$ | $f_{b}$
---|---|---|---|---|---|---|---|---
| (kpc) | ($10^{13}M_{\odot}$) | ($10^{12}M_{\odot}$) | ($10^{11}M_{\odot}$) | ($10^{14}M_{\odot}$) | ($10^{14}M_{\odot}$) |  |
A0085 | 1217 | 6.67 | 4.39 | 19.05 | 8.44 | 4.66 | 1.81 | 0.086
A0262 | 755 | 1.08 | 1.47 | 5.62 | 1.37 | 0.83 | 1.64 | 0.094
A0496 | 967 | 2.79 | 2.60 | 10.72 | 4.61 | 2.39 | 1.92 | 0.069
A0576 | 869 | 2.00 | 2.13 | 14.13 | 3.94 | 2.18 | 1.80 | 0.060
A1795 | 1085 | 4.95 | 3.68 | 3.89 | 4.59 | 4.41 | 1.04 | 0.117
A2029 | 1247 | 8.24 | 4.99 | 41.69 | 8.07 | 6.37 | 1.27 | 0.113
A2142 | 1371 | 13.40 | 6.68 | 12.02 | 8.05 | 7.47 | 1.08 | 0.176
A2589 | 837 | 1.77 | 1.98 | 10.72 | 3.06 | 1.81 | 1.69 | 0.068
A3158 | 1013 | 3.75 | 3.11 | 14.79 | 4.59 | 3.55 | 1.29 | 0.091
A3571 | 1133 | 5.16 | 3.77 | 18.62 | 7.36 | 5.02 | 1.47 | 0.078

As mentioned in Section 3.2, galaxy number distributions are described by
three Plummer spheres $\nu(r)$. For simplicity, we assume the satellite
galaxies share the same stellar mass. This should not affect our results as
long as the sample of galaxies is sufficiently large. With this assumption,
stellar mass distributions are simply galaxy number distributions up to a
proportional constant, which can be determined by the total stellar mass. Chiu
et al. (2018) reported the correlation between the gas mass and stellar mass
within $r_{500}$,

$M_{\star}(<r_{500})=4\times 10^{12}M_{\odot}\Big{(}\frac{M_{\rm
gas}(<r_{500})}{5.7\times 10^{13}M_{\odot}}\Big{)}^{0.6}.$ (23)

We take the values of $M_{\rm gas}(<r_{500})$ from Zhang et al. (2011) and
derive the enclosed stellar mass. The stellar mass distribution is then given
by

$M_{\star}(<r)=(M_{\star}(<r_{500})-M_{\rm BCG})\nu(r)/\nu(r_{500})+M_{\rm
BCG},$ (24)

which is plotted in green solid lines. Figure 3 shows the baryonic mass is
generally dominated by stellar mass in the inner regions of clusters.

###  5.2 Dynamical mass profiles

Similar to disturbed clusters, we derive the hydrostatic mass from the surface
brightness fits using the $\beta$ function (Chen et al., 2007) with a constant
temperature. For the five clusters in the X-COP, we overplot their results
from Ettori et al. (2019), which are supposed to be more robust because the
radial variations of the temperature are taken into account. Eckert et al.
(2022) also included the pressure measurements from the Planck Sunyaev-
Zel’dovich (SZ) survey (Planck Collaboration et al., 2014), so their dynamical
mass profiles are more extended than the gas distributions.

Interestingly, though the gas mass profiles from the X-COP and the $\beta$
function fits differ significantly at large radii, the total mass profiles are
quite consistent. This is partially because hydrostatic mass depends on the
density slope rather than the absolute mass density. The slope of the $\beta$
function is shallower than the true density distribution in the outer regions,
while the assumed constant temperature over-predicts the actual values. These
two deviations have opposite effects on the measurements of hydorstatic mass.
As a result, the hydrostatic mass from the $\beta$ function fits with a
constant temperature is fairly reliable in spite of the apparent technical
caveats.

In the inner regions, some clusters have a cool core within 0.1$R_{500}$. The
positive temperature gradient leads to a smaller pressure gradient, and
thereby the derived hydrostatic mass is smaller than that from a constant
temperature. On the other hand, the $\beta$ function is known to underestimate
the inner slope of the electron number density. The two opposite effects are
competing. As a result, the three cool-core clusters, A0085, A1795 and A2029
present roughly consistent measurements with constant and varying
temperatures. A2142 is identified as a non-cool-core cluster by Zhang et al.
(2011), but Eckert et al. (2022) show the core temperature is lower than the
peak value by 25%. The only non-cool-core cluster, A3158, shows that the
hydrostatic mass from the $\beta$ function fit is apparently lower than that
in X-COP at $r<200$ kpc, which might be due to the shallower inner slope of
the $\beta$ function. This does not affect the remaining five clusters for
which there are no X-COP results, since we mostly focus on the regions beyond
200 kpc.

In Figure 3 and 4, we also plot the dynamical mass profiles measured from the
line-of-sight velocity dispersion as discrete points with errors. The
corresponding total density profiles are presented in the appendix (Figure 7
and 8). Among those five clusters in X-COP, four clusters present roughly
consistent mass measurements. These are rich and massive clusters, so their
spatial dynamics can be well resolved by galaxy kinematics. The galaxy
distribution of A2029 is much more extended than its gas distribution. As
such, we can only compare our dynamical mass with the extrapolated hydrostatic
mass from the $\beta$ function fit at large radii. The good agreement suggests
that the $\beta$ function extrapolation might be a satisfactory approximation.
The galaxy kinematics of A0085 present consistent dynamical mass within
$r_{500}$, while beyond this radius, they overshoot the hydrostatic mass. This
indicates that although the X-ray gas is relaxed, the cluster galaxies in the
outskirts are not. It could be that relaxation is quicker at small radii but
takes longer at large radii. So after merging, the outskirt galaxies need more
time to reach equilibrium. One might think that A0085 is growing its structure
by accreting more galaxies. However, if these outskirt galaxies have higher
velocity dispersion than allowed by the total gravitational potential, they
cannot be gravitational bounded within the cluster. Therefore, it is more
likely that the hydrostatic mass is underestimated at large radii.

The remaining five clusters are low-mass and low-richness. Figures 3 and 4
show that the dynamical masses from galaxy kinematics for four clusters
overshoot their hydrostatic masses, while only A3571 presents an acceptable
agreement. The inconsistency may imply that these cluster galaxies are away
from dynamical equilibrium, though their X-ray images do not present abnormal
structures. If true, this would suggest that galaxies and intracluster medium
have different time scales for relaxation. However, our results are
qualitatively consistent with what Foëx et al. (2017) found. This indicates
that their hydrostatic masses are likely being underestimated.

In Table 1, we present the mass budget for each relaxed cluster. We also
compared the total masses measured from galaxy kinematics and hydrostatic
equilibrium. Our results show that the total mass from galaxy kinematics is
higher than hydrostatic mass by 4% to 92%. The mean percentage is $\sim$ 50%.
This could affect the $\sigma_{8}$ tension (Blanchard & Ilić, 2021). The
$\sigma_{8}$ parameter describes the magnitude of matter fluctuations in the
later Universe, which can be calculated from the CMB fluctuations via
extrapolations (Planck Collaboration et al., 2016a). It can also be directly
measured through cluster counts. The cluster sample for cosmological interests
is SZ selected, but the cluster mass is calibrated using hydrostatic mass
(Planck Collaboration et al., 2016b). So if the hydrostatic mass is
underestimated, the measured $\sigma_{8}$ would be lower as well. Blanchard &
Ilić (2021) pointed out the tension can be resolved if the hydrostatic mass is
scaled up by 40%, which is roughly consistent with our mass estimates.
However, our sample is quite small, so that it is unclear if the true mass
deviation measured from galaxy kinematics and hydrostatic equilibrium is
around 40% or not. To completely resolve the $\sigma_{8}$ tension, one would
also need to understand the cosmic shear analysis from weak lensing
measurements (Asgari et al., 2021), which also presents this tension with the
CMB results. Our results could be a step forward if statistically confirmed
with larger cluster samples.

Though promising for the $\sigma_{8}$ tension, our measurements result in the
baryonic fraction that is in tension with the cosmic value, $f_{b}=0.16$.
Table 1 shows that most of our clusters have a baryonic fraction much smaller
than 0.16. Therefore, the tension between the early-Universe and late-Universe
measurements is hardly resolved in the $\Lambda$CDM cosmology.

Figure 5: Radial distributions of the velocity anisotropy parameter $\beta$
for the ten undisturbed clusters. Light shadow regions show the 1$\sigma$
confidence intervals. The mean values of $\beta$ are consistent with an
isotropic velocity distribution within the errors at small radii, while at
large radii the velocity anisotropies are clearly presented in spite of the
large errors. The two big error bars at small and large radii represent the
allowed ranges set by the priors on $\beta_{0}$ and $\beta_{\infty}$,
respectively. The lower boundary is $-0.5$ at all radii.

###  5.3 Velocity anisotropy

Figure 5 plots the radial distributions of the velocity anisotropy parameter
for the ten undisturbed clusters. The derivation of their values relies on the
two virial shape parameters, $v_{s1}$ and $v_{s2}$, which help break the
$\rho-\beta$ degeneracy. However, $v_{s1}$ and $v_{s2}$ are not radially
dependent functions, but two single values. As such, they do not put strong
constraints. The measured velocity anisotropies hence present large
uncertainties. At small radii, the uncertainties are mostly around 0.1, while
they could be as large as 0.2 at large radii. The richest cluster, A2029, has
the smallest uncertainty, smaller than 0.1, given that the line-of-sight
velocity dispersion profile is better measured thanks to the larger sample of
galaxy tracers.

In spite of the large uncertainties, nine of the ten clusters present similar
trends for their mean velocity anisotropies: nearly isotropic in the inner
regions while apparently anisotropic in the outskirts. The boundary we impose
on the inner anisotropy is
$|\tilde{\beta_{0}}|=|\frac{\beta_{0}}{2-\beta_{0}}|<0.2$, which corresponds
to $-0.5<\beta_{0}<0.33$. Figure 5 shows all the inner anisotropies are well
within the set range. This justifies our priors and suggests that the isotropy
in the inner regions is truly physical rather than artificially imposed.
However, the anisotropies observed at large radii may not be solid given the
larger uncertainties, except for A2029. Our results are qualitatively
consistent with the expectations. Since the gravitational field is strong at
small radii, one expect galaxy motions are chaotic and hence present nearly
isotropic velocity distributions; while at large radii, the velocity
distribution could be more radially dominated. Whether our results are
quantitatively consistent with cosmological simulations is worth further
exploring.

###  5.4 Radial acceleration relation

The extended dynamical mass profiles probed by galaxy kinematics puts strong
constraints on dynamical relations. Here we test the radial acceleration
relation (RAR, McGaugh et al., 2016; Lelli et al., 2017b), a tight correlation
relating the observed acceleration from rotation curves $g_{\rm obs}$ and that
from the baryonic distributions $g_{\rm bar}$,

$g_{\rm obs}=\frac{g_{\rm bar}}{1-e^{-\sqrt{g_{\rm bar}/g_{\dagger}}}}.$ (25)

This relation was established statistically with 153 late-type galaxies (Lelli
et al., 2016a) and applied to early-type galaxies (Lelli et al., 2017a;
Shelest & Lelli, 2020). It also holds in individual galaxies once the
uncertainties on stellar mass-to-light ratio, galaxy distances and disc
inclinations are taken into account (Li et al., 2018). The empirical relation
is in line with modified Newtonian dynamics (MOND, Milgrom, 1983).

Figure 6: Radial acceleration relation of ten undisturbed galaxy clusters.
Black points show the results with both measured dynamical and baryonic
masses; while grey points show the radii where there are only measured
dynamical masses. The gas masses at these radii are derived from
extrapolations based on the surface brightness fits. The dotted line is the
line of unity, and the red line is the RAR established with late-type galaxies
(McGaugh et al., 2016; Lelli et al., 2017b). Light grey points represent the
SPARC galaxies from McGaugh et al. (2016). The total accelerations overshoot
the RAR in galaxies, suggesting a missing baryon problem in the inner regions
of galaxy clusters for MOND Milgrom (1983).

The RAR has also been explored on BCG-cluster scales. For example, Tian et al.
(2020) employed 20 massive clusters from the Cluster Lensing And Supernova
survey with Hubble (CLASH, Postman et al., 2012) with total mass profiles
measured from gravitational lensing (Umetsu et al., 2016). They found a
correlation between $g_{\rm obs}$ and $g_{\rm bar}$ that it is systematically
offset from the RAR observed in galaxies, with a significantly larger
acceleration scale. Eckert et al. (2022) investigated the RAR using 12 X-COP
clusters with dynamical mass profiles measured from hydrostatic equilibrium
and the SZ effect. They also reported higher dynamical accelerations than
expected from the RAR.

We calculate the total accelerations from the dynamical mass profiles measured
with galaxy kinematics,

$g_{\rm dyn}=\frac{GM_{\rm dyn}(<r)}{r^{2}}.$ (26)

Since galaxy distributions can be quite extended, their kinematics helps probe
lower accelerations than other approaches. The CLASH sample used by Tian et
al. (2020) probes accelerations as low as $10^{-10}$ m s-2. With the help of
the SZ technique, Eckert et al. (2022) extended the acceleration down to
slightly above $10^{-11}$ m s-2, similar to that observed in individual
clusters (Fig. 6).

It is challenging to derive the corresponding baryonic acceleration, which is
the summation of the gas and stellar contributions,

$g_{\rm bar}=\frac{G\big{(}M_{\rm gas}(<r)+M_{\star}(<r)\big{)}}{r^{2}}.$ (27)

In the outskirts, X-ray surface brightness becomes too low to be observable.
Therefore, we can only estimate the gas mass via extrapolations using the
$\beta$ function fits. For the five clusters in the X-COP, we adopt the
modified $\beta$ function from Vikhlinin et al. (2006),

$n^{2}_{e}=n^{2}_{0}\frac{(x/r_{c})^{\alpha}}{(1+x^{2}/r_{c}^{2})^{3\beta-\alpha/2}}\frac{1}{(1+x^{\gamma}/r^{\gamma}_{s})^{\epsilon/\gamma}},$
(28)

where $x=r/r_{500}$, and $\gamma=3$ is generally fixed. The free parameter
space is comprised of ($n_{0}$, $\alpha$, $\beta$, $\epsilon$, $r_{c}$,
$r_{s}$). Ghirardini et al. (2019) claimed that the density profile is
universal for the 12 X-COP clusters, so they can be described by a single set
of values: $\ln n_{0}=-4.44$, $\alpha=0.89$, $\beta=0.43$, $\epsilon=2.86$,
$\ln r_{c}=-2.99$, $\ln r_{s}=-0.29$. However, their Figure 3 shows the
electron number density profiles are quite diverse in the inner and outer
regions. To achieve more accurate extrapolations towards large radii, we only
adopt their values of $\alpha$, $\beta$, and $r_{c}$, but re-fit the outer
slope $\epsilon$, scale length $r_{s}$ and the overall normalization factor
$n_{0}$. The resultant fits properly describe the diversity of the density
profiles in the outskirts with greater accuracy, so we extrapolate the gas
mass distributions with these fits where X-ray data are absent.

Figure 6 plots the total acceleration measured with the line-of-sight velocity
dispersions of cluster member galaxies against the baryonic acceleration for
the ten relaxed clusters. Regions where the gas mass is estimated by
extrapolating the surface brightness fits should be given less credit. Similar
to previous works, we find that the total acceleration is systematically
higher than expected from the RAR defined by individual galaxies. This is
particularly true at intermediate radii. Black points roughly form a line,
parallel to the galactic RAR. A similar behavior has been reported by Tian et
al. (2020) with the CLASH sample. Tian et al. (2020) found that the higher
total acceleration imply a larger critical acceleration ($g_{\dagger}\sim
10^{-9}$ m s-2) in clusters than that in galaxies ($g_{\dagger}\sim 10^{-10}$
m s-2, McGaugh et al., 2016). The high critical acceleration is a reflection
of the missing baryon problem, a well-known problem for MOND in galaxy
clusters.

At larger radii and lower acceleration, the data appear to converge with the
galactic RAR. This indicates that the missing baryon problem that MOND suffers
in clusters is concentrated towards their centers (Sanders, 2003; Angus et
al., 2008), while at large radii the baryonic mass seems consistent with the
MOND prediction. However, we caution that the total baryonic mass in the
outskirts are extrapolations from surface brightness fits. As a result, the
total baryonic mass may be overestimated. Whether the missing baryon problem
appears in the outskirts of galaxy clusters remain an open question that will
be explored with more extended X-ray data in the future.

##  6 Discussion and conclusion

In this paper, we have applied the spherical Jeans equation to the motions of
galaxies in clusters to infer their orbital anisotropies and dynamical mass
profiles. We tested this approach with 16 HIFLUGCS clusters and derived the
dynamical mass profiles from the observed galaxy distributions and line-of-
sight velocity dispersion. We distinguish six clusters with indications of
disturbances in their X-ray images, and found their dynamical masses measured
from galaxy kinematics are systematically larger than the hydrostatic masses.
For the X-ray relaxed clusters, half of them present consistent dynamical mass
measurements from galaxy kinematics and hydrostatic equilibrium. These are
mostly massive and rich clusters. The rest of the clusters present higher
dynamical mass, which suggests that the hydrostatic mass might be
underestimated. This could be a step towards solving the $\sigma_{8}$ tension.
The latter requires a higher mass calibration than hydrostatic mass. But it
also causes difficulty in reproducing the cosmic baryonic fraction.
Alternatively, it could imply that these clusters are not in dynamical
equilibrium. Therefore, our approach can also be used to study the dynamical
state of galaxy clusters.

The galaxy samples for the clusters we employed in this paper are likely
incomplete, given we included only galaxies with a membership probability of
100%. This can affect our results in two ways: (1) it underestimates the
galaxy number density and surface density, which appear in equation 7 and 9;
(2) it affects the derivation of the line-of-sight velocity dispersion.
However, we do not expect the consequence of the incompleteness to be
significant. In equation 7 and 9, the galaxy number density appears in both
the numerator and denominator, so a constant fraction of incompleteness
cancels out. This largely reduces the effect, and only radially variable
incompleteness can affect our results. The line-of-sight velocity dispersion
is estimated with the Binulator which fits a group of galaxies to the
generalized Gaussian function. Missing a few galaxies does not seriously
affect the fits. To examine the effect of the incompleteness in detail, we
shall turn to more complete samples such as that from SDSS-V (Almeida et al.,
2023).

A big advantage of using galaxies to constrain cluster dynamics is that they
can trace the gravitational potential far beyond X-ray regions, where the
acceleration can be extremely low. In this paper, we showed the lowest
acceleration probed by the 10 HIFLUGCS clusters is slightly below $10^{-11}$ m
s-2. The low acceleration region is critical for testing some alternative
theories of gravity such as MOND, as well as different dark matter models like
fluid dark matter (Khoury, 2015). However, it also raises the concern that
galaxies in the outskirts of a cluster may be still on the way of collapsing
into the center, so they may have not reached the dynamical equilibrium yet.
Since there is no X-ray data in the outskirts, we cannot robustly test the
equilibrium, but only provided a rough investigation by extrapolating the
hydrostatic mass profiles with the fits using the $\beta$ function or modified
$\beta$ function when the actual data are available. The comparison supports
the application of galaxy kinematics in the outskirts, which makes our
approach promising. A robust calibration can only be carried out by cross-
checking the results from weak lensing measurements. Weak lensing measures the
projected mass along the line of sight, which includes any foreground or
background mass, whereas our dynamical analysis constrains the spherically
averaged 3D mass distribution. As such, a comparison between the two requires
both a projection of our 3D mass distribution on the sky and analysis of mock
data from full cosmological simulations to test the impact of line-of-sight
structures. For these reasons, we will consider detailed comparisons with weak
lensing studies in future work.

Our dynamical analysis refines and reinforces previous results that MOND
suffers a missing baryon problem in rich, $X$-ray emitting galaxy clusters.
This discrepancy is pronounced at intermediate radii where the dynamical
acceleration exceeds that predicted by the galaxy-calibrated RAR applied to
the observed baryons. Numerous solutions have been proposed for this, such as
unseen baryons (Milgrom, 2015), massive neutrinos (Sanders, 2003), sterile
neutrinos (Angus & Diaferio, 2011), mixed dark matter/alternative gravity
models (Berezhiani & Khoury, 2015), making the Lagrangian a function of
potential depth as well as acceleration (Hodson & Zhao, 2017), tinkering with
the MOND interpolation function (Zhao & Famaey, 2012), or some combination of
these. The cluster data converge towards the galaxy RAR at low accelerations
in their outskirts. Therefore, it seems the discrepancy is restricted to
cluster cores, while the total baryonic mass may be consistent with that
predicted by the Baryonic Tully-Fisher relation (McGaugh, 2005; Lelli et al.,
2016b; McGaugh et al., 2021). However, our baryonic mass estimates in the
outskirts are extrapolations and only three clusters extend to the low
acceleration region. More clusters with robust gas mass measurements in the
outskirts are necessary to further examine this behavior.

###### Acknowledgements.

P.L. is supported by the Alexander von Humboldt Foundation. Y.T. is supported
by Taiwan National Science and Technology Council NSTC 110-2112-M-008-015-MY3.
M.S.P. acknowledges funding of a Leibniz-Junior Research Group (project number
J94/2020) and a KT Boost Fund by the German Scholars Organization and Klaus
Tschira Stiftung. S.S.M. and J.M.S. are supported in part by NASA ADAP grant
80NSSC19k0570. S.S.M. also acknowledges support from NSF PHY-1911909. C.M.K.
is supported by Taiwan NSTC 111-2112-M-008-013.

## References

  * Almeida et al. (2023) Almeida, A., Anderson, S. F., Argudo-Fernández, M., et al. 2023, arXiv e-prints, arXiv:2301.07688
  * Andernach et al. (2005) Andernach, H., Tago, E., Einasto, M., Einasto, J., & Jaaniste, J. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 329, Nearby Large-Scale Structures and the Zone of Avoidance, ed. A. P. Fairall & P. A. Woudt, 283–287
  * Anders & Grevesse (1989) Anders, E. & Grevesse, N. 1989, Geochim. Cosmochim. Acta., 53, 197
  * Angus & Diaferio (2011) Angus, G. W. & Diaferio, A. 2011, MNRAS, 417, 941
  * Angus et al. (2008) Angus, G. W., Famaey, B., & Buote, D. A. 2008, MNRAS, 387, 1470
  * Asgari et al. (2021) Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104
  * Berezhiani & Khoury (2015) Berezhiani, L. & Khoury, J. 2015, Phys. Rev. D, 92, 103510
  * Binney & Mamon (1982) Binney, J. & Mamon, G. A. 1982, MNRAS, 200, 361
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Biviano et al. (2006) Biviano, A., Murante, G., Borgani, S., et al. 2006, A&A, 456, 23
  * Biviano et al. (2013) Biviano, A., Rosati, P., Balestra, I., et al. 2013, A&A, 558, A1
  * Blanchard & Ilić (2021) Blanchard, A. & Ilić, S. 2021, A&A, 656, A75
  * Böhringer et al. (2014) Böhringer, H., Chon, G., & Collins, C. A. 2014, A&A, 570, A31
  * Cappellari (2013) Cappellari, M. 2013, ApJ, 778, L2
  * Cappellari et al. (2013) Cappellari, M., Scott, N., Alatalo, K., et al. 2013, MNRAS, 432, 1709
  * Cavaliere & Fusco-Femiano (1976) Cavaliere, A. & Fusco-Femiano, R. 1976, A&A, 49, 137
  * Chen et al. (2007) Chen, Y., Reiprich, T. H., Böhringer, H., Ikebe, Y., & Zhang, Y. Y. 2007, A&A, 466, 805
  * Chiu et al. (2018) Chiu, I., Mohr, J. J., McDonald, M., et al. 2018, MNRAS, 478, 3072
  * Collins et al. (2021) Collins, M. L. M., Read, J. I., Ibata, R. A., et al. 2021, MNRAS, 505, 5686
  * Diaferio (1999) Diaferio, A. 1999, MNRAS, 309, 610
  * Diaferio & Geller (1997) Diaferio, A. & Geller, M. J. 1997, ApJ, 481, 633
  * Eckert et al. (2022) Eckert, D., Ettori, S., Pointecouteau, E., van der Burg, R. F. J., & Loubser, S. I. 2022, A&A, 662, A123
  * Eckert et al. (2019) Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40
  * Ettori et al. (2019) Ettori, S., Ghirardini, V., Eckert, D., et al. 2019, A&A, 621, A39
  * Foëx et al. (2017) Foëx, G., Böhringer, H., & Chon, G. 2017, A&A, 606, A122
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Ghirardini et al. (2019) Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41
  * Hodson & Zhao (2017) Hodson, A. O. & Zhao, H. 2017, A&A, 598, A127
  * Khoury (2015) Khoury, J. 2015, Phys. Rev. D, 91, 024022
  * Kravtsov et al. (2006) Kravtsov, A. V., Vikhlinin, A., & Nagai, D. 2006, ApJ, 650, 128
  * Lau et al. (2009) Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129
  * Lelli et al. (2016a) Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016a, AJ, 152, 157
  * Lelli et al. (2016b) Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016b, ApJ, 816, L14
  * Lelli et al. (2017a) Lelli, F., McGaugh, S. S., & Schombert, J. M. 2017a, MNRAS, 468, L68
  * Lelli et al. (2017b) Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017b, ApJ, 836, 152
  * Li et al. (2018) Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2018, A&A, 615, A3
  * Limber & Mathews (1960) Limber, D. N. & Mathews, W. G. 1960, ApJ, 132, 286
  * Liu et al. (2022) Liu, A., Bulbul, E., Ramos-Ceja, M. E., et al. 2022, arXiv e-prints, arXiv:2210.00633
  * Mamon & Łokas (2005) Mamon, G. A. & Łokas, E. L. 2005, MNRAS, 363, 705
  * Mathiesen & Evrard (2001) Mathiesen, B. F. & Evrard, A. E. 2001, ApJ, 546, 100
  * McGaugh (2005) McGaugh, S. S. 2005, ApJ, 632, 859
  * McGaugh et al. (2016) McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Physical Review Letters, 117, 201101
  * McGaugh et al. (2021) McGaugh, S. S., Lelli, F., Schombert, J. M., et al. 2021, AJ, 162, 202
  * Merrifield & Kent (1990) Merrifield, M. R. & Kent, S. M. 1990, AJ, 99, 1548
  * Milgrom (1983) Milgrom, M. 1983, ApJ, 270, 365
  * Milgrom (2015) Milgrom, M. 2015, MNRAS, 454, 3810
  * Nagai et al. (2007) Nagai, D., Vikhlinin, A., & Kravtsov, A. V. 2007, ApJ, 655, 98
  * Navarro et al. (1996) Navarro, J. F., Eke, V. R., & Frenk, C. S. 1996, MNRAS, 283, L72
  * Nelson et al. (2014) Nelson, K., Lau, E. T., & Nagai, D. 2014, ApJ, 792, 25
  * Planck Collaboration et al. (2014) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A29
  * Planck Collaboration et al. (2016a) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016a, A&A, 594, A13
  * Planck Collaboration et al. (2016b) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016b, A&A, 594, A24
  * Plummer (1911) Plummer, H. C. 1911, MNRAS, 71, 460
  * Postman et al. (2012) Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25
  * Read et al. (2016) Read, J. I., Agertz, O., & Collins, M. L. M. 2016, MNRAS, 459, 2573
  * Read & Steger (2017) Read, J. I. & Steger, P. 2017, MNRAS, 471, 4541
  * Read et al. (2018) Read, J. I., Walker, M. G., & Steger, P. 2018, MNRAS, 481, 860
  * Read et al. (2006) Read, J. I., Wilkinson, M. I., Evans, N. W., Gilmore, G., & Kleyna, J. T. 2006, MNRAS, 367, 387
  * Reiprich & Böhringer (2002) Reiprich, T. H. & Böhringer, H. 2002, ApJ, 567, 716
  * Sanders (2003) Sanders, R. H. 2003, MNRAS, 342, 901
  * Schombert et al. (2022) Schombert, J., McGaugh, S., & Lelli, F. 2022, AJ, 163, 154
  * Serra et al. (2011) Serra, A. L., Diaferio, A., Murante, G., & Borgani, S. 2011, MNRAS, 412, 800
  * Shelest & Lelli (2020) Shelest, A. & Lelli, F. 2020, A&A, 641, A31
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Tian et al. (2020) Tian, Y., Umetsu, K., Ko, C.-M., Donahue, M., & Chiu, I. N. 2020, ApJ, 896, 70
  * Tian et al. (2021) Tian, Y., Yu, P.-C., Li, P., McGaugh, S. S., & Ko, C.-M. 2021, ApJ, 910, 56
  * Tiret et al. (2007) Tiret, O., Combes, F., Angus, G. W., Famaey, B., & Zhao, H. S. 2007, A&A, 476, L1
  * Umetsu et al. (2016) Umetsu, K., Zitrin, A., Gruen, D., et al. 2016, ApJ, 821, 116
  * van der Marel (1994) van der Marel, R. P. 1994, MNRAS, 270, 271
  * Ventimiglia et al. (2008) Ventimiglia, D. A., Voit, G. M., Donahue, M., & Ameglio, S. 2008, ApJ, 685, 118
  * Vikhlinin et al. (2009) Vikhlinin, A., Burenin, R. A., Ebeling, H., et al. 2009, ApJ, 692, 1033
  * Vikhlinin et al. (2006) Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9
  * Zhang et al. (2011) Zhang, Y. Y., Andernach, H., Caretta, C. A., et al. 2011, A&A, 526, A105
  * Zhao & Famaey (2012) Zhao, H. & Famaey, B. 2012, Phys. Rev. D, 86, 067301
  * Zwicky (1933) Zwicky, F. 1933, Helvetica Physica Acta, 6, 110

##  Appendix A Total density profiles of ten relaxed clusters.

Figure 7: Density profiles corresponding to the total cumulative mass profiles
in Figure 3. Dark and light shadow area are 1$\sigma$ and 2$\sigma$ regions.
Note that these are not dark matter density profiles but the total density
profiles including baryonic contributions.

Figure 8: Continuation of Figure 7 for the remaining four clusters.