Issue 
A&A
Volume 571, November 2014
Planck 2013 results



Article Number  A18  
Number of page(s)  24  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201321540  
Published online  29 October 2014 
Planck 2013 results. XVIII. The gravitational lensinginfrared background correlation
^{1}
APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3,
CEA/lrfu, Observatoire de Paris, Sorbonne Paris Cité, 10 rue Alice Domon et Léonie Duquet,
75205
Paris Cedex 13,
France
^{2}
Aalto University Metsähovi Radio Observatory,
Metsähovintie 114, 02540
Kylmälä,
Finland
^{3}
African Institute for Mathematical Sciences,
68 Melrose Road, Muizenberg,
7701
Rondebosch Cape Town, South
Africa
^{4}
Agenzia Spaziale Italiana Science Data Center, via del Politecnico
snc, 00133
Roma,
Italy
^{5}
Agenzia Spaziale Italiana, viale Liegi 26, 00198
Roma,
Italy
^{6}
Astrophysics Group, Cavendish Laboratory, University
ofCambridge, J J Thomson
Avenue, Cambridge
CB3 0HE,
UK
^{7}
Astrophysics & Cosmology Research Unit, School of Mathematics,
Statistics & Computer Science, University of KwaZuluNatal,
Westville Campus, Private Bag
X54001, 4000
Durban, South
Africa
^{8}
Atacama Large Millimeter/submillimeter Array, ALMA Santiago Central
Offices, Alonso de Cordova 3107,
Vitacura, Casilla
763 0355
Santiago,
Chile
^{9}
CITA, University of Toronto, 60 St. George St., Toronto, ON
M5S 3H8,
Canada
^{10}
CNRS, IRAP, 9 Av.
colonel Roche, BP
44346, 31028
Toulouse Cedex 4,
France
^{11}
California Institute of Technology, Pasadena, California, USA
^{12}
Centre for Theoretical Cosmology, DAMTP, University of
Cambridge, Wilberforce
Road, Cambridge
CB3 0WA,
UK
^{13}
Centro de Estudios de Física del Cosmos de Aragón
(CEFCA), Plaza San Juan 1, planta
2, 44001
Teruel,
Spain
^{14}
Computational Cosmology Center, Lawrence Berkeley National
Laboratory, Berkeley,
California,
USA
^{15}
Consejo Superior de Investigaciones Científicas
(CSIC), 28006
Madrid,
Spain
^{16}
DSM/Irfu/SPP, CEASaclay, 91191
GifsurYvette Cedex,
France
^{17}
DTU Space, National Space Institute, Technical University of
Denmark, Elektrovej
327, 2800
Kgs. Lyngby,
Denmark
^{18}
Département de Physique Théorique, Université de
Genève, 24 Quai E.
Ansermet, 1211
Genève 4,
Switzerland
^{19}
Departamento de Física Fundamental, Facultad de Ciencias,
Universidad de Salamanca, 37008
Salamanca,
Spain
^{20}
Departamento de Física, Universidad de Oviedo,
Avda. Calvo Sotelo s/n,
93007
Oviedo,
Spain
^{21}
Department of Astronomy and Astrophysics, University of
Toronto, 50 Saint George Street,
Toronto, Ontario,
Canada
^{22}
Department of Astrophysics/IMAPP, Radboud University
Nijmegen, PO Box
9010, 6500 GL
Nijmegen, The
Netherlands
^{23}
Department of Electrical Engineering and Computer Sciences,
University of California, Berkeley, California, USA
^{24}
Department of Physics & Astronomy, University of British
Columbia, 6224 Agricultural Road,
Vancouver, British
Columbia, Canada
^{25}
Department of Physics and Astronomy, Dana and David Dornsife College
of Letter, Arts and Sciences, University of Southern California,
Los Angeles
CA
90089,
USA
^{26}
Department of Physics and Astronomy, University College
London, London
WC1E 6BT,
UK
^{27}
Department of Physics, Florida State University,
Keen Physics Building, 77 Chieftan
Way, Tallahassee,
Florida,
USA
^{28}
Department of Physics, Gustaf Hällströmin katu 2a, University of
Helsinki, Helsinki,
Finland
^{29}
Department of Physics, Princeton University,
Princeton, New Jersey, USA
^{30}
Department of Physics, University of California,
One Shields Avenue, Davis, California, USA
^{31}
Department of Physics, University of California,
Santa Barbara, California, USA
^{32}
Department of Physics, University of Illinois at
UrbanaChampaign, 1110 West Green
Street, Urbana,
Illinois,
USA
^{33}
Dipartimento di Fisica e Astronomia G. Galilei, Università degli
Studi di Padova, via Marzolo
8, 35131
Padova,
Italy
^{34}
Dipartimento di Fisica e Scienze della Terra, Università di
Ferrara, via Saragat
1, 44122
Ferrara,
Italy
^{35}
Dipartimento di Fisica, Università La Sapienza,
P.le A. Moro 2, 00185
Roma,
Italy
^{36}
Dipartimento di Fisica, Università degli Studi di
Milano, via Celoria
16, 20133
Milano,
Italy
^{37}
Dipartimento di Fisica, Università degli Studi di
Trieste, via A. Valerio
2, 34127
Trieste,
Italy
^{38}
Dipartimento di Fisica, Università di Roma Tor
Vergata, via della Ricerca Scientifica
1, 00133
Roma,
Italy
^{39}
Discovery Center, Niels Bohr Institute, Blegdamsvej 17, 2100
Copenhagen,
Denmark
^{40}
Dpto. Astrofísica, Universidad de La Laguna (ULL),
38206, La Laguna, Tenerife,
Spain
^{41}
European Southern Observatory, ESO Vitacura,
Alonso de Cordova 3107, Vitacura,
Casilla
19001
Santiago,
Chile
^{42}
European Space Agency, ESAC, Planck Science Office,
Camino bajo del Castillo s/n, Urbanización
Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain
^{43}
European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ
Noordwijk, The
Netherlands
^{44}
Finnish Centre for Astronomy with ESO (FINCA), University of
Turku, Väisäläntie
20, 21500
Piikkiö,
Finland
^{45}
Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University
of Helsinki, 00014
Helsinki,
Finland
^{46}
INAF – Osservatorio Astronomico di Padova,
Vicolo dell’Osservatorio 5,
35122
Padova,
Italy
^{47}
INAF – Osservatorio Astronomicodi Roma, via di Frascati 33, 00040
Monte Porzio Catone,
Italy
^{48}
INAF – Osservatorio Astronomico di Trieste,
via G.B. Tiepolo 11, 34131
Trieste,
Italy
^{49}
INAF Istituto di Radioastronomia, via P. Gobetti 101, 40129
Bologna,
Italy
^{50}
INAF/IASF Bologna, via Gobetti 101, 40129
Bologna,
Italy
^{51}
INAF/IASF Milano, via E. Bassini 15, 20133
Milano,
Italy
^{52}
INFN, Sezione di Bologna, via Irnerio 46,
40126
Bologna,
Italy
^{53}
INFN, Sezione di Roma 1, Università di Roma Sapienza,
Piazzale Aldo Moro 2,
00185
Roma,
Italy
^{54}
IPAG: Institut de Planétologie et d’Astrophysique de Grenoble,
Université Joseph Fourier, Grenoble 1/CNRSINSU, UMR 5274, 38041
Grenoble,
France
^{55}
IUCAA, Post Bag 4, Ganeshkhind, Pune University
Campus, 411 007
Pune,
India
^{56}
Imperial College London, Astrophysics group, Blackett
Laboratory, Prince Consort
Road, London,
SW7 2AZ,
UK
^{57}
Infrared Processing and Analysis Center, California Institute of
Technology, Pasadena
CA
91125,
USA
^{58}
Institut Néel, CNRS, Université Joseph Fourier Grenoble
I, 25 rue des
Martyrs, 38042
Grenoble,
France
^{59}
Institut Universitaire de France, 103 Bd SaintMichel, 75005
Paris,
France
^{60}
Institut d’Astrophysique Spatiale, CNRS (UMR 8617) Université
ParisSud 11, Bâtiment
121, 91405
Orsay,
France
^{61}
Institut d’Astrophysique de Paris, CNRS (UMR 7095),
98bis Bd Arago, 75014
Paris,
France
^{62}
Institute for Space Sciences, 077125
BucharestMagurale,
Romania
^{63}
Institute of Astronomy and Astrophysics, Academia
Sinica, 10617
Taipei,
Taiwan
^{64}
Institute of Astronomy, University of Cambridge,
Madingley Road, Cambridge
CB3 0HA,
UK
^{65}
Institute of Theoretical Astrophysics, University of
Oslo, Blindern,
0315
Oslo,
Norway
^{66}
Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, 38200, La
Laguna, Tenerife, Spain
^{67}
Instituto de Física de Cantabria (CSICUniversidad de
Cantabria), Avda. de los Castros
s/n, 39005
Santander,
Spain
^{68}
Jet Propulsion Laboratory, California Institute of
Technology, 4800 Oak Grove
Drive, Pasadena,
California,
USA
^{69}
Jodrell Bank Centre for Astrophysics, Alan Turing Building, School
of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13
9PL, UK
^{70}
Kavli Institute for Cosmology Cambridge,
Madingley Road, Cambridge, CB3 0HA, UK
^{71}
LAL, Université ParisSud, CNRS/IN2P3, Orsay, France
^{72}
LERMA, CNRS, Observatoire de Paris, 61 avenue de
l’Observatoire, 75014
Paris,
France
^{73}
Laboratoire AIM, IRFU/Service d’Astrophysique – CEA/DSM – CNRS –
Université Paris Diderot, Bât. 709,
CEASaclay, 91191
GifsurYvette Cedex,
France
^{74}
Laboratoire Traitement et Communication de l’Information, CNRS (UMR
5141) and Télécom ParisTech, 46 rue
Barrault, 75634
Paris Cedex 13,
France
^{75}
Laboratoire de Physique Subatomique et de Cosmologie, Université
Joseph Fourier Grenoble I, CNRS/IN2P3, Institut National Polytechnique de
Grenoble, 53 rue des
Martyrs, 38026
Grenoble Cedex,
France
^{76}
Laboratoire de Physique Théorique, Université ParisSud 11 &
CNRS, Bâtiment 210,
91405
Orsay,
France
^{77}
Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{78}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741
Garching,
Germany
^{79}
McGill Physics, Ernest Rutherford Physics Building, McGill
University, 3600 rue
University, Montréal,
QC, H3A 2T8, Canada
^{80}
MilliLab, VTT Technical Research Centre of Finland,
Tietotie 3, 05044
Espoo,
Finland
^{81}
Niels Bohr Institute, Blegdamsvej 17, 2100
Copenhagen,
Denmark
^{82}
Observational Cosmology, Mail Stop 36717, California Institute of
Technology, Pasadena
CA
91125,
USA
^{83}
Optical Science Laboratory, University College London,
Gower Street, London, UK
^{84}
SBITPLPPC, EPFL, 1015, Lausanne, Switzerland
^{85}
SISSA, Astrophysics Sector, via Bonomea 265, 34136
Trieste,
Italy
^{86}
School of Physics and Astronomy, Cardiff University,
Queens Buildings, The Parade,
Cardiff, CF24 3AA, UK
^{87}
Space Research Institute (IKI), Russian Academy of
Sciences, Profsoyuznaya Str,
84/32, 117997
Moscow,
Russia
^{88}
Space Sciences Laboratory, University of California,
Berkeley, California, USA
^{89}
Special Astrophysical Observatory, Russian Academy of
Sciences, Nizhnij Arkhyz,
Zelenchukskiy region, 369167
KarachaiCherkessian Republic,
Russia
^{90}
Stanford University, Dept of Physics, Varian Physics
Bldg, 382 via Pueblo
Mall, Stanford,
California,
USA
^{91}
SubDepartment of Astrophysics, University of Oxford,
Keble Road, Oxford
OX1 3RH,
UK
^{92}
Theory Division, PHTH, CERN, 1211, Geneva 23, Switzerland
^{93}
UPMC Univ. Paris 06, UMR7095, 98bis Boulevard Arago, 75014
Paris,
France
^{94}
Université de Toulouse, UPSOMP, IRAP, 31028
Toulouse Cedex 4,
France
^{95}
University of Granada, Departamento de Física Teórica y del Cosmos,
Facultad de Ciencias, 1807
Granada,
Spain
^{96}
Warsaw University Observatory, Aleje Ujazdowskie 4, 00478
Warszawa,
Poland
Received: 21 March 2013
Accepted: 18 March 2014
The multifrequency capability of the Planck satellite provides information both on the integrated history of star formation (via the cosmic infrared background, or CIB) and on the distribution of dark matter (via the lensing effect on the cosmic microwave background, or CMB). The conjunction of these two unique probes allows us to measure directly the connection between dark and luminous matter in the high redshift (1 ≤ z ≤ 3) Universe. We use a threepoint statistic optimized to detect the correlation between these two tracers, using lens reconstructions at 100, 143, and 217 GHz, together with CIB measurements at 100–857 GHz. Following a thorough discussion of possible contaminants and a suite of consistency tests, we report the first detection of the correlation between the CIB and CMB lensing. The well matched redshift distribution of these two signals leads to a detection significance with a peak value of 42/19σ (statistical/statistical + systematics) at 545 GHz and a correlation as high as 80% across these two tracers. Our full set of multifrequency measurements (both CIB auto and CIBlensing crossspectra) are consistent with a simple halobased model, with a characteristic mass scale for the halos hosting CIB sources of log_{10}(M/M_{⊙}) = 10.5 ± 0.6. Leveraging the frequency dependence of our signal, we isolate the high redshift contribution to the CIB, and constrain the star formation rate (SFR) density at z ≥ 1. We measure directly the SFR density with around 2σ significance for three redshift bins between z = 1 and 7, thus opening a new window into the study of the formation of stars at early times.
Key words: gravitational lensing: weak / cosmic background radiation / largescale structure of Universe / dark matter / galaxies: star formation
© ESO, 2014
1. Introduction
This paper, one of a set associated with the 2013 release of data from the Planck^{1} mission (Planck Collaboration I 2014), presents a first detection of a strong correlation between the infrared background anisotropies and a lensingderived 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 largescale 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 UVheated dust, enshrouding young stars. The CIB contains much of the energy from processes involved in structure formation. According to current models, the dusty starforming galaxies (DSFGs), which form the CIB have a redshift distribution peaked between z ∼ 1 and z ∼ 2, and tend to live in 10^{11}–10^{13} 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 arcminute 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 crosscorrelation 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 signaltonoise 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 crosscorrelation with measurements of the CIB in the Planck HFI bands to make the first detection of the lensinginfrared 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.
Fig. 1 Redshift and massintegrand for the CIB and CMB lensing potential power spectra at ℓ = 500, calculated using the CIB halo model of Planck Collaboration XVIII (2011), evaluated at 217 GHz. The good match between the redshift and halo mass distributions leads to an expected correlation of up to 80%. The sharper features in the CIB kernel are artefacts from the Béthermin et al. (2012c) model. We note that the low mass, highz behaviour of our calculation is limited by the accuracy of the mass function we use (Tinker & Wetzel 2010). All of our mass integrals use M_{min} = 10^{5} M_{⊙}. 
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 fullsky surveys. Here we give an overview of the HFI mapmaking process, with additional details given in Planck HFI Core Team (2011b); Planck Collaboration II (2014). The data are organized as timeordered 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) decorrelation 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 fullwidth at halfmaximum (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 intercalibration 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 “zodiremoved” 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 “halfring” (HR) maps to refer to maps made using the first and second half of the stable pointing period, “survey” for individual fullsky 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 highlatitude 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 crosscorrelation analysis, since it provides a method to reduce the largescale noise in our measurement, and the 3point 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 21cm 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^{20} 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^{2} regions. If the skewness is larger than a given value then the region is masked. This is an improvement over the usual cutoff 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^{2} maps performed by Land & Slosar (2007) and made available through the LAMBDA website^{3}. 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 brightnesstemperature 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.
Fig. 2 Combined Galactic, pointsource and H i masks used throughout the analysis, with unmasked sky fractions of 16% (red), 30% (red and orange), and 43% (red, orange and light blue). The dark blue area is removed for all of these masks. 
2.2.2. IRIS/IRAS maps
As a consistency test we will use an additional tracer of the CIB that derives from reprocessed IRAS maps at 60 and 100 μm. This new generation of IRAS maps, known as IRIS (MivilleDeschênes & Lagache 2005), benefits from improved zodiacal light subtraction, a calibration and zero level compatible with DIRBE, and an improved destriping procedure. IRAS made two fullsky maps (HCON1 and HCON2), as well as a final map that covers 75% of the sky (HCON3). The three maps had identical processing that included deglitching, checking of the zerolevel stability, visual examination for glitches and artefacts, and zodiacal light removal. The three HCONs were then coadded, taking into account the inhomogeneous sky coverage maps, to generate the average map (HCON0). 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 crosscorrelation 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. Crosscorrelation formalism and implementation
We now describe our statistical formalism and its implementation, with additional technical details given in the appendices. Our analysis consists of crosscorrelating a fullsky 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 , where 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 offdiagonal contributions to the CMB covariance matrix induced by lensing;

3.
subtract a “meanfield bias”, which corrects for known nonlensing 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, , which we use to weight our crosscorrelation 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 crosscorrelation 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 nonnegligible amount of CIB that would be removed this way is not easy, given the uncertainties in the crossfrequency 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 icorrelated dust component. As was performed in PER, we split the H i map into two velocity components: a lowvelocity local component (LC) typical of highlatitude H i emission, and a component of intermediatevelocity clouds (IVC). We found that the inclusion of a highvelocity component makes a negligible difference to the dustcleaned map. Unlike the dedicated highresolution H i observations used in PER and P2013 that only partially cover the sky, here we use the fullsky, 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 crosscorrelation measurement has high signaltonoise 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 dustH i correlation is constant. For the sake of convenience, we use regions of approximate size 13 (52) deg^{2} 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 CMBsubtracted maps, and use the two coefficients per region to assemble a fullsky (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 (CMBunsubtracted) Planck maps.
Fig. 3 Angular crossspectra between the reconstructed lensing map and the temperature map at the six HFI frequencies. The error bars correspond to the scatter within each band. The solid line is the expected result based on the PER model and is not a fit to these data (see Fig. 16 for an adjusted model), although it is already a satisfying model. In each panel we also show in grey the correlation between the lens reconstruction at 143 GHz and the 143 GHz temperature map. This is a simple illustration of the frequency scaling of our measured signal and also the strength of our signal as compared to possible intrafrequency systematic errors. 
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 dustfocused analysis detailed in Planck Collaboration V (2014). However, in the case of our crosscorrelation analysis the dustremoval 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 crossspectrum of removing this emission, and will place limits on the residual Galactic contamination in Sect. 5.3.5.
3.3. Measuring crosscorrelations
To estimate the crosscorrelation between the lensing potential and a tracer t, we calculate (1)Since the CIB has an approximately ℓ^{1} dependence and the lensing potential has an ℓ^{2} dependence, we multiply the crossspectra by ℓ^{3}, 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 meanfield 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 pseudoC_{ℓ} 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 bestfit models we have evaluated the correlation between different bins of the crosscorrelation signal for our nominal binning scheme. We find that the maskinduced 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 crossspectrum (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 crossspectrum we will assume that the noise correlation between bins can be neglected and that the bandpowers are flat.
3.4. Simulating crosscorrelations
In order to validate our measurement pipeline, and to confirm that our estimate of the crossspectrum 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 bestfit CIB model to generate CIB auto and CIBφ crossspectra, 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 crosspower 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 crossspectrum 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
Fig. 4 Temperature maps of size 1 deg^{2} at 545 and 857 GHz stacked on the 20 000 brightest peaks (left column), troughs (centre column) and random map locations (right column). The stacked (averaged) temperature maps is in kelvin. The arrows indicate the lensing deflection angle deduced from the gradient of the bandpassfiltered lensing potential map stacked on the same peaks. The longest arrow corresponds to a deflection of 6.3^{′′}, which is only a fraction of the total deflection angle because of our filtering. This stacking allows us to visualize in real space the lensing of the CMB by the galaxies that generate the CIB. A small and expected offset (around 1^{′}) was corrected when displaying the deflection field. 
We now describe the result of applying our pipeline to our nominal data set, i.e., the lens reconstruction at 143 GHz and the foregroundreduced 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 CIBlensing correlation predicted using the PER CIB model (the HOD parameters of the PER 217 GHz bestfit 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: (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 intrafrequency 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 CIBlensing signal. To get a better visualization of this detection, we show in Fig. 4 the realspace 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 bandpassfilter the lens map between ℓ = 400 and 600 to remove scales larger than our stacked map as well as smallscale noise. We stack a 1 deg^{2} 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 crosspower 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 noisefree 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^{2}). 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 starforming galaxies tracing the largescale 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 autospectra 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, , are calculated using the naive prescription (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, and are the spectra measured using the data, and is the model cross spectrum. This last term provides a negligible contribution, due to the large noise bias on , as we now describe.
The statistical error has two sources, instrumental and astrophysical. The measured autospectra in Eq. (3) contain a signal and noise contribution: . It is informative to split the right hand side of Eq. (3) into four pieces: (4)Here the first term is a signalonly piece, the second is a noiseonly 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: . 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, , 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 noiseonly term (in blue), which is proportional to the temperature noise spectrum. At 353 GHz and above the mixed signalnoise term (orange in Fig. 5) becomes important and is the largest contribution at 545 and 857 GHz at high ℓ.
Fig. 5 Naive analytical estimates of the contribution to the variance as a function of multipole and frequency, as given in Eq. (4). We assume the same bin sizes as in Fig. 3. The different lines are the contribution to the analytical error from: the signal only, (orange); noise only, (dark blue); and the mixed signal and noise terms, (orange) and (red). The total contribution is the solid black line, and the theory spectrum, , is the dashed black line. 
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 3point 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 3point measurement, and many effects that can lead to biases in the autospectrum 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 frequencydependent manner. As discussed in Planck Collaboration VIII (2014), uncertainties in the beam transfer (as well as the fiducial CMB power spectrum ) 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 smallscale hot and cold spots in the CMB. The quadratic estimator, which we use to reconstruct the lensing potential is optimized to measure the lensinginduced 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 “meanfield” 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 meanfield introduced by the application of a Galaxy and pointsource 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 meanfield 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 nonzero monopole, for example, this will correlate strongly with any error in the mask meanfield 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 χ^{2} 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 halfring (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 longterm 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 longterm (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 analoguetodigital 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 longterm gain evolution. Apart from this failure, which we take into account by increasing calibration uncertainties, the reasonable χ^{2}/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.
Fig. 6 Null tests at 545 GHz. Left: difference spectra obtained by nulling the signal in the halfring temperature map before correlating it with our nominal φ reconstruction. Middle: temperature signal nulled using different detectors at 545 GHz. Right: temperature signal nulled using the first and second survey maps. The black error bars correspond to the scatter measured within an ℓ bin, while the coloured bands correspond to the analytical estimate. Except for the survey null test (see text for details), these tests are passed satisfactorily, as illustrated by the quoted χ^{2} and N_{dof}, thus strengthening confidence in the cosmological nature of our signal. 
Fig. 7 Left: difference between the crossspectra measured using the 20% Galactic mask (20% is the unmasked sky fraction) from that measured with our default 40% Galactic mask. Middle: spectra obtained when differencing the 60% and 40% Galaxy mask measurements. For both left and middle panels and all Galactic masks the same point source and H i masks are used, which removes an additional fraction of the sky. Right: difference between the crossspectra calculated with the H icleaned temperature maps from those with no H i cleaning. This crossspectrum is thus the correlation between the H i template and the φ reconstruction. The error bars are calculated in the same way as in Fig. 6. Again, the null tests are passed with an acceptable χ^{2}. 
Fig. 8 Left: difference between crossspectra calculated using the lens reconstruction at 100 GHz with the nominal 143 GHz reconstruction. We see an overall shift, which leads to a high reduced χ^{2}. This shift can be explained by the expected overall normalization uncertainties of the 100 GHz and 143 GHz reconstructions. While this uncertainty is not included in the χ^{2} or the solid bars, it is included later in our analysis in Sect. 6.1. Middle: same as the left panel, but the 217 GHz reconstruction is used instead of the 100 GHz reconstruction. Right: difference between crossspectra when we consider the 143 GHz lens reconstruction calculated with a less restrictive Galaxy mask (that excludes only 20% of the sky) and the nominal reconstruction mask that excludes 40% of the sky. 
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, freefree 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 dustcleaning 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 ∝ ν^{2} 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 crossspectra 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 χ^{2} 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 icorrelated 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 noncleaned crossspectra, we observe that the largescale 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 crossspectrum, 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 crossspectra. From the three studies in this section we conclude that there are no obvious signs of Galactic foreground contamination in our crosscorrelation.
5.3.2. Point source contamination
We now discuss another wellknown 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 shotnoise 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 (5)where is the inversevariance filtered sky map. This estimator is simply the square of the inversevariance filtered sky map, which is a more sensitive probe of point sources than the standard lensing estimator.
In Fig. 10 we plot the crossspectrum of measured at 143 GHz and for the full set of HFI channels. This crossspectrum is probing the same point source contributions that could bias our estimates of , although with a greater signaltonoise 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 crossspectra is given by (6)
where , is the unlensed CMB spectrum, is the spectrum of , is as defined in Okamoto & Hu (2003), and we have used the fact that for our inversevariance filtering, . We have calculated this contribution using our measured and subtracted it from the data points of Fig. 10.
We can consider the effect of shot noise on this crossspectrum. With the shotnoise bispectrum defined by (7)the bias to the plotted crossspectrum is given as (8)This bias is plotted for bestfit values of ⟨ S^{3} ⟩ as the black lines in Fig. 10. To minimize systematic effects that might bias the ⟨S^{3}⟩ estimator, we have estimated S^{3} from the spectra of Fig. 10 for multipoles between ℓ = 500 and 2000. The fitted ⟨S^{3}⟩ 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 shotnoise bias. These conclusions have been verified using less restrictive point source masks that cover fewer sources.
With estimates of S^{3} in hand, we estimate a corresponding bias to , given by (9)where , with defined in Planck Collaboration VIII (2014). We show this contribution later in Fig. 12 as the dotted line. While nonzero, we see that the point source shotnoise 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 crossspectra).
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 socalled thermal SunyaevZeldovich (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 largescale matter potential that is traced by CMB lensing, we expect an SZinduced 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).
Point source estimator.
The frequency dependence of the SZ signal in our map depends on the detector bandpasses and is (10)where h(ν) is the detector bandpass and g(ν) is the SZ frequency dependence, which in the nonrelativistic limit is g(ν) = x(e^{x} + 1) / (e^{x}−1) − 4, with x = hν/k_{B}T_{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.
Fig. 9 Frequency spectrum of our crossspectra averaged within ℓ bins (black points with associated error bars). The light shaded regions correspond to the HFI frequency bands. The solid black curve corresponds to the bestfit CIB assuming a Gispert et al. (2000) spectrum, while the dotdashed line assumes a Fixsen et al. (1998) spectrum. The dashed black line corresponds to the bestfit model when allowing for an SZ component in addition to the Gispert et al. (2000) CIB shape. The blue dots correspond to the associated absolute value of the bestfit SZ component. We conclude from this plot that the SZ effect is not an important contaminant. 
To measure a contribution from the SZlensing 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 χ^{2} 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 bestfit SZlensing and CIBlensing spectra. For the CIBonly fit with the Gispert et al. (2000) frequency dependence we find a relatively poor fit in the lowest ℓ bin, χ^{2} (N_{d.o.f.}) = 15.5 (5), but an improved fit in the higher ℓ bin, χ^{2}(N_{d.o.f.}) = 4.15 (5). Including the SZ component gives an improvement in the χ^{2} 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 χ^{2}(N_{d.o.f.}) = 2.25 (5) and 5.49 (5) in the low and highℓ bins, respectively. Overall, the improvement in the χ^{2}/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 SZlensing 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 crossspectrum at 217–857 GHz we assume that the measured crossspectrum 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 crossspectrum at 143 GHz using only the SZlensing correlation requires an amplitude of (2.4 ± 1.6) times our calculated SZlensing crossspectrum. 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 SachsWolfe (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 ISWinduced 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 crossspectrum is consistent with zero.
Fig. 10 Results from the point source contamination estimator of Eq. (5). The bestfit crossspectra associated with shot noise are plotted in black. We note that the signaltonoise ratio at 545 and 857 GHz is particularly low. The grey line is a prediction for the bias from the CMB lensinginfrared correlation, and has been subtracted from the spectra (plotted as black points). We see that with the subtraction of the bias from CMB lensing, the measured bispectrumrelated spectrum is generally consistent either with zero or with the shape expected for shot noise. 
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 3point 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 nonGaussian 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 bestfit curve in Fig. 10, since the bestfit model that we used assumes only a Poissonian shotnoise contribution. We do not see any significant deviation in Fig. 10. Still, in order to isolate this effect we create crossspectra 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 crosscorrelation between the 100 GHz reconstruction and the 100 GHz temperature map does not show any large difference with the crossspectrum 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 crosscorrelation between the raw 545 GHz lens reconstruction and the 545 GHz temperature map. This crossspectrum is shown in Fig. 11 for three different Galaxy masks. The line shows the point source shotnoise template derived in Sect. 5.3.2, fit to the crossspectrum 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 ℓ.
Fig. 11 Crossspectrum of the 545 GHz lens reconstruction correlated with the 545 GHz temperature map, using different Galactic masks. The legend gives the visible sky fractions. The solid line represents the analytic unclustered shotnoise contribution fit to the f_{sky} = 0.09 points above ℓ = 1300. 
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 nonGalactic 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 CIBlensing 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 bandpowers in Table 2, along with our statistical and systematic errors. These bandpowers 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 bestfit CIB component with either a Fixsen et al. (1998) – dotdashed 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 crosscorrelation 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 highz 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 CIBlensing crossspectra measured in this paper and the CIB autospectra obtained from P2013. We use the CIBlensing crossspectra 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 crossspectra alone. As can be seen in P2013, the crossspectra of the CIB at different frequencies provide powerful constraints on the CIB emissivity.
Fig. 12 Foreground components at each frequency. The data points and error bars show our results. The dashed line is an estimated upper limit on the magnitude of the SZ contamination derived in Sect. 5.3.3. We show the absolute value of this contribution, which is negative at frequencies less than 217 GHz. The dotdashed line is the extragalactic point source contribution, with an amplitude measured from our data, as derived in Sect. 5.3.2. Again we show the absolute value, with the signal being negative below ℓ ∼ 1200. The oscillating solid line corresponds to the calculated ISW contamination. 
Crossspectrum detection bandpowers.
We use the loglikelihood, (11)where Ĉ and 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 crossspectrum as being uncorrelated. The small (≃2%) maskinduced modecoupling between neighbouring bins supports this approximation. However, calibration and beam errors (which are correlated between the auto and crossspectra 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 crosscorrelation, 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), (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
Fig. 13 Crosscorrelation coefficients calculated from the model φ spectrum and bestfit halo model at each frequency. The CIB is a spectacular tracer of CMB lensing, and viceversa. The data points represent the measured crosscorrelation divided by the bestfit autopower spectra models at 545 GHz. 
Fig. 14 Marginalized 2D distributions of z_{c} and σ_{z} for the linear bias model, fit to all frequencies simultaneously. The orange dots indicate the parameter values at the minimum χ^{2}. 
As a first pass at interpreting our measurement, we will consider a redshiftindependent linear bias model with a simple parametric SED. This model was found to provide a reasonable fit to the autospectra 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 (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: (14)Here b is the DSFG bias that we assume to be redshiftindependent, 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_{0} 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: (15)where we use a modified blackbody frequency dependence (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 crossspectra this approach would be equivalent to marginalizing over a frequencyindependent 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 (17)Using this equation and the measured FIRAS mean and uncertainty we calculate a χ^{2} value and add this to the χ^{2} in Eq. (11). Since there are no FIRAS constraints at 100 and 143 GHz, as well as no CIB autospectra measurements, and noisier crossspectra measurements at these frequencies, our constraints for the 100 and 143 GHz redshift parameters are weaker than for the other parameters.
Fig. 15 Measured crossspectra with the bestfit mean emissivity j reconstruction model fit to both the CIB auto and CIBlensing crossspectra (solid coloured), and the bestfit linear bias model (dashed coloured). The χ^{2} values quoted in each panel are the contribution to the global χ^{2} from the data in the panel for the halo model (with N_{d.o.f.} in brackets), and loosely indicate the goodness of fit (see text for details). The one and twohalo contributions are shown as the dotdashed and dashed thin black lines, respectively. A light dashed black horizontal line indicates the zero level. 
Fig. 16 P2013 (i.e., Planck Collaboration IX 2014) autospectra with the bestfit mean emissivity j reconstruction model fit for the CIB auto and CIBlensing crossspectra (solid coloured). The χ^{2} values are defined as in Fig. 15. The one and twohalo contributions are shown as the dotdashed and dashed thin black lines, respectively, while shot noise is the dotted thin black line. 
The linear bias model considers only linear clustering, and so when fitting the autospectra we restrict ourselves to ℓ < 500, where nonlinear contributions are negligible. Because we do not consider the highℓ modes, we also neglect the shotnoise contribution to the autospectra. The bestfit model is shown as the coloured dashed lines in Fig. 15, with χ^{2} values of 13.4, 16.8, 25.2, 21.8, 9.1, and 9.4, if we break up the χ^{2} contribution per frequency from 100 to 857 GHz, leading to an overall χ^{2} 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 twodimensional 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 χ^{2} 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 highredshift 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: (18)Here is the CIB emissivity at redshift z with mean , is the number density of DSFGs with mean , and is the DSFG overdensity, with power spectrum ⟨δ_{g}(k,z) δ_{g}(k^{′},z)^{∗}⟩ = (2π)^{3}δ(k − k^{′}) P_{gg}(k,z). We calculate this power spectrum, including the constituent 1 and 2halo 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 autospectra have a dependence, and so if the true value of has a strong z dependence within a bin then the bestfit emissivity in the bin will be difficult to interpret; the bestfit 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 bestfit 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 crossspectra are (Song et al. 2003): (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 autospectra and for the lensing crossspectra.
Overall, the halobased 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 crossspectra have a total of 120 ℓ bins, with four additional FIRAS data points. Solving for the likelihood described in Sect. 6.1, gives the bestfit models shown in Figs. 15 and 16 as the solid coloured lines. The χ^{2} values in each panel are the contribution to the total χ^{2} from the data within the panel. The combined reduced χ^{2} 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 bestfit parameters are log _{10}(M_{min}/M_{⊙}) = 12.2 and log _{10}(M_{sat}/M_{⊙}) = 12.8, which gives M_{sat}/M_{min} = 3.8. The mean parameter values are log _{10}(M_{min}/M_{⊙}) = 10.5 ± 0.6 and log _{10}(M_{sat}/M_{⊙}) = 10.8 ± 0.7. The bestfit 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: (20)Here S_{ν} is the flux density, and d^{2}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: (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) (22)The contribution to the infrared luminosity density from a given population is (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 (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, (25)Using Eq. (25) we can find the coefficients for each of the redshift bins and frequencies, which are presented in Table 3.
Fig. 17 Reconstructed mean emissivity, , for each frequency as a function of redshift. The solid line at low z and the dashed line at higher z correspond to the B11 model. The B11 emissivity model at z> 1 is not used, and is shown only for reference. The black error bars correspond to the 68% C.L. region, while the coloured shading displays the full posterior distribution. 
Fig. 18 Marginalized 1D distribution of the emissivity in the high redshift bin at 353 GHz with (black line) or without (blue line) including the CIBlensing correlation. Its inclusion helps to constrain the emissivity at high redshift, transforming an upper limit into a detection. 
Reconstructed emissivity as a function of redshift and associated star formation rate.
Fig. 19 Correlation between the lensing potential and the IRIS map at 100 μm, using our nominal lens reconstruction. We clearly see a correlation and estimate the significance to be 9σ, ignoring possible systematic effects. The solid line represents a simple reasonable prediction for this signal. 
6.4. Discussion and outlook
In the previous section we described a model that simultaneously fits the CIB autospectra and the CIBlensing crossspectra, 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 bestfit 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 _{10}(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 modeldependent. 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 crossspectrum measurement is worth quantifying. As an example, we show in Fig. 18 the highest redshift emissivity bin in the 353 GHz band. Adding the CIBlensing crossspectrum information tightens the constraint on the highredshift 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 crosscorrelation more sensitive to the highredshift CIB signal than the CIB autospectrum, 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 nonzero 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 nonGaussian 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 signaltonoise 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 CIBlensing 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 massindependent luminosity function. In addition the use of mapbased 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 datasets. 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 lowredshift 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 fullsky, multifrequency, 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.
Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states (in particular the lead countries France and Italy), with contributions from NASA (USA) and telescope reflectors provided by a collaboration between ESA and a scientific consortium led and funded by Denmark.
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/INSUIN2P3INP (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
 Addison, G. E.,Dunkley, J., &Bond, J. R. 2013, MNRAS, 436, 1896 [NASA ADS] [CrossRef] [Google Scholar]
 Amblard, A.,Cooray, A.,Serra, P., et al. 2011, Nature, 470, 510 [NASA ADS] [CrossRef] [Google Scholar]
 Arnal, E. M.,Bajaja, E.,Larrarte, J. J.,Morras, R., &Pöppel, W. G. L. 2000, A&AS, 142, 35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Babich, D., &Pierpaoli, E. 2008, Phys. Rev., 77, 3011 [NASA ADS] [Google Scholar]
 Bajaja, E.,Arnal, E. M.,Larrarte, J. J., et al. 2005, A&A, 440, 767 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bersanelli, M.,Mandolesi, N.,Butler, R. C., et al. 2010, A&A, 520, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M., Dole, H., Lagache, G., Le Borgne, D., &Pénin, A. 2011, A&A, 529, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M.,Daddi, E.,Magdis, G., et al. 2012a, ApJ, 757, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Béthermin, M., Floc’h, E. L.,Ilbert, O., et al. 2012b, A&A, 542, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M.,Doré, O., &Lagache, G. 2012c, A&A, 537, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blanchard, A., &Schneider, J. 1987, A&A, 184, 1 [NASA ADS] [Google Scholar]
 Blandford, R. D., &Jaroszynski, M. 1981, ApJ, 246, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Bleem, L. E., van Engelen, A.,Holder, G. P., et al. 2012, ApJ, 753, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Boulanger, F., &Perault, M. 1988, ApJ, 330, 964 [NASA ADS] [CrossRef] [Google Scholar]
 Boulanger, F.,Abergel, A.,Bernard, J. P., et al. 1996, A&A, 312, 256 [NASA ADS] [Google Scholar]
 Crawford, T., Schaffer, K., Bhattacharya, S., et al. 2013, ApJ, submitted [arXiv:1303.3535] [Google Scholar]
 Das, S.,Sherwin, B. D.,Aguirre, P., et al. 2011, Phys. Rev. Lett., 107, 1301 [NASA ADS] [Google Scholar]
 De Bernardis, F., &Cooray, A. 2012, ApJ, 760, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Dunkley, J.,Hlozek, R.,Sievers, J., et al. 2011, ApJ, 739, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Finkbeiner, D. P.,Davis, M., &Schlegel, D. J. 1999, ApJ, 524, 867 [NASA ADS] [CrossRef] [Google Scholar]
 Fixsen, D. J.,Dwek, E.,Mather, J. C.,Bennett, C. L., &Shafer, R. A. 1998, ApJ, 508, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Gispert, R.,Lagache, G., &Puget, J. L. 2000, A&A, 360, 1 [NASA ADS] [Google Scholar]
 Goldberg, D. M., &Spergel, D. N. 1999, Phys. Rev., D, 59, 3002 [Google Scholar]
 Górski, K. M.,Hivon, E.,Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Haiman, Z., &Knox, L. 2000, ApJ, 530, 124 [NASA ADS] [CrossRef] [Google Scholar]
 Hall, N. R.,Keisler, R.,Knox, L., et al. 2010, ApJ, 718, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Hand, N.,Addison, G. E.,Aubourg, E., et al. 2012, Phys. Rev. Lett., 109, 041101 [NASA ADS] [CrossRef] [Google Scholar]
 Hartmann, D., & Burton, W. B. 1997, Atlas of Galactic Neutral Hydrogen (Cambridge University Press) [Google Scholar]
 Hauser, M. G.,Arendt, R. G.,Kelsall, T., et al. 1998, ApJ, 508, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Hikage, C.,Mandelbaum, R.,Takada, M., &Spergel, D. N. 2013, MNRAS, 435, 2345 [NASA ADS] [CrossRef] [Google Scholar]
 Hirata, C. M.,Ho, S.,Padmanabhan, N.,Seljak, U., &Bahcall, N. A. 2008, Phys. Rev., D, 78, 3520 [Google Scholar]
 Hivon, E.,Gorski, K.,Netterfield, C., et al. 2002, ApJ, 567, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Hopkins, A. M., &Beacom, J. F. 2006, ApJ, 651, 142 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kalberla, P. M.,Burton, W.,Hartmann, D., et al. 2005, A&A, 440, 775 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kennicutt, R. C. 1998, ARA&A, 36, 189 [Google Scholar]
 Knox, L.,Christensen, N., &Skordis, C. 2001, ApJ, 563, L95 [NASA ADS] [CrossRef] [Google Scholar]
 Lacasa, F.,Aghanim, N.,Kunz, M., &Frommert, M. 2012, MNRAS, 421, 1982 [NASA ADS] [CrossRef] [Google Scholar]
 Lagache, G.,Abergel, A.,Boulanger, F., &Puget, J. L. 1998, A&A, 333, 709 [NASA ADS] [Google Scholar]
 Lagache, G.,Haffner, L. M.,Reynolds, R. J., &Tufte, S. L. 2000, A&A, 354, 247 [Google Scholar]
 Lagache, G.,Bavouzet, N.,FernandezConde, N., et al. 2007, ApJ, 665, L89 [NASA ADS] [CrossRef] [Google Scholar]
 Lamarre, J.,Puget, J.,Ade, P. A. R., et al. 2010, A&A, 520, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Land, K., &Slosar, A. 2007, Phys. Rev. D, 76, 7301 [Google Scholar]
 Leahy, J. P., Bersanelli, M., D’Arcangelo, O., et al. 2010, A&A, 520, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leauthaud, A.,George, M. R.,Behroozi, P. S., et al. 2012, ApJ, 746, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A., &Bridle, S. 2002, Phys. Rev., D, 66, 3511 [Google Scholar]
 Lewis, A., &Challinor, A. 2006, Phys. Rept., 429, 1 [Google Scholar]
 Lewis, A.,Challinor, A., &Lasenby, A. 2000, ApJ, 538, 473 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A.,Challinor, A., &Hanson, D. 2011, JCAP, 1103, 018 [NASA ADS] [CrossRef] [Google Scholar]
 Magdis, G. E.,Daddi, E.,Bethermin, M., et al. 2012, ApJ, 760, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Mandolesi, N.,Bersanelli, M.,Butler, R. C., et al. 2010, A&A, 520, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mennella, A.,Butler, R. C.,Curto, A., et al. 2011, A&A, 536, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MivilleDeschênes, M.A. &Lagache, G. 2005, ApJS, 157, 302 [NASA ADS] [CrossRef] [Google Scholar]
 Nozawa, S.,Itoh, N.,Kawana, Y., &Kohyama, Y. 2000, ApJ, 536, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Okamoto, T., &Hu, W. 2003, Phys. Rev., D, 67, 3002 [Google Scholar]
 Peebles, P. J. E. 1980, The largescale structure of the universe (Princeton university Press, Princeton, N.J.) [Google Scholar]
 Pénin, A.,Doré, O.,Lagache, G., &Béthermin, M. 2012, A&A, 537, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011, The Explanatory Supplement to the Planck Early Release Compact Source Catalogue (ESA) [Google Scholar]
 Planck Collaboration I. 2011, A&A, 536, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2011, A&A, 536, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2011, A&A, 536, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2011, A&A, 536, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2011, A&A, 536, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIV. 2011, A&A, 536, A24 [Google Scholar]
 Planck Collaboration I. 2014, A&A, 571, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2014, A&A, 571, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2014, A&A, 571, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2014, A&A, 571, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2014, A&A, 571, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2014, A&A, 571, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2014, A&A, 571, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2014, A&A, 571, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2014, A&A, 571, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2014, A&A, 571, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2014, A&A, 571, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2014, A&A, 571, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIV. 2014, A&A, 571, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XV. 2014, A&A, 571, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVII. 2014, A&A, 571, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2014, A&A, 571, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2014, A&A, 571, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XX. 2014, A&A, 571, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXI. 2014, A&A, 571, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2014, A&A, 571, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIII. 2014, A&A, 571, A23 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Planck Collaboration XXIV. 2014, A&A, 571, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXV. 2014, A&A, 571, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2014, A&A, 571, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2014, A&A, 571, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVIII. 2014, A&A, 571, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIX. 2014, A&A, 571, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXXI. 2014, A&A, 571, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck HFI Core Team 2011a, A&A, 536, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck HFI Core Team 2011b, A&A, 536, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puget, J. L.,Abergel, A.,Bernard, J. P., et al. 1996, A&A, 308, L5 [NASA ADS] [Google Scholar]
 Reichardt, C.,Shaw, L.,Zahn, O., et al. 2012, ApJ, 755, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Rosset, C.,Tristram, M.,Ponthieu, N., et al. 2010, A&A, 520, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sachs, R., &Wolfe, A. 1967, ApJ, 147, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Seljak, U., &Zaldarriaga, M. 1999, Phys. Rev., D, 60, 3504 [Google Scholar]
 Shang, C.,Haiman, Z.,Knox, L., &Oh, S. P. 2012, MNRAS, 421, 2832 [NASA ADS] [CrossRef] [Google Scholar]
 Sherwin, B. D.,Das, S.,Hajian, A., et al. 2012, Phys. Rev. D, 86, 3006 [Google Scholar]
 Smith, K. M.,Zahn, O., &Doré, O. 2007, Phys. Rev., D, 76, 3510 [Google Scholar]
 Smith, K. M.,Cooray, A.,Das, S., et al. 2009, AIP Conf. Proc., 1141, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Song, Y.S.,Cooray, A.,Knox, L., &Zaldarriaga, M. 2003, ApJ, 590, 664 [NASA ADS] [CrossRef] [Google Scholar]
 Sunyaev, R. A., &Zeldovich, Y. B. 1970, Ap&SS, 7, 3 [NASA ADS] [Google Scholar]
 Sunyaev, R. A., &Zeldovich, I. B. 1980, MNRAS, 190, 413 [NASA ADS] [CrossRef] [Google Scholar]
 Tauber, J. A.,Mandolesi, N.,Puget, J., et al. 2010, A&A, 520, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tinker, J. L., &Wetzel, A. R. 2010, ApJ, 719, 88 [NASA ADS] [CrossRef] [Google Scholar]
 van Engelen, A.,Keisler, R.,Zahn, O., et al. 2012, ApJ, 756, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Viero, M. P.,Ade, P. A. R.,Bock, J. J., et al. 2009, ApJ, 707, 1766 [NASA ADS] [CrossRef] [Google Scholar]
 Viero, M.,Wang, L.,Zemcov, M., et al. 2013, ApJ, 772, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Wilson, M. J.,Sherwin, B. D.,Hill, J. C., et al. 2012, Phys. Rev. D, 86, 2005 [Google Scholar]
 Zacchei, A.,Maino, D.,Baccigalupi, C., et al. 2011, A&A, 536, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Statistical errors
Fig. A.1 Ratio of various error estimation procedures to the errors obtained with the databased analytical estimate. At each frequency the numerator is given by: (i) the scatter within an ℓ bin in simulations (solid black line); (ii) the scatter within an ℓ bin in the data (solid dashed black line); (iii) the scatter of bins across simulated realizations (solid coloured line); (iv) the analytical errors calculated from the simulations (dashed coloured line); and (v) the scatter across realizations for the crosscorrelation between the simulated temperature map and the lensing potential reconstructed from the data (coloured dotdashed line). The grey envelope is the precision of the simulated errors expected from 100 simulations (shown as a spread around unity). 
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, , , and the model crossspectrum .

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 datadetermined crossspectrum.

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 crossspectra 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 crossspectrum.

6.
The error in the crossspectrum 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 crosscheck.
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 databased 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 databased analytical estimates at all frequencies, we use the databased analytical estimates as the basis for our statistical errors. We could alternatively scale the scatterdetermined 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 nonGaussian 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 crossspectrum of the reconstructed lensing potential in simulated maps with the datameasured temperature maps (Method 6, coloured dasheddotted 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
Reconstructed emissivity as a function of redshift and associated star formation rate.
All Figures
Fig. 1 Redshift and massintegrand for the CIB and CMB lensing potential power spectra at ℓ = 500, calculated using the CIB halo model of Planck Collaboration XVIII (2011), evaluated at 217 GHz. The good match between the redshift and halo mass distributions leads to an expected correlation of up to 80%. The sharper features in the CIB kernel are artefacts from the Béthermin et al. (2012c) model. We note that the low mass, highz behaviour of our calculation is limited by the accuracy of the mass function we use (Tinker & Wetzel 2010). All of our mass integrals use M_{min} = 10^{5} M_{⊙}. 

In the text 
Fig. 2 Combined Galactic, pointsource and H i masks used throughout the analysis, with unmasked sky fractions of 16% (red), 30% (red and orange), and 43% (red, orange and light blue). The dark blue area is removed for all of these masks. 

In the text 
Fig. 3 Angular crossspectra between the reconstructed lensing map and the temperature map at the six HFI frequencies. The error bars correspond to the scatter within each band. The solid line is the expected result based on the PER model and is not a fit to these data (see Fig. 16 for an adjusted model), although it is already a satisfying model. In each panel we also show in grey the correlation between the lens reconstruction at 143 GHz and the 143 GHz temperature map. This is a simple illustration of the frequency scaling of our measured signal and also the strength of our signal as compared to possible intrafrequency systematic errors. 

In the text 
Fig. 4 Temperature maps of size 1 deg^{2} at 545 and 857 GHz stacked on the 20 000 brightest peaks (left column), troughs (centre column) and random map locations (right column). The stacked (averaged) temperature maps is in kelvin. The arrows indicate the lensing deflection angle deduced from the gradient of the bandpassfiltered lensing potential map stacked on the same peaks. The longest arrow corresponds to a deflection of 6.3^{′′}, which is only a fraction of the total deflection angle because of our filtering. This stacking allows us to visualize in real space the lensing of the CMB by the galaxies that generate the CIB. A small and expected offset (around 1^{′}) was corrected when displaying the deflection field. 

In the text 
Fig. 5 Naive analytical estimates of the contribution to the variance as a function of multipole and frequency, as given in Eq. (4). We assume the same bin sizes as in Fig. 3. The different lines are the contribution to the analytical error from: the signal only, (orange); noise only, (dark blue); and the mixed signal and noise terms, (orange) and (red). The total contribution is the solid black line, and the theory spectrum, , is the dashed black line. 

In the text 
Fig. 6 Null tests at 545 GHz. Left: difference spectra obtained by nulling the signal in the halfring temperature map before correlating it with our nominal φ reconstruction. Middle: temperature signal nulled using different detectors at 545 GHz. Right: temperature signal nulled using the first and second survey maps. The black error bars correspond to the scatter measured within an ℓ bin, while the coloured bands correspond to the analytical estimate. Except for the survey null test (see text for details), these tests are passed satisfactorily, as illustrated by the quoted χ^{2} and N_{dof}, thus strengthening confidence in the cosmological nature of our signal. 

In the text 
Fig. 7 Left: difference between the crossspectra measured using the 20% Galactic mask (20% is the unmasked sky fraction) from that measured with our default 40% Galactic mask. Middle: spectra obtained when differencing the 60% and 40% Galaxy mask measurements. For both left and middle panels and all Galactic masks the same point source and H i masks are used, which removes an additional fraction of the sky. Right: difference between the crossspectra calculated with the H icleaned temperature maps from those with no H i cleaning. This crossspectrum is thus the correlation between the H i template and the φ reconstruction. The error bars are calculated in the same way as in Fig. 6. Again, the null tests are passed with an acceptable χ^{2}. 

In the text 
Fig. 8 Left: difference between crossspectra calculated using the lens reconstruction at 100 GHz with the nominal 143 GHz reconstruction. We see an overall shift, which leads to a high reduced χ^{2}. This shift can be explained by the expected overall normalization uncertainties of the 100 GHz and 143 GHz reconstructions. While this uncertainty is not included in the χ^{2} or the solid bars, it is included later in our analysis in Sect. 6.1. Middle: same as the left panel, but the 217 GHz reconstruction is used instead of the 100 GHz reconstruction. Right: difference between crossspectra when we consider the 143 GHz lens reconstruction calculated with a less restrictive Galaxy mask (that excludes only 20% of the sky) and the nominal reconstruction mask that excludes 40% of the sky. 

In the text 
Fig. 9 Frequency spectrum of our crossspectra averaged within ℓ bins (black points with associated error bars). The light shaded regions correspond to the HFI frequency bands. The solid black curve corresponds to the bestfit CIB assuming a Gispert et al. (2000) spectrum, while the dotdashed line assumes a Fixsen et al. (1998) spectrum. The dashed black line corresponds to the bestfit model when allowing for an SZ component in addition to the Gispert et al. (2000) CIB shape. The blue dots correspond to the associated absolute value of the bestfit SZ component. We conclude from this plot that the SZ effect is not an important contaminant. 

In the text 
Fig. 10 Results from the point source contamination estimator of Eq. (5). The bestfit crossspectra associated with shot noise are plotted in black. We note that the signaltonoise ratio at 545 and 857 GHz is particularly low. The grey line is a prediction for the bias from the CMB lensinginfrared correlation, and has been subtracted from the spectra (plotted as black points). We see that with the subtraction of the bias from CMB lensing, the measured bispectrumrelated spectrum is generally consistent either with zero or with the shape expected for shot noise. 

In the text 
Fig. 11 Crossspectrum of the 545 GHz lens reconstruction correlated with the 545 GHz temperature map, using different Galactic masks. The legend gives the visible sky fractions. The solid line represents the analytic unclustered shotnoise contribution fit to the f_{sky} = 0.09 points above ℓ = 1300. 

In the text 
Fig. 12 Foreground components at each frequency. The data points and error bars show our results. The dashed line is an estimated upper limit on the magnitude of the SZ contamination derived in Sect. 5.3.3. We show the absolute value of this contribution, which is negative at frequencies less than 217 GHz. The dotdashed line is the extragalactic point source contribution, with an amplitude measured from our data, as derived in Sect. 5.3.2. Again we show the absolute value, with the signal being negative below ℓ ∼ 1200. The oscillating solid line corresponds to the calculated ISW contamination. 

In the text 
Fig. 13 Crosscorrelation coefficients calculated from the model φ spectrum and bestfit halo model at each frequency. The CIB is a spectacular tracer of CMB lensing, and viceversa. The data points represent the measured crosscorrelation divided by the bestfit autopower spectra models at 545 GHz. 

In the text 
Fig. 14 Marginalized 2D distributions of z_{c} and σ_{z} for the linear bias model, fit to all frequencies simultaneously. The orange dots indicate the parameter values at the minimum χ^{2}. 

In the text 
Fig. 15 Measured crossspectra with the bestfit mean emissivity j reconstruction model fit to both the CIB auto and CIBlensing crossspectra (solid coloured), and the bestfit linear bias model (dashed coloured). The χ^{2} values quoted in each panel are the contribution to the global χ^{2} from the data in the panel for the halo model (with N_{d.o.f.} in brackets), and loosely indicate the goodness of fit (see text for details). The one and twohalo contributions are shown as the dotdashed and dashed thin black lines, respectively. A light dashed black horizontal line indicates the zero level. 

In the text 
Fig. 16 P2013 (i.e., Planck Collaboration IX 2014) autospectra with the bestfit mean emissivity j reconstruction model fit for the CIB auto and CIBlensing crossspectra (solid coloured). The χ^{2} values are defined as in Fig. 15. The one and twohalo contributions are shown as the dotdashed and dashed thin black lines, respectively, while shot noise is the dotted thin black line. 

In the text 
Fig. 17 Reconstructed mean emissivity, , for each frequency as a function of redshift. The solid line at low z and the dashed line at higher z correspond to the B11 model. The B11 emissivity model at z> 1 is not used, and is shown only for reference. The black error bars correspond to the 68% C.L. region, while the coloured shading displays the full posterior distribution. 

In the text 
Fig. 18 Marginalized 1D distribution of the emissivity in the high redshift bin at 353 GHz with (black line) or without (blue line) including the CIBlensing correlation. Its inclusion helps to constrain the emissivity at high redshift, transforming an upper limit into a detection. 

In the text 
Fig. 19 Correlation between the lensing potential and the IRIS map at 100 μm, using our nominal lens reconstruction. We clearly see a correlation and estimate the significance to be 9σ, ignoring possible systematic effects. The solid line represents a simple reasonable prediction for this signal. 

In the text 
Fig. A.1 Ratio of various error estimation procedures to the errors obtained with the databased analytical estimate. At each frequency the numerator is given by: (i) the scatter within an ℓ bin in simulations (solid black line); (ii) the scatter within an ℓ bin in the data (solid dashed black line); (iii) the scatter of bins across simulated realizations (solid coloured line); (iv) the analytical errors calculated from the simulations (dashed coloured line); and (v) the scatter across realizations for the crosscorrelation between the simulated temperature map and the lensing potential reconstructed from the data (coloured dotdashed line). The grey envelope is the precision of the simulated errors expected from 100 simulations (shown as a spread around unity). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.