© D. Eckert et al. 2022

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 so-called Navarro–Frenk–White (NFW) profile (Navarro et al. 1996, 1997), in which the matter density forms a central “cusp” with ρ(r)∝r⁻¹ and gradually steepens toward the outskirts as ρ(r)∝r⁻³, 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ρ/dlnr|r_s = −2). The characteristic shape is thought to arise from the universality of the structure formation process (Bullock et al. 2001; Macciò et al. 2008), and the distribution and redshift evolution of the halo structural parameters encode important cosmological information (e.g., Duffy et al. 2008; Dutton & Macciò 2014). As discussed in the original work by Navarro, Frenk, and White, the total mass at a given overdensity (typically 200 times the critical density of the universe enclosed in a sphere with radius R₂₀₀) correlates with the mass concentration of the halo, c = R₂₀₀/r_s. This correlation reflects the conditions in which a halo assembled, with halos formed earlier being more concentrated because the background density was larger back in cosmological time (Wechsler et al. 2002). More recently, several works (see, e.g., Dutton & Macciò 2014; Ludlow et al. 2016; Ragagnin et al. 2021) have modeled and quantified the effects of halo accretion histories and of the parameters that describe the cosmological background model on the halo mass concentration. Recent N-body simulations have also shown that the radial profiles of massive halos deviate from the NFW prediction (Navarro et al. 2004; Diemer & Kravtsov 2014; Klypin et al. 2016; Ludlow et al. 2016) and can be better fitted by functional forms that include a power-law dependence of the logarithmic radial slope, such as dlnρ/dlnr ∝ r^−α. The Einasto (1965) profile is an example of this revised functional form, where the Einasto index, α, can be related to the power-law index of the primordial power spectrum, n_s (Ludlow & Angulo 2017; Brown et al. 2020).

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

As the most massive collapsed structures in the Universe, galaxy clusters are privileged sites for studying the properties of the gravitational field over a wide range of densities. Given their deep potential well, gravitational processes largely dominate over feedback processes, such that galaxy clusters are essentially closed boxes for baryons (White et al. 1993; Evrard 1997; Eke et al. 1998; Kravtsov et al. 2005). Recent studies demonstrate that the baryon fraction enclosed within the virial radius of massive clusters is very close to the universal baryon fraction (Gonzalez et al. 2007; Laganá et al. 2013; Mantz et al. 2014; Chiu et al. 2018; Eckert et al. 2019), with the majority of the baryonic mass (∼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 (X-COP; Eckert et al. 2017), a very large program on XMM-Newton that provides complete X-ray mapping out to the virial radius for a sample of 12 nearby systems selected from the Planck Sunyaev–Zel’dovich (SZ) survey. We present a novel nonparametric technique for recovering hydrostatic mass profiles from joint X-ray and SZ data, which we include in a global framework for the estimation of hydrostatic mass profiles. In Sect. 4 we validate our reconstruction technique using mock XMM-Newton observations of a synthetic NFW cluster. We apply our technique to the 12 X-COP systems with the aim of comparing NFW and Einasto profiles. For a subset of systems, we use measurements of the gas mass as well as the stellar mass within the brightest cluster galaxy (BCG) and satellite galaxies to break down the gravitational field into its baryonic and DM components. Finally, we study the relation between the baryonic and total gravitational acceleration and compare it with the radial acceleration relation (RAR) of rotationally supported galaxies (McGaugh et al. 2016). In a companion paper (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 self-interacting DM.

Throughout the paper we assume a Planck 2015 ΛCDM cosmology (Planck Collaboration XIII 2016) with H₀ = 67.8 km s⁻¹ Mpc⁻¹ 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 easy-to-use Python package named hydromass. Our XMM-Newton mock data generator can be downloaded from GitHub.

2. Sample and data reduction

2.1. The X-COP sample

X-COP (Eckert et al. 2017) is a very large program on XMM-Newton, totaling 1.5 Ms of observing time. The original data set was awarded during XMM-Newton AO-13 (proposal ID 074441, PI: Eckert) and completed in 2015. The primary goal of the project is to use the combination of X-ray data from XMM-Newton with SZ data from Planck to probe the state of the ICM out to the virial radius. The sample was selected from the first Planck SZ catalog (PSZ1; Planck Collaboration XXIX 2014) according to the following criteria: (i) Planck PSZ1 signal-to-noise greater than 12, (ii) apparent size θ₅₀₀ > 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²¹ cm⁻².

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

2.2. XMM-Newton data analysis

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

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

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

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

2.3. Planck data analysis

The pressure profiles are recovered from the Planck survey data. We made use of a custom all-sky SZ map constructed in a very similar way as the public map of the full Planck survey delivered by the Planck Collaboration 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 unit-less Comptonization parameter, integrating the gas pressure along the line of sight,

(1)

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

y-profiles are computed over local maps of size 20 × θ₅₀₀ on a side, centered on the X-COP cluster position and extracted from the all-sky SZ map. The patches are over-sampled with respect to the Planck SZ map pixel. The induced correlation between the points of the y-profile is accounted for and propagated through the covariance matrix of the profile. We performed an azimuthal mean of the pixel fluxes falling in each radial bin, and subtracted a background residual offset estimated from radii larger than 5 × R₅₀₀ 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₅₀₀. 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₅₀₀. From the derived noise power spectrum, we draw a thousand realizations of the noise over each patch. Each realization is folded in the procedure to extract the y-profile and then used to compute the profile covariance matrix. The latter hence propagates the intrinsic properties of the noise in the SZ Planck map as well as the correlation introduced by our choice of pixel size.

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

2.4. Stellar mass

2.4.1. BCG

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

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

The surface brightness obtained from the MGE method was converted to surface brightness density, and the deprojection from surface density to intrinsic density was performed for the axisymmetric case using the Jeans anisotropic method (Cappellari 2008) assuming a constant stellar mass-to-light ratio, Υ_∗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₂₀₀ 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 mass-to-light 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 mass-to-light ratio to derive the enclosed stellar mass of the BCG within spheres of different radii.

For A644 and A2319, we did not have weak lensing masses available from Herbonnet et al. (2020) but used the average DM distribution (of the 23 other clusters in Loubser et al. 2020) to calculate the effect of a DM mass component on the best fitting stellar mass-to-light ratio (Υ_∗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⁻², 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₂₀₀ (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 r-band imaging data obtained with MegaCam at the CFHT, and additional u- and i-band photometry taken with the Wide Field Camera at the Isaac Newton Telescope. Aperture magnitude limits (5σ)are typically 24.3, 24.8, 24.2, and 23.3 in the ugri-filters, respectively, when measured on PSF-homogenized images using Gaussian weight functions. The stellar mass-to-light ratios that form the basis of the stellar mass measurements are obtained using spectral energy distribution fitting of these aperture fluxes. For this we assume the stellar population libraries from Bruzual & Charlot (2003), an exponentially declining star-formation history, and a Chabrier (2003) initial mass function. Further details are provided in van der Burg et al. (2015).

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

3. Mass modeling

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

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,

(2)

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 X-ray and SZ data, we performed joint fits to the two data sets by constructing 3D model profiles for the gas density, temperature, and pressure, which are related to one another through the ideal gas equation of state,

(3)

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

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

(4)

with the dictionary of input King functions (King 1962),

(5)

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 line-of-sight average of the 3D temperature profile. The average temperature is weighted by the local emissivity . On top of that, the temperatures obtained through single-temperature fits to multiphase spectra tend to be lower than mass-weighted temperatures because of the response of X-ray instruments. Mazzotta et al. (2004) showed that for temperatures of 3 keV and above, a condition that is satisfied by all X-COP clusters, the average spectroscopic temperature can be written as

(6)

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

3.2. Mass model

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

(7)

with M_mod(r) the model density profile integrated over the volume out to radius r, r₀ the outermost radius out to which the pressure can be measured, and P₀ = P(r₀) 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₀ is left free while fitting. We set a Gaussian prior on P₀ 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 X-COP sample and compared five different model profiles. We found that models with a central cusp are usually preferred by the data over cored models. For this reason, here we focus on the NFW model,

(8)

and the Einasto model,

(9)

Instead of the scale radius r_s and density normalization ρ_s we optimize for the overdensity radius R₂₀₀ and the concentration c₂₀₀ = R₂₀₀/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 mass-concentration relation.

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,

(10)

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),

(11)

with P₀, r_s, α, β, and γ as free parameters. This functional form was found to provide a good representation of gas pressure profiles determined from both X-ray and SZ observations (Arnaud et al. 2010; Planck Collaboration V 2013; Sayers et al. 2013; Bourdin et al. 2017; Ghirardini et al. 2019; Pointecouteau et al. 2021). The pressure gradient can be computed analytically from the reconstructed functional form to recover the mass profile from Eq. (2).

3.4. Nonparametric log-normal mixture reconstruction

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

(12)

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 log-normal 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 spectroscopic-like weights (Eq. (6)) and the normalizations can be fit to the spectroscopic X-ray temperatures and SZ pressure profile. The mass profile is then reconstructed by combining the 3D temperature and density profiles,

(13)

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

(14)

Similarly, since the King functions are analytical functions, the gas density gradient can be computed analytically from Eq. (4),

(15)

3.5. Optimization

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

For each cluster and model, we start by searching for the maximum likelihood using a gradient descent method, and then run four parallel HMC chains with 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 log-normal mixture method (hereafter NP). All the methods provide similar results across most of the radial range, although differences arise primarily in the innermost 100 kpc. Since the NP method is designed to trace the observed values as closely as possible, it cannot be extrapolated outside of the fitted range. Thus, for the remainder of the paper we evaluate it at the radii of the observed data and present them in the form of data points rather than continuous models. In the bottom-right panel we show the posterior distributions for the parameters of the NFW profile and the correlations between the parameters. Given the excellent statistical quality of the X-COP data, the posterior distributions are much narrower than the prior and clearly pull away from it, indicating that our data are able to break the degeneracy between NFW mass and concentration. To help convergence of the Einasto model, we optimize for the inverse of the Einasto index, μ = 1/α (see Paper II). The details of the Einasto fitting procedure, the adopted priors, and the resulting best-fit parameters are provided in Paper II.

Fig. 1. Mass reconstruction results for A1795. Top left: electron density profile reconstructed with the multi-scale method used here (Eckert et al. 2020) compared to the L1 regularization results presented in Ghirardini et al. (2019) and publicly released by the X-COP team. Top right: nonparametric reconstruction of the 3D temperature profile (green) compared to the spectroscopic X-ray measurements (red) and the 3D temperature profile obtained by dividing the SZ pressure by the X-ray density (orange). See Sect. 3.4 for details. The blue curve shows the projected, spectroscopically weighted, and PSF convolved temperature profile, to be compared to the red data points. Middle left: electron pressure profile measured by Planck and XMM-Newton compared to the reconstructions obtained from the NFW and Einasto mass models (Sect. 3.2), the parametric forward approach (labeled “Forward”; see Sect. 3.3), and the NP reconstruction (Sect. 3.4). The dotted vertical lines show the location of R500 and R200 from the NFW fit. Middle right: output mass profiles obtained from the four individual reconstructions (NFW, Einasto, Forward, and NP). The data points show the NP results evaluated at the radius of the X-ray spectroscopic measurements (labeled “NP Spectral bins,” blue) and the SZ pressure data (“NP SZ bins,” red). Bottom left: same as previous but for the gas mass fraction. The gray shaded area shows the cosmic baryon fraction Ω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₂₀₀ = 2000 kpc at a redshift z = 0.2, such that the simulated cluster fits well within the XMM-Newton FOV. The halo is filled with an ICM in HSE within the input NFW potential well. The input gas density profile follows the “universal” average X-COP gas density profile (see Fig. 3 and Table 3 of Ghirardini et al. 2019), which was fitted to a collection of individual nonparametric gas density profiles using a single common functional form adapted from Vikhlinin et al. (2006) and jointly accounting for the intrinsic scatter among the cluster population.

Given the input NFW mass profile and the universal gas density profile, we used Eq. (7) to compute the corresponding thermal pressure profile. The profiles were computed out to 5 × R₅₀₀, which is well beyond the range covered by our data to avoid uncertainties linked with the integration constant P₀. We then generated a 3D box inside which we set the gas density and temperature of the ICM inside each cell following spherical symmetry. The cell size was set to 2 kpc, which is substantially smaller than the XMM-Newton resolution. Inside each box cell, we used the APEC model (Smith et al. 2001) as implemented in 3ML (Vianello et al. 2015; Burgess et al. 2021)³ 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²⁰ cm⁻². Finally, the 3D model spectra were projected onto the line of sight. As a result, a spectrum made of the linear combination of individual APEC spectra is associated with each individual pixel of a projected 2D map.

4.2. Mock data creation

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

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

4.3. Code performance

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

The performance of our code on the mock data is shown in Fig. 2. In the left-hand panel we show the mean and dispersion of the ten deprojected 3D temperature profiles compared with the input profile. The nonparametric profile was evaluated at the average radii of the spectroscopic bins. We can see that all three methods are able to trace the true profile with minimal bias (i.e., less than 3% at each radius). In particular, the nonparametric method provides a good description of the profile with a scatter of ∼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 XMM-Newton observations of a synthetic NFW cluster (see Sect. 4). The curves show the mean and standard deviation of the ten realizations for the various reconstruction techniques (NFW, green; Einasto, cyan; and NP, red). The bottom panels show the ratio of reconstructed profiles to true input profiles.

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

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

We started by comparing the output of our new framework to the masses previously published in E19 in the case of the “standard” NFW reconstruction. As explained in Sect. 3.2, the main difference with respect to our previous analysis is the implementation of a modeling scheme for the PSF, both for the surface brightness and the temperature profile. Correcting for PSF smearing increases the concentration parameter of the NFW profile, which results in slightly higher masses in the core and lower masses in the outskirts. The average ratio between the E19 masses and the results presented here is given in Table 3 for several overdensities. Our new masses are slightly higher inside R₂₅₀₀ and slightly lower inside R₂₀₀; 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.

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₅₀₀ and R₅₀₀ values to visualize the differences in the shape of the profiles. The choice of the NFW results as scaling factors is motivated by the greater stability of the NFW masses compared to the other methods. Conversely, to study the radial shape of the profiles, the results obtained with the NP method are used as a benchmark of the local gravitational field because the method is fully data driven. In the left- and right-hand panels of the figure we show the reconstructed NFW and Einasto profiles, respectively, for the entire X-COP sample. We immediately notice that both parametric mass models trace the NP points closely, albeit with substantial scatter in the core and the outskirts. The additional degree of freedom afforded by the Einasto model allows it to trace the NP points more closely than the NFW, as there appears to be more diversity in the profile shapes than can be captured by the NFW model.

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

This visual impression is confirmed in a quantitative way in the left-hand panel of Fig. 4, where we show the sample median profiles and 68% percentiles for the various methods. The range of values in the cluster core is clearly broader for Einasto than NFW. The forward-fitting method also closely traces the NP profile, with the exception of the innermost regions, where the power-law shape of the gNFW pressure profile leads to excess pressure, and hence excess mass, compared to the NP data. In the left-hand panel of Fig. 4 we show the median NP-to-model ratio across the sample, as well as the dispersion of the model with respect to the reference NP value. The mass profiles estimated with all the methods are in remarkable agreement, with the median masses in the range [0.05 − 2]R₅₀₀ being within 3% of one another, even though the various methods make vastly different underlying assumptions. On average, we do not observe systematic deviations from an NFW shape, with the relative difference between NFW and NP points being less than 10% at all radii and consistent with the level of deviations found in the mock observations (see Sect. 4). Similarly, the Einasto profile accurately reproduces the average shape of the gravitational field in our systems. This study is therefore in agreement with that of E19, which concluded that the mass profiles of X-COP clusters are better represented by models that exhibit a rising density profile in their inner regions, whereas models that include large central cores (Burkert, isothermal sphere) are statistically disfavored. In all cases, the scatter is minimal in the radial range [0.2–1]R₅₀₀, 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₅₀₀.

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 mass-concentration relation for our sample. The small shaded ellipses indicate the 1σ contours of the fitted individual values for the 12 X-COP clusters individually, whereas the solid red and blue ellipses show the fitted sample mean and its 1σ and 2σ confidence regions, respectively. The curves show the predictions of N-body 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 mass-concentration relation. However, it is well suited to measure the mean and scatter of the relation around M₂₀₀ = 10¹⁵ M_⊙. To this end, we performed a Bayesian multivariate analysis of our sample, with the fractional log-normal scatter in mass and concentration as free parameters. Namely, we describe our sample as a set of values drawn from a distribution centered on the mean values for R₂₀₀ and c₂₀₀ with free intrinsic log-normal scatter on both axes. We set Gaussian priors on the two mean values, with the mean and standard deviation of the priors set to the mean and scatter of the sample values. A positive half-Cauchy prior with β = 0.5 is set on the fractional intrinsic scatter along both axes. Adopting instead a half-normal prior on the scatter has little impact on the results. To take the strong intrinsic correlation between mass and concentration in the NFW model into account, for each cluster we compute the covariance matrix between the two parameters from the output chains, and we set a global likelihood as the product of the Gaussian multivariate likelihoods of each object including the covariance matrix. We then sample the total likelihood using PyMC3.

The measured 1σ contours for all X-COP clusters in the M₂₀₀ − c₂₀₀ 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 R200 and c200 the mean values of R200 and NFW concentration in our sample and σR200 and σc the fractional log-normal intrinsic scatter of the points around the mean value.

We used the Python package Colossus (Diemer 2018) to compare our results with the predictions of various N-body simulation suites. The curves in the right-hand panel of Fig. 4 show the predictions of three different sets of simulations (Ludlow et al. 2016; Diemer & Joyce 2019; Ishiyama et al. 2021). The curves were calculated at the median redshift of the X-COP sample (z = 0.065) for a proper comparison with the data. We find a remarkably good agreement between our average values and the predictions of numerical simulations, with all three simulations considered here agreeing with our measurements within 1σ. On the other hand, numerical simulations also predict a large scatter in concentration at given mass, lower in more relaxed systems. Neto et al. (2007) find that a log-normal distribution represents the estimated concentrations in a given mass bin well, with a mean and a dispersion that decrease at higher masses, with the dispersion on the order of σ_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 X-COP data. Different halo collapse times, with higher concentration associated with a halo assembled earlier, can account for most of the measured scatter (e.g., Navarro et al. 1997; Ludlow et al. 2016).

All together, the shape parameters retrieved here agree well with the predictions of the ΛCDM framework, both for the average NFW concentration and the scatter of the c–M relation. A discussion of profile shape in the central regions of X-COP clusters as traced by the Einasto shape parameter is provided in Paper II, where we present our measurements of α 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,

(16)

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

The results of the mass reconstruction separating the baryonic and DM components are shown in Fig. 6. Here we focus on the Einasto reconstruction because of its greater flexibility; qualitatively similar results are obtained in the NFW case. The fitted NFW concentrations as well as the total masses obtained with this procedure are given in Table 2. In the left-hand panel of Fig. 6 we show the output differential mass density profiles for the BCG (green), gas mass (red), satellite galaxies (cyan) and DM (magenta). In the innermost regions (R < 0.02R₅₀₀) 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₅₀₀ as the dominant mass component. The NFW concentrations obtained for the fit to the DM component only are 10–20% smaller than what is obtained when fitting the total gravitational field (see Table 2), which results from the subtraction of the sharply peaked BCG profile in cluster cores. The baryon-driven adiabatic expansion of the DM halo under the influence of active galactic nucleus (AGN) feedback results in slightly lower concentrations (Duffy et al. 2010), in agreement with our findings. The DM profiles look remarkably similar in the radial range [0.2–0.6]R₅₀₀, 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₅₀₀ 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₅₀₀, such that the total mass of satellites within R₂₀₀ 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, dash-dotted cyan) and the gas mass (dashed red), whereas the fitted Einasto DM profiles are shown as the solid magenta curves. Right: total cumulative baryon fraction as a function of radius. The total baryon fraction and its uncertainty are shown as the solid curves and shaded areas, where the dashed and dotted lines show the relative contribution of gas and stars, respectively.

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 right-hand 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₅₀₀ and decreases steeply with radius. The integrated mass of the gas component dominates over that of stars beyond ∼0.1R₅₀₀, and inside R₅₀₀ 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₂₅₀₀, R₅₀₀, and R₂₀₀).

We observe an overall depletion of baryons in the 0.2 − 0.5R₅₀₀ range, which is compensated by the somewhat shallower slope of the gas density profile with respect to the DM. Inside R₅₀₀ 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₅₀₀ but slightly higher within R₂₀₀, 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₂₀₀ our baryon fractions exceed the cosmic value by 10–20%, potentially indicating a stronger impact of nonthermal pressure beyond R₅₀₀. 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 one-to-one relation when moving toward the external regions where the gravitational force is low. If this relation is found to be universal across all gravitationally bound halos, it could become the smoking gun for modified gravity theories such as MOND (Milgrom 1983) and its developments. The RAR was found to hold over a wide range in stellar mass (Lelli et al. 2017) with a low scatter of less than 0.13 dex and possibly even less than 0.06 dex (Li et al. 2018). Studies of the RAR in galaxy clusters found that the normalization of the observed acceleration is higher than what was found in rotationally supported galaxies (Chan & Del Popolo 2020; Tian et al. 2020; Pradyumna et al. 2021), although the uncertainties were too large to study in detail the shape of the RAR.

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

(17)

(18)

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

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,

(19)

with g_† = 1.2 × 10⁻¹⁰ m s⁻². At high acceleration (g > 10⁻¹⁰ m s⁻²), 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₅₀₀ 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⁻¹¹ m s⁻² we observe another change of regime, where both the observed and the baryonic components decrease with a slope that is almost parallel to the one-to-one relation. This corresponds to the regime where the baryonic component is dominated by the ICM mass and where the gas density profiles are slightly shallower than the DM profiles. Our systems eventually catch up with the McGaugh et al. (2016) RAR in the outermost regions, although the slope of the relation is widely different. The complex shape observed here indicates that the RAR in rotationally supported galaxy is not a universal property of gravity, since it fails to reproduce a gravitational field of similar strength in massive halos.

5.6. Modeling the galaxy cluster RAR

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

(20)

with g₁, g₂ the two characteristic acceleration scales describing the transition between the three regimes, and α₁, α₂ 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₂, transitioning to when g₂ < g_bar < g₁ and converging to 1 when g_bar ≫ g₁. Our empirical formula approaches the MOND formalism at low acceleration when α₁ = α₂ = 0.5 and g₁ is the classical MOND acceleration a₀ ≈ 10⁻¹⁰ m s⁻².

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 best-fit curve is shown in the left-hand panel of Fig. 8 where the functional form and its uncertainty are compared to the data. For the two acceleration scales we obtain m s⁻² and m s⁻², in qualitative agreement with the above discussion. The fitted slopes yield and . The scale g₁ is of the same order as the classical MOND acceleration, albeit slightly higher. However, the fitted slope α₁ is considerably larger than its corresponding MOND counterpart, which indicates that the transition away from the baryon-dominated regime is much faster than in disk galaxies. Below the second characteristic scale g₂ 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 MDM, 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 X-COP clusters within R₅₀₀ 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

(21)

We used the Einasto reconstructions to determine the mass profiles of M_DM, RAR. In the right-hand panel of Fig. 8 we show the enclosed mass profiles of the missing component assuming that Eq. (19) holds. The missing component dominates the mass profile beyond ∼0.05R₅₀₀, similar to the CDM case (Fig. 6). The required enclosed mass exceeds the gas mass and reaches a maximum around ∼0.5R₅₀₀ (∼600 − 700 kpc). This corresponds to the second transition scale g₂, 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₅₀₀ (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 X-COP hydrostatic masses are 10–15% lower than weak lensing estimates (Herbonnet et al. 2020). In Eckert et al. (2019) we used the integrated gas fraction as an anchor assuming that the true gas fraction can be robustly predicted by the ΛCDM framework to estimate the level of hydrostatic bias, which was found to be low (7% at R₅₀₀, 10% at R₂₀₀). 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₅₀₀ 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₅₀₀. The HSE assumption is likely to be a poorer approximation beyond this point.

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

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

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

6.2. Consistency with the ΛCDM framework

All of the results presented in this paper agree with the predictions of the ΛCDM framework and the bottom-up structure formation paradigm. The average mass profile of the X-COP clusters is well represented by an NFW model, with deviations of at most 10% over the entire radial range of interest ([0.01 − 2]R₅₀₀; see the left-hand panel of Fig. 4). The average NFW concentration of the sample matches perfectly the ΛCDM predictions (right-hand panel of Fig. 4), and N-body simulations also reproduce accurately the observed scatter in the mass-concentration relation at M₂₀₀ ∼ 10¹⁵ 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₂₀₀.

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₅₀₀) and outermost (R > R₅₀₀) regions covered in our study. The profiles estimated through our nonparametric technique appear to show a wider variety of shapes than can be described by the NFW parametric form (see Fig. 3). Conversely, the additional degree of freedom afforded by the Einasto parametric form allows the observed behavior to be reproduced more closely in each individual case. Some systems (e.g., A2255 and A644) show a larger amount of curvature in their mass profiles; the same systems were also found in E19 to be better described by models including a central core (Burkert, isothermal sphere). These systems exhibit a somewhat disturbed X-ray morphology as indicated by their large centroid shifts (see Paper II); thus, it is possible that the X-ray peak does not trace exactly the bottom of the potential well. This effect is likely to be present in A2255, which exhibits a very flat central gas density profile and a substantial offset between the X-ray peak and the position of the BCG (123 kpc, Rossetti et al. 2016). The apparent flat mass distribution in the core of this system may thus be explained by miscentering. On top of that, while the NFW profile is expected to provide an excellent description of the average mass profile extracted from CDM simulations, individual profiles show a scatter of ∼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 X-COP 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⁻¹⁰ m s⁻²), 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⁻¹¹ m s⁻², the relation steepens again and eventually catches up with the McGaugh et al. (2016) relation. The discrepancy cannot be attributed to potential systematics in our reconstruction such as hydrostatic bias (see Sect. 6.1) as they would most likely lead to an underestimation of the gravitational acceleration and increase the gap. Therefore, the behavior of the relation in X-COP clusters vastly differs from that in spiral galaxies, both qualitatively and quantitatively, indicating that the RAR is apparently not universal.

If the tight relation observed in spiral galaxies results from a modification of the theory of gravity at large scales (e.g., MOND), the observed behavior shows that an additional source of acceleration is required on top of the known baryonic components (see also Ettori et al. 2019). A possible way of alleviating the tensions with the MOND paradigm in the cluster regime is to invoke an additional matter component such as a light sterile neutrino with a mass in the range 1 − 10 eV (e.g., Sanders 2007; Angus et al. 2008; Nieuwenhuizen & Morandi 2013), which may be a natural candidate for a minimal extension of the standard model of particle physics. The missing component would dominate the matter content at intermediate radii where the discrepancy with MOND is large. Angus et al. (2010) argue that in the case of a sterile neutrino with a mass of 11 eV the gravitational field of galaxy clusters would be deep enough to retain the sterile neutrinos, while at galaxy scales the neutrinos would free stream from the halo. However, we note that invoking a substantial amount of hot DM affects the formation of structures in the Universe, which renders the expected halo mass function inconsistent with the observations (Angus et al. 2013). This scenario also invokes at the same time a modification of the gravitational law and a (less abundant) DM component, making it intrinsically more complex. Our analysis of the required missing mass in the standard RAR calibrated on disk galaxies indicates that this component should dominate the mass budget in the inner regions and its total mass should be comparable to the total baryonic mass (see Fig. 8). However, the mass distribution of the missing component should exhibit an unphysical decreasing trend beyond ∼0.5R₅₀₀, 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₂₀₀ 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₀ 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 X-COP clusters is quite complex and our model requires two characteristic acceleration scales, g₁ and g₂ (see Eq. (20)). The fitted value of g₁ is close to a₀, whereas the second characteristic scale g₂ is about an order of magnitude smaller. The g_bar − g_obs relation steepens again below g₂, which cannot be easily explained by a mass dependence of a₀.

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 X-COP galaxy clusters over the radial range [0.01 − 2]R₅₀₀ from a combination of deep XMM-Newton and Planck data. Our results can be summarized as follows:

• We introduced a novel framework for the reconstruction of hydrostatic mass profiles from X-ray and/or SZ data, which we distribute in the form of the public Python package hydromass. Hydrostatic equilibrium profiles can be reconstructed assuming a mass model, from a parametric pressure profile, or in a nonparametric way as a linear combination of log-normal functions. Our framework improves over our previous works (Ettori et al. 2019; the differences in the mass values are less than 10%; 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 XMM-Newton observations of a synthetic NFW cluster in HSE. Our mock observations include a wide range of instrumental effects, such as effective area, energy redistribution, vignetting, PSF convolution, and a sophisticated background model. The mock observations were analyzed using the same pipeline as the actual observations. Our code was found to reproduce the input profile very accurately, both for the 3D temperature profile and the hydrostatic mass (see Fig. 2).

• Applying our method to the 12 X-COP galaxy clusters, we find that the NFW and Einasto profiles both provide a very good representation of the average mass profiles over the entire radial range, in agreement with the predictions of the Λ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₂₀₀ = 3.69 ± 0.39 at M₂₀₀ ≈ 10¹⁵ M_⊙, in remarkable agreement with ΛCDM predictions. The intrinsic scatter of the mass-concentration relation was found to be , which again is consistent with the expectations of N-body simulations of structure formation.

• For a subset of five systems, we decomposed the gravitational field into its baryonic and DM components. We found that the stellar content of the BCG dominates the gravitational field in the innermost regions (< 0.02R₅₀₀), 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₅₀₀. The hot ICM dominates the baryonic mass beyond ∼0.1R₅₀₀ 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₂₀₀.

• 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⁻¹¹ m s⁻². We found that the RAR in galaxy clusters exhibits a complex shape with three distinct regimes delimited by two characteristic scales, which makes it difficult to reconcile with the MOND paradigm.

Application of the presented framework to larger samples (e.g., CHEX-MATE, The CHEX-MATE Collaboration 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 ASI-INAF Athena 2019-27-HH.0, “Attività di Studio per la comunità scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASI-INAF n. 2017-14-H.0), INAF mainstream project 1.05.01.86.10, and from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158). 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

Appendix A: Fitting results for individual clusters

Fig. A.1. Same as Fig. 1 but for A85.

Fig. A.2. Same as Fig. 1 but for A644.

Fig. A.3. Same as Fig. 1 but for A1644.

Fig. A.4. Same as Fig. 1 but for A2029.

Fig. A.5. Same as Fig. 1 but for A2142.

Fig. A.6. Same as Fig. 1 but for A2255.

Fig. A.7. Same as Fig. 1 but for A2319.

Fig. A.8. Same as Fig. 1 but for A3158.

Fig. A.9. Same as Fig. 1 but for A3266.

Fig. A.10. Same as Fig. 1 but for RXC1825.

Fig. A.11. Same as Fig. 1 but for Zw1215.

All Tables

All Figures

Fig. 1. Mass reconstruction results for A1795. Top left: electron density profile reconstructed with the multi-scale method used here (Eckert et al. 2020) compared to the L1 regularization results presented in Ghirardini et al. (2019) and publicly released by the X-COP team. Top right: nonparametric reconstruction of the 3D temperature profile (green) compared to the spectroscopic X-ray measurements (red) and the 3D temperature profile obtained by dividing the SZ pressure by the X-ray density (orange). See Sect. 3.4 for details. The blue curve shows the projected, spectroscopically weighted, and PSF convolved temperature profile, to be compared to the red data points. Middle left: electron pressure profile measured by Planck and XMM-Newton compared to the reconstructions obtained from the NFW and Einasto mass models (Sect. 3.2), the parametric forward approach (labeled “Forward”; see Sect. 3.3), and the NP reconstruction (Sect. 3.4). The dotted vertical lines show the location of R500 and R200 from the NFW fit. Middle right: output mass profiles obtained from the four individual reconstructions (NFW, Einasto, Forward, and NP). The data points show the NP results evaluated at the radius of the X-ray spectroscopic measurements (labeled “NP Spectral bins,” blue) and the SZ pressure data (“NP SZ bins,” red). Bottom left: same as previous but for the gas mass fraction. The gray shaded area shows the cosmic baryon fraction Ω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 XMM-Newton observations of a synthetic NFW cluster (see Sect. 4). The curves show the mean and standard deviation of the ten realizations for the various reconstruction techniques (NFW, green; Einasto, cyan; and NP, red). The bottom panels show the ratio of reconstructed profiles to true input profiles.

In the text

Fig. 3. Comparison between the mass profiles reconstructed using mass models (solid curves and shaded areas) and the nonparametric log-normal mixture method (data points). The results obtained for the NFW and Einasto profile are shown in the left and right panels, respectively. The data points show the posterior NP mass estimates at the radii of the X-ray spectral data (circles) and the SZ pressure data (asterisks). For clarity, all the profiles are scaled by the NFW values for R500 and M500. 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 mass-concentration relation for our sample. The small shaded ellipses indicate the 1σ contours of the fitted individual values for the 12 X-COP clusters individually, whereas the solid red and blue ellipses show the fitted sample mean and its 1σ and 2σ confidence regions, respectively. The curves show the predictions of N-body 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 R200 and c200 the mean values of R200 and NFW concentration in our sample and σR200 and σc the fractional log-normal 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, dash-dotted cyan) and the gas mass (dashed red), whereas the fitted Einasto DM profiles are shown as the solid magenta curves. Right: total cumulative baryon fraction as a function of radius. The total baryon fraction and its uncertainty are shown as the solid curves and shaded areas, where the dashed and dotted lines show the relative contribution of gas and stars, respectively.

In the text

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

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 MDM, 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