© ESO, 2014

1. Introduction

This paper, one of a set associated with the 2013 data release from the Planck¹ mission (Planck Collaboration I 2014), presents the first results on the integrated Sachs-Wolfe (ISW) effect using Planck data. The ISW effect (Sachs & Wolfe 1967; Rees & Sciama 1968; Martinez-Gonzalez et al. 1990; Hu & Sugiyama 1994) is a secondary anisotropy in the cosmic microwave background (CMB), caused by the interaction of CMB photons with the time-evolving potentials from large-scale structure (LSS). Photons follow a geodesic that is weakly perturbed by the Newtonian gravitational potential, Φ, and experience a fractional shift in their temperature given by $\mathrm{\Theta }\mathrm{=}\frac{\mathrm{\Delta }\mathit{T}}{{\mathit{T}}_{\mathrm{CMB}}}\mathrm{=}\frac{\mathrm{2}}{{\mathit{c}}^{\mathrm{3}}}{\mathrm{\int }}_{{\mathit{\eta }}_{\mathrm{\ast }}}^{{\mathit{\eta }}_{\mathrm{0}}}\mathrm{d}\mathit{\eta }\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial \eta }}\mathit{,}$(1)where the integral is expressed in terms of the conformal time η, defined differentially by dη/ da = 1 / (a²H(a)), with H(a) the Hubble function and a the scale factor. The integration limits here extend from the recombination time (η_∗) to the present time (η₀).

The sensitivity of the ISW effect to gravitational potentials (that can extend over Gpc scales) causes the power of the ISW to be concentrated on the largest scales. The largest scales for the CMB have been mapped out by the Wilkinson Microwave Anisotropy Probe (WMAP) to the statistical limit of cosmic variance. Some systematics (like foreground removal) can have an impact on the reconstruction of the CMB especially at the largest scales where our Galaxy can introduce significant residuals on the reconstructed CMB map. The superior sensitivity of Planck together with its better angular resolution and wider frequency coverage allows for a better understanding (and hence removal) of Galactic and extagalactic foregrounds, which reduces the possible negative impact of these residuals. Planck allows us to improve on previous measurements through a better systematic control, an improved removal of foregrounds (which permits us to explore the achromatic nature of the ISW signal on a wider frequency range), and a better understanding of systematics affecting tracer catalogues.

For cosmological models where Ω_m = 1, gravitational potentials remain constant during linear structure formation, and the ISW signal is negligible (to first order, although second-order nonlinear ISW is always expected around smaller over- and underdense regions). In the presence of dark energy, decaying potentials due to the accelerated expansion rate result in a net ISW effect that is positive when the CMB photons cross overdense regions and negative when the CMB photons cross underdense regions. Therefore, the ISW effect is an indicator of either nonzero curvature (Kamionkowski & Spergel 1994)², any form of dark energy, such as a cosmological constant Λ (Crittenden & Turok 1996), modified gravity (Hu 2002), or a combination of these possibilities. By measuring the rate at which gravitational potentials in the LSS decay (up to a redshift of about 2), the ISW effect can be used as an independent probe of cosmology and provides complementary and independent evidence for dark energy.

Detection of the ISW effect was first made possible with all-sky CMB maps from WMAP. Based on these data, many works can be found in the literature where the authors aim at making, and subsequectly improving, the measurement of the ISW effect through correlations with tracer catalogues: 2MASS (an infrared catalogue out to low redshifts around 0.1, Afshordi et al. 2004; Rassat et al. 2007; Francis & Peacock 2010b; Dupé et al. 2011), HEAO (an X-ray survey at low redshift, with the first positive claim of detection, Boughn & Crittenden 2004), Sloan Digital Sky Survey (SDSS, an optical survey at intermediate redshifts, Fosalba et al. 2003; Scranton et al. 2003; Fosalba & Gaztañaga 2004; Padmanabhan et al. 2005; Cabré et al. 2006; Giannantonio et al. 2006; Granett et al. 2009; Xia 2009; Bielby et al. 2010; López-Corredoira et al. 2010; Sawangwit et al. 2010), the NRAO VLA Sky Survey (NVSS, a radio catalogue with high-redshift sources, Boughn & Crittenden 2005; Vielva et al. 2006; Pietrobon et al. 2006a; McEwen et al. 2007; Raccanelli et al. 2008; Hernández-Monteagudo 2010; Massardi et al. 2010; Schiavon et al. 2012), and combined measurements with multiple tracers (Nolta et al. 2004; Ho et al. 2008; Corasaniti et al. 2005; Gaztañaga et al. 2006; Giannantonio et al. 2008, 2012). The significance of the ISW detections that can be found in the literature range between 0.9σ and 4.7σ. There are a number of peculiarities related to some of the detection claims, as noted by Hernández-Monteagudo (2010) and López-Corredoira et al. (2010). They both found lower significance levels than some previous studies and pointed out the absence of the signal at low multipoles where the ISW effect should be most prominent and the presence of point source emission on small scales for radio surveys.

The main result that is obtained from an ISW detection is a constraint on the cosmological constant (or dark energy), Ω_Λ. The general consensus from the variety of ISW analyses is for a value of Ω_Λ ≃ 0.75 with an error of about 20%, which provides independent evidence for the existence of dark energy (Fosalba et al. 2003; Fosalba & Gaztañaga 2004; Nolta et al. 2004; Corasaniti et al. 2005; Padmanabhan et al. 2005; Cabré et al. 2006; Giannantonio et al. 2006; Pietrobon et al. 2006b; Rassat et al. 2007; Vielva et al. 2006; McEwen et al. 2007; Ho et al. 2008; Schiavon et al. 2012). All tests on spatial flatness found an upper limit for Ω_K of a few percent (Nolta et al. 2004; Gaztañaga et al. 2006; Ho et al. 2008; Li & Xia 2010). Using a prior on spatial flatness, the dark energy equation of state parameter, w, was found to be close to −1 (Giannantonio et al. 2006; Vielva et al. 2006; Ho et al. 2008), and a strong time evolution has been excluded (Giannantonio et al. 2008; Li & Xia 2010).

The ISW effect is achromatic, that is, it conserves the Planck spectrum of the CMB. Nevertheless, it can be separated from other CMB fluctuations through cross-correlations with catalogues that trace the LSS gravitational potentials (see for instance Crittenden & Turok 1996). This cross-correlation can be studied in different ways: through angular cross-correlations in real space between the CMB and the catalogues tat trace the LSS; through the corresponding angular cross-power spectrum of the Fourier-transformed maps; or through the covariance of wavelet-filtered maps as a function of wavelet scale. The studies using WMAP data mentioned above follow this survey cross-correlation techique.

An alternative approach, similar to the angular cross-correlation in real space, consists of stacking CMB fields centred on known supersclusters or supervoids (Granett et al. 2008a,b; Pápai & Szapudi 2010). The advantage of this technique is that it allows for a detailed study of the profile of the CMB fluctuations caused by this secondary anisotropy.

A novel and powerful approach takes advantage of the fact that the same potentials that cause CMB photons to gain or lose energy along their path (ISW) create lensing distortions that can be directly measured from the CMB map (e.g., Hu & Okamoto 2002). The interplay between weak gravitational lensing and the ISW effect causes a non-Gaussian contribution to the CMB, which can be measured through the lensing-induced bispectrum between small and large angular scales. The measurement of the lensing potential requires a large number of modes that could not be measured before the arrival of Planck data.

This paper presents new measurements of the ISW effect carried out with Planck. Even although our detections are not in every case as strong as some previously claimed significance levels, we believe that our results are an improvement over earlier studies. This is because we can use the additional power enabled by the frequency coverage and sensitivity of Planck. To establish this we carried out a comprehensive study of all the main approaches that have previously been taken to estimate the ISW signal. We also present new results in relation to the non-Gaussian structure induced by the ISW effect.

The paper is organized as follows: in Sect. 2 we describe the data used in this work (both for the CMB and large-scale structure). The first-ever results of estimating the lensing-induced bispectrum are presented in Sect. 3. Cross-correlations with external surveys are investigated in Sect. 4, and in Sect. 5 we present the results for the stacking analysis on the temperature maps, and aperture photometry on superclusters and supervoids. The recovery of the ISW all-sky map is described in Sect. 6. Finally, we discuss our main results and their cosmological implications in Sect. 7.

2. Data description

In this section we describe the different data sets. This includes Planck data (the CMB temperature and lensing potential maps, see Planck Collaboration I 2014; Planck Collaboration II 2014; Planck Collaboration VI 2014; Planck Collaboration XII 2014; Planck Collaboration XVII 2014, for details), and external data sets (large-scale structure tracers) used in determining the ISW: the radio NVSS catalogue, optical luminous galaxies (CMASS/LOWZ) and the main galaxy sample from the Sloan Digital Sky Survey (SDSS), and several superstructure catalogues.

2.1. Planck data

Planck data and products are described in the following sections, in particular the foreground-cleaned CMB maps produced by the Planck component separation pipelines, and related products, such as dedicated component-separated frequency maps (Planck Collaboration XII 2014), and the Planck lensing map (Planck Collaboration XVII 2014).

2.1.1. CMB maps

We used the Planck foreground-cleaned CMB maps provided by the data processing centres (as described in the Planck component separation paper, Planck Collaboration XII 2014). In particular, to test their reliability, some of the results are presented for different cleaned CMB maps, which were constructed using four different component separation techniques: Commander-Ruler (C-R, which uses physical parametrization), NILC (an internal linear combination technique), SEVEM (a template-fitting method), and SMICA (which uses spectral matching). Since the contribution of the ISW signal is only significant on large scales, we used low-resolution maps, with HEALPix Górski et al. (2005) parameter N_side = 64 and a pixel size of about 55 arcmin for most of the analyses. One exception is the study of the correlation between the ISW and lensing signals, which requires the use of full-resolution maps (N_side = 2048, pixel size of 1.7 arcmin). The maps were degraded directly from the original full resolution to the corresponding N_side.

In addition, foreground-cleaned maps per frequency (from 44 to 353 GHz) at resolution N_side = 512 were used for the stacking analysis presented in Sect. 5. These maps were constructed by subtracting a linear combination of internal templates using SEVEM (for a detailed description of the method see the SEVEM appendix of Planck Collaboration XII 2014). As an example, the SEVEM CMB map is shown in Fig. 1 (left panel).

Finally, to minimize the foreground contamination in the maps, we used the official mask described in Planck Collaboration XII (2014), which excludes regions with stronger Galactic and point-source contamination (the U73 mask). This mask is given at the full Planck resolution and is downgraded to the required levels. The downgrading procedure consists of the following steps: the mask (originally a map with zero and one values) is convolved with a Gaussian beam of FWHM three times the characteristic pixel size of the final N_side resolution; this convolved map is then degraded to the required N_side, and, subsequently, a threshold of 0.75 is imposed (i.e., pixels with a value higher than this threshold are set to one, whereas the rest are set to zero).

2.1.2. Lensing potential map

Weak gravitational lensing distorts the CMB temperature anisotropy pattern. This effect is sensitive to the projected matter distribution in the large-scale structure at high redshifts, where structure growth is linear and the statistics are close to Gaussian. Weak lensing causes correlations between different multipoles, which are proportional to the lensing deflection field. These correlations can be exploited to reconstruct the density field and to measure its statistical properties (Hu & Okamoto 2002; Okamoto & Hu 2003). The lensing effect in the CMB can be estimated by this homogeneity breaking, and in this way, individual modes of the lensing potential at multipoles ℓ < 100 can be reconstructed with a significance of about 0.5σ, which shows that a statistical treatment is necessary. Nevertheless, the overall effect of the lensing is measured to better than 25σ (Planck Collaboration XVII 2014). The additional lensing effect in the temperature power spectrum is detectable with a significance of about 10σ (Planck Collaboration XV 2014).

With Planck data, we aim at detecting a correlation between the ISW effect and the lensing potential, where the latter is a tracer of the large-scale structure at high redshift. This correlation is restricted to 9σ, even in the ideal case, limited by cosmic variance and the weakness of the ISW effect in comparison to the primary CMB (Lewis et al. 2011). The data products used in this study are the Planck lensing potential reconstruction and specific lensing maps obtained from the component separation pipelines. The lensing potential is available as part of the first Planck data release. Its detailed development is described in the Planck lensing paper (Planck Collaboration XVII 2014). In Fig. 1 we reproduce (right panel) an optimally filtered version of the Planck lensing map, suitable for the ISW-lensing cross-correlation.

In addition to a direct correlation between the CMB sky and the reconstructed lensing map, we measured the bispectrum generated by weak lensing by applying a range of estimators: the KSW-bispectrum estimator, bispectra binned in multipole intervals and a modal decomposition of the bispectrum. This meassurement is made possible for the first time thanks to the Planck data. In addition, we used information from the lensing field as a tracer for an ISW map reconstruction at high redshift (see Sect. 6).

2.2. External data sets

As described in the introduction, the achromatic nature of the ISW effect requires a tracer of the gravitational potentials from the large-scale structure, so that by cross-correlating the CMB temperature map with that tracer distribution the fluctuations generated by the ISW effect are singled out. The prerequisites for a tracer catalogue to be used in ISW studies are a large survey volume, well-understood biasing properties, and low or at least well-modelled systematics. The radio NVSS catalogue and the optical luminous galaxies (SDSS-CMASS/LOWZ) and main photometric galaxy sample (SDSS-MphG) catalogues possess these qualities. Table 7 summarizes some basic properties of these catalogues. In addition, the redshift distributions of these catalogues are shown in Fig. 2. Note that NVSS presents the widest redshift coverage. The SDSS-CMASS/LOWZ sample peaks at about z ≈ 0.5, whereas the SDSS-MphG sample peaks at about z ≈ 0.3.

Figure 3 shows the all-sky density projection for these maps, where the grey area indicates regions not observed by these surveys (or that were discarded because they were contaminated or had a low galaxy number density, see next subsections for details). In Fig. 4 we give the angular power spectra (blue points) of the surveys (corrected for with a procedure similar to MASTER, e.g., Hivon et al. 2002), and the theoretical spectra (black lines) and their 1σ error bars (grey areas), estimated from the MASTER approach as well.

In addition to the cross-correlation between CMB and LSS tracers (Sect. 4), we present results from a different methodology in Sect. 5, where we used catalogues of superstructures to study the ISW through stacking of the CMB fluctuations on the positions of these superstructures. The relevant catalogues are described in Sect. 2.2.4.

2.2.1. NVSS radio catalogues

Luminous active galactic nuclei (hereafter AGN) are known to be powerful radio sources that are visible out to high redshifts. These sources are hence able to probe the cosmic density field during the entire redshift range from matter domination to accelerated expansion due to dark energy. If AGN are fair tracers of the underlying density field, these sources should likewise probe the spatial distribution of the large-scale potential wells that decay at late times after the accelerated expansion sets in and generates the ISW effect.

We focused on a single radio survey that has the level of sensitivity and sky coverage required for ISW studies, namely the NRAO VLA Sky Survey (hereafter NVSS, Condon et al. 1998). This survey was conducted using the Very Large Array (VLA) at 1.4 GHz, and covers up to an equatorial latitude of b_E = − 40deg, with an average noise level of 0.45 mJy beam^-1. It results in roughly 1.4 × 10⁶ sources above a flux threshold of 2.5 mJy. Figure 3 displays the number density map computed from the NVSS survey (top panel). The AGN population is known to be dominant in radio catalogues at 1.4 GHz in the high flux density regime. Condon et al. (1998) showed that at this frequency, star-forming galaxies (SFGs) contribute about 30% of the total number of weighted source counts above 1 mJy, but their presence decreases rapidly as higher flux thresholds are adopted. The NVSS SFGs are nearby sources (z< 0.01), and hence may distort the ability of our radio template to probe the intermediate and high-redshift density field.

We next address the systematic effects in the NVSS survey. Two different antenna configurations were used while conducting the NVSS survey: the D-configuration (for b_E ∈ [ − 10deg, 78deg ]), and the DnC-configuration for large zenith angles (b_E< − 10deg,b_E> 78deg). This change in the antenna configuration is known to introduce changes in the source number density above 2.5 mJy, as first pointed out by Blake & Wall (2002). The NVSS map at 2.5 mJy was corrected for this declination systematic using the following procedure: the sky was divided into equatorial strips and the mean number of sources in each strip was re-normalized to the full sky mean (see e.g., Vielva et al. 2006). With this procedure the average number of sources in the NVSS map is the same as before the correction, and hence the shot noise level does not change. The number of strips into which the map was divided is 70, but the results are independent of this choice.

Regarding the galaxy bias, we adopted the Gaussian bias-evolution model of Xia et al. (2011). If n(M,z) is the halo mass function and b(M,z) is the bias of halos with comoving mass M, then the bias of the survey is given by a mass-weighted integral, $\mathit{b}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\int }}_{{\mathit{M}}_{\min }}^{\mathrm{\infty }}\mathrm{d}\mathit{M}\mathit{b}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathit{M}\mathit{n}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}}{{\mathrm{\int }}_{{\mathit{M}}_{\min }}^{\mathrm{\infty }}\mathrm{d}\mathit{M}\mathit{M}\mathit{n}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}}\mathrm{·}$(2)This model depends on the minimum mass M_min of halos present in the survey. The upper limit in the mass is taken to be infinity because the effect of the high-mass end on the bias is negligible. Marcos-Caballero et al. (2013) proposed a theoretical model for the NVSS angular power spectrum, which also takes into account the information of the redshift distribution given by CENSORS data (Brookes et al. 2008). The redshift distribution is parametrized by $\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{z}}\mathrm{=}{\mathit{n}}_{\mathrm{0}}{\left(\frac{\mathit{z}}{{\mathit{z}}_{\mathrm{0}}}\right)}^{\mathit{\alpha }}{\mathrm{e}}^{\mathrm{-}\mathit{\alpha z}\mathit{/}{\mathit{z}}_{\mathrm{0}}}\mathit{,}$(3)where z₀ = 0.33 and α = 0.37. The parameter n₀ is a constant to obtain a distribution normalized to unity. This function is represented in Fig. 2. The bias follows the prescription of Eq. (2), with M_min equal to 10^12.67M_⊙, where the Sheth-Tormen (Sheth & Tormen 1999) mass function is adopted. Hereafter this model is our fiducial model for NVSS.

2.2.2. SDSS luminous galaxies

For this analysis we used the photometric luminous galaxy (LG) catalogue from the Baryonic Oscillation Spectroscopic Survey (BOSS) of the SDSS III. The data used consist of two sub-samples: CMASS and LOWZ. The two samples were combined to form a unique LG map (see Fig. 3, second panel). Hereafter, these samples are referred to as SDSS-CMASS, SDSS-LOWZ, and SDSS-CMASS/LOWZ for the combination.

SDSS-CMASS

We used the BOSS targets to obtain an approximately constant stellar mass; this is known as the photometric CMASS sample. This sample is mostly contained in the redshift range z = 0.4 − 0.7, with a galaxy number density close to 110 deg^-2, and was selected after applying the colour cuts explained in Ross et al. (2011).

While this colour selection yields a catalogue of about 1 600 000 galaxies, more cuts were needed to be applied to account for dust extinction (based on the maps by Schlegel et al. 1998 with the criterion E(B − V) < 0.08), for seeing in the r band (required to be <2.0′′), and for the presence of bright stars, similar to Ho et al. (2012). Finally, we neglected all pixels with a mask value inferred from the footprint below 0.9 on a HEALPix map of resolution N_side = 64. This procedure left about one million sources 10 500 deg². Photometric redshifts of this sample were calibrated using a selection of about 100 000 BOSS spectra as a training sample for the photometric catalogue. These LGs are among the most luminous galaxies in the Universe and therefore allow for a good sampling of the largest scales. Given the large number of the sources included in the sample, shot-noise does not dominate clustering errors. According to Ross et al. (2011), about 3.7% of these objects are either stars or quasars, which necessitates additional corrections, as explained at the end of this section.

SDSS-LOWZ

The photometric LOWZ sample is one of the two galaxy samples targeted by the BOSS of Sloan III. It selectd luminous, highly biased, mostly red galaxies, placed at an average redshift of ⟨ z ⟩ ~ 0.3 and below the redshifts of the CMASS sample (z< 0.4). Our selection criteria in terms of the Sloan five-model magnitudes ugriz followed those given in Sect. 2 of Parejko et al. (2013). With a total number of sources close to 600 000, this photometric sample contains a higher number density of galaxies in the southern part of the footprint than in the northern one (by more than 3%), which seems to be at odds with ΛCDM predictions. However, most of this effect vanishes after subtracting the dipole in the effective area under analysis in such a way that the low ℓ range of the auto power spectrum is consistent with a ΛCDM model and a constant bias b ≃ 2 (Hernández-Monteagudo et al. 2013).

The SDSS-CMASS and SDSS-LOWZ samples were also corrected for scaling introduced by possible systematics such as stars, mask value, seeing, sky emission, airmass, and dust extinction. Following exactly the same procedure as in Hernández-Monteagudo et al. (2013), we found that both LG samples were contaminated by stars in the sense that the galaxy number density decreased in areas with higher star density, since the latter tend to “blind” galaxy detection algorithms.

2.2.3. Main photometric SDSS galaxy sample

We used a sample of photometrically selected galaxies from the SDSS-DR8 catalogue, which covers a total sky area of 14 555 deg² (Aihara et al. 2011). The total number of objects labelled as galaxies in this data release is 208 million. From this catalogue, and following Cabré et al. (2006), we defined a subsample by selecting only objects within the range 18 <r< 21, where this r-band model magnitude corrected for extinction. Following Giannantonio et al. (2008), we also restricted our subsample to objects with redshifts 0.1 <z< 0.9 and with measured redshift errors such that σ_z< 0.5z. We relied on the photometric redshift estimates of the SDSS photo-z primary galaxy table, which were obtained through a kd-tree nearest-neighbour technique by fitting the spectroscopic objects observed with similar colour and inclination angle. The total number of galaxies in our final sample is about 42 million, with redshifts distributed around a median value of about 0.35, as shown in Fig. 2. To avoid possible errors introduced by singularities in the photometric redshift estimates, instead of using the real observed redshift distribution in our analysis, we resorted to the analytical function $\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{z}}\mathrm{=}\frac{\mathit{\beta }}{\mathrm{\Gamma }\left(\frac{\mathit{m}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}{\mathit{\beta }}\right)}\frac{{\mathit{z}}^{\mathit{m}}}{{\mathit{z}}_{\mathrm{0}}^{\mathit{m}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{z}\mathit{/}{\mathit{z}}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\beta }}}\mathit{,}$(4)which was fitted to the data with parameters m = 1.5, β = 2.3 and z₀ = 0.34, which are identical to those found by Giannantonio et al. (2012). For the galactic bias we used the value b = 1.2, which was found by Giannantonio et al. (2012) by fitting the ΛCDM prediction to the observed auto-correlation function of the galaxies, and we adopted the mask proposed by them.

2.2.4. SDSS superstructures

Granett et al. (2008b) produced a sample³ of 50 superclusters and 50 supervoids identified from the luminous red galaxies (LRGs) in the SDSS (sixth data release, DR6, Adelman-McCarthy et al. 2008), which covers an area of 7500 deg² on the sky. They used publicly available algorithms, based on the Voronoi tessellation, to find 2836 superclusters (using VOBOZ, VOronoi BOund Zones, Neyrinck et al. 2005) and 631 supervoids (using ZOBOV, ZOnes Bordering On Voidness, Neyrinck 2008) above a 2σ significance level (defined as the probability of obtaining, in a uniform Poisson point sample, the same density contrasts as those of clusters and voids). The 50 superclusters and 50 supervoids they published in their catalogue correspond to density contrasts of about 3σ and 3.3σ, respectively. They span a redshift range of 0.4 <z< 0.75, with a median of around 0.5, and inhabit a volume of around 5 h^-3 Gpc³. These superstructures can potentially produce measurable ISW signals, as suggested in Granett et al. (2008a,b). For each structure, the catalogue provides the position on the sky of the centre, the mean and maximum angular distance between the galaxies in the structure and its centre, the physical volume and three different measures of the density contrast (calculated from all its Voronoi cells, from only its over- or underdense cells, and from only its most over- or underdense cells). We concentrated here on using the supervoid catalogue by Granett et al. (2008b) because it can be compared with two other, more recent catalogues of voids.

The second catalogue of cosmic voids that we considered here was published by Pan et al. (2012)⁴. It has been built from the seventh data release (DR7) of the SDSS. Using the VoidFinder algorithm (Hoyle & Vogeley 2002), the authors identified 1055 voids with redshifts lower than z = 0.1. Each void is listed with its position on the sky, its physical radius (defined as the radius of the largest sphere enclosing the void), an effective radius defined as the radius of a sphere of the same volume, its physical distance to us, its volume, and its mean density contrast. The filling factor of the voids in the sample volume is 62%. The largest void is just over 47 Mpc in effective radius, while the median effective radius of the void sample is about 25 Mpc. Some of the voids are both very close to us and relatively large (larger than 30 Mpc in radius), which results in large angular sizes of up to 15deg.

The third void catalogue that we used has been released by Sutter et al. (2012) and was also made publicly available⁵. Note that it is being updated regularly, and the results reported here are based on the version of 21 February 2013 of the catalogue. Using their own modified version of ZOBOV, these authors identified 1495 voids distributed across the 0–0.44 redshift range. They subdivided their catalogue into six subsamples: dim1, dim2, bright1, and bright2, constructed from the main SDSS, and lrgdim and lrgbright built from the SDSS LRG sample. For each void, the information provided includes the position of the centre, the redshift, the volume, the effective radius, and the density contrast.

3. ISW-lensing bispectrum

There is an interesting interplay between gravitational lensing of the CMB and the ISW effect, which manifests itself as a non-Gaussian feature. CMB-lensing can be described by a convolution of the CMB-temperature map T with the weak-lensing potential φ (e.g., Lewis & Challinor 2006), $\mathit{T}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{\to }\mathit{T}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{-}\mathrm{\int }\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{\ell }}^{\mathrm{\prime }}}{\mathrm{2}\mathit{\pi }}{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{(}\mathit{\ell }\mathrm{-}{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{)}\mathit{\phi }\mathrm{(}\mathit{\ell }\mathrm{-}{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{)}\mathit{T}\mathrm{\left(}{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{\right)}\mathit{.}$(5)The CMB lensing can be measured by a direct estimate of the CMB bispectrum, because the bispectrum acquires first-order terms proportional to the product of two power spectra $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}\end{array}{\mathit{C}}_{\mathit{\ell }}^{\mathit{T\phi }}$, where $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}\end{array}$ is the lensed temperature power spectrum and ${\mathit{C}}_{\mathit{\ell }}^{\mathit{T\phi }}$ is the temperature-potential cross-power spectrum. The potential field φ and the temperature field T are correlated, because φ, which deflects the CMB photons, also gives rise to the ISW effect in T (Hu & Okamoto 2002; Seljak & Zaldarriaga 1999; Verde & Spergel 2002; Giovi et al. 2003). This secondary bispectrum contains new information about the cosmological redshift, because it is generated mainly at redshifts higher than unity, and biases measurements of the primordial bispectrum. The term ${\mathit{C}}_{\mathit{\ell }}^{\mathit{T\phi }}$ correlates the CMB temperature on small scales with the lensing potential on large scales, and causes the bispectrum to assume large amplitudes in the squeezed-triangles configuration (see e.g., Goldberg & Spergel 1999; Seljak & Zaldarriaga 1999; Hu 2000; Giovi et al. 2003; Okamoto & Hu 2003; Giovi & Baccigalupi 2005; Lewis & Challinor 2006; Serra & Cooray 2008; Mangilli & Verde 2009; Hanson et al. 2009, 2010; Smith & Zaldarriaga 2011; Lewis et al. 2011).

Owing to the rotational invariance of the sky, the CMB angular bispectrum ⟨ a_ℓ1m1a_ℓ2m2a_ℓ3m3 ⟩ can be factorized as follows: $\mathrm{⟨}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\mathrm{⟩}\mathrm{=}{\mathrm{𝒢}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathit{,}$(6)where ${\mathrm{𝒢}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}\mathrm{\equiv }{}^{\mathrm{\int }}\mathrm{d}^{\mathrm{\Omega }}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{1}}}^{{\mathit{m}}_{\mathrm{1}}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{m}}_{\mathrm{2}}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}{\mathit{Y}}_{{\mathit{l}}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{3}}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ is the Gaunt-integral and b_ℓ1ℓ2ℓ3 is the so-called reduced bispectrum. When the bispectral signal on the sky is generated by the ISW-lensing effect, ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{=}{\mathit{A}}^{\mathit{T\phi }}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{lens}\mathrm{-}\mathrm{ISW}}$, where A^Tφ parametrizes the amplitude of the effect and $\begin{array}{ccc}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{lens}\mathrm{-}\mathrm{ISW}}& \mathrm{=}& \{\frac{{\mathit{\ell }}_{\mathrm{1}}\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{-}{\mathit{\ell }}_{\mathrm{2}}\mathrm{(}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{+}{\mathit{\ell }}_{\mathrm{3}}\mathrm{(}{\mathit{\ell }}_{\mathrm{3}}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{2}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathit{T\phi }}\stackrel{}{{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}^{\mathit{TT}}}_{\mathrm{˜}}}\\ & & \end{array}$(7)A more general expression for intensity and polarization can be found in Lewis et al. (2011). Estimating the bispectrum then yields a measurement of A^Tφ.

We can also define an alternative rotationally invariant reduced bispectrum B_ℓ1ℓ2ℓ3 as ${\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{=}{\mathit{h}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{2}}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}$, where $\begin{array}{ccc}{\mathit{h}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{2}}& \mathrm{\equiv }& \sum _{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}{\mathrm{(}{{\mathrm{𝒢}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}}^{\mathrm{)}}}^{\mathrm{2}}\\ & \mathrm{=}& \frac{\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{3}}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}{\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)}^{\mathrm{2}}\mathit{.}\end{array}$(8)The interest in B_ℓ1ℓ2ℓ3 is that it can be directly estimated from the observed map using the expression ${\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{obs}}\mathrm{=}\mathrm{\int }\mathrm{d\Omega }{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{1}}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{2}}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{3}}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\mathit{,}$(9)where the filtered maps T_ℓ(Ω) are defined as ${\mathit{T}}_{\mathit{\ell }}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\mathrm{\equiv }\sum _{\mathit{m}}{\mathit{a}}_{\mathit{\ell m}}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\mathit{.}$(10)By basically combining the single-ℓ estimates B^obs/B^{lens − ISW} for A^Tφ using inverse variance weighting, the ISW-lensing bispectrum estimator can be written as (see Planck Collaboration XXIV 2014, for more details) ${\mathit{Â}}^{\mathit{T\phi }}\mathrm{=}\frac{\mathrm{⟨}{\mathit{B}}^{\mathrm{lens}\mathrm{-}\mathrm{ISW}}\mathit{,}\mathrm{(}{\mathit{B}}^{\mathrm{obs}}\mathrm{-}{\mathit{B}}^{\mathrm{lin}}\mathrm{)}\mathrm{⟩}}{\mathrm{⟨}{\mathit{B}}^{\mathrm{lens}\mathrm{-}\mathrm{ISW}}\mathit{,}{\mathit{B}}^{\mathrm{lens}\mathrm{-}\mathrm{ISW}}\mathrm{⟩}}\mathit{,}$(11)where the inner product is defined by $\mathrm{⟨}{\mathit{B}}^{\mathit{i}}\mathit{,}{\mathit{B}}^{\mathit{j}}\mathrm{⟩}\mathrm{\equiv }\sum _{{\mathit{\ell }}_{\mathrm{1}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{2}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{3}}}\frac{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathit{i}}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathit{j}}}{{\mathit{V}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}}\mathrm{·}$(12)Here, B^lin is a linear correction that has zero average but reduces the variance of the estimator in the presence of anisotropic noise and a mask. Furthermore, ${\mathit{V}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{=}{\mathit{g}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}{\mathit{h}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{2}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}$, with g being a simple permutation factor (g = 6 when all ℓ are equal, g = 2 when two ℓ are equal and g = 1 otherwise). As in all expressions in this section, we implicitly took the beam and noise of the experiment into accout, which means that C_ℓ should actually be ${\mathit{b}}_{\mathit{\ell }}^{\mathrm{2}}{\mathit{C}}_{\mathit{\ell }}\mathrm{+}{\mathit{N}}_{\mathit{\ell }}$, with b_ℓ the beam transfer function and N_ℓ the noise power spectrum.

In Eq. (11) we have also used the fact that, as discussed in detail in Planck Collaboration XXIV (2014), the full inverse covariance weighting can be replaced by a diagonal covariance term, (C^-1a)_ℓm → a_ℓm/C_ℓ, without loss of optimality, if the masked regions of the map are filled in with a simple diffusive inpainting scheme.

The normalization of the lensing-ISW estimator in the denominator of Eq. (11) can be replaced by (see e.g. Lewis et al. 2011) $\mathit{F}\mathrm{=}\sum _{\mathit{\ell }}({\mathit{F}}_{\mathit{\ell }}^{-1}\mathrm{+}\frac{\mathrm{1}\mathrm{+}{\mathit{r}}_{\mathit{\ell }}^{-2}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}{)}^{-1}\mathit{,}$(13)where ${\mathit{r}}_{\mathit{\ell }}\mathrm{\equiv }{\mathit{C}}_{\mathit{\ell }}^{\mathit{T\phi }}\mathit{/}\sqrt{\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}\end{array}{\mathit{C}}_{\mathit{\ell }}^{\mathit{\phi \phi }}}$ parameterizes the deviation from the Cauchy-Schwarz relation and F_ℓ is given in terms of the ISW-lensing bispectrum (see for example Lewis et al. 2011). The first term in Eq. (13) corresponds to the Fisher errors assuming Gaussian a_ℓm. However, contrary to the null hypothesis that is assumed, for example, in the primordial bispectra (Gaussianity), there is an actual non-Gaussian signal already present in the ISW-lensing bispectrum. This guarantees a larger variance for the estimators than are included in the additional terms present in the previous equations.

An important issue is the impact of the ISW-lensing bispectrum on estimates of the primordial non-Gaussianity. Assuming weak levels of non-Gaussianity and considering both the primordial bispectrum ${\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{prim}}$ and the ISW-lensing bispectrum ${\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{lens}\mathrm{-}\mathrm{ISW}}$, one can compute the expected bias Δ induced in the primordial bispectrum using the formula: ${\mathrm{\Delta }}^{\mathrm{prim}}\mathrm{=}\frac{\mathrm{⟨}{\mathit{B}}^{\mathrm{lens}\mathrm{-}\mathrm{ISW}}\mathit{,}{\mathit{B}}^{\mathrm{prim}}\mathrm{⟩}}{\mathrm{⟨}{\mathit{B}}^{\mathrm{prim}}\mathit{,}{\mathit{B}}^{\mathrm{prim}}\mathrm{⟩}}$(14)with the inner product defined in Eq. (12). Predictions of this bias on the primordial f_NL for the Planck resolution can be seen for example in Hanson et al. (2009), Mangilli & Verde (2009), Smith & Zaldarriaga (2011), and Lewis et al. (2011). The most important bias is introduced to the local shape and, considering ℓ_max ~ 2000, is expected to be Δ^local ~ 7 (Planck Collaboration XXIV 2014).

3.1. ISW-lensing estimators

There are several implementations of the optimal estimator given in Eq. (11). For a detailed description in the context of Planck see Planck Collaboration XXIV (2014) and Planck Collaboration XVII (2014). We applied four of these implementations to Planck data to constrain the ISW-lensing bispectrum. Three of them represent a direct bispectrum estimation: the KSW estimator (Komatsu et al. 2005; Creminelli et al. 2006), the binned bispectrum (Bucher et al. 2010), and the modal decomposition (Fergusson et al. 2010).The remaining approach is based on a previous estimation of the gravitational lensing potential field Lewis et al. (2011). These estimators differ in the implementation and approximations that are used to compute the expression given in Eq. (11), the direct computation of which is beyond current computing facilities. They are reviewed in the next subsections.

3.1.1. Lensing potential reconstruction

The estimator given in Eq. (6) can be written in terms of the lensing potential amplitude reconstruction as $\mathit{Â}\mathit{T\phi }\mathrm{\equiv }\mathit{Ŝ}\mathrm{=}\frac{\mathrm{1}}{\mathit{N}}\sum _{\mathit{\ell m}}{\mathit{C}}_{\mathit{\ell }}^{\mathit{T\phi }}\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathit{\ell m}}\end{array}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}\end{array}}\frac{\mathit{\phi ̂}\begin{array}{c}\mathrm{\ast }\\ \mathit{\ell m}\end{array}}{{\mathit{N}}_{\mathit{\ell }}^{\mathit{\phi \phi }}}\mathit{,}$(15)where $\mathit{\phi ̂}\begin{array}{c}\mathrm{\ast }\\ \mathit{\ell m}\end{array}$ can be estimated using a quadratic estimator (Okamoto & Hu 2003) and ${\mathit{N}}_{\mathit{\ell }}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ is given in terms of the ISW-lensing bispectrum (Lewis et al. 2011). Therefore, this estimator quantifies the amount of cross-correlation between the temperature map and the reconstruction of the lensing signal, and most of the correlation is found at multipoles below 100.

3.1.2. KSW-estimator

The KSW bispectrum estimator (Komatsu et al. 2005) for the ISW-lensing signal can be written as $\mathit{Â}\mathit{T\phi }\mathrm{=}\mathrm{\left(}{\mathit{F}}^{-1}\mathrm{\right)}\mathit{Ŝ}\mathit{,}$(16)where Ŝ can be computed from data as $\begin{array}{ccc}\mathit{Ŝ}& \mathrm{\equiv }& \frac{\mathrm{1}}{\mathrm{6}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\sum _{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\sum _{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}{\mathrm{𝒢}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{lens}\mathrm{-}\mathrm{ISW}}\\ & & \mathrm{\times }[\mathrm{\left(}{\mathit{C}}^{-1}\mathit{a}{\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\mathrm{\right(}{\mathit{C}}^{-1}\mathit{a}{\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{(}{\mathit{C}}^{-1}\mathit{a}{\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\\ & & \mathrm{-}\mathrm{3}\mathrm{\left(}{\mathit{C}}^{-1}{\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{\right(}{\mathit{C}}^{-1}\mathit{a}{\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}]\mathit{,}\end{array}$(17)and (F^-1) is the inverse of the ISW-lensing Fisher matrix F of Eq. (13). Details on the implementation of the KSW estimator for the ISW-lensing signal can be found in Mangilli et al. (2013). In particular, Eq. (16) takes the form $\mathit{Â}\mathit{T\phi }\mathrm{=}\mathrm{\left(}{\mathit{F}}^{-1}\mathrm{\right)}\mathrm{(}\mathit{Ŝ}\mathrm{cubic}\mathrm{+}\mathit{Ŝ}{\mathrm{linear}}^{\mathrm{)}}\mathit{,}$(18)where Ŝ_cubic is the term that extracts the amplitude information from the data contained in the bispectrum, while Ŝ_linear is a zero-mean term that reduces estimator variance when the experimental setup breaks rotational invariance, that is, in the presence of sky cut and anisotropic noise. To estimate Â^Tφ we used the KSW estimator with an implementation of the linear term truncated at ℓ_max as described in Munshi & Heavens (2010) and Planck Collaboration XXIV (2014).

3.1.3. Binned bispectrum

The binned bispectrum estimator (Bucher et al. 2010) achieves the required computational reduction in determining A^Tφ by binning Eq. (11). In particular, the maximally filtered maps in Eq. (10) are replaced by ${\mathit{T}}_{\mathit{i}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell }\mathrm{\in }{\mathrm{\Delta }}_{\mathit{i}}}\sum _{\mathit{m}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{\ell }}^{\mathrm{+}\mathit{\ell }}{\mathit{a}}_{\mathit{\ell m}}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\mathit{,}$(19)where the Δ_i are suitably chosen intervals (bins) of multipole values (chosen in such a way as to minimize the variance of the quantities to be estimated). These maps were then used in Eq. (9) to obtain the binned observed bispectrum, and analogously for B^lin. The bispectrum template B^{lens − ISW} and inverse-variance weights V were also binned by summing them over all ℓ values in the bin. Finally, these binned quantities were inserted into the general expression for A^Tφ (Eq. (11)), with the sum over ℓ replaced by a sum over bin indices i. Since most bispectrum shapes change rather slowly (with features on the scale of the acoustic peaks, like the power spectrum), the binned estimator works very well, increasing the variance only slightly, while achieving an enormous computational reduction (from about 2000 multipoles in each of the three directions to only about 50 bins).

3.1.4. Modal bispectra

Modal decomposition of bispectra has been introduced by Fergusson et al. (2010) as a way to compute reduced bispectra that uses a diagonalization ansatz such that the shape function in Fourier space can be separated, which reduces the dimensionality of the integration. At the same time, it greatly reduces the complexity of estimating bispectra from data. The separation of the bispectrum shape function into coefficients ${\mathit{q}}_{\mathit{p}}^{\mathit{\ell }}\mathrm{\left(}\mathit{x}\mathrm{\right)}$ allows the derivation of a filtered map , ${\mathit{M}}_{\mathit{p}}\mathrm{\left(}\mathit{n̂}\mathit{,x}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell m}}\frac{{\mathit{q}}_{\mathit{p}}^{\mathit{\ell }}\mathrm{\left(}\mathit{x}\mathrm{\right)}{\mathit{a}}_{\mathit{\ell m}}}{{\mathit{C}}_{\mathit{\ell }}}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathit{,}$(20)from the coefficients a_ℓm of the temperature map. With that expression, one can obtain a mode expansion coefficient β, ${\mathit{\beta }}_{\mathit{prs}}\mathrm{=}\mathrm{\int }\mathrm{d\Omega }\mathrm{\int }{\mathit{x}}^{\mathrm{2}}\mathrm{d}\mathit{x}{\mathit{M}}_{\mathit{p}}\mathrm{\left(}\mathit{n̂}\mathit{,x}\mathrm{\right)}{\mathit{M}}_{\mathit{r}}\mathrm{\left(}\mathit{n̂}\mathit{,x}\mathrm{\right)}{\mathit{M}}_{\mathit{s}}\mathrm{\left(}\mathit{n̂}\mathit{,x}\mathrm{\right)}\mathit{.}$(21)With that decomposition, the estimator of the bispectrum assumes a particularly simple diagonal shape, $\mathit{Ŝ}\mathrm{=}\frac{\mathrm{6}}{\mathit{N}}{\mathrm{\Delta }}_{\mathrm{\Phi }}^{\mathrm{2}}\sum _{\mathit{prs}}{\mathit{\alpha }}_{\mathit{prs}}{\mathit{\beta }}_{\mathit{prs}}\mathit{,}$(22)where the α_prs are the equivalent coefficients obtained by performing the modal decomposition of the theoretical bispectrum shape function. The relation between modal bispectra and wavelet bispectra was derived by Regan et al. (2013).

3.2. Results

The detection of the ISW effect via the non-Gaussian signal induced by the gravitational lensing potential is summarized in Table 2. We provide the estimates of the ISW-lensing amplitude A^Tφ, its uncertainty σ_A, and the signal-to-noise ratio obtained with the different estimator pipelines described in Sect. 3.1. The estimators were applied to the official Planck CMB maps made using C-R, NILC, SEVEM, and SMICA (Planck Collaboration XII 2014). The quantity σ_A was obtained from 200 simulations representative of the analysed CMB data maps. These Monte Carlo simulations (FFP-6, see Planck Collaboration I 2014) account for the expected non-Gaussian ISW-lensing signal, according to the Planck best-fit model, and were passed through the different component separation pipelines, as described in Planck Collaboration XII (2014). A more detailed description of the lensed simulations can be found in Planck Collaboration XVII (2014). The mask used in the analysis was the combined Galactic and point-source common mask (U73, Planck Collaboration I 2014) with sky fraction f_sky = 0.73.

The KSW and the Tφ estimators show similar sensitivity, finding A^Tφ = 0.81 ± 0.31 and A^Tφ = 0.70 ± 0.28, respectively, from the SMICA CMB map, which corresponds to a significance at about the 2.5σ level. The modal and binned estimators are slightly less optimal, but give consistent results, which is consistent with the imperfect overlap of the modal estimator templates with the ISW-lensing signal; the ISW-lensing bispectrum has a rapidly oscillating shape in the squeezed limit and both the binned and modal estimates are better suited (and originally implemented) to deal with smooth bispectra of the type predicted by primordial inflationary theories. Since the correlation coefficient of the binned and modal ISW-lensing templates relative to the actual ISW-lensing bispectrum (Eq. (8)) is generally 0.8 <r< 0.9 (to be compared with r = 0.99, achieved by both estimators for local, equilateral, and orthogonal primordial templates, Planck Collaboration XXIV 2014), the corresponding estimator’s weights are expected to be about 20 % suboptimal, consistent with observations.

The Tφ estimator was also applied to the specific Planck lensing baseline, that is, the MV map, which is a noise-weighted combination of the 217 GHz and 143 GHz channel maps, previously cleaned from infrared contamination through subtraction of the 857 GHz map, taken as a dust template. From this map the lensing potential was recovered and then correlated with that potential field to estimate the amplitude A^Tφ. The official baseline adopts a more conservative high-pass filtering, such that bacause only multipoles ℓ ≥ 10 are considered, and the mask with f_sky = 0.7 was used. In this case, the ISW-lensing estimate is 0.78 ± 0.32 (a 2.4σ detection, where the error bars are obtained from 1000 simulations), as reported in the first subrow for Tφ in Table 2. The full multipole range is considered in the second subrow, obtaining about 7% better sensitivity.

According to all the estimators, the C-R CMB map provides lower significance for ISW-lensing, since its resolution is slightly lower than that of the other maps. NILC and SMICA exhibit a somewhat larger detection of the ISW signal, since they are the least noisy maps.

To explore the agreement among the different estimators, we performed a validation test based on 200 lensed simulations processed through the SMICA pipeline. The results are summarized in Table 3. For each pair of statistics, we provide the difference in amplitudes estimated for the data (ΔA^Tφ), the dispersion of the difference of amplitudes obtained from the simulations (s_A), the ratio between this dispersion and the highest corresponding sensitivities (η, according to Table 2), and the correlation coefficient (ρ). As can be seen from the table, the agreement among estimators is good and the discrepancies are only around 0.5σ, which is the expected scatter, given the correlation between the weights of different estimators discussed above. Overall, the bispectrum estimators provide a higher value of the amplitude A^Tφ than the Tφ estimator.

We also explored the joint estimation of the two bispectra that are expected to be found in the data: the ISW-lensing and the residual point sources. A detailed description of the non-Gaussian signal coming from point sources can be found in Planck Collaboration XXIV (2014). The joint analysis of these two signals performed with the KSW estimator and the binned and modal estimators has shown that the ISW-lensing amplitude estimation can be considered almost completely independent of the non-Gaussian signal induced by the residual sources, and that the two bispectra are nearly perfectly uncorrelated.

There is no unique way of extracting a single signal-to-noise value from Table 2. However, all the estimators show evidence of ISW-lensing at about the 2.5σ level.

Finally, we estimate that the bias introduced by the ISW-lensing signal on the estimation of the primordial local shape bispectrum (Eq. (14)) is Δ^prim ≃ 7, corresponding to the theoretical expectation, as described in detail in Planck Collaboration XXIV (2014).

4. Cross-correlation with surveys

The ISW effect can be probed through several different approaches. Among those previously explored in the literature, the classical test is to study the cross-correlation of the CMB temperature fluctuations with a tracer of the matter distribution, typically a galaxy or cluster catalogue. As mentioned in the introduction, the correlation of the CMB with LSS tracers was first proposed by Crittenden & Turok (1996) as a natural way to amplify the ISW signal, which is otherwise very subdominant with respect to the primordial CMB fluctuations. Indeed, this technique led to the first reported detection of the ISW effect (Boughn & Crittenden 2004).

Several methods have been proposed in the literature to study statistically the cross-correlation of the CMB fluctuations with LSS tracers, and they can be divided into real-space statistics (e.g., the cross-correlation function, CCF), harmonic space statistics (e.g., the cross-angular power spectrum, CAPS), and wavelet space statistics (e.g., the covariance of the spherical mexican hat wavelet coefficients, SMHWcov). These statistics are equivalent (in the sense of the significance of the ISW detection) under ideal conditions. However, the ISW data analysis presents several problems (incomplete sky coverage, selection biases in the LSS catalogues, foreground residuals in the CMB map, etc.). Hence, the use of several different statistical approaches provides a more reliable framework for studying the ISW-LSS cross-correlation, since different statistics may have different sensitivity to these systematic effects. The individual methods are described in more detail in Sect. 4.1.

In addition to the choice of specific statistical tool, the ISW cross-correlation can be studied from two different (and complementary) perspectives. On the one hand, we can determine the amplitude of the ISW signal, as well as the corresponding signal-to-noise ratio, by comparing the observed cross-correlation with the expected one. On the other hand, we can postulate a null hypothesis (i.e., that there is no correlation between the CMB and the LSS tracer) and study the probability of obtaining the observed cross-correlation. Whereas the former answers a question regarding the compatibility of the data with the ISW hypothesis (and provides an estimate of the signal-to-noise ratio associated with the observed signal), the latter tells us how incompatible the measured signal is with the no-correlation hypothesis, that is, against the presence of dark energy (assuming that the Universe is spatially flat). Obviously, both approaches can be extended to account for the cross-correlation signal obtained from several surveys at the same time. These two complementary tests are described in detail in Sect. 4.2, with the results presented in Sect. 4.3.

4.1. Cross-correlation statistics

We denote the expected cross-correlation of two signals (x and y) by ${\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}$, where a stands for a distance measure (e.g., the angular distance θ between two points in the sky, the multipole ℓ of the harmonic transformation, or the wavelet scale R). For simplicity, we assume that the two signals are given in terms of a fluctuation field (i.e., with zero mean and dimensionless).

This cross-correlation could represent either the CCF, the CAPS, or the SMHWcov. It has to be understood as a vector of a_max components, where a_max is the maximum number of considered distances. Obviously, when x ≡ y, ${\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}$ represents an auto-correlation. The specific forms for ${\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}$ and C_ξ^xy for the different cross-correlation statistics (CAPS, CCF, and SMHWcov) are given below.

4.1.1. Angular cross-power spectrum

The angular cross-power spectrum (CAPS) is a natural tool for studying the cross-correlation of the CMB fluctuations and tracers of the LSS. Under certain conditions, it provides a statistical tool with uncorrelated (full-sky coverage) or nearly uncorrelated (binned spectrum for incomplete sky coverage) components. Even the unbinned CAPS, estimated on incomplete signals, can be easily worked out, since the correlations are mostly related to the geometry of the mask. This is the case for the CAPS obtained through MASTER approach (e.g., Hivon et al. 2002; Hinshaw et al. 2003). Another approach is to work in the map domain, making use of a quadratic maximum-likelihood (QML) estimator (Tegmark 1997) for the CAPS (Padmanabhan et al. 2005; Schiavon et al. 2012). This approach is optimal, that is, leads to unbiased estimates for the CAPS with minimum error bars.

Pseudo-angular power spectrum

We denote the CAPS between the CMB field T(p) and an LSS tracer G(p) map (where p = (θ,φ) represents a given pixel) as: ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TG}}$ (i.e., ${\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}\mathrm{\equiv }{\mathit{C}}_{\mathit{\ell }}^{\mathit{TG}}$ for this cross-correlation estimator). In the full-sky case, an optimal estimator of the CAPS is given by: $\mathit{Ĉ}\begin{array}{c}\mathit{TG}\\ \mathit{\ell }\end{array}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}\sum _{\mathit{m}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{\ell }}^{\mathrm{+}\mathit{\ell }}{\mathit{t}}_{\mathit{\ell m}}{\mathit{g}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathit{,}$(23)where t_ℓm and g_ℓm are the spherical harmonic coefficients of the CMB and the LSS maps, respectively. This CAPS can be seen as a vector with ℓ_max components, where ℓ_max is the highest multipole considered in the analysis. Here we adopted 3N_side − 1, which suffices for ISW analysis, since it is know that most of the ISW signal is contained within ℓ ≲ 80 (Afshordi 2004; Hernández-Monteagudo 2008). When a mask Π(p) is applied to the maps, it acts as a weight that modifies the underlying harmonic coefficients. Now, we have $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{t}}_{\mathit{\ell m}}\end{array}$ and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{g}}_{\mathit{\ell m}}\end{array}$, where $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ {\mathit{t}}_{\mathit{\ell m}}\end{array}& \mathrm{=}& \int \mathrm{\int }\mathrm{d}\mathrm{(}\cos {\mathit{\theta }}^{\mathrm{)}}\mathrm{d}\mathit{\phi }\hspace{0.17em}\mathit{T}\left(\mathit{\theta ,\phi }\right)\mathrm{\Pi }\left(\mathit{\theta ,\phi }\right){\mathit{Y}}_{\mathit{\ell m}}^{\mathrm{\ast }}\left(\mathit{\theta ,\phi }\right)\mathit{,}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{g}}_{\mathit{\ell m}}\end{array}& \mathrm{=}& \end{array}$(24)and Y_ℓm(θ,φ) are the spherical harmonic functions. In these circumstances, the estimator in Eq. (23) is no longer optimal and is referred to as pseudo-CAPS. A nearly optimal estimator is given by decoupling the masked CAPS (denoted by $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{TG}}\end{array}$) through the masking kernel B (e.g., Xia et al. 2011): $\mathit{Ĉ}\begin{array}{c}\mathit{TG}\\ \mathit{\ell }\end{array}\mathrm{=}{\mathrm{B}}^{-1}\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathit{TG}}}_{\mathrm{˜}}}\mathit{,}$(25)where ${\mathit{B}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,G}}\mathrm{=}\frac{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}{\mathrm{4}\mathit{\pi }}\sum _{{\mathit{\ell }}^{\mathrm{\prime \prime }}}{\mathit{J}}_{{\mathit{\ell }}^{\mathrm{\prime \prime }}}^{\mathit{G}}{\left(\begin{array}{c}\\ \mathit{\ell }& {\mathit{\ell }}^{\mathrm{\prime }}& {\mathit{\ell }}^{\mathrm{\prime \prime }}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)}^{\mathrm{2}}\mathit{,}$(26)with ${\mathit{J}}_{{\mathit{\ell }}^{\mathrm{\prime \prime }}}^{\mathit{G}}$ the cross-angular power spectrum of the T and G masks.

The estimator in Eq. (25) is nearly optimal because $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{TG}}\end{array}$ has to be understood as the mean value over an ensemble average of skies. When more than a single CAPS is considered, for instance, when one is interested in the cross-correlation of the Planck CMB map with more than one LSS tracer map, the CAPS estimator can be seen as a single vector with Nℓ_max components, with N being the number of surveys.

It can be shown that the element C_ℓℓ′,ij of the covariance matrix of the CAPS estimator in Eq. (25) (for the case of a masked sky and for N surveys) is given by ${\mathit{C}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,ij}}\mathrm{=}{\mathit{K}}_{\mathit{\ell ,ij}}{\mathit{K}}_{{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,ij}}\frac{{\mathrm{\left(}{\mathit{M}}_{\mathit{ij}}\mathrm{\right)}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}^{-1}}{\mathrm{(}\mathrm{2}{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{+}\mathrm{1}\mathrm{)}}\mathit{,}$(27)where ${\mathit{K}}_{\mathit{\ell ,ij}}\mathrm{=}{\left[{\mathit{C}}_{\mathit{\ell }}^{{\mathit{TG}}_{\mathit{i}}}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{TG}}_{\mathit{j}}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathit{T}}\left({\mathit{C}}_{\mathit{\ell }}^{{\mathit{G}}_{\mathit{i}}{\mathit{G}}_{\mathit{j}}}\mathrm{+}{\mathit{N}}_{\mathit{\ell }}^{{\mathit{G}}_{\mathit{i}}{\mathit{G}}_{\mathit{j}}}{\mathit{\delta }}_{\mathit{ij}}\right)\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(28)and ${{}^{\mathrm{(}}\mathit{M}_{\mathit{ij}}{\mathrm{)}}^{}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}^{-1}$ is the (ℓ,ℓ′) element of the inverse matrix of M for surveys i and j fixed, such as ${\mathit{M}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,ij}}\mathrm{=}\frac{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}{\mathrm{4}\mathit{\pi }}\sum _{{\mathit{\ell }}^{\mathrm{\prime \prime }}}{\mathit{H}}_{{\mathit{\ell }}^{\mathrm{\prime \prime }}}^{\mathit{ij}}{\left(\begin{array}{c}\\ \mathit{\ell }& {\mathit{\ell }}^{\mathrm{\prime }}& {\mathit{\ell }}^{\mathrm{\prime \prime }}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)}^{\mathrm{2}}\mathit{.}$(29)Here ${\mathit{H}}_{{\mathit{\ell }}^{\mathrm{\prime \prime }}}^{\mathit{ij}}$ is the angular cross-power spectrum of the two joint masks, that is, the masks resulting from the multiplication of the T with G_i and G_j, respectively. The quantities ${\mathit{C}}_{\mathit{\ell }}^{\mathit{xy}}$ are expected or fiducial spectra, ${\mathit{N}}_{\mathit{\ell }}^{\mathit{yy}}$ is the Poissonian noise of the y survey (deconvolved by any beam or pixel filter), and δ_ij is the Kronecker delta. In Eq. (28), the instrumental noise associated with the CMB data has been ignored, since the Planck sensitivity is such that the noise contribution on the scales of interest is negligible. When more than one survey has poor sky coverage, then the complexity of the correlations is not well reflected by the previous expression. Therefore, we computed C_ℓℓ′,ij from coherent CMB and LSS Monte Carlo simulations. For each simulation, we generated four independent, Gaussian, and white realizations (at N_nside = 64), which were afterwards properly correlated using the expected auto- and cross-correlations of the signals. Corresponding Poissonian shot-noise realizations were added to each survey map. The resulting four maps were masked with the corresponding masks (i.e., one for the CMB and one for mask for each survey).

The computation of the CAPS in Eq. (25) is extremely fast (especially for the resolutions used in the study of the ISW). However, as stated above, it is only a nearly optimal estimator of the underlying CAPS. Moreover, its departure from optimality is largest at the smallest multipoles (largest scales), where the ISW signal is more important.

The QML angular power spectrum

The QML method for the power spectrum estimation of temperature CMB anisotropies was introduced by Tegmark (1997) and later extended to polarization by Tegmark & de Oliveira-Costa (2001). For an application to temperature and polarization to WMAP data see Gruppuso et al. (2009) and Paci et al. (2013). The same method was employed to measure the cross-correlation between the CMB and LSS in Padmanabhan et al. (2005), Ho et al. (2008), and Schiavon et al. (2012). The QML method is usually stated to be optimal, since it provides unbiased estimates and the smallest error bars allowed by the Fisher-Cramér-Rao inequality. As a drawback from the computational point of view, the QML is a very expensive approach. We denote the QML estimator of the CAPS between the CMB map T and an LSS tracer G (at multipole ℓ) by ${\mathit{Q}}_{\mathit{\ell }}^{\mathit{TG}}$ (i.e., ${\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}\mathrm{\equiv }{\mathit{Q}}_{\mathit{\ell }}^{\mathit{TG}}$ for this cross-correlation estimator).

A detailed description of the algebra of the QML is given in Schiavon et al. (2012). We briefly recall here the basics of the CAPS estimator, which is given by $\mathit{Q̂}\begin{array}{c}\mathit{TG}\\ \mathit{\ell }\end{array}\mathrm{=}\sum _{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{X}}^{\mathrm{\prime }}}\mathrm{(}{\mathit{F}}^{-1}{\mathrm{)}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathit{TG}{\mathit{X}}^{\mathrm{\prime }}}\mathrm{[}{{x}}^{\mathit{T}}{\mathrm{E}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{{\mathit{X}}^{\mathrm{\prime }}}{x}\mathrm{-}\mathrm{Tr}\mathrm{(}\mathrm{N}{\mathrm{E}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{{\mathit{X}}^{\mathrm{\prime }}}{\mathrm{\left)}}^{\mathrm{\right]}}\mathit{,}$(30)where X represents any of the following spectra: X = {T,TG,G}. The vector x has 2N_pix elements (with N_pix being the total number of pixels allowed by the joint CMB and LSS mask): the first set of N_pix corresponds to the CMB map, and the second one accounts for the LSS map. The ${\mathit{F}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathit{TG}{\mathit{X}}^{\mathrm{\prime }}}$ is the Fisher matrix defined as ${\mathit{F}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathit{TG}{\mathit{X}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{Tr}\left[{\mathrm{C}}^{-1}\frac{\mathit{\partial }\mathrm{C}}{\mathit{\partial }{\mathit{C}}_{\mathit{\ell }}^{\mathit{TG}}}{\mathrm{C}}^{-1}\frac{\mathit{\partial }\mathrm{C}}{\mathit{\partial }{\mathit{C}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{{\mathit{X}}^{\mathrm{\prime }}}}\right]\mathit{,}$(31)and the E matrix is given by ${\mathrm{E}}_{\mathit{\ell }}^{\mathit{X}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{C}}^{-1}\frac{\mathit{\partial }\mathrm{C}}{\mathit{\partial }{\mathit{C}}_{\mathit{\ell }}^{\mathit{X}}}{\mathrm{C}}^{-1}\mathit{.}$(32)The object $\mathrm{C}\mathrm{=}\mathrm{S}\mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}^{\mathit{X}}\mathrm{\right)}\mathrm{+}\mathrm{N}$ is the global covariance matrix, including the signal S and noise N contributions, with ${\mathit{C}}_{\mathit{\ell }}^{\mathit{X}}$ being the fiducial theoretical angular power spectrum. The uncertainties in the QML estimates are given by the inverse of the Fisher matrix, which includes the correlation among different multipoles. The error associated with the shot noise of the galaxy surveys is modelled in the galaxy submatrix of N.

The results presented here are based on $\mathit{Ĉ}\begin{array}{c}\mathit{TG}\\ \mathit{\ell }\end{array}$, while $\mathit{Q̂}\begin{array}{c}\mathit{TG}\\ \mathit{\ell }\end{array}$ is used as a cross-check, applied to a lower-resolution version of the maps of N_side = 32. In addition, the maximum multipole considered in this case is ℓ_max = 2N_side, which has been already verified as a conservative limit for the QML.

4.1.2. Cross-correlation function

The cross-correlation function (CCF) is a suitable tool for studying the ISW effect via cross-correlation of the CMB fluctuations and tracers of the LSS, and it has been one of the most extensively used in this context (e.g., Boughn & Crittenden 2002; Giannantonio et al. 2008; Xia et al. 2009). On the one hand, the signal only appears at very large scales and, therefore, it is sufficient to work at resolutions at which the low performance of the CCF (in terms of computational time) is not a serious handicap. On the other hand, neither the CMB nor the LSS data are available with full sky coverage and, in some cases, the geometry of the masks is very complicated: the CCF adapts perfectly to the effects of partial sky coverage, since it is defined in real space. As a drawback, the Poissonian noise of the galaxy tracer appears at the smallest angular scales, where the signal-to-noise ratio of the ISW effect is higher for this estimator. Therefore, a proper characterization of the shot-noise is especially important for the CCF to obtain a good estimation of the uncertainties.

We denote the CCF between the CMB map T and an LSS tracer G (at an angular distance of θ) as C^TG(θ) (i.e., ${\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}\mathrm{\equiv }{\mathit{C}}^{\mathit{TG}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}$ for this cross-correlation statistic). The CCF estimator is defined as $\mathit{Ĉ}\mathit{TG}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}_{\mathit{\theta }}}\sum _{\mathit{i,j}}{\mathit{T}}_{\mathit{i}}{\mathit{G}}_{\mathit{j}}\mathit{,}$(33)where the sum runs over all pixels with a given angular separation. For each angular bin centred around θ, N_θ is the number of pixel pairs separated by an angle within the bin. Only the pixels allowed by the joint CMB and LSS mask are considered. The angular bins used here are θ₁ ∈ [0,1]°; θ₂ ∈ ⁽1,3^]°; θ₃ ∈ (3,5]°; ...; and θ₆₁ ∈ ⁽119,121^]°. The choice of binning does not affect the results significantly.

The covariance of the CCF estimator can be easily derived from the one already computed for the CAPS in Eq. (27). It is sufficient to know that the CCF can be expressed in terms of the CAPS as ${\mathit{C}}^{\mathit{TG}}\mathrm{(}{\mathit{\theta }}^{\mathrm{)}}\mathrm{=}\sum _{\mathit{\ell }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}^{{\mathit{\ell }}_{\max }}\frac{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathit{C}}_{\mathit{\ell }}^{\mathit{TG}}{\mathit{P}}_{\mathit{\ell }}\mathrm{(}\cos {\mathit{\theta }}^{\mathrm{)}}\mathit{,}$(34)where P_ℓ⁽cosθ⁾ are the Legendre polynomials. Hence, it is straightforward to prove that the covariance between the θ and θ′ components of the CCF for the surveys i and j, respectively, is given by ${\mathit{C}}_{\mathit{\theta }{\mathit{\theta }}^{\mathrm{\prime }}\mathit{,ij}}\mathrm{=}\sum _{\mathit{\ell }}\sum _{{\mathit{\ell }}^{\mathrm{\prime }}}\frac{\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}{\mathrm{1}}^{\mathrm{)}}}{\mathrm{4}\mathit{\pi }}\frac{\left(\mathrm{2}{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{+}\mathrm{1}\right)}{\mathrm{4}\mathit{\pi }}{\mathit{P}}_{\mathit{\ell }}\mathrm{(}\cos {\mathit{\theta }}^{\mathrm{)}}{\mathit{P}}_{{\mathit{\ell }}^{\mathrm{\prime }}}\left(\cos {\mathit{\theta }}^{\mathrm{\prime }}\right){\mathit{C}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,ij}}\mathit{.}$(35)

4.1.3. Wavelet covariance

Wavelets provide an interesting alternative to more traditional tools (e.g., CCF or CAPS) for studying the CMB-LSS correlation. They exploit the fact that the ISW signal is mostly concentrated at scales of a few degrees (e.g., Afshordi 2004). Wavelets are ideal kernels to enhance features with a characteristic size, since the wavelet analysis at an appropriate scale R amplifies those features over the background. Therefore, wavelets could provide most of the signal-to-noise ratio of the ISW effect by just analysing a narrow range of scales. They were first proposed for the ISW detection by Vielva et al. (2006), where the SMHW (Martínez-González et al. 2002) was proposed as the filtering kernel. The basic idea of this approach is to estimate the covariance of the SMHW coefficients (SMHWcov) as a function of the wavelet scale (see e.g., Vielva et al. 2006, for details). Other wavelet kernels can be considered, such as needlets (Pietrobon et al. 2006a), directional wavelets (McEwen et al. 2007), or steerable wavelets (McEwen et al. 2008).

We denote the SMHWcov between the CMB map T and a LSS tracer G (at a wavelet scale R) as ${\mathrm{\Omega }}^{\mathit{TG}}{}^{\mathrm{\right(}}\mathit{R}^{\mathrm{\left)}}$, i.e., ${\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}\mathrm{\equiv }{\mathrm{\Omega }}^{\mathit{TG}}{}^{\mathrm{\right(}}\mathit{R}^{\mathrm{\left)}}$ for this cross-correlation statistic). The SMHWcov estimator is defined as $\mathrm{\Omega ̂}\mathit{TG}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{pix}}}\sum _{\mathit{i}}{\mathit{\omega }}_{{\mathrm{T}}_{\mathit{i}}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}{\mathit{\omega }}_{{\mathrm{G}}_{\mathit{i}}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathit{,}$(36)where ${\mathit{\omega }}_{\mathit{T}}{}^{\mathrm{\right(}}\mathit{R}^{\mathrm{\left)}}$ and ${\mathit{\omega }}_{\mathit{G}}{}^{\mathrm{\right(}}\mathit{R}^{\mathrm{\left)}}$ are the SMHW coefficients for the CMB and the LSS at scale R, respectively (note that wavelet coefficients are forced to have zero mean on the observed sky). The scales considered in our study are R = { 60, 90, 120, 150, 200, 250, 300, 350, 400, 500, 600 } in arcminutes.

As for the CCF, the covariance of the SMHWcov estimator can be easily derived from the one already computed for the CAPS in Eq. (27). It is sufficient to know that the SMHWcov can be expressed in terms of the CAPS as ${\mathrm{\Omega }}^{\mathit{TG}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathrm{=}\sum _{\mathit{\ell }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}^{{\mathit{\ell }}_{\max }}\frac{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathit{C}}_{\mathit{\ell }}^{\mathit{TG}}{\mathit{\omega }}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathit{,}$(37)where ${\mathit{\omega }}_{\mathit{\ell }}{}^{\mathrm{\right(}}\mathit{R}^{\mathrm{\left)}}$ is the SMHW window function at the scale R. Hence, it is straightforward to prove that the covariance between the R and R′ components of the SMWHcov for surveys i and j, respectively, is given by ${\mathit{C}}_{\mathit{R}{\mathit{R}}^{\mathrm{\prime }}\mathit{,ij}}\mathrm{=}\sum _{\mathit{\ell }}\sum _{{\mathit{\ell }}^{\mathrm{\prime }}}\frac{\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}{\mathrm{1}}^{\mathrm{)}}}{\mathrm{4}\mathit{\pi }}\frac{\left(\mathrm{2}{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{+}\mathrm{1}\right)}{\mathrm{4}\mathit{\pi }}{\mathit{\omega }}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}{\mathit{\omega }}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathrm{2}}\left({\mathit{R}}^{\mathrm{\prime }}\right){\mathit{C}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,ij}}\mathit{.}$(38)

4.2. Cross-correlation tests

For any of the cross-correlation estimators described above, we performed two different statistical tests. First, if the observed cross-correlation is given by $\mathit{\xi ̂}\begin{array}{c}\mathit{xy}\\ \mathit{a}\end{array}$, a simple χ² can be proposed to estimate the amplitude A, such that $\mathit{A}{\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}$ is the closest solution to $\mathit{\xi ̂}\begin{array}{c}\mathit{xy}\\ \mathit{a}\end{array}$: ${\mathit{\chi }}^{\mathrm{2}}\mathrm{(}{\mathit{A}}^{\mathrm{)}}\mathrm{=}{\mathrm{[}\mathit{\xi ̂}\begin{array}{c}\mathit{xy}\\ \mathit{a}\end{array}\mathrm{-}\mathit{A}{{\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}}^{\mathrm{]}}}^{\mathit{T}}{\mathrm{C}}_{{\mathit{\xi }}^{\mathit{xy}}}^{-1}\mathrm{[}\mathit{\xi ̂}\begin{array}{c}\mathit{xy}\\ \mathit{a}\end{array}\mathrm{-}\mathit{A}{{\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}}^{\mathrm{]}}\mathit{,}$(39)where C_ξ^xy is the covariance matrix (of dimension a_max × a_max) of the expected cross-correlation ${\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}$, that is, ${{C}}_{{{\mathit{\xi }}^{\mathit{xy}}}_{\mathit{i,j}}}\mathrm{\equiv }{}^{\mathrm{⟨}}\mathit{\xi }_{{\mathit{a}}_{\mathit{i}}}^{\mathit{xy}}{{\mathit{\xi }}_{{\mathit{a}}_{\mathit{j}}}^{\mathit{xy}}}^{\mathrm{⟩}}$. It is straightforward to show that the best-fit amplitude A and its dispersion are given by

$\begin{array}{ccc}& & \mathit{A}\mathrm{=}{\mathrm{[}\mathit{\xi ̂}{\begin{array}{c}\mathit{xy}\\ \mathit{a}\end{array}}^{\mathrm{]}}}^{\mathit{t}}{\mathrm{C}}_{{\mathit{\xi }}^{\mathit{xy}}}^{-1}{\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}{\left[{\mathrm{[}{{\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}}^{\mathrm{]}}}^{\mathit{t}}{\mathrm{C}}_{{\mathit{\xi }}^{\mathit{xy}}}^{-1}{\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}\right]}^{-1}\mathit{,}\\ & & {\mathit{\sigma }}_{\mathit{A}}\mathrm{=}{\left[{\mathrm{[}{{\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}}^{\mathrm{]}}}^{\mathit{t}}{\mathrm{C}}_{{\mathit{\xi }}^{\mathit{xy}}}^{-1}{\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}\end{array}$(40)

An analogy with Eq. (39) can be defined for the null hypothesis case: ${{\mathit{\chi }}^{\mathrm{2}}}_{\mathrm{null}}\mathrm{=}{\mathrm{[}\mathit{\xi ̂}{\begin{array}{c}\mathit{xy}\\ \mathit{a}\end{array}}^{\mathrm{]}}}^{\mathit{T}}{\mathrm{D}}_{{\mathit{\xi }}^{\mathit{xy}}}^{-1}\mathit{\xi ̂}\begin{array}{c}\mathit{xy}\\ \mathit{a}\end{array}\mathit{,}$(41)where D_ξ^xy is the covariance of the cross-correlation of the two signals, in the absence of an intrinsic dependence, that is, when ${\mathit{\xi }}_{\mathit{a}}^{\mathit{xy}}\mathrm{\equiv }\mathrm{0}$. The ISW signal is very weak and, therefore it is a good approximation to assume that D_ξ^xy ≈ C_ξ^xy.

For Gaussian statistics, χ²_null already provides the direct probability of the observed cross-correlation $\mathit{\xi ̂}\begin{array}{c}\mathit{xy}\\ \mathit{a}\end{array}$ under the null hypothesis. However, several problems (sky coverage, systematics, foregrounds residuals, etc.) forced us to use alternative approaches to estimate the probability. One of the most common options is to perform the cross-correlation of survey signal y with realistic simulations of the CMB, x, and compute a joint statistics (e.g., χ²_null) for each simulation. The probabiltiy value associated with the data is then the fraction of simulations with a value of χ²_null equal to or higher than the one obtained for the data. Both C_ξ^xy and D_ξ^xy can be derived either analytically or numerically (via simulations).

The latter approach is computationally expensive, but in some cases could provide a more accurate defence against certain systematics, in particular, the incomplete sky coverage. There are several options to perform this type of simulations. The standard approach is the one mentioned above, that is, cross-correlation of the LSS map with CMB simulations. This is a very reliable approach, since it is usually hard to reproduce the systematics present in the LSS tracers, but incomplete because the LSS is fixed, resulting in a lack of randomness. An alternative method is to use a jack-knife test, which unfortunately tends to underestimate the errors. Finally, one can produce simulations of both the CMB and the LSS, assuming perfect knowledge of the properties of both signals, in particular of the LSS field (which, as mentioned above, is almost never the case). Comprehensive discussions of these approaches are given in Cabré et al. (2007) and Giannantonio et al. (2008).

4.3. Results

In this section we present the results obtained from the cross-correlation of the galaxy catalogues described in Sect. 2.2 (NVSS, SDSS-CMASS/LOWZand SDSS-MphG) with the four Planck CMB maps presented in Sect. 2.1.1 (C-R, NILC, SEVEM, and SMICA). All these maps were analysed at a HEALPix resolution of N_side = 64. The cross-correlation estimators described in the previous section were applied to all cases. This comprehensive analysis will help to achieve a reliable estimate of the ISW.

As already mentioned, the covariance among all the components of the estimators were obtained from coherent Gaussian simulations of the CMB and the three galaxy catalogues. Since we only considered large-scale effects (above about 1°), the same set of CMB simulations are equally valid for the four CMB maps, since they are nearly identical on such scales (see Planck Collaboration XII 2014, for details). We used 70 000 coherent Monte Carlo simulation sets (as described in Sect. 4.1.1) to compute the correlations; this is enough to characterize the covariance.

The expected signal-to-noise ratio for the ISW effect detection is summarized in Table 4. Values for all the cross-correlation estimators are given. We considered a survey-by-survey detection as well as the joint analysis of all the surveys. A signal-to-noise ratio of about 3σ is expected for the joint analysis, which is dominated by the NVSS cross-correlation. This is expected, since, firstly, NVSS covers a much larger portion of the sky compared to other surveys, and secondly, it extends over a redshift interval ideal for the detection of the ISW signal (e.g. Afshordi 2004).

The differences among estimators are not significant, indicating that none of them is clearly better than the others. To explore this agreement further, we analysed an extra set of 1000 CMB and LSS clustered simulations and compared, simulation by simulation, the ISW amplitude estimation derived for each cross-correlation estimator (${\mathit{C}}_{\mathit{\ell }}^{\mathit{TG}}$, ${\mathit{C}}^{\mathit{TG}}{}^{\mathrm{\right(}}\mathit{\theta }^{\mathrm{\left)}}$ and ${\mathrm{\Omega }}^{\mathit{TG}}{}^{\mathrm{\right(}}\mathit{R}^{\mathrm{\left)}}$). In Table 5 we summarize the comparison. We only report values for the joint fit to the ISW amplitude for the three surveys. Similar results are found survey by survey. For each pair of estimators, we provide the mean difference among the amplitude estimations (ΔA), the dispersion of these differences (s_A), the ratio (η) of this dispersion to the expected sensitivity (i.e., the inverse of the signal-to-noise ratio numbers given in the last column of Table 4), and the correlation coefficient (ρ). It is clear that the agreement between estimators is very high and that differences are, on average, lower than half the statistical uncertainty imposed by the sampling variance.

We fitted the observed cross-correlations to the expected ISW signal (${\mathit{C}}_{\mathit{\ell }}^{\mathit{TG}}$, ${\mathit{C}}^{\mathit{TG}}{}^{\mathrm{\right(}}\mathit{\theta }^{\mathrm{\left)}}$, and ${\mathrm{\Omega }}^{\mathit{TG}}{}^{\mathrm{\right(}}\mathit{R}^{\mathrm{\left)}}$, see Fig. 5), following Eq. (40), that is, allowing for a free amplitude of the expected signal. Results are summarized in Table 6. Note that the plots only show partial information on the total detection level of the ISW signal, specifically those from the diagonal of the covariance matrix C_ξ^xy. The level of correlation among the components (i.e., the multipoles, angles or scales) varies considerably for the different estimators, and, therefore the whole picture is only obtained when the full covariance is considered, as in Eq. (40).

Overall, the ISW detection is at about the 3σ level and, as expected, it is clearly dominated by the NVSS signal. There are only small differences among estimators and CMB maps (as expected from the above discussion), indicating that this is a reliable result. Note that all the estimated amplitudes are compatible with unity within the error bars (especially for NVSS and SDSS-CMASS/LOWZ). This is an additional validation of how CMB and LSS are modelled. Values of A that deviate significantly from unity would indicate some tension between the observed cross-correlation and the model (in particular on the LSS modelling, which is more complex). The CAPS-QML, applied to the SEVEM and NVSS (i.e., the survey with the highest signal-to-noise ratio), yields a value of A = 0.73 ± 0.33, which is compatible with the CAPS when applied to the same N_side = 32 and ℓ_max = 2N_side resolution (A = 0.84 ± 0.34). Preliminary tests indicated that running the CAPS-QML at N_side = 64 resolution may increase the sensitivity for detecting the ISW effect with NVSS by ≈20%.

Our results indicate a somewhat lower signal-to-noise ratio than in some previous analyses of WMAP data, where several (and in some case quite similar) surveys were also considered. For instance, Ho et al. (2008) and Giannantonio et al. (2012) found 3.7σ and 4.4σ detections, respectively. Compatibility with the former is below 1σ, while there is more tension (about 1.5σ) with the latter. A fraction of about 0.3σ of these differences might be explained in terms of the comological parameters adopted to defined the theoretical expectations. In particular, the lower values of H₀ and Ω_Λ found by Planck (Planck Collaboration XVI 2014) with respect to WMAP (e.g., Larson et al. 2011) imply a lower sensitivity for the ISW by ≈10%. The rest of the differences come either from the LSS side or from the error characterization, which depends on the presence of a correlated signal between CMB and LSS simulations (see for instance Cabré et al. 2007, for a discussion). Survey modelling is another important aspect: in addition to systematic errors associated with the galaxy identification and redshift estimation procedures, there are complicated aspects, such as the bias characterization. As was mentioned already, a strong point of our results is the excellent compatibility between the ISW amplitude estimates with respect to the expected value. While our estimation deviates by about 0.5σ from the expected value, the Giannantonio et al. (2012) result exceeds it by about 1σ and Ho et al. (2008) is higher by 2σ.

Nevertheless, the value of the ISW effect that we measured by means of NVSS (which is probably the best current catalogue for studying the ISW effect because of its large sky coverage, redshift range, and density of galaxies) is significant and agrees with previously published results using WMAP.

We also studied the ISW signal from the point of view of its compatibility with the null hypothesis. We considered in this analysis only the NVSS catalogue, since it provides the highest detection of the ISW effect and therefore is the best-suited of the existing surveys to challenge the null hypothesis. Probability values are summarized in Table 7. As mentioned in Sect. 4.2, there is no unique way of computing the null hypothesis. Our approach follows Eq. (41), where D_ξ^xy was computed out of 90 000 CMB simulations that were cross-correlated with the LSS data. This matrix was used to compute ${\mathit{\chi }}_{\mathrm{null}}^{\mathrm{2}}$ from the data. This value was then compared with its distribution for the null hypothesis, obtained from 1000 realistic CMB simulations (FFP-6) uncorrelated with NVSS, which were processed in the same way as the Planck data set. CAPS provides the lowest probability value, but the null hypothesis is rejected at about 10% only; this result is expected since the ISW effect is weak.

The fact that the CAPS statistic provides tighter limits than the CCF and SMHWcov could have been anticipated. In our implementation, the CAPS explores the maximum angular range allowed for a given map, while the CCF and the SMHWcov approaches are only evaluated at certain angles/scales. This limitation is not a problem in the analysis devoted to estimating the ISW amplitude, since these angles/scales are suitable for detecting the ISW. However, to discard the null hypothesis, the larger the number of “distances”, the better.

The previous approach is a frequentist one. However, there is an alternative way of addressing the null hypothesis compatibility in the framework of a hypotheses test. In particular, we can study the ratio of Bayesian evidence for both scenarios, comparing the alternative (there is ISW signal) and the null (there is no correlation between the CMB and the LSS) hypotheses. To compute the Bayesian evidence we just need the likelihoods and the priors associated with each hypothesis, which are already available. The likelihood for the alternative hypothesis is obtained from Eq. (39), peaked at the best-fit value for the ISW amplitude, and its prior can be described as a Gaussian probability peaked at A = 1 and with a dispersion given by σ_A in Eq. (40). Conversely, the likelihood for the null hypothesis is the function given by Eq. (41), and its prior would be like the one for the alternative hypothesis, but peaked at A ≡ 0 (notice that this is justified, since, as explained in the previous section, D_ξ^xy ≈ C_ξ^xy). Following this approach, we obtain a Bayesian evidence ratio of about 30, which provides strong support for the assumption that there is an ISW signal.

5. Stacking of large-scale structures

An alternative approach for measuring the ISW effect in Planck maps is to search for an ISW signal directly at the positions of positive and/or negative peaks in the potential. Since the expected (and observed) signal is very weak, a stacking technique needs to be applied for individual structures. Using the WMAP data, it has been shown that CMB maps show hot spots and cold spots in the direction of superclusters and supervoids, respectively (Granett et al. 2008a,b, GR08 hereafter), which appear to be barely consistent with the predictions of standard ΛCDM (see also Hernández-Monteagudo & Smith 2013). These structures, which are not yet virialized, are evolving while the CMB photons travel across them and this should contribute to the ISW effect. We applied here the same approach to the different Planck maps, using the catalogues of superstructures introduced in Sect. 2.2.4, and we tested the reliability of our findings. We first discuss our method and the results obtained using the catalogue provided by GR08 and then present the results obtained with the other catalogues.

5.1. Method

Our analysis was performed on the SMICA CMB map, although we verified that the results are compatible for the other three Planck maps. We also used the cleaned frequency maps from SEVEM (see Sect. 2.1.1) for some of the tests. We first removed the monopole and dipole of the maps (outside the U73 mask) and then applied a compact source mask based on the Planck Legacy Point Source Catalogue (Planck Collaboration XII 2014) to remove the contamination from individual point sources.

For the purpose of comparison with the results of GR08, we smoothed the CMB maps with a common Gaussian kernel of 30′ FWHM. We then projected them onto patches around each position in the supervoid and supercluster catalogue. The GR08 structures have a relatively small size on the sky (a few degrees), but the other two catalogues considered here contain many larger and closer voids, covering larger angular sizes. Therefore we worked with 30deg × 30deg CMB patches and chose the pixel size to be 6′, so that all voids considered are fully enclosed. We then co-added (stacked) the maps, taking into account the mask used. On the stacked images, we calculated both the radial temperature profile and the aperture photometry to characterize the signal around density structures. The temperature profile was obtained by computing the mean of the pixels in rings of fixed width and increasing angular radius; in practice, it was calculated for 150 radii between 0deg and 15deg, with a width of . The photometry profile was obtained by applying a compensated filter that subtracts the average temperature of a ring from the average temperature within the disk, whose radius θ is the inner radius ring, and where the outer radius was chosen to be $\mathit{\theta }\sqrt{\mathrm{2}}$, so that the disc and ring have the same area. This is expected to enhance fluctuations of typical angular size θ against fluctuations at smaller or larger scales. Aperture photometry results were also provided for at 150 angles, this time between 0deg and 15/$\sqrt{\mathrm{2}}\mathrm{\approx }\mathrm{10}.̊\mathrm{6}$. In addition to the monopole and dipole, we also removed from the CMB maps the contribution of large-scale angular modes, namely ℓ = 2 − 10. These modes correspond to angular scales much larger than those of the structures under investigation, and for our purposes their only effect is to introduce gradients in the stacked images; the high-pass filter essentially stops from gradients entering the stacked map (which is equivalent to removing gradients at the end). The contribution of the large-scale angular modes has no impact on the aperture photometry profiles and only introduces an offset in the temperature profiles (Ilić et al. 2013).

To estimate the significance of the resulting photometry and temperature profiles, we followed a Monte Carlo approach based on stacked CMB images chosen at random positions. In detail, we computed the photometry and the temperature profiles for 16 000 sets of 50 CMB patches randomly distributed over the SDSS area. We then compared the profiles obtained from the stacking at the location of the GR08 superstructures to these random profiles to compute their signal-to-noise ratio.

5.2. Results

We show in Fig. 6 the stacked images of the 50 supervoids and 50 superclusters of GR08 in the Planck map. The corresponding temperature and photometry profiles along with their significance levels are shown in Fig. 7. First, although the signatures are fairly weak, the sign of the effect certainly seems to be correct. Using the same catalogue and the Planck CMB map, we found a reasonable agreement with GR08. The highest photometric decrement −10.8 μK (essentially identical with the −11.3 μK found by GR08) induced by supervoids is obtained for a preferred scale of about (4deg in GR08) and a signal-to-noise ratio of 3.3 (3.7σ in GR08), as shown in Fig. 7. Superclusters produce a photometric increment of about 8.5 μK (slightly above the 7 μK in GR08), with a significance of 3.0σ (compared with 2.6σ in GR08) at a slightly larger angle of . Finally, the stack of the combined sample (clusters minus voids) gives a temperature deviation of 8.7 μK with a signal strength of 4.0σ at , which is consistent with the values reported in GR08. The values of statistical significance for our aperture photometry results are closely related to those for the temperature profiles. Indeed, as shown in the top panel of Fig. 7, the temperature profile for the void stack shows a deficit of about 2σ at small angular radii and an excess of about 2σ extending to large radii. Since the aperture photometry is essentially an integral of the temperature profile with a compensated filter, it picks up enhanced significance because of the shape of the temperature profile.

As noted previously by several authors (e.g., Hernández-Monteagudo & Smith 2013), the amplitude and shape of the photometric profile found for voids and clusters disagrees (around 2σ) with the values expected from pure ISW within ΛCDM. However, it is not straightforward to associate this entire signal with a pure ISW effect. As seen in the Fig. 6, many small-scale structures – both cold and hot – are present around the region delineated by the angular radius at which the signal-to-noise ratio of the aperture photometry is highest. This small-scale structure contributes to the amplitude of the photometric decrement, but at a few tens of μK, which is incompatible with the ΛCDM predictions for the ISW effect. These are rather simply background CMB fluctuations, and linger because of the small number (50) of patches that were used to produce the stack.

It is intriguing that the angular sizes of the catalogued superstructures are smaller than the angular sizes suggested by the photometry profiles. This result was more apparent when we repeated the stacking analysis after rescaling each CMB patch by the effective radius of the structure it contains. Since the voids and super clusters identified by GR08 are approximately the same size, the photometric results are similar after rescaling (−10 μK for voids and 7.9 μK for superclusters). However, the deviations have significance levels of 3.3σ and 2.7σ for supervoids and superclusters, respectively, at angular sizes of 0.9 (voids) and 2.6 (clusters) times the effective radius of the structures. This mismatch, especially for clusters, could be a result of underestimation of the structure extent of the ZOBOV and VOBOZ algorithms (as previously suggested by GR08) or because larger potential hills and valleys underlie the detected superstructures. Since structure in the potential is related to the density field through the Poisson equation, gravitational potential features are expected to cover larger scales than structures in the density field. Nevertheless, the factor of 2.6 for the case of superclusters seems high. It is also true that the GR08 superstructures were identified in the LRG subsample of the SDSS, and LRGs are known to be biased tracers of the matter density field (e.g., Tegmark et al. 2006). This bias might help to explain why structures are larger than the scales identified by the VOBOZ algorithm, although the argument would go in the opposite sense for the voids. Another way of stating this is (Hunt & Sarkar 2010) that the relatively strong effect decrement found for the GR08 voids can only be attributed to the ISW effect if the underdensities have been significantly underestimated or if the LRGs are under-biased.

It is therefore difficult to be confident that the signal is caused entirely by the ISW effect. We know, however, that the ISW signal generated by superstructures is expected to be achromatic, since it generates a fractional perturbation of the CMB temperature. To determine if the signal we measure is indeed independent of frequency, we applied the same technique to Planck individual SEVEM cleaned frequency maps from 44 to 353 GHz. Lower (higher) frequency maps may be contaminated by radio (IR) signals coming from our Galaxy and may thus introduce a bias in the measurement. Figure 8 shows the photometry profiles of supercluster- and supervoid-stacked maps at 44, 70, 100, 143, 217, and 353 GHz. The flux measured appears to be quite constant, which supports the idea that the signal is generated by the ISW effect induced by structures. In the remainder of this section, we therefore apply our analysis only to the SMICA CMB map.

5.3. Discussion and alternative catalogues

We that although the GR08 structures are considered to be good tracers of the cosmic matter distribution on scales larger than 10 h^-1 Mpc, they are also known for their sparsity at these redshifts (z ≈ 0.4–0.7). This sparsity could lead to biased estimates of the properties of the reconstructed voids, their sizes and depths in particular could be biased. Moreover, some of the structures overlap on the sky, which could lead to systematic effects in the stacking analysis.

We therefore turned to other samples, for example, to the catalogue of Pan et al. (2012) introduced in Sect. 2.2.4. The 1054 statistically significant voids it contains are larger than 10 h^-1 Mpc in radius and, with redshifts lower than 0.1, they are much closer to us than the structures of GR08. Direct stacking gives only a weak signal at about the 1σ level, which is difficult to reconcile with the previous results. This may be due to the inclusion of a large number of small voids that may dilute the signal. Moreover, unlike the voids of GR08, the voids of Pan et al. (2012) have a large scatter in angular sizes on the sky, from about 2deg to 20deg (e.g., Ilić et al. 2013). To try to enhance the signal, we repeated the stacking after rescaling the voids to their effective sizes. We also subdivided the catalogue into subsamples based on redshift, radius, and/or angular size. However, none of these attempts yielded any statistically significant result, in agreement with Ilić et al. (2013).

Finally, we applied our procedure to the catalogue of voids published by Sutter et al. (2012). These cover a rather extended range of angular scales (about 2–10deg), and accordingly we rescaled all the CMB patches by the effective radius of each void. Stacking subsample by subsample (dim1, dim2, bright1, bright2, lrgdim, lrgbright) did not yield any significant signal. Similarly, when stacked together, the entire catalogue does not yield an ISW detection. However, since the ISW signal is expected to be stronger for the largest and closest voids (e.g., Flender et al. 2013) we tried starting from the largest void and adding them one by one, looking for the optimal number of voids, that is, that for which the signal-to-noise ratio is highest (see Fig. 10). We found that stacking 231, 589, or 936 voids gives about the same signal-to-noise ratio (2.5σ, 2.0σ and 2.2σ, respectively). However, the more voids we stacked, the smaller was the amplitude of the photometry signal (see Fig. 9, this was about −2.0 μK for 936 voids, −2.1 μK for 589 voids and −4.1 μK for 231 voids, at an angular size of about 0.5 times the common rescaled radius. These amplitudes are lower than those found with the 50 GR08 voids, but higher than expected from numerical simulations (see e.g., Hernández-Monteagudo & Smith 2013, for higher redshift and larger voids), but see also Cai et al. (2013).

The apparent angular size detected (about 0.5 times the effective void radius) in the photometry profile is smaller than that for the GR08 voids, but agrees with expectations from simulations (Cai et al. 2013) and is consistent with the sizes detected using WMAP data (Ilić et al. 2013).

The profiles in Fig. 9 show intriguing hints (significance about 2σ) of a positive excess below about 0.2 times the effective void radius. This is somewhat surprising, since the centre of each void of the Sutter et al. (2012) catalogue has been defined as the volume-weighted barycentre of all the galaxies contained in the void volume. This definition in principle removes most of the bias that may arise in other methods, such as choosing the most underdense galaxy position, which would lead to a slight overdensity at the centre. Figure 9 also shows positive excess for larger apertures, partly caused by the large hot ring surrounding the cold feature in the stacked images, which raises the mean temperature of the stacked image for discs of radii of about 0.8−1.2 times the void radius.

6. ISW map recovery

In recent years, some effort has been invested not only to obtain the statistical cross-correlation signal between the CMB and LSS data, but also to recover a map of the ISW signal itself (Barreiro et al. 2008, 2013; Francis & Peacock 2010a; Dupé et al. 2011). In particular, assuming the existence of a correlation between the CMB and the gravitational potential, it is possible to recover a map of the ISW fluctuations using a filtering method when one has a tracer of the gravitational potential (e.g., the galaxy catalogues described in Sect. 2.2) and the CMB fluctuations. Given the weakness of the signal, the main objective of this section is to provide a qualitative image of the ISW fluctuations for visual inspection, and an additional consistency test of the validity of the assumed fiducial model, by comparing the statistical properties of the recovered and expected signals. In addition, this secondary anisotropy map can also be used to study the large-scale properties of the CMB and its possible relation to some possible large-angle anomalies found in the Planck data (Planck Collaboration XXIII 2014).

6.1. Method

We followed the methodology of Barreiro et al. (2008), which applies a linear filter to the CMB and to a gravitational potential tracer map to reconstruct an ISW map, assuming that the cross- and auto-spectra of the signals are known. This technique has recently been applied to reconstruct the ISW map from the WMAP data and NVSS galaxy map (Barreiro et al. 2013). The filter was implemented in harmonic space and the estimated ISW map ŝ_ℓm at each harmonic mode is given by (see Barreiro et al. 2008 for details) $\mathit{ŝ}\mathit{\ell m}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{12}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}}{{\mathit{L}}_{\mathrm{11}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}}{\mathit{g}}_{\mathit{\ell m}}\mathrm{+}\frac{{\mathit{L}}_{\mathrm{22}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}}{{\mathit{L}}_{\mathrm{22}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{n}}}\left({\mathit{d}}_{\mathit{\ell m}}\mathrm{-}\frac{{\mathit{L}}_{\mathrm{12}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}}{{\mathit{L}}_{\mathrm{11}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}}{\mathit{g}}_{\mathit{\ell m}}\right)\mathit{,}$(42)where L(ℓ) corresponds to the Cholesky decomposition of the covariance matrix between the considered tracer of the potential and the ISW signal at each multipole, which satisfies C(ℓ) = L(ℓ)L^T(ℓ). Here d_ℓm and g_ℓm are the CMB data and the gravitational potential tracer map, respectively, and ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{n}}$ is the power spectrum of the CMB signal excluding the ISW effect. If full-sky coverage is not available, the covariance matrix is obtained from the corresponding pseudo-spectra. It can be shown that the expected value of the power spectrum for the reconstructed signal is given by $\mathrm{⟨}{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{ŝ}}}^{\mathrm{⟩}}\mathrm{=}\frac{\mathrm{(}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{gs}}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\left|\mathit{C}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\right|\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{g}}{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{n}}}^{\mathrm{)}}\mathrm{+}{\left|\mathit{C}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\right|}^{\mathrm{2}}}{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{g}}\mathrm{(}\left|\mathit{C}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\right|\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{g}}{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{n}}}^{\mathrm{)}}}\mathit{,}$(43)where | C(ℓ) | is the determinant of the tracer-ISW covariance matrix at each multipole, and ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{gs}}$ and ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{g}}$ are the assumed cross-spectrum and gravitational potential tracer spectra, respectively. Note that the recovered ISW power spectrum will not contain the full ISW signal, since it can only account for the part of the ISW signal probed by the tracer considered. It is also worth noting that in detail the expected cross-correlation depends on the assumed model. However, in practice, given the weakness of the signal, it would be difficult to distinguish between two mild variants of the standard ΛCDM model. Nevertheless, this approach still provides a useful consistency check.

6.2. Results

We applied the filter described above to two different cases: combining information from the CMB and the NVSS galaxy catalogue, and applying the filter to the CMB and the recovered lensing potential map described in 2.1.2. Results were obtained for the four Planck maps, C-R, NILC, SEVEM, and SMICA. For simplicity, we show the reconstructions only for the SEVEM CMB map, since the four methods give very similar results. The resolution considered for both analyses is N_side = 64.

For the first case, we used the Planck fiducial model for the CMB and cross-power spectrum, while for the NVSS map we assumed the model described in Sect. 2.2.1. We also took into account the presence of Poissonian noise. We excluded the area obtained from combining the CMB mask at N_side = 64 (described in Sect. 2.1.1) and the area that has not been observed by NVSS. The final mask represents about 62% of the sky. Since the filter is constructed in harmonic space, we used an apodized version of the mask to reduce the mask-induced correlations. In any case, the degradation introduced by a mask is small (Barreiro et al. 2008).

For the second case, the lensing map involved applying a high-pass filter, which removed all multipoles with ℓ < 10. This filtering was made in harmonic space with a mask. To take this effect into account we used a direct estimate of the pseudo-power spectrum of these data for the power spectrum of the lensing map after applying the corresponding apodized mask. We used the Planck fiducial model for the other power spectra involved, but set the cross-power for ℓ < 10 to zero. A mask wasconstructed by combining the CMB mask plus that provided for the lensing potential map (described in Planck Collaboration XVII 2014), which represents about 67% of the sky. The corresponding apodized version of this mask was applied before reconstructing the ISW map. Note that the map given in Fig. 1 (right panel) corresponds, to a good approximation, to the first term of the right-hand side of Eq. (42).

Figure 11 shows the reconstructed ISW map using the Planck CMB map and NVSS (left panel) and that obtained combining the CMB with the lensing potential map (right panel). There are similar structures present in both maps, but they are not expected to look exactly the same, since each of them provides only a partial reconstruction of the ISW signal. This is because that the reconstruction accounts for the part of the ISW effect probed by the considered tracer, which is different (although correlated) for each case. Moreover, because of the high-pass filter applied to the lensing potential map, the power at ℓ < 10 for this case corresponds to the Wiener-filtered map of the CMB (to which the filter given by Eq. (42) defaults, if the cross-correlation is set to zero, as in this case), without additional information from the considered tracer.

For both cases, we tested that the power spectrum of the recovered ISW signal, as well as that of the cross-power between the reconstructed ISW and the considered gravitational potential tracer, are consistent with the corresponding expected values. This indicates the compatibility between the assumed fiducial model and the underlying statistical properties of the data.

7. Conclusions

We presented the first study of the ISW effect using Planck data. We derived results based on three different approaches: the detection of the interplay between weak lensing of the CMB and the ISW effect by looking at non-Gaussian signatures, the conventional cross-correlations with tracers of large-scale structure, and aperture photometry on stacks of the CMB field at the positions of known superstructures. A reconstruction of the ISW map inferred from the CMB and LSS tracers was also provided.

The correlation with lensing allows, for the first time, the detection of the ISW effect using only CMB data. This is an effective approach, because the gravitational potential responsible for deflecting CMB photons also generates ISW temperature perturbations. Using different estimators, we investigated the correlation of the Planck temperature map with a reconstruction of the lensing potential on the one hand, and the estimation of the ISW-lensing generated non-Gaussian signature on the other. We found that the signal strength is close to 2.5σ for several combinations of estimator implementation and foreground-cleaned CMB maps.

We computed cross-correlations between the Planck CMB temperature map and tracers of large-scale structure, namely the NVSS survey of radio sources; and the SDSS-CMASS/LOWZ, and SDSS-MphG galaxy samples. As estimators we considered the angular cross-correlation function, the angular cross-spectra, and the variance of wavelet coefficients as a function of angular scale. We performed a comparison on different component-separation maps, where we considered C-R, NILC, SEVEM, and SMICA, and found remarkable agreement between the results, indicating that the the low multipoles are reliably reconstructed. Covariance matrices between the cross-correlation quantities were estimated for a set of Gaussian realizations of the CMB for the Planck fiducial model. For the ISW effect, we report detection significance levels of 2.9σ (NVSS), 1.7σ (SDSS-CMASS/LOWZ), and 2.0σ (SDSS-MphG), which are consistent among the different estimators considered. Although these numbers are compatible with previous claims that used WMAP data, they are generally smaller. We believe that this discrepancy is mainly due to the different characterization of the surveys and treatment of uncertainties, since the measurement of the CMB fluctuations at the scales that contribute to the ISW detection are very similar for Planck and WMAP. Only a fraction of these differences (about 0.3σ) can be understood in terms of the different cosmological models used by each experiment – in particular, that the values of H₀ and Ω_Λ reported by Planck are lower than those of WMAP.

A strength of our new study lies in the fact that the amplitudes derived for the expected signals are largely consistent with unity (i.e., the model expectation), which indicates good modelling of the surveys. The CMB and LSS cross-correlation was also tested against the null hypothesis, that is, whether the observed signal was compatible with a null correlation. As expected for such a weak signal, there is no strong evidence of incompatibility with the lack of correlation. In this respect, the CAPS approach seems to provide better constraints than the other estimators investigated here (CCF and the SMHWcov).

We explored the aperture photometry of stacked CMB patches at the positions of superstructures identified in the SDSS galaxy distribution. Our analysis of the Granett et al. (2008a) catalogue (50 supervoids and 50 superstructures) reproduced previous results with similarly strong amplitude and significance levels (somewhat higher and lower than 3σ for voids and clusters, respectively). While the most plausible source of this signal is the ISW effect associated with these structures, it slightly disagrees with expectations, both in terms of amplitude and scale. The same type of analysis was carried on the latest and much larger void catalogues of Sutter et al. (2012) (about 1500 voids) and Pan et al. (2012) (about 1000 voids). The results range from negligible to evidence at the 2−2.5σ level, with a more moderate amplitude and a smaller scale, in better agreement with theoretical predictions found in the literature. The broad spectral coverage of Planck allows us to confirm the achromatic nature of these signals across the 44 to 353 GHz range, which supports the idea that they have a cosmological origin.

We reconstructed maps of the ISW effect using a linear filter by combining the Planck CMB and a gravitational potential tracers. In particular, we considered both the NVSS catalogue and the reconstructed CMB lensing map as LSS tracers. Again we found good agreement between different component separation methods and consistency between the expected and reconstructed auto- and cross-power spectra for the recovered ISW map.

We conclude that the ISW effect is present in Planck data at the level expected for ΛCDM-cosmologies, using a range of measurement methods, although there is a possible disagreement with the results from stacking of CMB fields centred on superstructures. Generally, our results are more conservative than previous claims using WMAP data, but the agreement with the expected signal is better. Future Planck data releases, including polarization information and improved understanding of foregrounds, may be able to improve on these results, in particular for ISW-lensing correlation and ISW-lensing map reconstruction.

---

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, 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); and PRACE (EU). A description of the Planck Collaboration and a list of its members, including the technical or scientific activities in which they have been involved, can be found at http://www.sciops.esa.int/index.php?project=planck&page=Planck_Collaboration. The modal and KSW bispectrum estimator analysis was performed on the COSMOS supercomputer, part of the STFC DiRAC HPC Facility. We acknowledge the computer resources, technical expertise and assistance provided by the Spanish Supercomputing Network (RES) node at Universidad de Cantabria, and the support provided by the Advanced Computing and e-Science team at IFCA.

References

All Tables

All Figures