© ESO, 2016

1. Introduction

This paper, one of a set associated with the 2015 release of data from the Planck¹ mission, describes the Planck Compton parameter (y) map, which is part of the Planck 2015 data release.

The thermal Sunyaev-Zeldovich (tSZ) effect (Sunyaev & Zeldovich 1972) is produced by the inverse Compton scattering of cosmic microwave background (CMB) photons by hot electrons along the line of sight and in particular in clusters of galaxies. The tSZ effect has proved to be a major tool for studying the physics of clusters of galaxies, and also structure formation in the Universe. Catalogues of clusters of galaxies selected via the tSZ effect have become available in the last few years, including, for example, those from the Planck satellite (Planck Collaboration VIII 2011; Planck Collaboration XXIX 2014), the Atacama Cosmology Telescope (ACT, Hasselfield et al. 2013), and the South Pole Telescope (SPT, Reichardt et al. 2013; Bleem et al. 2015). These catalogues and their associated sky surveys have been used to study the physics of clusters of galaxies (Planck Collaboration XII 2011; Planck Collaboration XI 2011; Planck Collaboration X 2011) and their cosmological implications (Planck Collaboration XX 2014; Benson et al. 2013; Das et al. 2013; Wilson et al. 2012; Mak & Pierpaoli 2012).

The study of cluster number counts and their evolution with redshift using tSZ selected catalogues is an important cosmological test (Carlstrom et al. 2002; Dunkley et al. 2013; Benson et al. 2013; Planck Collaboration XX 2014; Hou et al. 2014). This study is complemented by the measurement of the tSZ effect power spectrum (Komatsu & Seljak 2002), for which no explicit measurement of cluster masses is required. Low-mass clusters, thus fainter in tSZ, which can not be detected individually, also contribute statistically to the measured signal (Battaglia et al. 2010; Shaw et al. 2010). Another complementary approach (as pointed out in Rubiño-Martín & Sunyaev 2003) is to study the higher-order statistics of the tSZ signal, and in particular the skewness or, equivalently, the bispectrum. The bispectrum of the tSZ effect signal is dominated by massive clusters at intermediate redshifts (Bhattacharya et al. 2012), for which high-precision X-ray observations exist. This contrasts with the power spectrum, for which most of the signal comes from the lower mass and higher redshift groups and clusters (e.g. Trac et al. 2011). Thus, theoretical uncertainties in the tSZ bispectrum are expected to be significantly smaller than those in the estimation of the tSZ power spectrum. Therefore, combined measurements of the power spectrum and the bispectrum can be used to distinguish the contribution to the power spectrum from different cluster masses and redshift ranges. However, contamination from point sources (Rubiño-Martín & Sunyaev 2003; Taburet et al. 2010) and other foregrounds (Planck Collaboration XX 2014) needs to be dealt with carefully.

As shown in the Planck 2013 results (Planck Collaboration XX 2014), the all-sky coverage and unprecedented wide frequency range of Planck allowed us to produce all-sky tSZ Compton parameter maps, also referred to as y-maps. From these maps we can obtain an accurate measurement of the tSZ power spectrum at intermediate and large angular scales, for which the tSZ fluctuations are almost insensitive to the cluster core physics. Furthermore, the expected non-Gaussianity properties of the tSZ signal can be studied using higher order statistical estimators, such as the skewness, and the bispectrum. The y-map can be also used to extract the tSZ signal on regions centred at cluster positions and in particular to perform stacking analysis. Finally, this map may also be cross correlated with other cosmological probes (see for example Ma et al. 2015; Van Waerbeke et al. 2014; Hill & Spergel 2014) as a consistency test.

In this paper we construct revised tSZ all-sky maps from the individual Planck frequency maps. With respect to the Planck Compton parameter map in the Planck 2013 results (Planck Collaboration XX 2014), we have extended the analysis to the full mission data set and performed an in-depth characterization of its statistical properties. This extended analysis allows us to deliver this map as part of the Planck 2015 data release. The paper is structured as follows. Section 2 describes the Planck data used to compute the tSZ all-sky maps and the simulations used to characterize them. Section 3 presents the reconstruction of the Planck all-sky Compton parameter map. In Sect. 4 we discuss the validation of the reconstructed y-map in the pixel domain including signal and noise characterization. Section 5 describes the power spectrum analysis. Cross-checks using higher order statistics are presented in Sect. 6. The cosmological interpretation of the results is discussed in Sect. 7, and we finally present our conclusions in Sect. 8.

2. Data and simulations

2.1. The Planck data

This paper is based on the Planck’s full mission (Planck Collaboration I 2014; Planck Collaboration I 2016), corresponding to five and eight full-sky surveys for the High Frequency Instrument (HFI) and Low Frequency Instrument (LFI) data (Planck Collaboration II 2016; Planck Collaboration VII 2016; Planck Collaboration VIII 2016), respectively. The Planck channel maps are provided in HEALPix (Górski et al. 2005) pixelization scheme at N_side = 2048. A noise map is associated with each channel map. This map is obtained from the half difference of maps made from the first and second halves of each stable pointing period (also called “ring”). In the half-difference maps the astrophysical emission cancels out, which makes them a good representation of the statistical instrumental noise. These half-difference maps are used to estimate the noise in the final Compton parameter map. In addition, survey maps, which are also available for each channel, will be used to estimate possible residual systematic effects in the y-map.

For the purpose of this paper we approximate the Planck effective beams by circular Gaussians, the FWHM estimates of which are given in Table 1 for each frequency channel.

2.2. Simulations

We also use simulated Planck frequency maps obtained from the full focal plane (FFP) simulations (Planck Collaboration XII 2016), which are described in the Planck Explanatory supplement (Planck Collaboration 2013).

These simulated maps include the most relevant sky components at microwave and millimetre frequencies, based on foregrounds from the Planck sky model (PSM, Delabrouille et al. 2013), namely CMB anisotropies, thermal SZ effect, diffuse Galactic emissions (synchrotron, free-free, thermal and spinning dust and CO), radio and infrared point sources, and the clustered cosmic infrared background (CIB). The tSZ signal was constructed using hydrodynamical simulations of clusters of galaxies up to redshift 0.3. For higher redshifts pressure profile-based simulations of individual clusters of galaxies randomly drawn on the sky were added. The noise in the maps was obtained from realizations of Gaussian random noise in the time domain and therefore accounts for noise inhomogeneities in the maps.

The simulation set also includes Monte Carlo noise-only realizations for each Planck channel map. These will also be used to estimate the noise properties in the final y-map.

3. Reconstruction of the all-sky tSZ maps

The thermal SZ Compton parameter (Sunyaev & Zeldovich 1972) in a given direction, , can be expressed as (1)where k_B is the Boltzmann contsant, m_e the electron mass, σ_T the Thomson cross-section, ds the distance along the line of sight, , and n_e and T_e are the electron number density and temperature.

In CMB temperature units the tSZ effect contribution to the Planck maps for a given observational frequency ν is given by (2)Neglecting relativistic corrections, g(ν) = x coth(x/ 2)−4, with x = hν/ (k_BT_CMB). Table 1 presents the conversion factors for Compton parameter to CMB temperature, K_CMB, for each frequency channel, after integrating over the instrumental bandpass (Planck Collaboration IX 2014).

3.1. Reconstruction methods

The tSZ effect signal in the Planck frequency maps is subdominant with respect to the CMB and other foreground emissions. By contrast to the CMB emission, the tSZ effect from galaxy clusters is spatially localized and leads to a highly non-Gaussian signal. CMB-oriented component-separation methods (Planck Collaboration XII 2014) are not adequate to recover the tSZ signal. We therefore use specifically tailored component separation algorithms to reconstruct the tSZ signal from the Planck frequency channel maps, as in Planck Collaboration XX (2014). These algorithms rely on the spatial localization of the different astrophysical components and on their spectral diversity to separate them. As in Planck Collaboration XX (2014) we consider here the MILCA (Modified Internal Linear Combination Algorithm, Hurier et al. 2013) and NILC (Needlet Independent Linear Combination, Remazeilles et al. 2011) methods. Both are based on the Internal Linear Combination (ILC) approach that searches for the linear combination of the input maps that minimizes the variance of the final reconstructed map, under the constraint of offering unit gain to the component of interest (here the tSZ effect, whose frequency dependence is known). Both algorithms have been extensively tested on simulated Planck data.

For both methods, the Planck HFI maps from 100 to 857 GHz are used convolved to a common resolution of 10′. In the case of NILC we also use the LFI data at large angular scales (ℓ< 300). Similarly, for both methods the 857 GHz map, which traces the thermal dust emission on large angular scales, is only used for multipoles ℓ< 300 to minimize residuals from IR point sources and clustered CIB emission. The final y-maps have a resolution of 10′.

3.1.1. NILC

In the multi-component extensions of NILC (Delabrouille et al. 2009; Remazeilles et al. 2011), initially developed to extract the CMB, the weights for component separation (i.e., covariances) are computed independently in domains of a needlet decomposition (in the spherical wavelet frame). The needlet decomposition provides localization 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. We imposed constraints to remove the CMB contamination and preserve the tSZ effect. To avoid strong foreground effects, part of the Galactic plane (corresponding to about 2% of the sky) was masked before applying NILC to the Planck frequency maps.

The localization in multipole space is achieved by using ten Gaussian window functions { h_j(ℓ) } _{1 ≤ j ≤ 10} as bandpass filters², so-called “needlet bands”, allowing for smooth localization in ℓ (see Fig. 1, solid lines). NILC performs a weighted linear combination of the bandpass filtered Planck maps for each needlet scale independently. The localization in the spatial domain is achieved by defining scale-dependent zones over the sky on which the covariance matrices and ILC weights are computed. More precisely, the pixel domain on which the covariance is computed is defined from the smoothing of the product of the relevant needlet maps with a symmetric Gaussian window in pixel space whose FWHM depends on the needlet scale considered. This avoids artificial discontinuities at pixel edges. The localization reduces the number of modes on which the statistics is computed; this may be responsible for an “ILC bias” due to chance correlations between SZ and foregrounds (Delabrouille et al. 2009). At the coarsest scale in particular (first needlet band), the area of the spatial localization must be large enough to counterbalance the lack of modes in multipole space due to spectral localization. In practice, the zones for spatial localization are not pre-defined but their area is automatically adjusted to the needlet scale considered. The ILC bias b_ILC is related to the number of channels N_ch and to the number of modes N_mDelabrouille et al. (2009) as (3)This offers the possibility to adapt N_m and N_ch in order to control the ILC bias to a set threshold, as discussed in Remazeilles et al. (2013): given both the number of channels and the number of spectral modes in ℓ-space at the needlet scale considered, our NILC algorithm consistently computes the number of spatial modes (similarly, the zone area for spatial localization) required to control the ILC bias.

3.1.2. MILCA

MILCA (Hurier et al. 2013), when applied to the extraction of the tSZ signal, also uses two spectral constraints: preservation of the tSZ signal (the tSZ spectral signature discussed above is assumed); and removal of the CMB contamination in the final SZ map, making use of the well known spectrum of the CMB. In addition, to compute the weights of the linear combination, we have used the extra degrees of freedom in the linear system to minimize residuals from other components (two degrees of freedom) and from the noise (two additional degrees). The noise covariance matrix is estimated from the half-difference maps described in Sect. 2.1. To improve the efficiency of the MILCA algorithm, weights are allowed to vary as a function of multipole ℓ, and are computed independently on different sky regions. We have used 11 filters in ℓ space as shown in Fig. 1 (dashed lines). These filters have an overall transmission of unity, except for ℓ< 8. For these large angular scales we have used a Gaussian filter to reduce foreground contamination. To ensure sufficient spatial localization for each required resolution, the size of the independent sky regions was adapted to the multipole range. We used a minimum of 12 regions at low resolution and a maximum of 3072 regions at high resolution.

3.2. Reconstructed Compton parameter map

Figure 2 shows the reconstructed Planck all-sky Compton parameter map for NILC (top panel) and MILCA (bottom panel). For display purposes, the maps are filtered using the procedure described in the first paragraph of Sect. 6.1. Clusters appear as positive sources, the Coma cluster and Virgo supercluster clearly visible near the north Galactic pole. The Galactic plane is masked in both maps, leaving 67% of the sky. Residual Galactic contamination is also visible as diffuse positive structures in the MILCA y-map. We can also observe a granular structure in the NILC y-map that corresponds to an excess of noise at large angular scales, as discussed in Sect. 4.2.

Weaker and more compact clusters are visible in the zoomed region of the Southern cap, shown in Fig. 3. Strong Galactic and extragalactic radio sources show up as negative bright spots on the maps and were masked prior to any scientific analysis, as discussed below in Sect. 4.4.1. We can again observe residual Galactic contamination around the edges of the masked area, which is more important for MILCA. Finally, we note in the NILC and MILCA y-maps the presence of systematic residuals along the scanning direction that show up as stripes. These are the consequence of stripes in the original Planck maps (see Planck Collaboration VIII 2016, for a detailed discussion). We discuss the level of residual stripes in Sect. 4.1.

In addition to the full Compton parameter maps, we also produce the so-called first and last (“F” and “L” hereafter) Compton parameter maps from the first and second halves of the survey rings (i.e., pointing periods). These maps are used for power spectrum estimation in Sect. 5 as well as to estimate the noise properties in the y-maps (see Sect. 4.2).

3.3. Comparison to other Compton parameter maps in the literature

Van Waerbeke et al. (2014) and Hill & Spergel (2014) have used the Planck nominal data to reconstruct a Compton parameter map over a large fraction of the sky. As for the present paper, they use an ILC approach, imposing spectral constraints to have unitary gain for the SZ component and null contribution from CMB, but they do not use the spectral and spatial separation discussed above. Van Waerbeke et al. (2014) consider only the four HFI channels between 100 and 353 GHz and force a null contribution from dust emission. Hill & Spergel (2014) also include the 545 GHz map, while using the 857 GHz channel only as a template for dust emission: a flux cut is imposed in order to build a mask that keeps only the 30% of the sky. They also apply a point source mask (radio and IR) before computing the ILC coefficients, reducing the final sky fraction to about the 25%. Furthermore, Van Waerbeke et al. (2014) and Hill & Spergel (2014) aim mainly at studying the cross-correlation with the gravitational lensing mass map from the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS) and the publicly-released Planck CMB lensing potential map, respectively. Their y-map reconstruction methods are tailored to fulfil these objectives and as a consequence they give larger overall foreground contamination, and in the case of Hill & Spergel (2014) also a significantly smaller sky fraction.

4. Pixel space analysis

4.1. Stripes in the y-map

As discussed in Sect. 3.2, residual stripes are visible in the NILC and MILCA y-maps. These stripes are mainly due to in-scan direction residual systematics after subtraction of an offset for each Planck stable pointing period (Planck Collaboration VIII 2016). To study how these stripes contaminate the final y-maps we have decomposed them in their in-scan and cross-scan direction contributions. We first convert the y-maps from Galactic to Ecliptic coordinates for which the scan direction corresponds mainly to a fixed longitude value. Secondly, we apply a Galactic mask to the Ecliptic y-maps and decompose them into spherical harmonics. Third, we select the in-scan (cross-scan) direction components by nullifying the spherical harmonic coefficients, a_ℓ,m, for | m | >ℓ/ 2 (| m | < = ℓ/ 2). Finally, we construct maps from those transformed coefficients and convert back to Galactic coordinates. We have chosen to mask the brightest 40% of the sky in the 857 GHz Planck map to keep Galactic ringing residuals negligible and a sufficiently large fraction of the sky for the analysis.

Figure 4 shows (from left to the right) the original, in-scan, and cross-scan y-maps for NILC (top) and MILCA (bottom) in the southern sky region presented in Fig. 3. We note that the y-maps look noisier since they are not filtered using the procedure described in the first paragraph of Sect. 6.1. The stripes are apparent in both the original and in-scan y-maps, as expected. We find that the ratio of the rms of the in- and cross-scan maps is 1.16 (1.17) for NILC (MILCA), consistent with residual stripe contamination (for Gaussian noise only maps we find a ratio of 1). Here we have measured an overall increase of the rms of the maps due to stripe contamination; a more detailed estimate of the effect in terms of tSZ power spectrum is presented later in Sect. 5.2.1.

The in-scan and cross-scan decomposition method modifies the cluster signal significantly due to ringing effects – see positive and negative patterns around clusters in Fig. 4. To explore this effect we have also applied the in-scan and cross-scan decomposition method to the simulated Compton parameter map for the detected and confirmed clusters of galaxies in the Planck catalogue (Planck Collaboration XXIX 2014), which is presented in Sect. 5.3. We find that those negative and positive patterns are also present in the decomposed maps and the ratio between the rms of the in-scan and cross-scan maps is 1, as expected. We finally stress that the in-scan and cross-scan decomposition of the y-maps does not preserve the positiveness of the tSZ signal. As a consequence we can not use them to estimate the contamination of the stripe systematic effect in the analysis presented in Sect. 6.

4.2. Noise distribution on the y-map

The Planck maps present a highly non-homogeneous structure for the noise that is mainly due to the inhomogeneous scanning strategy. Such complicated structure is propagated into the y-map and needs to be considered for further analysis. Here, we have chosen to describe the noise structure in the y-map by: (A) a variance map, which will capture the inhomogeneity of the noise; and (B) the angular power spectrum of a homogenous noise map, which is obtained after correcting for inhomogeneity. The variance map and homogeneous noise power spectrum have been obtained using two different methods, which we now describe.

• 1.

  Half-difference. The half-difference of the half-ring y-maps is used to obtain an estimate of the noise in the NILC and MILCA y-maps. From this half-difference map we compute an estimate of the noise variance at lower resolution by squaring and downgrading it to Healpix N_side = 128. The homogenous noise contribution is obtained by dividing the half-difference map by the variance map after upgrading it to Healpix N_side = 2048.

• 2.

  Simulations-based. 100 noise only realizations from the simulations described in Sect. 2.2 are used to compute an estimate of the variance per pixel at full resolution, which is then averaged to Healpix N_side = 128 resolution. The homogenous noise contribution map is then computed as above.

The two methods give consistent results at high Galactic latitudes, but differ at low Galactic latitudes, where we expect larger foreground residuals that show up in the half-difference method. This is visible in the top panel of Fig. 5 where we display variance maps obtained from the half-difference method for NILC (left) and MILCA (right). For a wide range of analyses with the y-map (e.g., computing the radial profile, surface density or total flux of a cluster, stacking of faint sources, detection of shocks, and study of the ICM) an estimate of the overall uncertainties per pixel (including foreground contribution) is required. Thus, the half-difference maps are released along with the y-map. We also show in the bottom panel of Fig. 5 the power spectrum of the homogeneous noise contribution. We observe a significant large angular scale (low multipole) component in the homogeneous noise contribution. This low multipole component is significantly larger for the NILC y-map than for the MILCA one.

4.3. Foreground contamination and masking

As discussed above, a residual Galactic foreground contribution is observed at low Galactic latitudes in both the NILC and MILCA y-maps, while being more important for the MILCA y-map (see Sect. 5.2.1 for a more quantitative comparison). This contribution is mainly induced by thermal dust emission, as discussed in detail in Planck Collaboration XXIX (2014), and is therefore mainly associated with the Galactic plane. When dealing with individual objects, or with the stacking of faint objects, we recommend that it is account for by using the variance maps discussed above. For statistical analyses, such as those presented in Sects. 5 and 6, specific masks need to be defined and these are fully discussed in those sections.

4.4. tSZ signal from resolved sources

4.4.1. Point sources

Point source contamination is an important issue for the cosmological interpretation of the Planck Compton parameter map.

In the reconstructed tSZ maps radio sources will appear as negative peaks, while infrared sources will show up as positive peaks, mimicking the signal from clusters. To avoid contamination from these sources we introduce a point source mask (PSMASK, hereafter). This mask is the union of the individual frequency point-source masks discussed in Planck Collaboration XXVIII (2014). The reliability of this mask was verified by looking at the 1D Probability Density Function (PDF) of the pixel signal (as will be discussed in Sect. 6.1). We found that a more robust mask can be obtained by enlarging the mask size around the strongest radio sources in order to minimize their contribution.

Alternatively, we have also performed a blind search for negative sources in the y-maps using the MHW2 algorithm (López-Caniego et al. 2006). We have detected 997 and 992 negative sources for NILC and MILCA respectively. We find that 67 and 42 (for NILC and MILCA, respectively) of those detected sources are not masked by the PSMASK. However, most of them (54 for for NILC and 36 for MILCA) are false detections made by the algorithm in regions surrounding very strong positive sources (i.e., galaxy clusters). For those detected as additional negative sources, the PSMASK has been updated to account for them, even though the strongest are already excluded when applying the 50% Galactic mask used for the cosmological analysis below.

For infrared sources, estimating the efficiency of the masking is hampered by the tSZ signal itself. The residual contamination from point sources is discussed in Sects. 5.2 and 6. It is also important to note that the PSMASK may exclude some clusters of galaxies. This is particularly true in the case of clusters with strong central radio sources, such as the Perseus cluster (see Planck Collaboration XXIX 2014).

4.4.2. Blind search for clusters

Following Planck Collaboration XXI (2014) a blind search for tSZ (positive) sources has alsobeen performed on the all-sky NILC and MILCA y-maps. We use the IFCAMEX (MHW2, González-Nuevo et al. 2006; López-Caniego et al. 2006) and the single frequency matched filter (MF, Melin et al. 2006) methods. The detected sources with signal-to-noise ratio S/N> 4 have been matched to the Planck cluster catalogue (PSZ2, Planck Collaboration XXVII 2016). We have associated the detected sources to PSZ2 objects if the distance between their positions is smaller than 10′ (the resolution of the SZ all-sky maps). The MHW2 algorithm finds 1018 and 1522 candidates for NILC and MILCA, respectively, out of which 500 and 457 correspond to objects present within the PSZ2 catalogue. For NILC (MILCA) we have 41 (25) positions that are not masked by the PSMASK and that do not correspond to PSZ2 clusters. However 14 (5) are located in regions excluded by the mask used to build the PSZ2 catalogue. With the MF method we have 1472 and 1502 candidates for NILC and MILCA, out of which 1107 and 1096 correspond to objects present within the PSZ2 catalogue (867 and 835 corresponding to already validated objects, respectively). This implies good agreement between the Planck cluster sample and the detected sources in the y-maps.

This agreement is improved by taking the union of the MHW2 and MF catalogues. In this case, and considering only validated sources in the PSZ2 catalogue (1070 sources), we find 907 and 870 matches with blind detected cluster candidates for NILC and MILCA, respectively. We have also performed a visual inspection of the NILC and MILCA y-maps at the positions of validated PSZ2 sources for which we find no counterpart, and here we find evidence of an excess of signal. For most of these sources low S/N and/or foreground contamination can explain why they are not detected blindly in the y-map. This is consistent with the fact that we expect blind detection methods in the y-map to be less sensitive than multifrequency searches for compact sources (Melin et al. 2012). In some cases these non-detected sources are extended but with a relatively high S/N (between 6 and 12). We note that the MHW2 and MF methods are tuned to detect point-like and compact sources primarily, and thus they may fail for extended ones.

We present in the top panel of Fig. 6 a comparison of the flux measured for these common sources in the NILC y-map with that derived from the PSZ2. The quantity Y_5r500 corresponds to the integrated signal within a radius equal to 5 ×r₅₀₀³. We observe good agreement between the two. Most of the observed outliers (low y-map flux with respect to the PSZ2 one) correspond to cluster candidates detected close to masked point sources, or in regions with strong Galactic contamination. We also show in the bottom panel a comparison of the flux measured with the NILC and MILCA y-maps. We find good agreement between the two, with three outliers corresponding to sources lying either in a region badly affected by radio and/or infrared point sources or for which we observe large uncertainties in the estimated characteristic radius of the cluster candidate. For both maps (NILC and MILCA), the average S/N of the recovered clusters with the MHW2 algorithm is around 10, while S/N of PSZ2 clusters only found with the matched filter technique is, as expected, lower (about 6).

4.4.3. Maps of selected clusters

We have performed a more detailed analysis for some of the cluster candidates in the PSZ2 sample. In Fig. 7 we present Compton parameter maps centred at the positions of three of the newly discovered Planck clusters with S/N of 9.3 (left), 6.2 (middle) and 4.6 (right), both for MILCA (top row) and NILC (middle row). In the bottom row of the figure we present the radial profiles of these clusters as obtained from the MILCA (black lines) and NILC (blue) maps. We find that the profiles are consistent between MILCA and NILC. For comparison we also show the radial profile of a 10′ Gaussian beam. We observe that even for low S/N and compact clusters we are able to detect extended emission. Thus, these maps could be used to expand the pressure profile analysis presented in Planck Collaboration Int. V (2013) to fainter and higher redshift clusters.

One of the major outcomes of the NILC and MILCA y-maps is the possibility to study diffuse faint emission between clusters, as well as the emission from the cluster outskirts. We show in Fig. 8 some well-known merging systems including the Shapley supercluster, and the A3395-A3391 and the A339-A401 interacting systems (see Planck Collaboration Int. VIII 2013, for a detailed description). We observe that the clusters themselves, their outskirts, and the inter-cluster emission are well reconstructed in the NILC and MILCA y-maps (which show consistent results). The quality of these maps will permit analyses similar to what is presented in Planck Collaboration Int. VIII (2013), with a significant increase of the S/N.

4.5. SZ signal below the noise level

We use the catalog of 132 684 clusters of galaxies identified from the Sloan Digital Sky Survey III (Wen et al. 2012) in order to quantify the tSZ signal below the noise level in the y-map. This catalogue provides estimates of cluster redshift, richness (N₂₀₀), and characteristic radius (r₂₀₀).

We focus on unresolved groups and clusters for which the richness is 8 ≤ N₂₀₀ ≤ 100. These objects are expected to be at S/N well below 1 in the y-map and so their direct detection is not possible. However, their average tSZ signal can be detected using a stacking approach on the y-maps. Figure 9 shows the stacked maps (obtained with the IAS stacking library, Bavouzet 2008; Béthermin et al. 2010) for six richness intervals with N₂₀₀ varying from 8 to 100, on patches of 4° × 4°, for the full sky NILC y-map. For the stacking we exclude all the positions for which point sources are found within a radius of 10′ from the centre of the object. The noise of the stacked maps scales, as expected, with , and the average signal increases for increasing richness. Then, when considering a sufficiently large number of objects, we are able to significantly detect the average SZ signal, even for small groups (Fig. 9). We obtain consistent results for the MILCA and NILC y-maps.

For the N₂₀₀ intervals reported in Fig. 9 we have estimated the total stacked fluxes (Y_5r500) as the integrated signal within a radius 5 × r₅₀₀, with r₅₀₀ = 0.7 ⟨ r₂₀₀ ⟩. Here ⟨ r₂₀₀ ⟩ is the mean of the r₂₀₀ reported by Wen et al. (2012) for clusters belonging to each considered subsample. The average values and associated errors have been obtained using a bootstrap approach. For this, we have constructed and stacked cluster samples obtained by randomly replacing sources from the original sample, so that each of them contains a number of clusters equal to the initial one. By repeating the process several times, we have determined the statistical properties of the population being stacked. In Fig. 10 (left panel) we report Y_5r500 as a function of richness (N₂₀₀) for all the objects with z ≤ 0.42. In fact higher redshift clusters present in this catalogue may be biased to lower richness because of incompleteness of member galaxies, as detailed in Wen et al. (2012). Following Planck Collaboration XII (2011), we have fitted a power law of the form (4)with E(z)² = (H(z) /H₀)², D_A the angular diameter distance of the cluster, and Y₂₀ a normalization parameter.

By considering only N₂₀₀ ≤ 60 (for which the number of clusters in each richness bin is >1000), we find a slope α = 1.85 ± 0.20, which is consistent with the scaling obtained in Planck Collaboration XII (2011). We have verified that the scaling obtained is insensitive to whether we limit our cluster sample to high Galactic latitude objects or not, and that a change in the integration radius only affects the normalization Y₂₀. Using Y₅₀₀ = Y_5r500/ 1.79 (Planck Collaboration XXIX 2014), we obtain Y₂₀ = (16.5 ± 1.2) × 10^-5 arcmin². However, for independent cluster data sets, different choices are made for the fainter magnitude of the member galaxies contributing to N₂₀₀. This may affect the richness associated with an object of given mass and complicate the comparison of the constant term Y₂₀ for analyses that are based on different cluster data sets.

The mass-richness scaling relation, (5)can then be used to also derive the Y₅₀₀–M_tot scaling (see Fig. 10, right panel), by fitting a power law of the form (6)Ford et al. (2014) compares different results for the scaling of the mass with richness, in the form of Eq. (5). In particular they discuss the results obtained by Wen et al. (2012), for a subsample of clusters with already known masses, and the results of Covone et al. (2014), from weak lensing mass measurements of 1176 clusters of the catalogue used here. Assuming β = 1.2 ± 0.1 and M₀ = (1.1 ± 0.11) × 10¹⁴M_⊙ (Wen et al. 2012; Covone et al. 2014,and a 10% error) in Eq. (5), we compute via a Monte Carlo approach the average M₂₀₀ and associated uncertainties for each Y_5r500 bin. From this we find B = 1.92 ± 0.42 for the Y_5r500−M₂₀₀ scaling (Eq. (6)). This is consistent with the self-similar value (5/3) and with what is expected for this mass range according to simulations (Sembolini et al. 2014). Through a weak lensing analysis of more than 18 000 clusters in the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS), Ford et al. (2014) obtain a steeper calibration for the mass-richness relation, with β = 1.58 ± 0.34. For higher β, the Y_5r500−M₂₀₀ slope B becomes lower, but still consistent with the value obtained before. Again comparison of M₀ and Y₀ is complicated by the different choices of the fainter limit on galaxies contributing to N₂₀₀. By exploring the range 0.2 ≲ z ≲ 0.9, Ford et al. (2014) finds no evolution with redshift. Thus, we do not expect our result to be biased by the z ≤ 0.42 cut.

5. Angular power spectrum

We now consider the validation of the Compton parameter maps at the power spectrum level.

5.1. Power spectrum estimation

To estimate the power spectrum of the tSZ signal we use the XSPECT method (Tristram et al. 2005), initially developed for the cross-correlation of independent detector maps. XSPECT uses standard MASTER-like techniques (Hivon et al. 2002) to correct for the beam convolution and the pixelization, as well as the mode-coupling induced by masking foreground contaminated sky regions.

In the following, all the spectra will use a common multipole binning scheme, which was defined in order to minimize the correlation between adjacent bins at low multipoles and to increase the S/N at high multipole values. We will also only consider angular cross-power spectra between the y-maps reconstructed from the first (F) and second (L) halves of the data. This allows us to avoid the bias induced by the noise in the angular auto power spectrum, the correction of which would require a large number of noise simulations. The angular cross-spectra are named in the following “NILC F/L” and “MILCA F/L”. Error bars in the spectrum are computed analytically from the auto-power and cross-power spectra of the pairs of maps, as described in Tristram et al. (2005). We do not consider here higher order terms in the power spectrum variance. All of our Compton parameter maps assume a circular Gaussian beam of 10′ FWHM. The additional filtering at large angular scales in the MILCA Compton parameter maps is also accounted for and deconvolved when producing power spectra.

5.2. Foreground contamination

In Planck Collaboration XXI (2014) we identified the dominant foreground contributions to the angular power spectrum of the reconstructed y-maps using the FFP6 simulated data set. We have repeated here the same analysis on updated simulated data sets (see Sect. 2.2), for which the description of foreground components has been improved according to our current knowledge (Planck Collaboration XII 2016). We find no major difference between the two analyses and therefore we do not repeat the discussion here. At large angular scales the dominant foreground contribution is the Galactic thermal dust emission, as we have already discussed in Sect. 4.3. At intermediate and small angular scales the major contribution comes from the clustered CIB emission. Equally important at small angular scales are the residual contributions from radio and IR sources.

5.2.1. Low-multipole contribution

Assuming that at large angular scales the Compton parameter maps are mainly affected by diffuse Galactic dust emission, we have tested several Galactic masks by imposing flux cuts on the Planck 857 GHz channel intensity map. In particular we investigated masking out 30%, 40%, 50%, and 60% of the sky. The edges of these masks have been apodized to limit ringing effects on the reconstruction of the angular power spectrum. We stress that to compute power spectra we used the XSPECT algorithm (Tristram et al. 2005), which has been specifically developed to account for mode coupling. We have proved using simulations that we have no bias in the final power spectrum, and so we do not expect that any differences between MILCA and NILC could be induced by the different sizes of the masked area. Figure 11 shows the power spectra for NILC (left) and MILCA (right). We observe that the MILCA F/L large angular scale power decreases when imposing a more severe masking (larger fraction of the sky is masked out). A similar effect, but significantly smaller, is also observed for the MILCA F/L cross–power spectrum. The MILCA y-map is significantly more contaminated by thermal dust emission at multipoles below 30. However, we also observe for NILC F/L an extra noise contribution at large angular scales, as discussed in Sect. 4.2. These two problems limit the reliability of the results at multipoles below 30.

To extend the measurement of the angular power spectrum of the tSZ emission to multipoles below 30 we have considered two options: (1) as in Planck Collaboration XXI (2014) we apply a more severe Galactic mask (30% of the sky is masked) before the computation of the NILC weights to produce the y-map (NILC70); and (2) compute the cross-correlation of NILC and MILCA y-maps. Considering 50% of the sky, we show in Fig. 11 the angular power spectrum of the NILC70 y-map (black) and the cross-spectrum of the NILC first half and MILCA second half (NILC-MILCA F/L) y-maps. We observe that the two are compatible within error bars in the multipole range 10 <ℓ< 1500. As expected the cross spectrum NILC-MILCA F/L shows larger error bars.

Although the NILC70 y-map seems to be the best choice in terms of power spectrum estimation, it results in a significant reduction of the available sky area for other kind of studies as those presented in Sect. 4. Furthermore, it is difficult to accurately estimate the ultimate residual foreground contribution at very large angular scales. Because of this, and to preserve the coherence of the delivered products and the analysis presented in this paper, we have chosen the angular cross–power spectrum of NILC-MILCA F/L as a baseline. This is obviously a more conservative choice in terms of noise–induced uncertainties. The NILC-MILCA F/L angular cross–power spectrum bandpowers and uncertainties are further discussed in Sect. 7. Using the in-scan and cross-scan y-maps presented in Sect. 4.1 we find that stripe contamination accounts for less than 10% of the total signal in the NILC-MILCA F/L cross angular power spectrum. We have enlarged the error bars to account for this systematic effect.

5.2.2. High-multipole contribution

At small angular scales the measured tSZ power spectrum is affected by residual foreground contamination, coming from clustered CIB emission as well as radio and IR point sources. They show up in the MILCA-NILC F/L cross–power spectrum (see Fig. 11) as an excess of power at large multipoles.

To deal with those residuals we adopt the same strategy as in Planck Collaboration XXI (2014). We define physically motivated models of the angular power spectrum of the foreground components for each observation channel, including cross–correlations between channels. In contrast to Planck Collaboration XXI (2014) we also account for the cross correlation between the clustered CIB and the tSZ emission. A detailed description of this cross correlation, as well as of the clustered CIB model is presented in Planck Collaboration XXIII (2016). For the radio and IR point source models we refer to Planck Collaboration XXI (2014).

We use the outputs of Planck Collaboration XVIII (2011) and Planck Collaboration XXX (2014) for the clustered CIB modelling. For the six Planck HFI frequencies considered in this paper, the clustered CIB model consists of six auto-power spectra and 24 cross-power spectra. For frequencies above 217 GHz, these spectra are fitted in Planck Collaboration XXX (2014) to the measured CIB, consistently with Planck Collaboration XVIII (2011). The model is extrapolated at 100 and 143 GHz following Béthermin et al. (2012) and Planck Collaboration XVIII (2011). The uncertainties in the clustered-CIB model are mainly due to the cross-correlation coefficients that relate the cross-power spectra to the auto-power spectra. Following Planck Collaboration XXX (2014) we consider 5% global uncertainties on those coefficients.

We use the Béthermin et al. (2012) model to compute the star-forming dusty galaxy contribution. Finally, we use the Tucci et al. (2011) model, fitted to the Planck ERCSC (Planck Collaboration Int. VII 2013), for extragalactic radio sources. Notice that these models are also used for the study of the clustered CIB with Planck (Planck Collaboration XXX 2014).

Using the models described above we generate Gaussian realizations of the foreground contribution for each HFI frequency channel between 100 and 857 GHz. Note that the LFI channels are only used at large angular scales. We apply the MILCA or NILC weights to these simulated maps. From these simulations we find that the cross–correlation between the CIB and tSZ contribution can be neglected to first order with respect to the others and will therefore not be considered hereafter. Uncertainties on the parameters describing the foreground models were also propagated using simulations. We find that the clustered CIB model uncertainties might be as large as 50% in amplitude. In addition, we notice that the amplitude of the point source models can vary significantly depending on the point source mask applied. These uncertainties are taken into account hereafter. The amplitude of the residual foreground models are jointly fitted with the cosmology-dependent tSZ model, as detailed in Sect. 7.1.

5.3. Contribution of resolved clusters to the tSZ power spectrum

We simulate the expected Compton parameter map for the detected and confirmed clusters of galaxies in the Planck catalogue (Planck Collaboration XXIX 2014) from their measured integrated Compton parameter Y_5r500 and redshift z. We assume hydrostatic equilibrium and an Arnaud et al. (2010) pressure profile. The green solid line in Fig. 12 shows the power spectrum of this simulated map. Figure 12 also shows the cross-power spectra of the NILC and MILCA F/L maps (in black). We also compute the cross-power spectrum of the simulated cluster map and the Planck reconstructed NILC and MILCA Compton parameter maps. This is shown in red in the figure. Here again, the signal is consistent with the expected power spectrum of the confirmed Planck clusters of galaxies. These results show that the tSZ signal from clusters is preserved in the y-maps.

6. Higher–order statistics

The power spectrum analysis presented above only provides information on the 2-point statistics of the Compton parameter distribution over the sky. An extended characterization of the field can be performed by studying the higher-order moments in the 1D PDF of the map, or by measuring 3-point statistics, i.e., the bispectrum.

6.1. 1D PDF analysis

Following Planck Collaboration XXI (2014), we perform an analysis of the 1D PDF of the NILC and MILCA reconstructed Compton parameter maps. For the tSZ effect we expect an asymmetric distribution with a significantly positive tail (Rubiño-Martín & Sunyaev 2003). We thus focus on the asymmetry of the distribution and its unnormalized skewness. First, we filter the maps in order to enhance the tSZ signal with respect to foreground contamination and noise. We then apodize the combined PSMASK and Galactic mask to avoid residual point source ringing. We follow the approach of Wilson et al. (2012) and use a filter in harmonic space. For each multipole ℓ, the filter is constructed as the ratio between the angular power spectrum of the expected tSZ signal (, obtained from the simulations in Sect. 2.2) and the power spectrum of the half-difference y maps (, see Sect. 4.2), such that . We smooth this filter function in multipole space using a 21-point square kernel and normalize it to one. Notice that this filter only selects the multipole range for which the tSZ signal is large with respect to the noise, and thus, it does not modify significantly the non-Gaussianity properties of the y-maps. Furthermore, we have found that the filter used here behaves better than the more traditionally used Wiener filter, because it is less affected by point-source ringing. Following this procedure, the 1D PDF of the filtered Compton parameter maps, P(y), is computed from the histogram of the pixels. As for the power spectrum, several Galactic masks have been considered in order to tests the robustness of the results.

Figure 13 shows the 1D PDF for the NILC (left) and MILCA (right) Compton parameter map in red when masking 50% of the sky. They correspond to the convolution of the 1D PDF of the different components in the map: the tSZ effect; foregrounds; and the noise. The 1D PDFs of the NILC and MILCA y-maps clearly show three distinct contributions: a Gaussian central part that is slightly wider than the noise contribution, as expected from the half-difference map 1D PDF (black curve); a small negative tail, corresponding most likely to residual radio sources; and a positive tail corresponding mainly to the tSZ signal as observed for the PSZ2 cluster simulated map (dark blue). In comparison to Planck Collaboration XXI (2014) we now find a better agreement between the NILC and MILCA results. This is mainly due to the improved masking of radio sources presented in Sect. 4.4.1. Finally we also show the 1D PDF of the reconstructed y-maps after masking the PSZ2 clusters (orange). We observe that most of the tSZ is removed, indicating that the 1D PDF of the y-map is dominated by detected clusters.

A simple analysis of the measured 1D PDF can be performed by considering the asymmetry of the distribution: (7)where y_p is the peak value of the normalized PDF P(y). In addition, the non-Gaussianity of the positive tail can be quantified by (8)with G(y) the expected PDF if fluctuations were only due to noise. The latter can be obtained from the half-difference y-maps. See Hill et al. (2014) for a similar analysis conducted on ACT data.

Masking 60% of the sky, we find A = 0.218 and Δ = 0.10 for the NILC Compton parameter map. Equivalently, for the MILCA Compton parameter map we find A = 0.223 and Δ = 0.11. Uncertainties on A and Δ are mainly given by the accuracy to which measure the Gaussian noise distribution and are below 1 %. These results are in agreement and consistent with a positive tail in the 1D PDF, confirming the tSZ nature of the signal. Actually, we find similar values for the FFP6 simulations, A = 0.30 and Δ = 0.13 as presented in Planck Collaboration XX (2014). The agreement between the NILC and MILCA results worsens when reducing the masked area, as a consequence of a larger foreground contribution in the MILCA y-map.

6.2. Bispectrum

Since the SZ signal is non-Gaussian, significant statistical information is contained in the bispectrum, complementary to the power spectrum (Rubiño-Martín & Sunyaev 2003; Bhattacharya et al. 2012). Results on SPT data have been obtained by Crawford et al. (2014) and George et al. (2015). The bispectrum also provides an additional statistic assessing the compatibility of the NILC and MILCA reconstructed Compton parameter maps, as well as their reliability in terms of foreground contamination.

We therefore analyse the bispectra of the NILC and MILCA maps. The estimation method is essentially the same as for the 2013 results (Planck Collaboration XX 2014), briefly recapped here. We use the binned bispectrum estimator described in Bucher et al. (2010) and Lacasa et al. (2012), which is also used for the Planck primordial non-Gaussianity analysis (Planck Collaboration XXIV 2014; Planck Collaboration XVII 2016). We use a multipole bin size Δℓ = 64 and a maximum multipole ℓ_max = 2048 for the analysis, working at a resolution N_side = 1024 to reduce computing time. The NILC and MILCA maps are masked with the combination of PSMASK, described in Sect. 4.4.1, and a Galactic mask at 50, 60, or 70%, described in Sect. 5.2.1 (in the rest of this section we will simply denote as X% mask the combination of PSMASK and the Galactic mask at X%). The best-fit monopole and dipole outside the mask are finally removed before estimation.

An important part of the pipeline is then to correct for the bias introduced by masking. To this end, we compute the ratio of the full-sky and masked sky bispectra, on highly non-Gaussian simulations with a tSZ-like bispectrum and 10′ resolution. This ratio is used to correct the measured bispectra and flag unreliable (ℓ₁,ℓ₂,ℓ₃) configurations. Specifically we flag configurations where the ratio is different by more than 25% from the naive expectation f_skyB(ℓ₁) B(ℓ₂) B(ℓ₃), where B(ℓ) is the Gaussian 10′ beam.

For both NILC and MILCA, we find that the bispectra computed on the 50, 60 and 70% masks are consistent. This indicates that there is no detectable residual galactic contamination in these bispectra. However we did find Galactic contamination on less aggressive Galactic masks (when we mask less than 50% of the sky), specifically positive Galactic dust. As Galactic dust is highly non-Gaussian, we do not recommend the measurement of higher order statistics using Galactic masks smaller than 50%. In the following we adopt the 50% mask as a baseline, since it leaves the most sky available for estimation and minimizes masking effects in the measurement.

Figure 14 shows the obtained bispectra as a function of multipole for the NILC and MILCA Compton parameter maps. We observe good agreement between the bispectra of the two maps, and the bispectral behaviour is consistent with that expected from a tSZ signal (see, e.g., Lacasa 2014, chapter 5). We further compare these measurements with the bispectrum of the simulated map for the PSZ2 clusters, which is presented in blue in Fig. 14. We observe good agreement between the bispectra of the NILC and MILCA and that of the PSZ2 clusters. We therefore conclude that the observed bispectrum in the y-map is dominated by detected clusters.

Finally, in Fig. 14 are shown the ± 1σ uncertainties of the measurements, in black dotted lines. The error bars were computed in a similar manner to that of the 2013 results (Planck Collaboration XX 2014), and Appendix A.3 gives more detailed discussion.

With a detection per configuration at an average significance of 3.5σ, and a total significance of around 60σ, the Planck data thus provide a high quality measurement of the non-Gaussianity of the thermal Sunyaev-Zel’dovich signal, with contamination from foregrounds below the uncertainty levels when masking more than 50 % of the sky.

7. Cluster physics and cosmology

7.1. Power spectrum analysis

7.1.1. tSZ power spectrum modelling

As a measure of structure growth, the tSZ power spectrum can provide independent constraints on cosmological parameters. As shown by Komatsu & Seljak (2002), the power spectrum of the tSZ effect is highly sensitive to the normalization of the matter power spectrum, commonly parameterized by the rms of the z = 0 mass distribution on 8 h^-1 Mpc scales, σ₈, and to the total amount of matter Ω_m. We expect the tSZ power spectrum to also be sensitive to other cosmological parameters, e.g., the baryon density parameter Ω_b, the Hubble contant H₀, and the primordial special index n_s. For reasonable external priors on those parameters, however, the variations are expected to be negligible with respect to those introduced by changes in Ω_m and σ₈ and so they are not considered here.

Following Planck Collaboration XXI (2014) we consider here a two-halo model for the tSZ power spectrum, which is fully described in Appendix A.1. This model accounts for both intra-halo (1-halo term) and inter-halo (two-halos) correlations. Following Eq. (A.6), the tSZ spectrum is computed using the 2-halo model, the Tinker et al. (2008) mass function, and the Arnaud et al. (2010) pressure profile. In particular, we use the numerical implementation presented in Taburet et al. (2009, 2010, 2011), and integrate in redshift from 0 to 3 and in mass from 10¹³M_⊙ to 5 × 10¹⁵M_⊙. Our model allows us to compute the tSZ power spectrum at the largest angular scales and is consistent with the tSZ spectrum presented in Efstathiou & Migliaccio 2012 (EM12), which was used as a template in the CMB cosmological analysis in Planck Collaboration XVI (2014) and Planck Collaboration XI (2016). We also include the mass bias parameter, b, which accounts for bias between the observationally deduced (M^obs) and true (M^true) mass of the cluster (see Planck Collaboration XX 2014, for details), such that M^obs = (1−b)M^true.

Cosmological parameter results are very sensitive to the mass bias and in particular we expect σ₈ and Ω_m to be strongly degenerate with b (Planck Collaboration XXI 2014). By contrast to Planck Collaboration XXIV (2016), for which external priors in the mass bias have been used, we consider here only two distinct values, b = 0.2 and b = 0.4. The former corresponds to the average value that numerical simulations seem to indicate (see Fig. A.2 in Planck Collaboration XX 2014). On the other hand, the latter value for the bias alleviates the inconsistency with the constraints derived from the analysis based on primary CMB anisotropies (see Planck Collaboration XX 2014).

This value is, however, larger than that obtained by mass comparison of clusters present in both the Planck cosmology sample (Planck Collaboration XX 2014) and the Weighing the Giants (WtG von der Linden et al. 2014) project. Even if studies based on lensing mass measurements still provide different and inconsistent results for the cluster mass calibration, their number and their accuracy has improved dramatically in the recent past (e.g., Mahdavi et al. 2013, for the CCCP project, Umetsu et al. 2014, for CLASH; Israel et al. 2014, for 400d WL; Ford et al. 2014 for CFHTLenS; Gruen et al. 2014 for WISCy) and they are expected to provide more useful constraints in the near future.

7.1.2. Maximum likelihood analysis

As in Planck Collaboration XXI (2014), cosmological constraints are obtained from a fit to the tSZ power spectrum. Following the discussion of Sect. 5.2 we take the NILC-MILCA F/L cross-power spectrum (black dots in Fig. 15) as a reference and limit the analysis to 50% of the sky to minimize foreground residuals. In terms of astrophysical signal we consider a four-component model: tSZ; clustered CIB; radio point sources; and infrared point sources. We restrict the analysis to multipoles ℓ> 10 so that we can neglect the residual thermal dust contamination (see Sect. 5.2.1). For ℓ> 2000 the total signal in the tSZ map is dominated by correlated noise, which is also accounted for in the fit. Because of this correlated noise and the expected high value of foreground contamination at the smallest scales we limit the fit to multipoles ℓ< 1411. Finally, the observed y-map power spectrum, , is modelled as (9)Here is the tSZ power spectrum (in red in Fig. 15), is the clustered CIB power spectrum (in green), and and are the infrared and radio source power spectra, respectively (in cyan and in blue). is an empirical model for the high multipole correlated noise.

Foreground contamination is modelled following Sect. 5.2.2. As discussed there, the main uncertainties in the residual power spectrum translate into up to 50% uncertainty in the amplitude of the clustered CIB. We have not considered in this analysis the CIB-tSZ cross-correlation that was shown to be negligible in terms of the power spectrum. The amplitude of the IR and radio point-source contribution will depend on the exact Galactic mask used for the analysis. However, we expect the shape of the their power spectra to remain the same. We thus allow for a variation of the normalization amplitudes for the clustered CIB, A_CIB, and for the point sources, A_IR and A_rad.

We assume a Gaussian approximation for the likelihood function. Best-fit values and uncertainties are obtained using an adapted version of the CosmoMC algorithm (Lewis & Bridle 2002). Only σ₈ and Ω_m are allowed to vary. All other cosmological parameters are fixed to their best-fit values, as obtained in Table 2 of Planck Collaboration XVI (2014). The normalization amplitudes, A_CIB, A_rad, and A_IR, considered as nuisance parameters, are allowed to vary between 0 and 3. For the range of multipoles considered here, the tSZ angular power spectrum varies approximately as . The results are thus presented in terms of this parameter combination.

7.2. Best-fit parameters and tSZ power spectrum

Figure 16 presents the 2D and 1D likelihood distributions for the cosmological parameter combination σ₈(Ω_m/ 0.28)^{3 / 8} and for the foreground nuisance parameters. We present the results obtained assuming a mass bias of 0.2 (black) and 0.4 (red). We obtain very similar values for the nuisance parameters in both cases. In particular the best-fit values for a mass bias of 0.2 are , , and . However, there is a significant shift in the value of σ₈(Ω_m/ 0.28)^{3 / 8} as one would expect (Planck Collaboration XX 2014). In the case of a mass bias of 0.2 we have , while for a mass bias of 0.4 we have . Notice that these values are obtained in a specific framework, with all other cosmological parameters being fixed and a fiducial fixed model used for the signals. Relaxing this framework would likely weaken the constraints presented here, as discussed below.

Figure 15 shows the NILC-MILCA F/L angular cross-power spectrum before correcting (black dots) for foreground contribution. We also show the best-fit foreground models: clustered CIB (green line); and radio (blue line) and IR (cyan line) point sources. The statistical (thick line) and total (statistical plus foreground, thin line) are also shown. The best-fit tSZ power spectrum is presented as a solid red line. We conclude that the NILC-MILCA F/L angular cross-power spectrum is dominated by tSZ for multipoles ℓ< 700, and by foreground contribution for multipoles ℓ> 1200. We also note that for the best-fit model the radio point-sources contribution seems to be negligible with respect to the IR one. This is not a realistic result and it is most probably explained by the strong degeneracy observed between the radio and IR point-source amplitude (see Fig. 16).

Finally we present in Fig. 17 the NILC-MILCA F/L angular cross-power spectrum after correcting for foreground contributions. Uncertainties account for statistical and systematic errors, as well as for uncertainties in the foreground subtraction. The marginalized bandpowers and uncertainties are also presented in Table 2. We note that foreground-induced uncertainties dominate at multipoles ℓ> 100. Bandpowers for the best-fit model for the angular tSZ power spectrum are also given for comparison. We also show in Fig. 17 th tSZ power spectrum estimates at high multipoles obtained in CMB oriented analyses by the Atacama Cosmology Telescope (ACT) and the South Pole Telescope (SPT, George et al. 2015). The black line shows the tSZ power spectrum template (EM12, Efstathiou & Migliaccio 2012) used in the Planck CMB cosmological analysis (Planck Collaboration XVI 2014; Planck Collaboration XI 2016) assuming the best-fit amplitude A_tSZ and the grey region ±2σ uncertainties from Planck Collaboration XI (2016). A direct comparison of our results to the estimates of the SZ power spectrum amplitude by ACT, SPT and Planck Collaboration XI (2016) is difficult as we need on the one hand to account for foreground contamination in the Planck SZ data and on the other hand to marginalise over all uncertainties on the SZ power spectrum amplitude. Accounting for foreground uncertainties as in Table 2 we find that χ² of the best-fit value is not modified significantly by including the ACT and SPT data.

7.2.1. Dependence on cluster physics

As discussed in Planck Collaboration XXI (2014), we also expect the tSZ power spectrum amplitude to be sensitive to the physics of clusters of galaxies. To explore this dependence we have considered a set of predicted tSZ spectra for various physical models. In Fig. 18 we compare these models to the foreground-cleaned Planck tSZ power spectrum derived above, as well as to the Atacama Cosmology Telescope (ACT) and the South Pole Telescope (SPT, George et al. 2015) power spectrum estimates. We consider the predictions derived from hydrodynamical simulations (Battaglia et al. 2010, 2012), from N-body simulations plus semi-analytical models (Trac et al. 2011, TBO2) and from analytical calculations (Shaw et al. 2010). These models were originally computed for the set of cosmological parameters in Hinshaw et al. (2012) with σ₈ = 0.8 and have been rescaled in amplitude to our best-fit value for . We note that there is some dispersion in the predicted amplitudes and shapes of the tSZ power spectrum. These differences reflect the range of methodologies and assumptions used both in the physical properties of clusters and in the technical details of the computation. The latter includes differences in the redshift ranges and also in the mass intervals probed by the limited sizes of the simulation boxes of the hydrodynamical simulations. Analytical predictions are also sensitive to the model ingredients, such as the mass function, mass bias and scaling relations adopted.

We see from Fig. 18 that the models presented above (the tSZ template for CMB analyses, plus the Battaglia et al. 2012; Shaw et al. 2010; and TBO2 models) provide reasonable fits to the data for multipoles above 200. For lower multipoles the Shaw et al. 2010 and TBO2 models are not consistent with the data.

We have also performed a simplified likelihood analysis to evaluate the uncertainties on cosmological parameters induced by the uncertainties in the modelling of the cluster physics. We replace our own model of the tSZ power spectrum by the models discussed above and recompute σ₈(Ω_m/ 0.28)^{3 / 8}, A_CIB, A_rad and A_IR from a simple linear fit to the NILC-MILCA F/L cross-power spectrum. In the case of a mass bias of 0.2, we obtain values for σ₈(Ω_m/ 0.28)^{3 / 8} between 0.77 and 0.80, which lie within the 1σ uncertainties (0.03) presented above.

In the case of our fiducial model (see Appendix A.1), we can also consider uncertainties in the parameters describing the scaling relations allowing us to relate the observed tSZ flux to the mass of the cluster for a given redshift. Following Eq. (7) in Planck Collaboration XXVIII (2014) the main parameters to be considered are the mass bias b, the overall amplitude Y_∗ and the scaling slope β. As discussed above, the mass bias is fully degenerate with σ₈. Similar conclusions can be drawn for Y_∗, which is expected to be known at the percent level (see Table 1 in Planck Collaboration XXVIII 2014) and therefore it is subdominant with respect to the uncertainties in the mass bias. Although the uncertainties in the slope of the scaling relation are relatively large, we have checked that they lead to negligible uncertainties on cosmological parameters.

7.3. Higher-order statistics

7.3.1. Skewness measurements

The skewness of the 1D PDF distribution, / can also be used to derive constraints on cosmological parameters. Following Wilson et al. (2012) and Planck Collaboration XXI (2014) we have chosen a hybrid approach, by computing the skewness of the filtered Compton parameter maps outside the 50% sky mask. In particular, we have computed the skewness of the Planck data Compton parameter maps ⟨ y³ ⟩, and of the half-difference maps .

Using the models presented in Appendix A we can show that the unnormalized skewness of the tSZ fluctuation, scales approximately as , whereas the amplitude of the bispectrum scales as with α = 11–12, as shown by Bhattacharya et al. (2012). In the following we do not consider the dependency of the bispectrum and the unnormalized skewness on other cosmological parameters, since such dependencies are expected to be significantly lower than for σ₈ (Bhattacharya et al. 2012).

We derive constraints on σ₈ by comparing the measured unnormalized skewness and bispectrum amplitudes with those obtained from simulations of the tSZ effect. The tSZ contribution was obtained from a hybrid simulation including a hydrodynamic component for z< 0.3, plus extra individual clusters at z> 0.3, and with σ₈ = 0.789. This approach is strongly limited by systematic uncertainties and the details of the theoretical modelling (see Hill & Sherwin 2013). Uncertainties due to foreground contamination are computed using the simulations and are accounted for in the final error bars.

We obtain σ₈ = 0.77 for NILC and σ₈ = 0.78 for MILCA. Combining the two results and considering model and foreground uncertainties, we obtain σ₈ = 0.78 ± 0.02 (68% C.L.). Notice that the reported uncertainties are mainly dominated by foreground contamination. However the model uncertainties only account for the expected dependence of the unnormalized skewness upon σ₈, as shown in Appendix A. As was also the case in Wilson et al. (2012), we have neglected the dependence on other cosmological parameters. We have also not considered any uncertainties coming from the combination of the hydrodynamical and individual cluster simulations. Because of these limitations, our error bars might be underestimated.

7.3.2. Fit of the 1D PDF distribution

We also derive constraints on σ₈ by fitting the 1D PDF obtained in Sect. 6.1. Here, we follow the formalism described in Hill et al. (2014) that evaluates the tSZ 1D PDF theoretically integrating across individual cluster contributions. We specifically use the Tinker et al. (2008) mass function and the Arnaud et al. (2010) pressure profile. The latter is normalized following the Y–M scaling relations described in Planck Collaboration XX (2014), and considering a mass bias parameter of b = 0.2. All cosmological parameters are fixed to the Planck 2015 CMB analysis best-fit values, and we fit only for σ₈.

We note that the Hill et al. (2014) formalism explicitly neglects any effects due to overlapping clusters along the line of sight. For this reason, and given the uncertainties in the modelling of the foreground residuals, the best-fit solution to the observed 1D PDF is fitted in the region that is dominated by non-overlapping clusters and where the noise and foregrounds contributions are minimal. In our case, we use y> 4.5 × 10^-6.

The confidence limits on σ₈ are obtained from a maximum likelihood approach, in which the likelihood has a multivariate Gaussian shape with a covariance matrix which only depends on σ₈. This covariance matrix is evaluated numerically, accounting only for Poisson terms (both for pixel-to-pixel covariance and due to the correlations introduced by the cluster’s y-profile). We obtain σ₈ = 0.77 ± 0.02 (68% C.L.) both for NILC and MILCA maps, including statistical and systematic uncertainties.

7.4. Comparison to other Planck cosmological probes

We have shown in Fig. 17 that the amplitude of the tSZ power spectrum measured in this paper and from the Planck CMB analysis are in good agreement. However, the Planck 2013 results (Planck Collaboration XX 2014) have shown some tension between CMB and tSZ derived constraints on σ₈ for wide range of experiments including Planck. Figure 19 shows the marginalized likelihood distribution for σ₈(Ω_m/ 0.28)^{3 / 8}, as obtained from the combined Planck CMB and BAO analysis (Planck Collaboration XIII 2016) with and without lensing constraints. We also presents the results obtained from the Planck 2015 cluster number count analysis assuming a mass bias of 0.2 (Planck Collaboration XXIV 2016), and for the tSZ power spectrum analysis in this paper assuming a mass bias of 0.2 and 0.4. We observe that assuming a mass bias of 0.2 the two tSZ analyses are in good agreement, but the tension with the CMB measurements remains. Agreement would be better if we adopted a larger mass bias that increases the value of σ₈(Ω_m/ 0.28)^{3 / 8}. Furthermore, we find that including lensing constraints leads to smaller values of σ₈(Ω_m/ 0.28)^{3 / 8}, which as a consequence, are in better agreement with the tSZ results. From this we can conclude that the Large Scale Structure (clusters and lensing) data would have a marginal preference for small values of the mass bias.

8. Summary and conclusions

Because of its wide frequency coverage, from 30 to 857 GHz, the Planck satellite mission is particularly well suited for the measurement of the thermal Sunyaev-Zeldovich effect. Working with the Planck frequency channel maps from 30 to 857 GHz, we have reconstructed the tSZ signal over the full sky using tailored component separation methods.

We tested and validated the Plancky-maps extensively and characterized them in terms of noise and foreground contamination. As expected the noise in the y-map is inhomogeneous and can be characterized by pixel-dependent variance and homogeneous correlated Gaussian noise. Foreground contamination by thermal dust emission is found to be important at large angular scales. Additional foreground contamination is due to radio and IR point sources, for which a mask is provided. In terms of tSZ signal we find good agreement between the flux of blindly detected clusters in the y-map and that measured for clusters in the Planck cluster sample. Furthermore, we find that the sensitivity of the y-map is sufficient to detect faint and diffuse structures such as bridges between merging clusters. Moreover, we have proved, via a stacking analysis, that the very low signal-to-noise regions in the y-map preserve the tSZ signal even for small galaxy groups (tens of galaxies).

After accounting for foreground contribution, mainly thermal dust emission at large angular scales, and clustered CIB and point sources at small angular scales, we have derived from the y-map the tSZ angular power spectrum in the multipole range from 9 <ℓ< 1411. The results are compatible with previous Planck measurements (Planck Collaboration XXI 2014) and extend significantly the multiple range giving for the first time access to the 2-halo term contribution. The cosmological analysis of the tSZ power spectrum allows us to set constraints on cosmological parameters representing the matter content in the Universe, mainly σ₈ and Ω_m. These constraints are consistent with those obtained from cluster number counts (Planck Collaboration XXIV 2016) and in soft tension with those derived from CMB analysis (Planck Collaboration XIII 2016).

The analysis of the non-Gaussian properties of the y-map using the 1D PDF, the unnormalized skewness and the bispectrum of the map have confirmed the tSZ nature of the signal.

The Plancky-maps and additional ancillary data (noise variance maps, foreground masks and ILC weights) are made available to the public for the Planck 2015 release (see Appendix C for details). These y-maps are expected to be useful in a wide range of astrophysical and cosmological analyses with clusters. For any of these analyses, and depending on the scientific goals, the inhomogeneous properties of the noise, and the systematics and foreground contamination should be taken into account in different ways, as described in this paper. Regions masked by the point source mask should never be used. In the case of pixel-based analyses, quality flags can be defined by combining the information from the variance map and the various foreground masks. For power spectrum, cross correlation, and higher order statistic analyses, we note that the y-maps present significant foreground contamination that needs to be taken into account both by masking highly contaminated regions (mainly the Galactic plane region) and by using adequate foreground models to which the ILC weights are applied. Taking these necessary precautions, the Plancky-maps will prove a very useful tool for further studies.

---

Acknowledgments

The Planck Collaboration acknowledges the support of: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, J.A., and RES (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); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planck-collaboration.

References

Appendix A: Modelling the expected tSZ signal

Appendix A.1: tSZ power spectrum

The representation of the y-map in spherical harmonics, Y_ℓm, reads (A.1)Thus, its angular power spectrum is given by (A.2)Note that is a dimensionless quantity here, like y.

As in Planck Collaboration XX (2014), the tSZ power spectrum is modelled using a 2-halo model, to account both for intra-halo and inter-halo correlations: (A.3)Following Komatsu & Seljak (2002) the 1-halo term reads (A.4)where dV_c/ (dzdΩ) is the comoving volume per unit redshift and solid angle and n(M,z)dM dV_c/ (dzdΩ) is the probability of having a galaxy cluster of mass M at a redshift z in the direction dΩ. The quantity is the 2D Fourier transform on the sphere of the 3D radial profile of the Compton y-parameter of individual clusters, (A.5)where x = r/r_s, ℓ_s = D_A(z) /r_s, r_s is the scale radius of the 3D pressure profile, D_A(z) is the angular diameter distance to redshift z, and P_e is the electron pressure profile.

The 2-halo term (Komatsu & Kitayama 1999; Diego & Majumdar 2004; Taburet et al. 2011) is given by (A.6)where P(k,z) is the 3D matter power spectrum at redshift z. The quantity B(M,z) is the time-dependent linear bias factor that relates the matter power spectrum, P(k,z), to the power spectrum of the cluster correlation function. Following Komatsu & Kitayama (1999, see also Mo & White 1996 we adopt B(M,z) = 1 + (ν²(M,z)−1) /δ_c(z), where ν(M,z) = δ_c(M) /D(z)σ(M), σ(M) is the present-day rms mass fluctuation, D(z) is the linear growth factor, and δ_c(z) is the threshold over-density of spherical collapse.

Finally, we use the Tinker et al. (2008) mass function, dn(M,z) / dM, including an observed-to-true mass bias b, as discussed in detail in Planck Collaboration XX (2014), and we model the SZ Compton parameter using the pressure profile of Arnaud et al. (2010).

Appendix A.2: Nth moment of the tSZ field

Assuming a Poisson distribution (1-halo term) of the clusters on the sky, and neglecting clustering between clusters the Nth moment of the tSZ signal (Komatsu & Kitayama 1999; Wilson et al. 2012; Planck Collaboration XX 2014) reads (A.7)where y(θ,M,z) is the integrated Compton parameter along the line of sight for a cluster of mass M at redshift z.

Appendix A.3: Bispectrum

The angular bispectrum is given by (A.8)where the angle-averaged quantity in the full-sky limit can be written as (A.9)and satisfies the conditions m₁ + m₂ + m₃ = 0, ℓ₁ + ℓ₂ + ℓ₃ = even, and |ℓ_i−ℓ_j| ≤ ℓ_k ≤ ℓ_i + ℓ_j, for the Wigner 3j function in brackets. Assuming a Poissonian spatial distribution of the clusters, as above, the bispectrum reads (Bhattacharya et al. 2012)

Appendix B: Bispectrum cosmic variance

Following Lacasa (2014, chapter 2), the bispectrum cosmic variance is composed of a Gaussian term, a bispectrum × bispectrum term, a spectrum × trispectrum term and a connected 6-point term. Due to the lack of model or measurement of the trispectrum and 6-point function, we neglected the last two terms. Note that they are, however, expected to yield a subdominant contribution. Thus we have two contributions (for full sky)

• Gaussian cosmic variance: (B.1)where (B.2)and where C_ℓ is the auto power spectrum of the Compton parameter map, thus containing the noise contribution.

• Bispectrum × bispectrum cosmic variance (B.3)This is the only term which gives off-diagonal contributions to the covariance matrix.

For our purpose, this cosmic variance is multiplied by f_sky and binned to the appropriate binning scheme.

We also consider systematic errors induced by foreground residuals or masking effects. We estimate systematic errors due to component separation uncertainties from the half difference of the NILC and MILCA bispectra. Masking effects are normally corrected for using simulations, which may under- or overestimate leakage from large to small scales. We thus take a conservative ± 25% error on the debiasing ratio, consistent with the fact that the selected configurations have a ratio within ± 25 % of f_skyB(ℓ₁) B(ℓ₂) B(ℓ₃). This error is most likely a conservative overestimate.

Appendix C: Products

In the following we list the y-map related products delivered in the Planck 2015 data release⁴:

• full-sky MILCA and NILC y-maps for the full mission and for the first (F) and second (L) halves of Planck stable pointing period in Compton parameter units (see Sect. 3.2);

• ILC weights per filter and per frequency used for the reconstruction of the MILCA and NILC full mission y-maps, as described in Sect. 3.2;

• variance map accounting for the non-homogeneous coverage and power spectrum of the correlated homogeneous counterpart, , for the MILCA and NILC full mission y-maps in Compton parameter units (see Sect. 4.2);

• point source masks including known radio and IR sources, as described in Sect. 4.4.1;

• Galactic masks used in the analyses presented in Sects. 5 and 6;

All Tables

All Figures