11institutetext:  Department of Astronomy, University of Geneva, Ch. d’Ecogia
16, CH-1290 Versoix, Switzerland
11email<EMAIL_ADDRESS>22institutetext: INAF, Osservatorio di
Astrofisica e Scienza dello Spazio, via Piero Gobetti 93/3, 40129 Bologna,
Italy 33institutetext: INFN, Sezione di Bologna, viale Berti Pichat 6/2,
I-40127 Bologna, Italy 44institutetext: IRAP, Université de Toulouse, CNRS,
CNES, UPS, Toulouse, France 55institutetext: European Southern Observatory,
Karl-Schwarzschild-Str. 2, 85748, Garching, Germany 66institutetext: Centre
for Space Research, North-West University, Potchefstroom 2520, South Africa

# The gravitational field of X-COP galaxy clusters

D. Eckert 11 S. Ettori 2233 E. Pointecouteau 44 R. F. J. van der Burg 55 S. I.
Loubser 66

The mass profiles of massive dark matter halos are highly sensitive to the
nature of dark matter and potential modifications of the theory of gravity on
large scales. The $\Lambda$ cold dark matter (CDM) paradigm makes strong
predictions on the shape of dark matter halos and on the dependence of the
shape parameters on halo mass, such that any deviation from the predicted
universal shape would have important implications for the fundamental
properties of dark matter. Here we use a set of 12 galaxy clusters with
available deep X-ray and Sunyaev-Zel’dovich data to constrain the shape of the
gravitational field with an unprecedented level of precision over two decades
in radius. We introduce a nonparametric framework to reconstruct the shape of
the gravitational field under the assumption of hydrostatic equilibrium and
compare the resulting mass profiles to the expectations of Navarro-Frenk-White
(NFW) and Einasto parametric mass profiles. On average, we find that the NFW
profile provides an excellent description of the recovered mass profiles, with
deviations of less than 10% over a wide radial range. However, there appears
to be more diversity in the shape of individual profiles than can be captured
by the NFW model. The average NFW concentration and its scatter agree very
well with the prediction of the $\Lambda$CDM framework. For a subset of
systems, we disentangle the gravitational field into the contribution of
baryonic components (gas, brightest cluster galaxy, and satellite galaxies)
and that of dark matter. The stellar content dominates the gravitational field
inside $\sim 0.02R_{500}$ but is responsible for only 1-2% of the total
gravitational field inside $R_{200}$. The total baryon fraction reaches the
cosmic value at $R_{200}$ and slightly exceeds it beyond this point, possibly
indicating a mild level of nonthermal pressure support ($10-20\%$) in cluster
outskirts. Finally, the relation between observed and baryonic acceleration
exhibits a complex shape that strongly departs from the radial acceleration
relation in spiral galaxies, which shows that the aforementioned relation does
not hold at the galaxy-cluster scale.

###### Key Words.:

X-rays: galaxies: clusters - Galaxies: clusters: general - Galaxies: groups:
general - Galaxies: clusters: intracluster medium - cosmology: large-scale
structure

##  1 Introduction

a Numerical simulations of cosmic structure formation have long predicted that
the shape of collapsed halos is universal over a wide range of halo masses. In
the cold dark matter (CDM) paradigm, dark matter (DM) halos are expected to
follow the so-called Navarro-Frenk-White (NFW) profile (Navarro et al., 1996,
1997), in which the matter density forms a central “cusp” with $\rho(r)\propto
r^{-1}$ and gradually steepens toward the outskirts as $\rho(r)\propto
r^{-3}$, with a characteristic radius $r_{s}$ that defines the scale at which
the logarithmic slope of the density profile has an isothermal value of $-$2
(i.e., $d\ln\rho/d\ln r|r_{s}=-2$). The characteristic shape is thought to
arise from the universality of the structure formation process (Bullock et
al., 2001; Macciò et al., 2008), and the distribution and redshift evolution
of the halo structural parameters encode important cosmological information
(e.g., Duffy et al., 2008; Dutton & Macciò, 2014). As discussed in the
original work by Navarro, Frenk, and White, the total mass at a given
overdensity (typically 200 times the critical density of the universe enclosed
in a sphere with radius $R_{200}$) correlates with the mass concentration of
the halo, $c=R_{200}/r_{s}$. This correlation reflects the conditions in which
a halo assembled, with halos formed earlier being more concentrated because
the background density was larger back in cosmological time (Wechsler et al.,
2002). More recently, several works (see, e.g., Dutton & Macciò, 2014; Ludlow
et al., 2016; Ragagnin et al., 2021) have modeled and quantified the effects
of halo accretion histories and of the parameters that describe the
cosmological background model on the halo mass concentration. Recent N-body
simulations have also shown that the radial profiles of massive halos deviate
from the NFW prediction (Navarro et al., 2004; Diemer & Kravtsov, 2014; Klypin
et al., 2016; Ludlow et al., 2016) and can be better fitted by functional
forms that include a power-law dependence of the logarithmic radial slope,
such as $d\ln\rho/d\ln r\propto r^{-\alpha}$. The Einasto (1965) profile is an
example of this revised functional form, where the Einasto index, $\alpha$,
can be related to the power-law index of the primordial power spectrum,
$n_{s}$ (Ludlow & Angulo, 2017; Brown et al., 2020).

Additionally, deviations from the universal NFW shape can provide important
clues to the nature of DM. For instance, in cases where the DM self-
interaction (SIDM) cross section is non-negligible, high-density regions are
expected to be homogenized, leading to flatter cores compared to the cuspy
profiles expected in the CDM scenario (Robertson et al., 2021). “Warm” dark
matter (WDM) particle candidates with masses of a few keV would free-stream
across the high-density regions of halos, thereby flattening the observed
profiles (Macciò et al., 2013; Ludlow et al., 2016). Modified gravity theories
such as MOdified Newtonian Dynamics (MOND; Milgrom, 1983) and its relativistic
extensions (e.g., Bekenstein, 2004; Milgrom, 2009; Skordis & Złośnik, 2021)
can also be tested by comparing the observed gravitational field to
predictions (see Famaey & McGaugh, 2012, for a review).

As the most massive collapsed structures in the Universe, galaxy clusters are
privileged sites for studying the properties of the gravitational field over a
wide range of densities. Given their deep potential well, gravitational
processes largely dominate over feedback processes, such that galaxy clusters
are essentially closed boxes for baryons (White et al., 1993; Evrard, 1997;
Eke et al., 1998; Kravtsov et al., 2005). Recent studies demonstrate that the
baryon fraction enclosed within the virial radius of massive clusters is very
close to the universal baryon fraction (Gonzalez et al., 2007; Laganá et al.,
2013; Mantz et al., 2014; Chiu et al., 2018; Eckert et al., 2019), with the
majority of the baryonic mass ($\sim 90\%$) in the hot intracluster medium
(ICM) and the remaining 10% in stars. The total baryonic content (gas and
stars) of these structures is directly observable. Galaxy clusters thus
represent ideal structures for studying the properties of the gravitational
field in DM halos and its interplay with the various baryonic components.

The mass profiles of galaxy clusters have been the subject of numerous recent
studies (see Pratt et al., 2019, for a review). The shape of the mass
distribution can be estimated through weak and strong gravitational lensing
(Okabe et al., 2013; Umetsu et al., 2016; Umetsu & Diemer, 2017; Tam et al.,
2020), the dynamics of cluster galaxies (Biviano et al., 2013; Munari et al.,
2014; Geller et al., 2014), and the hydrostatic equilibrium (HSE) equation
(Pointecouteau et al., 2005; Vikhlinin et al., 2006; Ettori et al., 2010;
Bartalucci et al., 2019; Ettori et al., 2019). The NFW profile usually
provides an adequate description of the data at hand, irrespective of the
adopted method (Umetsu et al., 2016; Ettori et al., 2019). Models that include
a large central core, such as the Burkert (Salucci & Burkert, 2000) or
isothermal sphere models (King, 1962), provide a significantly worse
representation of the data (Umetsu & Diemer, 2017; Ettori et al., 2019).
However, the data quality is usually insufficient to distinguish between the
various models that include a central cusp, in particular between the NFW and
Einasto profiles. Finally, galaxy cluster mass profiles were also extensively
used to test modified gravity theories (e.g., Sanders, 1999, 2003;
Pointecouteau & Silk, 2005; Wilcox et al., 2015; Ettori et al., 2017; Mitchell
et al., 2018).

In this paper we investigate in detail the shape of the gravitational field in
local galaxy clusters under the HSE assumption. We make use of the data
acquired for the _XMM Cluster Outskirts Project_ (X-COP; Eckert et al., 2017),
a very large program on _XMM-Newton_ that provides complete X-ray mapping out
to the virial radius for a sample of 12 nearby systems selected from the
_Planck_ Sunyaev-Zel’dovich (SZ) survey. We present a novel nonparametric
technique for recovering hydrostatic mass profiles from joint X-ray and SZ
data, which we include in a global framework for the estimation of hydrostatic
mass profiles. In Sect. 4 we validate our reconstruction technique using mock
_XMM-Newton_ observations of a synthetic NFW cluster. We apply our technique
to the 12 X-COP systems with the aim of comparing NFW and Einasto profiles.
For a subset of systems, we use measurements of the gas mass as well as the
stellar mass within the brightest cluster galaxy (BCG) and satellite galaxies
to break down the gravitational field into its baryonic and DM components.
Finally, we study the relation between the baryonic and total gravitational
acceleration and compare it with the radial acceleration relation (RAR) of
rotationally supported galaxies (McGaugh et al., 2016). In a companion paper
(hereafter Paper II), we discuss our measurements of the Einasto shape
parameter and perform a detailed comparison to cosmological simulations that
include self-interacting DM.

Throughout the paper we assume a Planck 2015 $\Lambda$CDM cosmology (Planck
Collaboration XIII, 2016) with $H_{0}=67.8$ km/s/Mpc and $\Omega_{m}=0.308$.
The codes developed in the framework of this project are made publicly
available. In particular, we release our mass reconstruction code in the form
of an easy-to-use Python package named hydromass. Our _XMM-Newton_ mock data
generator can be downloaded from GitHub.

##  2 Sample and data reduction

###  2.1 The X-COP sample

X-COP (Eckert et al., 2017) is a very large program on _XMM-Newton,_ totaling
1.5 Ms of observing time. The original data set was awarded during _XMM-
Newton_ AO-13 (proposal ID 074441, PI: Eckert) and completed in 2015. The
primary goal of the project is to use the combination of X-ray data from _XMM-
Newton_ with SZ data from _Planck_ to probe the state of the ICM out to the
virial radius. The sample was selected from the first _Planck_ SZ catalog
(PSZ1; Planck Collaboration XXIX, 2014) according to the following criteria:
(i) _Planck_ PSZ1 signal-to-noise greater than 12, (ii) apparent size
$\theta_{500}>10^{\prime}$ to ensure that all systems are spatially resolved
by _Planck_ , (iii) redshift in the range $0.04<z<0.1$, and Galactic
absorption column density $N_{H}<10^{21}$ cm-2.

For more details about the sample we refer the reader to Eckert et al. (2017).
The full observation log was presented in Appendix F of Ghirardini et al.
(2019). For the present work, we add a set of 6 additional _XMM-Newton_
pointings that were obtained during AO-17 (observation ID 082365, PI: Eckert)
on two clusters (A3266 and A2029) to homogenize the data set.

###  2.2 XMM-Newton data analysis

With the exception of the two systems with new observations (A3266 and A2029),
we employ the publicly available data products presented in Ghirardini et al.
(2019)111https://dominiqueeckert.wixsite.com/xcop and already used to study
the gravitational field of galaxy clusters in several previous works (Ettori
et al., 2017, 2019; Pradyumna et al., 2021; Haridasu et al., 2021; Harikumar &
Biesiada, 2022). Here we briefly recall the data reduction chain developed for
X-COP. More details are available in Sect. 2 of Ghirardini et al. (2019). For
the remaining two systems, we applied the same analysis procedure to the newly
available data and jointly analyzed them together with the previously existing
data set.

The data were processed using the XMMSAS v13.5 package and the corresponding
calibration database. The standard event screening chains were run to extract
lists of valid events for the three detectors of the European Photon Imaging
Camera (EPIC). We used the mos-filter and pn-filter tasks to filter out time
periods affected by strong soft proton flares. We extracted count images from
the clean event files in the [0.7-1.2] keV band, which maximizes the source-
to-background ratio and allows us to minimize the systematics associated with
background subtraction. To model the high-energy particle background, we used
a large collection of filter-wheel-closed observations available in the
calibration database, and rescaled the normalization of the closed data such
that the high-energy count rate matches the count rate measured in the corners
of the detectors, which are located outside of the instrument’s field of view
(FOV) and thus do not record any incoming photon. We also attempted to model
the residual soft proton background by measuring the ratio of the high-energy
count rates inside and outside the FOV, which can be related to the expected
soft proton rate (Leccardi & Molendi, 2008). This relation was calibrated
using a large set of blank-sky pointings and presented in Appendix B of
Ghirardini et al. (2018).

For each X-COP object, the count images from the various individual pointings
were co-added to create a mosaic combining all the available data. The same
procedure was repeated for the individual exposure maps and non X-ray
background maps. We masked the detected point sources and used the azimuthal
median technique (Eckert et al., 2015) to remove the contribution of
unresolved gas clumps that might bias the measured X-ray density (Nagai & Lau,
2011). To compute the median surface brightness profile, we created adaptively
binned surface brightness maps using Voronoi tessellation (Cappellari & Copin,
2003), and measured the median surface brightness value in concentric annuli
centered on the surface brightness peak.

To determine the plasma temperature, we accumulated X-ray spectra in
concentric annuli and extracted the appropriate response files and non X-ray
background spectra using the Extended Source Analysis Software (ESAS) package
(Snowden et al., 2008) available within XMMSAS. We followed the XSPEC spectral
fitting procedure outlined in Eckert et al. (2014), in which the particle
background is fitted with a phenomenological model and adjusted to the corners
of the detectors, and the sky background is described by a three-component
model (cosmic X-ray background, Galactic halo, and local hot bubble; McCammon
et al., 2002). The free parameters of the sky background model are fitted to
the spectrum of a background region located far outside of the cluster
($>2R_{500}$) and then rescaled appropriately to the source region. The source
is described by a single-temperature thin-plasma emission model using the APEC
code (Smith et al., 2001) with the temperature, normalization and metal
abundance left free while fitting. Finally, the source model is modified by
the Galactic absorption along the line of sight, with the Galactic column
density fixed to the HI column density estimated by the LAB survey (Kalberla
et al., 2005).

###  2.3 Planck data analysis

The pressure profiles are recovered from the _Planck_ survey data. We made use
of a custom all-sky SZ map constructed in a very similar way as the public map
of the full _Planck_ survey delivered by the Planck Collaboration (2016). Both
maps were obtained from an internal linear combination of the _Planck_
frequency maps, using the modified internal linear combination algorithm
(MILCA) method (Hurier et al., 2013). Our SZ map has a spatial resolution of 7
arcmin full width half maximum and by construction bears intrinsic correlated
noise.

The SZ signal is provided in unit-less Comptonization parameter, integrating
the gas pressure along the line of sight,

$y=\frac{\sigma_{\rm T}}{m_{\rm e}c^{2}}\int_{\ell}P(\ell)\,d\ell.$ (1)

The $y$-profiles and 3D pressure profiles are obtained following the method
presented in Planck Collaboration V (2013), and subsequently used in studies
by the X-COP collaboration (e.g., Tchernin et al., 2016; Ghirardini et al.,
2019; Eckert et al., 2019). We recall here the main steps of the $y$ and
pressure profile computation.

$y$-profiles are computed over local maps of size
$20\times\theta_{\textrm{500}}$ on a side, centered on the X-COP cluster
position and extracted from the all-sky SZ map. The patches are over-sampled
with respect to the _Planck_ SZ map pixel. The induced correlation between the
points of the $y$-profile is accounted for and propagated through the
covariance matrix of the profile. We performed an azimuthal mean of the pixel
fluxes falling in each radial bin, and subtracted a background residual offset
estimated from radii larger than $5\times R_{\textrm{500}}$ safe from any
cluster signal. Obvious positive or negative (arising from the linear
combination of frequencies) point sources were manually masked. We then
performed a pixel clipping of the local map using a $2.5\sigma$ criterion with
respect to the mean of the flux beyond $5\times R_{\textrm{500}}$. The noise
of the _Planck_ SZ map can be considered as Gaussian and correlated. For each
cluster we characterized the statistical properties of the local noise over
the individual patches in regions excluding the cluster emission, that is,
$\theta>5\times R_{\textrm{500}}$. From the derived noise power spectrum, we
draw a thousand realizations of the noise over each patch. Each realization is
folded in the procedure to extract the $y$-profile and then used to compute
the profile covariance matrix. The latter hence propagates the intrinsic
properties of the noise in the SZ _Planck_ map as well as the correlation
introduced by our choice of pixel size.

To recover the 3D pressure profiles, we optimally binned the $y$-profiles. We
assumed spherical symmetry and applied a deconvolution from the Gaussian
_Planck_ point spread function (PSF) and a geometrical deprojection following
the method detailed by Croston et al. (2006). Uncertainties on the $y$-profile
are propagated through the covariance matrix by generating $10,000$
realizations of the noise profile. Each are co-added to the $y$-profile,
individually deconvolved and deprojected, and used to compute the covariance
matrix of the 3D pressure profile.

###  2.4 Stellar mass

####  2.4.1 BCG

For the BCGs of A644, A1795, A2029, A2142, and A2319, we used stellar mass
distributions from Loubser et al. (2020). The stellar masses were derived from
$r$-band imaging obtained on the Canada-France-Hawaii Telescope (CFHT; see
Sand et al. 2012). The observations and corrections (e.g., extinction) are
described in Loubser et al. (2020), and we only briefly summarize the BCG
stellar mass modeling here.

We used the multi-Gaussian expansion (MGE) method, as implemented by
Cappellari (2002), to obtain the stellar mass distributions from the $r$-band
imaging. The PSF-convolved MGE models extend beyond the effective radii
($R_{e}$) for the BCGs. In the case of A1795, A2142 and A2319, we masked some
foreground or background features on the outer edges of the images.

The surface brightness obtained from the MGE method was converted to surface
brightness density, and the deprojection from surface density to intrinsic
density was performed for the axisymmetric case using the Jeans anisotropic
method (Cappellari 2008) assuming a constant stellar mass-to-light ratio,
$\Upsilon_{*\rm DYN}$. We assumed an inclination of $i$ = 90$\degr$ (see Smith
et al. 2017; Fasano et al. 2010; van der Marel 1991). In addition to the
stellar mass component, the BCG mass models include a central mass component
for a supermassive black hole, and a DM halo that follows a spherically
symmetric NFW profile. The DM mass ($M_{\rm DM}$) within $R_{200}$ was
approximated from weak lensing results (Herbonnet et al., 2020), resulting in
only two free parameters: (i) the stellar velocity anisotropy ($\beta_{z}$)
and (ii) the stellar mass-to-light ratio ($\Upsilon_{*\rm DYN}$). The BCG mass
model was compared to the observed kinematic profiles and the two parameters
$\beta_{z}$ and $\Upsilon_{*\rm DYN}$ were adjusted until the predicted
stellar kinematics best matched the spectroscopic observations (as described
in Loubser et al. 2020), placing constraints on those parameters. We use this
best fitting stellar mass-to-light ratio to derive the enclosed stellar mass
of the BCG within spheres of different radii.

For A644 and A2319, we did not have weak lensing masses available from
Herbonnet et al. (2020) but used the average DM distribution (of the 23 other
clusters in Loubser et al. 2020) to calculate the effect of a DM mass
component on the best fitting stellar mass-to-light ratio ($\Upsilon_{*\rm
DYN}$). When a DM mass component is included in the dynamical mass models for
the BCGs, the best fitting $\Upsilon_{*\rm DYN}$ decrease by 8.3 $\pm$ 2.9 per
cent on average.

Our surface brightness measurement are accurate at the level of $<0.25$ mag
arcsec-2, and our stellar kinematic measurement errors are $\pm 4.7$ per cent
on average (see Loubser et al. 2018). We also find that our dynamical modeling
is robust against the DM distribution profile that we assume, as well as
against assumptions regarding the black hole mass component. However, we find
that the observational errors on $M_{\rm DM}$ and $R_{200}$ (obtained from the
weak lensing results by Herbonnet et al. 2020) are by far the dominant
uncertainty. We derived the errors on the best fitting parameters by
incorporating the 1$\sigma$ errors on the weak lensing masses as well as a
$\pm$ 10 per cent uncertainty on the calculated value for the concentration
parameter. The uncertainties on the stellar mass profiles are therefore
constant with radius.

####  2.4.2 Satellite galaxies

The stellar mass content contained in the satellite population of the same
five clusters (A644, A1795, A2029, A2142, and A2319), is probed using $g$\-
and $r$-band imaging data obtained with MegaCam at the CFHT, and additional
$u$\- and $i$-band photometry taken with the Wide Field Camera at the Isaac
Newton Telescope. Aperture magnitude limits (5-$\sigma$)are typically 24.3,
24.8, 24.2, and 23.3 in the $ugri$-filters, respectively, when measured on
PSF-homogenized images using Gaussian weight functions. The stellar mass-to-
light ratios that form the basis of the stellar mass measurements are obtained
using spectral energy distribution fitting of these aperture fluxes. For this
we assume the stellar population libraries from Bruzual & Charlot (2003), an
exponentially declining star-formation history, and a Chabrier (2003) initial
mass function. Further details are provided in van der Burg et al. (2015).

While initially all galaxies are assumed to be members of the cluster, we
clean the galaxy distribution from interlopers (i.e., foreground and
background galaxies) by performing a statistical subtraction of sources
observed in the COSMOS field (Muzzin et al., 2013). Details are provided in
van der Burg et al. (2015) and Ghizzardi et al. (2021). In summary, to make a
fair accounting, we redo the photometry in the COSMOS field, only utilizing
the $ugri$-filters. We also accounted for field-to-field (often called
“cosmic”) variance, given that the COSMOS field is relatively small. We
estimate this uncertainty using Moster et al. (2011), as detailed in van der
Burg et al. (2015).

##  3 Mass modeling

Here we present the scheme that we developed to recover HSE masses and the
various methods used to reconstruct the mass profiles. The framework presented
here implements several methods to estimate the underlying mass profile:
parametric mass models (Sect. 3.2), forward fitting with parametric functions
(Sect. 3.3), and nonparametric log-normal mixture reconstruction (Sect. 3.4).
We distribute our code jointly with this paper in the form of the publicly
available Python package hydromass 222https://github.com/domeckert/hydromass.

###  3.1 Framework

The gravitational field of galaxy clusters can be estimated from the measured
thermodynamic profiles under the assumption that the gas is in HSE within the
potential well of the underlying halo. We assume that the ICM is fully
thermalized and that the gas pressure locally balances the gravitational
force,

$\frac{dP_{\rm gas}}{dr}=-\rho_{\rm gas}\frac{GM_{\rm tot}(<r)}{r^{2}},$ (2)

with $\rho_{\rm gas},P_{\rm gas}$ the gas density and pressure, respectively,
and $M_{\rm tot}(<r)$ the total (Newtonian) mass enclosed within a radius $r$.
To take advantage of the sum of information provided by the combination of
X-ray and SZ data, we performed joint fits to the two data sets by
constructing 3D model profiles for the gas density, temperature, and pressure,
which are related to one another through the ideal gas equation of state,

$P_{\rm gas}=\frac{k}{\mu m_{p}}\rho_{\rm gas}T.$ (3)

Solving the HSE equation requires a good control over the gradient of the
thermodynamic quantities, which can be highly sensitive to statistical
fluctuations. We use three complementary approaches to model the radial
temperature and pressure profiles and reconstruct $M_{HSE}$, which we describe
in Sects. 3.2 through 3.4. We integrate the three methods within a common
Bayesian fitting framework. The X-ray (surface brightness and spectroscopic
temperature) and SZ (electron pressure) observables are fitted jointly to
optimize the model parameters. We adopt a forward modeling approach in which
the model 3D quantities are projected along the line of sight and compared to
the projected data.

The gas density profile is modeled using the nonparametric multi-scale
approach implemented in Eckert et al. (2020). Namely, the 3D X-ray emissivity
profile is described as the linear combination of a large number $N_{K}$ of
King profiles,

$\epsilon(r)=\sum_{k=1}^{N_{K}}\alpha_{k}\Phi_{k}(r),$ (4)

with $\\{\Phi_{k}\\}_{k=1}^{N_{K}}$ the dictionary of input King functions
(King, 1962),

$\Phi_{k}(r)=\left(1+\frac{r}{r_{c,k}}^{2}\right)^{-3\beta_{k}},$ (5)

and $\\{\alpha_{k}\\}_{k=1}^{N_{K}}$ their multiplicative coefficients.
Following Eckert et al. (2020) the King function parameters are set adaptively
on a grid of values for the parameters $r_{c}$ and $\beta$, with one value of
$r_{c}$ for every set of 4 surface brightness points and 10 values of $\beta$
sampling the range 0.6-3. Our final dictionary contains $N_{K}=5N/2$
functions, with $N$ the number of surface brightness points. The 3D model is
projected onto the line of sight and convolved with the instrumental PSF to
create a model for the observed surface brightness. The coefficients
$\alpha_{k}$ are then fitted to the data. More details on the method and an
extensive validation with numerical simulations are provided in Eckert et al.
(2020).

The observed spectroscopic temperature is the line-of-sight average of the 3D
temperature profile. The average temperature is weighted by the local
emissivity $\epsilon\propto n_{e}^{2}$. On top of that, the temperatures
obtained through single-temperature fits to multiphase spectra tend to be
lower than mass-weighted temperatures because of the response of X-ray
instruments. Mazzotta et al. (2004) showed that for temperatures of 3 keV and
above, a condition that is satisfied by all X-COP clusters, the average
spectroscopic temperature can be written as

$T_{spec}(r)=\frac{\int T_{3D}\,n_{e}^{2}T_{3D}^{-3/4}d\ell}{\int
n_{e}^{2}T_{3D}^{-3/4}d\ell}.$ (6)

For a given model 3D temperature profile, the profile is projected onto the
line of sight using Eq. 6, convolved with the instrumental PSF and adjusted
onto the measured temperatures. The model 3D pressure profile is fitted
jointly to the SZ pressure profile, taking into account the correlated noise
of the _Planck_ SZ signal (Planck Collaboration V, 2013). Therefore, our joint
fitting procedure takes advantage of both the X-ray and SZ information by
constructing a joint likelihood comparing the 3D model pressure to the SZ
pressure and the spectroscopic-like temperature profile to the X-ray
temperatures (see Eq. 2 of Ettori et al., 2019).

###  3.2 Mass model

If a parametric form for the mass model can be defined, the HSE equation (Eq.
2) can be integrated to predict the expected pressure and temperature
profiles,

$P_{3D}(r)=P_{0}+\int_{r}^{r_{0}}\frac{\rho_{\rm gas}GM_{\rm
mod}(r^{\prime})}{r^{\prime 2}}\,dr^{\prime},$ (7)

with $M_{\rm mod}(r)$ the model density profile integrated over the volume out
to radius $r$, $r_{0}$ the outermost radius out to which the pressure can be
measured, and $P_{0}=P(r_{0})$ the integration constant at the outer boundary.
This method is based on the “backward” approach introduced by Ettori et al.
(2010), with the notable difference that the integration constant $P_{0}$ is
left free while fitting. We set a Gaussian prior on $P_{0}$ with mean and
standard deviation set to the outermost pressure value and its 1-sigma
uncertainty.

In Ettori et al. (2019) (hereafter E19) we applied a similar technique to the
X-COP sample and compared five different model profiles. We found that models
with a central cusp are usually preferred by the data over cored models. For
this reason, here we focus on the NFW model,

$\rho_{NFW}(r)=\frac{\rho_{0}}{\left(\frac{r}{r_{s}}\right)\left(1+\frac{r}{r_{s}}\right)^{2}},$
(8)

and the Einasto model,

$\rho_{Einasto}(r)=\rho_{s}\exp\left[-\frac{2}{\alpha}\left(\left(\frac{r}{r_{s}}\right)^{\alpha}-1\right)\right].$
(9)

Instead of the scale radius $r_{s}$ and density normalization $\rho_{s}$ we
optimize for the overdensity radius $R_{200}$ and the concentration
$c_{200}=R_{200}/r_{s}$. We expand on the work of E19 in several ways. First,
by correcting for the telescope’s PSF (see Sect. 3.5), which allows us to
accurately trace the DM profile in the inner regions. Second, in the Einasto
case we leave the Einasto index $\alpha$ free, whereas previously its value
was fixed to $1/\alpha\equiv\mu=5$ (Navarro et al., 2004; Mamon & Łokas,
2005). We set weakly informative independent Gaussian priors on the model
parameters of the NFW and Einasto models to make sure the fitting procedure
quickly converges toward a reasonable solution whilst not biasing the
resulting posteriors. In Table 1 we give the details of the priors on the NFW
parameters; a similar description of the priors on the Einasto model is
provided in Paper II. In case the available data are of lower quality, our
framework also allows the user to set informative priors on the model
parameters, for example by including prior information on the mass-
concentration relation.

Table 1: Normal priors on the NFW fit parameters. Here $P_{m}$ and $dP_{m}$ denote the outermost SZ pressure value and its error. Parameter | Mean | $\sigma$ | Min | Max
---|---|---|---|---
$R_{200}$ [kpc] | 2000 | 1000 | 300 | 4000
$c_{200}$ | 4 | 4 | 0 | 15
$P_{0}$ | $P_{m}$ | $dP_{m}$ | $P_{m}-2dP_{m}$ | $P_{m}+2dP_{m}$

Whenever information on the baryonic mass profile (stellar and gas) is
available, we perform an additional reconstruction in which we attempt to
model the gravitational field as the sum of baryonic and DM profiles,

$\rho_{\rm tot}(r)=\rho_{DM}+\rho_{\rm
gas}+\rho_{\star,BCG}+\rho_{\star,sat}.$ (10)

In this case, the DM profile $\rho_{DM}$ is described by a mass model (NFW or
Einasto) and the contribution of the baryonic components is provided as input.

###  3.3 Parametric forward model

In this case, we introduce a parametric form for the pressure profile and
attempt to reconstruct the shape of the mass profile without prior assumption
on its shape. This method is similar to the approach introduced by Vikhlinin
et al. (2006). We model the pressure profile as a generalized NFW (gNFW)
profile (Nagai et al., 2007),

$P_{forw}(r)=\frac{P_{0}}{(r/r_{s})^{\gamma}\left(1+(r/r_{s})^{\alpha}\right)^{\frac{\beta-\gamma}{\alpha}}},$
(11)

with $P_{0}$, $r_{s}$, $\alpha$, $\beta,$ and $\gamma$ as free parameters.
This functional form was found to provide a good representation of gas
pressure profiles determined from both X-ray and SZ observations (Arnaud et
al., 2010; Planck Collaboration V, 2013; Sayers et al., 2013; Bourdin et al.,
2017; Ghirardini et al., 2019; Pointecouteau et al., 2021). The pressure
gradient can be computed analytically from the reconstructed functional form
to recover the mass profile from Eq. 2.

###  3.4 Nonparametric log-normal mixture reconstruction

Here we introduce a nonparametric method for temperature deprojection and HSE
mass reconstruction. Our method attempts at the same time to make minimal
assumptions on the shape of the mass profile and limit the fluctuations
introduced by the deprojection process and the computation of the pressure
gradient. To this aim, we describe the 3D temperature profile as a linear
combination of a large number $N_{g}$ of log-normal functions,

$T_{NP}(r)=\sum_{i=1}^{N_{g}}G_{i}\frac{1}{\sqrt{2\pi\sigma_{i}^{2}}}\exp\left(-\frac{(\ln(r)-\ln(\mu_{i}))^{2}}{2\sigma_{i}^{2}}\right).$
(12)

The mean values $\\{\mu_{i}\\}_{i=1}^{N_{g}}$ and standard deviations
$\\{\sigma_{i}\\}_{i=1}^{N_{g}}$ of the Gaussians are set a priori to allow as
much freedom as possible to the fitted function, whilst at the same time
preventing small, unphysical fluctuations induced by the statistical nature of
the problem. For large values of $N_{g}$ the model is essentially independent
of the choice of $N_{g}$ and $\mu_{i}$, whereas the standard deviations
$\sigma_{i}$ act as effective smoothing scales. We implement the model with
$N_{g}=200$ and the values of $\mu_{i}$ logarithmically spaced from the
innermost to the outermost data points. We set the values of $\sigma_{i}$ to
the width of the nearest spectroscopic annulus or SZ radial bin, such that
fluctuations on scales smaller than the bin size are suppressed and the model
follows closely the data on larger scales.

Once the number of log-normal functions $N_{g}$, the mean values $\mu_{i}$ and
the standard deviations $\sigma_{i}$ are adaptively set, the model can be
projected onto the line of sight including spectroscopic-like weights (Eq. 6)
and the normalizations $\\{G_{i}\\}_{i=1}^{N_{g}}$ can be fit to the
spectroscopic X-ray temperatures and SZ pressure profile. The mass profile is
then reconstructed by combining the 3D temperature and density profiles,

$M(<r)=-\frac{rkT_{3D}(r)}{G\mu m_{p}}\left(\frac{\partial\log
T_{3D}}{\partial\log r}+\frac{\partial\log n_{\rm gas}}{\partial\log
r}\right).$ (13)

The temperature gradient can be computed analytically from Eq. 12,

$\frac{\partial T_{3D}}{\partial
r}=\sum_{i=1}^{N_{g}}G_{i}\frac{1}{\sqrt{2\pi}\sigma_{i}^{3}}\exp\left(-\frac{(\ln(r)-\ln(\mu_{i}))^{2}}{2\sigma_{i}^{2}}\right)\frac{\ln(\mu_{i})-\ln(r)}{r}.$
(14)

Similarly, since the King functions $\\{\Phi_{k}\\}_{k=1}^{N_{K}}$ are
analytical functions, the gas density gradient can be computed analytically
from Eq. 4,

$\frac{\partial n_{\rm gas}}{\partial
r}=\sum_{k=1}^{N_{K}}\alpha_{k}\frac{\partial\Phi_{k}}{\partial r}.$ (15)

###  3.5 Optimization

Since our models can contain several hundred parameters, we require the use of
a statistical sampler that is suitable for high-dimensional optimization
problems. We use the No U-Turn Sampler (NUTS; Hoffman & Gelman, 2014), a
Hamiltonian Monte Carlo (HMC) algorithm implemented in the PyMC3 package
(Salvatier et al., 2016). The convergence of HMC algorithms is typically much
faster than that of traditional Markov chain Monte Carlo samplers thanks to
the use of gradient information. At the same time, NUTS includes advanced
self-tuning features, such that the sampler requires minimal input from the
user.

For each cluster and model, we start by searching for the maximum likelihood
using a gradient descent method, and then run four parallel HMC chains with
2,000 tuning steps and 2,000 samples each. The four chains are then combined
to measure the posterior distributions of the parameters. Running multiple
chains in parallel improves the computational efficiency of the code and
reduces the risk that the posterior be drawn from a single chain that has not
fully converged.

In Fig. 1 we show an example of reconstruction for A1795 ($z=0.0622$). The
results obtained from the NFW and Einasto mass models are compared to the
results of the parametric forward approach (hereafter labeled as “Forward”)
and the nonparametric log-normal mixture method (hereafter NP). All the
methods provide similar results across most of the radial range, although
differences arise primarily in the innermost 100 kpc. Since the NP method is
designed to trace the observed values as closely as possible, it cannot be
extrapolated outside of the fitted range. Thus, for the remainder of the paper
we evaluate it at the radii of the observed data and present them in the form
of data points rather than continuous models. In the bottom-right panel we
show the posterior distributions for the parameters of the NFW profile and the
correlations between the parameters. Given the excellent statistical quality
of the X-COP data, the posterior distributions are much narrower than the
prior and clearly pull away from it, indicating that our data are able to
break the degeneracy between NFW mass and concentration. To help convergence
of the Einasto model, we optimize for the inverse of the Einasto index,
$\mu=1/\alpha$ (see Paper II). The details of the Einasto fitting procedure,
the adopted priors, and the resulting best-fit parameters are provided in
Paper II.

Figure 1:  Mass reconstruction results for A1795. _Top left:_ Electron density
profile reconstructed with the multi-scale method used here (Eckert et al.,
2020) compared to the L1 regularization results presented in Ghirardini et al.
(2019) and publicly released by the X-COP team. _Top right:_ Nonparametric
reconstruction of the 3D temperature profile (green) compared to the
spectroscopic X-ray measurements (red) and the 3D temperature profile obtained
by dividing the SZ pressure by the X-ray density (orange). See Sect. 3.4 for
details. The blue curve shows the projected, spectroscopically weighted, and
PSF convolved temperature profile, to be compared to the red data points.
_Middle left:_ Electron pressure profile measured by _Planck_ and _XMM-Newton_
compared to the reconstructions obtained from the NFW and Einasto mass models
(Sect. 3.2), the parametric forward approach (labeled “Forward”; see Sect.
3.3), and the NP reconstruction (Sect. 3.4). The dotted vertical lines show
the location of $R_{500}$ and $R_{200}$ from the NFW fit. _Middle right:_
Output mass profiles obtained from the four individual reconstructions (NFW,
Einasto, Forward, and NP). The data points show the NP results evaluated at
the radius of the X-ray spectroscopic measurements (labeled “NP Spectral
bins,” blue) and the SZ pressure data (“NP SZ bins,” red). _Bottom left:_ Same
as previous but for the gas mass fraction. The gray shaded area shows the
cosmic baryon fraction $\Omega_{b}/\Omega_{m}$ in Planck Collaboration XIII
(2016) cosmology. _Bottom right:_ Posterior distributions and correlations for
the parameters of the NFW profile.

##  4 Validation using mock data

To validate our numerical framework, we created synthetic observations of an
idealized cluster in HSE and ran our entire pipeline on the mock data. In the
following, we describe the various steps of this validation process and
present the accuracy of our reconstruction technique.

Figure 2: Result of temperature deprojection (left) and HSE mass
reconstruction (right) reconstructed from a set of ten mock _XMM-Newton_
observations of a synthetic NFW cluster (see Sect. 4). The curves show the
mean and standard deviation of the ten realizations for the various
reconstruction techniques (NFW, green; Einasto, cyan; and NP, red). The bottom
panels show the ratio of reconstructed profiles to true input profiles.

###  4.1 Synthetic galaxy cluster simulations

We consider a fiducial DM halo with a mass profile described by a generic NFW
profile (Eq. 8) with a concentration $c_{NFW}=4$ and an overdensity radius
$R_{200}=2000$ kpc at a redshift $z=0.2$, such that the simulated cluster fits
well within the _XMM-Newton_ FOV. The halo is filled with an ICM in HSE within
the input NFW potential well. The input gas density profile follows the
“universal” average X-COP gas density profile (see Fig. 3 and Table 3 of
Ghirardini et al., 2019), which was fitted to a collection of individual
nonparametric gas density profiles using a single common functional form
adapted from Vikhlinin et al. (2006) and jointly accounting for the intrinsic
scatter among the cluster population.

Given the input NFW mass profile and the universal gas density profile, we
used Eq. 7 to compute the corresponding thermal pressure profile. The profiles
were computed out to $5\times R_{500}$, which is well beyond the range covered
by our data to avoid uncertainties linked with the integration constant
$P_{0}$. We then generated a 3D box inside which we set the gas density and
temperature of the ICM inside each cell following spherical symmetry. The cell
size was set to 2 kpc, which is substantially smaller than the _XMM-Newton_
resolution. Inside each box cell, we used the APEC model (Smith et al., 2001)
as implemented in 3ML (Vianello et al., 2015; Burgess et al.,
2021)333https://threeml.readthedocs.io/en/latest/notebooks/APEC_doc.html to
generate model spectra in each cell. A constant metal abundance of
$0.3Z_{\odot}$ relative to the Anders & Grevesse (1989) Solar abundance table
was used in all cells. The spectra were absorbed by Galactic absorption with
$N_{H}=5\times 10^{20}$ cm-2. Finally, the 3D model spectra were projected
onto the line of sight. As a result, a spectrum made of the linear combination
of individual APEC spectra is associated with each individual pixel of a
projected 2D map.

###  4.2 Mock data creation

From the projected model spectra and images, we generated mock _XMM-Newton_
observations of our synthetic cluster including a wide variety of instrumental
effects: effective area, energy redistribution, vignetting, PSF smearing, and
a realistic background model. Our mock data generation tool is in the form of
a Python package, which is made publicly available together with this
paper444https://github.com/domeckert/xmm_simulator. The tool reads the on-axis
effective area of each of the three _XMM-Newton_ instruments from the
calibration database, and combines the on-axis area with the telescope’s
energy-dependent vignetting curve to generate a local effective area curve at
each position in the FOV. The projected model spectra described above are
multiplied by the local effective area curve and convolved with the
instrument’s redistribution matrix. The spectra are then convolved with the
instrumental PSF. Finally, we generate a mock source event list from Poisson
realizations of the model spectra.

On top of the source spectrum, we add a realistic background spectrum made of
a vignetted component for the sky background and an un-vignetted component for
the particle-induced background. The un-vignetted component is drawn out of
the collection of filter-wheel-closed observations available in the
calibration database. The vignetted component is made of a three-component
model (see Sect. 2.2) consisting of an absorbed power law for the cosmic X-ray
background and two APEC models describing the foregrounds. The normalization
of these components is set to be representative of a typical extragalactic
field. A Poisson realization of the sky background is then generated and
merged with the source event list to create a total event file including our
source, the astrophysical background and the cosmic-ray induced background.
Images and spectra can then be extracted from the final event list and
analyzed in the same way as real observations.

###  4.3 Code performance

We performed a set of ten individual 50 ks mock observations of our synthetic
cluster using the procedure described above and analyzed the simulated data
using the exact same analysis pipeline as for the actual data. Namely, for
each realization, an image of the synthetic cluster in the [0.7-1.2] keV band
is extracted from the simulated events. From the resulting image, we extract a
background subtracted surface brightness profile in the same way as for real
observations. We extract mock X-ray spectra in ten logarithmically spaced
annuli and a background region located at $R>2R_{500}$ from the cluster
center, and we fit the spectra using a single-temperature APEC model and the
sum of vignetted and un-vignetted background components. We then run our mass
reconstruction code on the surface brightness and spectroscopic temperature
profiles. In each case, we fit the data using the NFW and Einasto mass models
as well as the nonparametric log-normal mixture model.

The performance of our code on the mock data is shown in Fig. 2. In the left-
hand panel we show the mean and dispersion of the ten deprojected 3D
temperature profiles compared with the input profile. The nonparametric
profile was evaluated at the average radii of the spectroscopic bins. We can
see that all three methods are able to trace the true profile with minimal
bias (i.e., less than 3% at each radius). In particular, the nonparametric
method provides a good description of the profile with a scatter of $\sim
5\%$, even though it makes no assumption on the shape of the mass profile. The
NFW model traces the true profile very closely, which is expected since the
true input mass profile was assumed to be an NFW. The Einasto profile also
provides an excellent description of the data, although it overestimates the
temperature by 3% on average.

The comparison between the mean of the fitted mass profiles and the true
profile can be seen in the right-hand panel of Fig. 2, with the residuals
shown in the bottom panel. Again, as expected the NFW reconstruction matches
very well the input profile, with a bias of at most 3% at all radii but
consistent within the uncertainties. The Einasto and nonparametric profiles
closely trace the true profile, with the Einasto exceeding the true mass by
about 6% on average, which can be explained by the application of a different
parametric model compared to the true one. The nonparametric technique follows
the true profile with minimal bias (¡5%) but somewhat higher uncertainties
given the freedom allowed to the model. The uncertainties are higher at the
inner and outer edges of the profile, since the temperature gradient cannot be
directly constrained there.

Overall, the tests presented here show that our analysis pipeline and HSE
reconstruction framework are able to reconstruct the properties of the
gravitational field in our synthetic cluster with a bias of at most 3% when
adopting the NFW method. The uncertainties associated with the reconstruction
technique are therefore subdominant with respect to other sources of
systematic uncertainty such as hydrostatic bias (10-20%; e.g., Rasia et al.,
2006; Lau et al., 2009), temperature inhomogeneities ($\sim 10\%$, Rasia et
al., 2014) and effective area calibration ($15\%$, Schellenberger et al.,
2015).

##  5 Results

Table 2:  Sample properties and results of NFW fits. The columns labeled “DM” refer to the NFW fits to the DM only for the five systems with complete information on the baryonic components (see Sect. 5.4), whereas the “TOT” subscripts refer to the fits of the total gravitational field with an NFW model (see Fig. 3). In this case, the quoted total mass is computed from the sum of DM and baryonic density profiles. Cluster | z | $N_{H}$ | $R_{500,TOT}$ | $M_{500,TOT}$ | $M_{200,TOT}$ | $c_{200,TOT}$ | $M_{200,DM}$ | $c_{200,DM}$
---|---|---|---|---|---|---|---|---
|  |  $10^{20}$ cm-2 | kpc | $10^{14}M_{\odot}$ | $10^{14}M_{\odot}$ |  | $10^{14}M_{\odot}$ |
A85 | 0.0555 | 2.8 | $1214_{-4}^{+4}$ | $5.03_{-0.05}^{+0.05}$ | $7.04_{-0.07}^{+0.08}$ | $4.708_{-0.044}^{+0.042}$ | - | -
A644 | 0.0704 | 7.5 | $1398_{-16}^{+16}$ | $7.80_{-0.27}^{+0.26}$ | $10.97_{-0.43}^{+0.45}$ | $4.597_{-0.156}^{+0.163}$ | $10.98_{-0.41}^{+0.45}$ | $4.458_{-0.169}^{+0.171}$
A1644 | 0.0473 | 4.1 | $1031_{-10}^{+11}$ | $3.06_{-0.09}^{+0.09}$ | $4.93_{-0.20}^{+0.21}$ | $2.548_{-0.135}^{+0.140}$ | - | -
A1795 | 0.0622 | 1.2 | $1160_{-6}^{+6}$ | $4.42_{-0.07}^{+0.07}$ | $6.04_{-0.11}^{+0.11}$ | $5.454_{-0.076}^{+0.072}$ | $6.82_{-0.13}^{+0.13}$ | $4.510_{-0.072}^{+0.074}$
A2029 | 0.0766 | 3.2 | $1340_{-9}^{+9}$ | $6.91_{-0.14}^{+0.14}$ | $9.29_{-0.21}^{+0.22}$ | $5.929_{-0.101}^{+0.102}$ | $11.67_{-0.28}^{+0.29}$ | $4.185_{-0.091}^{+0.093}$
A2142 | 0.09 | 3.8 | $1453_{-9}^{+9}$ | $8.93_{-0.16}^{+0.17}$ | $13.33_{-0.30}^{+0.32}$ | $3.440_{-0.070}^{+0.076}$ | $13.94_{-0.31}^{+0.30}$ | $3.046_{-0.074}^{+0.076}$
A2255 | 0.0809 | 2.5 | $1202_{-16}^{+17}$ | $5.01_{-0.20}^{+0.22}$ | $8.94_{-0.56}^{+0.66}$ | $1.800_{-0.180}^{+0.193}$ | - | -
A2319 | 0.0557 | 3.2 | $1424_{-4}^{+5}$ | $8.12_{-0.08}^{+0.08}$ | $11.61_{-0.14}^{+0.14}$ | $4.233_{-0.070}^{+0.068}$ | $12.16_{-0.14}^{+0.15}$ | $4.098_{-0.079}^{+0.079}$
A3158 | 0.059 | 1.4 | $1146_{-11}^{+12}$ | $4.25_{-0.13}^{+0.14}$ | $6.46_{-0.25}^{+0.28}$ | $3.178_{-0.146}^{+0.145}$ | - | -
A3266 | 0.0589 | 1.6 | $1381_{-10}^{+10}$ | $7.43_{-0.15}^{+0.16}$ | $12.22_{-0.36}^{+0.41}$ | $2.352_{-0.098}^{+0.097}$ | - | -
RXC1825 | 0.065 | 9.4 | $1109_{-8}^{+9}$ | $3.88_{-0.08}^{+0.09}$ | $5.66_{-0.17}^{+0.18}$ | $3.807_{-0.159}^{+0.165}$ | - | -
Zw1215 | 0.0766 | 1.7 | $1346_{-21}^{+21}$ | $7.00_{-0.32}^{+0.35}$ | $11.44_{-0.70}^{+0.76}$ | $2.413_{-0.148}^{+0.155}$ | - | -

###  5.1 Total density profiles

We used our hydromass Python package introduced in Sect. 3 to analyze the
gravitational field of the 12 X-COP clusters. In each case, we ran the
reconstruction with the two mass models (NFW and Einasto), the parametric
forward method, and the nonparametric log-normal mixture technique, similar to
the example shown in Fig. 1 for A1795. In the case of the NP reconstruction,
the flexibility of the model makes it highly uncertain away from the fitted
data points. Therefore, we restrict specifically to the radii corresponding to
the existing spectral and SZ data points to avoid extrapolating. The masses
estimated from the NFW fits as well as the basic properties of the sample are
shown in Table 2. The results obtained with the Einasto model are presented
and discussed in Paper II.

We started by comparing the output of our new framework to the masses
previously published in E19 in the case of the “standard” NFW reconstruction.
As explained in Sect. 3.2, the main difference with respect to our previous
analysis is the implementation of a modeling scheme for the PSF, both for the
surface brightness and the temperature profile. Correcting for PSF smearing
increases the concentration parameter of the NFW profile, which results in
slightly higher masses in the core and lower masses in the outskirts. The
average ratio between the E19 masses and the results presented here is given
in Table 3 for several overdensities. Our new masses are slightly higher
inside $R_{2500}$ and slightly lower inside $R_{200}$; the difference,
however, is always within 10%. The validation on mock data presented in Sect.
4 shows that our framework is able to recover the NFW shape very accurately,
and thus we conclude that the NFW concentrations given in E19 were slightly
underestimated.

Table 3: Ratio between the masses published in E19 ($M_{\rm E19}$) and our new results (labeled $M_{\rm hydromass}$) for several overdensity factors. For consistency, the masses estimated here were converted into the same cosmology as that used in E19. Overdensity | 2500 | 1000 | 500 | 200
---|---|---|---|---
mean($M_{\rm hydromass}/M_{\rm E19}$) | 1.08 | 1.01 | 0.96 | 0.93
median($M_{\rm hydromass}/M_{\rm E19}$) | 1.05 | 1.02 | 0.95 | 0.90

Figure 3: Comparison between the mass profiles reconstructed using mass models
(solid curves and shaded areas) and the nonparametric log-normal mixture
method (data points). The results obtained for the NFW and Einasto profile are
shown in the left and right panels, respectively. The data points show the
posterior NP mass estimates at the radii of the X-ray spectral data (circles)
and the SZ pressure data (asterisks). For clarity, all the profiles are scaled
by the NFW values for $R_{500}$ and $M_{500}$. The bottom panels show the
local ratio of NP to model.

In Fig. 3 we present a detailed comparison between the NFW and Einasto mass
models and the nonparametric reconstruction. All the masses and radii were
scaled by the NFW $M_{500}$ and $R_{500}$ values to visualize the differences
in the shape of the profiles. The choice of the NFW results as scaling factors
is motivated by the greater stability of the NFW masses compared to the other
methods. Conversely, to study the radial shape of the profiles, the results
obtained with the NP method are used as a benchmark of the local gravitational
field because the method is fully data driven. In the left- and right-hand
panels of the figure we show the reconstructed NFW and Einasto profiles,
respectively, for the entire X-COP sample. We immediately notice that both
parametric mass models trace the NP points closely, albeit with substantial
scatter in the core and the outskirts. The additional degree of freedom
afforded by the Einasto model allows it to trace the NP points more closely
than the NFW, as there appears to be more diversity in the profile shapes than
can be captured by the NFW model.

Figure 4: Average profile shapes. _Left:_ Median scaled mass profiles for the
four methods investigated here. The bottom panel shows the median ratio of NP
to model, with the error bars indicating the 68% interval. The ratios are
evaluated at the same radius but are slightly shifted for better readability.
_Right:_ NFW mass-concentration relation for our sample. The small shaded
ellipses indicate the $1\sigma$ contours of the fitted individual values for
the 12 X-COP clusters individually, whereas the solid red and blue ellipses
show the fitted sample mean and its $1\sigma$ and $2\sigma$ confidence
regions, respectively. The curves show the predictions of N-body simulations
(Ludlow et al., 2016; Diemer & Joyce, 2019; Ishiyama et al., 2021).

This visual impression is confirmed in a quantitative way in the left-hand
panel of Fig. 4, where we show the sample median profiles and 68% percentiles
for the various methods. The range of values in the cluster core is clearly
broader for Einasto than NFW. The forward-fitting method also closely traces
the NP profile, with the exception of the innermost regions, where the power-
law shape of the gNFW pressure profile leads to excess pressure, and hence
excess mass, compared to the NP data. In the left-hand panel of Fig. 4 we show
the median NP-to-model ratio across the sample, as well as the dispersion of
the model with respect to the reference NP value. The mass profiles estimated
with all the methods are in remarkable agreement, with the median masses in
the range $[0.05-2]R_{500}$ being within 3% of one another, even though the
various methods make vastly different underlying assumptions. On average, we
do not observe systematic deviations from an NFW shape, with the relative
difference between NFW and NP points being less than 10% at all radii and
consistent with the level of deviations found in the mock observations (see
Sect. 4). Similarly, the Einasto profile accurately reproduces the average
shape of the gravitational field in our systems. This study is therefore in
agreement with that of E19, which concluded that the mass profiles of X-COP
clusters are better represented by models that exhibit a rising density
profile in their inner regions, whereas models that include large central
cores (Burkert, isothermal sphere) are statistically disfavored. In all cases,
the scatter is minimal in the radial range [0.2-1]$R_{500}$, where all the
methods converge within less than 10% of one another. The difference between
the models is largest in the outer regions, where the scatter of the NFW
values relative to the NP points is nearly twice as large as for the other two
methods. This result shows that the properties of the gas in cluster outskirts
require more diversity in the DM profile shape than can be described by the
NFW model, although the validity of the HSE assumption at large radii is
unclear (see the discussion in Sect. 6.1). Our results agree with the
prediction of Neto et al. (2007), who showed that while the NFW profile
provides a good description of $\Lambda$CDM density profiles on average,
deviations from an NFW shape can be substantial in individual systems, with
average deviations on the order of $0.1$ dex. Conversely, the Einasto profile
traces the measured ICM properties closely even beyond $R_{500}$.

###  5.2 Concentration-mass relation

The relation between NFW concentration and mass is a clear prediction of the
$\Lambda$CDM paradigm (e.g., Duffy et al., 2008; Bhattacharya et al., 2013;
Dutton & Macciò, 2014; Diemer & Kravtsov, 2014). It arises from the
universality of the structure formation process through mergers and accretion.
Our precise measurements of halo concentrations therefore allow us to test the
$\Lambda$CDM framework.

Given that our sample spans a narrow mass range (only a factor of $\sim 2$),
it is not sensitive to the shallow slope of the NFW mass-concentration
relation. However, it is well suited to measure the mean and scatter of the
relation around $M_{200}=10^{15}M_{\odot}$. To this end, we performed a
Bayesian multivariate analysis of our sample, with the fractional log-normal
scatter in mass and concentration as free parameters. Namely, we describe our
sample as a set of values drawn from a distribution centered on the mean
values for $R_{200}$ and $c_{200}$ with free intrinsic log-normal scatter on
both axes. We set Gaussian priors on the two mean values, with the mean and
standard deviation of the priors set to the mean and scatter of the sample
values. A positive half-Cauchy prior with $\beta=0.5$ is set on the fractional
intrinsic scatter along both axes. Adopting instead a half-normal prior on the
scatter has little impact on the results. To take the strong intrinsic
correlation between mass and concentration in the NFW model into account, for
each cluster we compute the covariance matrix between the two parameters from
the output chains, and we set a global likelihood as the product of the
Gaussian multivariate likelihoods of each object including the covariance
matrix. We then sample the total likelihood using PyMC3.

The measured $1\sigma$ contours for all X-COP clusters in the
$M_{200}-c_{200}$ plane are shown in Fig. 4 together with the fitted mean
values. The strong intrinsic anticorrelation between mass and concentration
can be clearly seen on the plot. The posterior distributions for the
parameters of interest are shown in Fig. 5. The mass-concentration relation
for our sample is centered on $M_{200}=8.64_{-0.81}^{+0.91}\times
10^{14}M_{\odot}$ and $c_{200}=3.69_{-0.36}^{+0.39}$ with an intrinsic scatter
$\sigma_{\ln c_{200}}=0.37_{-0.07}^{+0.11}$.

Figure 5: Posterior distributions for the fitted concentration-mass relation,
with $R_{200}$ and $c_{200}$ the mean values of $R_{200}$ and NFW
concentration in our sample and $\sigma_{R_{200}}$ and $\sigma_{c}$ the
fractional log-normal intrinsic scatter of the points around the mean value.

We used the Python package Colossus (Diemer, 2018) to compare our results with
the predictions of various N-body simulation suites. The curves in the right-
hand panel of Fig. 4 show the predictions of three different sets of
simulations (Ludlow et al., 2016; Diemer & Joyce, 2019; Ishiyama et al.,
2021). The curves were calculated at the median redshift of the X-COP sample
($z=0.065$) for a proper comparison with the data. We find a remarkably good
agreement between our average values and the predictions of numerical
simulations, with all three simulations considered here agreeing with our
measurements within $1\sigma$. On the other hand, numerical simulations also
predict a large scatter in concentration at given mass, lower in more relaxed
systems. Neto et al. (2007) find that a log-normal distribution represents the
estimated concentrations in a given mass bin well, with a mean and a
dispersion that decrease at higher masses, with the dispersion on the order of
$\sigma_{\ln c}\sim 0.21$ for more relaxed systems. For the most massive
clusters, Bhattacharya et al. (2013) estimate $\sigma_{\ln c}=0.33$, in
excellent agreement with our measurements. Diemer & Kravtsov (2015) measure a
scatter of about 0.16 dex ($\sigma_{\ln c}=0.37$) at all redshifts, masses,
and for all mass definition, when the entire ensemble of halos is considered,
again in agreement with the X-COP data. Different halo collapse times, with
higher concentration associated with a halo assembled earlier, can account for
most of the measured scatter (e.g., Navarro et al., 1997; Ludlow et al.,
2016).

All together, the shape parameters retrieved here agree well with the
predictions of the $\Lambda$CDM framework, both for the average NFW
concentration and the scatter of the c-M relation. A discussion of profile
shape in the central regions of X-COP clusters as traced by the Einasto shape
parameter is provided in Paper II, where we present our measurements of
$\alpha$ and discuss in detail the implications of our measurements.

###  5.3 Breaking down into baryonic and DM components

Figure 6:  Decomposition of the gravitational field into DM and baryonic
components. _Left:_ Breakdown of the total differential density profile into
baryonic and DM components for the five systems for which stellar mass
profiles are available for both the BCG and the satellite galaxies. The
baryonic components are split into the stellar mass (BCG, dashed green;
satellites, dash-dotted cyan) and the gas mass (dashed red), whereas the
fitted Einasto DM profiles are shown as the solid magenta curves. _Right:_
Total cumulative baryon fraction as a function of radius. The total baryon
fraction and its uncertainty are shown as the solid curves and shaded areas,
where the dashed and dotted lines show the relative contribution of gas and
stars, respectively.

As stated in Sect. 2, for a subset of systems (A644, A1795, A2029, A2142, and
A2319) we have access to precise measurements of the stellar mass profiles of
the cluster’s BCG (Loubser et al., 2020) and of the cumulative stellar mass of
satellite galaxies (van der Burg et al., 2015). Therefore, with the exception
of a potential contribution of intracluster light (ICL) we have direct
measurements of the mass profiles of all the baryonic components (gas, BCG,
and satellites). For these systems, we conducted an additional analysis where
the cumulative mass of the baryonic components is added as a fixed component
to the mass profile at every radii. The remaining mass is then modeled as an
additional model component (NFW or Einasto) to constrain the shape of the
gravitational field induced by DM only. This procedure is described in detail
in Sect. 3.2.

Given that the stellar mass profile of the BCG can be estimated only within a
limited radial range out to $\sim 50$ kpc, we fitted the BCG stellar mass
using a Sérsic functional form, which was found to provide a good description
of the data,

$\rho_{BCG}(r)=\rho_{0}\exp\left\\{-b_{n}\left((r/r_{e})^{1/n}-1\right)\right\\}.$
(16)

Similarly, the cumulative stellar mass profiles of the satellite galaxies are
quite noisy and not smooth, which is difficult to treat in the context of this
analysis. van der Burg et al. (2015) found that the stellar mass density
profiles of satellite galaxies can be well described by a generalized NFW
profile (see Eq. 11). Therefore, we fitted the stellar mass profiles beyond
the BCG ($R>100$ kpc) with a gNFW profile with shape parameters fixed to van
der Burg et al. (2015), and used the resulting functional form as input for
our mass reconstruction code.

The results of the mass reconstruction separating the baryonic and DM
components are shown in Fig. 6. Here we focus on the Einasto reconstruction
because of its greater flexibility; qualitatively similar results are obtained
in the NFW case. The fitted NFW concentrations as well as the total masses
obtained with this procedure are given in Table 2. In the left-hand panel of
Fig. 6 we show the output differential mass density profiles for the BCG
(green), gas mass (red), satellite galaxies (cyan) and DM (magenta). In the
innermost regions ($R<0.02R_{500}$) the stellar mass of the BCG dominates
because of the gravitational collapse of cold baryons at the bottom of the
potential well. The DM component takes over beyond $\sim 0.02R_{500}$ as the
dominant mass component. The NFW concentrations obtained for the fit to the DM
component only are 10-20% smaller than what is obtained when fitting the total
gravitational field (see Table 2), which results from the subtraction of the
sharply peaked BCG profile in cluster cores. The baryon-driven adiabatic
expansion of the DM halo under the influence of active galactic nucleus (AGN)
feedback results in slightly lower concentrations (Duffy et al., 2010), in
agreement with our findings. The DM profiles look remarkably similar in the
radial range [0.2-0.6]$R_{500}$, where our method is most accurate (see Sect.
6.1). The shape of the hot gas component is similar to that of the DM. It
dominates the baryonic content of clusters beyond $\sim 0.05R_{500}$ and even
slightly “catches up” with the DM in the outskirts. In the case of A2319 and
A644, the mass profile was found to be better described by a Burkert profile,
which leads to the curved shape of the DM profiles observed here.
Nongravitational heating processes such as AGN feedback inject energy into the
gaseous atmosphere, making its profile somewhat shallower than that of the DM.
Finally, satellite galaxies dominate the stellar mass content beyond $\sim
0.1R_{500}$, such that the total mass of satellites within $R_{200}$ exceeds
that of the BCG by more than an order of magnitude. However, their
contribution to the total gravitational field never exceeds a few per cent.
Intracluster light, which is not included in our analysis, may slightly change
the shape of the total stellar component, as its distribution is thought to be
centrally concentrated but somewhat shallower than the BCG (Montes & Trujillo,
2014); however, the contribution of ICL to the total baryon budget is expected
to be small ($<5\%$; e.g., Gonzalez et al., 2007).

###  5.4 Cumulative baryon fraction

Table 4:  Relative contribution of gas and stars as well as the total baryonic fraction at several overdensity radii for the five objects with available complete stellar mass information. All the values are given as percentages. Cluster | $f_{\rm gas,2500}$ | $f_{\star,2500}$ | $f_{\rm bar,2500}$ | $f_{\rm gas,500}$ | $f_{\star,500}$ | $f_{\rm bar,500}$ | $f_{\rm gas,200}$ | $f_{\star,200}$ | $f_{\rm bar,200}$
---|---|---|---|---|---|---|---|---|---
A644 | $7.7\pm 0.1$ | $1.2\pm 0.4$ | $8.9\pm 0.5$ | $12.8\pm 0.6$ | $1.3\pm 0.5$ | $14.1\pm 1.1$ | $16.9\pm 1.2$ | $1.5\pm 0.6$ | $18.4\pm 1.7$
A1795 | $9.6\pm 0.1$ | $1.7\pm 0.6$ | $11.3\pm 0.6$ | $13.4\pm 0.3$ | $1.3\pm 0.5$ | $14.7\pm 0.7$ | $15.8\pm 0.5$ | $1.3\pm 0.4$ | $17.1\pm 1.0$
A2029 | $11.1\pm 0.1$ | $2.8\pm 0.9$ | $13.8\pm 1.1$ | $14.4\pm 0.3$ | $2.1\pm 0.7$ | $16.5\pm 1.0$ | $17.1\pm 0.7$ | $2.0\pm 0.7$ | $19.1\pm 1.4$
A2142 | $12.0\pm 0.1$ | $2.3\pm 0.8$ | $14.3\pm 0.9$ | $14.1\pm 0.2$ | $1.7\pm 0.6$ | $15.8\pm 0.8$ | $14.7\pm 0.5$ | $1.5\pm 0.5$ | $16.2\pm 1.0$
A2319 | $9.2\pm 0.1$ | $2.0\pm 0.7$ | $11.3\pm 0.8$ | $17.0\pm 0.2$ | $1.5\pm 0.5$ | $18.5\pm 0.7$ | $24.9\pm 0.5$ | $1.5\pm 0.5$ | $26.4\pm 1.0$

The data presented here also allow us to study the total cumulative baryon
fraction in our systems as a function of radius. The right-hand panel of Fig.
6 shows the integrated baryon fraction $f_{\rm bar}(r)=M_{\rm bar}(<r)/M_{\rm
tot}(<r)$ from the core to the outskirts, as well as the relative contribution
of gas and stars. As already highlighted above, the stellar content of the BCG
dominates inside a few per cent of $R_{500}$ and decreases steeply with
radius. The integrated mass of the gas component dominates over that of stars
beyond $\sim 0.1R_{500}$, and inside $R_{500}$ the ICM accounts for about 90%
of the total baryonic content of clusters. In Table 4 we quote the measured
baryon fraction values for the five systems studied here at three different
overdensity radii ($R_{2500}$, $R_{500}$, and $R_{200}$).

We observe an overall depletion of baryons in the $0.2-0.5R_{500}$ range,
which is compensated by the somewhat shallower slope of the gas density
profile with respect to the DM. Inside $R_{500}$ our estimated baryon fraction
is equal to the cosmic value, and even slightly exceeds it beyond this point.
Such an excess relative to the cosmic value is unlikely to be a real effect,
as clusters are expected to be closed boxes for baryons within their virial
radius (e.g., White et al., 1993). An excess baryon fraction is more likely a
trace of deviations from the HSE assumption, which are expected to bias HSE
masses toward low values (e.g., Rasia et al., 2006; Nelson et al., 2014). This
effect is particularly striking in the case of A2319, for which HSE masses are
known to be biased low (see the detailed discussion in Ghirardini et al.,
2018). In Eckert et al. (2019) we related the cumulative hydrostatic gas
fraction to the level of hydrostatic bias by comparing the measured gas
fractions computed with the NFW model with the predictions of numerical
simulations, which expect that at large radii the baryon fraction should match
the cosmic value with a small depletion of $5-10\%$ (Planelles et al., 2013).
We found that deviations from HSE are low, at the level of $6-10\%$, with the
notable exception of A2319. However, the stellar content had to be assumed.
Our measurements of the total baryon fraction are consistent with the results
of Eckert et al. (2019) within $R_{500}$ but slightly higher within $R_{200}$,
which can be explained by the different model used here. Indeed, the study
presented in Eckert et al. (2019) was based on NFW fits to the total density
profile, whereas here we consider the baryonic and DM components separately
and model the DM profile with the more flexible Einasto functional form. At
$R_{200}$ our baryon fractions exceed the cosmic value by 10-20%, potentially
indicating a stronger impact of nonthermal pressure beyond $R_{500}$. This
analysis thus shows that our mass profiles are accurate at the $\lesssim 20\%$
level throughout the range considered here, which has little impact on our
main conclusions.

###  5.5 The galaxy cluster radial acceleration relation

Figure 7: Radial acceleration relation for X-COP galaxy clusters. The solid
colored curves show the result of the Einasto fit to the baryons+DM model,
whereas the individual data points indicate the accelerations computed using
the nonparametric log-normal mixture model (NP). For comparison, the solid
black curve and shaded area show the RAR obtained by McGaugh et al. (2016) for
rotationally supported galaxies and the associated 0.06 dex scatter. The
dashed line indicates the one-to-one relation.

The availability of precise measurements of the gas, stellar, and DM mass
components allows us to study the relation between the total, observed
gravitational acceleration $g_{\rm obs}$ and the acceleration that would be
expected in the absence of DM from the sum of the detected baryonic
components, $g_{\rm bar}$. Studying the rotation curves of rotationally
supported galaxies, McGaugh et al. (2016) showed that the relation between
$g_{\rm bar}$ and $g_{\rm obs}$ is nearly universal across the considered
sample and smoothly deviates from the one-to-one relation when moving toward
the external regions where the gravitational force is low. If this relation is
found to be universal across all gravitationally bound halos, it could become
the smoking gun for modified gravity theories such as MOND (Milgrom, 1983) and
its developments. The RAR was found to hold over a wide range in stellar mass
(Lelli et al., 2017) with a low scatter of less than 0.13 dex and possibly
even less than 0.06 dex (Li et al., 2018). Studies of the RAR in galaxy
clusters found that the normalization of the observed acceleration is higher
than what was found in rotationally supported galaxies (Chan & Del Popolo,
2020; Tian et al., 2020; Pradyumna et al., 2021), although the uncertainties
were too large to study in detail the shape of the RAR.

In Fig. 7 we show the RAR in X-COP galaxy clusters. The baryonic and observed
acceleration were computed locally from the reconstructed profiles,

$g_{\rm obs}(r)=\frac{GM_{\rm HSE}(<r)}{r^{2}},$ (17) $g_{\rm
bar}(r)=\frac{G(M_{\rm
gas}(<r)+M_{\rm\star,BCG}(<r)+M_{\rm\star,sat}(<r))}{r^{2}}.$ (18)

For the observed acceleration we consider both the nonparametric
reconstruction and the results of the Einasto fit to the DM component only. As
can be seen in Fig. 7, the two methods provide consistent results throughout
the entire range. The data are compared with the RAR relation determined by
McGaugh et al. (2016) in rotationally supported galaxies,

$g_{\rm obs}=\frac{g_{\rm bar}}{1-\exp(-\sqrt{g_{\rm
bar}/g_{\dagger}})}\equiv\mathcal{F}(g_{\rm bar}),$ (19)

with $g_{\dagger}=1.2\times 10^{-10}$ m s-2. At high acceleration
($g>10^{-10}$ m s-2), the baryonic acceleration induced by the stellar mass of
the BCG dominates, as expected from the breakdown of the density profiles
shown in Fig. 6. Beyond this point, the observed gravitational acceleration
starts to exceed the expectation of the RAR in spiral galaxies. In this
regime, the baryonic acceleration rapidly decreases, whereas the total
acceleration remains almost flat. This corresponds to the $[0.05-0.2]R_{500}$
radial range in Fig. 6 where the stellar mass of the BCG falls off rapidly and
the DM component largely dominates. The behavior of the relation is
inconsistent with the McGaugh et al. (2016) relation at high statistical
significance, even when considering the nonparametric reconstruction only. The
observed acceleration is up to 5 times higher than expected, which cannot be
explained by a potential bias in our mass reconstruction, since the HSE method
rather tends to underestimate the true mass (see Sect. 6.1 for a full
discussion).

Around $g_{\rm bar}\sim 2\times 10^{-11}$ m s-2 we observe another change of
regime, where both the observed and the baryonic components decrease with a
slope that is almost parallel to the one-to-one relation. This corresponds to
the regime where the baryonic component is dominated by the ICM mass and where
the gas density profiles are slightly shallower than the DM profiles. Our
systems eventually catch up with the McGaugh et al. (2016) RAR in the
outermost regions, although the slope of the relation is widely different. The
complex shape observed here indicates that the RAR in rotationally supported
galaxy is not a universal property of gravity, since it fails to reproduce a
gravitational field of similar strength in massive halos.

###  5.6 Modeling the galaxy cluster RAR

Figure 8:  Modeling of the RAR in galaxy clusters. _Left:_ Same as in Fig. 7.
The thick dashed orange curve and shaded area show the best fit to the NP data
with the model that describes two transitions (Eq. 20) and its corresponding
$1\sigma$ error envelope. The dotted line shows the classical MOND relation.
_Right:_ Breakdown of the various mass components in the RAR paradigm for the
five systems with complete descriptions of the baryonic components assuming
that the McGaugh et al. (2016) relation is universal. The curves show the
cumulative radial profiles of the gas mass (dashed green), stellar mass
(dotted cyan), and the required missing $M_{DM,RAR}$ computed using Eq. 21
(solid magenta).

While our data unequivocally confirm that galaxy clusters do not follow an
acceleration relation similar to that of disk galaxies, Fig. 7 shows that
galaxy clusters define their own relation in the $g_{\rm obs}-g_{\rm bar}$
plane, which we can attempt to quantify. The complex shape of the relation
implies that it cannot be described simply with a single acceleration scale.
To model the observed relation, we introduce a modification of the MOND law
with a second acceleration scale as

$g_{\rm obs}=g_{\rm bar}\left(1+\frac{g_{1}}{g_{\rm
bar}}\right)^{\alpha_{1}}\left(1+\frac{g_{2}}{g_{\rm
bar}}\right)^{\alpha_{2}-\alpha_{1}},$ (20)

with $g_{1},g_{2}$ the two characteristic acceleration scales describing the
transition between the three regimes, and $\alpha_{1},\alpha_{2}$ the slopes
of the correction term in the regime where the corresponding term dominates.
This functional forms models a correction that is proportional to $g_{\rm
bar}^{-\alpha_{2}}$ when $g_{\rm bar}\ll g_{2}$, transitioning to $g_{\rm
bar}^{-\alpha_{1}}$ when $g_{2}<g_{\rm bar}<g_{1}$ and converging to 1 when
$g_{\rm bar}\gg g_{1}$. Our empirical formula approaches the MOND formalism at
low acceleration when $\alpha_{1}=\alpha_{2}=0.5$ and $g_{1}$ is the classical
MOND acceleration $a_{0}\approx 10^{-10}$ m/s2.

We fitted the nonparametric points on the $g_{\rm obs}-g_{\rm bar}$ plane with
the model described in Eq. 20 and attempted to quantify the parameters of the
model and the intrinsic scatter. The model provides an adequate description of
the data at hand, with an intrinsic scatter $\sigma_{\ln g_{\rm obs}}=0.21\pm
0.02$. The best-fit curve is shown in the left-hand panel of Fig. 8 where the
functional form and its uncertainty are compared to the data. For the two
acceleration scales we obtain $g_{1}=2.3_{-0.8}^{+1.2}\times 10^{-10}$ m s-2
and $g_{2}=1.2_{-0.7}^{+0.9}\times 10^{-11}$ m s-2, in qualitative agreement
with the above discussion. The fitted slopes yield
$\alpha_{1}=1.4_{-0.5}^{+0.7}$ and $\alpha_{2}=-1.2_{-0.7}^{+0.9}$. The scale
$g_{1}$ is of the same order as the classical MOND acceleration, albeit
slightly higher. However, the fitted slope $\alpha_{1}$ is considerably larger
than its corresponding MOND counterpart, which indicates that the transition
away from the baryon-dominated regime is much faster than in disk galaxies.
Below the second characteristic scale $g_{2}$ the $g_{\rm obs}-g_{\rm bar}$
relation is much steeper than expected in MOND.

As previously recognized by several authors (e.g., Sanders, 1999, 2003), an
extra mass component is required to reconcile the gravitational field of
galaxy clusters with the MOND paradigm. In E19 we showed that the total
baryonic mass in X-COP clusters within $R_{500}$ should exceed the measured
baryonic mass by about a factor of 2. Our data can be used to determine the
mass profile of the missing component assuming that the McGaugh et al. (2016)
relation (Eq. 19) applies. Specifically, for each value of $g_{\rm obs}$ we
can determine the missing enclosed mass $M_{DM,RAR}(<r)$ such that the
expected acceleration matches the observed one. In other terms, for any given
value of $g_{\rm obs}$ and $g_{\rm bar}$ we numerically search for the value
of $M_{DM,RAR}$ that satisfies the condition

$g_{\rm obs}=\mathcal{F}\left(g_{\rm
bar}+\frac{GM_{DM,RAR}(<r)}{r^{2}}\right).$ (21)

We used the Einasto reconstructions to determine the mass profiles of
$M_{DM,RAR}$. In the right-hand panel of Fig. 8 we show the enclosed mass
profiles of the missing component assuming that Eq. 19 holds. The missing
component dominates the mass profile beyond $\sim 0.05R_{500}$, similar to the
CDM case (Fig. 6). The required enclosed mass exceeds the gas mass and reaches
a maximum around $\sim 0.5R_{500}$ ($\sim 600-700$ kpc). This corresponds to
the second transition scale $g_{2}$, beyond which the slope of the relation is
steeper than the expectation. As described in Sect. 5.5, the galaxy cluster
data eventually catch up with the McGaugh et al. (2016) relation in the
outskirts, implying that the required correction is negligible and the
enclosed mass profile decreases. Since the cumulative mass cannot decrease,
this analysis shows that the required mass distribution of the missing mass in
MOND is unphysical, unless our analysis is hampered by substantial systematic
uncertainties (see Sect. 6.1).

##  6 Discussion

###  6.1 Systematic uncertainties

A key issue to be addressed is the potential impact of systematic
uncertainties. While our tests with mock data demonstrate that our analysis
pipeline introduces minimal biases, from spectral fitting to hydrostatic mass
reconstruction (see Sect. 4), a number of other sources of systematic
uncertainty need to be addressed. Here we discuss the main sources of
uncertainty: hydrostatic bias, gas inhomogeneities, and effective area
calibration.

Hydrostatic bias: In the presence of mergers and nongravitational energy
input, such as AGN feedback, residual gas motions can act as an additional
source of pressure support on top of thermal pressure (e.g., Rasia et al.,
2004; Lau et al., 2009; Nelson et al., 2014; Biffi et al., 2016). In the
presence of a substantial level of residual gas motions, mass estimates
derived under the hydrostatic assumption are known to be biased toward low
values, an effect that various numerical simulations predict to be in the
range $5-30\%$ at $R_{500}$ (Rasia et al., 2004; Nagai et al., 2007;
Angelinelli et al., 2020; Barnes et al., 2021; Bennett & Sijacki, 2021).
Through a direct comparison between hydrostatic masses and masses obtained
using alternative methods (Ettori et al., 2019), we found that X-COP
hydrostatic masses are 10-15% lower than weak lensing estimates (Herbonnet et
al., 2020). In Eckert et al. (2019) we used the integrated gas fraction as an
anchor assuming that the true gas fraction can be robustly predicted by the
$\Lambda$CDM framework to estimate the level of hydrostatic bias, which was
found to be low (7% at $R_{500}$, 10% at $R_{200}$). We find similar results
in this analysis, as described in Sect. 5.4, where we can see that, with the
exception of a single system that is known to be strongly biased (A2319;
Ghirardini et al., 2018), the measured baryon fractions agree with the
expectations at $R_{500}$ and slightly exceed the cosmic baryon fraction
beyond this point, indicating a mild level of hydrostatic bias. Since the
relative level of nonthermal pressure support is expected to increase with
distance to the cluster core (e.g., Nelson et al., 2014), we conclude that our
measurements are mildly biased ($<10\%$) out to $R_{500}$. The HSE assumption
is likely to be a poorer approximation beyond this point.

Gas inhomogeneities: Inhomogeneities in the gas distribution can potentially
impact the recovered thermodynamic profiles, both for the gas density
(Mathiesen et al., 1999; Nagai & Lau, 2011; Vazza et al., 2013) and the
temperature profiles (Rasia et al., 2014). In the presence of overdense, cool
substructures that have not yet mixed with the surrounding plasma, the
recovered X-ray densities and temperatures are expected to be biased toward
high and low values, respectively, compared to the mass-weighted quantities.
The multi-temperature structure may also be important in the innermost regions
($<10$ kpc) where the hot ICM may mix with the cooler gas content of the BCG.
Several works based on numerical simulations (e.g., Rasia et al., 2012; Pearce
et al., 2020) predict that temperature inhomogeneities would bias the X-ray
derived temperatures, which might greatly increase the bias in the recovered
HSE masses. For instance, Pearce et al. (2020) predict that masses based on
X-ray data only may be largely biased low (up to 50% in the X-COP mass range),
whereas HSE masses using SZ-derived pressure profiles should be closer to the
true, mass-weighted hydrostatic masses. However, we note that our gas density
profiles are estimated using the azimuthal median technique (Eckert et al.,
2015), which excises the regions of enhanced surface brightness and returns
gas density profiles that are free from the clumping effect. The regions
corresponding with obvious substructures are masked during the spectral
extraction procedure, such that their impact on the recovered temperatures
should be limited. Our masses estimated from X-ray data only are consistent
with the results obtained in combination with SZ data (Ettori et al., 2019)
and the pressure profiles determined independently from X-ray and SZ data are
in excellent agreement (Ghirardini et al., 2019). Therefore, we conclude that
the impact of gas inhomogeneities is small ($<10\%$) and certainly not as
severe as predicted by, for example, Pearce et al. (2020). The difference
could be explained either by the different procedures used in the simulations
and observations or by an incomplete mixing in the corresponding simulations.
To investigate the origin of the difference, it would be necessary to apply
our pipeline on mock observations of simulated clusters, similar to the
analysis presented in Sect. 4, which is beyond the scope of this paper.

Effective area calibration: The calibration of the _XMM-Newton_ effective area
affects the measured spectroscopic temperatures, as the spectral model needs
to be folded through the instrumental response to fit the X-ray spectra. There
is a known inconsistency between the temperatures determined by _XMM-Newton_
/EPIC and by _Chandra_ /ACIS (Nevalainen et al., 2010; Schellenberger et al.,
2015), with ACIS returning systematically higher temperatures than EPIC when
fitting over the entire 0.5-10 keV energy range. The discrepancy increases
with plasma temperature, and reaches $\sim 20\%$ for 10 keV plasma. It is
still unclear whether the _XMM-Newton_ temperatures derived here are
accurately reproducing the true spectroscopic temperatures. If the _Chandra_
temperatures are correct, our masses at fixed radii would be underestimated by
a similar amount, given that to first order the HSE mass is proportional to
the fitted temperature. Our reconstruction of the gravitational field is based
on a combination of _XMM-Newton_ and _Planck_ data, the latter providing the
dominant contribution to the temperature profiles at large radii ($\geq
R_{500}$). Therefore, the uncertainty in the effective area calibration mostly
affects the regions located inside $0.5R_{500}$, where our temperatures are
constrained primarily by the spectroscopic X-ray measurements.

Overall, we estimate that the systematic uncertainties associated with each of
the effects discussed here amount to $10-20\%$, which should not affect the
conclusions of this paper with the exception of the results on the baryon
fraction (Table 4).

###  6.2 Consistency with the $\Lambda$CDM framework

All of the results presented in this paper agree with the predictions of the
$\Lambda$CDM framework and the bottom-up structure formation paradigm. The
average mass profile of the X-COP clusters is well represented by an NFW
model, with deviations of at most 10% over the entire radial range of interest
($[0.01-2]R_{500}$; see the left-hand panel of Fig. 4). The average NFW
concentration of the sample matches perfectly the $\Lambda$CDM predictions
(right-hand panel of Fig. 4), and N-body simulations also reproduce accurately
the observed scatter in the mass-concentration relation at $M_{200}\sim
10^{15}M_{\odot}$. Finally, with the exception of a system known for its high
level of nonthermal pressure support (A2319), the total cumulative baryon
budget (gas + stars) is close to the cosmic baryon fraction
$\Omega_{b}/\Omega_{m}$, with a slight excess of 10-20% that can most likely
be attributed to a hydrostatic bias (see Sect. 5.4). Numerical simulations
including different gas physics and hydrodynamic solvers generically predict
that in the most massive halos the gravitational binding energy largely
exceeds nongravitational energy input, for example from AGN and supernovae. As
a result, the baryon fraction enclosed within the virial radius should closely
match the universal value (White et al., 1993; Eke et al., 1998; Ettori et
al., 2009; Planelles et al., 2013; Eckert et al., 2019; Mantz et al., 2022),
in agreement with the findings presented here. We note however that given the
uncertain impact of hydrostatic bias, our baryon fractions only provide a weak
test of $\Lambda$CDM, and more accurate mass reconstructions combined, for
example, with weak lensing data are required to determine the exact baryon
fraction of galaxy clusters within $R_{200}$.

While the NFW model provides an adequate representation of the average mass
profiles, we do observe deviations from the predicted universal shape, both in
the innermost ($R<0.1R_{500}$) and outermost ($R>R_{500}$) regions covered in
our study. The profiles estimated through our nonparametric technique appear
to show a wider variety of shapes than can be described by the NFW parametric
form (see Fig. 3). Conversely, the additional degree of freedom afforded by
the Einasto parametric form allows the observed behavior to be reproduced more
closely in each individual case. Some systems (e.g., A2255 and A644) show a
larger amount of curvature in their mass profiles; the same systems were also
found in E19 to be better described by models including a central core
(Burkert, isothermal sphere). These systems exhibit a somewhat disturbed X-ray
morphology as indicated by their large centroid shifts (see Paper II); thus,
it is possible that the X-ray peak does not trace exactly the bottom of the
potential well. This effect is likely to be present in A2255, which exhibits a
very flat central gas density profile and a substantial offset between the
X-ray peak and the position of the BCG (123 kpc, Rossetti et al., 2016). The
apparent flat mass distribution in the core of this system may thus be
explained by miscentering. On top of that, while the NFW profile is expected
to provide an excellent description of the average mass profile extracted from
CDM simulations, individual profiles show a scatter of $\sim 0.1$ dex around
the mean (see, e.g., Bhattacharya et al., 2013; Ludlow et al., 2016), in
agreement with our observations (see Fig. 4). In conclusion, the study
presented here does not show any significant tension with $\Lambda$CDM
predictions, and extensions of the paradigm, either in the form of a different
DM candidate (e.g., WDM or SIDM) or modifications of the theory of gravity,
are not required by our data.

###  6.3 The RAR is not universal

Through our decomposition of the measured gravitational field into DM and
baryonic components (gas, BCG, and satellite galaxies), we provided a detailed
reconstruction of the RAR in a subset of five systems (see Fig. 7). Our study
improves over previous attempts to test the RAR in the galaxy cluster regime
(Chan & Del Popolo, 2020; Tian et al., 2020; Pradyumna et al., 2021) as it
combines our measurements of the gravitational field over two decades in
radius with a full characterization of the various baryonic components (see
Fig. 6). The relation derived by McGaugh et al. (2016) for disk galaxies fails
to reproduce the relation between observed and baryonic accelerations in X-COP
clusters at very high significance. As described in Sect. 5.5, we observe
three distinct regimes in the $g_{\rm obs}-g_{\rm bar}$ plane: at high
acceleration ($g_{\rm bar}>10^{-10}$ m s-2), the gravitational field is
dominated by the stellar content of the BCG. Below this threshold, the
relation departs from the McGaugh et al. (2016) relation and the observed
acceleration becomes up to 5 times higher than that expected from the RAR in
spiral galaxies. Below a second characteristic scale, $g_{\rm bar}\sim 2\times
10^{-11}$ m s-2, the relation steepens again and eventually catches up with
the McGaugh et al. (2016) relation. The discrepancy cannot be attributed to
potential systematics in our reconstruction such as hydrostatic bias (see
Sect. 6.1) as they would most likely lead to an underestimation of the
gravitational acceleration and increase the gap. Therefore, the behavior of
the relation in X-COP clusters vastly differs from that in spiral galaxies,
both qualitatively and quantitatively, indicating that the RAR is apparently
not universal.

If the tight relation observed in spiral galaxies results from a modification
of the theory of gravity at large scales (e.g., MOND), the observed behavior
shows that an additional source of acceleration is required on top of the
known baryonic components (see also Ettori et al., 2019). A possible way of
alleviating the tensions with the MOND paradigm in the cluster regime is to
invoke an additional matter component such as a light sterile neutrino with a
mass in the range $1-10$ eV (e.g., Sanders, 2007; Angus et al., 2008;
Nieuwenhuizen & Morandi, 2013), which may be a natural candidate for a minimal
extension of the standard model of particle physics. The missing component
would dominate the matter content at intermediate radii where the discrepancy
with MOND is large. Angus et al. (2010) argue that in the case of a sterile
neutrino with a mass of 11 eV the gravitational field of galaxy clusters would
be deep enough to retain the sterile neutrinos, while at galaxy scales the
neutrinos would free stream from the halo. However, we note that invoking a
substantial amount of hot DM affects the formation of structures in the
Universe, which renders the expected halo mass function inconsistent with the
observations (Angus et al., 2013). This scenario also invokes at the same time
a modification of the gravitational law and a (less abundant) DM component,
making it intrinsically more complex. Our analysis of the required missing
mass in the standard RAR calibrated on disk galaxies indicates that this
component should dominate the mass budget in the inner regions and its total
mass should be comparable to the total baryonic mass (see Fig. 8). However,
the mass distribution of the missing component should exhibit an unphysical
decreasing trend beyond $\sim 0.5R_{500}$, since our data at large radii
eventually catch up with the McGaugh et al. (2016) relation. If the mass
density of the MOND DM is set to exactly zero beyond the radius where the
missing component peaks, our observed gravitational acceleration would need to
be biased low by about a factor of 2 at $R_{200}$ to match the predicted
acceleration. While not totally impossible, this discrepancy is quite a bit
larger than our assessment of systematic uncertainties (see Sect. 6.1) and
thus this scenario is disfavored by our data.

Alternatively, it has been claimed that the MOND characteristic scale $a_{0}$
may be mass dependent (Hodson & Zhao, 2017). While our clusters indeed define
their own RAR with a scatter $\sigma_{\ln g}=0.21$, the behavior of the
recovered RAR in X-COP clusters is quite complex and our model requires two
characteristic acceleration scales, $g_{1}$ and $g_{2}$ (see Eq. 20). The
fitted value of $g_{1}$ is close to $a_{0}$, whereas the second characteristic
scale $g_{2}$ is about an order of magnitude smaller. The $g_{\rm bar}-g_{\rm
obs}$ relation steepens again below $g_{2}$, which cannot be easily explained
by a mass dependence of $a_{0}$.

Overall, our results can be more easily explained in the $\Lambda$CDM
scenario. The decomposition of the gravitational field into its various
components can be well explained: the prevalence of the stellar content in the
innermost regions reflects the collapse of cold baryons to the bottom of the
potential well, whereas the somewhat flatter distribution of hot gas with
respect to the DM profile can be explained by AGN feedback, which is known to
inject nongravitational energy into the surrounding medium (e.g., Le Brun et
al., 2014). Moreover, a few galaxy evolution models are also able to reproduce
the shape of the RAR in rotationally supported galaxies in the $\Lambda$CDM
context (Navarro et al., 2017; Ludlow et al., 2017; Dutton et al., 2019).
Therefore, it is conceivable that the RAR may arise as a byproduct of the
galaxy formation process as an interplay between baryonic and total
acceleration.

##  7 Conclusions

In this paper we have performed a detailed analysis of the gravitational field
of X-COP galaxy clusters over the radial range $[0.01-2]R_{500}$ from a
combination of deep _XMM-Newton_ and _Planck_ data. Our results can be
summarized as follows:

  * •

We introduced a novel framework for the reconstruction of hydrostatic mass
profiles from X-ray and/or SZ data, which we distribute in the form of the
public Python package hydromass. Hydrostatic equilibrium profiles can be
reconstructed assuming a mass model, from a parametric pressure profile, or in
a nonparametric way as a linear combination of log-normal functions. Our
framework improves over our previous works (Ettori et al. 2019; the
differences in the mass values are less than 10 per cent; see Table 3) as it
implements PSF deconvolution, affecting in particular the inner regions, and
fits with the Einasto mass model, where the parameter $\alpha$ is left free to
vary. The various methods discussed here are integrated within a common
framework that includes an efficient Bayesian optimization scheme.

  * •

We validated our method extensively using a set of mock _XMM-Newton_
observations of a synthetic NFW cluster in HSE. Our mock observations include
a wide range of instrumental effects, such as effective area, energy
redistribution, vignetting, PSF convolution, and a sophisticated background
model. The mock observations were analyzed using the same pipeline as the
actual observations. Our code was found to reproduce the input profile very
accurately, both for the 3D temperature profile and the hydrostatic mass (see
Fig. 2).

  * •

Applying our method to the 12 X-COP galaxy clusters, we find that the NFW and
Einasto profiles both provide a very good representation of the average mass
profiles over the entire radial range, in agreement with the predictions of
the $\Lambda$CDM paradigm. The additional flexibility afforded by the Einasto
model allows it to trace the individual profiles more closely than the NFW in
the innermost and outermost regions.

  * •

The very high statistical quality of our data allows us to measure the NFW
concentration of individual systems with typical uncertainties of a few per
cent (see Fig. 4). Modeling the relation between concentration and mass in our
sample, we find an average concentration $c_{200}=3.69\pm 0.39$ at
$M_{200}\approx 10^{15}M_{\odot}$, in remarkable agreement with $\Lambda$CDM
predictions. The intrinsic scatter of the mass-concentration relation was
found to be $\sigma_{\ln c_{200}}=0.37_{-0.07}^{+0.11}$, which again is
consistent with the expectations of N-body simulations of structure formation.

  * •

For a subset of five systems, we decomposed the gravitational field into its
baryonic and DM components. We found that the stellar content of the BCG
dominates the gravitational field in the innermost regions ($<0.02R_{500}$),
which can be explained by the collapse of cold baryons to the bottom of the
potential well. The contribution of the stellar mass decreases sharply with
radius, such that we observe a depletion of baryons in the range
$[0.1-0.5]R_{500}$. The hot ICM dominates the baryonic mass beyond $\sim
0.1R_{500}$ and exhibits a somewhat shallower profile than the DM. Our systems
reach the cosmic baryon fraction at the virial radius. We note a slight excess
at large radii with respect to the cosmic baryon fraction, which we interpret
as evidence for a mild contribution of nonthermal pressure ($10-20\%$) at
$R_{200}$.

  * •

We studied in detail the relation between observed and baryonic acceleration
(RAR) in the five systems for which the baryonic content can be completely
characterized (Fig. 7). The recovered relation deviates from the relation
observed in spiral galaxies (McGaugh et al., 2016) at high significance, with
the observed acceleration exceeding the expected value by a factor of $\sim 5$
around $g_{\rm bar}=2\times 10^{-11}$ m s-2. We found that the RAR in galaxy
clusters exhibits a complex shape with three distinct regimes delimited by two
characteristic scales, which makes it difficult to reconcile with the MOND
paradigm.

Application of the presented framework to larger samples (e.g., CHEX-MATE, The
CHEX-MATE Collaboration et al., 2020) will allow us to study in detail the
shape of galaxy cluster mass profiles and test the validity of the HSE
assumption.

###### Acknowledgements.

SE acknowledges financial contribution from the contracts ASI-INAF Athena
2019-27-HH.0, “Attività di Studio per la comunità scientifica di Astrofisica
delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASI-INAF n.
2017-14-H.0), INAF mainstream project 1.05.01.86.10, and from the European
Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n.
871158). SIL is supported in part by the National Research Foundation of South
Africa (NRF Grant Number: 120850). Opinions, findings and conclusions or
recommendations expressed in this publication is that of the author(s), and
the NRF accepts no liability whatsoever in this regard.

## References

  * Anders & Grevesse (1989) Anders, E. & Grevesse, N. 1989, Geochim. Cosmochim. Acta., 53, 197
  * Angelinelli et al. (2020) Angelinelli, M., Vazza, F., Giocoli, C., et al. 2020, MNRAS, 495, 864
  * Angus et al. (2013) Angus, G. W., Diaferio, A., Famaey, B., & van der Heyden, K. J. 2013, MNRAS, 436, 202
  * Angus et al. (2008) Angus, G. W., Famaey, B., & Buote, D. A. 2008, MNRAS, 387, 1470
  * Angus et al. (2010) Angus, G. W., Famaey, B., & Diaferio, A. 2010, MNRAS, 402, 395
  * Arnaud et al. (2010) Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92
  * Barnes et al. (2021) Barnes, D. J., Vogelsberger, M., Pearce, F. A., et al. 2021, MNRAS, 506, 2533
  * Bartalucci et al. (2019) Bartalucci, I., Arnaud, M., Pratt, G. W., Démoclès, J., & Lovisari, L. 2019, A&A, 628, A86
  * Bekenstein (2004) Bekenstein, J. D. 2004, Phys. Rev. D, 70, 083509
  * Bennett & Sijacki (2021) Bennett, J. S. & Sijacki, D. 2021, arXiv e-prints, arXiv:2110.07326
  * Bhattacharya et al. (2013) Bhattacharya, S., Habib, S., Heitmann, K., & Vikhlinin, A. 2013, ApJ, 766, 32
  * Biffi et al. (2016) Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112
  * Biviano et al. (2013) Biviano, A., Rosati, P., Balestra, I., et al. 2013, A&A, 558, A1
  * Bourdin et al. (2017) Bourdin, H., Mazzotta, P., Kozmanyan, A., Jones, C., & Vikhlinin, A. 2017, ApJ, 843, 72
  * Brown et al. (2020) Brown, S. T., McCarthy, I. G., Diemer, B., et al. 2020, MNRAS, 495, 4994
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Bullock et al. (2001) Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559
  * Burgess et al. (2021) Burgess, J. M., Fleischhack, H., Vianello, G., et al. 2021, The Multi-Mission Maximum Likelihood framework (3ML), Zenodo
  * Cappellari (2002) Cappellari, M. 2002, MNRAS, 333, 400
  * Cappellari (2008) Cappellari, M. 2008, MNRAS, 390, 71
  * Cappellari & Copin (2003) Cappellari, M. & Copin, Y. 2003, MNRAS, 342, 345
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Chan & Del Popolo (2020) Chan, M. H. & Del Popolo, A. 2020, MNRAS, 492, 5865
  * Chiu et al. (2018) Chiu, I., Mohr, J. J., McDonald, M., et al. 2018, MNRAS, 478, 3072
  * Croston et al. (2006) Croston, J. H., Arnaud, M., Pointecouteau, E., & Pratt, G. W. 2006, A&A, 459, 1007
  * Diemer (2018) Diemer, B. 2018, ApJS, 239, 35
  * Diemer & Joyce (2019) Diemer, B. & Joyce, M. 2019, ApJ, 871, 168
  * Diemer & Kravtsov (2014) Diemer, B. & Kravtsov, A. V. 2014, ApJ, 789, 1
  * Diemer & Kravtsov (2015) Diemer, B. & Kravtsov, A. V. 2015, ApJ, 799, 108
  * Duffy et al. (2008) Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008, MNRAS, 390, L64
  * Duffy et al. (2010) Duffy, A. R., Schaye, J., Kay, S. T., et al. 2010, MNRAS, 405, 2161
  * Dutton & Macciò (2014) Dutton, A. A. & Macciò, A. V. 2014, MNRAS, 441, 3359
  * Dutton et al. (2019) Dutton, A. A., Macciò, A. V., Obreja, A., & Buck, T. 2019, MNRAS, 485, 1886
  * Eckert et al. (2017) Eckert, D., Ettori, S., Pointecouteau, E., et al. 2017, Astronomische Nachrichten, 338, 293
  * Eckert et al. (2020) Eckert, D., Finoguenov, A., Ghirardini, V., et al. 2020, The Open Journal of Astrophysics, 3, 12
  * Eckert et al. (2019) Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40
  * Eckert et al. (2014) Eckert, D., Molendi, S., Owers, M., et al. 2014, A&A, 570, A119
  * Eckert et al. (2015) Eckert, D., Roncarelli, M., Ettori, S., et al. 2015, MNRAS, 447, 2198
  * Einasto (1965) Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87
  * Eke et al. (1998) Eke, V. R., Navarro, J. F., & Frenk, C. S. 1998, ApJ, 503, 569
  * Ettori et al. (2010) Ettori, S., Gastaldello, F., Leccardi, A., et al. 2010, A&A, 524, A68
  * Ettori et al. (2017) Ettori, S., Ghirardini, V., Eckert, D., Dubath, F., & Pointecouteau, E. 2017, MNRAS, 470, L29
  * Ettori et al. (2019) Ettori, S., Ghirardini, V., Eckert, D., et al. 2019, A&A, 621, A39
  * Ettori et al. (2009) Ettori, S., Morandi, A., Tozzi, P., et al. 2009, A&A, 501, 61
  * Evrard (1997) Evrard, A. E. 1997, MNRAS, 292, 289
  * Famaey & McGaugh (2012) Famaey, B. & McGaugh, S. S. 2012, Living Reviews in Relativity, 15, 10
  * Fasano et al. (2010) Fasano, G., Bettoni, D., Ascaso, B., et al. 2010, MNRAS, 404, 1490
  * Geller et al. (2014) Geller, M. J., Hwang, H. S., Diaferio, A., et al. 2014, ApJ, 783, 52
  * Ghirardini et al. (2019) Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41
  * Ghirardini et al. (2018) Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7
  * Ghizzardi et al. (2021) Ghizzardi, S., Molendi, S., van der Burg, R., et al. 2021, A&A, 646, A92
  * Gonzalez et al. (2007) Gonzalez, A. H., Zaritsky, D., & Zabludoff, A. I. 2007, ApJ, 666, 147
  * Haridasu et al. (2021) Haridasu, B. S., Karmakar, P., De Petris, M., Cardone, V. F., & Maoli, R. 2021, arXiv e-prints, arXiv:2111.01101
  * Harikumar & Biesiada (2022) Harikumar, S. & Biesiada, M. 2022, European Physical Journal C, 82, 241
  * Herbonnet et al. (2020) Herbonnet, R., Sifón, C., Hoekstra, H., et al. 2020, MNRAS, 497, 4684
  * Hodson & Zhao (2017) Hodson, A. O. & Zhao, H. 2017, A&A, 598, A127
  * Hoffman & Gelman (2014) Hoffman, M. D. & Gelman, A. 2014, Journal of Machine Learning Research, 15, 1593
  * Hurier et al. (2013) Hurier, G., Macías-Pérez, J. F., & Hildebrandt, S. 2013, A&A, 558, A118
  * Ishiyama et al. (2021) Ishiyama, T., Prada, F., Klypin, A. A., et al. 2021, MNRAS, 506, 4210
  * Kalberla et al. (2005) Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775
  * King (1962) King, I. 1962, AJ, 67, 471
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340
  * Kravtsov et al. (2005) Kravtsov, A. V., Nagai, D., & Vikhlinin, A. A. 2005, ApJ, 625, 588
  * Laganá et al. (2013) Laganá, T. F., Martinet, N., Durret, F., et al. 2013, A&A, 555, A66
  * Lau et al. (2009) Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129
  * Le Brun et al. (2014) Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2014, MNRAS, 441, 1270
  * Leccardi & Molendi (2008) Leccardi, A. & Molendi, S. 2008, A&A, 486, 359
  * Lelli et al. (2017) Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017, ApJ, 836, 152
  * Li et al. (2018) Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2018, A&A, 615, A3
  * Loubser et al. (2020) Loubser, S. I., Babul, A., Hoekstra, H., et al. 2020, MNRAS, 496, 1857
  * Loubser et al. (2018) Loubser, S. I., Hoekstra, H., Babul, A., & O’Sullivan, E. 2018, MNRAS, 477, 335
  * Ludlow & Angulo (2017) Ludlow, A. D. & Angulo, R. E. 2017, MNRAS, 465, L84
  * Ludlow et al. (2017) Ludlow, A. D., Benítez-Llambay, A., Schaller, M., et al. 2017, Phys. Rev. Lett., 118, 161103
  * Ludlow et al. (2016) Ludlow, A. D., Bose, S., Angulo, R. E., et al. 2016, MNRAS, 460, 1214
  * Macciò et al. (2008) Macciò, A. V., Dutton, A. A., & van den Bosch, F. C. 2008, MNRAS, 391, 1940
  * Macciò et al. (2013) Macciò, A. V., Ruchayskiy, O., Boyarsky, A., & Muñoz-Cuartas, J. C. 2013, MNRAS, 428, 882
  * Mamon & Łokas (2005) Mamon, G. A. & Łokas, E. L. 2005, MNRAS, 362, 95
  * Mantz et al. (2014) Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2014, MNRAS, 440, 2077
  * Mantz et al. (2022) Mantz, A. B., Morris, R. G., Allen, S. W., et al. 2022, MNRAS, 510, 131
  * Mathiesen et al. (1999) Mathiesen, B., Evrard, A. E., & Mohr, J. J. 1999, ApJ, 520, L21
  * Mazzotta et al. (2004) Mazzotta, P., Rasia, E., Moscardini, L., & Tormen, G. 2004, MNRAS, 354, 10
  * McCammon et al. (2002) McCammon, D., Almy, R., Apodaca, E., et al. 2002, ApJ, 576, 188
  * McGaugh et al. (2016) McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Phys. Rev. Lett., 117, 201101
  * Milgrom (1983) Milgrom, M. 1983, ApJ, 270, 384
  * Milgrom (2009) Milgrom, M. 2009, Phys. Rev. D, 80, 123536
  * Mitchell et al. (2018) Mitchell, M. A., He, J.-h., Arnold, C., & Li, B. 2018, MNRAS, 477, 1133
  * Montes & Trujillo (2014) Montes, M. & Trujillo, I. 2014, ApJ, 794, 137
  * Moster et al. (2011) Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113
  * Munari et al. (2014) Munari, E., Biviano, A., & Mamon, G. A. 2014, A&A, 566, A68
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJS, 206, 8
  * Nagai & Lau (2011) Nagai, D. & Lau, E. T. 2011, ApJ, 731, L10+
  * Nagai et al. (2007) Nagai, D., Vikhlinin, A., & Kravtsov, A. V. 2007, ApJ, 655, 98
  * Navarro et al. (2017) Navarro, J. F., Benítez-Llambay, A., Fattahi, A., et al. 2017, MNRAS, 471, 1841
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Navarro et al. (2004) Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039
  * Nelson et al. (2014) Nelson, K., Lau, E. T., & Nagai, D. 2014, ApJ, 792, 25
  * Neto et al. (2007) Neto, A. F., Gao, L., Bett, P., et al. 2007, MNRAS, 381, 1450
  * Nevalainen et al. (2010) Nevalainen, J., David, L., & Guainazzi, M. 2010, A&A, 523, A22
  * Nieuwenhuizen & Morandi (2013) Nieuwenhuizen, T. M. & Morandi, A. 2013, MNRAS, 434, 2679
  * Okabe et al. (2013) Okabe, N., Smith, G. P., Umetsu, K., Takada, M., & Futamase, T. 2013, ApJ, 769, L35
  * Pearce et al. (2020) Pearce, F. A., Kay, S. T., Barnes, D. J., Bower, R. G., & Schaller, M. 2020, MNRAS, 491, 1622
  * Planck Collaboration (2016) Planck Collaboration. 2016, A&A, 594, A22
  * Planck Collaboration V (2013) Planck Collaboration V. 2013, A&A, 550, A131
  * Planck Collaboration XIII (2016) Planck Collaboration XIII. 2016, A&A, 594, A13
  * Planck Collaboration XXIX (2014) Planck Collaboration XXIX. 2014, A&A, 571, A29
  * Planelles et al. (2013) Planelles, S., Borgani, S., Dolag, K., et al. 2013, MNRAS, 431, 1487
  * Pointecouteau et al. (2005) Pointecouteau, E., Arnaud, M., & Pratt, G. W. 2005, A&A, 435, 1
  * Pointecouteau et al. (2021) Pointecouteau, E., Santiago-Bautista, I., Douspis, M., et al. 2021, A&A, 651, A73
  * Pointecouteau & Silk (2005) Pointecouteau, E. & Silk, J. 2005, MNRAS, 364, 654
  * Pradyumna et al. (2021) Pradyumna, S., Gupta, S., Seeram, S., & Desai, S. 2021, Physics of the Dark Universe, 31, 100765
  * Pratt et al. (2019) Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25
  * Ragagnin et al. (2021) Ragagnin, A., Saro, A., Singh, P., & Dolag, K. 2021, MNRAS, 500, 5056
  * Rasia et al. (2006) Rasia, E., Ettori, S., Moscardini, L., et al. 2006, MNRAS, 369, 2013
  * Rasia et al. (2014) Rasia, E., Lau, E. T., Borgani, S., et al. 2014, ApJ, 791, 96
  * Rasia et al. (2012) Rasia, E., Meneghetti, M., Martino, R., et al. 2012, New Journal of Physics, 14, 055018
  * Rasia et al. (2004) Rasia, E., Tormen, G., & Moscardini, L. 2004, MNRAS, 351, 237
  * Robertson et al. (2021) Robertson, A., Massey, R., Eke, V., Schaye, J., & Theuns, T. 2021, MNRAS, 501, 4610
  * Rossetti et al. (2016) Rossetti, M., Gastaldello, F., Ferioli, G., et al. 2016, MNRAS, 457, 4515
  * Salucci & Burkert (2000) Salucci, P. & Burkert, A. 2000, ApJ, 537, L9
  * Salvatier et al. (2016) Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Computer Science, 2:e55
  * Sand et al. (2012) Sand, D. J., Graham, M. L., Bildfell, C., et al. 2012, ApJ, 746, 163
  * Sanders (1999) Sanders, R. H. 1999, ApJ, 512, L23
  * Sanders (2003) Sanders, R. H. 2003, MNRAS, 342, 901
  * Sanders (2007) Sanders, R. H. 2007, MNRAS, 380, 331
  * Sayers et al. (2013) Sayers, J., Czakon, N. G., Mantz, A., et al. 2013, ApJ, 768, 177
  * Schellenberger et al. (2015) Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30
  * Skordis & Złośnik (2021) Skordis, C. & Złośnik, T. 2021, Phys. Rev. Lett., 127, 161302
  * Smith et al. (2017) Smith, R. J., Lucey, J. R., & Edge, A. C. 2017, MNRAS, 471, 383
  * Smith et al. (2001) Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91
  * Snowden et al. (2008) Snowden, S. L., Mushotzky, R. F., Kuntz, K. D., & Davis, D. S. 2008, A&A, 478, 615
  * Tam et al. (2020) Tam, S.-I., Jauzac, M., Massey, R., et al. 2020, MNRAS, 496, 4032
  * Tchernin et al. (2016) Tchernin, C., Eckert, D., Ettori, S., et al. 2016, A&A, 595, A42
  * The CHEX-MATE Collaboration et al. (2020) The CHEX-MATE Collaboration, :, Arnaud, M., et al. 2020, arXiv e-prints, arXiv:2010.11972
  * Tian et al. (2020) Tian, Y., Umetsu, K., Ko, C.-M., Donahue, M., & Chiu, I. N. 2020, ApJ, 896, 70
  * Umetsu & Diemer (2017) Umetsu, K. & Diemer, B. 2017, ApJ, 836, 231
  * Umetsu et al. (2016) Umetsu, K., Zitrin, A., Gruen, D., et al. 2016, ApJ, 821, 116
  * van der Burg et al. (2015) van der Burg, R. F. J., Hoekstra, H., Muzzin, A., et al. 2015, A&A, 577, A19
  * van der Marel (1991) van der Marel, R. P. 1991, MNRAS, 253, 710
  * Vazza et al. (2013) Vazza, F., Eckert, D., Simionescu, A., Brüggen, M., & Ettori, S. 2013, MNRAS, 429, 799
  * Vianello et al. (2015) Vianello, G., Lauer, R. J., Younk, P., et al. 2015, arXiv e-prints, arXiv:1507.08343
  * Vikhlinin et al. (2006) Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691
  * Wechsler et al. (2002) Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., & Dekel, A. 2002, ApJ, 568, 52
  * White et al. (1993) White, S. D. M., Navarro, J. F., Evrard, A. E., & Frenk, C. S. 1993, Nature, 366, 429
  * Wilcox et al. (2015) Wilcox, H., Bacon, D., Nichol, R. C., et al. 2015, MNRAS, 452, 1171

##  Appendix A Fitting results for individual clusters

Figure 9: Same as Fig. 1 but for A85.

Figure 10: Same as Fig. 1 but for A644.

Figure 11: Same as Fig. 1 but for A1644.

Figure 12: Same as Fig. 1 but for A2029.

Figure 13: Same as Fig. 1 but for A2142.

Figure 14: Same as Fig. 1 but for A2255.

Figure 15: Same as Fig. 1 but for A2319.

Figure 16: Same as Fig. 1 but for A3158.

Figure 17: Same as Fig. 1 but for A3266.

Figure 18: Same as Fig. 1 but for RXC1825.

Figure 19: Same as Fig. 1 but for Zw1215.