© ESO, 2013

1. Introduction

In a pure hierarchical gravitational collapse scenario in the concordance ΛCDM cosmology, concentrations of matter (“halos”) are fully characterised by their redshift and mean matter density, which in turn are related to the power spectrum of initial density fluctuations (Peebles 1980). The scale-free dark matter (DM) collapse drives the evolution of halo concentration across cosmic time (see e.g., recent work by Bhattacharya et al. 2011; Gao et al. 2012), and the ensuing similarity yields both a universal dark matter distribution and simple global scaling relations that should describe the entire halo population (Bertschinger 1985; Kaiser et al. 1995; Navarro et al. 1995; Evrard et al. 1996; Navarro et al. 1997; Voit 2005; Arnaud et al. 2005). However, the observable properties of clusters are determined by the visible baryonic component, which is subject to more complex physical processes related to galaxy formation and feedback. As the main baryonic reservoir in massive halos, the hot gas in the intra-cluster medium (ICM) is the natural target for studying the physical processes at play and their link to the underlying cluster DM content. Modelling and understanding the baryon physics and disentangling the effect of various feedback processes is a very challenging task (see e.g., Borgani & Kravtsov 2011, for a review), underlining the need for appropriate observational constraints.

The ICM attains X-ray emitting temperatures due to gravitational heating within the halo potential well. X-ray emission is proportional to the square of the gas density, thus it probes the denser regions of the hot gas. The thermal Sunyaev-Zeldovich (SZ) effect (Sunyaev & Zeldovich 1972), due to inverse Compton scattering of CMB photons by the ICM, is proportional to the thermal gas pressure integrated along the line of sight. It is sensitive to density and/or temperature variations, such as shocks and compression. These two independent observational probes are thus complementary and allow us to further constrain the physics of the ICM.

Recent X-ray observations based primarily on representative samples have returned a consistent picture of the scaling and structural properties of halos, from high mass clusters down to the low mass group regime (see, e.g., Böhringer et al. 2007; Vikhlinin et al. 2006; Croston et al. 2008; Pratt et al. 2009; Arnaud et al. 2010; Sun et al. 2011; Sun 2012). In parallel, SZ observations with instruments such as the Sunyaev-Zeldovich Array (SZA, Muchovej et al. 2007), the South Pole Telescope (SPT, Carlstrom et al. 2011), the Atacama Cosmology Telescope (ACT, Swetz et al. 2011), and Planck (Tauber et al. 2010)¹ have recently started to deliver on the promise of SZ observations for cluster studies.

Building on earlier works on smaller cluster samples (Benson et al. 2004; Andersson et al. 2011), recent results from Planck have underlined the consistency between the X-ray and SZ view of the ICM within R₅₀₀². These constraints were achieved using three different approaches. The first, detailed in Planck Collaboration (2011e), involved bin-averaging of an X-ray selected sample from the Meta-Catalogue of X-ray detected Clusters of galaxies (MCXC, Piffaretti et al. 2011). The second, described in Planck Collaboration (2011f), concerned the comparison of SZ measurements of 62 local clusters from the Early Release Compact Source Catalogue (ERCSC, Planck Collaboration 2011b,c) with good quality XMM-Newton archive data. Finally Planck Collaboration (2011d) examined the scaling properties of 21 newly-detected Planck clusters confirmed by XMM-Newton. These studies have provided well-constrained scaling relations (e.g., Y₅₀₀ − M₅₀₀, Y₅₀₀ − L_X,500) in the local Universe, including the first measurement of their intrinsic scatter, to be used as a reference for future evolution and cosmology studies. Other dedicated investigations of smaller samples (e.g., Andersson et al. 2011; Sifon et al. 2012) have reached conclusions similar to Planck early results regarding the agreement between the X-ray and SZ view of the ICM.

Beyond global properties, similarity of shape in cluster radial quantities has long been shown in X-ray observations (e.g., Arnaud et al. 2001; Pointecouteau et al. 2005; Vikhlinin et al. 2006; Croston et al. 2008; Arnaud et al. 2010), in the optical (e.g., Rines & Diaferio 2006; Wojtak & Lokas 2010), and with weak-lensing (e.g., Postman et al. 2012). For SZ studies, the radial distribution of the thermal ICM pressure in massive halos is of particular interest. Recent X-ray observations of REXCESS, a representative sample of nearby objects, have shown that the ICM pressure distribution, when scaled appropriately, follows an approximately “universal” shape (Arnaud et al. 2010, A10 hereafter) up to a radius of R₅₀₀. The small differences in shape for the pressure profiles between relaxed and unrelaxed clusters, especially in the central regions, does not seem to have an impact on the integrated thermal content within R₅₀₀, or on other global properties, as noted by Planck Collaboration (2011f).

A10 also compared their X-ray pressure profiles to those predicted from numerical simulations. Although limited to R₅₀₀, good agreement was seen between radial pressure profile observations and predictions from various numerical simulations within this radius. However, there are presently few observational constraints beyond R₅₀₀, and consequently in this region the shape of the “universal pressure profile” was extrapolated according to the predictions from numerical simulations (Borgani et al. 2004; Nagai et al. 2007b; Piffaretti & Valdarnini 2008).

The thermodynamical state of the gas beyond R₅₀₀ bears the signature of the complex physics taking place in the outer parts of the clusters. Characterisation of the gas in cluster outskirts is also necessary to unveil the level of thermal pressure in the cluster periphery, thus constraining their dynamical state. These constraints are crucial for our understanding of the formation and evolution of massive halos. In the cluster surroundings, accretion along filaments is three-dimensional and non-spherical (Tozzi & Norman 2001) and the initial conditions of the accretion shock are driven by the thermodynamical state of the (pre-shocked) in-falling material (Voit et al. 2002, 2003). The physical origin of any possible pre-heating is still unclear, but it may smooth the continuous accretion of substructures. Thus the entropy production at the accretion shock may be boosted (Voit & Ponman 2003; Borgani et al. 2005). Feedback may also improve the degree of thermalisation of the ICM after the shock (with a residual kinetic energy of about 10% beyond R₅₀₀, Kay et al. 2004).

Observations with X-ray telescopes have only recently started to provide insight into the physical properties of the gas in the cluster outskirts beyond R₅₀₀ (George et al. 2009; Reiprich et al. 2009; Urban et al. 2011). The steep decline of the X-ray emission with cluster-centric radius makes such observations extremely challenging with current instruments. A deep X-ray study of the Perseus cluster with Suzaku by Simionescu et al. (2011) indicated that the gas mass fraction in the cluster outskirts was overestimated, possibly due to the effect of gas clumping and the non-virialised state of the ICM in these regions, as suggested in some numerical simulations (Nagai et al. 2007b).

The different sensitivity of the SZ effect to the radial ICM distribution means that SZ observations have the potential to contribute greatly to the discussion on cluster outskirts. The radial pressure distribution of the first SZ cluster samples have recently been presented based on observations by SPT (Plagge et al. 2010), ACT (Sehgal et al. 2011) and SZA/CARMA (Bonamente et al. 2012). These studies confirmed that the ICM properties as seen by SZ and X-ray observations are consistent at least out to R₅₀₀. Beyond this radius the SZ effect offers the interesting possibility to further constrain the thermal pressure support. Planck is highly competitive in this regard. It is the only SZ experiment with a full sky coverage, able to map even nearby clusters to their outermost radii and offering the possibility of an in-depth statistical study through the combination of many observations. In this paper we present constraints on the thermal pressure support derived using SZ observations from the Planck survey. Following our previous methodology for scaling relations, we investigate these issues from a statistical point of view, working with data from 62 local clusters selected from the Planck ESZ sample (Planck Collaboration 2011c), for which there are good quality XMM-Newton archival data.

The paper is organised as follows. In the next section we describe the Planck mission and data, together with the sample of galaxy clusters used in this study. In Sect. 3 we recall the basic formalism on the SZ effect and the parametrisaton of the cluster pressure profile we use throughout the paper. Section 4 is devoted to a detailed description of the processing involved in the reconstruction of the SZ and pressure profiles from Planck data. In Sects. 5 and 6 we present the stacked profile of our sample and its best analytical representation. We discuss the comparison with other observational and theoretical constraints in Sect. 7, before presenting our conclusions in Sect. 8.

Throughout the paper we adopt a ΛCDM cosmology with H₀ = 70  km s^-1 Mpc^-1, Ω_M = 0.3 and Ω_Λ = 0.7 . The quantity E(z) is the ratio of the Hubble parameter at redshift z to its present value, H₀, i.e., E(z)² = Ω_M(1 + z)³ + Ω_Λ.

2. Data

2.1. Planck data

Planck (Tauber et al. 2010; Planck Collaboration 2011a) is the third generation space mission to measure the anisotropy of the CMB. Planck observes the sky in nine frequency bands covering 30–857 GHz, with high sensitivity and angular resolution from 31′ to 5′. The Low Frequency Instrument (LFI; Mandolesi et al. 2010; Bersanelli et al. 2010; Mennella et al. 2011) covers the 30, 44, and 70 GHz bands with amplifiers cooled to 20 K. The High Frequency Instrument (HFI; Lamarre et al. 2010; Planck HFI Core Team 2011a) covers the 100, 143, 217, 353, 545, and 857 GHz bands with bolometers cooled to 0.1 K. Early astrophysics results, based on data taken between 13 August, 2009 and 7 June, 2010 (Planck HFI Core Team 2011b; Zacchei et al. 2011), are given in Planck Collaboration 2011c,d,e,f,g. Intermediate astrophysics results are now based on data taken between 13 August, 2009 and 27 November, 2010.

We used the full sky maps in the nine Planck frequency bands provided in healpix (Górski et al. 2005) N_side = 2048 full resolution. An error map associated to each frequency band is obtained from the difference of the first half and second half of the Planck rings for a given position of the satellite. The resulting jack-knife maps are basically free from astrophysical emission, whilst being a good representation of the statistical instrumental noise and systematic error. We adopted a circular Gaussian as the beam pattern for each frequency, as described in Planck HFI Core Team (2011b) and Zacchei et al. (2011). Uncertainties in flux measurements due to beam corrections, map calibrations and uncertainties in bandpasses are expected to be small, as discussed extensively in Planck Collaboration (2011c), Planck Collaboration (2011d) and Planck Collaboration (2011e).

2.2. Cluster sample

The first all-sky coverage by the Planck satellite (Planck Collaboration 2011a) has allowed the detection of dozens of clusters via their SZ signature on the CMB. The Early Release SZ sample (ESZ, Planck Collaboration 2011c) comprises 189 clusters of galaxies characterised over the sub-millimetre to centimetre wavelength range. The sample consists of SZ clusters and candidates detected with signal-to-noise ratios (S/N) spanning from 6 to 29 in the first all-sky survey. The sample was thoroughly validated by Planck internal quality assessment, external X-ray and optical data cross-correlations, and a multi-frequency follow-up programme for confirmation. The ESZ-cluster sample spans over a decade in mass, from 0.9 to 15 × 10¹⁴ M_⊙, which is essentially the full galaxy cluster mass range.

In the following, we focus on a sub-sample of 62 clusters from the ESZ catalogue. This sub-sample is defined and extensively described in Planck Collaboration (2011f), where it was used to calibrate scaling relations between the SZ and X-ray cluster properties. All 62 clusters were already known in X-rays (found in the MCXC, Piffaretti et al. 2011), and have good XMM-Newton archive data, allowing for a high quality X-ray data analysis. Masses and radii for the sample were estimated using the M₅₀₀–Y_X relation of A10 (see also Pratt et al. 2010), assuming standard evolution (see Eq. (9)).

The radius R₅₀₀ was calculated iteratively as described in Kravtsov et al. (2006). While this sample is neither representative nor complete, it is the largest, highest-quality SZ-X-ray data set currently available. The majority of objects lie at a redshift lower than 0.3 (and all have z < 0.5) and cover approximately a decade in mass. In angular size, i.e., θ₅₀₀, they range between 3.7 and 22.8 arcmin, with a median value of 7.6 arcmin. These extended, SZ-bright clusters are thus ideal targets to investigate the spatial distribution of ICM thermal pressure support in clusters of galaxies by means of the spatial distribution of their SZ signal.

3. Basic formalism

3.1. The Sunyaev-Zeldovich effect

The inverse Compton effect for thermalised electrons under a blackbody radiation field is treated in the Kompaneets (1957) equation. It was studied and characterised for the case of thermal electrons in clusters of galaxies by Sunyaev & Zeldovich (1970, 1972). The ensuing spectral distortion of the CMB spectrum in the direction of clusters is named the Sunyaev-Zeldovich (SZ) effect.

In the following, we leave aside the kinetic SZ effect, a doppler effect, resulting from the cluster peculiar motion within the comoving reference frame of the Hubble flow. Also, the relativistic corrections on the SZ spectrum, of the order of kT_e/m_ec², are not relevant for our study and are thus neglected (see, e.g., Pointecouteau et al. 1998; Challinor & Lasenby 1998; Sazonov & Sunyaev 1998).

The intensity of the SZ effect is characterised by the dimensionless Comptonisation parameter y, the product of the average fractional energy transferred per collision and the average number of collisions: (1)where σ_T is the Thomson cross-section, m_e the mass of the electron and c the speed of light. P is the pressure produced by the plasma of thermal electrons along the line of sight.

Assuming the clusters are spherically symmetric, we can express the profile in the Comptonisation parameter as a geometrical projection of the spherical pressure profile along the axis of a cylinder: (2)As y is dimensionless we have y(θ) ≡ y(r). Actual SZ measurements derive from the convolution of the y-profile on the sky with the instrument spatial response, f_PSF: (3)If the SZ brightness, which is proportional to y, is indeed independent of redshift, the SZ flux is proportional to the integrated Comptonisation parameter and thus depends on the source distance via (4)where D_A is the angular diameter distance.

In the following, the observed value of the integrated Comptonisation parameter is given: in units of Mpc² when expressed in the source intrinsic reference frame, i.e., Y(R); and in arcmin² when expressed in the Planck satellite reference frame, i.e. Y(θ). Analogously we can define the total integrated pressure within a sphere of radius r, and express it in Y units. We define this “pseudo” or “spherical” integrated Comptonisation parameter as (5)

3.2. The scaled pressure profile

As mentioned in Sect. 1, in a hierarchical scenario of structure formation the halo population is self-similar in scale and structure. Profiles of physical quantities are universal once scaled according to their radius and reference quantities defined at a given density contrast, δ. Within this self-similar framework, we adopt throughout the paper the value of δ = 500. The scaled pressure profile thus reads (6)We adopt for the pressure profile the analytical formulation given by Nagai et al. (2007a) for the generalised Navarro-Frenk White (Navarro et al. 1997; Nagai et al. 2007a, GNFW) profile (7)where x = r/R₅₀₀, and the model is defined by the following parameters: P₀, normalisation; c₅₀₀, concentration parameter defined at the characteristic radius R₅₀₀; and the slopes in the central (x ≪ 1/c₅₀₀), intermediate (x ~ 1/c₅₀₀) and outer regions (x ≫ 1/c₅₀₀), given by γ, α and β, respectively.

A10 have fitted this analytical profile to a combination of observed pressure profiles derived for the REXCESS sample (Böhringer et al. 2007) from high quality XMM-Newton data, together with three sets of predicted profiles from numerical simulations of structure formation implementing DM and baryon physics (i.e., radiative cooling and recipes for feedback). Measurements of the gas mass (integrating the density profile) and the X-ray spectroscopic temperature lead to the quantity Y_X,500 = M_gas,500T_X,500. This parameter links to the actual integrated parameter as seen from the SZ, i.e., Y₅₀₀, as (8)where μ_e = 1.148, is the mean molecular weight of electrons for a 0.3 solar abundance plasma and m_p is the proton mass. Here A_XSZ = 0.924 ± 0.004 from Eq. (19) of A10 derived from the REXCESS sample, and A_XSZ = 0.95 ± 0.04 from the best fit SZ scaling relation between Y₅₀₀ and Y_X,500 from Planck Collaboration (2011f) using the same sample as in this work. For consistency purposes we use this latest value in the following. The predicted and measured values for Y₅₀₀ are compatible within the  ± 1σ limits.

Y_X,500 is a good proxy for the cluster total mass (Kravtsov et al. 2006; Nagai et al. 2007a; Arnaud et al. 2007). The non-standard scaling relation fitted against the REXCESS data is provided by their Eq. (2): (9)In turn the characteristic pressure P₅₀₀ scales with the cluster total mass, reflecting the mass variation expected in the standard self-similar model, purely based on gravitation (Eq. (5) in A10): (10)The M − Y_X relation given in Eq. (9) deviates from the standard self-similar case (which has a slope of 3/5). From the definition of P₅₀₀, any deviation from the standard self-similar scaling will appear as a variation of the scaled pressure profiles. As shown by A10 in their Eq. (9), this variation can be expressed as a function of the total mass. At δ = 500, this mass dependence is almost constant with radius, and it can be approximated by (11)

3.3. Conversions and normalisations

In the following, we measure SZ profiles from the Planck data (Sect. 4.2). All through this work, we compare these observed profiles with the predictions that arise from the X-ray constraints. For all clusters we assume the universal pressure profile shape as derived from A10. We parametrise it according to the quantities derived from the X-ray analysis R₅₀₀ and Y_X,500 (see previous section). We can then derive predicted pressure profiles (see Eq. (7)), and Comptonisation parameter profiles (Eqs. (2) and (4)). The later convolved by the instrument beam (Eq. (3)) is directly comparable to our observed Planck profiles.

We apply a statistical approach to the observed SZ profile and to the pressure profile, averaging individual profiles in our sample once scaled. Profiles are scaled in radius according to R₅₀₀. The observed SZ profiles (or predicted SZ profiles from X-ray constraints) are dimensionless, and we therefore normalise them by the quantity (12)In the following we used the value of Y₅₀₀ as defined in Eq. (5).

To translate our measured SZ profiles, we deconvolve and deproject them for each cluster following Eqs. (2) and (3). In practice this operation is performed as described in Sect. 6.1. The resulting profiles are converted to pressure such that (13)

4. Reconstruction of the SZ profile

4.1. SZ signal extraction methods

The thermal SZ maps were recovered from a combination of Planck channels, making use of either real or Fourier space methods. We applied three methods based on internal linear combination (ILC) algorithms to the Planck data: (1) the Modified Internal Linear Combination Algorithm (MILCA, Hurier et al. 2010); (2) Needlet Internal Linear Combination (NILC, Delabrouille et al. 2009); and (3) the Generalized Morphological Component Analysis (GMCA Bobin et al. 2008). A description of each method is given in Appendix A.

To optimise the reconstructed SZ map we have to take into account that: (i) the brightness of the SZ effect (as an increment or decrement) is maximum in the sub-millimetre to millimetre range; (ii) the typical angular size of the cluster is a few arcmin; and (iii) the final resolution of the reconstructed SZ map is determined by the lowest resolution of the combined frequency maps. We therefore restricted ourselves to the use of the six Planck-HFI channels in the SZ map reconstruction process. The final resolution of the reconstructed SZ maps is that of the lowest-resolution 100 GHz channel, or 10 arcmin. SZ maps are in units of the (dimensionless) Comptonisation parameter. Integrated Comptonisation parameters are expressed in arcmin² for the observed values and in Mpc² for the intrinsic value in the source reference frame. The associated SZ noise maps are built in a fully consistent way from the individual frequency error maps (see Sect. 2).

4.2. SZ profile computation

4.2.1. Profile extraction

Individual cluster profiles were computed from the reconstructed all-sky SZ maps. Uncertainties were obtained from the all-sky reconstructed SZ error maps, which were derived by applying the same reconstruction methods to the frequency error maps (see Sect. 2). For instance, in the case of a linear method (e.g., ILC based), the coefficients used for the linear combination of the frequency maps were propagated to the combination of the frequency error maps.

We extracted a square patch of side 20 × θ₅₀₀ around each cluster position from the all-sky SZ map. Patches were projected from healpix (Górski et al. 2005) to a tangential projection. For each patch, the pixel size was adapted so that it was constant in scaled units of θ₅₀₀ over the full cluster sample. This unavoidably leads to an oversampling of the pixels or, more precisely, to a redundancy of the original all-sky map pixels in the reprojected patch. We produced an associated patch that tracks this redundancy. We also extracted equivalent patches from the associated all-sky SZ error map, and from the all-sky variance map.

We then computed a profile from each SZ map patch. These were calculated on a regular radial grid with bins of width Δr/R₅₀₀ = 0.25, allowing us to sample each cluster profile with four points within R₅₀₀. The y value of a bin was defined as the mean of the values of the pixels falling in each annulus. We subtracted a background offset from the maps prior to the profile computation. The offset value was estimated from the surrounding region of each cluster where r > 7 R₅₀₀. The uncertainty associated with this baseline offset subtraction was propagated into the uncertainty of each bin of the radial SZ profile.

4.3. Expected SZ profile

We used the GNFW pressure profile shape characterised by A10 from the REXCESS sample and numerical simulations to generate a pressure model for each cluster of our sample. We kept the A10 best fit values of c₅₀₀ = 1.18 and the three slopes, α = 1.05, β = 5.49, γ = 0.31. We derived the normalisation P₀ from the observed value of Y_X,500 (see Eq. (5)), using a conversion of 0.95 between Y_X and Y_SZ as described in Sect. 3.2. Finally R₅₀₀ was fixed to the values derived from the XMM-Newton analysis. Both Y_X and R₅₀₀ are reported in Planck Collaboration (2011f). We computed the projection matrix and the PSF redistribution matrix (for a Gaussian beam of 10 arcmin FWHM) as expressed in Eqs. (2) and (3), respectively, to multiply each pressure profile model. The derived y profile model is directly comparable to the Planck measurements. The average profile model across the sample was derived similarly to the observed stacked profile (see above), and is used to compare to the observed stacked profile.

4.3.1. Accounting for correlation between points

A certain level of correlation is introduced between the points of a radial profile derived in this manner. We account for this correlation in the covariance matrix of each profile, which is computed as follows.

For each map patch, we masked radii θ < 7 θ₅₀₀ centred on the cluster and computed the power spectrum of the noise on the remaining area. The use of the surrounding regions of the cluster to characterise the noise properties allows us to account for the effect of astrophysical contamination in the cluster vicinity, as well as the instrument noise and systematics. These sources of contamination were, by construction, excluded from the all-sky frequency error maps (see Sect. 2). The power spectrum drawn from the SZ error map of each patch is therefore, systematically, slightly lower than the one computed over the cluster surroundings.

We then simulated m = 500 realisations of the noise patch based on this power spectrum, accounting for the variance map and sky pixel redundancy and assuming an inhomogeneous correlated Gaussian noise. We then extracted a noise profile from each realisation, reproducing the baseline background subtraction used for the observed profile. The covariance matrix was built from all the simulated noise profiles, i.e., , where P_n is an npoints × m matrix of simulated noise profiles. Considering two points as correlated when their correlation coefficient is larger than 0.3, the typical level of point-to-point correlation in our profiles is about 20% (16%, 21% and 28% for MILCA, NILC and GMCA, respectively).

Fig. 1 Average scaled SZ profiles for each of the nine Planck frequency bands (decreasing frequency from left to right and from top to bottom). The points within each individual profile are correlated at about the 20% level (see Sect. 4.3.1). At each frequency, the stacked radial profile is obtained from the average of the 62 individual cluster profiles scaled by R500 and Φi/ ⟨ Φi ⟩ , in the x- and y-axis directions, respectively (see Sect. 5.1). The light-red shaded area in each panel indicates the statistical uncertainty about the average observed profile. The solid black line corresponds to the prediction of the universal pressure profile (A10) for our cluster sample in each band, and the two dashed lines indicate its associated dispersion.

4.3.2. Stacking procedure

Given that all profiles are computed on the same grid in scaled radius, they can easily be stacked together. Each profile y_i was rescaled by R_500,i and Φ_i (as given in Eq. (12)), respectively in the X and Y axis directions. The stacked profile and associated covariance matrix are given by (14)For the computation of Φ, we used Y_sph,500 in Mpc² as given in Eq. (5). For simulated data (Appendix B), Φ was derived from the Planck M₅₀₀–Y₅₀₀ relation (Planck Collaboration 2011d); for observed data (Sect. 5.3), Φ was measured directly from the Planck data.

When for a given bin, values for all clusters are dominated by the signal, we assumed a log-normal distribution of their scatter and stacked them in logarithmic space. The logarithm of the stacked value and associated error were then translated back into linear space. In other terms, for bin j, if all the clusters satisfied the condition , we stacked all the measurements in the logarithmic space; otherwise, we did it in the linear space.

We computed the statistical and weighted (i.e., 1/σ²) average profiles, as well as the median profile of the sample, and checked their consistency. Over a radial range of 3 × R₅₀₀, the average of the relative error of the weighted mean and median profiles with respect to the statistical average profile is smaller than 10% and 5%, respectively. Based on this agreement, we used the statistical average to compute the stacked profile throughout this study.

Of the 62 clusters in our ESZ-XMM sample, only two are spatially close and thus are potentially physically connected: A3528 and A3532. These two clusters are members of the Shapley supercluster. Each of those clusters was masked over an area of radius 3 R₅₀₀ when processing the other. Thus, we consider that their profiles and covariance matrices are independent.

At the end of our data processing, the stacked covariance matrix encompasses the statistical errors due to instrumental noise, astrophysical fluctuations at the cluster location and systematic effects (e.g., instrumental and arising from data processing). For the purposes of our study we have also computed the dispersion across our sample for each position within our profiles. Both uncertainties are propagated throughout the analysis, and are cross-compared in the following sections. In Figs. 3, 4 and 6 the error bars shown on the Planck data points are purely statistical and correspond to the square root of the diagonal elements of the covariance matrix (i.e., Eq. (14)). The rms scatter of individual profiles around the mean is indicated by a coloured band.

5. Planck galaxy cluster SZ profile

5.1. Frequency stacked profiles

We first looked for the raw cluster signature in each of the nine Planck frequencies. We extracted a brightness profile at each frequency from the raw Planck maps for each of the 62 clusters in the ESZ-XMM sample. These were rescaled in units of R₅₀₀ in the x-axis direction and in units of Φ_i/ ⟨ Φ_i ⟩  in the y-axis direction (see Sect. 3.3), and then stacked. The resulting average brightness profiles are shown in Fig. 1, where the shaded area in each panel depicts the statistical error about the mean flux value in each bin.

The cluster signal appears clearly in most of the nine frequencies, following the SZ thermal spectral signature. The signal is positive at 545 and 353 GHz, compatible with zero at 217 GHz, and then negative down to 30 GHz. The 62 clusters of our sample are strong SZ sources, as expected for objects detected at S/N > 6 (Planck Collaboration 2011c). The stacking procedure mostly averages out the effect of foreground and background contaminants; however, they contribute to the final dispersion in the stacked profiles. Owing to their Gaussian nature, CMB fluctuations are more easily washed out in the stacking procedure than the Galactic and extragalactic dust emission in the high frequency range, and Galactic and point source radio emission (i.e., free-free and synchrotron) in the lower frequency range.

We also plotted for each frequency the expected profile computed as described in Sect. 4.3. This comparison shows a first order agreement between the Planck measurements and the expected SZ profiles derived from X-ray constraints assuming a GNFW pressure profile shape.

5.2. Stacked SZ image

Figure 2 presents the stacked average image of the 62 clusters of our ESZ-XMM sample, plotted on a logarithmic scale. Before averaging, each individual SZ map was normalised identically to the profiles (see Eq. (12)) and randomly rotated by 0°, 90°, 180° or 270°. The stacked image was then renormalised by  ⟨ Φ_i ⟩  in order to be expressed in y units. The image shows that the SZ signal is detected out to  ~ 3 × R₅₀₀, marked by the white circle in the figure. The extent of existing observational constraints on the gas distribution in clusters in SZ and X-rays, i.e., R₅₀₀, is marked by the black circle. We also generated a jackknife-type map by averaging the scaled individual SZ maps, m_i, as  ∑ ( − 1)ⁱm_i. The rms values of both the stacked SZ and on-off maps beyond 5 × R₅₀₀ are compatible: 4.6 × 10^-7 and 4.8 × 10^-7, and, in turn, compatible with the rms value below this radius on the jackknife stacked map: 3.5 × 10^-7. This check demonstrates that our average map is not strongly affected by residuals or biases. The radial profile computed from the stacked map is fully compatible with the stacked profile discussed below.

Fig. 2 Stacked scaled image (size is 10 R500 on a side) using a logarithmic stretch in y. Individual maps are rescaled by Φi before averaging, and then multiplied by  ⟨ Φi ⟩ . The black and white circles mark the loci of 1 and 3 × R500, respectively.

5.3. Observed stacked SZ profile

We derived the observed stacked Planck SZ profiles for the ESZ-XMM sample as described in Sect. 4.2. Y₅₀₀ is obtained from the algorithms used for blind detection of SZ clusters in the Planck survey (Planck Collaboration 2011c), namely the Powell Snakes (PwS, Carvalho et al. 2009, 2011) and multi-frequency matched filter (MMF3, Melin et al. 2006) algorithms. In both cases the algorithms implement a universal pressure profile shape (A10) with the position of the cluster fixed to the XMM-Newton coordinates and the size of the universal pressure profile fixed to R₅₀₀. The fluxes from both methods (i.e., PwS and MMF3) are consistent over the whole sample. The median value for the ratio PwS/MMF3 is 0.96 ± 0.05 (see Appendix B.2). To further validate the above, and as a consistency check with previous Planck results, we fitted each individual SZ profile with a projected, PSF-convolved universal pressure profile. We fixed R₅₀₀ to the best fitting X-ray value from Planck Collaboration (2011f) and only fitted the normalisation, Y₅₀₀. Given the agreement, we used the MMF3 values to compute Φ for each cluster (Eq. (12)).

Appendix B presents a detailed investigation of the convergence between stacked profiles derived from each of the three different methods of SZ signal reconstruction. We find that all three methods agree remarkably well over the entire radial range, both for simulated and observed SZ profiles. The three SZ reconstruction methods lead to profiles fully compatible with each other. Across the range of radii over which the profiles have been computed (i.e., [0–10] × R₅₀₀), accounting for the correlated errors of each profile, the reduced χ² of NILC and GMCA with respect to MILCA are 0.48 and 0.62, respectively³.

The left panel of Fig. 3 shows the stacked MILCA profile for the 62 clusters in the ESZ-XMM sample compared to the individual profiles (see Sect. 4.2). The SZ signal is statistically detected over more than two decades in intensity, and out to a remarkably large radius of  ~3 R₅₀₀, reaching far into the cluster outskirts. Assuming that the virial theorem can be applied (which at such large radius might be breaking down), the outer radius of our statistical detection corresponds to a density contrast of δ  ~  50. More conservatively, we can assume that we are statistically probing the average SZ and pressure distribution down to regions of density contrast of δ ~ 50–100. The dispersion about the mean profile dominates the statistical uncertainties. It is minimal (by construction) at  ~20% around R₅₀₀, but increases towards the centre and the outskirts to  ~50% and  ~65% at 0.3 and 2 × R₅₀₀, respectively. At the most external radius of our statistical detection, i.e.,  ~3 × R₅₀₀, the dispersion of the individual profiles about the mean is more than 100%, as at these large radii, the individual SZ profiles are fully dominated by noise.

The right panel of Fig. 3 shows the correlation coefficient matrix for the stacked profiles presented on the left panel, and illustrates the degree of correlation between the points in the profiles.

5.4. Comparison with expectations from universal profile

The stacked model used for comparison was computed as described in Sect. 4.3. When considering the statistical errors only, the measured stacked profile is significantly above the model. Taking into account the error on the model and the correlated errors between points of the reconstructed SZ profile, we obtained a reduced χ² value of 3.53 within 3 R₅₀₀. If we omit the error on the model, this value becomes 7.48 As our tests on simulations show a very good agreement between the input model and the output reconstructed profiles (see Appendix B), the difference observed here between the measured SZ profiles and the predicted model from X-ray constraints is not an artefact of the method, but an intrinsic difference. This difference is significant at a 2–3σ level from around R₅₀₀ out to 3 R₅₀₀.

6. The galaxy cluster pressure profile from Planck

Fig. 3 Left: individual profiles (grey lines) and scaled stacked profile (red points) for our sample of 62 clusters. The light-red shaded area marks the dispersion about the average stacked profile, with its upper and lower limits highlighted by the two dashed red lines. The points within each individual profile are correlated at about the 20% level (see Sect. 4.3.1). The error bars on the Planck data points are purely statistical and correspond to the square root of the diagonal elements of the covariance matrix. The solid black line (labeled in the legend as “Expected from PA10”) is the average stacked profile obtained from the expected SZ individual profiles drawn from the universal pressure profile of A10, parameterised according R500 and Y500 derived from the XMM-Newton data analysis presented in Planck Collaboration (2011f). The two dotted black lines indicate the associated dispersion about this model profile. The bottom panel shows the value of χ at each point with respect to the expectations from the universal profile taking into account the measured errors. Right: matrix of the correlation coefficient for the observed radial stacked profiles. Both panels are illustrated here for the MILCA results.

Fig. 4 Left: Planck pressure profile obtained from the average of the individual pressure profiles across our sample (red points) shown together with the stacked pressure profile derived from the XMM-Newton data for the same sample (purple points). The dispersions about the SZ and X-ray profiles are depicted, respectively, by the red and purple shaded areas. Our best fit GNFW profile is shown as a solid black line and that of A10 as a solid green line. The lower panel shows the χ profile of these two best models taking into account the statistical errors and dispersion about the observed profile. Right: stacked profile from Planck and XMM-Newton for the sub-samples of cool-core (open circles) and non cool-core (open stars) clusters within the ESZ-XMM sample. The shaded areas are identical to the one shown in the left panel. Our best fit models for each sub-sample are shown as black solid and dashed lines for the cool-core and non cool-core clusters, respectively (see Table 1). The best A10 fit for cool-core and non cool-core clusters are shown as green solid and dashed lines, respectively. The lower panel gives the χ profiles of our best fit models with respect to the measured profiles and associated errors (including dispersion). On both panels the Planck data points are correlated at about the 20% level (see Sect. 4.3.1). As for the stacked SZ profile, the error bars on the Planck points are purely statistical and correspond to the square root of the diagonal elements of the covariance matrix. For the XMM-Newton points they correspond to the statistical error on mean.

6.1. PSF deconvolution and deprojection

The pressure profile is derived for each cluster in our sample by applying a deconvolution and deprojection to the observed SZ profiles. As already mentioned, we assume spherical symmetry. The convolution by the instrumental beam and the geometrical projection of a spherical pressure profile into a Comptonisation parameter profile are expressed in Eqs. (3) and (2), respectively. We have applied a deconvolution plus deprojection algorithm adapted from the method described by Croston et al. (2006). This method allows us to perform a straightforward real space deconvolution and deprojection using a regularisation procedure originally applied to X-ray surface brightness profiles to derive the underlying density profile. The method is adapted to lower statistics (i.e., lower number of bins per radial profile). In the present case, the PSF response matrix was computed for a circular Gaussian beam with FWHM = 10 arcmin, i.e., the angular resolution of our SZ maps. SZ cylindrical profiles deconvolved from the PSF were then deprojected into 3D profiles (assuming spherical symmetry for the clusters), and converted to pressure following Eq. (13).

We used a Monte Carlo (MC) method to propagate the measured error over the SZ profiles onto the PSF-corrected and deprojected profiles, accounting for the covariance between points expressed in the form of the matrix C_i for the cluster i. For each cluster i, its covariance matrix was Cholesky decomposed (i.e., ). We assumed correlated Gaussian noise, and generated m = 10 000 realisations of the SZ profile, P_i, with respect to this decomposition (i.e., , where x is an array of random numbers following a normal distribution). With this method we make sure we properly sample the noise properties of each profile. Each realisation was then rebinned logarithmically in radius and run through our deconvolution and deprojection algorithm. We conservatively chose to feed to the algorithm values of the weight per point equal to the inverse of the square of the flux dispersion in this bin (i.e., ). The weights remained the same for all realisations. These m realisations of the pressure profile were used to compute the covariance matrix of the cluster pressure profile similarly to what is described in Sect. 4.3.1. With this MC process, we ensure a proper propagation of the correlated noise through the PSF-correction and deprojection of the SZ profile.

6.2. The stacked pressure profile

As we did for the SZ profiles, the 62 pressure profiles were rescaled in the x- and y-axis directions according to R₅₀₀ and P₅₀₀, and then stacked together.

To further compare with A10’s results, we accounted for the slight difference in mass range between the REXCESS sample and ours via the factor (Eq. (11), and see also Sun et al. 2011). We divided our stacked profile by the average value across our ESZ-XMM sample, i.e.,  ⟨ f(M) ⟩  = 1.09.

Our three reconstruction methods lead to compatible stacked pressure profiles. However, we have accounted for their fluctuations by adding to the diagonal terms of the covariance matrix of the stacked profile the maximum point-to-point difference between the MILCA profile and the two other methods. In parallel to our Planck SZ analysis, we have also derived the pressure profiles of our 62 clusters from the analysis of XMM-Newton data. We followed the method applied to the REXCESS sample and presented in A10. Details of the XMM-Newton analysis and results are provided in Démocles et al. (in prep.).

The Planck and XMM-Newton stacked pressure profiles derived from our sample of 62 ESZ-XMM clusters, are displayed in Fig. 4. The two sets of data agree remarkably well. They are fully compatible within their respective dispersions, and they overlap in the radial range  ~(0.3−1) R₅₀₀. We recall that the points from the Planck data suffer a degree of correlation at about the 20% level, whereas all the points in the XMM-Newton profile can be considered as uncorrelated. These two independently derived pressure profiles are fully complementary. Indeed, in the central region the pressure profiles are very well constrained from the X-ray data, whilst the Planck measurements hardly reach down to 0.1 R₅₀₀ due to the moderate spatial resolution of our SZ maps. At larger radii, the extent of the X-ray observations is limited to radii smaller than  ~(0.7−1.0) R₅₀₀ because of the XMM-Newton field of view and of the quick drop in X-ray surface brightness with increasing radius. The Planck profile extends far beyond this radius.

With this joint constraint, we bring for the first time a comprehensive observational view of the distribution of the average thermal pressure distribution in clusters of galaxies out to a density contrast of δ ~ 50–100.

Following A10, we have investigated two sub-samples of our ESZ-XMM clusters. Keeping the cool-core (CC) versus non cool-core (non-CC) classification as provided in Planck Collaboration (2011f), we computed the stacked profiles for the 22 CC versus 40 non-CC clusters. Both stacked profiles are displayed in the right panel of Fig. 4. They are, as expected, different in the central parts, with a more peaked profile for CC systems and a shallower one for non-CC clusters. However, in the outer parts, our observed profiles for the two subsamples have very similar slopes.

Fig. 5 Marginalised posterior likelihood distribution for our best MCMC fit to a GNFW pressure profile with free parameters [P0,c500,α,β] and γ = 0.31. The best fit values of the parameters are marked as white crosses and solid red lines, respectively, on the 2D and 1D posterior distributions. On the 2D distributions, the white triangles mark the value of the A10 parameters, while the red, purple and green solid lines define, respectively, the 68, 95 and 99% confidence levels.

6.3. Best-fit to the generalised NFW profile

We have combined the Planck and XMM-Newton pressure profiles into a joint fit to the generalised NFW model (Eq. (7)). The pressure profiles derived over our sample from the Planck and XMM-Newton data cover a radial range from 0.02 × R₅₀₀ to 3 × R₅₀₀, with an overlap in the radial interval of (0.1−1) × R₅₀₀. In order to find the best parametric representation of our observed stacked profile, we ran a Monte Carlo Markov chain analysis to find the maximum likelihood solution parametrising the GNFW profile. The Planck and XMM-Newton data are two different and fully independent observational data sets. We therefore computed the likelihood as the product of the two independent likelihoods: ℒ = ℒ_X × ℒ_SZ. For the Planck data, we have accounted for the correlation between points through the use of the covariance matrix (see Sect. 6.1). We also accounted for the dispersion about the average X-ray and SZ profiles across the sample by quadratically adding it to the errors bars of the XMM-Newton points and to the diagonal elements of the covariance matrix for the Planck profile.

We performed our MCMC analysis in log–log space for various combinations of free parameters. Fixed parameters were assigned the values from the best fit provided by A10. Our best fit parameters are reported in Table 1. In order to find the best possible analytical representation of our data, we assessed the quality of the fit via the value of the reduced χ² using only the statistical errors, . The configuration with four and five free parameters are of equivalent quality, the former being slightly better (i.e., and 1.0, respectively). Unsurprisingly the degeneracy between free parameters increases with their number, translating into a drastic increase in the parameter covariance. We therefore adopted the four free parameter approach as our best parametrisation of a GNFW profile of our data, i.e., [P₀,c₅₀₀,α,β] free, with γ fixed to 0.31. The associated values are marked in bold face in Table 1, and are given purposefully with high precision with respect to the errors on parameters. These uncertainties are illustrated by the corresponding marginalised posterior likelihood distributions that are shown for individual (1D) and pairs of parameters (2D) in Fig. 5.

In parallel, following exactly the same procedure, we also fitted the stacked profiles for the sub-samples of cool-core and non cool-core clusters. The best fit parameters are also reported in Table 1. Both the stacked profiles and best fit model are shown on the right panel of Fig. 4.

We further adopted this four free parameter configuration to fit each of our measured individual pressure profiles, combining XMM-Newton and Planck pressure data. As describe in the above section for the stacked profile, we accounted for the point-to-point differences between our three reconstruction methods (see Appendix A) in the error budget. Results are reported in Appendix C, Fig. C.1 and Table C.1.

7. Discussion

7.1. The observed pressure profile

7.1.1. The core of clusters

Our present results are in very good agreement with the universal profile derived by A10 outside the core and down to R₅₀₀. Within the core, i.e., R < 0.15 R₅₀₀, our observed profile lies significantly below the A10 profile. In comparison to the REXCESS sample which is X-ray selected, our ESZ-XMM sample is SZ selected (except for the intersection with the XMM-Newton archives) and is thus, a priori, closer to being a mass-selected sample. It therefore contains more dynamically disturbed clusters (e.g., A2163, RXJ2228+2036, etc.). As has also been indicated by results in the validation follow-up of Planck clusters with XMM-Newton (Planck Collaboration 2011d, 2012, 2013b), this seems to indicate that X-ray selection under-samples the population of morphologically disturbed clusters. When splitting our sample into CC and non-CC clusters, the average profile for the CC clusters (and its associated best fit) is in better agreement with the universal profile given in A10. However, it is still slightly lower in the central parts than the A10 profile for CC clusters (see right panel of Fig. 4). There might thus still be differences between the population sampled by our CC clusters and those from A10. On the contrary, our non-CC cluster profile agrees well with the A10 non-CC profile, except for the very central parts (R < 0.04 R₅₀₀).

7.1.2. The inner profile

Observationally, it is hard to accurately determine the distribution of the SZ signal and, thus, of the underlying thermal pressure. Early single target SZ studies were limited to the clusters’ inner regions (i.e., within R₅₀₀), and were not really competitive with the X-ray measurements in terms of resolution and sensitivity (see e.g., Pointecouteau et al. 2001; Halverson et al. 2009; Korngut et al. 2011). However, the combination of SZ and X-ray tracers has already shown its potential in terms of structural studies of clusters (Pointecouteau et al. 2002; Kitayama et al. 2004; Jia et al. 2008; Basu et al. 2010). The first studies working on SZ cluster samples were presented by the SPT and ACT collaborations (Plagge et al. 2010; Sehgal et al. 2011), with 15 and 9 high significance clusters, respectively. The SPT team found good agreement between the X-ray predicted signal and their SZ measurement within R₅₀₀. Similarly, in a recent study, Bonamente et al. (2012) used a sample of 25 clusters observed with SZA/CARMA (Muchovej et al. 2007) and found an excellent agreement when modelling the SZ emission over their sample, either with the universal pressure profile from A10 or with the model from Bulbul et al. (2010). Our SZ and X-ray data are in excellent agreement over the radial range (0.1–1) R₅₀₀, providing tighter observational constraints.

Our results, together with the SPT, and ACT results, have to be compared with other works that reported a difference between SZ and X-ray measurements. Two papers based on WMAP data reported a measured SZ signal lower than the expected signal from their X-ray properties. Based on a WMAP (one year) analysis for 31 randomly selected nearby clusters, Lieu et al. (2006) found a ratio of one fourth between the measured SZ signal and the expected signal from X-ray constraints. An analysis based on WMAP-7 data for 49 clusters (Komatsu et al. 2011) found this ratio to be 0.5–0.7. However, these claims did not agree with other independent WMAP-5 statistical analysis for 893 clusters (Melin et al. 2011). Neither do we, obtaining instead an excellent agreement between X-ray and SZ properties within R₅₀₀, in particular, when using the Planck frequencies overlapping those of WMAP (see profiles for 30, 44 and 70 GHz in Fig. 1). The WMAP fluxes derived by Lieu et al. (2006) and Komatsu et al. (2011) are also discrepant with fluxes from OVRO/BIMA (Bonamente et al. 2006). These differences could be due, for instance, to complex effects involving large-scale beams and/or non-linear gains (e.g. Whitbourn et al. 2011, and references therein).

7.1.3. The outskirts

Beyond R₅₀₀, the observational situation is even more wanting. Only a few X-ray observations with XMM-Newton, Chandra and/or Suzaku constrain the density, temperature or gas fraction profiles out to  ~R₂₀₀ ≃ 1.4 R₅₀₀ (George et al. 2009; Reiprich et al. 2009; Urban et al. 2011; Simionescu et al. 2011; Walker et al. 2012). This type of X-ray detection remains very challenging and requires very long exposure times, as at larger radii the X-ray emission is extremely faint. We recall that beyond R₅₀₀ the universal pressure profile from A10 was constrained by predictions from numerical simulations, not by observations. Probing the gas with the SZ effect is therefore a powerful alternative, as shown, for instance, by the SPT average emission over 15 clusters (Plagge et al. 2010) or the WMAP statistical analysis over  ~700 clusters by Atrio-Barandela et al. (2008), though the latter is affected by the limited resolution and sensitivity of the survey. Conversely, from the tentative analysis of WMAP-3 data for 193 clusters with X-ray temperature above 3 keV, Afshordi et al. (2007) provided constraints for the cluster pressure profile (significant out to  ~1.3 R₅₀₀). In agreement with these earlier works, the Planck measurements of the present study are the first to allow a precise description of the thermal pressure distribution out to the cluster outskirts.

As shown in Appendix C, Planck also resolves some individual pressure profiles. The work presented by the Planck Collaboration (2013c) on the Coma cluster is even more striking, as the SZ profile detection reaches beyond 3 × R₅₀₀, i.e., as far as the statistical measurement presented in this work. The derived best fit model for the Planck observations of Coma is in full agreement with our stacked result.

While in the inner parts the profiles of CC and non-CC clusters are significantly different, and although the CC profiles lie just below the non-CC beyond  ~1 R₅₀₀, the two are compatible within our statistical errors. The same is true when comparing our best-fit models and A10’s for CC and non-CC clusters. This suggests that across our sample the average differences between the outskirts of the two types of clusters are smaller than the scatter between clusters.

We have also investigated the change in integrated Comptonisation parameter derived by comparing our best fit of a GNFW model to that of A10. Within a fixed aperture R₅₀₀, known from ancillary data, we have computed Y₅₀₀ for each cluster in our sample and for both sets of parameters (see Eq. (B.2)). The average ratio across the sample between the value of Y₅₀₀ from our best profile and that of A10’s is 1.02 ± 0.03. The difference is thus marginal as the two profiles are alike, and consistent with the ratio of 1.01 between the two GNFW parametrisations at R₅₀₀. This consistency demonstrates the robustness of the previously-published Planck SZ analysis and scaling relations where the A10 profile was adopted as a fiducial model for known R₅₀₀ values (Planck Collaboration 2011d,e,f, 2012, 2013a,b). The ratio for the values of Y for the two profiles is 1.13 and 1.19 within fixed apertures of 3 R₅₀₀ and 5 R₅₀₀ respectively, as expected given the slightly flatter outer slope derived in this work with respect to that of A10.

7.2. Comparison with theoretical predictions

7.2.1. Sets of numerical data

We have investigated three sets of simulated clusters in order to compare to our combined Planck and XMM-Newton pressure profile. All are taken from numerical simulations of structure formation in a ΛCDM cosmology:

• 1.

  We first compared to the set of combined simulations used byA10 to derive their universal pressure profile together with theREXCESS data. They comprise 93, 14, and 88simulated clusters withM > 10¹⁴ M_⊙, from Borgani et al. (2004), Nagai et al. (2007a) and Piffaretti & Valdarnini (2008), respectively. We refer to this set of simulations as B04+N07+P08 hereafter.

• 2.

  The second set is built from the pressure profiles of 64 massive clusters with M > 3 × 10¹⁴ M_⊙ from a simulation by Dolag (in prep.).

• 3.

  The third comprises the 40 most massive clusters (i.e., M > 3 × 10¹⁴ M_⊙) provided from the numerical simulations described in Battaglia et al. (2010, 2012).

All simulations include treatment of radiative cooling, star formation and energy feedback from supernova explosions. Simulations by Dolag (in prep.) and Battaglia et al. (2010, 2012) have prescriptions for AGN feedback. The different simulation sets use different techniques and different implementations of the physical processes. This ensures a fair description of up-to-date theoretical predictions, hopefully bracketing the plausible range on the thermal pressure profiles distribution in clusters of galaxy.

We kept the B04+N07+P08 sample as used by A10, even though its mass limit is slightly smaller than our ESZ-XMM sample. The higher mass threshold for the two other samples enforce the best possible match to our data. As for the observed profiles (see Sect. 6.2), the simulated cluster profiles were renormalised by the factor  ⟨ f(M) ⟩ . The average value for this factor is 1.0, 1.08 and 1.03 for the B04+N07+P08, Dolag (in prep.) and Battaglia et al. (2010) simulations, respectively. Finally, we accounted for the differences in the definition of R₅₀₀ and M₅₀₀ between observations and simulations. We corrected from the differences between the true mass and the hydrostatic masses which impact the estimation of both R₅₀₀ and P₅₀₀ (see Eq. (10)). We refer to Sect. 4.1 of A10 for a complete discussion. Assuming a hydrostatic bias r_cor = M_true/M_HE = 1.15 (Kay et al. 2004; Piffaretti & Valdarnini 2008), we divided the values of R₅₀₀ and P₅₀₀ derived from the simulations by and , respectively.

All three sets of data are reported in Fig. 6. It is well beyond the scope of this paper to extensively discuss the comparison between theoretical predictions. Nevertheless, we note that they agree within their respective dispersions across the whole radial range. The Dolag (in prep.) and Battaglia et al. (2010, 2012) profiles best agree within the central part, and are flatter than the B04+N07+P08 profile. This is likely due to the implementation of AGN feedback, which triggers energy injection at the clusters centre, balancing radiative cooling and thus stopping the gas cooling. In the outer parts where cooling is negligible, the B04+N07+P08 and Dolag (in prep.) profiles are in perfect agreement. The Battaglia et al. (2010) profile is slightly higher, but still compatible within its dispersion with the two other sets. Here again it might be due to the specific implementation of the simulatons.

Fig. 6 Measured pressure profiles from Planck and XMM-Newton, displayed as in Fig. 4. The three shaded areas mark the dispersions about the average profiles for three samples of simulated clusters: in blue, the simulations from Borgani et al. (2004), Nagai et al. (2007a) and Piffaretti & Valdarnini (2008), which were used in A10 to derive the universal pressure profile together with the REXCESS data; in green, the simulated sample of clusters from Dolag (in prep.); and in brown the simulated clusters from Battaglia et al. (2012). The corresponding average profiles for each set of simulations are plotted as solid lines using the same colour scheme. In the bottom panel, we present the χ profiles between the observed profile and the simulated average profiles, taking into account their associated dispersion. See Sect. 6.2 for more details.

Fig. 7 Gas mass fraction profile derived from the combined Planck and XMM-Newton pressure profile, assuming for the temperature profile: (H1) the average best fit model across the sample from X-ray spectroscopy (red line and striped area); or (H2) the same but extrapolating beyond R500 to a constant value equal to the average temperature measured in the last radial bin across the sample (blue line and striped area). The green dashed curve marks the expected gas fraction profile assuming the A10 pressure profile and kT(r) as in hypothesis (1). The star gives fgas(r = R500) for REXCESS clusters with M500 > 5 × 1014 M⊙ (Pratt et al. 2009). In maroon we reproduce the gas mass fraction profile derived from Suzaku measurements for the Perseus cluster (Simionescu et al. 2011). The solid and dashed black lines mark the cosmic baryon fraction expected from CMB measurements (Komatsu et al. 2011) and the expected gas fraction, assuming that 12% of baryons are in stars, respectively. The shaded blue and red areas translate the dispersion in the pressure profile across the ESZ-XMM sample as shown in Fig. 4, assuming hypothesis (H1) for kT(r). See Sect. 7.3.

7.2.2. Comparison with the Planck pressure profile

Figure 6 compares our SZ+X-ray observational constraints to the theoretical predictions described above. Overall, our observed pressure profile lies within the scatter across the profiles from these various sets of simulated clusters.

In the central parts, the data points lie on the lower envelope of the scatter of the B04+N07+P08 simulations, similarly to the REXCESS sample profiles (left panel of Figs. 4 and 7 in A10). This flatter inner slope is more compatible with that of the clusters from simulations which include the effect of AGN feedback (Dolag, in prep.; Battaglia et al. 2012), although our measured points lie above them.

With increasing radius (i.e., R > R₅₀₀) both SZ and X-ray profiles are marginally compatible with the stacked profiles for simulated clusters. Interestingly our profile lies above the B04+N07+P08 and Dolag (in prep.) simulated clusters, and below the Battaglia et al. (2012) clusters. In the outer parts of halos, the predicted pressure in numerical simulations is essentially sustained by gravitational heating. The general agreement (within the dispersion) with our measured pressure profile favours the equilibrium between ions and electrons in the outskirts, i.e., both populations have an equivalent temperature. Indeed, a difference would result in a drop in thermal pressure (Rudd & Nagai 2009). This may suggest that we have the correct global understanding of the intra-cluster thermal pressure support. However, the discrepancies between simulations and the average observed profile in the clusters’ outer parts calls for a deeper investigation of the gas distribution and structure in clusters. Thus we need more detailed modelling of the baryon physics in cluster outskirts, where incomplete virialisation of the intra-cluster medium affects the thermal pressure. Even if a universal behaviour of cluster quantities is anticipated from the simplest model of gravitational collapse, dispersion is also expected, even without considering the effects of non-gravitational physics, as shown by numerical simulations of halo formation (e.g., Navarro et al. 1997; Reed et al. 2011; Bhattacharya et al. 2011; Gao et al. 2012). This is plainly reflected in the large level of dispersion that we observe across our sample of simulated clusters.

On the observational side, our investigation of the sub-sample of CC and non-CC clusters did not show differences in the outer shape of the pressure profile (CC and non-CC cluster profiles are compatible within the measured statistical error bars). This might indicate that the aforementioned dispersion is not mainly driven by the dynamical state of clusters, which likely affects primarily the inner parts of halos. If our Planck constraints over the whole ESZ-XMM sample give a flatter outer shape with respect to the average simulated pressure profiles by B04+N07+P08 and Dolag (in prep.), we need nonetheless to keep in mind that our current sample may be affected by selection biases, as discussed in Planck Collaboration (2011f).

7.3. Constraint on the gas mass fraction

As already mentioned, X-ray measurements hardly reach density contrasts of  ~ 200 (e.g., George et al. 2009; Reiprich et al. 2009; Urban et al. 2011; Simionescu et al. 2011), so modelling beyond R₂₀₀ only relies on predictions from numerical simulations. We used our Planck and XMM-Newton derived pressure profile to investigate the average gas mass fraction distribution across our sample. Assuming the ICM to be a perfect gas, the thermal pressure is the genuine product of the temperature and the density, i.e., P(r) ∝ n_e(r) × kT(r). Following previous X-ray works (Pratt et al. 2007; Arnaud et al. 2010; Démoclès et al. 2010), we used the best fit to the spectroscopic temperature profile (Démocles et al., in prep.) with the analytical function proposed by Vikhlinin et al. (2006) for each of our clusters. Beyond the reach of the X-ray observations (i.e., [0.5–1.2] R₅₀₀), we have worked with two hypotheses for the temperature: (H1) the extrapolation of the best fit model (note tha ); or (H2) a constant temperature fixed to the last radial bin value. We have then derived the density profiles, and integrated them over the cluster volume to obtain the gas mass radial distribution. Meanwhile, we modelled each total mass profile with an NFW profile with the mass, M₅₀₀, derived from Eq. (9), and the concentration estimated from the c- − M relation by Bhattacharya et al. (2011). We checked that the average of the individual NFW mass models agrees with the NFW model derived from the average values of M₅₀₀ and c₅₀₀ across the sample. The average gas fraction profile is computed as the ratio of the average gas mass model and the total mass model profiles.

Figure 7 shows the resulting gas mass fraction distribution for the two hypotheses on the temperature at large radii (red and blue solid lines). The red and blue striped areas picture the statistical errors measured from the SZ+X-ray pressure profile. The red and blue shaded areas overplot the XMM-Newton and Planck dispersion onto the pressure profile (see Fig. 4) using hypothesis (H1) for the temperature, and provide an illustration of the cluster-to-cluster variation of f_gas within our sample.

From the stacked extrapolated temperature models, the average temperature is 1.3 keV at 3 R₅₀₀, which leads to f_gas ~ 0.2. An underestimation of this temperature by a third (i.e., kT ~ 1 keV) would increase the gas fraction likewise (i.e., f_gas  ~  0.3). For the lower bound, hypothesis (H2) fixes the temperature of the averaged profile beyond R₅₀₀ to  ~  4 keV, which likely overestimates the true average value of the temperature across our sample in the (2–3) R₅₀₀ radial range. Therefore we are confident that both hypotheses fairly bracket the range for the gas fraction distribution in massive clusters out to 3 R₅₀₀. In other words, our constraints on the pressure profile imply that the temperature required to flatten the gas fraction profile at 3 R₅₀₀ to the expected cosmic value lies between 1 and 4 keV and is likely closer to the lower value.With hypothesis (H1) for kT(r), the green dashed line gives the gas fraction for the A10 pressure profile. As for our SZ+X-ray pressure profile, it leads to f_gas values above the expected gas fraction out to 3 R₅₀₀.

At R₅₀₀, the measurement of f_gas for massive clusters presented by Pratt et al. (2009) (see also Giodini et al. 2009, for massive clusters with (4–10) × 10¹⁴ M_⊙) is lower by a factor 1.3 with respect to our estimate (although compatible within  ± 2σ). We can also compare our result to the gas fraction profile reported by Simionescu et al. (2011) using Suzaku measurements in the direction of the Perseus cluster. The derived gas fraction profile is in very good agreement with ours within R₅₀₀. Beyond this radius, it rises within our plausible dispersion envelope (red and blue shaded areas in Fig. 7). Note however, that such a time consuming X-ray observation sampled less than 5% of the sky area covered by the cluster and as a specific measurement is neither representative nor out of the allowed range for the gas fraction in cluster outskirts.

The constrained interval for the gas fraction from our observed pressure profile encompasses the cosmic gas fraction derived from CMB measurements (solid black line in Fig. 7, from WMAP-7 results by Komatsu et al. 2011) and the expected hot gas fraction (dashed black line) assuming that 12% of cluster baryons are in stars (Gonzalez et al. 2007; Giodini et al. 2009).

8. Conclusion

From the Planck nominal mission (i.e., 14 months of survey), making use of its full spectral coverage of the SZ spectrum, we have extracted and reconstructed the SZ signal distribution in clusters of galaxies for a sample of 62 massive nearby clusters. All are individually detected in the survey with high significance (Planck Collaboration 2011c) and were previously used to investigate the total integrated SZ flux and the SZ scaling relations (Planck Collaboration 2011f). We have scaled and averaged together all the SZ profiles in the sample into a stacked profile. We have statistically detected the SZ signal out to 3 × R₅₀₀, providing, for the first time, stringent observational constraints on the ICM gas beyond R₅₀₀ and out to a density contrast of δ  ~  50–100.

From the Planck SZ raw profiles, we have derived the underlying thermal pressure profiles of the ICM. Together with the pressure profile derived from the XMM-Newton data, we have provided for the first time a comprehensive observational view of the distribution of thermal pressure support in clusters from 0.01 out 3 × R₅₀₀. We have fitted these unique measurements to a generalised NFW profile. Our best analytical representation over this wide radial range is given by the parameters  [P₀,c₅₀₀,γ,α,β ]  =  [6.41,1.81,0.31,1.33,4.13 ] .

Our observational measurements further confirm the agreement of the SZ and X-ray constraints on the intra-cluster gas properties within the inner part of the clusters commonly considered to be virialised, i.e., R₅₀₀, as found in Planck Collaboration (2011e) and Planck Collaboration (2011f). Overall it also agrees with a wide range of simulated clusters through the whole radial range, although in the central regions it matches best the numerical predictions that implement prescriptions for AGN feedback. The statistical nature of our stacked detection at large radii provides the average trend for the thermal pressure distribution in cluster outskirts out to 3 × R₅₀₀, which is slightly flatter than most theoretical predictions.

In conjunction with X-ray constraints on the temperature profile of our clusters, we have derived the profile of the gas mass fraction out to the cluster outskirts. From reasonable hypotheses on the gas temperature to account for the lack of constraints beyond R₅₀₀, we have bracketed a range for the gas fraction in the cluster outer regions, which is compatible with the cosmic baryon fraction and the expected gas fraction in halos.

The processes governing the thermodynamical state of the outer regions in clusters still need to be understood from the theoretical and observational points of view (see Kravtsov & Borgani 2012, for a review and references). Issues such as gas clumping will affect the pressure estimation (i.e., effect on the X-ray surface brightness, see Roncarelli et al. 2006), departures from hydrostatic equilibrium and contribution from non-thermal pressure (due to magnetic fields and/or cosmic rays) will modify the gas fraction. In this regard SZ observations provide a straightforward description of the thermal pressure distribution in massive halos, and a clear path to the gas fraction determination. Future SZ instruments with increased spatial resolution and sensitivity, and retaining the ability to map clusters of galaxies out to large radii will certainly provide us with further details and insight. In the meantime, with Planck we bring unique observational constraints that are extremely valuable to further test and understand the physics at play in the outskirts of clusters.

Online material

Appendix A: SZ map reconstruction methods

In this Appendix we discuss the different SZ map reconstruction methods applied to the Planck data in Sect. 4.1.

• 1.

  MILCA (Hurier et al. 2010):The thermal SZ signal reconstruction is performed on thesix Planck all-sky maps from 100 GHzto 857 GHz. MILCA (Modified InternalLinear Combination Algorithm) is a component separa-tion approach aiming at extracting a chosen component(here the thermal SZ signal) from a multi-channel set ofinput maps. It is mainly based on the well known InternalLinear Combination approach (e.g., Eriksenet al. 2004), that searches for thelinear combination of the input maps that minimises the varianceof the final reconstructed map by imposing spectral constraints.In this paper, we applied MILCA using two constraints: preserva-tion of the thermal SZ signal (knowing the SZ spectral signature);and removal of the CMB contamination in the final SZ map (alsomaking use of the well known spectrum of the CMB). In addition,to compute the weights of the linear combination, we have usedthe extra degrees of freedom to minimise residuals from othercomponents (2 degrees) and from the noise (2 degrees). The noisecovariance matrix was estimated from the frequency error maps(see Sect. 2).

• 2.

  NILC: Needlet ILC performs a linear combination of the observed maps which has minimum variance under the constraint of offering unit response to the component of interest (here the thermal SZ, whose frequency scaling is known). The weights of the ILC depend on the covariance between the various observations, and can be computed, for example, on domains of the observed pixels (pixel space) or angular scales (spherical harmonics). In the case of NILC, covariances (and hence weights for component separation) are computed independently in domains of a needlet decomposition (spherical wavelet frame). The needlet decomposition provides localisation of the ILC filters both in pixel and in multipole space, allowing us to deal with local contamination conditions varying both in position and in scale. NILC was developed to extract a CMB map from WMAP data (Delabrouille et al. 2009) and was also tested for SZ effect extraction in Leach et al. (2008). Multi-component extensions have been investigated by Remazeilles et al. (2011).

• 3.

  GMCA (Bobin et al. 2008) is a blind source separation method developed for separating sources from instantaneous linear mixtures. The components are assumed to be sparsely represented (i.e., have a few significant samples in a specific basis) in a so-called sparse representation Ψ (typically wavelets). The assumption that the components have a sparse representation in the wavelet domain is equivalent to assuming that most components have a certain spatial regularity. These components and their spectral signatures are then recovered by minimising the number of significant coefficients in Ψ. Recently, L-GMCA has been further introduced to analyse the CMB data in a local and multi-scale manner (Bobin et al. 2008). More precisely, the multi-channel data are analysed in four frequency bands in spherical harmonics (i.e., wavelet bands). In each wavelet band, GMCA is applied locally on small patches with band-dependent sizes. In each band, the observations are analysed at the same band-dependent resolution. In Bobin et al. (2008), it was shown this local analysis approach enhances the separation quality. The spectral signatures of CMB, free-free and SZ are assumed to be known.

Fig. A.1 Comparison between the SZ profile reconstruction from MILCA, NILC and GMCA The comparison is presented in the form of the averaged stacked radial profile from our 62 clusters of the ESZ-XMM sample, individually scaled respectively in the x- and y-axis directions according to R500 and Φ (Eq. (12)). The light-red shaded area marks the dispersion about the average stacked profile for the MILCA method. The points within each individual profile are correlated at about the  ~ 20% level (see Sect. 4.2) and the plotted errors correspond to the square root of the diagonal elements of the covariance matrix of the profile. The solid black line (labelled in the legend as “Expected from PA10”) is the average stacked profile obtained from the expected SZ individual profiles drawn from the universal pressure profile by A10, parameterised according R500 and Y500 derived from the XMM-Newton data analysis presented in Planck Collaboration (2011f). The two dotted lines indicate the associated dispersion across the sample. The bottom panel shows the value of χ at each point of the observed profiles with respect to the expectations from the universal profile of A10.

Appendix B: Validation on simulations

B.1. Profiles

We have validated our reconstruction of the SZ signal and profile computation methods with simulations. We used the Planck sky model (PSM, Delabrouille 2012) to simulate the whole sky as seen by Planck. The PSM includes most astrophysical components acting as foreground or background contamination with respect to the SZ signal (i.e., Galactic dust, free-free and synchrotron emission, extragalactic sources in the infrared and radio parts of the spectrum and CMB). For the instrumental noise we added to the simulation the all-sky error map drawn from the Planck jack-knife maps at each frequency (see Sect. 2), thus mimicking the noise properties of the actual survey (noise and systematic effects). For the SZ component, we built a full sky template for all the 1743 clusters from the MCXC meta-catalogue (Piffaretti et al. 2011), which includes all 62 clusters in our sample. Each cluster was modelled assuming spherical symmetry and a thermal pressure distribution following the universal profile from A10. We recall that the universal pressure profile provided by these authors is the best fit of a generalised NFW profile (Nagai et al. 2007a) to the median profile derived from the REXCESS pressure profile and the predictions from three different sets of numerical simulations (i.e., from Borgani et al. 2004; Nagai et al. 2007a,b; Piffaretti & Valdarnini 2008). In our case, the universal pressure profile is parametrised using the values of R₅₀₀ and M₅₀₀ from the MCXC, together with the M–Y_X scaling relation from Arnaud et al. (2007, 2010) to estimate Y₅₀₀.

We processed the simulated Planck sky maps through each of our SZ map reconstruction methods. For our sample of 62 clusters, we extracted the patches and computed the individual profiles as described in Sect. 4.2. We stacked the simulated profiles for the three methods after renormalisation, making use of the aforementioned R₅₀₀ and Y₅₀₀ values from the MCXC prescription.

The resulting stacked profiles are displayed in Fig. A.1. All three tested methods agree remarkably well over the whole radial range. Taking into account the correlated errors of each profile, the χ² over the radial range (0–10) × R₅₀₀ for NILC and GMCA with respect to MILCA are 0.78 and 1.54, respectively. The reconstructed profiles are also in good agreement with the reconstructed stacked model. Accounting for the errors on the model and the correlated errors between points on the reconstructed profiles, the reduced χ² for MILCA, NILC and GMCA with respect to the stacked model profile are 0.80, 0.66 and 1.03 within 3 × R₅₀₀. This excellent agreement also translates into a small relative error below 15% in the same range of radii. Beyond R ~ 2 R₅₀₀, the three reconstructed profiles drop slightly below the input model, even though they still agree within their dispersions and errors. This drop is caused by an intrinsic over-prescription of point sources in the PSM simulations with respect to real data. Despite our careful masking, the unidentified or unresolved extra sources in the PSM affect the power spectrum of the noise, increasing the background level and making it more difficult to detect the SZ signal further out in the simulations than in real data. We stress that this is an effect limited to our test simulations, and not relevant for the real data.

B.2. Consistency between Y₅₀₀ measurements

Fig. B.1 Comparison of Y500 values from the Planck cluster detection method (MMF3) against our best fit for each of our 62 clusters (from NILC in this case). In both cases the distribution of thermal pressure is assumed to follow a universal pressure profile (A10) for which the scaled radius R500 is fixed to the XMM-Newton best fit value. See Sect. 6.3 for details.

To further validate the above, and as a consistency check with previous Planck results, we fitted each individual SZ profile with a projected, PSF-convolved universal pressure profile. We fixed R₅₀₀ to the best fitting X-ray value from Planck Collaboration (2011f) and only fitted the normalisation, Y₅₀₀. For a given profile, P, with an associated covariance matrix, C, and a model M, the χ² statistic can be expressed as (B.1)The solution that minimises χ², Y₅₀₀, and its associated uncertainty, σ_Y500, are analytic: We compared the resulting value of Y₅₀₀ to those obtained from the detection algorithms PwS and MMF3. Both detection algorithms were run under the same assumptions, i.e., universal pressure profile with the “non-standard” parameterisation and a fixed size (R₅₀₀). The fluxes were derived from the Planck nominal mission survey maps. These fluxes are in agreement with those extracted from the Planck first year survey (Planck Collaboration 2011f). For instance, for MMF3, the average ratio is 1.04 ± 0.09, compatible within the average relative errors with the fluxes extracted from the nominal mission.

The relation between the integrated Comptonisation parameters derived from the SZ profiles and from the detection methods is displayed for MMF3 vs NILC fuxes in Fig. B.1. A BCES linear fit accounting for the errors in both x- and y-axis directions (Akritas & Bershady 1996) yields a best fitting slope and normalisation of 1.00 ± 0.02 and 0.01 ± 0.01 (), respectively. There is 4% dispersion about this relation. Furthermore, the median and associated deviation for the one-to-one ratios of PwS and MMF3 fluxes compared to the fitted NILC profile values are 1.03 ± 0.11 and 0.98 ± 0.08, respectively. The median ratios for the errors are 1.10 ± 0.36 are 1.64 ± 0.34. The errors associated with our best fit values are in very good agreement with those from PwS; however, they are smaller than those derived from the MMF3 method. This difference with MMF3 might be explained by the different treatment of the noise. Whilst we estimated the noise from the surroundings of the clusters, defined according to the scaled radius of each target, and PwS also used a local computation of the noise, the MMF3 algorithm used a larger area to characterise the noise properties. Our tests against the PwS and MMF3 fluxes emphasise our grasp on the characterisation of the local noise for each target, and further assess the reliability of our estimates of the statistical errors in our stacked profiles.

Such agreement between the fluxes derived from our fit over the cluster profiles and those from the detection algorithm, strongly confirms the reliability of our SZ flux estimate and further strengthen the overall result stated in the Planck collaboration early papers that the SZ and the X-ray fluxes within R₅₀₀ are fully consistent (Planck Collaboration 2011c,d,e,f).

Appendix C: Pressure profiles of individual clusters

In this appendix we present the individual pressure profiles obtained from Planck and XMM-Newton data for each of the 62 clusters in our sample (see Fig. C.1). We show the agreement with the universal pressure profile (A10) and provide the best fitting model for each cluster to the GNFW profile in Table C.1. The following hard limit priors were adopted during the fitting procedure: 0 < P₀ < 100; 0 < c₅₀₀ < 10; 0 < α < 10; and 0 < β < 15.

Fig. C.1 Measured individual pressure profiles for the ESZ-XMM sample using the Planck and XMM-Newton data. As for the stacked pressure profile, the maximum difference point-to-point between the three SZ reconstruction methods have been taken into account in the Planck measurement error budget (Sect. 6.3). The best fit model on each individual profile is shown as a blue solid line. The black and green solid curves mark the best fit model to the stacked SZ+X-ray pressure profile (Sect. 6.2) and the A10 best fit model, respectively.

Fig. C.1 continued.

Fig. C.1 continued.

Fig. C.1 continued.

Acknowledgments

The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN and JA (Spain); Tekes, AoF and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and DEISA (EU). We acknowledge the use of the Healpix software.

References

All Tables

All Figures

Fig. 1 Average scaled SZ profiles for each of the nine Planck frequency bands (decreasing frequency from left to right and from top to bottom). The points within each individual profile are correlated at about the 20% level (see Sect. 4.3.1). At each frequency, the stacked radial profile is obtained from the average of the 62 individual cluster profiles scaled by R500 and Φi/ ⟨ Φi ⟩ , in the x- and y-axis directions, respectively (see Sect. 5.1). The light-red shaded area in each panel indicates the statistical uncertainty about the average observed profile. The solid black line corresponds to the prediction of the universal pressure profile (A10) for our cluster sample in each band, and the two dashed lines indicate its associated dispersion.

In the text

	Fig. 2 Stacked scaled image (size is 10 R500 on a side) using a logarithmic stretch in y. Individual maps are rescaled by Φi before averaging, and then multiplied by  ⟨ Φi ⟩ . The black and white circles mark the loci of 1 and 3 × R500, respectively.
In the text

Fig. 3 Left: individual profiles (grey lines) and scaled stacked profile (red points) for our sample of 62 clusters. The light-red shaded area marks the dispersion about the average stacked profile, with its upper and lower limits highlighted by the two dashed red lines. The points within each individual profile are correlated at about the 20% level (see Sect. 4.3.1). The error bars on the Planck data points are purely statistical and correspond to the square root of the diagonal elements of the covariance matrix. The solid black line (labeled in the legend as “Expected from PA10”) is the average stacked profile obtained from the expected SZ individual profiles drawn from the universal pressure profile of A10, parameterised according R500 and Y500 derived from the XMM-Newton data analysis presented in Planck Collaboration (2011f). The two dotted black lines indicate the associated dispersion about this model profile. The bottom panel shows the value of χ at each point with respect to the expectations from the universal profile taking into account the measured errors. Right: matrix of the correlation coefficient for the observed radial stacked profiles. Both panels are illustrated here for the MILCA results.

In the text

Fig. 4 Left: Planck pressure profile obtained from the average of the individual pressure profiles across our sample (red points) shown together with the stacked pressure profile derived from the XMM-Newton data for the same sample (purple points). The dispersions about the SZ and X-ray profiles are depicted, respectively, by the red and purple shaded areas. Our best fit GNFW profile is shown as a solid black line and that of A10 as a solid green line. The lower panel shows the χ profile of these two best models taking into account the statistical errors and dispersion about the observed profile. Right: stacked profile from Planck and XMM-Newton for the sub-samples of cool-core (open circles) and non cool-core (open stars) clusters within the ESZ-XMM sample. The shaded areas are identical to the one shown in the left panel. Our best fit models for each sub-sample are shown as black solid and dashed lines for the cool-core and non cool-core clusters, respectively (see Table 1). The best A10 fit for cool-core and non cool-core clusters are shown as green solid and dashed lines, respectively. The lower panel gives the χ profiles of our best fit models with respect to the measured profiles and associated errors (including dispersion). On both panels the Planck data points are correlated at about the 20% level (see Sect. 4.3.1). As for the stacked SZ profile, the error bars on the Planck points are purely statistical and correspond to the square root of the diagonal elements of the covariance matrix. For the XMM-Newton points they correspond to the statistical error on mean.

In the text

Fig. 5 Marginalised posterior likelihood distribution for our best MCMC fit to a GNFW pressure profile with free parameters [P0,c500,α,β] and γ = 0.31. The best fit values of the parameters are marked as white crosses and solid red lines, respectively, on the 2D and 1D posterior distributions. On the 2D distributions, the white triangles mark the value of the A10 parameters, while the red, purple and green solid lines define, respectively, the 68, 95 and 99% confidence levels.

In the text

Fig. 6 Measured pressure profiles from Planck and XMM-Newton, displayed as in Fig. 4. The three shaded areas mark the dispersions about the average profiles for three samples of simulated clusters: in blue, the simulations from Borgani et al. (2004), Nagai et al. (2007a) and Piffaretti & Valdarnini (2008), which were used in A10 to derive the universal pressure profile together with the REXCESS data; in green, the simulated sample of clusters from Dolag (in prep.); and in brown the simulated clusters from Battaglia et al. (2012). The corresponding average profiles for each set of simulations are plotted as solid lines using the same colour scheme. In the bottom panel, we present the χ profiles between the observed profile and the simulated average profiles, taking into account their associated dispersion. See Sect. 6.2 for more details.

In the text

Fig. 7 Gas mass fraction profile derived from the combined Planck and XMM-Newton pressure profile, assuming for the temperature profile: (H1) the average best fit model across the sample from X-ray spectroscopy (red line and striped area); or (H2) the same but extrapolating beyond R500 to a constant value equal to the average temperature measured in the last radial bin across the sample (blue line and striped area). The green dashed curve marks the expected gas fraction profile assuming the A10 pressure profile and kT(r) as in hypothesis (1). The star gives fgas(r = R500) for REXCESS clusters with M500 > 5 × 1014 M⊙ (Pratt et al. 2009). In maroon we reproduce the gas mass fraction profile derived from Suzaku measurements for the Perseus cluster (Simionescu et al. 2011). The solid and dashed black lines mark the cosmic baryon fraction expected from CMB measurements (Komatsu et al. 2011) and the expected gas fraction, assuming that 12% of baryons are in stars, respectively. The shaded blue and red areas translate the dispersion in the pressure profile across the ESZ-XMM sample as shown in Fig. 4, assuming hypothesis (H1) for kT(r). See Sect. 7.3.

In the text

Fig. A.1 Comparison between the SZ profile reconstruction from MILCA, NILC and GMCA The comparison is presented in the form of the averaged stacked radial profile from our 62 clusters of the ESZ-XMM sample, individually scaled respectively in the x- and y-axis directions according to R500 and Φ (Eq. (12)). The light-red shaded area marks the dispersion about the average stacked profile for the MILCA method. The points within each individual profile are correlated at about the  ~ 20% level (see Sect. 4.2) and the plotted errors correspond to the square root of the diagonal elements of the covariance matrix of the profile. The solid black line (labelled in the legend as “Expected from PA10”) is the average stacked profile obtained from the expected SZ individual profiles drawn from the universal pressure profile by A10, parameterised according R500 and Y500 derived from the XMM-Newton data analysis presented in Planck Collaboration (2011f). The two dotted lines indicate the associated dispersion across the sample. The bottom panel shows the value of χ at each point of the observed profiles with respect to the expectations from the universal profile of A10.

In the text

	Fig. B.1 Comparison of Y500 values from the Planck cluster detection method (MMF3) against our best fit for each of our 62 clusters (from NILC in this case). In both cases the distribution of thermal pressure is assumed to follow a universal pressure profile (A10) for which the scaled radius R500 is fixed to the XMM-Newton best fit value. See Sect. 6.3 for details.
In the text

Fig. C.1 Measured individual pressure profiles for the ESZ-XMM sample using the Planck and XMM-Newton data. As for the stacked pressure profile, the maximum difference point-to-point between the three SZ reconstruction methods have been taken into account in the Planck measurement error budget (Sect. 6.3). The best fit model on each individual profile is shown as a blue solid line. The black and green solid curves mark the best fit model to the stacked SZ+X-ray pressure profile (Sect. 6.2) and the A10 best fit model, respectively.

In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text