Highlight
Free Access
Issue
A&A
Volume 621, January 2019
Article Number A41
Number of page(s) 29
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/201833325
Published online 03 January 2019

© ESO 2019

1. Introduction

Galaxy clusters are the largest bound structures in the Universe. In the hierarchical structure formation scenario, they grow from primordial density fluctuations to form the massive structures we observe today (e.g., Kravtsov & Borgani 2012). Pristine gas falls into the dark matter potential and is progressively heated up to temperatures of 107 − 108 K, such that the majority of the baryons end up in the form of a fully ionized plasma, the intracluster medium (ICM), which produces X-ray emission through bremsstrahlung radiation and line emission. In addition, photons of the cosmic microwave background (CMB) crossing galaxy clusters are subject to inverse Compton scattering off the hot ICM electrons. This produces a spectral shift in the CMB signal, the Sunyaev-Zel’dovich (SZ) effect (Sunyaev & Zeldovich 1972), which is detectable at microwave wavelengths. X-ray emission and the SZ effect thus provide highly complementary diagnostics of the state of the ICM.

The thermodynamical properties of the ICM encode valuable information on the processes governing the formation and evolution of galaxy clusters. At first approximation, the state of the gas is determined by the properties of the gravitational potential and the merging history of the host halo alone. Gravitational collapse implies the existence of tight scaling laws between ICM properties and cluster mass (the self-similar model, Kaiser 1986; Bryan & Norman 1998) and non-radiative cosmological simulations predict that the scaled thermodynamical profiles of massive systems are nearly universal (e.g., Frenk et al. 1999; Borgani et al. 2005; Voit 2005). Thus, the distribution of the thermodynamical properties across the ICM is a powerful tool for probing the formation of galaxy clusters. Deviations from gravitational collapse predictions can be used as a way to quantify the impact of non-gravitational physics such as gas cooling and feedback from supernovae and active galactic nuclei (AGN; Tozzi & Norman 2001; Kay et al. 2002; Borgani et al. 2004; Kravtsov et al. 2005; Gaspari et al. 2014) and to probe hydrodynamical phenomena induced by structure formation such as shocks, turbulence and bulk motions (Dolag et al. 2005; Rasia et al. 2006; Vazza et al. 2009; Burns et al. 2010), cosmic rays (Pfrommer et al. 2007), and magnetic fields (Dolag et al. 1999). In recent years, hydrodynamical simulations of increasing complexity have attempted to model simultaneously all known non-gravitational effects and to determine their impact on the structural properties of galaxy clusters (e.g., Le Brun et al. 2014; Hahn et al. 2017; Planelles et al. 2017; Barnes et al. 2018; Gaspari et al. 2018). These simulations highlight the role played by AGN feedback in maintaining the cooling/heating balance and regulating the star formation efficiency.

In the past decade, X-ray observations have provided constraints of increasing quality on the distribution of gas density (Croston et al. 2006; Eckert et al. 2012), temperature (De Grandi & Molendi 2002; Vikhlinin et al. 2006; Pratt et al. 2007; Leccardi & Molendi 2008), pressure (Arnaud et al. 2010), and entropy (Cavagnolo et al. 2009; Pratt et al. 2010) throughout the ICM. However, these observations sampled a radial range limited to less than half of the virial radius, and thus the majority of the cluster volume, and the entirety of the accretion region, remained unexplored. Thanks to its low instrumental background, the Suzaku satellite allowed us to extend the accessible radial range to the virial radius in a dozen clusters (e.g., Reiprich et al. 2013; Akamatsu et al. 2011; Simionescu et al. 2011, 2017; Walker et al. 2012a,b; Urban et al. 2014). Somewhat surprisingly, most studies find steeply decreasing temperature profiles and a flattening of the entropy at R500 and beyond, at odds with the predictions of gravitational collapse. Possible explanations for these behaviors include biases in X-ray measurements caused by gas clumping (Nagai & Lau 2011; Simionescu et al. 2011; Vazza et al. 2013; Roncarelli et al. 2013), non-equilibrium between electrons and ions (Hoshino et al. 2010; Avestruz et al. 2015), breakdown of hydrostatic equilibrium caused by turbulence and bulk motions (Lapi et al. 2010; Okabe et al. 2014; Avestruz et al. 2015; Khatri & Gaspari 2016), or as-yet-unknown systematics in the subtraction of the Suzaku background. Clearly, accurate observational data out to the virial radius are required to understand the origin of these deviations.

Since the first X-ray observations of galaxy clusters it has become evident that clusters can be divided into two categories differing only in their central properties (Jones & Forman 1984). The origin of these two populations is still unclear; however, several recent works (Rossetti et al. 2017; Lovisari et al. 2017; Andrade-Santos et al. 2017) point out that there is a clear difference between X-ray and SZ selected samples, also generating possible biases in the measured thermodynamic quantities.

In this paper, we present the universal thermodynamical properties of the galaxy clusters in the XMM Cluster Outskirts Project (X-COP), a sample of 12 massive, local galaxy clusters (z <  0.1) with deep XMM-Newton and Planck data covering the entire azimuth out to the virial radius. Unlike previous studies that exclusively utilized the X-ray signal, we take advantage of the high signal-to-noise ratio (S/N) of our clusters in the Planck survey (Planck Collaboration XXIX 2014), and combine X-ray and SZ data to increase the precision of our measurements while keeping a good control of systematic errors. This method was applied to reconstruct the thermodynamical properties of a few clusters (Basu et al. 2010; Eckert et al. 2013; Adam et al. 2015; Ruppin et al. 2017), and we demonstrated the ability of XMM-Newton and Planck to measure accurately the state of the gas out to the virial radius in two pilot studies (Tchernin et al. 2016; Ghirardini et al. 2018).

This paper is the first in a series presenting the results of the analysis of the full X-COP sample. Here we consider the global properties of the X-COP sample and present universal functional forms describing the structural properties of the ICM over more than two decades in radius ([0.01 − 2]R500). We also determine the slope and intrinsic scatters of the various quantities as a function of radius. In Ettori et al. (2019) we present our reconstruction of the mass profiles of our clusters, whereas in Eckert et al. (2019) we estimate the amount of non-thermal pressure support affecting our sample.

The paper is organized as follows. In Sect. 2 we describe the available dataset and the analysis procedures. In Sect. 3 we present our results on the universal thermodynamic profiles, slopes, and intrinsic scatter. Our findings are discussed in Sect. 4. Throughout the paper, we assume a flat Λ cold dark matter (CDM) cosmology with Ωm = 0.3, ΩΛ = 0.7, and H0 = 70 km s−1 Mpc−1. All our fitting is performed using the Bayesian nested sampling algorithm MultiNest (Feroz et al. 2009) unless otherwise stated.

2. Dataset and analysis procedures

2.1. The X-COP project

X-COP (Eckert et al. 2017) is a very large program (VLP) on XMM-Newton whose aim is to advance our understanding of the virialization region of galaxy clusters. The strategy adopted for the project is to target the most significant Sunyaev-Zel’dovich (SZ) sources in the Planck survey in order to combine X-ray and SZ information out to the virial radius. This strategy was already applied to the cases of A2142 (Tchernin et al. 2016) and A2319 (Ghirardini et al. 2018) and it was shown to allow the detection and characterization of the ICM out to R200 and even beyond. The sample was selected from the first Planck SZ catalogue (PSZ1; Planck Collaboration XXIX 2014) according to the following criteria:

  • PSZ1 S/N > 12. We select only the sources with the strongest SZ signal to ensure detection at high significance beyond R500;

  • 0.04 < z < 0.1. This criterion allows us to select objects for which the region of interest can be covered with five XMM-Newton pointings (one central and four offset);

  • θ5001 > 10 arcmin. Given the 7 arcmin Planck beam, this criterion ensures that the clusters are well resolved by Planck;

  • NH2 <  1021 cm2. Since we use a soft energy band ([0.7–1.2] keV, see below) to extract surface brightness and gas density profiles, this cut avoids objects for which the X-ray flux is strongly suppressed below ∼1 keV.

These criteria lead to a set of 16 clusters from the PSZ1 catalogue. Of these systems, we excluded four: two (A754 and A3667) because they exhibit strongly disturbed morphologies, rendering an azimuthally averaged analysis difficult; one (A2256) because of its bad visibility for XMM-Newton; and one (A3827) because its apparent size θ500 is very close to our cut and so may not be properly resolved by Planck. The final list of clusters is given in Table 1, together with some basic properties of these systems and classification into the cool-core and non-cool-core classes, using the value of K0 from Cavagnolo et al. (2009) as discerning value (see also Sect. 3.7 for further details). Like other Planck-selected samples (Rossetti et al. 2017; Andrade-Santos et al. 2017; Lovisari et al. 2017), the X-COP sample is dominated by non-cool-core systems (8 out of 12).

Table 1.

Basic properties of the X-COP sample.

The nominal XMM-Newton exposure time was set to 25 ks per pointing, which allows us to reach a limiting surface brightness of 3 × 10−16 ergs cm−2 s−1 arcmin−2 in the [0.5–2.0] keV band. Including central and archival pointings, the total observing time is about 2 Ms. The observation log, observation IDs, and observing time, after applying flare filtering, are given in Table F.1.

2.2. XMM-Newton data analysis

Here we describe in detail the data analysis pipeline that we set up for the analysis of the XMM-Newton data. A flow chart describing the main steps of the analysis is presented in Fig. 1.

thumbnail Fig. 1.

Flow chart of the XMM-Newton data analysis pipeline. The steps of the analysis are shown in red, the main intermediate and final products are described in the black boxes, and the procedures are shown in green italics.

2.2.1. Data reduction

We reduced all the data using XMMSAS v13.5 and the Extended Source Analysis Software (ESAS) data reduction scheme (Snowden et al. 2008). To perform basic data reduction, we used the emchain and epchain pipelines to extract calibrated event files from the observations, and we reran epchain in out-of-time mode to create event files for pn out-of-time events. To filter out time periods affected by soft proton flares we ran the mos-filter and pn-filter executables, which extract the light curve of each observation in the hard band, and applied a sigma-clipping technique to exclude time intervals with enhanced background. We used the unexposed corners of the MOS detectors to monitor the particle background level during each observation, and measured the count rates in the high-energy band ([7.5–11.8] keV) of the MOS from the regions located inside and outside the field of view (FOV) of the telescope (hereafter IN and OUT). The comparison between IN and OUT count rates was then used to estimate the contamination of residual soft protons to the spectrum (De Luca & Molendi 2004; Leccardi & Molendi 2008; Salvetti et al. 2017).

2.2.2. Image extraction and preparation

We extracted photon count images from the three EPIC detectors in the [0.7–1.2] keV band. This narrow band maximizes the ratio between source and background emission and is thus best suited to minimize the systematics in the subtraction of the EPIC background (Ettori et al. 2010). We used eexpmap to compute exposure maps taking vignetting effects into account for all three detectors independently. To create total EPIC images, we summed the count maps of the three detectors and repeated the same operation with the exposure maps, multiplying the pn exposure maps by a factor of 3.44 representing the ratio of pn to MOS effective areas in our band of interest. The high-energy particle background are modeled and subtracted by simple subtraction of the rescaled background images, taken from the unexposed corners of the detectors.

Even after cleaning the light curves from soft-proton flares, it is known that a fraction of residual soft proton contamination remains within the datasets. This component can introduce systematics at the level of ∼20% in the subtraction of the EPIC background in our band of choice (Tchernin et al. 2016). To model the contribution of residual soft protons, we follow the method outlined in Ghirardini et al. (2018), which was calibrated using a large set of ∼500 blank-sky pointings. Namely, we measure the high-energy ([7.5–11.85] keV) MOS count rates in the exposed (IN) and unexposed (OUT) parts of the FOV, and we use the difference between IN and OUT count rates as an indicator of the contamination of each observation by residual soft protons (see Salvetti et al. 2017, for a detailed overview of this approach). We use our large blank field dataset to calibrate an empirical relation between the IN–OUT indicator and the required intensity of the soft proton component (see Appendix A of Ghirardini et al. 2018), and we use this relation to create a 2D soft proton model. This procedure was shown to bring the systematics in the subtraction of the EPIC background to an accuracy better than 5%.

For each cluster, we combine the resulting EPIC count maps in the [0.7–1.2] keV from each individual observation (central or offset) to create a mosaic image. We apply the same procedure to the combined EPIC exposure maps and to the background maps, summing the non-X-ray background components (quiescent particle background and residual soft protons). We then obtain mosaicked photon maps, exposure maps, and non-X-ray background maps for each source. In Fig. 2 we show the resulting background-subtracted and exposure-corrected mosaics, which we adaptively smoothed using the asmooth code (Ebeling et al. 2006).

thumbnail Fig. 2.

Adaptively smoothed and exposure corrected XMM-Newton mosaic images in the [0.7–1.2] keV energy band for all X-COP clusters. The superimposed white contours represent the Planck SZ S/N maps between 70 and 857 GHz. The contour levels correspond to 1, 3, 5, 7, 10, 15, 20, 30, 40, and 50σ. The spatial scale is given by the red thick lines, which have a common length of 15 arcmin. The color scale is given below the images in units of cts s−1 pixel−1, and is the same for all clusters.

2.2.3. Point source subtraction

To detect point sources present within the field, we extract photon count maps from the three detectors in a soft ([0.5–2] keV) and a hard ([2–7] keV) band, and we use the XMMSAS tool ewavelet with wavelet scales in the range 1–8 pixels and S/N threshold of 5.0. We then cross-match the soft and hard band detections between the multiple (central and offset) observations of each cluster to create a global point source list per cluster. Since the vignetting and the point spread function of the XMM-Newton telescopes depend on off-axis angle, the sensitivity threshold for source detection depends on the position of a source on the detector. At a fixed observing time, XMM-Newton thus detects point sources down to lower fluxes near the aim point than close to the edge of the FOV, and the fraction of the cosmic X-ray background (CXB) that is resolved by the instrument is spatially dependent. To correct for this effect, we draw the distribution of measured count rates from the detected sources and we determine the count rate at which the distribution peaks. Since the logN–logS of distant sources contributing to the CXB is a monotonically decreasing function (e.g., Moretti et al. 2003) the peak in the count rate distribution of our observation roughly corresponds to the threshold down to which our source detection is complete. We then excise only the sources with a measured count rate greater than our threshold and leave the fainter sources to enforce a constant flux threshold across the FOV and avoid biasing local measurements of the CXB intensity.

2.2.4. Surface brightness profiles

An important complication with the analysis of X-ray data of cluster outskirts lies in the presence of accreting structures and inhomogeneities in the gas distribution, which are expected to contribute substantially to the measured X-ray flux beyond ∼R500 (e.g., Nagai & Lau 2011; Vazza et al. 2013; Roncarelli et al. 2013). Since the X-ray emissivity is proportional to the square of the gas density, overdense regions contribute predominantly to the measured X-ray flux and can bias high the recovered gas density. To avoid this problem we apply the azimuthal median method outlined in Eckert et al. (2015). Numerical simulations show that the median surface brightness in concentric annuli is robust against the presence of outliers in the gas density distribution (Zhuravleva et al. 2013). We thus construct background-subtracted and exposure-corrected surface brightness maps using a Voronoi tessellation technique (Cappellari & Copin 2003) with a target number of 20 counts per bin.

The intensity of the sky background is determined by averaging the surface brightness distribution in the regions with R >  2R500, where we assume the cluster emission to be negligible. A systematic uncertainty of 5% of the background level (see Sect. 2.2.2) was added in quadrature to the error budget of the surface brightness profiles. This procedure was applied to all clusters except A3266, for which the current mosaic does not extend out to 2R500. In this case, we estimate the sky background intensity from the ROSAT all-sky survey background tool3 and included a systematic error of 30% in quadrature to the full error budget.

To compute the surface brightness profiles, we draw the surface brightness distribution from the Voronoi-binned images in concentric radial bins starting from the X-ray peak and then choosing the annuli such that the emissivity in each bin is almost constant. The errors on the azimuthal median are estimated from 104 bootstrap resampling of the pixel distribution. Circular regions of 30 arcsec radius are excised around the positions of point sources selected through the procedure described in Sect. 2.2.3, corresponding to an encircled energy fraction of 90% of the point source flux.

2.3. Deprojection and gas density profiles

To extract gas density profiles, we take advantage of the fact that the X-ray surface brightness in our energy band of choice is proportional to the squared gas density integrated along the line of sight. To convert surface brightness profiles into emission measure, we describe the emissivity of the source with a thin-plasma model absorbed by the Galactic NH and folded through the on-axis EPIC/MOS effective area. This approach allows us to calculate the conversion between the observed count rate in MOS units and the normalization of the APEC model, which is related to the plasma emission measure as

(1)

where dA is the angular diameter distance of the source, and ne and nH are the electron and ion number densities in units of cm−3, with ne = 1.17nH in a fully ionized plasma (Anders & Grevesse 1989). Since we are using a soft energy band, the conversion between count rate and emission measure shows little dependence on the temperature as long as the temperature exceeds ∼1.5 keV, which is the case in all X-COP systems. The resulting emission measure profiles can then be deprojected under the assumption of spherical symmetry by computing the projected volumes V of each spherical shell onto each 2D annulus. To recover the 3D emissivity and density profiles from the projected data, we apply two different deprojection methods that we briefly outline here.

  • L1 regularization. This method builds on the non-parametric regularization approaches developed by Croston et al. (2006) and Ameglio et al. (2007), introducing a penalty term on the modulus of the second derivative of the 3D density profile to kill spurious small-scale fluctuations introduced by the random nature of the data (Diaz-Rodriguez et al. 2017). Given an observed 2D emission measure profile EM = (EM1EMn) and corresponding uncertainties σEM = (σEM, 1σEM, n), the values of the 3D emissivity profile ϵ = (ϵ1ϵn) are obtained by maximizing the following likelihood function

    (2)

    where Vi, j is the geometrical matrix volume of the jth shell intercepted by the ith annulus, # is the symbol for matrix product, and the sum is performed along all the annuli. Moreover the second derivative of the emissivity is computed as a numerical derivative of the ϵ(r) vector. The parameter λ controls the degree of regularization of the profile. To maximize the likelihood function described in Eq. (2), we use the Markov chain Monte Carlo (MCMC) tool emcee (Foreman-Mackey et al. 2013), leaving the value of the 3D density profile at each radius as a free parameter and setting a logarithmic prior (i.e., uniform prior in logarithmic space) on each parameter value to enforce positivity of the resulting profile. The value of the parameter λ is chosen such that the log-likelihood is about 1 per data point, to allow for typical statistical deviations of 1σ. We note that λ = 0 is equivalent to using the onion-peeling technique directly (see Kriss et al. 1983; Ettori et al. 2002, 2010).

  • Multiscale fitting. This method follows the technique developed in Eckert et al. (2016), whereby the projected emission measure profile is decomposed into a sum of analytical multiscale functions which can be individually deprojected. Following Eckert et al. (2016) we write the observed 2D profile as a sum of N King functions with fixed core radii and normalizations and slopes allowed to vary while fitting, choosing N = Npoints/4; i.e., one model component is added for every set of four data points, fixing a core radius to the mean radial value of these four data points. Since the projection kernel is linear, each King function can be individually deprojected and the 3D profile can be analytically reconstructed from the fit to the projected data. As above, we use emcee to optimize for the parameters and reconstruct the error envelope around the best fitting curve.

In the top left panel of Fig. 3 we compare the density profiles reconstructed with the two methods and find a remarkable agreement between them, with an average scatter < 5% at each radius. By construction, the profiles reconstructed with the L1 regularization method shows more pronounced features as the method imposes fewer constraints on the shape of the profile, whereas the profiles obtained with the multiscale method are smoother. Thus, we adopt the results of the L1 regularization when attempting to determine the exact shape of our profiles, whereas the multiscale technique is preferred when reconstructing hydrostatic mass profiles to provide better control over the gradient.

thumbnail Fig. 3.

Top left panel: density profiles for all X-COP clusters obtained with two different deprojection methods: L1 regularization (data points) and multiscale fitting (solid lines). The magenta shaded area represents the scatter of the median profile in Eckert et al. (2012). Bottom part: ratio of the two methods for each individual system. Top right panel: joint fit to all the density profiles using piecewise power laws in several radial ranges (color-coded). The best fits and 1σ error envelope are shown by the solid lines, while the dashed lines represent the intrinsic scatter. Bottom left panel: joint fit to the density profiles using the functional form introduced by Vikhlinin et al. (2006), in red, with the shaded area indicating the 1σ error envelope around the best fit. The dashed lines represent the intrinsic scatter in the functional form as a function of radius. Bottom right panel: slope of the density profiles as a function of radius. The green data points show the results of the piecewise power law fits, whereas the red curve indicates the fit to the entire radial range using the Vikhlinin et al. (2006) functional form. In all panels, the vertical dashed and dotted lines represent the location of R500 and R200, respectively.

2.4. Spectral analysis

2.4.1. Spectral extraction

We extract spectra in concentric annuli around the X-ray peak covering approximately the radial range [0 − 1]R500, removing the point sources which contaminates the spectra (see Sect. 2.2.3). In the spectral analysis, differently from the imaging case in Sect. 2.2.2, we compute models of the high-energy particle background using the ESAS tools mos-spectra and pn-spectra in imaging mode. To this end, we select the filter-wheel-closed observations recorded at the nearest possible time to the observation, and we use the spectra of the unexposed corners of the detectors to rescale the filter-wheel-closed observations. This procedure is performed individually for all CCDs, and the CCDs operating in anomalous mode are discarded. We then extract an image from the rescaled filter-wheel-closed data in the [0.7–1.2] keV band using the mos-back and pn-back executables, which we use as our model for the high-energy particle background. In the case of the pn, we repeat the operation with the out-of-time event files and create a model for the intensity and spatial distribution of out-of-time events.

We selected the binning such that the width of the bins increases exponentially, but choosing a minimum width of 0.5 arcmin for the innermost bins such that the instrumental PSF does not contribute much to the photon in each bin. We group the output spectra with a minimum of five counts per bin to ensure stable fitting results, and discard the data below 0.5 keV where the EPIC calibration is uncertain. We then use XSPEC v12.9 and the C-statistic (Cash 1979) to fit the spectra and determine the plasma parameters (see Sect. 2.4).

2.4.2. Spectral modeling

To extract spectral diagnostics from the observed spectra (see Sect. 2.4.1), we proceeded using a full spectral modeling approach following the method described in detail in Eckert et al. (2014). Here we describe our approach to model all the individual background components and the source spectra.

  • High-energy particle background. We use the rescaled filter-wheel-closed spectra to determine the intensity and spectral shape of the particle background. We fit the filter-wheel-closed spectra using a diagonal response matrix and a phenomenological model including a broken power law and several Gaussians to reproduce the shape of the continuum and fluorescence lines. We then apply the fitted model to the source spectrum, leaving the normalization free to vary within ±10% to account for possible systematics in the scaling of the filter-wheel-closed data.

  • Sky background. We model the X-ray background and foreground emission as the sum of three components: (i) an absorbed power law with a photon index fixed to 1.46 to describe the residual CXB (De Luca & Molendi 2004); (ii) an absorbed APEC thermal plasma model with a temperature allowed to vary in the range [0.15–0.6] keV to model the Galactic halo emission (McCammon et al. 2002); and (iii) an unabsorbed APEC model with a temperature fixed to 0.11 keV to represent the local hot bubble. The Galactic hydrogen column density NH was fixed to the LAB value (Kalberla et al. 2005). Similarly to what was done for the imaging case, the parameters of the sky emission model are fitted to the spectra of background regions located at R >  2R500 from the cluster core. Again the exception to this procedure is A3266, for which we use the ROSAT all-sky survey background tool to determine the sky background parameters. The best fit model is then applied to the source spectra, rescaling the intensity of the components according to the area of each region.

  • Residual soft protons. In cases where our IN–OUT indicator of soft proton contamination is found to be high (IN–OUT > 0.1 counts s−1), we include an additional model component to the particle background model to take soft protons into account. We model the soft proton component as a broken power law with fixed spectral shape (slopes 0.4 and 0.8 and break energy 5 keV, Leccardi & Molendi 2008) and leave the normalization of this component free to vary in the overall fitting procedure.

  • Source. We model the source emission in each annulus as an absorbed single-temperature APEC model with temperature, emission measure, and metal abundance free to vary. In cases where multiple observations were available for the same regions, we fit all the available spectra jointly, tying the source parameters between the different spectra. The solar abundance table is set to Anders & Grevesse (1989). Since our objects are nearby and extended on scales much larger than the XMM-Newton PSF, we neglect the potential cross-talk between the various annuli.

All the spectra were fitted in the energy range [0.5–12] keV using XSPEC v12.9, ATOMDB v3.0.7, and the C-statistic (Cash 1979). When several observations of the same region were available, we extracted the spectra from each individual pointing and fit them jointly. We ignored the energy ranges [1.2–1.9] keV (MOS) and [1.2–1.7] keV, [7.0–9.2] keV (pn) where bright and time-variable fluorescence lines were present. We then constructed projected gas temperature profiles from the best fit results. We also deprojected our 2D temperature profiles using the projection matrix V and the emissivity in each annulus as weights, adopting the spectroscopic-like temperature scaling of Mazzotta et al. (2004).

2.5. Planck data analysis

The signal was recovered from the Planck survey (Tauber et al. 2010; Planck Collaboration I 2016) making use of the six frequency maps provided by the High Frequency Instrument (HFI; Lamarre et al. 2010; Planck HFI Core Team 2011). They were combined with a modified internal linear combination algorithm (MILCA; Hurier et al. 2013) to produce Comptonization parameter maps, i.e., y-maps tracing the intensity of the SZ effect. The maps used for the X-COP project are provided with a resolution of 7 arcmin FWHM.

The intra-cluster gas thermal pressure integrated is recovered from the y parameter (Sunyaev & Zeldovich 1972) as

(3)

where the integral is computed along the line of sight, θ is the angular distance to the cluster core, σT is the Thomson cross-section, me is the mass of the electron, and c is the speed of light.

As illustrated for the cases of A2142 (Tchernin et al. 2016) and A2319 (Ghirardini et al. 2018), we followed for the whole X-COP sample the methodology presented and used in Planck Collaboration Int. V (2013). Assuming azimuthal symmetry of the cluster, we computed the y radial profile for each cluster over a regular grid scaled in units of R500 with radial bins of size Δθ/θ500 = 0.2 out to 10 × R500. The local background is assumed to be flat and constant and is computed from the area beyond 5 × θ500. A covariance matrix is computed for each profile to account for the correlation between points due to the profile binning, and intrinsic noise correlation introduced from the y map construction.

We corrected from the Planck beam redistribution through a real space deconvolution of the instrument PSF. Further assuming the spherical symmetry of the source, we reconstructed each pressure profile through a geometrical deprojection. The two steps follow the method initially presented in Croston et al. (2006). The associated covariance matrix for the pressure profile is obtained via a Monte Carlo procedure by randomizing over the initial y profile covariance matrix. In the following we ignore the innermost three Planck data points because of the difficulty deconvolving from the large Planck beam.

As an alternative method, we also extracted Planck pressure profiles using the forward-modeling approach of Bourdin et al. (2017). In this case, a spectral model for the relevant components is constructed (CMB, dust, synchrotron, and thermal SZ). The model is folded through the Planck response and fitted to the multi-frequency data points (see Bourdin et al. 2017, for details). In Appendix B we compare the results obtained with this approach to the results of the MILCA component separation method and show their consistency. For the remainder of the paper, we use the MILCA pressure profiles as our default choice.

2.6. Mass estimates

To estimate scale radii and self-similar scaling quantities, we use the high-precision hydrostatic mass reconstructions presented in Ettori et al. (2019). The mass models were obtained by combining X-ray and SZ information for each individual system and solving the hydrostatic equilibrium equation. For the present work, we adopt as our reference mass model the backward NFW results, which were obtained by assuming that the mass profile follows a Navarro-Frenk-White shape (Navarro et al. 1996) with scale radius and concentration c200 as free parameters. This method was shown to provide the best representation of the data (Ettori et al. 2019) and at R500 it matches the results obtained without assuming a functional form for the mass profile with an accuracy of ∼5%. Comparison of our mass reconstruction with weak lensing and SZ estimates (Ettori et al. 2019) and convergence toward the expected universal gas fraction (Ettori et al. 2019) show that our masses and the corresponding values of R500 are accurate at the 10% and 3% level, respectively. For more details on the reconstruction of hydrostatic masses and estimates of systematic uncertainties we refer to Ettori et al. (2019) and Eckert et al. (2019).

3. Thermodynamic properties

In the following section we describe how we derived the universal profiles, slopes, and intrinsic scatter of all our thermodynamic variables. We then present our main results and provide best fitting functional forms describing the X-COP cluster population.

3.1. Fitting procedure

We adopt two different approaches to fit the thermodynamic properties:

  • Piecewise power law fits. In this case, we split our data into several radial ranges as a fraction of R500 and we approximate the global behavior of the population in each range as a power law with free log-normal intrinsic scatter

    (4)

    with x = R/R500; Q/Q500 the rescaled thermodynamic quantity at an overdensity of 500; and A, B, and σint the normalization, slope, and intrinsic scatter in the radial range of interest. The values of Q500 are computed adopting the virial theorem as in Voit et al. (2005), and are shown in Eqs. (8), (10), and (12). The fitting procedure thus has three free parameters (A, B, and σint) in each of the chosen radial ranges. We note that this procedure provides model independent measurements of the slope and intrinsic scatter at different radii.

  • Global functional forms. In this case, we describe the radial dependence of the thermodynamic quantity of interest throughout the entire radial range with a parametric functional form found in the literature

    (5)

    with Q/Q500 the rescaled thermodynamic quantity at an overdensity of 500, as in Eq. (4), and f(x) the chosen functional form for the thermodynamic quantity Q. In this case, since we model the whole radial range covered by our measurements, we allow the intrinsic scatter to vary with radius following a quadratic functional form to model the radial dependence of the intrinsic scatter

    (6)

    with σ1 the width of the log-parabola, and x0 and σ0 respectively the location and the intercept of the minimum of the log-parabola. A total of n + 3 parameters are allowed to vary during the fitting procedure, with n the number of parameters of the adopted functional form f(x). Optimizing jointly for the parameters of the intrinsic scatter profile allows us to determine the shape of σint(x).

In Figs. 36 we show our rescaled thermodynamic quantities for the gas density, pressure, temperature, and entropy. The best fits with piecewise power laws and using functional forms are presented as well, together with the slopes for the parametric and non-parametric cases.

thumbnail Fig. 4.

Same as Fig. 3 for the pressure profiles rescaled by the self-similar quantity P500 (Eq. (8)). The squares indicate data points obtained from the deprojection of the SZ signal, while the filled circles are computed by combining the X-ray gas density profiles with the spectroscopic temperature. The solid red curve in the bottom panels shows the joint best fit to the data with the generalized NFW functional form (Nagai et al. 2007, see Eq. (9)). In all these plots the dotted and dashed-dotted lines represent the result of Planck Collaboration Int. V (2013) and Arnaud et al. (2010) respectively. The shadow areas represent the envelope obtained by Planck Collaboration Int. V (2013) and the Early release SZ sample (XMM-ESZ; Planck Collaboration VIII 2011).

thumbnail Fig. 5.

Same as Fig. 3 for the projected temperature profiles rescaled by the self-similar quantity T500 (Eq. (10)). The filled circles show the measurements of the X-ray spectroscopic temperature (see Sect. 2.4.2 ), whereas the filled squares indicate the data points obtained by combining the SZ pressure with the gas density, projected along the line of sight assuming the spectroscopic-like scaling of Mazzotta et al. (2004). The solid red curves in the bottom panels show the best fit to the joint dataset with the functional form introduced by Vikhlinin et al. (2006; see Eq. (11)).

thumbnail Fig. 6.

Same as Fig. 3 for the entropy profiles rescaled by the self-similar quantity K500 (Eq. (12)). The filled circles show the measurements obtained from the X-ray spectroscopic temperature as , whereas the filled squares indicate the data points obtained by combining the SZ pressure with the gas density as . The solid red curves in the bottom panels indicate the best fit to the entire population with the functional form presented in Eq. (13), whereas the solid black line shows the prediction of pure gravitational collapse (Voit et al. 2005).

3.2. Density

We rescaled our density profiles by the self similar quantities, E2(z)=Ωm(1 + z)3 + ΩΛ and R500 for density and radius, respectively. In Fig. 3 we compare our scaled gas density profiles with the “universal” density profiles from Eckert et al. (2012) from a sample of 31 clusters with available ROSAT/PSPC pointed data. We observe excellent agreement with their results; A2319 is the only exception that deviates at large radii, as shown in Ghirardini et al. (2018), because of its large non-thermal pressure support.

In the top right panel of Fig. 3 we show our density profiles fitted with piecewise power laws in several radial ranges, we show the best fitting parameters in Table 2. We parametrize the behavior of our density profiles over the whole radial range by adopting the Vikhlinin et al. (2006) functional form

(7)

Table 2.

Results of the piecewise power law fits (normalizations, slopes, and intrinsic scatter; see Eq. (4)) for the various thermodynamic quantities in several radial ranges, marked by the inner and outer rescaled radii xin and xout.

with x = R/R500 and γ = 3 fixed. The form thus has six free parameters (n0, rc, α, β, rs, and ϵ) and is able to reproduce both the core and the outer parts of the density profile. We apply flat priors in logarithmic space to n0, rc, and rs and flat priors in linear space to the remaining parameters, constraining ϵ <  5 (we specify the priors adopted in Table 3). We show the posterior distributions of the parameters of this functional form in Fig. E.1 and the covariance between them.

Table 3.

Best fit parameters of the functional forms describing the universal thermodynamic quantities.

The resulting profile is consistent at all radii with the universal envelope computed by Eckert et al. (2012). We can see in all panels how the profiles become progressively less scattered going toward the outskirts. While the core is affected by a large scatter likely caused by cooling, AGN feedback and different merger states, the profiles show a high degree of self-similarity in the radial range [0.3 − 1]R500. Then in the outskirts the scatter increases again, likely caused by different accretion rates from one system to another. We note from the plot in bottom right panel in Fig. 3 that the slope of the density profiles steepens steadily with radius. The slopes computed from the piecewise power law fits and from the global fit with Eq. (7) agree within 1σ at all points. Again, this result agrees with the findings of Eckert et al. (2012) and Morandi et al. (2015), but is at variance with the relatively flat slopes reported in several clusters observed by Suzaku. For instance, several papers report density slopes flatter than −2.0, for example −1.7 in the outskirts of the Perseus cluster (Urban et al. 2014) or even as low as −1.2 in A1689 (Kawaharada et al. 2010 and Virgo Simionescu et al. 2017). These measurements are clearly in tension with the slope of −2.5 at R200 measured here for the X-COP cluster population. It must be noted, however, that thanks to the azimuthal median technique, our gas density profiles are essentially free of the clumping effect, whereas the Suzaku studies could not properly excise overdense regions because of the lower resolution of the instrument and/or observations performed along narrow arms.

3.3. Pressure

Pressure in galaxy clusters is usually the smoothest thermodynamic quantity along the azimuth, if the cluster is not affected by an ongoing merger.

We recover the gas pressure both through the combination of X-ray gas density and spectral temperature and through the direct deprojection of the SZ effect. In the former case, we deproject the spectral X-ray temperature (e.g., Mazzotta et al. 2004, see Sect. 2.4.2) and combine the deprojected temperature with the gas density interpolated on the same grid to infer the pressure PX = kBTX × ne. In the latter we recover the pressure directly from the Planck data by deprojecting the measured y profiles (see Sect. 2.5) from which we exclude the first three points from the analysis. We thus combine the higher resolution and precision of XMM-Newton in the inner region with the high quality of the Planck data at R500 and beyond, which allows us to constrain the shape and intrinsic scatter of the pressure profiles in the radial range [0.01 − 2.5]R500. We rescale our pressure profiles by the self-similar quantities at an overdensity of 500,

(8)

where fb is the Universal gas fraction, which we take to be Ωbm = 0.16 (Planck Collaboration XIII 2016), rounding the number to 2 significant figures), μ and μe are the mean molecular weight per particle and mean molecular weight per electron for which we adopt the values measured by Anders & Grevesse (1989). In Fig. 4 we show the scaled pressure profiles of our 12 objects obtained through X-ray and SZ measurements of the ICM. We note that our profiles agree with the results obtained by the Planck Collaboration for a sample of 62 clusters (Planck Collaboration Int. V 2013), falling well within the two envelopes, with the exception, as in the case of the density, of A2319 (see the discussion in Ghirardini et al. 2018).

Similarly to the density, we fitted the profiles using piecewise power laws in several radial ranges, also obtaining in this case a scatter that decreases with radius out to R500 and then (as for the density) that increases in the outskirts; these profiles become progressively steeper with radius, see Table 2. Our profiles in the outskirts are compatible with the results of the Planck Collaboration, both for the central value and the slope. We also fitted our data using the generalized NFW functional form introduced by Nagai et al. (2007)

(9)

where x = R/R500, with five free parameters P0, c500, and three slopes, γ, α, and β representing respectively the inner, intermediate, and outer slopes (we specify the priors adopted in Table 3). Since the parameters are strongly degenerate, it was advised to fix at least one of the slopes (Arnaud et al. 2010); therefore, we fixed the central slope α to the best fit value of 1.3 estimated by the Planck Collaboration (Planck Collaboration Int. V 2013). The resulting best fit, the intrinsic scatter around the median profile, and the slope computed from the fit are shown in Fig. 4. The posterior distributions of the parameters and the covariance between them are shown in Fig. E.2.

Similar to the case of the gas density, we find that the slope of the profiles steepens steadily with radius, as expected for a gas in hydrostatic equilibrium within a NFW potential. The best fit with the generalized NFW functional form does an excellent job of reproducing the slopes estimated from the piecewise power law fits.

In the range where pressure measurements are available both from XMM-Newton and Planck, we find excellent agreement between the two (see Appendix A) even though the pressure profiles were obtained using completely independent probes. This shows that X-ray and SZ observations provide a consistent picture of the state of the ICM and gives us confidence that systematics in our measurements are small and under control.

3.4. Temperature

Temperature profiles in X-ray studies are usually obtained by performing spectral fitting in concentric annuli (see Sect. 2.4.2), called the spectroscopic temperature. In addition, we also use our Planck SZ pressure profiles and combine them with the X-ray density profiles to obtain gas-mass-weighted temperatures Tgmw = PSZ/ne, which are then projected (using the X-ray emissivity as weight, as in Mazzotta et al. 2004) on the plane of the sky and overplotted on the spectroscopic temperatures. Our X-ray and SZ measurements of pressure and temperature cover different radial ranges. X-ray spectroscopy probes the temperature of the gas within R500, while SZ probes temperatures from 0.7R500 to 2R500 (excluding the first three SZ data points), which highlights the complementarity of the two ICM diagnostics. In Fig. 5 we show our 2D spectral temperature profiles rescaled by T500, defined as

(10)

While density and pressure change by three to four orders of magnitude going from the center to the outskirts of the cluster, temperature shows much milder variations. In particular, it is almost constant out to ∼0.5R500, and then declines beyond this point.

In Fig. 5 we show the results of the piecewise power law fits in several radial ranges. We perform a global fit to the temperature profiles with the functional form described in Vikhlinin et al. (2006), which is able to describe the temperature profiles of both the core and the outer parts of galaxy clusters

(11)

with x = R/R500, and six free parameters: T0, Tmin, rcool, acool, rt, and c (we specify the priors adopted in Table 3). The posterior distribution of these parameters and the covariances are shown in Fig. E.3. This functional form provides an accurate description of the shape of the temperature profiles. The slopes estimated from the global fit follow the slopes measured from the piecewise power law fits at different radial ranges within 1σ at every radius.

The average slope of the temperature profiles is slightly positive in the central regions because of effects due to cooling, especially in cool-core clusters. Beyond ∼0.5R500 the slope remains relatively flat at a value of −0.3. This value is consistent with the slopes measured inside R500 from XMM-Newton and Chandra data (Leccardi & Molendi 2008; Pratt et al. 2007), but it is flatter than the typical slopes measured in Suzaku data. From a collection of a dozen clusters with published Suzaku temperature profiles, Reiprich et al. (2013) report that the universal shape of the profiles can be well described by the form T/⟨T⟩=1.19 − 0.84(R/R200), i.e., the data are consistent with a slope of −1.0 in the outskirts, which is much steeper than the results presented here. Again, the low angular resolution of Suzaku may have prevented the authors from removing cool, overdense regions that could bias low the measured spectroscopic temperature. On the other hand, SZ pressure profiles are much less sensitive to gas clumping (e.g., Khedekar et al. 2013; Roncarelli et al. 2013) and our density profiles were corrected for the statistical effect of gas clumping (see Sect. 2.2.2), thus our X/SZ temperatures are closer to gas-mass-weighted temperatures (see the discussion in Adam et al. 2017).

3.5. Entropy

Entropy traces the thermal evolution of the ICM plasma, which can be altered via cooling/heating, mixing, and convection. Simulations predict that in the presence of non-radiative processes only, entropy increases steadily with radius out to ∼2 × R200, following a power law with a slope of 1.1 (Tozzi & Norman 2001; Voit et al. 2005; Borgani et al. 2005; Lau et al. 2015). The entropy profiles of the cluster population should scale self-similarly when rescaled by the quantity

(12)

Similarly to pressure and temperature, entropy can be recovered from X-ray-only quantities as (using the deprojected temperature; see Sect. 2.4) or by combining SZ pressure with X-ray density as (ignoring the first three Planck points). We show our scaled entropy profiles in Fig. 6. Our profiles very closely match the predicted power law model beyond 0.3R500, with just A2319 showing a significant flattening not compatible within the error bars. In the central regions our profiles flatten, with non-cool-core clusters flattening more than cool-core clusters, as already noted in numerous studies (e.g., Pratt et al. 2010; Cavagnolo et al. 2009), and references therein).

By fitting the profile using piecewise power laws we observe a gradual steepening of the entropy slope, which becomes consistent with the predictions of gravitational collapse (Voit et al. 2005) beyond ∼0.5R500, i.e., from a slope of ∼0.6 in the core to ∼1.1 in the outskirts. As for the previous cases, we fitted our profiles throughout the entire radial range with the functional form introduced by Cavagnolo et al. (2009), which consists of a power law with a constant entropy floor

(13)

with x = R/R500, and three free parameters K0, K1, and α (we specify the priors adopted in Table 3). The posterior distributions of the parameters and the covariances are shown in Fig. E.4. We note that this functional form does not provide an accurate description of the data in the outer parts of the profiles where it is not able to follow the observed gradual change in slope throughout the radial range covered. At large radii, the best fitting slope using Eq. (13) reads α = 0.84 ± 0.04 (see Table 3), whereas the data prefer a slope consistent with the self-similar prediction of 1.1 beyond 0.6R500 (see Table 2).

3.6. Scatter

The high data quality of X-COP allows us to probe the intrinsic scatter of our profiles as a function of radius. The piecewise fit using power laws allows us to measure the scatter in a nearly model independent way, whereas the global fit with functional forms and scatter described as a log-parabola provides a consistent description of both the profile shape and the intrinsic scatter throughout the entire radial range. In Fig. 7 we show the scatter of all our thermodynamic quantities obtained in both cases. We recall that our definition of the intrinsic scatter is relative: a value of 0.1 on the y-axis indicates that the considered quantity is intrinsically scattered by 10% of its value.

thumbnail Fig. 7.

Measured intrinsic scatter of all our thermodynamic quantities, density (red), pressure (green), temperature (blue), and entropy (black). The data points indicate the results of piecewise power law fits on several radial ranges, whereas the dashed lines and shaded areas show the intrinsic scatter described as a log-parabola (Eq. (6)) around the best fitting functional forms.

We notice that our thermodynamic profiles generally exhibit a high scatter in the central parts of the profile. The scatter decreases toward the outskirts, reaching a minimum in the range [0.2 − 0.8]R500, and increases slightly beyond this point. We find that temperature is the least scattered thermodynamic quantity, with intrinsic scatters ranging from 10% to 20%. On the contrary, and surprisingly, pressure is the most scattered quantity at all radii (looking at the scatter reconstructed from the piecewise power law fits), ranging from 25% to 60%.

In all cases, we note that our profiles present a high degree of self-similarity in the radial range [0.2 − 0.8]R500, with a typical intrinsic scatter less than 0.3 (∼0.1 dex) in all the measured quantities. This radial range corresponds to the region where gravity dominates and baryonic physics (gas cooling, AGN, and supernova feedback) is relatively unimportant, whereas gas accretion still plays a subdominant role. This is consistent with tightly self-regulated mechanical AGN feedback (e.g., via chaotic cold accretion, Gaspari et al. 2012), which can only affect the region < 0.1R500, with predicted moderate scatter in T/ne as similarly retrieved here.

We note that the intrinsic scatter profiles shown in Fig. 7 include the scatter that is induced by uncertainties on the cluster mass, hence on the self-similar scaling quantities. In Appendix C we estimate numerically the residual scatter coming from uncertainties in the self-similar scaling on the various thermodynamic quantities. We found that the scaled pressure is the quantity that is most strongly affected by mass uncertainties, which introduce a scatter of ∼11% at R500, compared to 6% for the temperature, 5% for the density, and 3% for the entropy. However, this effect appears insufficient to fully explain the difference in scatter between density and pressure at 0.5R500. We also checked whether the higher scatter in pressure could be explained by intrinsic differences between X-ray and SZ pressure profiles (see Appendix A). However, we find no statistically significant differences between the pressure profiles measured with the two methods, and the scatter in pressure remains the same when considering X-ray and SZ data separately.

3.7. CC versus NCC

We divide our cluster sample into two populations based on the central entropy value, shown as the last column in Table 1. We use as an indicator of dynamical state the central entropy of our clusters as measured by Chandra (Cavagnolo et al. 2009), which has a better spatial resolution than XMM-Newton, and therefore is able to trace more accurately the behavior of the entropy profiles in the inner regions. Using this indicator we identify four clusters as cool-core (CC) and eight as non-cool-core (NCC) using the value of K0 = 30 keV cm2 as the cutoff value.

In Fig. 8 we show the data split into the CC and NCC populations, together with the fit using the functional forms used above, Eq. (7) for density, Eq. (9) for pressure, Eq. (11) for temperature, and Eq. (13) for entropy. The best fitting functional forms for the CC and NCC classes separately are provided in Table 3, and the results of piecewise power law fits to the two populations individually are given in Table D.1 and Fig D.1. We note that in the core, the CC and NCC systems separate out. However, we do not observe any significant differences between CC and NCC systems outside the core: the properties of our SZ selected clusters beyond 0.3R500 are not influenced by the properties of the core. We remark that in the case of the temperature there is a slight difference between the two best fits, with NCC having steeper temperature profiles, however well within the 1σ error envelope. The only marginally significant difference is found in the entropy profiles, which appear slightly flatter in the outskirts of NCC clusters. As shown in Table 3, we measure an outer slope αCC = 0.95 ± 0.03 for the CC populations, as opposed to αNCC = 0.85 ± 0.07. However, we note that this difference can be an artifact of the poor fit to the data obtained with a simple power law with an entropy floor (Eq. (13)). Indeed, similar to the case of the fit to the overall population, we find a steeper slope at large radii when fitting the data points for the two populations with a piecewise power law (αCC = 1.23 ± 0.14, αNCC = 0.94 ± 0.14, see Table D.1), which is consistent with the self-similar slope of 1.1 within 1σ. Thus, the evidence for a flatter entropy slope beyond R500 in the NCC population is marginal.

thumbnail Fig. 8.

Thermodynamic quantities of the X-COP clusters dividing the clusters into cool-core (CC, in blue) and non-cool-core (NCC, in red) populations, compared with the entire population (ALL, in green). Top left panel: density profiles fitted using the functional form presented in Eq. (7), overplotted on the data and on the “universal” density profile of Eckert et al. (2012, pink shaded area). Top right panel: pressure profiles fitted using the functional form presented in Eq. (9), overplotted on the data and compared to the Planck results (Planck Collaboration Int. V 2013), dotted black line) and the universal pressure profile of Arnaud et al. (2010, dash-dotted black line). Bottom left panel: temperature profiles fitted using the functional form presented in Eq. (7) and overplotted on the data. Bottom right panel: entropy profiles fitted using the functional form presented in Eq. (13), overplotted on the data and on the gravitational collapse predictions (Voit et al. 2005, solid black line).

4. Discussion

4.1. Systematic uncertainties

In this section we describe the potential systematic errors affecting our analysis.

  • Gas density. As described in Sect. 2.2.2, we paid special attention to the minimization of the systematics in the subtraction of the XMM-Newton background. The method that we used to model the contribution of each individual background component was calibrated using a large set of ∼500 blank-sky pointings and leads to residual systematics on the order of 3% on the subtraction of the local background (see Appendix A in Ghirardini et al. 2018). For the present work, we conservatively increased the level of systematics to 5% to include potential uncertainties associated with the application of the method to a cluster field instead of a blank field. A systematic error of 5% of the background value was thus added in quadrature to all our surface brightness measurements. We note that the systematic uncertainty becomes comparable to the statistical errors only beyond ∼2 × R500. At R200 the systematic uncertainty is typically 20% or less of the measured signal. Further improvements in the modeling of the XMM-Newton background could allow us in the future to provide information beyond the current limiting radii since in many cases our SZ pressure profiles extend beyond 2 × R500.

  • Pressure profiles. A possible source of systematics on the reconstruction of SZ pressure profiles is the relativistic corrections to the SZ effect (Itoh et al. 1998), which reduce the amplitude of the SZ increment in the high-frequency part of the CMB spectrum. Several recent works claimed a detection of the relativistic SZ corrections on stacked Planck data (Hurier 2016; Erler et al. 2018). In particular, Erler et al. (2018) noted that the relativistic corrections could lead to an underestimate of the integrated SZ signal up to 15% for the hottest clusters, which could thus affect our pressure profiles too. However, we note that the gas temperature decreases by a factor of 2 − 2.5 from the core to the outskirts, such that the impact of SZ corrections should be limited to the central regions, where spectroscopic X-ray measurements are preferred because of their higher S/N and resolution. For typical temperatures of ∼5 keV at R500 and beyond the expected effect is less than 5% (Erler et al. 2018). For more discussion on the impact of systematic uncertainties we refer to Planck Collaboration Int. V (2013).

  • Spectroscopic temperatures. Although our modeling of the XMM-Newton spectra is fairly sophisticated (see Sect. 2.4.2), uncertainties in the subtraction of the XMM-Newton background can lead to systematics in our spectral measurements in the outermost regions considered. Following Leccardi & Molendi (2008) we do not attempt to perform spectral measurements in the regions where our signal is less than 60% of the background intensity to avoid introducing biases. Another potential source of systematics is the calibration of the telescope’s effective area. Schellenberger et al. (2015) reported systematic differences at the level of 15% between XMM-Newton and Chandra temperature measurements for the same regions, Chandra returning systematically higher temperatures than XMM-Newton. As shown in Fig. 4 and demonstrated in Appendix A, we observe a very good agreement between XMM-Newton and Planck pressure profiles; the only exception is ZwCl 1215, for which the pressure measured by XMM-Newton actually exceeds the SZ pressure by ∼20%, which could be explained by orientation effects since the X-ray and SZ signals have different line-of-sight dependencies. Since our X-ray and SZ pressure profiles are obtained in an independent way from different instruments and different techniques, we conclude that our spectral measurements are robust.

  • Self-similar scaling. Given that the scaling quantities depend on the measured mass, and that we use as our reference mass model the backward NFW mass model (Ettori et al. 2010; Ghirardini et al. 2018), uncertainties on the mass measurements should be taken into account. In Eckert et al. (2019) we discuss the accuracy of our mass models by testing our mass measurements using various methods (forward fitting, Gaussian processes, and several functional forms for the mass model). We find that all the methods agree with the NFW mass reconstruction, with the mass profiles scattered by less than 5% at a fixed radius of 1.5 Mpc. The uncertainty in our scaling is therefore less than 3% on P500 and K500, and less than 2% on R500. In Eckert et al. (2019) we also assess the level of non-thermal pressure support by comparing the X-COP gas fraction profiles with the expected universal gas fraction. We find that the bias in our mass measurements at R500 is just 6% on average, again resulting in very small corrections in the self-similar quantities.

4.2. Regular outskirts

The wide radial range accessible with the X-COP data allows us to study the properties of the gas at R500 and beyond and to constrain the shape of the universal thermodynamic profiles throughout the entire cluster volume for the first time. Compared to previous works addressing the state of the gas in cluster outskirts (e.g., with Suzaku data) the study presented here constitutes a substantial improvement in several ways: (i) the ability of our azimuthal median method to excise overdense regions down to scales of 10–20 kpc depending on the cluster redshift (Eckert et al. 2015), which allows us to measure gas density profiles that are free of the effects of gas clumping on the scales we are able to resolve (typically 30 kpc); (ii) a nearly uniform azimuthal coverage for all our clusters out to R200, which guarantees that our measurements are representative of the global behavior and were not obtained along preferential directions; and (iii) an exquisite control of systematic uncertainties even in the faint cluster outskirts regime (see above).

As described in Sect. 3, our reconstruction of clumping-free thermodynamic quantities leads to results that differ substantially from the typical results obtained with Suzaku. We recall that several studies found relatively flat density profiles, steep temperature profiles, and entropy profiles that fall below the prediction of gravitational collapse and sometimes even roll over (e.g., Kawaharada et al. 2010; Walker et al. 2012a,b; Urban et al. 2014; Simionescu et al. 2017). Conversely, our clumping-corrected reconstruction yields density and pressure profiles that steepen steadily with radius (see Figs. 3 and 4), temperature profile decreasing with a mild slope of −0.3 that is consistent with the slopes observed inside R500 by XMM-Newton and Chandra (Leccardi & Molendi 2008; Pratt et al. 2007; Vikhlinin et al. 2006), and entropy profiles rising with a slope that is consistent with the self-similar slope of 1.1 beyond 0.6R500 and all the way out to the largest radii considered (2 × R500).

All the results presented here point to gas clumping as the primary origin for the deviations from the predictions reported so far by Suzaku, in agreement with the results presented in Tchernin et al. (2016) for the case of Abell 2142. The low resolution of Suzaku (∼2 arcmin) indeed prevented the authors from excising cool, overdense structures that would bias at the same time the gas density toward high values and the spectroscopic temperature toward low values, resulting in underestimated values for the entropy that are not representative of the bulk of the ICM. If the gas in such infalling structures is in pressure equilibrium with its environment, as usually predicted (e.g., Roncarelli et al. 2013; Planelles et al. 2017), pressure profiles reconstructed from the SZ effect are mildly affected by such inhomogeneities and the combination of SZ pressure and clumping-free gas density is representative of the state of the ICM well beyond R500.

The only exception to this scenario is the case of Abell 2319 (Ghirardini et al. 2018), which deviates systematically from the measured universal profiles even when the profiles are corrected for clumping. In Ghirardini et al. (2018) we showed that the deviations from self-similarity cannot be explained by azimuthal variations, but rather that the ongoing merging activity causes a high level of non-thermal pressure support. This conclusion is supported by the high hydrostatic gas fraction of this system and a clear deficit of entropy beyond R500, even after excising clumps. Abell 2319 is the only system within the X-COP sample that exhibits such a behavior (see also Ettori et al. 2019), which suggests that this system is currently experiencing a transient phase of high non-thermal pressure induced by a violent merger with a mass ratio of 3 to 1 (Oegerle et al. 1995).

Overall, the results presented here establish that in the majority of cases, the bulk of the ICM is virialized and follows the predictions of gravitational collapse out to 2 ×  R500 ≈  R100. Accretion shocks that are expected to raise the entropy level of the smooth infalling gas should be located approximately at 3 − 4 × R500 (Lau et al. 2015), and we would expect the entropy of the ICM to turn over around this radius. These radii should also correspond to the approximate location of the splashback radius (Diemer & Kravtsov 2014; Diemer et al. 2017), which represents a natural boundary of dark matter halos. Future X-ray and SZ facilities such as ATHENA (Nandra et al. 2013) and CMB-S4 (Abazajian et al. 2016) will attempt to detect the ICM at the cluster boundary to constrain the location of accretion shocks and the accretion rate. The results presented here highlight the need for relatively high angular resolution experiments with a low and highly reproducible background to reach these goals.

4.3. Self-similarity of the profiles

Our analysis shows that the thermodynamic profiles exhibit a high level of similarity once the profiles are rescaled according to the self-similar model (Kaiser 1986). The level of self-similarity is particularly remarkable beyond the core (R >  0.3R500) and it reaches a maximum in the radial range [0.2 − 0.8]R500. As already discussed in Sect. 3.6, the region of minimum scatter observed in this study corresponds to the region where the gas is highly virialized and baryonic effects are negligible. In the central regions (R <  0.3R500) baryonic effects (cooling, AGN feedback) lead to a substantial scatter within the cluster population. Beyond ∼R500, we again observe an increase in the measured scatter, which might be explained by different accretion rates from one system to another. Importantly, our study shows that the properties of the X-COP cluster population beyond 0.3R500 are not correlated with the core state (CC or NCC). While the core state probably retains memory of past major mergers, it does not trace the accretion rate on large scales at the present epoch. This result agrees with the predictions of Planelles et al. (2017), which did not find differences in the accretion rate of simulated CC and NCC systems. For instance, the case of A2029 is striking. This cluster hosts a strong cool-core and it is one of the most regular in our sample. However, our large-scale mosaic reveals that it is located within a chain of at least three X-ray detected structures (see Fig. 2) with overlapping R200, and the optical information shows that this system is part of a larger filamentary structure extending over ∼20 Mpc (Smith et al. 2012).

Another important result of our study is that beyond the central regions pressure is the most scattered thermodynamic quantity (see Fig. 7). The scatter in Pe = TX × ne is about 50% larger than the scatter in either TX or ne, which is expected when fluctuations in temperature and density are uncorrelated. This result is opposite to the widely accepted view that temperature and density variations are anti-correlated, which has lead people to postulate that the quantity YX = Mgas × TX has the lowest scatter at fixed mass (Kravtsov et al. 2006). Our results disagree with this conclusion and imply that the scatter in Mgas is less than the scatter in YX at fixed mass. These results are consistent with the recent predictions of Truong et al. (2018), which found that in the simulation runs including gas cooling and subgrid thermal AGN feedback, temperature and density are essentially uncorrelated (see their Fig. 10), implying that the scatter in Mgas and TX is less than that in YX. Beyond the core, X-ray observables appear to behave self-similarly to a high level of precision. In the case a selection based on the integrated gas mass or the core-excised X-ray luminosity can be achieved, future X-ray surveys such as eROSITA (Predehl et al. 2010) will yield large cluster catalogues and low-scatter mass proxies, even in comparison to SZ surveys (Mantz et al. 2018).

5. Conclusion

In this paper, we presented the universal thermodynamic properties of the ICM for 12 SZ selected galaxy clusters observed with XMM-Newton and Planck. Our observational strategy allowed us to construct radial profiles of gas density, pressure, temperature, and entropy over an unprecedentedly wide radial range from 0.01R500 to 2 × R500, i.e., covering the entire cluster volume. We fitted our self-similar scaled profiles with universal functional forms and provided estimates of the radial dependence of the slope and intrinsic scatter. Our findings can be summarized as follows:

  • Our gas density and pressure profiles are in excellent agreement with previous determinations of the universal density (Eckert et al. 2012) and pressure profiles (Planck Collaboration Int. V 2013). The typical uncertainties in the gas density and pressure at R200 are at the level of 10%, allowing us to perform a detailed analysis of the shape and intrinsic scatter.

  • The logarithmic slope of the density and pressure profiles steepens steadily with radius, reaching a value of −2.5 and −3.0 at R200 for density and pressure, respectively. These results are consistent with the expectations for an ideal gas in hydrostatic equilibrium within a NFW potential well.

  • Beyond ∼0.3R500 the temperature profiles decrease steadily with radius with a logarithmic slope of −0.3, which is somewhat shallower than the slope of ∼ − 1.0 observed in the outer regions of several systems with Suzaku (Reiprich et al. 2013).

  • With the exception of one system, beyond ∼0.5R500 all clusters follow the gravitational collapse prediction for entropy generation in galaxy clusters (Voit et al. 2005) out to the largest radii considered (2 × R500). This result is at odds with the conclusions usually reached from Suzaku observation, which often show a deficit of entropy beyond R500. The difference is explained by the steep slope of the Suzaku temperature profiles compared to ours and by our treatment of gas clumping. We postulate that the impossibility of properly excising clumps in low-resolution Suzaku data is responsible for biasing the observed temperatures low and gas densities high.

  • The outer regions of galaxy clusters exhibit a high level of self-similarity. Beyond ∼0.3R500 we find no significant difference between the cool-core and non-cool-core cluster populations in any of the quantities of interest. This result implies that the core properties are determined by the merging history of a system but do not trace efficiently the current accretion rate, which determines the state of the gas in the outskirts.

  • We determined for the first time the scatter of each thermodynamical quantity within the cluster population as a function of radius. The scatter of all quantities is maximum in the core and reaches a minimum in the radial range [0.2 − 0.8]R500 (see Table 2 and Fig. 7). We find that the gas temperature is the least scattered quantity at nearly all radii.

A recently accepted XMM-Newton program will extend the X-COP sample to objects that were initially excluded (A754, A3667, and A3827), which will allow us to perform a similar analysis on a statistically complete SZ-selected sample. Furthermore, since pressure profiles extend beyond 2 × R500, a further reduction of the systematics on the surface brightness profile would be useful to extend the thermodynamic profiles beyond the current limits, provided that observations with higher statistical quality can be performed.


1

Where θ500 refers to the apparent angular size of R500.

2

Where NH refers to the hydrogen column density along the line of sight.

Acknowledgments

X-COP data products are available for download at https://www.astro.unige.ch/xcop. We thank Alexis Finoguenov, Kaustuv Basu, Jens Erler, and Gus Evrard for useful discussions. The research leading to these results has received funding from the European Union’s Horizon 2020 Programme under the AHEAD project (grant agreement n. 654215). Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. S.E. acknowledges financial contribution from the contracts NARO15 ASI-INAF I/037/12/0, ASI 2015-046-R.0 and ASI-INAF n.2017-14-H.0. F.V. acknowledges financial support from the ERC Starting Grant “MAGCOW”, no. 714196. M.G. is supported by NASA through Einstein Postdoctoral Fellowship Award Number PF5-160137 issued by the Chandra X-ray Observatory Center, which is operated by the SAO for and on behalf of NASA under contract NAS8-03060. Support for this work was also provided by Chandra grant GO7-18121X. E.R. acknowledges the ExaNeSt and Euro Exa projects, funded by the European Union’s Horizon 2020 research and innovation program under grant agreement No. 671553 and No. 754337 and financial contribution from the agreement ASI-INAF n.2017-14-H.0. H.B. and P.M. acknowledge financial support by ASI Grant 2016-24-H.0.

References

  1. Abazajian, K. N., Adshead, P., & Ahmed, Z. 2016, ArXiv e-prints [arXiv:1610.02743] [Google Scholar]
  2. Adam, R., Comis, B., Macías-Pérez, J.-F., et al. 2015, A&A, 576, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Adam, R., Arnaud, M., Bartalucci, I., et al. 2017, A&A, 606, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Akamatsu, H., Hoshino, A., Ishisaki, Y., et al. 2011, PASJ, 63, 1019 [NASA ADS] [Google Scholar]
  5. Ameglio, S., Borgani, S., Pierpaoli, E., & Dolag, K. 2007, MNRAS, 382, 397 [NASA ADS] [CrossRef] [Google Scholar]
  6. Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]
  7. Andrade-Santos, F., Jones, C., Forman, W. R., et al. 2017, ApJ, 843, 76 [NASA ADS] [CrossRef] [Google Scholar]
  8. Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Avestruz, C., Nagai, D., Lau, E. T., & Nelson, K. 2015, ApJ, 808, 176 [NASA ADS] [CrossRef] [Google Scholar]
  10. Barnes, D. J., Vogelsberger, M., & Kannan, R. 2018, MNRAS, 481, 1809 [NASA ADS] [CrossRef] [Google Scholar]
  11. Basu, K., Zhang, Y.-Y., Sommer, M. W., et al. 2010, A&A, 519, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Borgani, S., Murante, G., Springel, V., et al. 2004, MNRAS, 348, 1078 [NASA ADS] [CrossRef] [Google Scholar]
  13. Borgani, S., Finoguenov, A., Kay, S. T., et al. 2005, MNRAS, 361, 233 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bourdin, H., Mazzotta, P., Kozmanyan, A., Jones, C., & Vikhlinin, A. 2017, ApJ, 843, 72 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]
  16. Burns, J. O., Skillman, S. W., & O’Shea, B. W. 2010, ApJ, 721, 1105 [Google Scholar]
  17. Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cash, W. 1979, ApJ, 228, 939 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cavagnolo, K. W., Donahue, M., Voit, G. M., & Sun, M. 2009, ApJS, 182, 12 [NASA ADS] [CrossRef] [Google Scholar]
  20. Croston, J. H., Arnaud, M., Pointecouteau, E., & Pratt, G. W. 2006, A&A, 459, 1007 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. De Grandi, S., & Molendi, S. 2002, ApJ, 567, 163 [NASA ADS] [CrossRef] [Google Scholar]
  22. De Luca, A., & Molendi, S. 2004, A&A, 419, 837 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Diaz-Rodriguez, J., Eckert, D., Monajemi, H., Paltani, M., & Sardy, S. 2017, ArXiv e-prints [arXiv:1703.00654] [Google Scholar]
  24. Diemer, B., & Kravtsov, A. V. 2014, ApJ, 789, 1 [NASA ADS] [CrossRef] [Google Scholar]
  25. Diemer, B., Mansfield, P., Kravtsov, A. V., & More, S. 2017, ApJ, 843, 140 [NASA ADS] [CrossRef] [Google Scholar]
  26. Dolag, K., Bartelmann, M., & Lesch, H. 1999, A&A, 348, 351 [NASA ADS] [Google Scholar]
  27. Dolag, K., Vazza, F., Brunetti, G., & Tormen, G. 2005, MNRAS, 364, 753 [NASA ADS] [CrossRef] [Google Scholar]
  28. Ebeling, H., White, D. A., & Rangarajan, F. V. N. 2006, MNRAS, 368, 65 [NASA ADS] [Google Scholar]
  29. Eckert, D., Vazza, F., Ettori, S., et al. 2012, A&A, 541, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Eckert, D., Molendi, S., Vazza, F., Ettori, S., & Paltani, S. 2013, A&A, 551, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Eckert, D., Molendi, S., Owers, M., et al. 2014, A&A, 570, A119 [Google Scholar]
  32. Eckert, D., Roncarelli, M., Ettori, S., et al. 2015, MNRAS, 447, 2198 [NASA ADS] [CrossRef] [Google Scholar]
  33. Eckert, D., Ettori, S., Coupon, J., et al. 2016, A&A, 592, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Eckert, D., Ettori, S., Pointecouteau, E., et al. 2017, Astron. Nachr., 338, 293 [NASA ADS] [CrossRef] [Google Scholar]
  35. Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Erler, J., Basu, K., Chluba, J., & Bertoldi, F. 2018, MNRAS, 476, 3360 [NASA ADS] [CrossRef] [Google Scholar]
  37. Ettori, S., De Grandi, S., & Molendi, S. 2002, A&A, 391, 841 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Ettori, S., Gastaldello, F., Leccardi, A., et al. 2010, A&A, 524, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Ettori, S., Ghirardini, V., Eckert, D., et al. 2019, A&A, 621, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
  41. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
  42. Frenk, C. S., White, S. D. M., Bode, P., et al. 1999, ApJ, 525, 554 [NASA ADS] [CrossRef] [Google Scholar]
  43. Gaspari, M., Ruszkowski, M., & Sharma, P. 2012, ApJ, 746, 94 [Google Scholar]
  44. Gaspari, M., Brighenti, F., Temi, P., & Ettori, S. 2014, ApJ, 783, L10 [NASA ADS] [CrossRef] [Google Scholar]
  45. Gaspari, M., McDonald, M., Hamer, S. L., et al. 2018, ApJ, 854, 167 [NASA ADS] [CrossRef] [Google Scholar]
  46. Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Hahn, O., Martizzi, D., Wu, H.-Y., et al. 2017, MNRAS, 470, 166 [NASA ADS] [Google Scholar]
  48. Hoshino, A., Patrick Henry, J., Sato, K., et al. 2010, PASJ, 62, 371 [NASA ADS] [CrossRef] [Google Scholar]
  49. Hurier, G. 2016, A&A, 596, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Hurier, G., Macías-Pérez, J. F., & Hildebrandt, S. 2013, A&A, 558, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Itoh, N., Kohyama, Y., & Nozawa, S. 1998, ApJ, 502, 7 [NASA ADS] [CrossRef] [Google Scholar]
  52. Jones, C., & Forman, W. 1984, ApJ, 276, 38 [NASA ADS] [CrossRef] [Google Scholar]
  53. Kaiser, N. 1986, MNRAS, 222, 323 [NASA ADS] [CrossRef] [Google Scholar]
  54. Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Kawaharada, M., Okabe, N., Umetsu, K., et al. 2010, ApJ, 714, 423 [NASA ADS] [CrossRef] [Google Scholar]
  56. Kay, S. T., Pearce, F. R., Frenk, C. S., & Jenkins, A. 2002, MNRAS, 330, 113 [NASA ADS] [CrossRef] [Google Scholar]
  57. Khatri, R., & Gaspari, M. 2016, MNRAS, 463, 655 [NASA ADS] [CrossRef] [Google Scholar]
  58. Khedekar, S., Churazov, E., Kravtsov, A., et al. 2013, MNRAS, 431, 954 [NASA ADS] [CrossRef] [Google Scholar]
  59. Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353 [NASA ADS] [CrossRef] [Google Scholar]
  60. Kravtsov, A. V., Nagai, D., & Vikhlinin, A. A. 2005, ApJ, 625, 588 [NASA ADS] [CrossRef] [Google Scholar]
  61. Kravtsov, A. V., Vikhlinin, A., & Nagai, D. 2006, ApJ, 650, 128 [NASA ADS] [CrossRef] [Google Scholar]
  62. Kriss, G. A., Cioffi, D. F., & Canizares, C. R. 1983, ApJ, 272, 439 [NASA ADS] [CrossRef] [Google Scholar]
  63. Lamarre, J.-M., Puget, J.-L., Ade, P. A. R., et al. 2010, A&A, 520, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Lapi, A., Fusco-Femiano, R., & Cavaliere, A. 2010, A&A, 516, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Lau, E. T., Nagai, D., Avestruz, C., Nelson, K., & Vikhlinin, A. 2015, ApJ, 806, 68 [NASA ADS] [CrossRef] [Google Scholar]
  66. Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2014, MNRAS, 441, 1270 [NASA ADS] [CrossRef] [Google Scholar]
  67. Leccardi, A., & Molendi, S. 2008, A&A, 486, 359 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Lovisari, L., Forman, W. R., Jones, C., et al. 2017, ApJ, 846, 51 [NASA ADS] [CrossRef] [Google Scholar]
  69. Mantz, A. B., Allen, S. W., Morris, R. G., & von der Linden, A. 2018, MNRAS, 473, 3072 [NASA ADS] [CrossRef] [Google Scholar]
  70. Mazzotta, P., Rasia, E., Moscardini, L., & Tormen, G. 2004, MNRAS, 354, 10 [NASA ADS] [CrossRef] [Google Scholar]
  71. McCammon, D., Almy, R., Apodaca, E., et al. 2002, ApJ, 576, 188 [Google Scholar]
  72. Morandi, A., Sun, M., Forman, W., & Jones, C. 2015, MNRAS, 450, 2261 [CrossRef] [Google Scholar]
  73. Moretti, A., Campana, S., Lazzati, D., & Tagliaferri, G. 2003, ApJ, 588, 696 [NASA ADS] [CrossRef] [Google Scholar]
  74. Nagai, D., & Lau, E. T. 2011, ApJ, 731, L10 [NASA ADS] [CrossRef] [Google Scholar]
  75. Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1 [NASA ADS] [CrossRef] [Google Scholar]
  76. Nandra, K., Barret, D., & Barcons, X. 2013, ArXiv e-prints [arXiv:1306.2307] [Google Scholar]
  77. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [NASA ADS] [CrossRef] [Google Scholar]
  78. Oegerle, W. R., Hill, J. M., & Fitchett, M. J. 1995, AJ, 110, 32 [NASA ADS] [CrossRef] [Google Scholar]
  79. Okabe, N., Umetsu, K., Tamura, T., et al. 2014, PASJ, 66, 99 [NASA ADS] [Google Scholar]
  80. Pfrommer, C., Springel, V., Jubelgas, M., & Ensslin, T. A. 2007, in Cosmic Frontiers, eds. N. Metcalfe, & T. Shanks, ASP Conf. Ser., 379, 221 [NASA ADS] [Google Scholar]
  81. Planck Collaboration VIII. 2011, A&A, 536, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Planck Collaboration XXIX. 2014, A&A, 571, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Planck Collaboration I. 2016, A&A, 594, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Planck Collaboration VII. 2016, A&A, 594, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Planck Collaboration Int. V. 2013, A&A, 550, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Planck HFI Core Team 2011, A&A, 536, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Planelles, S., Fabjan, D., Borgani, S., et al. 2017, MNRAS, 467, 3827 [NASA ADS] [CrossRef] [Google Scholar]
  89. Pratt, G. W., Böhringer, H., Croston, J. H., et al. 2007, A&A, 461, 71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Pratt, G. W., Arnaud, M., Piffaretti, R., et al. 2010, A&A, 511, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Predehl, P., Andritschke, R., & Böhringer, H. 2010, Proc. SPIE, 7732, 77320U [Google Scholar]
  92. Rasia, E., Ettori, S., Moscardini, L., et al. 2006, MNRAS, 369, 2013 [NASA ADS] [CrossRef] [Google Scholar]
  93. Reiprich, T. H., Basu, K., Ettori, S., et al. 2013, Space Sci. Rev., 177, 195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Roncarelli, M., Ettori, S., Borgani, S., et al. 2013, MNRAS, 432, 3030 [NASA ADS] [CrossRef] [Google Scholar]
  95. Rossetti, M., Gastaldello, F., Eckert, D., et al. 2017, MNRAS, 468, 1917 [NASA ADS] [CrossRef] [Google Scholar]
  96. Ruppin, F., Adam, R., Comis, B., et al. 2017, A&A, 597, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Salvetti, D., Marelli, M., Gastaldello, F., et al. 2017, Exp. Astron., 44, 309 [NASA ADS] [CrossRef] [Google Scholar]
  98. Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Simionescu, A., Allen, S. W., Mantz, A., et al. 2011, Science, 331, 1576 [NASA ADS] [CrossRef] [Google Scholar]
  100. Simionescu, A., Werner, N., Mantz, A., Allen, S. W., & Urban, O. 2017, MNRAS, 469, 1476 [NASA ADS] [CrossRef] [Google Scholar]
  101. Smith, A. G., Hopkins, A. M., Hunstead, R. W., & Pimbblet, K. A. 2012, MNRAS, 422, 25 [NASA ADS] [CrossRef] [Google Scholar]
  102. Snowden, S. L., Mushotzky, R. F., Kuntz, K. D., & Davis, D. S. 2008, A&A, 478, 615 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comm. Astrophys. Space Phys., 4, 173 [Google Scholar]
  104. Tauber, J. A., Mandolesi, N., Puget, J.-L., et al. 2010, A&A, 520, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Tchernin, C., Eckert, D., Ettori, S., et al. 2016, A&A, 595, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Tozzi, P., & Norman, C. 2001, ApJ, 546, 63 [NASA ADS] [CrossRef] [Google Scholar]
  107. Truong, N., Rasia, E., Mazzotta, P., et al. 2018, MNRAS, 474, 4089 [NASA ADS] [CrossRef] [Google Scholar]
  108. Urban, O., Simionescu, A., Werner, N., et al. 2014, MNRAS, 437, 3939 [NASA ADS] [CrossRef] [Google Scholar]
  109. Vazza, F., Brunetti, G., Kritsuk, A., et al. 2009, A&A, 504, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Vazza, F., Eckert, D., Simionescu, A., Brüggen, M., & Ettori, S. 2013, MNRAS, 429, 799 [NASA ADS] [CrossRef] [Google Scholar]
  111. Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691 [NASA ADS] [CrossRef] [Google Scholar]
  112. Voit, G. M. 2005, Rev. Mod. Phys., 77, 207 [NASA ADS] [CrossRef] [Google Scholar]
  113. Voit, G. M., Kay, S. T., & Bryan, G. L. 2005, MNRAS, 364, 909 [Google Scholar]
  114. Walker, S. A., Fabian, A. C., Sanders, J. S., & George, M. R. 2012a, MNRAS, 424, 1826 [NASA ADS] [CrossRef] [Google Scholar]
  115. Walker, S. A., Fabian, A. C., Sanders, J. S., George, M. R., & Tawara, Y. 2012b, MNRAS, 422, 3503 [NASA ADS] [CrossRef] [Google Scholar]
  116. Zhuravleva, I., Churazov, E., Kravtsov, A., et al. 2013, MNRAS, 428, 3274 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Consistency between X-ray and SZ pressure measurements

thumbnail Fig. A.1.

Posterior distribution of the parameter ηSZ ≈  PSZ/PX. The shaded blue area indicates the region containing 68% of the posterior distribution.

We checked the consistency between X-ray and SZ pressure profiles to test how our results are affected by discrepancies between the two measurements. To perform this check we introduced a parameter ηSZ, which is the ratio between the SZ and X-ray pressure profile, and proceeded with a joint fit, allowing the scatter on X-ray and SZ data to be independent. Mathematically we can write the following system of equations:

(A.1)

where ηSZ, σint, SZ, and σint, X are free parameters. For Pmodel we checked both the piecewise powerlaw fit case in the radial range where we have both X-ray and SZ measurements, and the global functional form on the entire radial range. In both cases the measured scatters are in good agreement, and are compatible with the scatter shown in Fig. 7. More importantly, the parameter ηSZ shows a distribution consistent with unity (see Fig. A.1), indicating a very good general agreement between X-ray and SZ pressure measurements.

Appendix B: Comparison between MILCA and forward-modeling pressure profiles

To test the robustness of our pressure profile measurements, we compared the MILCA results with the pressure measured using an alternative technique (see Bourdin et al. 2017). Following this technique, Planck-HFI frequency maps are first wavelet cleaned for CMB and thermal dust anisotropies, a parametric pressure template (Nagai et al. 2007) is subsequently projected onto the sky plane, convolved with the frequency dependent HFI beams and fitted to the CMB and dust cleaned maps. Being fully parametric, this technique allow us to take advantage of the frequency dependent angular resolution of each HFI channel during the template fitting. This angular resolution is about 9.7 and 7.3 arcmin at 100 and 143 GHz, respectively, but reaches about 5 arcmin in the energy range [217, 857] GHz (Planck Collaboration VII 2016).

In Fig. B.1 we compare the resulting best fitting pressure profile from the above procedure with the MILCA maps (see Sect. 2.5). The residuals are shown in the left panel and show the nice agreement, within the statistical uncertainty, of the two different methods applied. In the right panel all the residuals are grouped together to create a distribution, which is compared with a Gaussian centered at 0 and width 1, showing that the residuals follow this distribution very closely, indicating that statistically the two pressure profiles are in very good agreement.

thumbnail Fig. B.1.

Left panel: residuals of the comparison of the two different methods we used to estimate the SZ pressure profile. A remarkable good agreement within the statistical uncertainties is reached at all radii, especially excluding the first 3 MILCA points which are the most affected by the Planck PSF. Right panel: distribution of the residuals compared with the statistical prediction of a set of residuals: a Gaussian centered in zero and width one (red line).

Appendix C: Mass-induced scatter in thermodynamic profiles

Since the scaling of our thermodynamic quantities depends on the cluster mass both through the scale radius R500 and the self-similar quantities Q500 (see Sect. 3), the measured scatter profiles presented in Fig. 7 depend on the accuracy of the adopted masses. Both statistical and systematic fluctuations of the measured mass around the true mass will induce fluctuations of the scaling quantities, thus introducing an irreducible source of scatter originating from the limited accuracy of our mass calibration.

To take this effect into account, we estimated numerically the scatter in each quantity that is induced by uncertainties in our mass scaling. We started from scatter-free universal profiles following the measurements provided in Table 3, and perturbed the mass scaling according to the known statistical uncertainties and biases in our mass scaling. Namely, for each X-COP cluster we randomly drew new values of the observed mass Mobs as

(C.1)

with Mtrue the assumed true mass, b the hydrostatic mass bias, and ΔM the statistical uncertainty in our hydrostatic masses (see Ettori et al. 2019 for details). For the hydrostatic mass bias, we used the distribution of non-thermal pressure values determined in Eckert et al. (2019) from the measured gas fraction. We then scaled the scatter-free profiles for each quantity Q by the perturbed values of R500 and Q500.

thumbnail Fig. C.1.

Scatter in the various thermodynamic profiles induced by uncertainties in the mass calibration. The shaded areas show the range of intrinsic scatter obtained from 1000 simulations.

We applied this procedure to each X-COP cluster and computed the resulting scatter as a function of radius. We repeated this procedure 1000 times to get an idea of the uncertainty introduced by sample variance. In Fig. C.1 we show the resulting mass-induced scatter for the scaled pressure, density, temperature, and entropy. We can see that the effect of the mass scaling is largest on the pressure and ranges between 6% and 12% as a function of radius. Conversely, the effect on the entropy is minimal (2%–3%). The scatter in temperature and density induced by uncertainties in the mass scaling lies somewhere in between. Pressure is more affected than the other thermodynamic quantities simply because its slope is the most steep among the quantities.

Appendix D: Piecewise power law fits for CC and NCC clusters separately

In Fig. D.1 and Table D.1 we show the results of piecewise power law fits to the X-COP clusters divided into the CC and NCC populations. In Fig. D.2 we also show the best fit scatter of the populations as a function of radius, split into CC and NCC clusters and compared with the full population. In this case, we caution that the number of systems in each category is small (four CC and eight NCC systems) and the measurements of the scatter may be unreliable.

thumbnail Fig. D.1.

Same as Fig. 8, but using the piecewise fit instead of global functional form. The cluster population of X-COP is divided into cool-core (CC, in blue) and non-cool-core clusters (NCC, in red), compared with the entire population (ALL, in green).

thumbnail Fig. D.2.

Scatter of our thermodynamic quantities dividing our cluster sample in CC and NCC.

Table D.1.

Same as Table 2, but discriminating between cool-core and non-cool-core clusters.

Table D.2.

Stacked thermodynamic profiles.

Table D.2.

continued.

Appendix E: Marginalized posterior likelihood

thumbnail Fig. E.1.

Density – marginal probability Vikhlinin+06. Parameter distribution for the best fit on the density of all clusters using the Vikhlinin et al. (2006) functional form, Eq. (7). The priors on the parameters are shown in Table 3.

thumbnail Fig. E.2.

Pressure – marginal probability Nagai+07 – α fixed. Parameter distribution for the best fit on the pressure of all clusters using the Nagai et al. (2007) gNFW functional form, Eq. (9), fixing the intermediate slope α. The priors on the parameters are shown in Table 3.

thumbnail Fig. E.3.

Temperature – marginal probability Vikhlinin+06. Parameter distribution for the best fit on the temperature of all clusters using the Vikhlinin et al. (2006) functional form, Eq. (11). The priors on the parameters are shown in Table 3.

thumbnail Fig. E.4.

Entropy – marginal probability powerlaw + const. Parameter distribution for the best fit on the entropy of all clusters using a power law plus constant functional form, Eq. (13). The priors on the parameters are shown in Table 3.

Appendix F: Log of scientific observations

Table F.1.

Log of X-COP observations.

All Tables

Table 1.

Basic properties of the X-COP sample.

Table 2.

Results of the piecewise power law fits (normalizations, slopes, and intrinsic scatter; see Eq. (4)) for the various thermodynamic quantities in several radial ranges, marked by the inner and outer rescaled radii xin and xout.

Table 3.

Best fit parameters of the functional forms describing the universal thermodynamic quantities.

Table D.1.

Same as Table 2, but discriminating between cool-core and non-cool-core clusters.

Table D.2.

Stacked thermodynamic profiles.

Table D.2.

continued.

Table F.1.

Log of X-COP observations.

All Figures

thumbnail Fig. 1.

Flow chart of the XMM-Newton data analysis pipeline. The steps of the analysis are shown in red, the main intermediate and final products are described in the black boxes, and the procedures are shown in green italics.

In the text
thumbnail Fig. 2.

Adaptively smoothed and exposure corrected XMM-Newton mosaic images in the [0.7–1.2] keV energy band for all X-COP clusters. The superimposed white contours represent the Planck SZ S/N maps between 70 and 857 GHz. The contour levels correspond to 1, 3, 5, 7, 10, 15, 20, 30, 40, and 50σ. The spatial scale is given by the red thick lines, which have a common length of 15 arcmin. The color scale is given below the images in units of cts s−1 pixel−1, and is the same for all clusters.

In the text
thumbnail Fig. 3.

Top left panel: density profiles for all X-COP clusters obtained with two different deprojection methods: L1 regularization (data points) and multiscale fitting (solid lines). The magenta shaded area represents the scatter of the median profile in Eckert et al. (2012). Bottom part: ratio of the two methods for each individual system. Top right panel: joint fit to all the density profiles using piecewise power laws in several radial ranges (color-coded). The best fits and 1σ error envelope are shown by the solid lines, while the dashed lines represent the intrinsic scatter. Bottom left panel: joint fit to the density profiles using the functional form introduced by Vikhlinin et al. (2006), in red, with the shaded area indicating the 1σ error envelope around the best fit. The dashed lines represent the intrinsic scatter in the functional form as a function of radius. Bottom right panel: slope of the density profiles as a function of radius. The green data points show the results of the piecewise power law fits, whereas the red curve indicates the fit to the entire radial range using the Vikhlinin et al. (2006) functional form. In all panels, the vertical dashed and dotted lines represent the location of R500 and R200, respectively.

In the text
thumbnail Fig. 4.

Same as Fig. 3 for the pressure profiles rescaled by the self-similar quantity P500 (Eq. (8)). The squares indicate data points obtained from the deprojection of the SZ signal, while the filled circles are computed by combining the X-ray gas density profiles with the spectroscopic temperature. The solid red curve in the bottom panels shows the joint best fit to the data with the generalized NFW functional form (Nagai et al. 2007, see Eq. (9)). In all these plots the dotted and dashed-dotted lines represent the result of Planck Collaboration Int. V (2013) and Arnaud et al. (2010) respectively. The shadow areas represent the envelope obtained by Planck Collaboration Int. V (2013) and the Early release SZ sample (XMM-ESZ; Planck Collaboration VIII 2011).

In the text
thumbnail Fig. 5.

Same as Fig. 3 for the projected temperature profiles rescaled by the self-similar quantity T500 (Eq. (10)). The filled circles show the measurements of the X-ray spectroscopic temperature (see Sect. 2.4.2 ), whereas the filled squares indicate the data points obtained by combining the SZ pressure with the gas density, projected along the line of sight assuming the spectroscopic-like scaling of Mazzotta et al. (2004). The solid red curves in the bottom panels show the best fit to the joint dataset with the functional form introduced by Vikhlinin et al. (2006; see Eq. (11)).

In the text
thumbnail Fig. 6.

Same as Fig. 3 for the entropy profiles rescaled by the self-similar quantity K500 (Eq. (12)). The filled circles show the measurements obtained from the X-ray spectroscopic temperature as , whereas the filled squares indicate the data points obtained by combining the SZ pressure with the gas density as . The solid red curves in the bottom panels indicate the best fit to the entire population with the functional form presented in Eq. (13), whereas the solid black line shows the prediction of pure gravitational collapse (Voit et al. 2005).

In the text
thumbnail Fig. 7.

Measured intrinsic scatter of all our thermodynamic quantities, density (red), pressure (green), temperature (blue), and entropy (black). The data points indicate the results of piecewise power law fits on several radial ranges, whereas the dashed lines and shaded areas show the intrinsic scatter described as a log-parabola (Eq. (6)) around the best fitting functional forms.

In the text
thumbnail Fig. 8.

Thermodynamic quantities of the X-COP clusters dividing the clusters into cool-core (CC, in blue) and non-cool-core (NCC, in red) populations, compared with the entire population (ALL, in green). Top left panel: density profiles fitted using the functional form presented in Eq. (7), overplotted on the data and on the “universal” density profile of Eckert et al. (2012, pink shaded area). Top right panel: pressure profiles fitted using the functional form presented in Eq. (9), overplotted on the data and compared to the Planck results (Planck Collaboration Int. V 2013), dotted black line) and the universal pressure profile of Arnaud et al. (2010, dash-dotted black line). Bottom left panel: temperature profiles fitted using the functional form presented in Eq. (7) and overplotted on the data. Bottom right panel: entropy profiles fitted using the functional form presented in Eq. (13), overplotted on the data and on the gravitational collapse predictions (Voit et al. 2005, solid black line).

In the text
thumbnail Fig. A.1.

Posterior distribution of the parameter ηSZ ≈  PSZ/PX. The shaded blue area indicates the region containing 68% of the posterior distribution.

In the text
thumbnail Fig. B.1.

Left panel: residuals of the comparison of the two different methods we used to estimate the SZ pressure profile. A remarkable good agreement within the statistical uncertainties is reached at all radii, especially excluding the first 3 MILCA points which are the most affected by the Planck PSF. Right panel: distribution of the residuals compared with the statistical prediction of a set of residuals: a Gaussian centered in zero and width one (red line).

In the text
thumbnail Fig. C.1.

Scatter in the various thermodynamic profiles induced by uncertainties in the mass calibration. The shaded areas show the range of intrinsic scatter obtained from 1000 simulations.

In the text
thumbnail Fig. D.1.

Same as Fig. 8, but using the piecewise fit instead of global functional form. The cluster population of X-COP is divided into cool-core (CC, in blue) and non-cool-core clusters (NCC, in red), compared with the entire population (ALL, in green).

In the text
thumbnail Fig. D.2.

Scatter of our thermodynamic quantities dividing our cluster sample in CC and NCC.

In the text
thumbnail Fig. E.1.

Density – marginal probability Vikhlinin+06. Parameter distribution for the best fit on the density of all clusters using the Vikhlinin et al. (2006) functional form, Eq. (7). The priors on the parameters are shown in Table 3.

In the text
thumbnail Fig. E.2.

Pressure – marginal probability Nagai+07 – α fixed. Parameter distribution for the best fit on the pressure of all clusters using the Nagai et al. (2007) gNFW functional form, Eq. (9), fixing the intermediate slope α. The priors on the parameters are shown in Table 3.

In the text
thumbnail Fig. E.3.

Temperature – marginal probability Vikhlinin+06. Parameter distribution for the best fit on the temperature of all clusters using the Vikhlinin et al. (2006) functional form, Eq. (11). The priors on the parameters are shown in Table 3.

In the text
thumbnail Fig. E.4.

Entropy – marginal probability powerlaw + const. Parameter distribution for the best fit on the entropy of all clusters using a power law plus constant functional form, Eq. (13). The priors on the parameters are shown in Table 3.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.