© ESO, 2014

1. Introduction

This paper, one of a set associated with the 2013 release of data from the Planck¹ mission (Planck Collaboration I 2014), presents a first detection of a strong correlation between the infrared background anisotropies and a lensing-derived projected mass map. The broad frequency coverage of the Planck satellite provides two important probes of the high redshift Universe. In the central frequency bands of Planck (70, 100, 143, and 217 GHz), cosmic microwave background (CMB) fluctuations dominate over most of the sky. Gravitational lensing by large-scale structure produces weak shear and magnification effects in the observed fluctuations, which can be exploited to reconstruct an integrated measure of the gravitational potential along the line of sight Okamoto & Hu (2003). This “CMB lensing potential” is sourced primarily by dark matter halos located at 1 ≲ z ≲ 3, halfway between ourselves and the last scattering surface (see Blandford & Jaroszynski 1981; Blanchard & Schneider 1987; or Lewis & Challinor 2006, for a review). In the upper frequency bands (353, 545, and 857 GHz), the dominant extragalactic signal is not the CMB, but the cosmic infrared background (CIB), composed of redshifted thermal radiation from UV-heated dust, enshrouding young stars. The CIB contains much of the energy from processes involved in structure formation. According to current models, the dusty star-forming galaxies (DSFGs), which form the CIB have a redshift distribution peaked between z ∼ 1 and z ∼ 2, and tend to live in 10¹¹–10¹³ M_⊙ dark matter halos (see, e.g., Béthermin et al. 2012a, and references therein).

As first pointed out by Song et al. (2003), the halo mass and redshift dependence of the CMB lensing potential and the CIB fluctuations are well matched, and as such a significant correlation between the two is expected. This point is illustrated quantitatively in Fig. 1, where we plot estimates for the redshift – and mass – kernels of the two tracers. In this paper we report on the first detection of this correlation.

Measurements of both CMB lensing and CIB fluctuations are currently undergoing a period of rapid development. While the CIB mean was first detected using the FIRAS and DIRBE instruments aboard COBE (Puget et al. 1996; Fixsen et al. 1998; Hauser et al. 1998), CIB fluctuations were only later detected by the Spitzer Space Telescope (Lagache et al. 2007) and then by the BLAST balloon experiment (Viero et al. 2009) and the Herschel Space Observatory (Amblard et al. 2011; Viero et al. 2013), as well as the new generation of CMB experiments, including Planck, which have extended these measurements to longer wavelengths (Hall et al. 2010; Dunkley et al. 2011; Planck Collaboration XVIII 2011; Reichardt et al. 2012). The Planck early results paper (Planck Collaboration XVIII 2011, henceforth referred to as PER) presented measurements of the angular power spectra of CIB anisotropies from arc-minute to degree scales at 217, 353, 545, and 857 GHz, establishing Planck as a potent probe of the clustering of the CIB, both in the linear and nonlinear regimes. A substantial extension of PER is presented in a companion paper to this work (Planck Collaboration IX 2014, henceforth referred to as P2013).

The CMB lensing potential, on the other hand, which was first detected statistically through cross-correlation with galaxy surveys (Smith et al. 2007; Hirata et al. 2008; and more recently Bleem et al. 2012; Sherwin et al. 2012), has now been observed directly in CMB maps by the Atacama Cosmology Telescope (Das et al. 2011) and the South Pole Telescope (van Engelen et al. 2012).

Planck’s frequency coverage, sensitivity and survey area, allow high signal-to-noise measurements of both the CIB and the CMB lensing potential. Accompanying the release of this paper, Planck Collaboration VIII (2014) reports the first measurement and characterization of the CMB lensing potential with the Planck data; this has several times more statistical power than previous measurements, over a large fraction (approximately 70%) of the sky. We will use this measurement of the lensing potential in cross-correlation with measurements of the CIB in the Planck HFI bands to make the first detection of the lensing-infrared background correlation. In addition to our measurement, we discuss the implications for models of the CIB fluctuations. The outline of this paper is as follows. In Sect. 2 we describe the data we will use, followed by a description of our pipeline for correlating the CIB and lensing signals in Sect. 3. Our main result is presented in Sect. 4, with a description of our error budget, consistency tests and an array of systematic tests in Sect. 5. We discuss the implications of the measured correlation for CIB modelling in Sect. 6.

2. Data sets

2.1. Planck maps

Planck (Tauber et al. 2010; Planck Collaboration I 2011) is the third generation space mission to measure the anisotropy of the CMB. It observes the sky with high sensitivity in nine frequency bands covering 30–857 GHz at an angular resolution from 31^′ to 5^′. The Low Frequency Instrument (LFI; Mandolesi et al. 2010; Bersanelli et al. 2010; Mennella et al. 2011) covers the 30, 44, and 70 GHz bands with radiometers that incorporate amplifiers cooled to 20 K. The High Frequency Instrument (HFI; Lamarre et al. 2010; Planck HFI Core Team 2011a) covers the 100, 143, 217, 353, 545, and 857 GHz bands with bolometers cooled to 0.1 K. Polarization is measured in all but the highest two bands (Leahy et al. 2010; Rosset et al. 2010). A combination of radiative cooling and three mechanical coolers produces the temperatures needed for the detectors and optics (Planck Collaboration II 2011). Two data processing centres (DPCs) check and calibrate the data and make maps of the sky (Planck HFI Core Team 2011b; Zacchei et al. 2011). Planck’s sensitivity, angular resolution, and frequency coverage make it a powerful instrument for Galactic and extragalactic astrophysics as well as cosmology. Early astrophysics results are given in Planck Collaboration VIII–XXVI 2011, based on data taken between 13 August 2009 and 7 June 2010. Intermediate astrophysics results are now being presented in a series of papers based on data taken between 13 August 2009 and 27 November 2010. This paper uses data corresponding to the second Planck data release, with data acquired in the period up to 27 November 2010 with improved processing compared to the first release.

We use the Planck HFI temperature maps at all six frequencies, i.e., 100, 143, 217, 353, 545, and 857 GHz. The maps at each frequency were created using almost three full-sky surveys. Here we give an overview of the HFI map-making process, with additional details given in Planck HFI Core Team (2011b); Planck Collaboration II (2014). The data are organized as time-ordered information, hereafter TOI. The attitude of the satellite as a function of time is provided by two star trackers on the spacecraft. The pointing for each bolometer is computed by combining the attitude with the location of the bolometer in the focal plane, as determined by planet observations. The raw bolometer TOI for each channel is first processed to produce cleaned timelines and to set flags that mark bad data (for example data immediately following a cosmic ray strike on the detector). This TOI processing includes: (1) signal demodulation and filtering; (2) deglitching, which flags the strong part of any glitch and subtracts the tails; (3) conversion from instrumental units (volts) to physical units (watts of absorbed power, after a correction for the weak nonlinearity of the response); (4) de-correlation of thermal stage fluctuations; (5) removal of the systematic effects induced by 4 K cooler mechanical vibrations; and (6) deconvolution of the bolometer time response. Focal plane reconstruction and beam shape estimation is made using observations of Mars. The simplest description of the beams, as elliptical Gaussians, leads to full-width at half-maximum (FWHM) values of 9.́65, 7.́25, 4.́99, 4.́82, 4.́68, and 4.́33, as given in Table 4 of Planck Collaboration II (2014). As explained in this paper, the inter-calibration accuracy between channels is better than the absolute calibration. This leads us to adopt conservative absolute calibration uncertainties of 0.64%, 0.53%, 0.69%, 2.53%, 10.0%, and 10.0%, at 100, 143, 217, 353, 545, and 857 GHz, respectively. We convert between emissivities given in  MJy sr^-1 (with the photometric convention νI_ν = const.) and temperatures in  μK, using the measured bandpass filters (see PER and P2013 for details). Note that at 545 and 857 GHz, an extra step is also applied to reduce the zodiacal emission. We are using “zodi-removed” maps throughout this paper.

For the sake of consistency testing (presented in particular in Sect. 5), we will sometimes use temperature maps where only a fraction of the TOI is used to generate the sky map. In particular, throughout this paper we use the terminology “half-ring” (HR) maps to refer to maps made using the first and second half of the stable pointing period, “survey” for individual full-sky survey maps (note that the Survey 3 is incomplete for the particular data release used in the intermediate papers), and “detset” for maps made using two independent sets of detectors per frequency (for details see Planck HFI Core Team 2011b).

We create three masks to exclude regions with bright Galactic and extragalactic foreground emission. The first mask accounts for diffuse Galactic emission as observed in the Planck data. In order to test for the effects of residual Galactic emission on our results we create three different versions of this mask, each with a different masked area, such that 20%, 40%, or 60% of the sky is unmasked. Each version of this mask is created directly from the Planck 353 GHz map, from which we remove the CMB using the 143 GHz channel as a CMB template before smoothing by a Gaussian with FWHM of 5^◦. The map is then thresholded so that the mask has the required sky fraction. Although the Galactic emission is stronger at 857 GHz, we expect the 353 GHz mask to better trace dust emission at the lower frequencies that we use. The mask therefore accounts for Galactic dust and Galactic CO emission, as explained in Planck Collaboration XII (2014). We will be ignoring synchrotron emission, which is important at low frequencies, since its contribution at 100 GHz and at high Galactic latitudes is small, and, as with the dust component, will be uncorrelated with the lensing potential. The second mask covers bright point sources. This mask is created using algorithms tailored to detect point sources in the Planck data and is optimized for each frequency, as detailed in Planck Collaboration VII (2011) and Planck Collaboration (2011). The third mask is designed to remove extended high-latitude Galactic dust emission (“cirrus”), which is traced by external H i data, as we will describe in Sect. 2.2.1. While the first two masks are described in Planck Collaboration XII (2014), the third is specific to our cross-correlation analysis, since it provides a method to reduce the large-scale noise in our measurement, and the 3-point nature of our estimate ensures that it will not introduce a bias (although we test for this in Sect. 5). Ultimately, when we combine the three masks, we obtain an effective sky fraction of 16, 30 and 43% for the 20%, 40%, and 60% Galactic masks, respectively, but we will still refer to them as the 20%, 40%, and 60% masks for the sake of simplicity.

2.2. External data sets

2.2.1. H i maps

We use measurements of 21-cm emission from Galactic neutral hydrogen (H i) as a cirrus monitor. Outside of our Galactic and point source masks we use the H i data to construct a template of the dust emission in regions where the H i column density is low (less than N_HI ≤ 2 × 10²⁰ cm^-2), and we mask regions where it is high, since in these regions the H i and dust emission are not well correlated (Boulanger et al. 1996; Boulanger & Perault 1988, PER). The masking procedure that we use is described in detail in Planck Collaboration XXIV (2011); it consists of subtracting the H i dust template from the Planck temperature map at 857 GHz and calculating the skewness of the residuals in 5 deg² regions. If the skewness is larger than a given value then the region is masked. This is an improvement over the usual cut-off in H i column density. We use the latest release from the Leiden/Argentina/Bonn (LAB) survey (Kalberla et al. 2005), which consists of the Leiden/Dwingeloo Survey (LDS) (Hartmann & Burton 1997) north of − 30^° declination, combined with the Instituto Argentino de Radioastronomia Survey (Arnal et al. 2000; Bajaja et al. 2005) south of − 25^° declination. The angular resolution of the combined map is approximately 0.6^° FWHM. The LAB Survey is the most sensitive Milky Way H i survey to date, with the greatest coverage both spatially and kinematically. We make use of projections of the LAB maps onto N_side = 512HEALPix² maps performed by Land & Slosar (2007) and made available through the LAMBDA website³. The local standard of rest velocity coverage spans the interval − 450 km s^-1 to + 400 km s^-1, at a resolution of 1.3 km s^-1, with an rms brightness-temperature noise of 0.07–0.09 K, and with additional errors due to defects in the correction for stray radiation that are less than 20–40 mK for most of the data.

2.2.2. IRIS/IRAS maps

As a consistency test we will use an additional tracer of the CIB that derives from re-processed IRAS maps at 60 and 100 μm. This new generation of IRAS maps, known as IRIS (Miville-Deschênes & Lagache 2005), benefits from improved zodiacal light subtraction, a calibration and zero level compatible with DIRBE, and an improved de-striping procedure. IRAS made two full-sky maps (HCON-1 and HCON-2), as well as a final map that covers 75% of the sky (HCON-3). The three maps had identical processing that included deglitching, checking of the zero-level stability, visual examination for glitches and artefacts, and zodiacal light removal. The three HCONs were then co-added, taking into account the inhomogeneous sky coverage maps, to generate the average map (HCON-0). Note that the Finkbeiner et al. (1999) maps are also constructed from the IRAS 100 μm data, and as such we will not investigate their cross-correlation properties, since the IRIS map contains the same information. For simplicity we will assume that the effective IRIS beam is uniform across the sky and described by a Gaussian with FHWM of 4.́3.

3. Cross-correlation formalism and implementation

We now describe our statistical formalism and its implementation, with additional technical details given in the appendices. Our analysis consists of cross-correlating a full-sky reconstruction of the CMB lensing potential with a temperature map.

3.1. Reconstructing the CMB lensing potential

The CMB is lensed by the gravitational potential of all matter along the photon trajectory from the last scattering surface to us. The lensed CMB is a remapping of the unlensed CMB with the lensed temperature equal to $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathrm{\Theta }\end{array}\mathrm{\left(}\widehat{{n}\hspace{0.17em}}\mathrm{\right)}\mathrm{=}\mathrm{\Theta }\mathrm{(}\widehat{{n}\hspace{0.17em}}\mathrm{+}\mathrm{\nabla }\mathit{\phi }\mathrm{)}$, where $\mathrm{\Theta }\mathrm{\left(}\widehat{{n}\hspace{0.17em}}\mathrm{\right)}$ is the unlensed CMB temperature and φ is the lensing potential. We use the methodology described in Planck Collaboration VIII (2014) to obtain estimates of the lensing potential in harmonic space, using the standard Okamoto & Hu (2003) quadratic estimator.

Complete details of the lens reconstruction procedure that we use are given in Planck Collaboration VIII (2014), although we review it briefly in point form here. Our estimates of are obtained by the following set of steps:

• 1.

  inverse variance filter the CMB map;

• 2.

  use the filtered CMB map as the input to a quadratic lensing estimator, which is designed to extract the off-diagonal contributions to the CMB covariance matrix induced by lensing;

• 3.

  subtract a “mean-field bias”, which corrects for known non-lensing contributions to the covariance matrix, including instrumental noise inhomogeneity, beam asymmetry, and the Galaxy+point source mask.

The output from this pipeline is an estimate of the lensing potential in harmonic space, , and an associated noise power spectrum, ${\mathit{N}}_{\mathit{L}}^{\mathit{\phi \phi }}$, which we use to weight our cross-correlation estimates. We also produce a set of simulated lens reconstructions, which we use to establish our statistical error bars.

Our nominal lens reconstructions use the 143 GHz channel. However, there is almost equivalent power to measure lensing using the 217 GHz channel. Combining both channels would reduce the noise power spectrum of our lens reconstruction by approximately 25%, compared with using either individually (the improvement is significantly less than 50% because a significant portion of the lens reconstruction noise is due to the finite number of CMB modes that we are able to observe, and is correlated between the two channels). We choose to focus on 143 GHz here because it is significantly less susceptible to CIB contamination. We will use lens reconstructions based on the 100 and 217 GHz data for consistency tests.

3.2. Decreasing the foreground noise

An important source of noise (but, as we will explain below, not bias) in our cross-correlation measurement is Galactic foreground emission. Dust emission is the dominant Galactic component at HFI frequencies above 217 GHz (see Sect. 5.1 for a quantitative discussion). In order to reduce the Galactic dust emission we create a dust template and subtract it from the temperature maps described in Sect. 2.1. At 100 and 143 GHz the CMB signal is significantly brighter than the dust emission outside the Galactic mask. We therefore do not create and subtract a dust template at these frequencies. Note that while we could use other frequency maps to trace the CMB and remove it, to quantify the non-negligible amount of CIB that would be removed this way is not easy, given the uncertainties in the cross-frequency CIB correlation structure.

We rely on the well documented (but complex) correlation between Galactic H i and dust (e.g., Boulanger & Perault 1988; Boulanger et al. 1996; Lagache et al. 1998, PER) to reduce the contamination by subtracting the H i-correlated dust component. As was performed in PER, we split the H i map into two velocity components: a low-velocity local component (LC) typical of high-latitude H i emission, and a component of intermediate-velocity clouds (IVC). We found that the inclusion of a high-velocity component makes a negligible difference to the dust-cleaned map. Unlike the dedicated high-resolution H i observations used in PER and P2013 that only partially cover the sky, here we use the full-sky, low resolution LAB survey introduced in Sect. 2.2.1 as our H i tracer. Although the resolution of this survey is lower than the Planck resolution, it allows us to perform dust cleaning on large scales, where our cross-correlation measurement has high signal-to-noise ratio. The emissivity of the dust varies across the sky, and so the correlation between the dust and H i emission is expected to vary. To account for this we divide the sky into regions within which we assume that the dust-H i correlation is constant. For the sake of convenience, we use regions of approximate size 13 (52) deg² defined by the Healpix pixels at resolution N_side = 16 (8) that are outside the Galactic mask. We test that our conclusions do not depend on this resolution.

The details of our procedure are as follows. We subtract the 143 GHz Planck temperature map from each of the 217–857 GHz temperature maps to remove the CMB signal (this CMB subtraction is only done for the purposes of creating the dust template). We upgrade each of the N_side = 512 LAB maps compiled in Land & Slosar (2007) to the Planck map resolution of N_side = 2048. Within each region we then simultaneously fit for the amplitude of each H i velocity component in the CMB-subtracted maps, and use the two coefficients per region to assemble a full-sky (minus the mask) dust template for each of the 217–857 GHz channels. We smooth each template with a Gaussian of FWHM 10^′ to remove the discontinuity at the patch boundaries, and then subtract the template from the original (CMB-unsubtracted) Planck maps.

We note that the accuracy of this procedure would be difficult to evaluate for all possible uses of the map, in other words whether it constitutes a robust component separation method remains to be demonstrated. Our approach and results are however consistent with the dust-focused analysis detailed in Planck Collaboration V (2014). However, in the case of our cross-correlation analysis the dust-removal requirements are less severe, since the dust emission only contributes to our measurement as a noise source. We will describe later in Sect. 5.2 the effect on the cross-spectrum of removing this emission, and will place limits on the residual Galactic contamination in Sect. 5.3.5.

3.3. Measuring cross-correlations

To estimate the cross-correlation between the lensing potential and a tracer t, we calculate $\begin{array}{ccc}\mathit{Ĉ}\begin{array}{c}\mathit{t\phi }\\ \mathit{\ell }\end{array}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}\sum _{\mathit{m}}\mathit{t̂}\mathit{\ell m}\mathit{\phi ̂}\begin{array}{c}\mathrm{\ast }\\ \mathit{\ell m}\end{array}\mathit{.}& & \end{array}$(1)Since the CIB has an approximately ℓ^-1 dependence and the lensing potential has an ℓ^-2 dependence, we multiply the cross-spectra by ℓ³, and then bin it in 15 linearly spaced bins between ℓ = 100 and 2000. As we will discuss in Sect. 5, modes with ℓ < 100 are not considered, due to possible lens mean-field systematic effects, and modes with ℓ > 2000 are also removed, due to possible extragalactic foreground contamination. We have tested that our results are robust to an increase or decrease in the number of ℓ bins.

We expect the error bars to be correlated across bins to some extent, due to pseudo-C_ℓ mixing induced by the mask, and between frequencies, because the lens reconstruction noise is common. In addition, any foregrounds that are present in multiple channels will introduce correlated noise. The foreground mask will also induce a coupling between different modes of the unmasked map. This extra coupling can be calculated explicitly using the mixing matrix formalism introduced in Hivon et al. (2002). Using this formalism and our best-fit models we have evaluated the correlation between different bins of the cross-correlation signal for our nominal binning scheme. We find that the mask-induced correlation is less than 2% across all bins at all frequencies. We will thus neglect it in our analysis. For this reason, and based on the results we obtain from simulations, we do not attempt to “deconvolve” the mask from the cross-spectrum (see e.g., Hivon et al. 2002) and instead correct for the power lost through masking the maps by a simple sky fraction, f_sky, ignoring the mode coupling.

As will be discussed later in Sect. 6.1, when we fit models to the cross-spectrum we will assume that the noise correlation between bins can be neglected and that the band-powers are flat.

3.4. Simulating cross-correlations

In order to validate our measurement pipeline, and to confirm that our estimate of the cross-spectrum is unbiased, we create simulated maps of the lensed CMB and CIB that have the expected statistical properties.

Using the Planck only favoured ΛCDM cosmology as described in Planck Collaboration VII (2014) we generate a theoretical prediction of the lensing potential spectrum using CAMB (Lewis et al. 2000). We use this to generate 300 maps of φ that are used to lens 300 CMB realizations using the approach described in Planck Collaboration VIII (2014). We then use the PER best-fit CIB model to generate CIB auto- and CIB-φ cross-spectra, from which we create CIB simulations that are correctly correlated with φ in each HFI band. The PER model that we use describes the CIB clustering at HFI frequencies through a halo approach, and simultaneously reproduces known number count and luminosity function measurements. At each frequency we add a lensed CMB realization to each of the CIB simulations and then smooth the maps using a symmetric beam with the same FWHM as the beam described in Sect. 2.1. Once this set of realizations has been generated we apply the reconstruction procedure described above to produce an estimate of the lensing potential map, and then calculate the cross-power spectrum using our measurement pipeline.

These simulations will miss some complexities inherent in the Planck mission. They do not take into account inhomogeneous and correlated noise, and we do not simulate asymmetric beam effects. In addition, we do not simulate any foreground components, and we instead take a different approach to determine their contribution. While simplistic, we believe that our simulations are good enough for the purposes of this particular measurement. In Sect. 5 we will discuss possible limitations, as well as how we test for systematic effects that are not included in the simulations.

We use the simulated maps to check that our pipeline correctly recovers the cross-spectrum that was used to generate the simulations. For the ℓ bins used in our analysis, we find that the recovered spectrum is unbiased (to within the precision achievable with 300 simulations), and with a noise level consistent with expectations. The noise in the recovered spectrum is discussed in Sect. 5.1.

4. A strong signal using Planck HFI data

We now describe the result of applying our pipeline to our nominal data set, i.e., the lens reconstruction at 143 GHz and the foreground-reduced Planck HFI temperature maps with a 40% Galactic mask (which when combined with the point source mask and H i mask leaves 30.4% of the sky unmasked). The results are presented in Fig. 3. The error bars correspond to the naive scatter measured within each bin. The thin black line corresponds to the expected CIB-lensing correlation predicted using the PER CIB model (the HOD parameters of the PER 217 GHz best-fit model were used at 100 and 143 GHz, since no CIB clustering measurement at these frequencies is available). As can be seen from these plots, the noise is strongly correlated across frequencies, especially at the lowest frequencies where the CMB dominates the error budget. A detailed analysis of the uncertainties and potential systematic errors attached to this result is presented in Sect. 5.

As clearly visible in Fig. 3, a strong signal is detected. To set a reference point and naively quantify its statistical significance when taken at face value, we define a detection significance as follows. We count the number of standard deviations as the quadrature sum of the significance in the different multipole bins: $\begin{array}{ccc}{\mathit{s}}_{\mathit{\nu }}\mathrm{=}\sqrt{\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{15}}{\left(\frac{{\mathit{C}}_{\mathit{i}}^{\mathrm{T}\mathit{\phi }}}{\mathrm{\Delta }{\mathit{C}}_{\mathit{i}}^{\mathrm{T}\mathit{\phi }}}\right)}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(2)For our nominal parameters this gives us 3.6σ, 4.3σ, 8.3σ, 31σ, 42σ, and 32σ, at 100, 143, 217, 343, 545, and 857 GHz, respectively. Note that these numbers include an additional 20% contribution to the statistical error to account for mode correlations (which we discuss in Sect. 5.1), but do not include systematic errors or our point source correction. As a comparison, in each panel we plot the correlation between the lens reconstruction at 143 GHz and the 143 GHz map in grey. This shows the frequency scaling of our measured signal and also the strength of the signal, as compared to possible intra-frequency systematic effects. This will be studied in depth in Sect. 5.

This first pass on our raw data demonstrates a strong detection that is in good agreement with the expected CIB-lensing signal. To get a better visualization of this detection, we show in Fig. 4 the real-space correlation between the observed temperature and the lens deflection angles. This figure allows us to visualize the correlation between the CIB and the CMB lensing deflection angles for the first time. These images were generated using the following stacking technique. We first mask the 545 and 857 GHz temperature maps with our combined mask (that includes the 20% Galaxy mask), and identify 20 000 local maxima and minima in these maps. We also select 20 000 random locations outside the masked region to use in a null test. We then bandpass-filter the lens map between ℓ = 400 and 600 to remove scales larger than our stacked map as well as small-scale noise. We stack a 1 deg² region around each point in both the filtered temperature map and lensing potential map, to generate stacked CIB and stacked lensing potential images. We take the gradient of the stacked lensing potential to calculate the deflection angles, which we display in Fig. 4 as arrows. The result of the stacking over the maxima, minima and random points is displayed from left to right in Fig. 4. The strong correlation seen already in the cross-power spectrum is clearly visible in both the 545 and 857 GHz extrema, while the stacking on random locations leads to a lensing signal consistent with noise. From simulations, we expect a small offset (≃1′) in the deflection field. This offset was corrected for in this plot. We have verified in simulations that this is due to noise in the stacked lensing potential map that creates a random miscentring, even after stacking 20 000 points. This effect is not present when we consider noise-free simulations. It wouldt thus disappear were we to increase this number, but it is obviously not possible given the size of our patch (1 deg²). As expected, we see that the temperature maxima of the CIB, which contain a larger than average number of galaxies, deflect light inward, i.e., they correspond to gravitational potential wells, while temperature minima trace regions with fewer galaxies and deflect light outward, i.e., they correspond to gravitational potential hills.

5. Statistical and systematic error budget

The first pass of our pipeline suggests a strong correlation of the CIB with the CMB lensing potential. We now turn to investigate the strength and the origin of this signal. We will first discuss the different contributions to the statistical error budget in Sect. 5.1, and then possible systematic effects in Sect. 5.2. Although the most straightforward interpretation of the signal is that it arises from dusty star-forming galaxies tracing the large-scale mass distribution, in Sect. 5.3 we consider other potential astrophysical origins for the observed correlation.

5.1. Statistical error budget

In this section we discuss any noise contribution that does not lead to a bias in our measurement. The prescription adopted throughout this paper is to obtain the error estimates from the naive Gaussian analytical error bars calculated using the measured auto-spectra of the CIB and lensing potential. We find that these errors are approximately equal to 1.2 times the naive scatter within an ℓ bin, and we will sometimes use this prescription for convenience, where appropriate (as will be stated in the text). This is justified in Appendix A where we consider six different methods of quantifying the statistical errors using both simulations and data. The Gaussian analytical errors, $\mathrm{\Delta }\mathit{Ĉ}\begin{array}{c}\mathrm{T}\mathit{\phi }\\ \mathit{\ell }\end{array}$, are calculated using the naive prescription $\begin{array}{ccc}{\mathit{f}}_{\mathrm{sky}}\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}\mathrm{\Delta }\mathit{\ell }\hspace{0.17em}{\mathrm{(}\mathrm{\Delta }\mathit{Ĉ}{\begin{array}{c}\mathrm{T}\mathit{\phi }\\ \mathit{\ell }\end{array}}^{\mathrm{)}}}^{\mathrm{2}}& \mathrm{=}\mathit{Ĉ}\begin{array}{c}\mathrm{TT}\\ \mathit{\ell }\end{array}\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{\ell }\end{array}\mathrm{+}{\mathrm{(}{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{T}\mathit{\phi }}}^{\mathrm{)}}}^{\mathrm{2}}\mathit{,}& \end{array}$(3)where, as before, f_sky is the fraction of the sky that is unmasked, Δℓ = 126 for our 15 linear bins between ℓ = 100 and ℓ = 2000, $\mathit{Ĉ}\begin{array}{c}\mathrm{TT}\\ \mathit{\ell }\end{array}$ and $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{\ell }\end{array}$ are the spectra measured using the data, and ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{T}\mathit{\phi }}$ is the model cross spectrum. This last term provides a negligible contribution, due to the large noise bias on $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{\ell }\end{array}$, as we now describe.

The statistical error has two sources, instrumental and astrophysical. The measured auto-spectra in Eq. (3) contain a signal and noise contribution: $\mathit{Ĉ}\begin{array}{c}\mathit{XX}\\ \mathit{\ell }\end{array}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathit{XX}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathit{XX,N}}$. It is informative to split the right hand side of Eq. (3) into four pieces: $\begin{array}{ccc}[{\mathit{C}}_{\mathit{\ell }}^{\mathit{\phi \phi }}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}}& \mathrm{+}& ({\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}\mathit{\phi }}{)}^{\mathrm{2}}]\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathit{\phi \phi ,}\mathrm{N}}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}\mathit{,}\mathrm{N}}\\ & \mathrm{+}& {\mathit{C}}_{\mathit{\ell }}^{\mathit{\phi \phi }}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}\mathit{,}\mathrm{N}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathit{\phi \phi ,}\mathrm{N}}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}}\mathit{.}\end{array}$(4)Here the first term is a signal-only piece, the second is a noise-only piece, and the remaining two terms are mixed signal and noise pieces. To discuss the relative importance of these terms, we will use for the signal terms the model spectra, and for the noise terms we subtract the model spectra from the measured spectra: $\mathit{Ĉ}\begin{array}{c}\mathit{XX,N}\\ \mathit{\ell }\end{array}\mathrm{=}\mathit{Ĉ}\begin{array}{c}\mathit{XX}\\ \mathit{\ell }\end{array}\mathrm{-}{\mathit{C}}_{\mathit{\ell }}^{\mathit{XX}}$. With this definition, the noise will contain the instrumental contribution, as well as other astrophysical signals including the CMB, which we do not remove from our data for reasons previously mentioned. We show the different terms in Fig. 5. Up to 353 GHz the measured temperature spectrum, $\mathit{Ĉ}\begin{array}{c}\mathrm{TT}\\ \mathit{\ell }\end{array}$, is dominated by the CMB at low ℓ and the instrumental noise at high ℓ. At higher frequencies Galactic emission dominates at low ℓ and the CIB at high ℓ. For all frequencies up to 353 GHz the dominant contribution to the errors in our signal comes from the noise-only term (in blue), which is proportional to the temperature noise spectrum. At 353 GHz and above the mixed signal-noise term ${\mathit{C}}_{\mathit{\ell }}^{\mathit{\phi \phi ,}\mathrm{N}}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}}$ (orange in Fig. 5) becomes important and is the largest contribution at 545 and 857 GHz at high ℓ.

5.2. Instrumental and observational systematic effects

A number of systematic errors affect the Planck HFI analysis and we briefly discuss some of them here. A more complete discussion can be found in Planck HFI Core Team (2011b). We will illustrate how the very nature of our measurement, a 3-point function, makes it particularly robust to many systematic effects, and we will check for their signatures using null tests. For example, there is no noise bias in the 3-point measurement, and many effects that can lead to biases in the auto-spectrum of φ do not affect us.

5.2.1. Potential sources of systematic error

We begin by describing our knowledge of known systematic effects, before discussing others that could bias our result. To account for an error in the calibration of the temperature maps, we simply add in quadrature a calibration uncertainty to our error bars. In Sect. 5.2.2 we use null tests to check that these errors are consistent with the data. In addition we use the null tests to search for evidence that the calibration has changed between surveys, for example due to gain drifts. We account for beam errors in a conservative manner by using a constant uncertainty at each frequency equal to the maximum error in the beam multipoles, B_ℓ, at any ℓ (see P2013 for details). The B_ℓ uncertainties are larger at high ℓ, but still remains small. For ℓ = 1500–2000 they are of order 1% at 100 GHz and below 0.5% at 143 and above.

The calibration error is therefore larger than the beam error at all ℓs between 217 and 857 GHz, but smaller at high ℓ in the 100 and 143 GHz channels. We add the beam error in quadrature in a multipole- and frequency-dependent manner. As discussed in Planck Collaboration VIII (2014), uncertainties in the beam transfer (as well as the fiducial CMB power spectrum ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{TT}}$) propagate directly to a normalization uncertainty in the lens reconstruction. Based on the beam eigenmodes of Planck Collaboration III (2014), it is estimated in Planck Collaboration VIII (2014) that beam uncertainty leads to an effective normalization uncertainty of approximately 0.2% at 143 and 217 GHz, and 0.8% at 100 GHz. To be conservative, on top of the calibration and beam error we will add in quadrature a 2% uncertainty on the overall lens normalization.

CMB lens reconstruction recovers modes of the lensing potential through their anisotropic distorting effect on small-scale hot and cold spots in the CMB. The quadratic estimator, which we use to reconstruct the lensing potential is optimized to measure the lensing-induced statistical anisotropy in CMB maps. However, other sources of statistical anisotropy, such as the sky mask, inhomogeneous noise, and beam asymmetries, produce signals that can potentially overlap with lensing. These introduce a “mean-field” bias, which we estimate using Monte Carlo simulations and then subtract from our lensing estimates. Innaccuracies in the simulation procedure will lead to errors in this correction, particularly if the correction is large. The mean-field introduced by the application of a Galaxy and point-source mask, for example, which can be several orders of magnitude larger than the lensing signal at ℓ < 100; this is discussed further in Appendix B of Planck Collaboration VIII (2014). The mask mean-field is a particular concern for us because it has the same phases as the harmonic transform of the mask. If our masked CIB maps have a non-zero monopole, for example, this will correlate strongly with any error in the mask mean-field correction. For this reason we do not use any data below ℓ = 100 in our analysis.

To summarize, we do not expect these known systematic effects to be present at a significant level. Nevertheless, we still perform a set of consistency tests that would be sensitive to them or other unexpected effects.

5.2.2. Null tests

The Planck scanning strategy, its multiple frequency bands, and its numerous detectors per frequency, offer many opportunities to test the robustness of our signal (see Sect. 2.1). We focus on such tests in this section. Our aim is to reveal any systematic effects that could lead to a spurious correlation. For all of the tests presented, we will quote a χ² value, as well as the number of degrees of freedom (N_d.o.f.) as a measure of the consistency with the expected null result. Throughout this section, black error bars in plots will correspond to the measured scatter within an ℓ bin, multiplied by 1.2 as was justified in Sect. 5.1 and Appendix A, and will also include a CIB calibration error and a beam error, while the coloured boxes correspond to the analytical errors of the corresponding signal (i.e., not the difference corresponding to the null test). Plotting these two error bars illustrates how important any deviation could be to our signal. Throughout this section, we will illustrate our findings with the 545 GHz correlation, since it is our prime band for this measurement, but our conclusions hold at other frequencies.

The first test we conduct is to take the temperature difference between the two half-ring (HR) maps to cancel any signal contribution, and therefore investigate the consistency of our measurements with our statistical errors over all time scales. We null the temperature maps and correlate with our nominal lensing map. The results are shown in the left panel of Fig. 6. No significant deviation is found.

The second test uses multiple detectors at a given frequency that occupy different parts of the focal plane. These detector sets are used to construct the “detset” maps that were described in Sect. 2.1. The two “detset” maps are subtracted and then correlated with our nominal lens reconstruction. This test is particularly sensitive to long-term noise properties or gain variations, as we do not expect these to be correlated from detector to detector. Since this detector division breaks the focal plane symmetry, it is also a good check for beam asymmetry effects. Here again, we do not find any significant deviation, as illustrated in the middle panel of Fig. 6.

The third test we conduct makes use of the redundant sky coverage, using multiple surveys to cancel the signal. As above, we null the temperature signal and correlate with the nominal lens reconstruction. This test is particularly sensitive to any long-term (i.e., month timescale) drifts that could affect our measurement. It is also a good test for any beam asymmetry effects, as individual pixels are observed with a different set of orientations in each survey. Since only the first two surveys are complete for this particular data release, we only use the two full survey maps, to avoid complications with the partially completed third survey. The results are displayed in the right panel of Fig. 6. We see a significant deviation from null. This particular failure can probably be explained by apparent gain drifts due to nonlinearity in the analogue-to-digital conversion (Planck Collaboration II 2014; Planck Collaboration IV 2014), not yet corrected at this frequency. Note, however, that the predicted variation is about 1%, while the deviation from null would call for a variation of 1.5–2%. But in any case, its amplitude is too small to significantly affect our measurement.

To conclude, this first set of stringent consistency tests have shown that there is no obvious contamination of our measurements due to instrumental effects, except for the effect of known long-term gain evolution. Apart from this failure, which we take into account by increasing calibration uncertainties, the reasonable χ²/N_d.o.f. obtained gives us confidence in our statistical noise evaluation. Although we measure the noise directly from the data, this success was not guaranteed.

5.3. Astrophysical contamination

We now turn to possible astrophysical biases to our measurement. We will discuss successively known astrophysical contaminants that have either Galactic or extragalactic origin. Once again, besides our knowledge of these signals, we will rely heavily on consistency tests made possible by having multiple full sky frequency maps.

5.3.1. Galactic foregrounds

Galactic foregrounds produce two possible effects on our measurement. The first is the introduction of an extra source of noise. The second is that contamination of the lensing reconstruction by any Galactic signal, e.g., synchrotron, free-free or dust, which could then correlate with foreground emission present in the temperature maps, remains a source of bias that has to be investigated. We will show that this bias is small. To do so, we take three approaches. We first investigate several Galactic masks, then perform the lensing reconstruction at various frequencies, and finally investigate the effect of a dust-cleaning procedure.

First, we consider two additional masks, either more aggressive or more conservative than our fiducial one. Both were introduced in Sect. 2.1. The first one leaves approximately 20% of the sky unmasked, while the second one leaves approximately 60% of the sky. Given the strong dependence of Galactic foregrounds on Galactic latitude, any Galactic contamination should vary strongly when we switch between masks. Comparing the measurements using these masks with our fiducial 40% mask in the left and centre panels of Fig. 7, we do not see any substantial deviation from our fiducial measurements. This excludes strong Galactic contamination of our results.

Second, we perform the lens reconstruction at 100 and 217 GHz, different from the fiducial frequency of 143 GHz, and compare their correlation with the temperature maps. Given the strong dependence of the Galactic emission with frequency, T ∝ ν^-3 for synchrotron and T ∝ ν² for dust in this frequency range, any contamination of our signal would have a strong frequency dependence. The comparison with the 100 GHz (217 GHz) reconstruction is presented in the left (centre) panel of Fig. 8. The right panel shows the difference of the cross-spectra calculated using the 143 GHz reconstruction with a more aggressive Galaxy mask (20% instead of 40%, to reduce possible Galactic contaminants in the reconstruction), and the nominal reconstruction. The three differences are consistent with null, as demonstrated by the quoted χ² and N_d.o.f..

Third, we investigate more specifically how cirrus, the dominant Galactic contaminant for our higher frequency channels, affects our measurements. We rely on the dust cleaning procedure detailed in Sect. 3.2 that aims to remove the H i-correlated dust component. This procedure leads to a decrease in the variance measured outside the mask of 22, 65, 73, and 73% in the 217, 353, 545, and 857 GHz maps, respectively. This frequency dependence is expected, given the dust scaling. However, in Fig. 7, where we show the differences between the cleaned and non-cleaned cross-spectra, we observe that the large-scale H i cleaning, even though it has a substantial impact on the power within our map, only makes a small change at low ℓ in the cross-spectrum, as well as reducing the noise at all multipoles. If we quantify the effect of our “local” H i cleaning on the detection significance level computed using only statistical errors, we find that the significance is increased by 4, 4, 28, and 36% at 217, 353, 545, and 857 GHz, respectively. Also, not surprisingly, we observe that for frequencies up to 353 GHz, where the statistical errors are dominated by the CMB, the H i cleaning has almost no effect on the cross-spectra. From the three studies in this section we conclude that there are no obvious signs of Galactic foreground contamination in our cross-correlation.

5.3.2. Point source contamination

We now discuss another well-known potential source of contamination, namely the contribution of unresolved point sources visible either through their radio or dust emission. Our concern is that a correlation between a spurious lens reconstruction caused by unresolved point sources can correlate with sources in the temperature map. Although in Sect. 5.3.1 our null test using lens reconstructions at different frequencies suggests that unresolved point sources are not an obvious contaminant, we will now perform a more detailed test designed specifically to search for point source contamination. Following Smith et al. (2007) and Osborne et al. (in prep.), we will construct a point source estimator designed to be more sensitive than the lensing estimator to point source contamination. Our focus here will be on possible contamination from the point source shot-noise bispectrum. In Sect. 5.3.5 we will discuss contamination from a scale dependent bispectrum.

Our (unnormalized) quadratic estimator, which is designed to detect point source contributions, is given by $\begin{array}{ccc}{\mathrm{[}\mathrm{\Theta ̅}\mathrm{143}\mathrm{(}\mathit{n̂}{\mathrm{)}}^{\mathrm{]}}}_{\mathit{LM}}^{\mathrm{2}}\mathrm{\equiv }\sum _{\mathit{LM}}{\mathit{Y}}_{\mathit{LM}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\hspace{0.17em}{\mathrm{[}\mathrm{\Theta ̅}\mathrm{143}\mathrm{(}\mathit{n̂}{\mathrm{)}}^{\mathrm{]}}}^{\mathrm{2}}\mathit{,}& & \end{array}$(5)where is the inverse-variance filtered sky map. This estimator is simply the square of the inverse-variance filtered sky map, which is a more sensitive probe of point sources than the standard lensing estimator.

In Fig. 10 we plot the cross-spectrum of $[\mathrm{\Theta ̅}{}^{\mathrm{143}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}{]}_{\mathit{LM}}^{\mathrm{2}}$ measured at 143 GHz and $\mathrm{\Theta ̅}\begin{array}{c}\mathit{\nu }\\ \mathit{LM}\end{array}$ for the full set of HFI channels. This cross-spectrum is probing the same point source contributions that could bias our estimates of ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{T}\mathit{\phi }}$, although with a greater signal-to-noise ratio.

There is one complication here, which is that just as lens reconstruction may be biased by point source contributions, the point source estimator is correspondingly biased by lensing. The bias to the plotted cross-spectra is given by $\begin{array}{ccc}& & {⟨{\mathrm{[}\mathrm{\Theta ̅}\mathrm{143}\mathrm{(}\mathit{n̂}{\mathrm{)}}^{\mathrm{]}}}_{\mathit{LM}}^{\mathrm{2}}\mathrm{\Theta ̅}\begin{array}{c}\mathit{\nu }\mathrm{\ast }\\ \mathit{LM}\end{array}⟩}_{\mathit{\phi }}\mathrm{=}\frac{{\mathit{C}}_{\mathit{L}}^{\mathit{\nu \phi }}}{{\mathit{C}}_{\mathit{L}}^{\mathit{\nu \nu }}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\sum _{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\frac{{\mathit{G}}_{\mathit{L}{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{-}\mathit{M}{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}}}{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{tot}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{tot}}}\\ & & \mathrm{\times }\mathrm{[}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{{\mathit{m}}_{\mathrm{1}}}{\mathit{I}}_{{\mathit{\ell }}_{\mathrm{2}}\mathit{L}{\mathit{\ell }}_{\mathrm{1}}}^{{\mathit{m}}_{\mathrm{2}}\mathrm{-}\mathit{M}\mathrm{-}{\mathit{m}}_{\mathrm{1}}}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}\end{array}\mathrm{+}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{{\mathit{m}}_{\mathrm{2}}}{\mathit{I}}_{{\mathit{\ell }}_{\mathrm{1}}\mathit{L}{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{m}}_{\mathrm{1}}\mathrm{-}\mathit{M}\mathrm{-}{\mathit{m}}_{\mathrm{2}}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}\end{array}}^{\mathrm{]}}\mathit{,}\end{array}$(6)

where ${\mathit{G}}_{{\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{\int }}\mathrm{d}\mathit{n̂}\hspace{0.17em}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}$, $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$ is the unlensed CMB spectrum, ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{tot}}$ is the spectrum of , ${\mathit{I}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}$ is as defined in Okamoto & Hu (2003), and we have used the fact that for our inverse-variance filtering, $\mathrm{\Theta ̅}{}_{\mathit{\ell m}}\mathrm{\approx }{\mathrm{\Theta }}_{\mathit{\ell m}}\mathit{/}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{tot}}$. We have calculated this contribution using our measured ${\mathit{C}}_{\mathit{L}}^{\mathrm{T}}\mathit{\phi }$ and subtracted it from the data points of Fig. 10.

We can consider the effect of shot noise on this cross-spectrum. With the shot-noise bispectrum defined by $\begin{array}{ccc}{\mathrm{⟨}{\mathrm{\Theta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}{\mathrm{\Theta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}{{\mathrm{\Theta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}^{\mathit{\nu }}}^{\mathrm{⟩}}}_{{\mathit{S}}^{\mathrm{3}}}\mathrm{=}{\mathit{G}}_{{\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{⟨}{\mathit{S}}^{\mathrm{3}}\mathrm{⟩}& & \end{array}$(7)the bias to the plotted cross-spectrum is given as $\begin{array}{ccc}{⟨{\mathrm{[}\mathrm{\Theta ̅}\mathrm{143}\mathrm{(}\mathit{n̂}{\mathrm{)}}^{\mathrm{]}}}_{\mathit{LM}}^{\mathrm{2}}\mathrm{\Theta ̅}\begin{array}{c}\mathit{\nu }\\ \mathit{LM}\end{array}⟩}_{{\mathit{S}}^{\mathrm{3}}}\mathrm{=}\mathrm{⟨}{\mathit{S}}^{\mathrm{3}}\mathrm{⟩}\frac{\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{M}}}{{\mathit{C}}_{\mathit{L}}^{\mathit{\nu \nu }}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\sum _{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\frac{{\mathit{G}}_{\mathit{L}{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{-}\mathit{M}{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}}{\mathit{G}}_{\mathit{L}{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{M}{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}}}{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{tot}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{tot}}}\mathrm{·}& & \\ & & \end{array}$(8)This bias is plotted for best-fit values of ⟨ S³ ⟩ as the black lines in Fig. 10. To minimize systematic effects that might bias the ⟨S³⟩ estimator, we have estimated S³ from the spectra of Fig. 10 for multipoles between ℓ = 500 and 2000. The fitted ⟨S³⟩ amplitudes are given in Table 1.

These amplitudes match our expectations (for example see Planck Collaboration XIII 2011). We observe a decrease in the amplitude of the point source contribution going from 100 to 217 GHz, which corresponds to a dominant contribution from radio point sources. We do not see any evidence of a dusty galaxy contribution to the shot-noise bias. These conclusions have been verified using less restrictive point source masks that cover fewer sources.

With estimates of S³ in hand, we estimate a corresponding bias to ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{T}\mathit{\phi }}$, given by $\begin{array}{ccc}{⟨\mathit{\phi ̂}\mathit{LM}\mathrm{\Theta ̅}\begin{array}{c}\mathit{\nu }\mathrm{\ast }\\ \mathit{LM}\end{array}⟩}_{{\mathit{S}}^{\mathrm{3}}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{M}}\frac{\mathrm{⟨}{\mathit{S}}^{\mathrm{3}}\mathrm{⟩}}{{\mathit{C}}_{\mathit{L}}^{\mathit{\nu \nu }}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\sum _{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\frac{{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}\mathit{M}}}{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{tot}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{tot}}}{\mathit{G}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}\mathrm{-}\mathit{M}}\mathit{,}& & \end{array}$(9)where ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}\mathit{M}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{{\mathit{m}}_{\mathrm{2}}}\mathrm{(}{\mathit{I}}_{{\mathit{\ell }}_{\mathrm{2}}\mathit{L}{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{-}{\mathit{m}}_{\mathrm{2}}\mathit{M}{\mathit{m}}_{\mathrm{1}}}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}\end{array}\mathrm{+}{\mathit{I}}_{{\mathit{\ell }}_{\mathrm{1}}\mathit{L}{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{m}}_{\mathrm{1}}\mathit{M}\mathrm{-}{\mathit{m}}_{\mathrm{2}}}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}\end{array}\mathrm{)}\mathit{/}\mathrm{2}{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{\phi \phi }}$, with ${\mathrm{ℛ}}_{\mathit{L}}^{\mathit{\phi \phi }}$ defined in Planck Collaboration VIII (2014). We show this contribution later in Fig. 12 as the dotted line. While non-zero, we see that the point source shot-noise contribution is always negligible in the ℓ range we consider, except at lower frequencies, where the radio point sources are important (but still not strong enough to lead to any clear signal in the cross-spectra).

5.3.3. SZ contamination

A fraction of CMB photons travelling from the surface of recombination are scattered by hot electrons in galaxy clusters. In the most massive clusters approximately 1% of CMB photons passing through them get scattered. On average, their energy will be increased, which leads to a measurable spectral distortion. This is the so-called thermal Sunyaev-Zeldovich (SZ) effect (Sunyaev & Zeldovich 1970). At the location of a galaxy cluster the CMB appears colder at frequencies below about 220 GHz and hotter at higher frequencies, with a temperature change proportional to the cluster optical depth for Compton scattering and to the electron temperature. Since hot electrons in clusters also trace the large-scale matter potential that is traced by CMB lensing, we expect an SZ-induced contamination in our measurement at some level. We will show below that the level of contamination is negligible. In these calculations we ignore the small relativistic corrections to the thermal SZ spectrum (e.g., Nozawa et al. 2000). We also ignore the kinetic SZ signal coming from the bulk motion of hot electrons in clusters, since it is subdominant to the thermal signal (Sunyaev & Zeldovich 1980; Reichardt et al. 2012; Hand et al. 2012).

The frequency dependence of the SZ signal in our map depends on the detector bandpasses and is $\begin{array}{ccc}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\frac{\mathrm{\int }\mathrm{d}\mathit{\nu h}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{g}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}{\mathrm{\int }\mathrm{d}\mathit{\nu h}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}\mathit{,}& & \end{array}$(10)where h(ν) is the detector bandpass and g(ν) is the SZ frequency dependence, which in the non-relativistic limit is g(ν) = x(e^x + 1) / (e^x−1) − 4, with x = hν/k_BT_CMB. The effect of the bandpass only makes a large difference at 217 GHz near the null of the SZ signal. The thermal SZ signal affects our measurement in two ways. First, since the SZ emission in our maps is not a Gaussian random field (e.g., Wilson et al. 2012), it introduces a spurious signal into our lens reconstruction that will correlate with the SZ signal in our CIB map. As shown in Osborne et al. (in prep.), this is well approximated by a Poisson noise term and is thus already addressed by our treatment of point sources in Sect. 5.3.2. The spurious lensing signal will also correlate with other components in our map, such as the CIB. However, we ignore these terms, since they will be smaller than those that correlate directly with the SZ emission. Additionally, a contribution comes from SZ emission in our CIB map that correlates with the lensing potential itself. The latter is the dominant term and we discuss it in this section.

To measure a contribution from the SZ-lensing correlation we attempt to separate the SZ and CIB emission based on their differing spectral shapes. We consider all frequencies from 100 to 857 GHz, but we will illustrate this procedure by considering only two ℓ bands: ℓ = 300–450; and 1200–1450. The first is well inside the linear regime, while the second receives a more important nonlinear contribution. However, we have checked that if we consider different ℓ bins we obtain similar conclusions. We model the signal within each ℓ band as s_ℓ(ν) = a_1,ℓc(ν) + a_2,ℓf(ν), where c(ν) and f(ν) are, respectively, the CIB frequency dependence (as proposed in Fixsen et al. 1998; or Gispert et al. 2000) and the SZ frequency dependence obtained from Eq. (10). For each ℓ band, we will solve for a_1,ℓ and a_2,ℓ, minimizing the associated χ² while forcing both amplitudes to be positive. As an approximation to the error in each multipole band we calculate the scatter of the signal within the band and multiply it by 1.2, as discussed in Sect. 5.1.

In Fig. 9 we show the measured frequency spectrum within each ℓ band, along with the best-fit SZ-lensing and CIB-lensing spectra. For the CIB-only fit with the Gispert et al. (2000) frequency dependence we find a relatively poor fit in the lowest ℓ bin, χ² (N_d.o.f.) = 15.5 (5), but an improved fit in the higher ℓ bin, χ²(N_d.o.f.) = 4.15 (5). Including the SZ component gives an improvement in the χ² of 0.52 and 1.34 in the low- and high-ℓ bins for one extra degree of freedom. When we use the Fixsen et al. (1998) frequency dependence we find an improved fit, with χ²(N_d.o.f.) = 2.25 (5) and 5.49 (5) in the low- and high-ℓ bins, respectively. Overall, the improvement in the χ²/N_d.o.f. when including the SZ component does not justify inclusion of the SZ component in the model, with the poor fit driven by the lowest frequency bands, where the CIB scaling is rather unconstrained. In fact, our measurements might constitute the first constraints to date on this scaling. From these results we conclude that including the SZ-lensing correlation in our data does not improve the fit in the ℓ range of interest to us and thus we do not consider it necessary to make a correction.

As an extra validation of this result, we now verify its consistency with current models of the CIB and SZ emission. For this purpose, we use the calculation of the correlation from Osborne et al. (in prep.), based on Babich & Pierpaoli (2008), which models the SZ emission as a statistically isotropic signal modulated by a biased density contrast, where the bias depends on the cluster mass and redshift. To obtain an estimate of the contribution to the cross-spectrum at 217–857 GHz we assume that the measured cross-spectrum at 143 GHz is entirely due to thermal SZ emission (note that we do this to find what we believe to be an upper limit on the SZ contribution at 217–857 GHz; for the reasons stated above we do not expect the 143 GHz correlation to be due to SZ). Since the SZ signal at 143 GHz gives a decrement in the CMB, and the CIB emission gives an enhanced signal, it is possible that this approach could still underestimate the SZ signal. We find that in order to fit the cross-spectrum at 143 GHz using only the SZ-lensing correlation requires an amplitude of (2.4 ± 1.6) times our calculated SZ-lensing cross-spectrum. In Fig. 12 the dashed line shows the magnitude of this SZ signal scaled to each frequency using Eq. (10). The small contribution it makes at 217–857 GHz further suggests that we can neglect this component. At 217 GHz the signal is negative, while at higher frequencies it is positive.

5.3.4. ISW contamination

The integrated Sachs-Wolfe (ISW) effect describes the redshifting (blueshifting) of photons travelling through gravitational potential wells (hills) that decay as the photons travel through them (Sachs & Wolfe 1967). The induced modulation of the CMB mean by the gravitational potential generates CMB fluctuations that correlate with the lensing potential, which also traces out the gravitational potential perturbations (Seljak & Zaldarriaga 1999; Goldberg & Spergel 1999; Lewis et al. 2011). Note that because the mean of the CIB is relatively much smaller than its fluctuation, the ISW-induced CIB fluctuations make a negligible change to the total CIB anisotropy. The CMB ISW induced signal has the same frequency dependence as the CMB and so is only a significant contaminant for us at low frequencies. We evaluate this signal using a theoretical calculation performed in CAMB (Lewis et al. 2000). The results are shown as the solid line in Fig. 12; it is a negligible contribution at all frequencies, except in the lowest ℓ bin for the lowest frequencies, where the measured cross-spectrum is consistent with zero.

5.3.5. CIB bispectrum

Having calculated the bias from the point source shot noise in Sect. 5.3.2, we now discuss a more complicated form of the unresolved point source 3-point function that could be present in our data, namely the clustering contribution. While its amplitude is essentially unknown (although the first detection was recently reported in Crawford et al. 2013), the CIB bispectrum is potentially a direct contaminant to our measurement. Because of the nonlinear clustering of DSFGs (PER), it certainly has to exist. But because of the very large redshift kernel that characterizes the CIB, this non-Gaussian effect will be washed out, reducing its importance. Nevertheless, we ought to study it carefully.

If important, this effect would show up as a departure of the data from the best-fit curve in Fig. 10, since the best-fit model that we used assumes only a Poissonian shot-noise contribution. We do not see any significant deviation in Fig. 10. Still, in order to isolate this effect we create cross-spectra with increased sensitivity to the clustered point source signal. We do this by calculating the lens reconstruction at 100 GHz and 545 GHz, where, respectively, the radio and dusty point source contribution is stronger. The 545 GHz map has a much larger Galactic dust signal than our nominal 143 GHz map. However, unlike in our fiducial estimates, here we do not attempt to project out dust contamination from the map used to perform our lens reconstruction, because this would also remove some of the CIB signal in the bispectrum. As was found in Sect. 5.3.1, the cross-correlation between the 100 GHz reconstruction and the 100 GHz temperature map does not show any large difference with the cross-spectrum obtained using the 143 GHz signal. We are thus not sensitive enough to detect a bias from the clustering of radio sources using this method. However, we do detect a strong cross-correlation between the raw 545 GHz lens reconstruction and the 545 GHz temperature map. This cross-spectrum is shown in Fig. 11 for three different Galaxy masks. The line shows the point source shot-noise template derived in Sect. 5.3.2, fit to the cross-spectrum with the 10% Galaxy mask at ℓ > 1300. If the signal were entirely due to extragalactic point sources, then the signal would be independent of masking, and we do see a convergence of the signal at high ℓ as the size of the Galactic mask is increased. At low ℓ, however, there is a large Galactic contribution and the convergence with the reduced mask size is less clear. We thus conclude that a strong contribution from Galactic dust is present in this measurement, at all ℓ.

We do not attempt here to calculate accurately the shape of the clustering contribution to the CIB bispectrum, since it is beyond the scope of this work, even though a simple prescription for it has recently been proposed in Lacasa et al. (2012). To separate the Galactic from non-Galactic contributions in our bispectrum measurement is difficult, even if a strong Galactic signal is clearly present, given the strong dependence of the signal on variations of the Galactic mask in Fig. 11. However, the combination of dust cleaning that we perform in our nominal pipeline, coupled with the fact that our nominal pipeline uses the 143 GHz map for lens reconstruction, means that we do not observe any dependence with masking in our measurement, as seen in Fig. 7. Because of this, the CIB bispectrum is unlikely to make a large contribution to our measurement. Furthermore, even if we were to assume that all of the signal seen in Fig. 11 was extragalactic in nature, using the Gispert et al. (2000) frequency scaling for the CIB (also appropriate for the Galactic dust in fact, Planck Collaboration XXIV 2011), the roughly − 1700 μK.sr observed at ℓ = 400 for the 40% Galactic mask would only lead to a − 0.02 μK.sr signal in Fig. 12, which is an order of magnitude smaller than our measured signal. To conclude, although our analysis does not lead to a clean measurement of the CIB bispectrum, we can safely assume that it is not a contaminant to our CIB-lensing measurement.

5.4. Final statistical and systematic error budget

Throughout the suite of tests for instrumental and observational systematic errors presented in Sect. 5.2, as well as the suite of tests for possible astrophysical contaminants presented in Sect. 5.3, we have established the robustness of our measurement. The fact that our consistency tests do not lead to any significant deviation gives us confidence in our error budget. As described in Sect. 5.2 we add to them an overall calibration uncertainty, beam uncertainty, and lens normalization uncertainty, consistent with the Planck data processing analysis (Planck HFI Core Team 2011b). We gather the measured band-powers in Table 2, along with our statistical and systematic errors. These band-powers have been corrected for the point source component measured in Sect. 5.3.2, whose amplitude is also given in Table 2.

Once all systematic effects are factored in, we claim a detection significance of 3.6 (3.5), 4.3 (4.2), 8.3 (7.9), 31 (24), 42 (19), and 32 (16)σ statistical (statistical and systematic) at 100, 143, 217, 353, 545, and 857 GHz, respectively.

6. Interpretation and discussion

The correlation we have investigated leads to a very strong signal at most HFI frequencies. After a thorough examination of possible instrumental and astrophysical origins, we interpret it as originating from the spatial correlation between the sources of the CIB and the matter responsible for the gravitational deflection of CMB photons. In this section, we build on this result and interpret the measurement using both angular and frequency information.

Before doing so, we highlight the spectral information contained in the signal. It was shown in PER that both the frequency spectrum of the CIB mean and fluctuation rms are well approximated by the two modified blackbody spectra proposed by Fixsen et al. (1998) and Gispert et al. (2000), with a slight preference for the latter. We expect our measurements to follow the same spectral energy distribution (SED). Following the procedure outlined in Sect. 5.3.3, we plot in Fig. 9 the best-fit CIB component with either a Fixsen et al. (1998) – dot-dashed black line – or Gispert et al. (2000) – solid black line – SED. We do so for two ℓ bins. We can see that, indeed, for a given ℓ bin, our measurements qualitatively follow the expected CIB spectrum. Unlike PER, we do not find a preference for the Gispert et al. (2000) frequency shape in our low ℓ bin and only a slight one in the high-ℓ bin. It is worth emphasizing that by carrying out a cross-correlation measurement we can obtain constraints at the lowest frequency, which is usually heavily contaminated by the CMB (see P2013 for discussion). This is particularly interesting, because these measurements are simultaneously the most sensitive to high-z star formation processes and the most discrepant with either of the SEDs, i.e., they are both systematically low by about 0.5σ. The models presented in this section will allow us to use both the spectral dependence and the relative amplitudes of the ℓ bins that was lost in Fig. 9. We now describe the general methodology we will use, before describing our models in detail.

6.1. Model comparison methodology

For the purpose of model fitting, we will utilize both the CIB-lensing cross-spectra measured in this paper and the CIB auto-spectra obtained from P2013. We use the CIB-lensing cross-spectra for two purposes: to improve constraints on the model parameters; and to provide a consistency test of models fit to the CIB auto- and frequency cross-spectra alone. As can be seen in P2013, the cross-spectra of the CIB at different frequencies provide powerful constraints on the CIB emissivity.

We use the log-likelihood, $\begin{array}{ccc}\ln \mathrm{ℒ}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{\ell }}\sum _{\mathit{ij}}\mathrm{[}\mathit{Ĉ}\begin{array}{c}\mathit{i}\\ \mathit{\ell }\end{array}\mathrm{-}\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathit{i}}}_{\mathrm{˜}}}\mathrm{(}\mathit{p}{\mathrm{)}}^{\mathrm{]}}{\mathrm{(}{{\mathrm{N}}_{\mathit{\ell }}^{\mathit{ij}}}^{\mathrm{)}}}^{-1}{\mathrm{[}\mathit{Ĉ}\begin{array}{c}\mathit{j}\\ \mathit{\ell }\end{array}\mathrm{-}\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathit{j}}}_{\mathrm{˜}}}\mathrm{(}\mathit{p}{\mathrm{)}}^{\mathrm{]}}}^{†}\mathit{,}& & \end{array}$(11)where Ĉ and $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{C}\end{array}$ are the data and theory spectra with parameters p, the i and j indices denote the type of spectra (e.g., 100 GHz × φ or 100 GHz × 100 GHz), and N is the covariance matrix that includes both statistical and systematic errors. We make the approximation that the covariance matrix is diagonal, i.e., we treat the errors for different bins of each auto- and cross-spectrum as being uncorrelated. The small (≃2%) mask-induced mode-coupling between neighbouring bins supports this approximation. However, calibration and beam errors (which are correlated between the auto- and cross-spectra at a given frequency), as well as the lens normalization error (which is also correlated across spectra), are not accounted for in this approximation. In addition, the lens reconstruction has some sensitivity to all modes of the temperature maps, and so different φ modes are correlated to some degree. We also neglect the fact that the contribution to the error from the CIB signal itself (the orange line in Fig. 5) is also substantially correlated from frequency to frequency. However, our evaluation using simulations suggests that these effects are too small to significantly affect our procedure. We thus resort to simply adding the beam, calibration, and normalization uncertainties in quadrature to the statistical errors. The posterior probability distributions of model parameters are determined using now standard Markov chain Monte Carlo techniques (e.g., Knox et al. 2001; Lewis & Bridle 2002).

6.2. Two modelling approaches

The strength of the correlation signal should come as no surprise, given our current knowledge of CMB lensing and the CIB. The PER model predicts a high correlation between the CIB and the lensing potential. As clearly illustrated in Fig. 1, the broad overlap of the redshift distributions of the CIB with the lensing kernel, peaking at z ≈ 2–3, leads to a correlation of 60–80%. This is comparable to our measurements at all of the HFI frequencies, as illustrated in Fig. 13,

In models of the cross-correlation, the underlying properties we can probe come from a combination of the lensing potential and the characteristics of the DSFGs, in particular their emissivity and clustering properties. Mostly driven by linear physics, the former is well understood theoretically, as confirmed by recent observations (Smith et al. 2009; Hirata et al. 2008; Das et al. 2011; van Engelen et al. 2012). Assuming the currently favoured ΛCDM cosmology, we can consider it to be known to better than 10% in the multipole range of interest to us, an error dominated by the uncertainty in the normalization of the primordial power spectrum. Given that this is much smaller than the theoretical uncertainties related to DSFGs, we will fix the cosmology to the currently favoured Planck alone flat ΛCDM model in Planck Collaboration VII (2014), and will focus our analysis on the modelling of the DSFGs.

At a given redshift we model the fluctuations in the mean CIB emission, , as being proportional to the fluctuations in the number of galaxies, n_g (Haiman & Knox 2000), $\begin{array}{ccc}\mathit{\delta j}\mathrm{\propto }\mathit{j̅}\frac{\mathit{\delta }{\mathit{n}}_{\mathrm{g}}}{{\mathit{n}}_{\mathrm{g}}}\mathrm{·}& & \end{array}$(12)With this hypothesis, the goal of the CIB modelling becomes twofold: first, to better understand the statistical properties of the dusty galaxy number density, δn_g; and second, to reconstruct the mean emissivity of the CIB as a function of redshift.

In this paper we will use two different models of the CIB emission. The first model (described in Sect. 6.2.1 and inspired by Hall et al. 2010) uses a single bias parameter at all frequencies, with the mean CIB emissivity modelled as a two parameter Gaussian. This model is not designed to be physically realistic, and furthermore we will marginalize over an arbitrary amplitude in this case. Nevertheless, we present this simple model to show that our measurements are quite consistent with broad expectations of the CIB. The second model, described in Sect. 6.2.2, is a natural extension of the halo occupation density (HOD) approach used in PER (see also Pénin et al. 2012, and references therein). But unlike the results obtained in PER we now use a single HOD to describe the spectra at all frequencies. This is possible by allowing for deviations from the Béthermin et al. (2011) model (hereafter B11) that was used to fix the emissivity.

6.2.1. Linear bias model

As a first pass at interpreting our measurement, we will consider a redshift-independent linear bias model with a simple parametric SED. This model was found to provide a reasonable fit to the auto-spectra in the linear regime in PER. Throughout this paper we use the Limber approximation (see, e.g., Peebles 1980), and in this section, since we are using a linear model, we write the relevant angular spectra as $\begin{array}{ccc}{\mathit{C}}_{\mathit{\ell }}^{\mathit{XY}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{\chi }}_{\mathrm{\ast }}}\mathrm{d}\mathit{\chi }{\mathit{W}}^{\mathit{X}}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}{\mathit{W}}^{\mathit{Y}}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}{\mathit{P}}_{\mathit{\delta \delta }}\mathrm{(}\mathit{k}\mathrm{=}\mathit{\ell }\mathit{/}\mathit{\chi ,\chi }\mathrm{)}\mathit{,}& & \end{array}$(13)where X and Y are either the CIB at frequency ν or the lensing potential φ, the integral is over χ, the comoving distance along the line of sight, χ_∗ is the comoving distance to the last scattering surface, P_δδ(k,χ) is the matter power spectrum at distance χ, and the W^X functions are the redshift weights for each of the signals X: $\begin{array}{ccc}{\mathit{W}}^{\mathit{\nu }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}& \mathrm{=}& \mathit{b}\frac{\mathit{a}\mathit{j̅}\mathit{\nu }\mathrm{\left(}\mathit{\chi }\mathrm{\right)}}{\mathit{\chi }}\mathrm{;}\\ {\mathit{W}}^{\mathit{\phi }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{3}}{{\mathit{\ell }}^{\mathrm{2}}}\hspace{0.17em}{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}\hspace{0.17em}\frac{\mathit{\chi }}{\mathit{a}}\left(\frac{{\mathit{\chi }}_{\mathrm{\ast }}\mathrm{-}\mathit{\chi }}{{\mathit{\chi }}_{\mathrm{\ast }}\mathit{\chi }}\right)\mathrm{·}\end{array}$(14)Here b is the DSFG bias that we assume to be redshift-independent, a is the scale factor, is the mean CIB emissivity at frequency ν, as defined in PER, Ω_m is the matter density today in critical density units, and H₀ is the Hubble parameter today. We use the Hall et al. (2010) model for the CIB kernel, which treats the CIB emissivity as a Gaussian in redshift: $\begin{array}{ccc}\mathit{j̅}\mathit{\nu }\mathrm{\left(}\mathit{\chi }\mathrm{\right)}\mathrm{\propto }\mathit{a}{\mathit{\chi }}^{\mathrm{2}}\exp \mathrm{[}\mathrm{-}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{c}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}\mathrm{2}{{\mathit{\sigma }}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{]}}{\mathit{f}}_{\mathit{\nu }\mathrm{(}\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{z}\mathrm{)}}\mathit{,}& & \end{array}$(15)where we use a modified blackbody frequency dependence $\begin{array}{ccc}{\mathit{f}}_{\mathit{\nu }\mathrm{(}\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{z}\mathrm{)}}\mathrm{\propto }{\mathit{\nu }}^{\hspace{0.17em}\mathit{\beta }}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}\mathit{.}& & \end{array}$(16)We fix the dust temperature to T_d = 34 K, the spectral index to β = 2 (Hall et al. 2010), and assume a constant bias b. We include a single normalization parameter for j, which we marginalize over. Since the normalization and bias parameters are degenerate in Eq. (13), if we were to only use the measured auto- and cross-spectra this approach would be equivalent to marginalizing over a frequency-independent bias parameter. However, we will further constrain our model using the FIRAS data, which breaks this degeneracy. We constrain the z_c and σ_z parameters at each frequency, giving us a total of 13 free parameters.

For 217–857 GHz, we use the FIRAS measurements of the CIB mean intensity from Lagache et al. (2000) as an additional constraint on our model. The mean intensity is simply $\begin{array}{ccc}{\mathit{I}}_{\mathit{\nu }}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{\chi }}_{\mathrm{\ast }}}\mathrm{d}\mathit{\chi a}\mathit{j̅}\mathit{\nu }\mathrm{\left(}\mathit{\chi }\mathrm{\right)}\mathit{.}& & \end{array}$(17)Using this equation and the measured FIRAS mean and uncertainty we calculate a χ² value and add this to the χ² in Eq. (11). Since there are no FIRAS constraints at 100 and 143 GHz, as well as no CIB auto-spectra measurements, and noisier cross-spectra measurements at these frequencies, our constraints for the 100 and 143 GHz redshift parameters are weaker than for the other parameters.

The linear bias model considers only linear clustering, and so when fitting the auto-spectra we restrict ourselves to ℓ < 500, where nonlinear contributions are negligible. Because we do not consider the high-ℓ modes, we also neglect the shot-noise contribution to the auto-spectra. The best-fit model is shown as the coloured dashed lines in Fig. 15, with χ² values of 13.4, 16.8, 25.2, 21.8, 9.1, and 9.4, if we break up the χ² contribution per frequency from 100 to 857 GHz, leading to an overall χ² of 95.7 for N_d.o.f. = 59. We see that the model captures some characteristics of the data, but we also have evidence it is significantly missing some features as well. This is perhaps not surprising, given the simplicity of the model. The two-dimensional marginal distributions of the z_c and σ_z parameters are shown in Fig. 14. Although we allowed for these parameters to be frequency dependent we note that the point z_c = 1 and σ_z = 2.2 is in a region of high probability at all frequencies, and gives a redshift distribution for the emissivity density roughly consistent with our expectations, rising towards z = 1 due to the χ² term and then only slowly falling off towards higher redshifts.

Althought it is useful to see the extent to which such a simple model can explain our data, we now turn to make a stronger connection between the properties of the infrared light and the distribution of the underlying dark matter applicable into the nonlinear regime.

6.2.2. An extended halo model based analysis

In this section we use a description of the CIB motivated by the halo model, which has been used successfully to describe the transition between the linear and nonlinear clustering regimes for optical galaxies. We use the halo model to attempt to reconstruct the CIB emissivity as a function of redshift. This is an extension of the approach taken in PER, where the modelled CIB emissivity at high redshift was treated as a single bin, with amplitude constrained by the data. The goal of this approach is to isolate the high-redshift contribution to the CIB, which is difficult to probe using observations of individual galaxies, due to their low brightness. The power of such an approach is further demonstrated in P2013.

We replace the linear bias used in Sect. 6.2.1 with a halo model and an HOD prescription that assigns galaxies to host dark matter halos (see PER for references and definitions). This allows a consistent description of the linear and nonlinear parts of the galaxy power spectrum and its redshift evolution. Because it is built on the clustering of dark matter halos, the halo model allows us to describe in a consistent way the clustering of DSFGs and the gravitational lensing caused by the halos. However, it is important to realise that the HOD prescription was developed to describe stellar mass within dark matter halos – an application for which it has been thoroughly tested – while here we are applying it to star formation within halos. The accuracy of this approach needs to be further quantified. However, it provides a good phenomenological description of our data and other CIB measurements, as well as other astrophysical probes of the relation between dark matter and light (e.g., Leauthaud et al. 2012; Hikage et al. 2013).

Unlike the model presented in PER we use a single HOD to describe our data at all frequencies. This is made possible by allowing for a deviation from the B11 emissivity model. Note however that we will still consider the CIB emissivity to depend only on redshift and not on the galaxy host halo mass, a simplification highlighted in Shang et al. (2012) that will be relaxed in the P2013 model. The emissivity of the CIB is inhomogeneous, due to spatial variations in the number density of galaxies: $\begin{array}{ccc}\frac{\mathit{\delta }{\mathit{j}}_{\mathit{\nu }}}{\mathit{j̅}\mathit{\nu }}\mathrm{\left(}\widehat{{n}\hspace{0.17em}}\mathit{,z}\mathrm{\right)}\mathrm{=}\frac{\mathit{\delta }{\mathit{n}}_{\mathrm{g}}}{\mathit{n̅}\mathrm{g}}\mathrm{\left(}\widehat{{n}\hspace{0.17em}}\mathit{,z}\mathrm{\right)}\mathrm{\equiv }{\mathit{\delta }}_{\mathrm{g}}\mathrm{\left(}\widehat{{n}\hspace{0.17em}}\mathit{,z}\mathrm{\right)}\mathit{.}& & \end{array}$(18)Here $\mathit{j}\mathrm{\left(}\widehat{{n}\hspace{0.17em}}\mathit{,z}\mathrm{\right)}$ is the CIB emissivity at redshift z with mean , ${\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\widehat{{n}\hspace{0.17em}}\mathit{,z}\mathrm{\right)}$ is the number density of DSFGs with mean , and ${\mathit{\delta }}_{\mathrm{g}}\mathrm{\left(}\widehat{{n}\hspace{0.17em}}\mathit{,z}\mathrm{\right)}$ is the DSFG overdensity, with power spectrum ⟨δ_g(k,z) δ_g(k^′,z)^∗⟩ = (2π)³δ(k − k^′) P_gg(k,z). We calculate this power spectrum, including the constituent 1- and 2-halo terms, using the procedure described in Appendix C of PER, with the constraint α_sat = 1, a theoretically favoured value (Tinker & Wetzel 2010). We remove the relationship between M_sat, a characteristic satellite mass scale, and M_min, the halo mass at which a halo has a 50% probability of containing a central galaxy that was enforced in PER (i.e., M_sat = 3.3M_min), and allow both M_sat and M_min to vary independently.

At redshift z< 1 we fix the emissivity to the B11 value, but at higher redshift we assume that the emissivity is constant within z bins and solve for the amplitude of the bins. Two factors affect the number of bins that we choose. The auto-spectra have a dependence, and so if the true value of has a strong z dependence within a bin then the best-fit emissivity in the bin will be difficult to interpret; the best-fit bin values could be significantly different from those that would be obtained by binning the true emissivity. However, as more bins are used and the number of parameters increases, it becomes more difficult to determine the best-fit parameters and the parameters will be highly correlated. After investigation using simulations, we found that three bins was a good compromise, given the expected slow redshift evolution. The bins are defined by: 1 <z ≤ 1.5; 1.5 <z ≤ 3; and 3 <z ≤ 7. As in Sect. 6.2.1 we use the FIRAS results at 217–857 GHz to add an integral constraint on the emissivity. The CIB auto- and lensing cross-spectra are (Song et al. 2003): $\begin{array}{ccc}{\mathit{C}}_{\mathit{\ell }}^{\mathit{\nu }{\mathit{\nu }}^{\mathrm{\prime }}}& \mathrm{=}& \int \mathrm{d}\mathit{\chi }{\mathit{W}}^{\mathit{\nu }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}{\mathit{W}}^{{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}{\mathit{P}}_{\mathrm{gg}}\mathrm{(}\mathit{k}\mathrm{=}\mathit{\ell }\mathit{/}\mathit{\chi ,\chi }\mathrm{)}\mathrm{;}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{\nu \phi }}& \mathrm{=}& \mathrm{\int }\mathrm{d}\mathit{\chi }{\mathit{W}}^{\mathit{\nu }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}{\mathit{W}}^{\mathit{\phi }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}{\mathit{P}}_{\mathit{\delta }\mathrm{g}}\mathrm{(}\mathit{k}\mathrm{=}\mathit{\ell }\mathit{/}\mathit{\chi ,\chi }\mathrm{)}\mathit{.}\end{array}$(19)Since we fix at z< 1, the model spectra consist of a low redshift part that is independent of the emissivity parameters, and a contribution from z> 1 that is proportional to for the auto-spectra and for the lensing cross-spectra.

Overall, the halo-based model contains two halo parameters that describe the galaxy clustering and are independent of frequency, and three j amplitudes at each frequency, giving a total of 20 parameters for the six frequencies of interest to us. The auto- and cross-spectra have a total of 120 ℓ bins, with four additional FIRAS data points. Solving for the likelihood described in Sect. 6.1, gives the best-fit models shown in Figs. 15 and 16 as the solid coloured lines. The χ² values in each panel are the contribution to the total χ² from the data within the panel. The combined reduced χ² is 102.1 for N_d.o.f. = 104, indicating a good fit. We force M_sat ≥ M_min in the MCMC fitting procedure. The best-fit parameters are log ₁₀(M_min/M_⊙) = 12.2 and log ₁₀(M_sat/M_⊙) = 12.8, which gives M_sat/M_min = 3.8. The mean parameter values are log ₁₀(M_min/M_⊙) = 10.5 ± 0.6 and log ₁₀(M_sat/M_⊙) = 10.8 ± 0.7. The best-fit value of M_min is consistent with that derived in PER at multiple frequencies, where a similar model was considered (even though we now set α_sat = 1 and reconstruct the mean emissivity as a function of redshift). The associated mean emissivity parameters are given in Table 3 and displayed in Fig. 17, where we also plot the B11 model for reference. As can be seen in Fig. 17, we remain consistent with the B11 model in most redshift bins.

6.3. Interpreting the reconstructed emissivities

We now illustrate one interesting consequence of this measurement and show how using the constrained emissivities, j_ν(z), we can estimate the star formation rate (SFR) density at different redshifts. Following Pénin et al. (2012), we begin by writing the emissivity as an integral over the galaxy flux densities: $\begin{array}{ccc}{\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\left(\mathit{a}\frac{\mathrm{d}\mathit{\chi }}{\mathrm{d}\mathit{z}}\right)}^{-1}\mathrm{\int }{\mathit{S}}_{\mathit{\nu }}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{N}}{\mathrm{d}{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{z}}\mathrm{d}{\mathit{S}}_{\mathit{\nu }}\mathit{.}& & \end{array}$(20)Here S_ν is the flux density, and d²N/ dS_ν dz is the number of galaxies per flux density element and redshift interval. The galaxies contributing to the CIB can be divided into various populations (labelled as p) based on the galaxy SED, e.g., according to galaxy type or dust temperature: $\begin{array}{ccc}{\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\left(\mathit{a}\frac{\mathrm{d}\mathit{\chi }}{\mathrm{d}\mathit{z}}\right)}^{-1}\sum _{\mathit{p}}\mathrm{\int }{\mathit{S}}_{\mathit{\nu }}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{N}}_{\mathit{p}}}{\mathrm{d}{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{z}}\mathrm{d}{\mathit{S}}_{\mathit{\nu }}\mathit{.}& & \end{array}$(21)If we define s_ν as the flux density of an L_IR = L_⊙ source with the SED of a given population, i.e., S_ν = s_νL_IR (with L_IR in units of L_⊙), then we can write Eq. (21) as (Pénin et al. 2012) $\begin{array}{ccc}{\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\left(\mathit{a}\frac{\mathrm{d}\mathit{\chi }}{\mathrm{d}\mathit{z}}\right)}^{-1}\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{z}}\hspace{0.17em}\sum _{\mathit{p}}{\mathit{s}}_{\mathit{\nu }}\mathrm{\int }{\mathit{L}}_{\mathrm{IR}}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{N}}_{\mathit{p}}}{\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\hspace{0.17em}\mathrm{d}\mathit{V}}\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\mathit{.}& & \end{array}$(22)The contribution to the infrared luminosity density from a given population is $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{IR}\mathit{,p}}\mathrm{=}\mathrm{\int }{\mathit{L}}_{\mathrm{IR}}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{N}}_{\mathit{p}}}{\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\hspace{0.17em}\mathrm{d}\mathit{V}}\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\mathit{.}& & \end{array}$(23)We assume a simple conversion between L_IR and the star formation rate density, ρ_SFR, using the Kennicutt constant K (Kennicutt 1998). Since by definition ρ_SFR = K ∑ _pρ_IR,p, we can rewrite Eq. (22) as $\begin{array}{ccc}{\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\hspace{0.17em}{\mathit{\chi }}^{\mathrm{2}}\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{SFR}}}{\mathit{K}}\left(\frac{{\sum }_{\mathit{p}}{\mathit{s}}_{\mathit{\nu }}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{IR}\mathit{,p}}}{{\sum }_{\mathit{p}}{\mathit{\rho }}_{\mathrm{IR}\mathit{,p}}}\right)\mathit{,}& & \end{array}$(24)where the final term in brackets is the effective SED of infrared galaxies, which we write as s_ν,eff. We derive these SEDs following the evolution model of Béthermin et al. (2012a) using Magdis et al. (2012) templates. The construction of these effective SEDs will be explained in detail in future work. Finally, we obtain the conversion factor between mean emissivity and SFR density, $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{SFR}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathit{K}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\hspace{0.17em}{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\hspace{0.17em}{\mathit{S}}_{\mathit{\nu ,}\mathrm{eff}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\hspace{0.17em}{\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{.}& & \end{array}$(25)Using Eq. (25) we can find the coefficients for each of the redshift bins and frequencies, which are presented in Table 3.

6.4. Discussion and outlook

In the previous section we described a model that simultaneously fits the CIB auto-spectra and the CIB-lensing cross-spectra, at all frequencies and with a single HOD prescription. Given that we use an emissivity function that is close to the B11 emissivities (to within our uncertainties), we expect predictions of the galaxy number counts derived from our best-fit emissivity to agree with current estimates (Béthermin et al. 2012b). The fact that our measurement is consistent with previous models of the CIB lends support to our current understanding of its origin. For example, the characteristic mass scale at which halos host galaxies, M_min, is consistent with values derived in PER, and is consistent with, but slightly higher than, the value derived more recently in Viero et al. (2013), log ₁₀(M_min/M_⊙) = 10.5 ± 0.5 (although a direct comparison could be misleading given the different model assumptions). In particular, it is clear that our model has limitations, some of which have been partially addressed in recent work (Shang et al. 2012; Béthermin et al. 2012c; De Bernardis & Cooray 2012; Béthermin et al. 2012a; Viero et al. 2013; Addison et al. 2013) and are points of focus in P2013, among them the mass independence of the emissivity and the redshift evolution of the SED. In P2013 it is also shown that the value of M_min is in fact poorly constrained in this more complex approach. As such, although reasonable, we interpret the constraints as model-dependent. Another question worth further investigation is the dependence of our results on the binning scheme chosen for the emissivity, which will be addressed in a future paper.

Given the consistency of our model with the PER results, the information added by our cross-spectrum measurement is worth quantifying. As an example, we show in Fig. 18 the highest redshift emissivity bin in the 353 GHz band. Adding the CIB-lensing cross-spectrum information tightens the constraint on the high-redshift part of the emissivity. This statement also holds for the other frequencies and stems from the fact that the CMB lensing kernel peaks at relatively high redshift, making the cross-correlation more sensitive to the high-redshift CIB signal than the CIB auto-spectrum, as is illustrated in Fig. 1. Although this gain does not translate into a substantial improvement in the determination of M_min, it leads to interesting constraints on the SFR density, as can be seen in Table 3.

The results at frequencies above 217 GHz each lead to around 2σ evidence for a non-zero SFR density for 1.5 <z< 3 and for 3 <z< 7. The values inferred are consistent with other probes of the SFR in these redshift ranges, as compiled for example in Fig. 1 of Hopkins & Beacom (2006). Assuming that each frequency is independent, we obtain SFR densities for the three redshift bins of (0.42 ± 0.12), (0.29 ± 0.14), and (0.23 ± 0.10) M_⊙ Mpc^-3 yr^-1, respectively, where the errors are 68% C.L. We note that the j distributions are rather non-Gaussian so that the 95% C.L. become 0.23, 0.25, and 0.19, respectively. This roughly 2σ detection per bin compares very favourably with other published measurements. The constraints clearly illustrate how this particular correlation can be used to better isolate the high redshift component of the CIB and improve our understanding of the star formation rate at high redshift. Such constrains will improve with future measurements, in particular if we can increase the signal-to-noise ratio in our lower frequency channels, where the high redshift contribution is the greatest. This will likely require an accurate removal of the CMB, our dominant source of noise at low frequencies. A more thorough discussion of this possibility is given in P2013.

To fully utilize the richness of the CIB-lensing correlation will require more studies. Future work could involve using more sophisticated halo models, specifically designed to model star formation within halos, as well as relaxing some of the assumptions made here, such as the mass-independent luminosity function. In addition the use of map-based methods, which enable estimates of the galaxy host halo mass by stacking the lensing potential maps, is worth pursuing, as is the extension to other data-sets. For illustration purposes, we show in Fig. 19 raw measurements of the correlation between our lensing potential map and the IRIS map at 100 μm that was introduced in Sect. 2.2.2. We use our nominal mask and lens reconstruction, with no H i cleaning performed on the IRIS map. We clearly see a strong correlation, whose significance we estimate to be 9σ, ignoring any possible systematic effects. To guide the eye we have added a prediction (not a fit) based on the HOD model presented in Pénin et al. (2012). The full analysis of this signal is beyond the scope of this paper, but it illustrates possible future uses of the lensing potential map. In this case the IRIS wavelength range will help us to isolate the low-redshift contribution to the CIB.

To conclude, we have presented the first measurement of the correlation between lensing of the CMB and the CIB. Planck’s unique full-sky, multi-frequency, deep survey enables us to make an internal measurement of this correlation. Measurements with high statistical significance are obtained, even after accounting for possible systematic errors. The high degree of correlation that is measured (around 80%) allows for unprecedented visualization of lensing of the CMB and holds great promise for new CIB and CMB focused science. CMB lensing appears promising as a probe of the origin of the CIB, while the CIB is now established as an ideal tracer of CMB lensing.

---

Acknowledgments

Based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN and JA (Spain); Tekes, AoF and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and 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. We acknowledge the use of the HEALPix package, and the LAMBDA archive (http://lambda.gsfc.nasa.gov).

References

Appendix A: Statistical errors

In this section we compare six different methods to estimate our statistical errors. This comparison is used to validate our claim (presented in the main text) that we can obtain our uncertainties from the naive analytical errors calculated from the data, i.e., Method 1 below. The six different methods we compare are as follows.

• 1.

  The naive, Gaussian, analytical errors estimated from the datathrough the measured total power of the T and φ fields, respectively, $\mathit{Ĉ}\begin{array}{c}\mathrm{TT}\\ \mathit{\ell }\end{array}$, $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{\ell }\end{array}$, and the model cross-spectrum ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{T}\mathit{\phi }}$.

• 2.

  As above, but instead of the data maps we use one of our simulations of the CIB and CMB temperature maps described in Sect. 3.4. The lens reconstruction is obtained from the simulated maps using the same procedure that we use for the data.

• 3.

  The scatter directly within individual ℓ bins in the data-determined cross-spectrum.

• 4.

  As above, but the scatter is measured in each ℓ bin for each simulation realization, and the errors are averaged over 100 realizations.

• 5.

  The scatter of the bins is calculated using the cross-spectra measured from our simulated maps. This is a direct measurement of the statistical error we require (to the extent that our simulations are realistic), and will differ from the scatter within the ℓ bins, for example due to noise correlations between different multipoles of the cross-spectrum.

• 6.

  The error in the cross-spectrum of the reconstructed lensing potential in simulated maps with the measured temperature maps. This will only give part of the contribution to the error, since the temperature maps are fixed, but it is still a useful cross-check.

In Fig. A.1 we show a comparison of the errors found from our six measurement methods. The precision achievable with 100 simulations is indicated by the grey envelope. We show the errors in each ℓ bin from the different methods divided by the data-based analytical estimate. To discuss the implications of these results we shall focus on the 100 GHz panel first. The scatter measured within an ℓ bin is fairly consistent in the simulations (Method 4, black solid line) and in the data (Method 3, black dashed line) giving us confidence in our simple simulation procedure. This rules out important systematic contributions and shows that our signal is mostly Gaussian, as expected. Note that the consistency with the simulations is not surprising, since at low frequencies we are dominated by CMB and instrumental noise, which are well understood. In addition, the fact that the analytical errors calculated on the simulations (Method 2, coloured dashed line) are mostly within the grey shaded region, and are therefore close to the analytical errors calculated from the data (Method 1), gives us further confidence in the simulations. To the extent that the simulations are accurate, the scatter of the ℓ bins in simulations (solid coloured line) is the error that we require. The fact that it is essentially all within the shaded envelope means that this method gives errors that are consistent with the analytical errors measured using the data, justifying our nominal choice for calculating the errors at low frequencies.

However, comparing the black lines with the coloured lines clearly indicates that, in both the data and simulations, obtaining the error bars by measuring the scatter within the ℓ bins leads to an underestimation of the errors by approximately 20%. Given the fact that this difference is observed in both the data and the simulations, we exclude any instrumental systematic effect as its cause and explain it as being due to noise correlations within the ℓ bins. Such a correlation is expected, since most of the lens reconstruction signal in the ℓ-range of interest to us comes from modes in the CMB map within a relatively narrow range at ℓ ≃ 1500, and so multipoles in the lens reconstruction are correlated. We have also checked that the mask induced ℓ bin coupling is negligible, given the bin width we have chosen, and is always smaller than 2%.

All of these conclusions remain valid up to 353 GHz. However, at 545 and 857 GHz, we see, by looking at the errors measured using simulations (solid black for Method 4, solid coloured for Method 5, and dashed coloured for Method 2), that the errors deduced from the analytical estimates measured from the data are substantially higher than those we measure in the simulations. This is easily explained through the fact that we are omitting any foreground emission in our simulations. The relative contribution of Galactic foreground emission is higher at low ℓ, which is expected because the Galactic cirrus emission has a steeply declining power spectrum. Overall, the amplitude of this contribution is also consistent with what is seen in Fig. 5.

Since the scatter within the ℓ bins measured in simulations (black dashed line) is about 20% lower than the data-based analytical estimates at all frequencies, we use the data-based analytical estimates as the basis for our statistical errors. We could alternatively scale the scatter-determined errors by 20% and obtain consistent results. This approach accounts for the foreground emission seen at 545 GHz and 857 GHz, but will in practice neglect the contribution to the errors from the non-Gaussian part of the foregrounds. However, we show in Sect. 5.3.5 that this contribution is small enough that we can ignore it.

The remaining method to discuss is obtained from the cross-spectrum of the reconstructed lensing potential in simulated maps with the data-measured temperature maps (Method 6, coloured dashed-dotted line in Fig. A.1). At 545 and 857 GHz the CIB signal is dominant over a large ℓ range, and so the error obtained from this method is equal to the “signal” terms in Eq. (4), which are the orange and green lines in Fig. 5. These two lines make up a significant fraction of the total error and provide a reasonable approximation to the true error at high-ℓ. However, at low-ℓ, where Galactic emission is important, and at 100–353 GHz where the CMB and instrumental noise are the largest components, the orange and green curves do not accurately describe the total error. We can see from Fig. A.1 that the errors obtained using this method are close to the errors measured using the other techniques. However, this method will underestimate the true errors, since the scatter in the CMB and noise components is neglected.

Note that the results presented in Fig. A.1 are all computed using the 40% Galaxy mask, but we have checked that they also hold when using the 20% and 60% Galaxy masks (which are discussed in detail in Sect. 5.2) and that the results show the appropriate f_sky scaling.

All Tables

All Figures