Issue 
A&A
Volume 662, June 2022



Article Number  A123  
Number of page(s)  30  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202142507  
Published online  30 June 2022 
The gravitational field of XCOP galaxy clusters
^{1}
Department of Astronomy, University of Geneva, Ch. d’Ecogia 16, 1290 Versoix, Switzerland
email: Dominique.Eckert@unige.ch
^{2}
INAF, Osservatorio di Astrofisica e Scienza dello Spazio, Via Piero Gobetti 93/3, 40129 Bologna, Italy
^{3}
INFN, Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
^{4}
IRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France
^{5}
European Southern Observatory, KarlSchwarzschildStr. 2, 85748 Garching, Germany
^{6}
Centre for Space Research, NorthWest University, Potchefstroom 2520, South Africa
Received:
22
October
2021
Accepted:
25
April
2022
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 Λ 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 Xray 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 Λ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 ∼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 galaxycluster scale.
Key words: galaxies: clusters: general / dark matter / galaxies: clusters: intracluster medium / Xrays: galaxies: clusters / gravitation
© D. Eckert et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
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 socalled Navarro–Frenk–White (NFW) profile (Navarro et al. 1996, 1997), in which the matter density forms a central “cusp” with ρ(r)∝r^{−1} and gradually steepens toward the outskirts as ρ(r)∝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., dlnρ/dlnrr_{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 Nbody 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 powerlaw dependence of the logarithmic radial slope, such as dlnρ/dlnr ∝ r^{−α}. The Einasto (1965) profile is an example of this revised functional form, where the Einasto index, α, can be related to the powerlaw 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 selfinteraction (SIDM) cross section is nonnegligible, highdensity 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 freestream across the highdensity 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 (∼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, 2019; Bartalucci 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 (XCOP; Eckert et al. 2017), a very large program on XMMNewton that provides complete Xray 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 Xray 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 XMMNewton observations of a synthetic NFW cluster. We apply our technique to the 12 XCOP 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 (Eckert et al. 2022, hereafter Paper II), we discuss our measurements of the Einasto shape parameter and perform a detailed comparison to cosmological simulations that include selfinteracting DM.
Throughout the paper we assume a Planck 2015 ΛCDM cosmology (Planck Collaboration XIII 2016) with H_{0} = 67.8 km s^{−1} Mpc^{−1} and Ω_{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 easytouse Python package named hydromass. Our XMMNewton mock data generator can be downloaded from GitHub.
2. Sample and data reduction
2.1. The XCOP sample
XCOP (Eckert et al. 2017) is a very large program on XMMNewton, totaling 1.5 Ms of observing time. The original data set was awarded during XMMNewton AO13 (proposal ID 074441, PI: Eckert) and completed in 2015. The primary goal of the project is to use the combination of Xray data from XMMNewton 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 signaltonoise greater than 12, (ii) apparent size θ_{500} > 10′ 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 XMMNewton pointings that were obtained during AO17 (observation ID 082365, PI: Eckert) on two clusters (A3266 and A2029) to homogenize the data set.
2.2. XMMNewton 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)^{1} 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 XCOP. 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 mosfilter and pnfilter 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 sourcetobackground ratio and allows us to minimize the systematics associated with background subtraction. To model the highenergy particle background, we used a large collection of filterwheelclosed observations available in the calibration database, and rescaled the normalization of the closed data such that the highenergy 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 highenergy 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 blanksky pointings and presented in Appendix B of Ghirardini et al. (2018).
For each XCOP object, the count images from the various individual pointings were coadded to create a mosaic combining all the available data. The same procedure was repeated for the individual exposure maps and non Xray 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 Xray 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 Xray spectra in concentric annuli and extracted the appropriate response files and non Xray 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 threecomponent model (cosmic Xray 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 singletemperature thinplasma 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 allsky SZ map constructed in a very similar way as the public map of the full Planck survey delivered by the Planck Collaboration XXII (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 unitless Comptonization parameter, integrating the gas pressure along the line of sight,
The yprofiles and 3D pressure profiles are obtained following the method presented in Planck Collaboration V (2013), and subsequently used in studies by the XCOP 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.
yprofiles are computed over local maps of size 20 × θ_{500} on a side, centered on the XCOP cluster position and extracted from the allsky SZ map. The patches are oversampled with respect to the Planck SZ map pixel. The induced correlation between the points of the yprofile 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 × R_{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σ criterion with respect to the mean of the flux beyond 5 × R_{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, θ > 5 × R_{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 yprofile 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 yprofiles. 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 yprofile are propagated through the covariance matrix by generating 10 000 realizations of the noise profile. Each are coadded to the yprofile, 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 rband imaging obtained on the CanadaFranceHawaii 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 multiGaussian expansion (MGE) method, as implemented by Cappellari (2002), to obtain the stellar mass distributions from the rband imaging. The PSFconvolved 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 masstolight ratio, Υ_{∗DYN}. We assumed an inclination of i = 90° (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_{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 (β_{z}) and (ii) the stellar masstolight ratio (Υ_{∗DYN}). The BCG mass model was compared to the observed kinematic profiles and the two parameters β_{z} and Υ_{∗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 masstolight 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 masstolight ratio (Υ_{∗DYN}). When a DM mass component is included in the dynamical mass models for the BCGs, the best fitting Υ_{∗DYN} decrease by 8.3 ± 2.9% on average.
Our surface brightness measurement are accurate at the level of < 0.25 mag arcsec^{−2}, and our stellar kinematic measurement errors are ±4.7% 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_{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σ errors on the weak lensing masses as well as a ±10% 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 rband imaging data obtained with MegaCam at the CFHT, and additional u and iband photometry taken with the Wide Field Camera at the Isaac Newton Telescope. Aperture magnitude limits (5σ)are typically 24.3, 24.8, 24.2, and 23.3 in the ugrifilters, respectively, when measured on PSFhomogenized images using Gaussian weight functions. The stellar masstolight 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 starformation 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 ugrifilters. We also accounted for fieldtofield (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 lognormal mixture reconstruction (Sect. 3.4). We distribute our code jointly with this paper in the form of the publicly available Python package hydromass^{2}.
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,
with ρ_{gas}, P_{gas} the gas density and pressure, respectively, and M_{tot}(< r) the total (Newtonian) mass enclosed within a radius r. To take advantage of the sum of information provided by the combination of Xray 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,
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 Xray (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 multiscale approach implemented in Eckert et al. (2020). Namely, the 3D Xray emissivity profile is described as the linear combination of a large number N_{K} of King profiles,
with the dictionary of input King functions (King 1962),
and 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 β, with one value of r_{c} for every set of 4 surface brightness points and 10 values of β 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 α_{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 lineofsight average of the 3D temperature profile. The average temperature is weighted by the local emissivity . On top of that, the temperatures obtained through singletemperature fits to multiphase spectra tend to be lower than massweighted temperatures because of the response of Xray instruments. Mazzotta et al. (2004) showed that for temperatures of 3 keV and above, a condition that is satisfied by all XCOP clusters, the average spectroscopic temperature can be written as
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 Xray and SZ information by constructing a joint likelihood comparing the 3D model pressure to the SZ pressure and the spectroscopiclike temperature profile to the Xray 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,
with M_{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σ uncertainty.
In Ettori et al. (2019, hereafter E19) we applied a similar technique to the XCOP 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,
and the Einasto model,
Instead of the scale radius r_{s} and density normalization ρ_{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 α free, whereas previously its value was fixed to 1/α ≡ μ = 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 massconcentration relation.
Normal priors on the NFW fit parameters. Here P_{m} and dP_{m} denote the outermost SZ pressure value and its error.
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,
In this case, the DM profile ρ_{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),
with P_{0}, r_{s}, α, β, and γ as free parameters. This functional form was found to provide a good representation of gas pressure profiles determined from both Xray 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 lognormal 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 lognormal functions,
The mean values and standard deviations 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 μ_{i}, whereas the standard deviations σ_{i} act as effective smoothing scales. We implement the model with N_{g} = 200 and the values of μ_{i} logarithmically spaced from the innermost to the outermost data points. We set the values of σ_{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 lognormal functions N_{g}, the mean values μ_{i} and the standard deviations σ_{i} are adaptively set, the model can be projected onto the line of sight including spectroscopiclike weights (Eq. (6)) and the normalizations can be fit to the spectroscopic Xray temperatures and SZ pressure profile. The mass profile is then reconstructed by combining the 3D temperature and density profiles,
The temperature gradient can be computed analytically from Eq. (12),
Similarly, since the King functions are analytical functions, the gas density gradient can be computed analytically from Eq. (4),
3.5. Optimization
Since our models can contain several hundred parameters, we require the use of a statistical sampler that is suitable for highdimensional optimization problems. We use the No UTurn 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 selftuning 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 2000 tuning steps and 2000 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 lognormal 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 bottomright 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 XCOP 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, μ = 1/α (see Paper II). The details of the Einasto fitting procedure, the adopted priors, and the resulting bestfit parameters are provided in Paper II.
Fig. 1. Mass reconstruction results for A1795. Top left: electron density profile reconstructed with the multiscale method used here (Eckert et al. 2020) compared to the L1 regularization results presented in Ghirardini et al. (2019) and publicly released by the XCOP team. Top right: nonparametric reconstruction of the 3D temperature profile (green) compared to the spectroscopic Xray measurements (red) and the 3D temperature profile obtained by dividing the SZ pressure by the Xray 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 XMMNewton 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 Xray 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 Ω_{b}/Ω_{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.
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 XMMNewton 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 XCOP 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 × 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 XMMNewton 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)^{3} to generate model spectra in each cell. A constant metal abundance of 0.3 Z_{⊙} 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 × 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 XMMNewton 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 paper^{4}. The tool reads the onaxis effective area of each of the three XMMNewton instruments from the calibration database, and combines the onaxis area with the telescope’s energydependent 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 unvignetted component for the particleinduced background. The unvignetted component is drawn out of the collection of filterwheelclosed observations available in the calibration database. The vignetted component is made of a threecomponent model (see Sect. 2.2) consisting of an absorbed power law for the cosmic Xray 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 cosmicray 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 Xray 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 singletemperature APEC model and the sum of vignetted and unvignetted 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 lognormal mixture model.
The performance of our code on the mock data is shown in Fig. 2. In the lefthand 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 ∼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.
Fig. 2. Result of temperature deprojection (left) and HSE mass reconstruction (right) reconstructed from a set of ten mock XMMNewton 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. 
The comparison between the mean of the fitted mass profiles and the true profile can be seen in the righthand 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.
Overa ll, 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 (∼10%, Rasia et al. 2014) and effective area calibration (15%, Schellenberger et al. 2015).
5. Results
5.1. Total density profiles
We used our hydromass Python package introduced in Sect. 3 to analyze the gravitational field of the 12 XCOP clusters. In each case, we ran the reconstruction with the two mass models (NFW and Einasto), the parametric forward method, and the nonparametric lognormal 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.
Sample properties and results of NFW fits.
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.
Ratio between the masses published in E19 (M_{E19}) and our new results (labeled M_{hydromass}) for several overdensity factors.
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 righthand panels of the figure we show the reconstructed NFW and Einasto profiles, respectively, for the entire XCOP 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.
Fig. 3. Comparison between the mass profiles reconstructed using mass models (solid curves and shaded areas) and the nonparametric lognormal 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 Xray 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. 
This visual impression is confirmed in a quantitative way in the lefthand 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 forwardfitting method also closely traces the NP profile, with the exception of the innermost regions, where the powerlaw shape of the gNFW pressure profile leads to excess pressure, and hence excess mass, compared to the NP data. In the lefthand panel of Fig. 4 we show the median NPtomodel 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 XCOP 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 Λ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}.
Fig. 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 massconcentration relation for our sample. The small shaded ellipses indicate the 1σ contours of the fitted individual values for the 12 XCOP clusters individually, whereas the solid red and blue ellipses show the fitted sample mean and its 1σ and 2σ confidence regions, respectively. The curves show the predictions of Nbody simulations (Ludlow et al. 2016; Diemer & Joyce 2019; Ishiyama et al. 2021). 
5.2. Concentration–mass relation
The relation between NFW concentration and mass is a clear prediction of the Λ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 ΛCDM framework.
Given that our sample spans a narrow mass range (only a factor of ∼2), it is not sensitive to the shallow slope of the NFW massconcentration relation. However, it is well suited to measure the mean and scatter of the relation around M_{200} = 10^{15} M_{⊙}. To this end, we performed a Bayesian multivariate analysis of our sample, with the fractional lognormal 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 lognormal 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 halfCauchy prior with β = 0.5 is set on the fractional intrinsic scatter along both axes. Adopting instead a halfnormal 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σ contours for all XCOP 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 and with an intrinsic scatter .
Fig. 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 σ_{R200} and σ_{c} the fractional lognormal 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 Nbody simulation suites. The curves in the righthand 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 XCOP 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σ. 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 lognormal 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 σ_{lnc} ∼ 0.21 for more relaxed systems. For the most massive clusters, Bhattacharya et al. (2013) estimate σ_{lnc} = 0.33, in excellent agreement with our measurements. Diemer & Kravtsov (2015) measure a scatter of about 0.16 dex (σ_{lnc} = 0.37) at all redshifts, masses, and for all mass definition, when the entire ensemble of halos is considered, again in agreement with the XCOP 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 Λ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 XCOP clusters as traced by the Einasto shape parameter is provided in Paper II, where we present our measurements of α and discuss in detail the implications of our measurements.
5.3. Breaking down into baryonic and DM components
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 ∼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,
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 lefthand 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 ∼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 baryondriven 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 ∼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 ∼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).
Fig. 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, dashdotted 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. 
5.4. Cumulative baryon fraction
The data presented here also allow us to study the total cumulative baryon fraction in our systems as a function of radius. The righthand panel of Fig. 6 shows the integrated baryon fraction f_{bar}(r) = M_{bar}(< r)/M_{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 ∼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}).
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.
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 ≲20% level throughout the range considered here, which has little impact on our main conclusions.
5.5. The galaxy cluster radial acceleration 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_{obs} and the acceleration that would be expected in the absence of DM from the sum of the detected baryonic components, g_{bar}. Studying the rotation curves of rotationally supported galaxies, McGaugh et al. (2016) showed that the relation between g_{bar} and g_{obs} is nearly universal across the considered sample and smoothly deviates from the onetoone 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 XCOP galaxy clusters. The baryonic and observed acceleration were computed locally from the reconstructed profiles,
Fig. 7. Radial acceleration relation for XCOP 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 lognormal 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 onetoone relation. 
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,
with g_{†} = 1.2 × 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_{bar} ∼ 2 × 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 onetoone 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
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_{obs} − g_{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
with g_{1}, g_{2} the two characteristic acceleration scales describing the transition between the three regimes, and α_{1}, α_{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 when g_{bar} ≪ g_{2}, transitioning to when g_{2} < g_{bar} < g_{1} and converging to 1 when g_{bar} ≫ g_{1}. Our empirical formula approaches the MOND formalism at low acceleration when α_{1} = α_{2} = 0.5 and g_{1} is the classical MOND acceleration a_{0} ≈ 10^{−10} m s^{−2}.
We fitted the nonparametric points on the g_{obs} − g_{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 σ_{lngobs} = 0.21 ± 0.02. The bestfit curve is shown in the lefthand panel of Fig. 8 where the functional form and its uncertainty are compared to the data. For the two acceleration scales we obtain m s^{−2} and m s^{−2}, in qualitative agreement with the above discussion. The fitted slopes yield and . The scale g_{1} is of the same order as the classical MOND acceleration, albeit slightly higher. However, the fitted slope α_{1} is considerably larger than its corresponding MOND counterpart, which indicates that the transition away from the baryondominated regime is much faster than in disk galaxies. Below the second characteristic scale g_{2} the g_{obs} − g_{bar} relation is much steeper than expected in MOND.
Fig. 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σ 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). 
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 XCOP 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_{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_{obs} and g_{bar} we numerically search for the value of M_{DM, RAR} that satisfies the condition
We used the Einasto reconstructions to determine the mass profiles of M_{DM, RAR}. In the righthand 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 ∼0.05R_{500}, similar to the CDM case (Fig. 6). The required enclosed mass exceeds the gas mass and reaches a maximum around ∼0.5R_{500} (∼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 2022). Through a direct comparison between hydrostatic masses and masses obtained using alternative methods (Ettori et al. 2019), we found that XCOP 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 Λ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 Xray densities and temperatures are expected to be biased toward high and low values, respectively, compared to the massweighted quantities. The multitemperature 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 Xray derived temperatures, which might greatly increase the bias in the recovered HSE masses. For instance, Pearce et al. (2020) predict that masses based on Xray data only may be largely biased low (up to 50% in the XCOP mass range), whereas HSE masses using SZderived pressure profiles should be closer to the true, massweighted 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 Xray data only are consistent with the results obtained in combination with SZ data (Ettori et al. 2019) and the pressure profiles determined independently from Xray 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 XMMNewton effective area affects the measured spectroscopic temperatures, as the spectral model needs to be folded through the instrumental response to fit the Xray spectra. There is a known inconsistency between the temperatures determined by XMMNewton/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 ∼20% for 10 keV plasma. It is still unclear whether the XMMNewton 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 XMMNewton and Planck data, the latter providing the dominant contribution to the temperature profiles at large radii (≥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 Xray 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 ΛCDM framework
All of the results presented in this paper agree with the predictions of the ΛCDM framework and the bottomup structure formation paradigm. The average mass profile of the XCOP 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 lefthand panel of Fig. 4). The average NFW concentration of the sample matches perfectly the ΛCDM predictions (righthand panel of Fig. 4), and Nbody simulations also reproduce accurately the observed scatter in the massconcentration relation at M_{200} ∼ 10^{15} M_{⊙}. 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 Ω_{b}/Ω_{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 Λ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 Xray morphology as indicated by their large centroid shifts (see Paper II); thus, it is possible that the Xray 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 Xray 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 ∼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 Λ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 XCOP clusters at very high significance. As described in Sect. 5.5, we observe three distinct regimes in the g_{obs} − g_{bar} plane: at high acceleration (g_{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_{bar} ∼ 2 × 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 XCOP 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 ∼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 σ_{lng} = 0.21, the behavior of the recovered RAR in XCOP 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_{bar} − g_{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 Λ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 Λ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 XCOP galaxy clusters over the radial range [0.01 − 2]R_{500} from a combination of deep XMMNewton and Planck data. Our results can be summarized as follows:

We introduced a novel framework for the reconstruction of hydrostatic mass profiles from Xray 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 lognormal functions. Our framework improves over our previous works (Ettori et al. 2019; the differences in the mass values are less than 10%; see Table 3) as it implements PSF deconvolution, affecting in particular the inner regions, and fits with the Einasto mass model, where the parameter α 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 XMMNewton 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 XCOP 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 Λ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 ± 0.39 at M_{200} ≈ 10^{15} M_{⊙}, in remarkable agreement with ΛCDM predictions. The intrinsic scatter of the massconcentration relation was found to be , which again is consistent with the expectations of Nbody 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 ∼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 ∼5 around g_{bar} = 2 × 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., CHEXMATE, The CHEXMATE Collaboration 2020) will allow us to study in detail the shape of galaxy cluster mass profiles and test the validity of the HSE assumption.
Acknowledgments
S.E. acknowledges financial contribution from the contracts ASIINAF Athena 201927HH.0, “Attività di Studio per la comunità scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASIINAF n. 201714H.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). S.I.L. 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, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]
 Angelinelli, M., Vazza, F., Giocoli, C., et al. 2020, MNRAS, 495, 864 [NASA ADS] [CrossRef] [Google Scholar]
 Angus, G. W., Famaey, B., & Buote, D. A. 2008, MNRAS, 387, 1470 [CrossRef] [Google Scholar]
 Angus, G. W., Famaey, B., & Diaferio, A. 2010, MNRAS, 402, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Angus, G. W., Diaferio, A., Famaey, B., & van der Heyden, K. J. 2013, MNRAS, 436, 202 [NASA ADS] [CrossRef] [Google Scholar]
 Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92 [CrossRef] [EDP Sciences] [Google Scholar]
 Barnes, D. J., Vogelsberger, M., Pearce, F. A., et al. 2021, MNRAS, 506, 2533 [NASA ADS] [CrossRef] [Google Scholar]
 Bartalucci, I., Arnaud, M., Pratt, G. W., Démoclès, J., & Lovisari, L. 2019, A&A, 628, A86 [EDP Sciences] [Google Scholar]
 Bekenstein, J. D. 2004, Phys. Rev. D, 70, 083509 [Google Scholar]
 Bennett, J. S., & Sijacki, D. 2022, MNRAS, 514, 313 [CrossRef] [Google Scholar]
 Bhattacharya, S., Habib, S., Heitmann, K., & Vikhlinin, A. 2013, ApJ, 766, 32 [Google Scholar]
 Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Biviano, A., Rosati, P., Balestra, I., et al. 2013, A&A, 558, A1 [Google Scholar]
 Bourdin, H., Mazzotta, P., Kozmanyan, A., Jones, C., & Vikhlinin, A. 2017, ApJ, 843, 72 [Google Scholar]
 Brown, S. T., McCarthy, I. G., Diemer, B., et al. 2020, MNRAS, 495, 4994 [NASA ADS] [CrossRef] [Google Scholar]
 Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
 Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559 [Google Scholar]
 Burgess, J. M., Fleischhack, H., Vianello, G., et al. 2021, https://doi.org/10.5281/zenodo.5646954 [Google Scholar]
 Cappellari, M. 2002, MNRAS, 333, 400 [NASA ADS] [CrossRef] [Google Scholar]
 Cappellari, M. 2008, MNRAS, 390, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345 [Google Scholar]
 Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
 Chan, M. H., & Del Popolo, A. 2020, MNRAS, 492, 5865 [NASA ADS] [CrossRef] [Google Scholar]
 Chiu, I., Mohr, J. J., McDonald, M., et al. 2018, MNRAS, 478, 3072 [Google Scholar]
 Croston, J. H., Arnaud, M., Pointecouteau, E., & Pratt, G. W. 2006, A&A, 459, 1007 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Diemer, B. 2018, ApJS, 239, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Diemer, B., & Joyce, M. 2019, ApJ, 871, 168 [NASA ADS] [CrossRef] [Google Scholar]
 Diemer, B., & Kravtsov, A. V. 2014, ApJ, 789, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Diemer, B., & Kravtsov, A. V. 2015, ApJ, 799, 108 [Google Scholar]
 Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008, MNRAS, 390, L64 [Google Scholar]
 Duffy, A. R., Schaye, J., Kay, S. T., et al. 2010, MNRAS, 405, 2161 [NASA ADS] [Google Scholar]
 Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359 [Google Scholar]
 Dutton, A. A., Macciò, A. V., Obreja, A., & Buck, T. 2019, MNRAS, 485, 1886 [NASA ADS] [CrossRef] [Google Scholar]
 Eckert, D., Molendi, S., Owers, M., et al. 2014, A&A, 570, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eckert, D., Roncarelli, M., Ettori, S., et al. 2015, MNRAS, 447, 2198 [Google Scholar]
 Eckert, D., Ettori, S., Pointecouteau, E., et al. 2017, Astron. Nach., 338, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eckert, D., Finoguenov, A., Ghirardini, V., et al. 2020, The Open J. Astrophys., 3, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Eckert, D., Ettori, S., Robertson, A., et al. 2022, A&A, submitted (Paper II) [Google Scholar]
 Einasto, J. 1965, Trudy Astrofizicheskogo Instituta AlmaAta, 5, 87 [NASA ADS] [Google Scholar]
 Eke, V. R., Navarro, J. F., & Frenk, C. S. 1998, ApJ, 503, 569 [NASA ADS] [CrossRef] [Google Scholar]
 Ettori, S., Morandi, A., Tozzi, P., et al. 2009, A&A, 501, 61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ettori, S., Gastaldello, F., Leccardi, A., et al. 2010, A&A, 524, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ettori, S., Ghirardini, V., Eckert, D., Dubath, F., & Pointecouteau, E. 2017, MNRAS, 470, L29 [CrossRef] [Google Scholar]
 Ettori, S., Ghirardini, V., Eckert, D., et al. 2019, A&A, 621, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Evrard, A. E. 1997, MNRAS, 292, 289 [CrossRef] [Google Scholar]
 Famaey, B., & McGaugh, S. S. 2012, Liv. Rev. Relat., 15, 10 [Google Scholar]
 Fasano, G., Bettoni, D., Ascaso, B., et al. 2010, MNRAS, 404, 1490 [NASA ADS] [Google Scholar]
 Geller, M. J., Hwang, H. S., Diaferio, A., et al. 2014, ApJ, 783, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ghizzardi, S., Molendi, S., van der Burg, R., et al. 2021, A&A, 646, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gonzalez, A. H., Zaritsky, D., & Zabludoff, A. I. 2007, ApJ, 666, 147 [Google Scholar]
 Haridasu, B. S., Karmakar, P., De Petris, M., Cardone, V. F., & Maoli, R. 2021, ArXiv eprints [arXiv:2111.01101] [Google Scholar]
 Harikumar, S., & Biesiada, M. 2022, Eur. Phys. J. C, 82, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Herbonnet, R., Sifón, C., Hoekstra, H., et al. 2020, MNRAS, 497, 4684 [NASA ADS] [CrossRef] [Google Scholar]
 Hodson, A. O., & Zhao, H. 2017, A&A, 598, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hoffman, M. D., & Gelman, A. 2014, J. Mach. Learn. Res., 15, 1593 [Google Scholar]
 Hurier, G., MacíasPérez, J. F., & Hildebrandt, S. 2013, A&A, 558, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ishiyama, T., Prada, F., Klypin, A. A., et al. 2021, MNRAS, 506, 4210 [NASA ADS] [CrossRef] [Google Scholar]
 Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 King, I. 1962, AJ, 67, 471 [Google Scholar]
 Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340 [Google Scholar]
 Kravtsov, A. V., Nagai, D., & Vikhlinin, A. A. 2005, ApJ, 625, 588 [Google Scholar]
 Laganá, T. F., Martinet, N., Durret, F., et al. 2013, A&A, 555, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129 [NASA ADS] [CrossRef] [Google Scholar]
 Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2014, MNRAS, 441, 1270 [NASA ADS] [CrossRef] [Google Scholar]
 Leccardi, A., & Molendi, S. 2008, A&A, 486, 359 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017, ApJ, 836, 152 [Google Scholar]
 Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2018, A&A, 615, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Loubser, S. I., Hoekstra, H., Babul, A., & O’Sullivan, E. 2018, MNRAS, 477, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Loubser, S. I., Babul, A., Hoekstra, H., et al. 2020, MNRAS, 496, 1857 [Google Scholar]
 Ludlow, A. D., & Angulo, R. E. 2017, MNRAS, 465, L84 [Google Scholar]
 Ludlow, A. D., Bose, S., Angulo, R. E., et al. 2016, MNRAS, 460, 1214 [Google Scholar]
 Ludlow, A. D., BenítezLlambay, A., Schaller, M., et al. 2017, Phys. Rev. Lett., 118, 161103 [NASA ADS] [CrossRef] [Google Scholar]
 Macciò, A. V., Dutton, A. A., & van den Bosch, F. C. 2008, MNRAS, 391, 1940 [Google Scholar]
 Macciò, A. V., Ruchayskiy, O., Boyarsky, A., & MuñozCuartas, J. C. 2013, MNRAS, 428, 882 [CrossRef] [Google Scholar]
 Mamon, G. A., & Łokas, E. L. 2005, MNRAS, 362, 95 [Google Scholar]
 Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2014, MNRAS, 440, 2077 [NASA ADS] [CrossRef] [Google Scholar]
 Mantz, A. B., Morris, R. G., Allen, S. W., et al. 2022, MNRAS, 510, 131 [Google Scholar]
 Mathiesen, B., Evrard, A. E., & Mohr, J. J. 1999, ApJ, 520, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Mazzotta, P., Rasia, E., Moscardini, L., & Tormen, G. 2004, MNRAS, 354, 10 [NASA ADS] [CrossRef] [Google Scholar]
 McCammon, D., Almy, R., Apodaca, E., et al. 2002, ApJ, 576, 188 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Phys. Rev. Lett., 117, 201101 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 1983, ApJ, 270, 384 [Google Scholar]
 Milgrom, M. 2009, Phys. Rev. D, 80, 123536 [NASA ADS] [CrossRef] [Google Scholar]
 Mitchell, M. A., He, J.H., Arnold, C., & Li, B. 2018, MNRAS, 477, 1133 [NASA ADS] [CrossRef] [Google Scholar]
 Montes, M., & Trujillo, I. 2014, ApJ, 794, 137 [Google Scholar]
 Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.W. 2011, ApJ, 731, 113 [Google Scholar]
 Munari, E., Biviano, A., & Mamon, G. A. 2014, A&A, 566, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJS, 206, 8 [Google Scholar]
 Nagai, D., & Lau, E. T. 2011, ApJ, 731, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Nagai, D., Vikhlinin, A., & Kravtsov, A. V. 2007, ApJ, 655, 98 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]
 Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039 [Google Scholar]
 Navarro, J. F., BenítezLlambay, A., Fattahi, A., et al. 2017, MNRAS, 471, 1841 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, K., Lau, E. T., & Nagai, D. 2014, ApJ, 792, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Neto, A. F., Gao, L., Bett, P., et al. 2007, MNRAS, 381, 1450 [NASA ADS] [CrossRef] [Google Scholar]
 Nevalainen, J., David, L., & Guainazzi, M. 2010, A&A, 523, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nieuwenhuizen, T. M., & Morandi, A. 2013, MNRAS, 434, 2679 [NASA ADS] [CrossRef] [Google Scholar]
 Okabe, N., Smith, G. P., Umetsu, K., Takada, M., & Futamase, T. 2013, ApJ, 769, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Pearce, F. A., Kay, S. T., Barnes, D. J., Bower, R. G., & Schaller, M. 2020, MNRAS, 491, 1622 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration V. 2013, A&A, 550, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIX. 2014, A&A, 571, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2016, A&A, 594, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planelles, S., Borgani, S., Dolag, K., et al. 2013, MNRAS, 431, 1487 [Google Scholar]
 Pointecouteau, E., & Silk, J. 2005, MNRAS, 364, 654 [NASA ADS] [CrossRef] [Google Scholar]
 Pointecouteau, E., Arnaud, M., & Pratt, G. W. 2005, A&A, 435, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pointecouteau, E., SantiagoBautista, I., Douspis, M., et al. 2021, A&A, 651, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pradyumna, S., Gupta, S., Seeram, S., & Desai, S. 2021, Phys. Dark Universe, 31 [Google Scholar]
 Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25 [Google Scholar]
 Ragagnin, A., Saro, A., Singh, P., & Dolag, K. 2021, MNRAS, 500, 5056 [Google Scholar]
 Rasia, E., Tormen, G., & Moscardini, L. 2004, MNRAS, 351, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Rasia, E., Ettori, S., Moscardini, L., et al. 2006, MNRAS, 369, 2013 [CrossRef] [Google Scholar]
 Rasia, E., Meneghetti, M., Martino, R., et al. 2012, New J. Phys., 14, 055018 [Google Scholar]
 Rasia, E., Lau, E. T., Borgani, S., et al. 2014, ApJ, 791, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Robertson, A., Massey, R., Eke, V., Schaye, J., & Theuns, T. 2021, MNRAS, 501, 4610 [NASA ADS] [CrossRef] [Google Scholar]
 Rossetti, M., Gastaldello, F., Ferioli, G., et al. 2016, MNRAS, 457, 4515 [Google Scholar]
 Salucci, P., & Burkert, A. 2000, ApJ, 537, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Comput. Sci., 2 [Google Scholar]
 Sand, D. J., Graham, M. L., Bildfell, C., et al. 2012, ApJ, 746, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H. 1999, ApJ, 512, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H. 2003, MNRAS, 342, 901 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H. 2007, MNRAS, 380, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Sayers, J., Czakon, N. G., Mantz, A., et al. 2013, ApJ, 768, 177 [Google Scholar]
 Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Skordis, C., & Złośnik, T. 2021, Phys. Rev. Lett., 127, 161302 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91 [Google Scholar]
 Smith, R. J., Lucey, J. R., & Edge, A. C. 2017, MNRAS, 471, 383 [NASA ADS] [CrossRef] [Google Scholar]
 Snowden, S. L., Mushotzky, R. F., Kuntz, K. D., & Davis, D. S. 2008, A&A, 478, 615 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tam, S.I., Jauzac, M., Massey, R., et al. 2020, MNRAS, 496, 4032 [NASA ADS] [CrossRef] [Google Scholar]
 Tchernin, C., Eckert, D., Ettori, S., et al. 2016, A&A, 595, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 The CHEXMATE Collaboration (Arnaud, M., et al.) 2020, A&A, 650, A104 [Google Scholar]
 Tian, Y., Umetsu, K., Ko, C.M., Donahue, M., & Chiu, I. N. 2020, ApJ, 896, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Umetsu, K., & Diemer, B. 2017, ApJ, 836, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Umetsu, K., Zitrin, A., Gruen, D., et al. 2016, ApJ, 821, 116 [Google Scholar]
 van der Burg, R. F. J., Hoekstra, H., Muzzin, A., et al. 2015, A&A, 577, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van der Marel, R. P. 1991, MNRAS, 253, 710 [NASA ADS] [CrossRef] [Google Scholar]
 Vazza, F., Eckert, D., Simionescu, A., Brüggen, M., & Ettori, S. 2013, MNRAS, 429, 799 [Google Scholar]
 Vianello, G., Lauer, R. J., Younk, P., et al. 2015, ArXiv eprints [arXiv:1507.08343] [Google Scholar]
 Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691 [Google Scholar]
 Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., & Dekel, A. 2002, ApJ, 568, 52 [NASA ADS] [CrossRef] [Google Scholar]
 White, S. D. M., Navarro, J. F., Evrard, A. E., & Frenk, C. S. 1993, Nature, 366, 429 [Google Scholar]
 Wilcox, H., Bacon, D., Nichol, R. C., et al. 2015, MNRAS, 452, 1171 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Fitting results for individual clusters
Fig. A.10. Same as Fig. 1 but for RXC1825. 
Fig. A.11. Same as Fig. 1 but for Zw1215. 
All Tables
Normal priors on the NFW fit parameters. Here P_{m} and dP_{m} denote the outermost SZ pressure value and its error.
Ratio between the masses published in E19 (M_{E19}) and our new results (labeled M_{hydromass}) for several overdensity factors.
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 Figures
Fig. 1. Mass reconstruction results for A1795. Top left: electron density profile reconstructed with the multiscale method used here (Eckert et al. 2020) compared to the L1 regularization results presented in Ghirardini et al. (2019) and publicly released by the XCOP team. Top right: nonparametric reconstruction of the 3D temperature profile (green) compared to the spectroscopic Xray measurements (red) and the 3D temperature profile obtained by dividing the SZ pressure by the Xray 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 XMMNewton 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 Xray 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 Ω_{b}/Ω_{m} in Planck Collaboration XIII (2016) cosmology. Bottom right: posterior distributions and correlations for the parameters of the NFW profile. 

In the text 
Fig. 2. Result of temperature deprojection (left) and HSE mass reconstruction (right) reconstructed from a set of ten mock XMMNewton 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. 

In the text 
Fig. 3. Comparison between the mass profiles reconstructed using mass models (solid curves and shaded areas) and the nonparametric lognormal 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 Xray 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 the text 
Fig. 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 massconcentration relation for our sample. The small shaded ellipses indicate the 1σ contours of the fitted individual values for the 12 XCOP clusters individually, whereas the solid red and blue ellipses show the fitted sample mean and its 1σ and 2σ confidence regions, respectively. The curves show the predictions of Nbody simulations (Ludlow et al. 2016; Diemer & Joyce 2019; Ishiyama et al. 2021). 

In the text 
Fig. 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 σ_{R200} and σ_{c} the fractional lognormal intrinsic scatter of the points around the mean value. 

In the text 
Fig. 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, dashdotted 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. 

In the text 
Fig. 7. Radial acceleration relation for XCOP 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 lognormal 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 onetoone relation. 

In the text 
Fig. 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σ 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). 

In the text 
Fig. A.1. Same as Fig. 1 but for A85. 

In the text 
Fig. A.2. Same as Fig. 1 but for A644. 

In the text 
Fig. A.3. Same as Fig. 1 but for A1644. 

In the text 
Fig. A.4. Same as Fig. 1 but for A2029. 

In the text 
Fig. A.5. Same as Fig. 1 but for A2142. 

In the text 
Fig. A.6. Same as Fig. 1 but for A2255. 

In the text 
Fig. A.7. Same as Fig. 1 but for A2319. 

In the text 
Fig. A.8. Same as Fig. 1 but for A3158. 

In the text 
Fig. A.9. Same as Fig. 1 but for A3266. 

In the text 
Fig. A.10. Same as Fig. 1 but for RXC1825. 

In the text 
Fig. A.11. Same as Fig. 1 but for Zw1215. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.