Issue 
A&A
Volume 536, December 2011
Planck early results



Article Number  A18  
Number of page(s)  30  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201116461  
Published online  01 December 2011 
Planck early results. XVIII. The power spectrum of cosmic infrared background anisotropies^{⋆}
^{1}
Aalto University Metsähovi Radio Observatory, Metsähovintie 114, 02540 Kylmälä, Finland
^{2}
Agenzia Spaziale Italiana Science Data Center, c/o ESRIN, via Galileo Galilei, Frascati, Italy
^{3}
Astroparticule et Cosmologie, CNRS (UMR7164), Université Denis Diderot Paris 7, Bâtiment Condorcet, 10 rue A. Domon et Léonie Duquet, Paris, France
^{4}
Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{5}
Atacama Large Millimeter/submillimeter Array, ALMA Santiago Central Offices, Alonso de Cordova 3107, Vitacura, Casilla 763 0355, Santiago, Chile
^{6}
Australia Telescope National Facility, CSIRO, P.O. Box 76, Epping, NSW 1710, Australia
^{7}
CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada
^{8}
CNRS, IRAP, 9 Av. Colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
^{9}
California Institute of Technology, Pasadena, California, USA
^{10} Centre of Mathematics for Applications, University of Oslo, Blindern, Oslo, Norway
^{11}
DAMTP, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, UK
^{12}
DSM/Irfu/SPP, CEASaclay, 91191 GifsurYvette Cedex, France
^{13}
DTU Space, National Space Institute, Juliane Mariesvej 30, Copenhagen, Denmark
^{14}
Departamento de Física, Universidad de Oviedo, Avda. Calvo Sotelo s/n, Oviedo, Spain
^{15}
Department of Astronomy and Astrophysics, University of Toronto, 50 Saint George Street, Toronto, Ontario, Canada
^{16}
Department of Physics & Astronomy, University of BritishColumbia, 6224 Agricultural Road, Vancouver, British Columbia, Canada
^{17}
Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK
^{18}
Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, Helsinki, Finland
^{19}
Department of Physics, Princeton University, Princeton, New Jersey, USA
^{20}
Department of Physics, Purdue University, 525 Northwestern Avenue, West Lafayette, Indiana, USA
^{21}
Department of Physics, University of California, Berkeley, California, USA
^{22}
Department of Physics, University of California, One Shields Avenue, Davis, California, USA
^{23}
Department of Physics, University of California, Santa Barbara, California, USA
^{24}
Department of Physics, University of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, Illinois, USA
^{25}
Dipartimento di Fisica G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131 Padova, Italy
^{26}
Dipartimento di Fisica, Università La Sapienza, P. le A. Moro 2, Roma, Italy
^{27}
Dipartimento di Fisica, Università degli Studi di Milano, via Celoria, 16, Milano, Italy
^{28}
Dipartimento di Fisica, Università degli Studi di Trieste, via A. Valerio 2, Trieste, Italy
^{29}
Dipartimento di Fisica, Università di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{30}
Dipartimento di Fisica, Università di Roma Tor Vergata, via della Ricerca Scientifica, 1, Roma, Italy
^{31}
Discovery Center, Niels Bohr Institute, Blegdamsvej 17, Copenhagen, Denmark
^{32}
Dpto. Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{33}
European Southern Observatory, ESO Vitacura, Alonso de Cordova 3107, Vitacura, Casilla 19001, Santiago, Chile
^{34}
European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, Madrid, Spain
^{35}
European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{36}
Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland
^{37}
INAF – Osservatorio Astrofisico di Catania, via S. Sofia 78, Catania, Italy
^{38}
INAF – Osservatorio Astronomico di Padova, vicolo dell’Osservatorio 5, Padova, Italy
^{39}
INAF – Osservatorio Astronomico di Roma, via di Frascati 33, Monte Porzio Catone, Italy
^{40}
INAF – Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, Trieste, Italy
^{41}
INAF/IASF Bologna, via Gobetti 101, Bologna, Italy
^{42}
INAF/IASF Milano, via E. Bassini 15, Milano, Italy
^{43}
INRIA, Laboratoire de Recherche en Informatique, Université ParisSud 11, Bâtiment 490, 91405 Orsay Cedex, France
^{44}
IPAG: Institut de Planétologie et d’Astrophysique de Grenoble, Université Joseph Fourier, Grenoble 1/CNRSINSU, UMR 5274, 38041 Grenoble, France
^{45}
Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK
^{46}
Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA
^{47}
Institut Néel, CNRS, Université Joseph Fourier Grenoble I, 25 rue des Martyrs, Grenoble, France
^{48}
Institut d’Astrophysique Spatiale, CNRS (UMR8617) Université ParisSud 11, Bâtiment 121, Orsay, France
^{49}
Institut d’Astrophysique de Paris, CNRS UMR7095, Université Pierre & Marie Curie, 98bis boulevard Arago, Paris, France
^{50}
Institut de Ciències de l’Espai, CSIC/IEEC, Facultat de Ciències, Campus UAB, Torre C5 par2, Bellaterra 08193, Spain
^{51}
Institute of Astronomy and Astrophysics, Academia Sinica, Taipei, Taiwan
^{52}
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{53} Institute of Theoretical Astrophysics, University of Oslo, Blindern, Oslo, Norway
^{54}
Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, La Laguna, Tenerife, Spain
^{55}
Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, Santander, Spain
^{56}
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, USA
^{57}
Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13 9PL, UK
^{58}
Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
^{59}
LERMA, CNRS, Observatoire de Paris, 61 avenue de l’Observatoire, Paris, France
^{60}
Laboratoire AIM, IRFU/Service d’Astrophysique – CEA/DSM – CNRS – Université Paris Diderot, Bât. 709, CEASaclay, 91191 GifsurYvette Cedex, France
^{61}
Laboratoire Traitement et Communication de l’Information, CNRS (UMR 5141) and Télécom ParisTech, 46 rue Barrault, 75634 Paris Cedex 13, France
^{62}
Laboratoire de Physique Subatomique et de Cosmologie, CNRS/IN2P3, Université Joseph Fourier Grenoble I, Institut National Polytechnique de Grenoble, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{63}
Laboratoire de l’Accélérateur Linéaire, Université ParisSud 11, CNRS/IN2P3, Orsay, France
^{64}
Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{65}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{66}
MilliLab, VTT Technical Research Centre of Finland, Tietotie 3, Espoo, Finland
^{67}
NRAO, P.O. Box 2, Rt 28/92, Green Bank, WV 249440002, USA
^{68}
National University of Ireland, Department of Experimental Physics, Maynooth, Co. Kildare, Ireland
^{69}
Niels Bohr Institute, Blegdamsvej 17, Copenhagen, Denmark
^{70}
Observational Cosmology, Mail Stop 36717, California Institute of Technology, Pasadena, CA, 91125, USA
^{71}
Optical Science Laboratory, University College London, Gower Street, London, UK
^{72}
SISSA, Astrophysics Sector, via Bonomea 265, 34136, Trieste, Italy
^{73}
SUPA, Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
^{74}
School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff, CF24 3AA, UK
^{75}
Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, Moscow, 117997, Russia
^{76}
Space Sciences Laboratory, University of California, Berkeley, California, USA
^{77}
Stanford University, Dept of Physics, Varian Physics Bldg, 382 via Pueblo Mall, Stanford, California, USA
^{78}
Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex 4, France
^{79}
Universities Space Research Association, Stratospheric Observatory for Infrared Astronomy, MS 2113, Moffett Field, CA 94035, USA
^{80}
University of Granada, Departamento de Física Teórica y del Cosmos, Facultad de Ciencias, Granada, Spain
^{81}
University of Miami, Knight Physics Building, 1320 Campo Sano Dr., Coral Gables, Florida, USA
^{82}
Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
Received: 7 January 2011
Accepted: 17 June 2011
Using Planck maps of six regions of low Galactic dust emission with a total area of about 140 deg^{2}, we determine the angular power spectra of cosmic infrared background (CIB) anisotropies from multipole ℓ = 200 to ℓ = 2000 at 217, 353, 545 and 857 GHz. We use 21cm observations of Hi as a tracer of thermal dust emission to reduce the already low level of Galactic dust emission and use the 143 GHz Planck maps in these fields to clean out cosmic microwave background anisotropies. Both of these cleaning processes are necessary to avoid significant contamination of the CIB signal. We measure correlated CIB structure across frequencies. As expected, the correlation decreases with increasing frequency separation, because the contribution of highredshift galaxies to CIB anisotropies increases with wavelengths. We find no significant difference between the frequency spectrum of the CIB anisotropies and the CIB mean, with ΔI / I = 15% from 217 to 857 GHz. In terms of clustering properties, the Planck data alone rule out the linear scale and redshiftindependent bias model. Nonlinear corrections are significant. Consequently, we develop an alternative model that couples a dusty galaxy, parametric evolution model with a simple halomodel approach. It provides an excellent fit to the measured anisotropy angular power spectra and suggests that a different halo occupation distribution is required at each frequency, which is consistent with our expectation that each frequency is dominated by contributions from different redshifts. In our bestfit model, half of the anisotropy power at ℓ = 2000 comes from redshifts z < 0.8 at 857 GHz and z < 1.5 at 545 GHz, while about 90% come from redshifts z > 2 at 353 and 217 GHz, respectively.
Key words: diffuse radiation / submillimeter: diffuse background / submillimeter: galaxies / cosmology: observations
© ESO, 2011
1. Introduction
In addition to instrument noise, deep cosmological surveys in the farinfrared to millimeter spectral range are limited in depth by confusion from extragalactic sources (e.g. Blain et al. 1998; Lagache et al. 2003; Dole et al. 2004; FernandezConde et al. 2008; Nguyen et al. 2010). This limitation arises from the high density of faint, distant galaxies that produce signal fluctuations within the telescope beam. As a consequence, the cosmic infrared background (CIB), which records much of the radiant energy released by processes of structure formation that have occurred since the decoupling of matter and radiation following the Big Bang (Puget et al. 1996; Hauser & Dwek 2001; Dole et al. 2006), is barely resolved into its constituents. Indeed, less than 10% of the CIB is resolved by Spitzer at 160 μm (Béthermin et al. 2010a), about 10% by Herschel at 350 μm (Oliver et al. 2010) and a negligible fraction is resolved by Planck^{1} (FernandezConde et al. 2008). Thus, in the absence of foreground (Galactic dust) and background (cosmic microwave background, CMB) emissions, and when the instrument noise is subdominant, maps of the diffuse emission at the angular resolution probed by the current surveys reveal a web of structures, characteristic of CIB anisotropies. With the advent of large area farinfrared to millimeter surveys (Herschel, Planck, SPT, and ACT), CIB anisotropies constitute a new tool for structure formation and evolution study.
Cosmic infrared background anisotropies are expected to trace largescale structures and probe the clustering properties of galaxies, which in turn are linked to those of their hosting dark matter halos. Because the clustering of dark matter is well understood, observations of anisotropies in the CIB constrain the relationship between dusty, starforming galaxies and the dark matter distribution. The connection between a population of galaxies and dark matter halos can be described by its halo occupation distribution (HOD; Peacock & Smith 2000; Seljak 2000; Benson et al. 2000; White et al. 2001; Berlind & Weinberg 2002; Cooray & Sheth 2002), which specifies the probability distribution of the number of objects with a given property (e.g., luminosity, stellar mass, or starformation rate) within a dark matter halo of a given mass and their radial distribution within the halo. The HOD and the halo model provide a powerful theoretical framework for describing the connection between galaxies and dark matter halos. Once decisions are made about which properties of the halos and their environment the HOD depends upon, what the moments of the HOD are and what the radial profile of objects within halos is, the halo model can be used to predict any clusteringrelated observable. In particular, the halo model predicts that the bias, describing the clustering of galaxies in relation to the dark matter, becomes scaleindependent at large scales. This assumption of a scaleindependent bias is often made in modelling the CIB.
The way galaxies populate dark matter halos is not the only ingredient that enters into the CIB anisotropy modelling. Correlated anisotropies also depend on the mean emissivity per comoving unit volume of dusty, starforming galaxies, that results from dusty galaxies evolution models. Such models are more and more constrained thanks to the increasing number of observations (mainly galaxies number counts and luminosity functions), but remain largely empirical. So far, CIB anisotropy models have combined (i) a scaleindependent bias clustering with a very simple emissivity model based on the CIB mean (Knox et al. 2001; Hall et al. 2010) or an empirical model of dusty galaxy evolution (Lagache et al. 2007) or the predictions of the physical model by Granato et al. (2004) for the formation and evolution of spheroidal galaxies (Negrello et al. 2007); (ii) a HOD with the Lagache et al. (2003) dusty galaxies evolution model (Amblard & Cooray 2007; Viero et al. 2009); and (iii) a merger model of dark matter halos with a very simple dust evolution model (Righi et al. 2008).
The angular power spectrum of CIB anisotropies has two contributions, a whitenoise component caused by shot noise and an additional component caused by spatial correlations between the sources of the CIB. Correlated CIB anisotropies have been measured at 3330 GHz by AKARI (Matsuura et al. 2011), 3000 GHz by IRAS/IRIS (Pénin et al. 2011b), 1875 GHz by Spitzer (Lagache et al. 2007; Grossan & Smoot 2007), 1200, 857, 600 GHz by BLAST and SPIRE (Viero et al. 2009; Amblard et al. 2011), and 220 GHz by SPT (Hall et al. 2010) and ACT (Dunkley et al. 2011). Depending on the frequency, the angular resolution and size of the survey these measurements can probe two different clustering regimes. On small angular scales (ℓ ≥ 2000), they measure the clustering within a single dark matter halo and accordingly the physics governing how dusty, starforming galaxies form within a halo. On large angular scales, CIB anisotropies measure clustering between galaxies in different dark matter halos. These measurements primarily constrain the largescale, linear bias, b, of dusty galaxies, which is usually assumed to be scaleindependent over the relevant range. Given their limited dynamic range in scale, current measurements are equally consistent with an HOD model, a powerlaw correlation function or a scaleindependent, linear bias. All models return a value for the largescale bias that is 2–4 times higher than that measured for local, dusty, starforming galaxies (where b ≃ 1).
Owing to its frequency coverage from 100 to 857 GHz, the High Frequency Instrument (HFI) onboard Planck is ideally suited to probe the dark matter – starformation connection. Planck (Tauber et al. 2010; Planck Collaboration 2011a) is the thirdgeneration space mission to measure the anisotropy of the CMB. It observes the sky in nine frequency bands covering 30–857 GHz with high sensitivity and angular resolution from 31′ to 5′. The Low Frequency Instrument (LFI; Mandolesi et al. 2010; Bersanelli et al. 2010; Mennella et al. 2011) covers the 30, 44, and 70 GHz bands with amplifiers cooled to 20K. The 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.1K. 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 2011b). 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 results are given in Planck Collaboration (2011c–w).
The primary objective of this paper is to measure with Planck HFI the CIB anisotropies caused by the clustering of starforming galaxies. To achieve this, we analyse small regions of sky with a total area of about 140deg^{2}, where we are able to cleanly separate the foreground (Galactic cirrus) and background (CMB) components from the signal. Unlike previous CIB anisotropy studies (but see Pénin et al. 2011b), we do not remove the cirrus by fitting a powerlaw power spectra at large scales, but use an independent, external tracer of diffuse dust emission (the Hi gas). We accurately measure the instrumental contributions (noise, beam) to the power spectra of CIB anisotropies and use a dedicated optimal method to measure power spectra (Ponthieu et al. 2011). All these steps allow us to recover for the first time the power spectra of CIB anisotropies from 200 ≤ ℓ ≤ 2000 at four frequencies simultaneously: 217, 353, 545 and 857 GHz.
The paper is organized as follows: in Sect. 2 we present the data we are using, the field selection and the removal of foreground and background components (CMB, Galactic cirrus, bright point sources) from the CIB. In Sect. 3 we discuss the different contributions to the power spectra of the residual maps. Section 4 describes how we estimated the power spectrum, its bias, and errors. Our main results are presented in Sect. 5. This section also describes our modelling and discusses the clustering of highredshift, dusty galaxies. We conclude in Sect. 6. In the appendices we show two flow charts summarizing the data processing and cleaning, and the power spectra measurements (Appendix A), and we give some details about the dusty starforming galaxy evolution model (Appendix B) and the halo model (Appendix C) we are using. Throughout the paper we use the WMAP7 cosmological parameters for standard ΛCDM cosmology (Larson et al. 2011).
2. Selected fields and data cleaning
2.1. Planck data
We used Planck channel maps of the HFI at 5 frequencies: 143, 217, 353, 545 and 857 GHz. Their characteristics and how they were created is described in detail in the companion paper on HFI early processing (Planck HFI Core Team 2011b). In summary, the channel maps correspond to temperature observations for the two first sky surveys by Planck. The data are organized as timeordered information (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 (see below). Timeordered informations of raw bolometer data are first processed to produce cleaned timelines and to set flags to mark data we do not currently fit. 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 4K cooler mechanical vibrations; and (6) deconvolution of the bolometer time constant. Focal plane reconstruction and beam shape estimation is made using observations of Mars. The simplest description of the beams, an elliptical Gaussian, leads to fullwidth at halfmaximum (FWHM) values, θ_{S}, given in Table 3 of Planck HFI Core Team (2011b) (i.e., 7.08′,4.71′,4.50′,4.72′ and 4.42′ from 143 to 857 GHz, with an uncertainty between 0.12′ and 0.28′). From the cleaned TOI and the pointing, channel maps were computed using all the bolometers at a given frequency. The path from TOI to maps in the HFI DPC is schematically divided into three steps: ringmaking, ring offset estimation, and mapmaking. The first step combines the data within a stable pointing period, during which the same circle on the sky is scanned repeatedly to create rings with higher signaltonoise ratio, taking full advantage of the redundancy of observations provided by the Planck scanning strategy. The lowfrequency component of the noise is accounted for in a second step by using a destriping technique that models this component as an offset of the ring values. Finally, cleaned maps are produced by coadding the offsetcorrected rings. The maps are produced in Galactic coordinates, using the HEALPix pixelisation scheme (see http://healpix.jpl.nasa.gov and Górski et al. 2005). Photometric calibration is performed either at ring level (using the CMB dipole) for the lower frequency channels or at the map level (using FIRAS data) for the higher frequency channels (545 and 857 GHz). The absolute gain calibration of the HFI Planck maps is known to better than 2% for the lower frequencies (143, 217 and 353 GHz) and 7% for the higher frequencies (545 and 857 GHz), as summarised in Planck HFI Core Team (2011b) Table 3. Intercalibration accuracy between channels is better than absolute calibration.
We made use of the socalled DX4 HFI data release, a dataset from which the CMB has not been removed. We used the 217, 353, 545, and 857 GHz channels for the CIB analysis and the 143 GHz channel for CMB removal. Maps are given either in MJysr^{1} (with the photometric convention νI_{ν} = const.^{2}) or μK_{CMB}, the conversion between the two was exactly computed using the bandpass filters (see Planck HFI Core Team 2011b).
CIB field description: centre (Galactic coordinates), size, mean and dispersion of Hi column density.
Fig. 1 From left to right and top to bottom: N1, AG, SP, LH2 and bootes fields overlaid on IRIS 100μm map (MivilleDeschênes & Lagache 2005). Fields Bootes 1 and 2 are both included in the large rectangle. All IRIS images have the same dynamic range, with a linear colour scale ranging from dark red to white from 0 to 2 MJysr^{1}. 
2.2. Extragalactic fields with high angular resolution HI data
Although Planck is an allsky survey, we restricted our first CIB anisotropy measurements to a few fields at high Galactic latitude to minimize the Galactic dust contamination. The choice of the fields was driven by the availability of Hi data at an angular resolution close to HFI.
The 21cm Hi spectra used here were obtained with the 100m Green Bank Telescope (GBT) over the period 2005 to 2010. Details of this highlatitude survey are presented by Martin et al. (in prep.). The total area mapped is about 825deg^{2}.
The spectra were taken with onthefly mapping. The primary beam of the GBT at 21cm has a FWHM of 9.1′, and the integration time (4s) and telescope scan rate were chosen to sample every 3.5′, more finely than the Nyquist interval, 3.86′. The beam is only slightly broadened to 9.4′ in the inscan direction. Scans were made moving the telescope in one direction (galactic longitude or right ascension), with steps of 3.5′ in the orthogonal coordinate direction before the subsequent reverse scan.
Data were recorded with the GBT spectrometer by inband frequency switching, yielding spectra with a velocity coverage −450 ≤ V_{LSR} ≤ +355 kms^{1} at a resolution of 0.80 kms^{1}. Spectra were calibrated, corrected for stray radiation, and placed on a brightness temperature (T_{b}) scale as described in Blagrave et al. (2010), Boothroyd et al. (in prep.), and Martin et al. (in prep.). A thirdorder polynomial was fitted to the emissionfree regions of the spectra to remove any residual instrumental baseline. The spectra were gridded on the natural GLS (SFL) projection to produce a data cube. Some regions were mapped two or three times. With the broad spectral coverage, all Hi components from local gas to highvelocity clouds are accessible.
We selected from this GBT cirrus survey the six faintest fields in terms of Hi column densities. Their main characteristics are given in Table 1 and the IRAS 100μm maps are shown in Fig. 1. They all have very low dust contamination and consequently Hi column densities, including the faintest allsky sight line (referenced as LH2 in the Table). The field areas are between 16 to 25deg^{2} for a total coverage of about 140deg^{2}. Going to higher average Hi column densities (N(Hi) > 2.5 × 10^{20} cm^{2}) is not recommended because dust emission associated with molecular gas starts to contaminate the signal (see Fig. 9 of Planck Collaboration 2011t) and Hi is no longer a good tracer of dust emission.
The HEALPix HFI maps were reprojected onto the small Hi maps by binning the original HEALPix data (sampled with HEALPix Nside of 2048, corresponding to a pixel size of 1.72′) into Hi map pixels (pixel size 3.5′ for all fields). An average of slightly more than four HEALPix pixels were averaged for each small map pixel.
2.3. Removing the bright sources from HFI maps
We removed from the maps all sources listed in the Planck Early Release Compact Source Catalog (ERCSC) (Planck Collaboration 2011c). This represents only a few sources per field (if any), but the bright source removal is important for both power spectrum analysis and CMB map construction. It is also important to know the flux limit to compute the radio and dusty galaxy shotnoise contribution to the power spectra. Since our fields have roughly the same (very low) dust contamination, source detection is not limited by cirrus. Indeed, the flux cut is set by extragalactic source confusion at high frequencies and CMB contamination at low frequencies. The same flux cut can therefore be applied to all our fields. We took the minimum ERCSC flux densities in our fields as the flux cuts. They are given in Table 3.
In practice, point source removal is performed in the original HFI HEALPix data prior to reprojection. For each source, a disc of size equal to the FWHM of the beam centred on the source position is blanked. Holes caused by missing data are then filled by a gapfilling process, which interpolates/extrapolates into the hole the values of neighbouring pixels.
2.4. Removing the CMB contamination from HFI maps
Cosmic microwave background anisotropies contribute significantly to the total HFI map variance in all channels at frequencies up to and including 353 GHz. The detection and characterisation of CIB anisotropies at these frequencies requires the separation of the contribution from the CMB.
The present work focuses on very clean regions of the sky, for which Galactic foregrounds are very faint and are monitored using ancillary Hi observations. To remove the CMB in the fields retained for our analysis, we used a simple subtraction technique. While this simple method could be improved in future, it enables us to reliably evaluate CMB residuals, noise contamination, and to propagate errors due to imperfect instrumental knowledge. It also guarantees that highfrequency CIB anisotropy signals will not leak into lower frequency, CMBfree maps.
We removed the CMB contamination in the 217 and 353 GHz channels by subtracting a CMB template obtained from the lower frequency data. We modelled the data for each frequency ν as (1)where x_{ℓm}(ν) represents the channel data at frequency ν, the CMB map, the CIB map and n_{ℓm}(ν) is a noise term comprising (if needed) any other astrophysical contaminant. The effect of the beam was accounted for with a (channel dependent) multiplicative factor, b_{ℓ}(ν) (see Sect. 3.2). For the purpose of CMB removal, b_{ℓ}(ν) was obtained from the Gaussian bestfit to the effective HFI beam of the channel maps.
Fig. 2 Wiener filter applied to the 143 GHz map for CMB subtraction. The filter essentially cuts out high multipoles where the CMBtonoise ratio of the 143 GHz map is low. Whereas this filter has to be known for estimating and subtracting the contribution of the residual CMB and 143 GHz instrument noise to the power spectrum of CMBcleaned channels, the exact value of the filter is not really critical. 
Fig. 3 Power spectra of the different components for field SP (the figure is similar for the other fields). Power spectra of the 217, 353, 545, and 857 GHz Planck maps (continuous black line) are compared to the noise power spectra (diamonds), to the CMBcleaned power spectra (red), and to the CMB and interstellar dustcleaned power spectra (green). In this plot signal power spectra have not been corrected for the beam window function. Noise power spectra are computed using halfpointing period maps, as explained in Sect. 3.3. 
At 100 and 143 GHz, we assumed that in the fields of interest only CMB and noise is present. Cosmic infrared background anisotropies are very small, and in the selected fields the contamination by other sources (e.g., cirrus) is negligible (see Sect. 3.4). In principle, both channels can be used to make a template of CMB emission. However, the 100 GHz channel is significantly less sensitive than 143 GHz and has an angular resolution two times worse than the 217 and 353 GHz channels. Therefore, we only used the 143 GHz channel as a CMB template. We corrected the 217 GHz maps for CMB contamination as follows (1)where w_{ℓ} is a Wiener filter, designed to minimize the total contamination of the new map, y(ν_{217}), by CMB and instrument noise. The 353 GHz map is corrected from CMB contamination in a similar way. Note that this cleaning was performed on a large region comprising all small fields used in the present analysis. The Wiener filter is obtained as (3)where is the current bestfit CMB model spectrum, and Y_{ℓ}(ν_{143}) is the power spectrum of the 143 GHz map. The Wiener filter w_{ℓ} is close to 1 at low ℓ, and close to 0 at large ℓ (see Fig. 2). Note that errors on the beam estimate, on the assumed CMB power spectrum, or on the estimation of the 143 GHz power spectrum would result in suboptimal filtering rather than in biases. We checked that the CMB remaining in the CMBcleaned maps does not change significantly with different assumptions leading to different w_{ℓ}.
Errors in photometric calibration between channels are a problem. Although these errors are estimated to be small (2% at 143, 217, and 353 GHz), they may result in residual CMB at low ℓ. They are accounted for in the processing, as detailed in Sect. 4.2.1.
Figure 3 shows the HFI power spectra of the raw and CMBcleaned maps for one of the six fields. The CMB correction is very large at 217 GHz: the residual is a factor ~100 below the raw power spectrum at ℓ ≃ 430 (it is a factor ~2 below at 353 GHz). Note that whereas this illustrates the effectiveness of CMB removal, it is also a source of worry about the impact of relative calibration errors for the 217 GHz channel. However, the power spectrum after CMB cleaning is ~1% of the original map power spectrum only for ℓ ≤ 600. Cosmic microwave backgroundcleaned maps are shown in Fig 4. We see that the CMB has been efficiently removed.
Finally we remark that an alternative method of removing CMB contamination, based on an internal linear combination of frequency maps and a needlet analysis (Delabrouille et al. 2009), was extensively studied and used in some of the Planck early papers, but it was not well suited to our purposes. The method tended to perform well over large patches of sky but left visible, largescale residuals in the sky patches of interest, and had leakage between the faint CIB and the CMB when other components (noise and Galactic cirrus) are present.
2.5. Removing the cirrus contamination from HFI maps
From 100 μm to 1mm, at high Galactic latitude and outside molecular clouds a tight correlation is observed between farinfrared emission from dust and the 21cm emission from gas^{3} (e.g. Boulanger et al. 1996; Lagache et al. 1998). Hi can thus be used as a tracer of cirrus emission in our fields, and indeed it is the best tracer of diffuse interstellar dust emission.
Fig. 4 Maps of the 26 deg^{2} of the N1 field, from left to right: 217, 353, 545 and 857 GHz. From top to bottom: raw HFI maps; CMB and ERCSC sourcecleaned maps; residual maps (CMB, sources, and cirruscleaned); residual maps smoothed at 10′ to highlight the CIB anisotropies. The joint structures clearly visible (bottom row) correspond to the anisotropies of CIB. Residual point sources are also visible. They have fluxes lower than the fluxes of the ERCSC removed sources. They have no impact on our analysis. 
Fig. 5 Hi and dust maps for two fields: SP (top)and AG (bottom). The first two maps on the left show the Hi components (Local and IVC for SP, IVC and HVC for AG), the third maps show the 857 GHz emission associated with Hi () and the maps on the right side show the HFI 857 GHz maps. Those HFI maps have been convolved by the GBT beam to allow a better comparison by eye. Hi maps are given in units of 10^{20} atomscm^{2}. Note the correlation of the dust emission with the different Hi velocity components and its variation from field to field. 
Hi components – The Hi data in each field show different velocity components: a local component, typical of highlatitude Hi emission, intermediatevelocity clouds (IVCs) and highvelocity clouds (HVCs). These clouds are typically defined as concentrations of neutral hydrogen at velocities inconsistent with a simple model of differential Galactic rotation. The distinction between IVCs and HVCs is loosely based on the observed radial velocities of the clouds; IVCs have radial velocities with respect to the local standard of rest (LSR) of 30 ≤ V_{LSR} ≤ 90km s^{1}, while HVCs have velocities V_{LSR} > 90kms^{1}. Highvelocity clouds might be infalling clouds fueling the Galaxy with lowmetallicity gas, whereas IVCs might have a Galactic origin (e.g. Richter et al. 2001). For each field, we constructed integrated Hi emission maps of the three velocity components. The selection of the velocity range for each component was based on inspection of the median 21cm spectrum and of the rms 21cm spectrum (i.e., the standard deviation of each channel map). It is fully described in the Planck Collaboration (2011t). The Hi maps were then converted to Hi column density using the optically thin approximation: (4)where T_{b} is the 21cm brightness temperature and v the velocity. Corrections have been applied for opacity (see Planck Collaboration 2011t), they are lower than 5% for our CIB fields. As illustrated in Fig. 5, the different fields have clearly distinct Hi contributions, with, e.g., no local component in the direction of the AG field.
HIdust correlation – To remove the cirrus contamination from the HFI maps, we need to determine the farIR to mm emission of the different Hi components identified with the 21cm observations. We assumed that HFI maps, I_{ν}(x,y), at frequency ν can be represented by the following model (5)where is the column density of the ith Hi component, is the farIR to mm – Hi correlation coefficient of component i at frequency ν and C_{ν}(x,y) is a residual. The correlation coefficients (often called emissivities) were estimated using a χ^{2} minimization given the Hi and HFI data and the model (Eq. (5)). Although the Hi column densities of the different velocity components are quite similar (see Fig. 5), the emissivities may vary by factors of more than 10 between the local/IVC and HVC (see Planck Collaboration 2011t), so it is important to consider them separately. The emissivities can be used to characterise the opacity and temperature of the dust emission in the different components. This is beyond the scope of this paper, but it is extensively discussed in the Planck Collaboration (2011t).
Cirrus contamination removal – We removed from the HFI maps the Hi velocity maps multiplied by the correlation coefficients. Maps are shown in the last two columns of Fig. 5 for two fields. The removal was made at the HFI angular resolution, even though the Hi map is of lower resolution. This is not a problem because cirrus, with a k^{3} powerlaw power spectrum (MivilleDeschênes et al. 2007), has negligible power between the GBT and HFI angular resolutions, in comparison to the power in the CIB.
Residual maps and power spectra – Figure 3 shows the HFI power spectra before and after the dust removal in the SP field. Cirrus removal has more impact for the two highfrequency channels. At 217 GHz, the correction is very small (13% at ℓ = 500). This method of using Hi data to remove the cirrus contamination from power spectra has also been successfully applied by Pénin et al. (2011b) at higher frequencies than ours, where the cirrus contamination is higher. The authors have been able to isolate precisely the CIB anisotropies power spectra at 1875 and 3000 GHz with Spitzer and IRAS/IRIS, in the N1 field.
The residual maps at the HFI angular resolution are shown in Fig. 4 for the N1 field. We clearly see that the cirrus has been efficiently removed. The bottom row shows the residual maps, smoothed at 10′. Common structures, corresponding to the CIB anisotropies, are clearly visible at the four frequencies. Table 2 gives the Pearson correlation coefficients between the CIB anisotropy maps. They are about 0.9 between the 545 and 857 GHz maps and 0.5 between the 217 and 857 GHz CIB maps. The decrease when the frequency difference between the maps is larger is expected because the contribution of highredshift galaxies to the CIB (and its anisotropies) increases with wavelength. This is illustrated in Fig. 6, extracted from Béthermin et al. (2011), where we show the redshift distribution of the CIB. The redshift distribution of correlated CIB anisotropies is discussed in FernandezConde et al. (2008, 2010) and Pénin et al. (2011a).
Pearson correlation coefficient between CIB anisotropy maps.
3. Astrophysical and instrumental components of residual HFI maps power spectra
Once the CMB and cirrus have been removed, there are three main astrophysical contributors to the power spectrum at the HFI frequencies: two from dusty starforming galaxies (with both shot noise, , and clustering, , components), and one from radio galaxies (with only a shotnoise component, , the clustering of radio sources being negligible, see Hall et al. (2010)). If the instrument noise and the signal are not correlated, the measured power spectrum C_{ℓ}(ν) is (6)where b_{ℓ}(ν) is the beam window function, and N_{ℓ}(ν) the power spectrum of the instrument noise. Note that we neglect here the SunyaevZeldovich (SZ; Sunyaev & Zeldovich 1980) contribution to the power spectra. Extrapolation of SPT (Lueker et al. 2010) and ACT (Dunkley et al. 2011) constraints show that SZ is negligible compared to CIB anisotropies at ν ≥ 217GHz. Our goal is to accurately measure , which we present and extensively discuss in Sect. 5. We begin by discussing all other components of Eq. (6) in this section.
3.1. Shot noise
The shot noise arises from sampling of a background composed of a finite number of sources. We assumed the distribution is Poisson, so that its power spectrum is independent of ℓ. If we identify and remove all sources brighter than S_{cut}, the shot noise from the remaining sources fainter than S_{cut} is given by (e.g., Scott & White 1999) (7)where S is the source flux and dN/dS the differential number counts. These counts can be directly measured or derived from evolution models of the relevant population of galaxies (dusty, starforming and radio galaxies in our case).
Fig. 6 Contribution to the CIB per redshift slice, extracted from Béthermin et al. (2011). The black solid line is the CIB spectrum predicted by the model. The contribution to the CIB from 0 < z < 0.3, 0.3 < z < 1, 1 < z < 2 and z > 2 galaxies is given by the red shortdashed, green dotdashed, blue three dotdashed and purple longdashed lines, respectively. Lower limits coming from the stacking analysis at 100 μm, 160 μm (Berta et al. 2010), 250 μm, 350 μm, 500 μm (Marsden et al. 2009), 850 μm (Greve et al. 2010) and 1.1 mm (Scott et al. 2010) are shown as black arrows. The black diamonds give the Matsuura et al. (2011) absolute measurements with AKARI. The black square the Lagache et al. (2000) absolute measurements with DIRBE/WHAM and the cyan line the Lagache et al. (2000) FIRAS measurement. 
3.1.1. Starforming, dusty galaxy shot noise,
We used the recent model of Béthermin et al. (2011) to compute the IR galaxy shotnoise power. This is an updated version of the Lagache et al. (2004) model that better reproduces new observational constraints (e.g., from Herschel). This new, empirical model uses the same galaxy spectral energy distribution (SED) templates as Lagache et al. (2004), but a fully parametric evolution of the luminosity function. The parameters of the model were determined by fitting the infrared/submm number counts, and some midIR luminosity functions, with a MonteCarlo Markov Chain (MCMC). More details on the model are given in Appendix B. The derived shotnoise power is given in Table 3, with uncertainties computed from the MCMC. The quoted numbers include statistical and photometric calibration uncertainties. This model has less energy output at high redshift (z ≃ 2) and consequently lower shotnoise power at long wavelength than the Lagache et al. (2004) model. The shot noise levels depend on the flux cut, which itself has an uncertainty linked to the flux uncertainty in the ERCSC. If we change the flux cut S_{cut} by 30% in Eq. (7) based on the uncertainty in ERCSC fluxes, the power spectra change by less than 5% at all frequencies (and less than 1% at 217 GHz).
As we will discuss in Sect. 5, the dusty galaxy shotnoise level will be a major factor in the interpretation of CIB anisotropy power spectra. Because we are obtaining this value from a model, not measuring it directly in this paper (see Sect. 5), we briefly discuss here the constraints on the model and the plausible range of values using the 857 GHz channel as an example (the same conclusions are reached for the other Planck channels). Fig. 7 shows a compilation of models from the literature superimposed on the latest number counts observed by BLAST and Herschel and the expected shot noise as a function of S_{cut}. First we see, as stated above, that a small variation in S_{cut} leads to only a small variation in shotnoise power. Second, we see that the highest shotnoise level is around 13500Jy^{2}sr^{1}, a factor ~2.3 above our nominal value, but it comes from a model that overestimates the observed number counts by a large factor (3–4 for 50 ≤ S ≤ 300mJy). Models that agree reasonably well with the number counts have a shotnoise level below 8000Jy^{2}sr^{1}. The Béthermin et al. (2011) model has the lowest shot noise. However, it is currently the model that best reproduces all available constraints, from the midinfrared to the millimeter, including the differential contribution of the S_{24} ≥ 80μJy sources to the CIB as a function of redshift, which is a difficult observation to predict. Eventually, the shot noise derived from this model agrees well also with the Herschel/SPIRE measurements in the LockmanSWIRE field, when none of the point sources is removed (Amblard, priv. comm.), as detailed in Sect. 5.3.
Fig. 7 A number of recent models of dustygalaxy evolution and their associated shot noise for different flux cuts at 857 GHz. Top: comparison of the models with the Herschel and BLAST differential numbers counts. Models are from Lagache et al. (2004), Negrello et al. (2007), Le Borgne et al. (2009), Patanchon et al. (2009), Pearson & Khan (2009), Valiante et al. (2009), Béthermin et al. (2011), Franceschini et al. (2010), Lacey et al. (2010), Marsden et al. (2011), RowanRobinson (2009), Wilman et al. (2010). Data points are from Oliver et al. (2010), Béthermin et al. (2010b), Glenn et al. (2010). Bottom: shotnoise level as a function of the flux cut for the same models (same colour and line coding between the two figures). The vertical and horizontal continuous dark lines show the Planck flux cut and shotnoise level from Table 3, respectively. The Béthermin et al. (2011) model is shown by the continuous dark line. This figure shows that models predicting a very high shot noise (e.g. continuous and dashed lightblue, reddashed, continuous and dashed darkblue lines) are incompatible with the measured number counts. 
Flux cut from the ERCSC for our six fields, and the shotnoise power for dusty and radio galaxies appropriate to those cuts (see text).
3.1.2. Radio galaxy shot noise,
The shotnoise power from radio galaxies is subdominant to that from dusty sources at the frequencies relevant to CIB anisotropy analysis. The radio galaxy shotnoise power can be estimated from the model of de Zotti et al. (2005). At frequencies ≤ 100 GHz, the model agrees with the source counts computed using the extragalactic radio sources from the ERCSC. At 143 and 217 GHz, and for fluxes below 300mJy (i.e., the case listed in Table 3) the de Zotti et al. (2005) model agrees with the source counts of Vieira et al. (2010). At higher fluxes the model needs to be scaled to reproduce the number counts obtained using the ERCSC. The estimated scaling factors are 2.03 and 2.65 at 143 and 217 GHz, respectively (see Planck Collaboration 2011i). At even higher frequencies the number counts by the ERCSC are no longer complete. We therefore use the 217 GHz scaling factor to set upper limits for the shot noise. It is negligible compared to at these frequencies (see Table 3). Changing the flux cut by 30% affects the shot noise by 30%, but because the radio contribution is subdominant at the frequencies relevant for CIB anisotropy analysis, it has little impact on our results.
3.2. The beam window function, b_{ℓ}(ν)
Because the HFI beams are not azimuthally symmetric, the scanning strategy has to be taken into account in modelling the effective beam response. We used two different methods to compute the effective beam: FEBeCoP and FICSBell. With FEBeCoP, we computed one effective beam per field, with FICSBell, one effective beam for the entire sky.
FICSBell – The FICSBell method (Hivon et al, in prep.) generalizes the approach of Hinshaw et al. (2007) and Smith et al. (2007) to polarization and to include other sources of systematics. The different steps of the method used for this study can be summarized as follows:

1.
The scanningrelated information (i.e., statistics of the orientation of each detector within each pixel) is computed first, and only once for a given observation campaign. The hit moments are only computed up to degree 4, for reasons described below.

2.
The (Marsbased) beam map or beam model of each detector, d, is decomposed into its spherical harmonic coefficients (8)where B_{d}(r) is the beam map centred on the North pole, and Y_{ℓs}(r) is a spherical harmonic. Higher s indices describes higher degrees of departure from azimuthal symmetry and, for HFI beams, the coefficients are decreasing functions of s at most ℓ considered. It also appears that for ℓ < 3000, the coefficients with s > 4 account for ≤ 1% of the beam throughput. For this reason, only modes with s ≤ 4 are considered in the present analysis (ArmitageCaplan & Wandelt (2009) reached a similar conclusion in their analysis of PlanckLFI beams).

3.
The coefficients computed above are used to generate sspin weighted maps for a given CMB sky realisation.

4.
The spinweighted maps and hit moments of the same order, s, are combined for all detectors involved, to provide an “observed” map.

5.
The power spectrum of this map can then be computed, and compared to the input CMB power spectrum to estimate the effective beam window function over the whole sky, or over a given region of the sky.
MonteCarlo (MC) simulations in which the sky realisations are changed can be performed by repeating steps 3, 4, and 5. The impact of beam model uncertainties can be studied by including step 2 into the MC simulations.
FEBeCoP – As mentioned in Sect. 2.1, map making reduces timeordered data to pixelised maps. Each pixel of a map represents a convolution of the true sky with the combined effect of scanning beam and scan pattern. FEBeCoP computes this combination of beams and scans – the effective beams – as is, in the pixel space. The FEBeCoP methodology and algorithm has been described in Mitra et al. (2011), and Planck HFI Core Team (2011b). Below, we list for completeness the essential steps made in computation of the beam window functions:

1.
For each pixel i in the map (or CIB field) we computed the FourierLegendre transform, B_{ℓ}, of the pixel space effective beams using the formula (9)where is the direction vector of the centre of the ith pixel on the sky, P_{ℓ} represents Legendre polynomials of order ℓ and the integration is performed over the (small) solid angle ΔΩ_{i}, outside which the beam can be taken as zero. This formula can be readily transformed to a discretised form with a careful correction for the “pixel window function” as (10)where the summation is over pixels that fall inside the beam solid angle ΔΩ_{i}, Ω_{pix} is the area of each (equal area) pixel and is the pixel window function that compensates for the systematic error that is introduced when integration over a pixel is replaced by the value of the integrand at the pixel centre times the area of the pixel.

2.
We then computed b_{ℓ} at uniformly sampled directions in each field to find the average window functions. The samples were chosen as the HEALPix pixel centres at a coarser resolution (N_{side} = 128) to ensure uniform sampling. Thus we obtained the average window functions for each frequency and field.

3.
To validate the average window functions obtained using the above prescription, we performed MonteCarlo simulations separately for each field and each frequency. We simulated 16 realisations of the sky starting from a ∝ ℓ^{2} angular power spectrum, which are convolved in two ways – (1) with FEBeCoPgenerated effective beams in pixel space and; (2) with analytical Gaussian beam in harmonic space for a beam size appropriate for the given frequency channel. The convolved maps were then “masked” using a function that is unity in the given field and smoothly (in ~25% of field radius) goes to zero outside the field. Finally, we computed the ratio of the angular power spectra of these two maps, multiplied the ratio by the theoretical window function for the same beam size and averaged over the realisations. Though these “transfer functions” suffer from ringing effects often seen in Fourier transforms of a narrow function, they wiggle around the average window functions, confirming the validity of the latter.
Figure 8 shows the FICSBell and FEBeCoP effective beams at 545 GHz. Also shown is the Gaussian beam with a FWHM of 4.72′ ± 0.21′. This is the average FWHM of the scanning beam, determined on Mars obtained by unweighted averaging the individual detectors FWHM. Each FWHM is that of the Gaussian beam, which would have the same solid angle as that determined by using a full GaussHermite expansion on destriped data (see Planck HFI Core Team 2011b, for more details). We see a quite good agreement between the FICSBell, allsky and FEBeCoP, smallfield effective window functions, with a 2% difference at ℓ ≃ 2000, the highest ℓ that will be considered for CIB anisotropy analysis (see Sect. 4). We also see from the figure that the error on the input scanning beam is larger than this difference and will dominate the uncertainties at high ℓ and high frequency (Sect. 4.2.2). Below we will use the FEBeCoP window functions because they are exactly computed for each of our fields.
Fig. 8 Effective beam window functions (b_{ℓ}) from FICSBell (black) and FEBeCoP (red) at 545 GHz (see Sect. 3.2 for more details). The six FEBeCoP beam window functions from each field are superimposed (red lines). Also shown for comparison is the Gaussian beam with a FWHM of 4.72′ ± 0.2′ (green lines), which is the equivalent FWHM of the beam determined on Mars. 
3.3. Instrument noise, N_{ℓ}(ν)
We can use three different jackknife difference maps to derive noise power spectra: maps made from the first and second halves of each pointing period (a halfpointing period is of the order of 20 min), maps made using half of the focal plane array, and maps using the two different surveys (surveys I and II). In each case the noise power spectrum, N_{ℓ}, is obtained by measuring the power spectrum of the difference maps. The three methods give similar N_{ℓ}, as is illustrated for one frequency and one field in Fig. 9. We chose, however, to use the halfpointing period maps because (1) the two survey maps are only fully covered for the LH2 and SP fields; and (2) there are only three bolometers at 545 and 857 GHz, making halffocal plane maps less accurate. We also computed the noise power spectrum from the difference between the auto and cross power spectrum of the two half maps. In the range of interest, 1500 ≤ ℓ ≤ 2100, where the contribution from the noise becomes important, they agree at better than 0.5, 1, 3, and 4% at 217, 353, 545 and 857 GHz, respectively. Figure 10 shows the noise power spectra for all fields. They are nearly flat, the deviation from flatness is caused by the effect of deconvolution from the instrumental response at high frequency and residual lowfrequency noise. Removing the ERCSC sources has no impact on the noise determination.
Fig. 9 Three independent noisepowerspectrum measurements in the SP field at 353 GHz: red continuous line, half pointing period; green dashed, surveys I and II; black dotdashed, half focal plane array). 
Fig. 10 Instrument noise power spectra of the six fields obtained using halfpointing period maps. From top to bottom: 217, 353, 545 and 857 GHz (continuous: N1, dotted: AG, dashed: SP, dashdotted: Bootes 1, longdash: Bootes 2, dash3 dotted: LH2). 
Figure 3 shows the noise power spectra compared to the HFI map power spectra for one illustrative field. We see that we have a very high signaltonoise ratio. At 545 and 857 GHz, the signal is dominating even at the highest spatial frequencies. At 217 and 353 GHz, the residual signal (i.e., CMB and cirruscleaned) is comparable to the noise at high ℓ (ℓ ≥ 2000 − 2500 depending on the field).
3.4. Additional corrections
Two additional corrections linked to the CMB cleaning were made for the power spectra. First we removed the extra instrument noise that has been introduced by CMB removal: (11)with ν equal to 217 or 353 GHz. N_{ℓ}(ν_{143}) is the noise power spectrum of the 143 GHz map. It is computed as the noise in the other frequency channels, using the halfpointing period maps, following Sect. 3.3.
Second, owing to the lower angular resolution of the 143 GHz channel compared to the 217 and 353 GHz, we also had to remove the CMB contribution that is left close to the angular resolution of the 217 and 353 GHz channels: (12)with F_{p} the pixel and reprojection transfer function (detailed in Sect. 4.1).
Finally, we had to assess the level of the astrophysical components that were removed from or added to the 217 and 353 GHz channels, using the filtered 143 GHz channel as a CMB template. Cirrus emission is highly correlated between 143, 217 and 353 GHz channels. Consequently, filtered cirrus emission was removed from the 217 and 353 GHz. This has no impact on our CIB anisotropy analysis because this extra cirrus removal only modifies the emissivities, with no consequence on our residual maps (it should be understood for a further interpretation of the Hicorrelated dust emission, which is not the goal of this paper). We expect the shotnoise powers to be quite decorrelated for the (143, 217) and (143, 353) sets of maps because the 143 GHz shot noise is dominated by radio sources, whereas the 217 and 353 GHz shot noise is dominated by dusty galaxies (see Table 3). To have an idea of the maximal effect (i.e., perfect decorrelation between shot noise at 143, and 217, and 353 GHz) we computed the contamination by the 143 GHz shot noise, summing the contribution of the radio and dusty galaxies and following (13)The last term accounts for the filtering and “rebeaming” of the 143 GHz map. The contamination is the highest in the 217 GHz channel. It is a factor 1.2 and 120 smaller than the sum of the predicted radio and dusty galaxies shotnoise powers at 217 GHz at ℓ ≃ 200, and 2000, respectively, but is equivalent at ℓ ≃ 1000. Anyway, it is smaller by factors of 20, 2.9 and 325 than the CIB anisotropies at 217 GHz, at ℓ= 200, 1000, and 2000, respectively. Because this is the maximal contamination and because it is quite low (and completely negligible at high ℓ), we did not apply any correction to the CIB anisotropy power spectra.
We still have to consider the case of CIBcorrelated anisotropies at 143 GHz. They have been marginally constrained at 150 GHz by SPT and ACT at high ℓ. The power is <5.2 × 10^{6}μK^{2} and <9.8 × 10^{6}μK^{2} at ℓ = 3000 in Dunkley et al. (2011) and Hall et al. (2010), respectively. This contribution is also completely negligible compared to the signal at 217 GHz.
In conclusion, we can ignore the CIB and cirrus components that are left in the CMB maps.
4. Angular power spectrum estimation
The angular power spectrum estimator used in this work is POKER (Ponthieu et al. 2011), which is an adaptation to the flat sky of the pseudospectrum technique developed for CMB analysis (see e.g. MASTER, Hivon et al. 2002). In brief, POKER computes the angular power spectrum of the masked data (a.k.a., the pseudopower spectrum) and deconvolves it from the power spectrum of the mask to obtain an unbiased estimate of the binned signal angular power spectrum. We summarize the main features of POKER in the following section and then detail how it was used to produce the final estimate of the CIB anisotropy power spectrum and its associated error bars.
In the following, the power spectrum associated to CIB anisotropies will be denoted C_{ℓ} and its unbiased estimator in the flatsky approximation P(ℓ). As already suggested, this final estimate makes use of the power spectrum of the masked data. This socalled pseudopower spectrum will be denoted . In the flatsky approximation, the standard angular frequencies labelled by their zenithal and azimuthal numbers, usually called ℓ and m respectively, are replaced by an “angular” wavevector ℓ; its norm ℓ is equal to the zenithal number (see e.g., the appendix of White et al. 1999). Finally, we will assume that CIB anisotropies arise from a statistically isotropic process. As is the case for the CMB, the CIB fluctuations are viewed as isotropic and homogeneous stochastic variables on the celestial sphere, leading to (14)with a(ℓ) the Fourier coefficients of CIB anisotropies. This assumption is theoretically reasonable, moreover, we checked that a(ℓ)^{2} computed from our CIB maps does not depend on the direction of ℓ.
4.1. POKER
The POKER implementation of the pseudospectrum approach uses the discrete Fourier transform (DFT). For a map of scalar quantity D_{jk} (j,k denote pixel indices), it is defined as where D_{mn} is the set of discrete Fourier coefficients of D_{jk}. For a given wavevector ℓ, labelled by the m and n indices, its corresponding norm is denoted by with m′ = m (respectively n′ and n) if m ≤ N_{x} / 2 and m′ = N_{x} − m if m > N_{x} / 2. The power spectrum of the map is defined as the squaremodulus of its Fourier coefficients, i.e., P(ℓ_{mn}) = D_{mn}^{2}.
The direct DFT of the masked data relates the true Fourier coefficients to the pseudoFourier coefficients of the signal (17)in which is a convolution kernel that depends only on the mask DFT coefficients. Replacing D_{mn} by in the definition of the power spectrum of a given map leads to the power spectrum of the masked data (a.k.a. the pseudopower spectrum). For a signal T plus noise N map, the ensemble averaged of the pseudospectrum tracing a statistically isotropic process, reads (18)where we have introduced the total transfer function F_{m′n′} accounting for the beam, the “mapmaking” pixelisation effects and reprojection from curved, HEALPix maps to flat, square maps. The beam transfer function is given by the beam power spectrum described in Sect. 3.2. The “mapmaking” pixelisation effects are described by the power spectrum of the pixel window function for fullsky maps provided by the HEALPix package (the initial HEALPix maps are built with N_{side} = 2048 corresponding to a pixel size of 1.7′). As explained in Planck HFI Core Team (2011a), timedomain filtering is included as part of the scanning beam, such that any timedomain filtering effects end up in the estimate of the beam instead of as part of F_{m′n′}. Finally, each curved map with a 1.7′ resolution is reprojected onto its tangent, flat space with a pixel size of 3.5′. This induces first a repixelisation effect because the output map is less resolved than the input one and second, a slight displacement of the pixel centres. The cumulative impact of “image deformation” and “repixelisation” is estimated via MonteCarlo: we first generated a set of fullsky maps and computed the MC average of their pseudospectra. This set of maps was then reprojected onto flat maps for which MC average of their pseudospectra in the flatsky approximation were computed. The ratio of the flatsky pseudospectrum divided by the fullsky pseudospectrum gives a measurement of reprojection effect. Note that those simulations have been performed assuming different shapes for the input angular power spectra. The derived reprojection transfer functions agreed perfectly, which underlines the robustness of the approach.
An unbiased estimate of C_{ℓ} is obtained by first subtracting the noise contribution and then deconvolving the mask and beam effects encoded in the convolution kernel . For the sky coverage of the considered fields, the rapid oscillations of the convolution kernels introduce strong correlations between spatial frequencies and make its inversion numerically intractable. (Pseudo)Power spectra are therefore estimated in frequency bands (labelled b hereafter). The binning operator is (19)where Δ_{b} is the number of wave vectors ℓ_{mn} that fall into the bin b. The reciprocal operator that relates the theoretical value of the onedimensional binned power spectrum P_{b} to its value at ℓ_{mn} is (20)For optimal results, the spectral index β should be chosen to get ℓ^{β}C_{ℓ} as flat as possible. For the CMB, β ≃ 2 is the equivalent of the standard ℓ(ℓ + 1) prefactor. For the CIB anisotropies C_{ℓ} scales roughly as ℓ^{1}, and we therefore adopted a binning with β = 1. Nevertheless, we checked that our results were robust against the choice of β: we simulated a power spectrum scaling as ℓ^{1} but reconstructed it assuming β = 0, 1 and 2 in POKER. For each choice of β, the estimated power spectrum perfectly agreed with the input one (for a more complete discussion, see Ponthieu et al. 2011).
The binned pseudopower spectra is (21)and the CIB power spectrum is related to its binned value, C_{b}, via (22)With these binned quantities, Eq. (18) can be rewritten as (23)with (24)An unbiased estimate of the binned angular power spectrum of the signal is thus given by (25)It is easily checked that ⟨ P_{b} ⟩ = C_{b}.
The complete recovery of the CIB anisotropy angular power spectra is therefore made in three steps:

A.
Given the mask, W, associated to our CIB map and the transfer function, F_{m,n}, compute and invert M_{bb′} as given by Eq. (24). The different fields (with a size at most 5.1° × 5.1°) are systematically embedded in a 10° × 10° square map. We use binary masks, i.e., W = 1 for observed pixels (i.e., those kept in the analysis), and W = 0 for unobserved pixels. The estimated power spectra are binned with a bandwidth of Δ_{b} = 200 and the first bin starting at ℓ = 80^{4}.

B.
Derive the noise bias , given by first the instrument noise described in Sect. 3.3 and second the additional corrections given in Sect. 3.4.

C.
Compute the final estimate of C_{b} from Eq. (25) and C_{b} = ⟨ P_{b} ⟩ .
Uncertainties on P_{b} come from sampling variance, noise variance, astrophysical contaminants, and systematic effects. In the following section, we present how we estimated each of these contributions from MonteCarlo simulations.
4.2. Error bar estimation
4.2.1. Statistical uncertainties
As presented in Sect. 4.1, the uncertainties on P_{b} come from signal sampling variance, noise, and uncertainties on the subtraction of CMB and Galactic dust. The first two are described by stochastic processes with known power spectra, the last two come from uncertainties in the weights applied to templates in the subtraction process at the map level.
We developed the tools necessary to simulate maps given any input angular power spectrum for each field and frequency. The simulation pipeline consists of simulating maps given an input power spectrum (in the case of CIB, noise and CMB residual) and maps of template residuals (conservative Gaussian random fractions of the templates). These maps are then combined and analysed by the power spectrum estimator. Each realisation provides an estimated power spectrum with the same statistical properties as our estimate on the data, and alltogether these simulations provide the uncertainties on our estimate. The covariance matrix of P_{b} is (26)with ⟨ · ⟩ _{MC} standing for MonteCarlo averaging. The error bar on each P_{b} is (27)and the binbin correlation matrix is given by its standard definition (28)
Simulation pipeline – The simulated maps are 10° × 10°. They contain six components, accounting for the different ingredients supposedly present in the actual data maps:

1.
A CIB anisotropy component obtained from a random, Gaussian realisation of the CIB anisotropy power spectrum. As a model for such a spectrum, we used a fit to P_{CIB}, estimated from the data additionally multiplied by the power spectrum of the beam, pixel, and reprojection transfer function;

2.
a residual CMB component derived from a random, Gaussian realisation of the power spectrum given in Eq. (12) using the known Wiener filter, beam differences between 143 GHz and other channels, and the WMAP best fit CMB temperature power spectrum;

3.
the instrument noise as derived in Sect. 3.3. Since the noise is slightly coloured, we simulate it using a fit of its measured power spectrum;

4.
extra instrument noise incurred by CMB removal using Eq. (11) as a model of its power spectrum.
In addition to those four ingredients standing for signal and noise (a CMB residual being viewed as an extrasource of “noise” from the viewpoint of CIB), we added the two foreground templates that are removed, with some uncertainties:

5.
A CMB map with a Gaussian uncertainty distributed with 2% and 3% standard deviation at 217 and 353 GHz, respectively (the CMB is negligible at higher frequencies compared to CIB and dust). The 2% and 3% are justified in Sect. 4.3;

6.
an Hi map with a 5%, 10% and 10% standard deviation for its emissivity (local, IVC and HVC components respectively), consistent with both the intercalibration errors (see Sect. 4.3) and the emissivity errors computed by the Planck Collaboration (2011t) using Monte Carlo simulations.
The analysis pipeline – The analysis pipeline works in four steps:

1.
A first set of 1000 MC simulations of CMB residual and noise only is performed to assess first, the pseudopower spectrum of the instrument noise and CMB residual used to debias the simulated data pseudospectrum and, second, the noise variance given by ;

2.
a second set of 1000 MC simulations, including all the components, is performed. For a given simulated map, CMB and dust templates are removed assuming the estimated emissivities of Sect. 2.5 for dust;

3.
the POKER algorithm is then applied to these “foregroundcleaned” maps to obtain a final estimate of the CIB angular power spectrum. In this step, the bias involved in Eq. (25) contains the pseudopower spectrum of the instrument noise model and of the CMB residual model as described in the simulation pipeline;

4.
the total error bars and binbin correlation matrix on P_{b} are obtained as the RMS of 1000 MonteCarlo realisations of the simulation pipeline, as described in the previous paragraph and using Eqs. (27) and (28).
The statistical uncertainties contain:

A.
sampling variance from CIB anisotropies and residual CMB as modelled in Eq. (12);

B.
noise variance from instrument noise and extra noise given by Eq. (11);

C.
uncertainties on the CMB and Hi template subtraction.
In this set of simulations and analysis, we assumed the beam profile as described in Sect. 3.2 and ignored potential beam uncertainties (see the next section for a discussion of this systematic effect). Below, we present the results obtained using the FEBeCoPderived beam profiles, but the estimated power spectra using either the FICSBellderived or the FEBeCoPderived beam agree very well (because the difference between the two beams is small, as shown in Fig. 8). Figure 12 shows the results for all fields and frequencies. We also display in Fig. 13 the binbin correlation matrix, showing that two bins are not correlated by more than 10%.
Fig. 11 Contribution of residuals to the final CIB anisotropy estimate (illustrated here with the SP field at 353 GHz). The bias induced by each dust and CMB component is negligible compared to both the CIB anisotropy signal (black dots) and the statistical noise (in green, including cosmic variance on the noise estimate itself). 
Fig. 12 CIB anisotropy power spectra of the six fields and their combined spectrum. 
Fig. 13 Modulus of the bintobin correlation matrix derived from the simulation pipeline for the SP field at 353 GHz. Spatial frequency bins are not correlated by more than ≃ 10%. 
4.2.2. Systematic errors
Our estimate of each power spectrum is affected by different systematic errors that must be accounted for separately from the statistical errors derived in the previous section. These systematic uncertainties may introduce a bias in the final estimate and/or bias our MonteCarlo estimate of the statistical uncertainties presented above. We review here the different sources of these systematics and evaluate their level.
Mask impact – Our power spectrum estimation is performed on a limited sky patch. This induces power aliasing from angular scales larger than the size of the patch, an effect that increases as the signal power spectrum steepens. POKER is designed to account for this effect (as well as the extra aliasing induced by holes in the map, if present) because the convolution kernel, M_{bb′}, contains the information on mode coupling to the larger scales. We ran POKER on data maps that were embedded in larger regions that are zeropadded. There is no general prescription as to the size of the zeropadding that one should use, but fortunately the results are very insensitive to the particular choice and whether or not the mask is apodized. By comparing different choices we found uncertainties at the 2% level, well below the statistical uncertainties.
Template subtraction impact – Imperfect template subtraction will also lead to “foreground” residuals that slightly bias our final estimate, , of the CIB anisotropy angular power spectrum. This residual level is given at first order by (29)with the pseudospectrum of the template and δα the error on the global amplitude of this template. Figure 11 illustrates the level of these residuals for the particular case of the SP field. Although negligible compared to the statistical errors, they are accounted for in the error budget.
Beam uncertainties – Uncertainty in the beam will also bias the estimate of the power spectrum because:

1.
The beam window function enters the computation of the convolution kernel M_{bb′}. Any beam error biases our estimate of M_{bb′} and thus our final result.

2.
Beam uncertainties will translate into slight misestimation of and , potentially biasing our final estimate.
Moreover, any beam uncertainty will also affect our estimation of the statistical error bars. For example, the contribution of noise to the error bars scales roughly as (see noise curve in Fig. 11). As a consequence, any beam misestimation couples to the noise variance and leads to additional uncertainties on the final power spectrum estimate. This additional error is given at first order by 29where σ_{Nb} is the noise error and δb_{b}(ν) is the beam uncertainty. This may be important at small angular scales, depending on δb_{b}(ν). A similar argument applies equally to sample variance.
In this work we used the current, best determination of the beams. As explained in Sect. 3.2, the uncertainties on the scanning beam dominate these uncertainties. They are ~ 3% to ~ 6% of the FWHM, depending on the frequency (Planck HFI Core Team 2011b).
We simulated the impact of such an error using the {simulation+analysis} pipeline presented in Sect. 4.2.1, assuming the appropriate frequencydependent discrepancy between the beam used to simulate the maps and the beam used to analyse the maps. The bias induced by these uncertainties and their impact on the statistical error bars are derived by comparing the estimated power spectra and their MCvariance with and without the beam discrepancy. The bias so induced is the dominant uncertainty, larger even than the statistical error bars at small angular scales for measurements at 545 and 857 GHz. We took this bias into account in the modelling (Sect. 5.5).
4.3. Robustness
Many additional tests have been done to test the robustness of CIB anisotropy power spectra. First, instead of removing the CMB contribution and fitting only the Hi correlation coefficients, we searched for the best fit simultaneously using the Hi and CMB templates (31)This allowed us to take into account the photometric intercalibration uncertainties, which are about 2% for the CMB channels (Planck HFI Core Team 2011b). We also fitted the lowfrequency channels for only CMB (32)We found that β_{1} and β_{2} agree at the ~0.05% level. The difference between the two coefficients and unity is less than 1% (3%) at 217 GHz (353 GHz). The determination of the β_{i} is however very noisy because of the small area of our fields. They are compatible with β_{1,2} = 1 at 217 GHz, though at 353 GHz they fall below unity by 2–3%. We did not correct for these intercalibration coefficients and took conservative errors on the residual CMB contamination in our error pipeline (2% at 217 GHz and 3% at 353 GHz, see Sect. 4.2). Fitting for both CMB and Hi or just Hi on CMBcleaned maps changes the dustHi emissivities by less than 2%. This was again taken into account in the error simulation pipeline.
To test the reliability of the noise power spectrum measurement, we recomputed the CIB anisotropy power spectrum using the crosscorrelation between halfpointing period maps instead of the autocorrelation. We removed from each half map the same CMB and Hi (with emissivities taken as those of the total map). On average, for the range of ℓ considered in this paper, the two methods agree at better than the 1% level (1% and 0.05% at 217 and 857 GHz, respectively).
Finally, one way to asses the robustness of our determination is to compare the CIB anisotropy power spectra for our different fields. The fields all have different noise, cirrus, and CMB contributions and they are all independent. The comparison made in Fig. 12 shows that they are all compatible within error bars.
4.4. Power spectra from combined fields
Our final estimate of the CIB anisotropy angular power spectra at different frequencies was derived by combining each power spectrum estimated on the six different fields (33)with A an index running over fields and an appropriate weight, to be defined below. The same binning was adopted for each field. The binbin covariance of the “combinedfield” estimator, , is a function of the binbin covariance of each “singlefield” estimator, , and the covariance between two “singlefield” estimators, and follows (34)The error bars on are simply given by and optimal weights are derived by searching for those minimizing these.
Our MC simulations showed that different bins are not correlated by more than 10%, and our fields were widely enough separated that individual measurements are uncorrelated. Neglecting fieldtofield and bintobin covariances, the optimal weights become the usual inverse variance weighting: (35)where is the statistical error bars derived from the MC simulations as previously described.
The final CIB anisotropy power spectra estimates were therefore computed by inserting the inverse variance weights in Eq. (33). The total statistical uncertainties were obtained by inserting those weights in Eq. (34), where stands for the statistical uncertainties only (i.e. the standard quadratic summation). The total systematic errors were obtained by linearly summing the weighted systematic error on each field (any covariance between fields is neglected). We stress that irrespectively of the way weights are derived the full binbin covariance matrix could be computed^{5}. However, the forthcoming cosmological interpretation of the derived power spectra assumes a zero bintobin correlation, and we only provide here the diagonal elements, i.e., error bars, of the final covariance.
Our results are displayed in Fig. 14 and the data points are given in Table 4. Though our weighting is slightly suboptimal, the final angular power spectra are measured with high signaltonoise compared to the singlefield measurements.
Fig. 14 Fieldcombined CIB anisotropy power spectra at 217, 353, 545, and 857 GHz. The dashed line shows the expected sum of the dusty and radio galaxy shotnoise power (from Table 3). The power spectra at 217, 353, and 545 GHz were arbitrary scaled to allow for a better comparison between frequencies (they were multiplied by 2 × 10^{6}, 10^{5} and 10^{3}, respectively). 
CIB anisotropy C_{ℓ}, at 217, 353, 545, and 857 GHz in sr.
5. Overview and comparison with previous work
The measured power and the shot noise predicted by the model discussed previously is shown in Fig. 14. Our measurements do not allow us to detect the shotnoise component, which will dominate on smaller scales than those probed here. However, the predictions are quite close to the measured power at the highest ℓ, indicating that further analysis of the CIB in Planck up to ℓ ~ 3000 might allow a measurement of the shotnoise component. The figure also reveals that the shape of the power spectrum is remarkably similar at the four frequencies, being identical within the 1σ statistical errors for all data points but the last two at 217 GHz (i.e., ℓ ≃ 1717 and 2060, which are 1.6σ from the points at the other frequencies). This suggests that over the range of frequency and ℓ probed here the clustering properties do not evolve much and/or the galaxy populations responsible for CIB anisotropies are the same. We will return to this point in Sect. 5.5. We start in this section by analysing the frequency dependence of the CIB anisotropies and CIB mean, and comparing our anisotropy measurements with previous measurements at the same (or nearby) frequencies.
5.1. Comparing the CIB mean and anisotropy SEDs
The rms fluctuation in the CIB is related to the anisotropy power spectrum as (36)Table 5 gives approximations to this integral using the measured values of ℓC_{ℓ} over the range 200 < ℓ < 2000. Statistical error bars on σ are computed with Monte Carlo simulations using the statistical errors on the power spectra. The table also gives systematic errors (the second error term) corresponding to the photometric calibration uncertainties. Those values can be compared to the CIB absolute level. Cosmic infrared background determinations based on FIRAS data from two groups can be used:

Fixsen et al. (1998) used threedifferent methods to obtain the CIB mean. They average the threespectra to obtain one CIB mean spectrum, and then fit it by amodified black body. They find a dust temperatureT = 18.5 ± 1.2 K, an optical depth τ = (1.3 ± 0.4) × 10^{5} and an emissivity index β = 0.64 ± 0.12. The FIRAS spectrum is quite noisy so that the uncertainties on the parameters are large and the three parameters are highly degenerate.

Lagache et al. (1999) made two determinations of the CIB mean spectrum, using different methods to remove the foregrounds (Hi and Hii) than Fixsen et al. (1998). One CIB mean spectrum is obtained in the Lockman Hole and one on 51% of the sky (to test isotropy). The two spectra agree very well (see Fig. 6 of Lagache et al. 1999). There is a refinement of the measurement in Lagache et al. (2000), which agrees within errors with the previous measurement. Gispert et al. (2000) fit the Lockman Hole spectrum with a modified black body to derive the CIB mean values and uncertainties at certain wavelengths (their Fig. 5). The best fit has T = 13.6 ± 1.5 K, τ = (8.9 ± 2.9) × 10^{5} and β = 1.4 ± 0.2. These parameters are very different to those found by Fixsen et al. (1998), but because of their degeneracy, the spectrum is quite close to that of Fixsen et al. (1998) in the relevant frequency range (200–1500 GHz).
We integrated the two CIB mean fits through the HFI bandpass filters to obtain the values for the absolute signal given in Table 5. The two determinations agree to better than 10%, 1%, and 20% at 857, 545, and 353 GHz, respectively. At 217 GHz they differ by 60%, but are compatible within errors. It is unknown which of the determinations is the best, because the errors on the spectrum are dominated by systematic effects linked to foreground removal that are difficult to estimate. For the uncertainties listed in Table 5 we fixed T and β to their bestfit values and considered only the uncertainty on τ (since the errors on the three parameters are highly correlated). Comparison between the CIB mean and anisotropy SED is shown in Fig. 15. We see that the CIB anisotropy SED is well described by the CIB mean spectrum of Gispert et al. (2000) with the amplitude scaled by 0.15. The CIB mean spectrum of Fixsen et al. (1998) is flatter, with an 857/353 colour of 4.1 compared to 5.3 and 5.4 and an 857/217 colour of 12, compared to 27 and 21, for CIB anisotropies and the Gispert et al. (2000) CIB mean, respectively. A steeper rise in the CIB anisotropy SED compared to the Fixsen et al. (1998) CIB mean has also been seen in Hall et al. (2010). These authors combine SPT and CIB anisotropy data from the literature and obtain an 857/217 colour of ≃ 25 at ℓ = 3000. This compares very well to the CIB anisotropy colour obtained with Planck, integrated over 200 < ℓ < 2000.
In conclusion, we do not see any evidence for a different CIB mean and anisotropy SED, which is consistent with the galaxies that dominate the CIB mean being those responsible for CIB anisotropies. The Gispert et al. (2000) fit of the Lagache et al. (1999) CIB mean spectrum (200 ≤ ν ≤ 1500 GHz) describes both the CIB mean and the CIB anisotropy SED equally well.
Fig. 15 Comparison of the observed CIB mean and anisotropy SED. The CIB measurements are from Lagache et al. (1999) (FIRAS spectrum in black) and Pénin et al. (2011b) (Spitzer and IRIS, pink diamond data points). The green and blue continuous (dashed) lines are the CIB fits from Gispert et al. (2000) and Fixsen et al. (1998) (multiplied by 0.15). The rms fluctuations of the CIB anisotropies, measured for 200 < ℓ < 2000, are shown with the red dots. Their error bars include both statistical and photometric calibration systematic errors (linearly added), as given in Table 5. This figure shows that the CIB anisotropy SED is steeper than the Fixsen et al. (1998) best fit but very close to the Gispert et al. (2000) best fit. We see no evidence for different CIB mean and anisotropy SED. 
Fig. 16 Comparison of SPT (Hall et al. 2010, dark open diamonds) and HFI measurements (red dots) at 217 GHz. The green dashed line corresponds to the SPT shot noise and the green dotdashed line to the clustering model of Hall et al. (2010), the sum of the two is the continuous green line. The clustering model overpredicts by a factor ≃ 2.4 the HFI power at ℓ ~ 800. The blue dashdotted line shows the clustering model divided by this factor. The clustering+shot noise (blue continuous line) now underpredicts the SPT data points, which may be the signature of nonlinear contributions. 
5.2. Comparison with SPT and ACT measurements
SPT and ACT both measured the amplitude of CIB anisotropies, though at higher multipoles than those presented in this paper. The ACT measurement is on scales too small to be directly compared with our measurements (the amplitude is given at ℓ = 3000 while our last data point is at ℓ = 2060). The SPT team computed the residual bandpowers after subtracting the tSZ, kSZ, CMB, and cirrus model components, quoting data points from ℓ = 2000 to 10000. Figure 16 shows the comparison of the HFI measurement at 217 GHz and the SPT measurements at 220 GHz. Because the bandpass filters are not exactly the same, we applied a multiplicative correction factor (colour correction) to overplot the SPT CIB anisotropy power spectrum on that of the HFI. This factor, the square of the HFI/SPT colour, is computed using the CIB SED of Gispert et al. (2000) convolved with the bandpass filters and is equal to 1.04.
To interpret their data, Hall et al. (2010) used a phenomenological model of CIB sources that assumed each galaxy has the same, nonevolving, modified blackbody SED and that their light was a biased tracer of the mass fluctuations, calculated in linear perturbation theory. Moreover, the redshift distribution of the luminosity density was set by two parameters that fix the width and peak redshift.
The green curve of Fig. 16 shows the Hall et al. (2010) model, normalized to the SPT bandpowers. We see that with this normalization, the power at large angular scales is larger by more than a factor of 2 than the HFI data. We also show as the blue curve the same model, except with amplitude adjusted to better agree with the HFI data. This downward adjustment of amplitude could arise from either a reduction in bias or in the amplitude of the mean CIB. This correction, of course, shifts the discrepancy to the smallerscale SPT data. Since we expect the linear theory assumption will be better at large scales than at small, the discrepancy between model and data at small scales may be signaling the importance of nonlinear corrections.
Fig. 17 Comparison of BLAST and HFI measurements at 545 and 857 GHz. HFI data points are the red circles; BLAST data points are the black diamonds. They were colourcorrected for the comparison (the colour was computed using the CIB SED of Gispert et al. (2000), integrated through the BLAST and HFI bandpass filters). The dashed line is the BLAST shot noise (also colourcorrected). Also shown is the BLAST bestfit clustering model (black dashdotted line) and the total contribution (shot noise plus clustering; continuous green line). It provides a good fit to the Planck data. Finally, we report in this figure a revised version of the SPIRE data points from Amblard et al. (2011) (blue triangles from Fig. 18, see text for more details). 
5.3. Comparison with BLAST and SPIRE measurements
Viero et al. (2009) presented BLAST power spectrum measurements at 1200, 857, and 600 GHz in the GOODSSouth field. They detect CIB anisotropy and shotnoise power in the range 940 ≤ ℓ ≤ 10800. The measured correlations are well fitted by a powerlaw over scales of 5–25′, with ΔI / I = 15.1% ± 1.7%. This level with respect to the CIB is the same as that found at the four HFI frequencies (see Sect. 5.1 and Fig. 15). Fitting to a linear theory power spectrum, they find that the BLAST galaxies responsible for the CIB fluctuations have bias parameters, b = 3.9 ± 0.6 and b = 4.4 ± 0.7 at 857 and 600 GHz, respectively. They further interpret their results using the halo model and find that the simplest prescription does not fit very well. One way to improve the fit is to increase the radius at which dark matter halos are truncated in the model (the virial radius) and thereby distribute satellite galaxies over a larger volume. They interpret this as being equivalent to having some starforming galaxies at z ≥ 1 located in the outskirts of groups and clusters.
We show in Fig. 17 the comparison between the BLAST and HFI measurements at 857 and 545 GHz. Because the bandpass filters are quite different (particularly the 600 and 545 GHz BLAST and HFI channels), we applied a colour correction as explained in Sect. 5.2, multiplying the BLAST CIB anisotropy power spectra by 0.7 and 1.05 at 545 and 857 GHz, respectively. We see from Fig. 17 that the BLAST power spectra agree quite well with those from HFI, except that their largestscale data points are systematically higher. This may be caused by contamination by residual Galactic cirrus emission in the BLAST power spectra. Also shown in the figure are shotnoise powers measured by BLAST. Once the colour corrections are applied, they are 1843 and 7326Jy^{2}sr^{1} at 545 and 857 GHz, respectively. Their flux cuts are comparable to ours (they removed two sources above 400mJy at 857 GHz, and no sources at 600 GHz). The measured shot noise levels are 1.6 and 1.2 times higher than the model predictions shown in Table 3 at 545 and 857 GHz, respectively. We also plot their bestfit halo model, which has a minimum halo mass required to host a galaxy of , and an effective bias b_{eff} ≃ 2.4. We see from Fig. 17 that their model is a very good fit to the HFI data points. Indeed, it provides a much better fit of the HFI data points than the BLAST data points!
Fig. 18 Comparison of SPIRE and HFI measurements at 545 and 857 GHz in the overlapping multipole range. HFI data points are the red circles; SPIRE data points from Amblard et al. (2011) are the black triangles. For these data points sources down to 50 mJy have been masked, there is consequently less power compared to HFI. The green triangles (Amblard, priv. comm.) show the SPIRE CIB measurements identical to Amblard et al. (2011), but without a flux cut applied and thus they are directly comparable to the HFI measurement. They should agree with HFI, but are a factor ~1.7 and ~1.2 below the HFI CIB data points for 400 < ℓ < 1500 at 857 and 545 GHZ, respectively. Indeed, they suffer from an overestimate of the cirrus contamination (by a factor 2). Moreover, preliminary crosscalibration between SPIRE and HFI is increasing the Amblard et al. (2011) SPIRE power spectra by 10 and 20% at 857 and 545 GHz, respectively (see Sect. 5.3 for more details). When corrected from these too factors (cirrus and crosscalibration), the SPIRE (blue triangles) and HFI measurements agree well. For this figure, all SPIRE data points were colourcorrected (colours were computed using the CIB SED of Gispert et al. (2000), integrated through the SPIRE and HFI bandpass filters). Error bars include only statistical errors (for SPIRE, error bars are only shown for the green triangles for sake of clarity). 
We now compare our results with the Herschel/SPIRE measurements (Amblard et al. 2011). This comparison has to be made with caution, because the sources are masked down to a flux cut of 50 mJy in SPIRE, which is much smaller than the 540 and 710 mJy flux cut used in HFI at 545 and 857 GHz, respectively. This large difference in the source removal will affect both the shot noise and the correlated components. We thus compare in Fig. 18 our HFI data points with a SPIRE measurement identical to the Amblard et al. (2011), but with no flux cut applied (Amblard, priv. comm.). In this figure, we only show the SPIRE measurements over the multipole range of 200 to 3000 overlapping with Planck. With higher angular resolution maps, SPIRE CIB anisotropy measurements extend down to subarcminute angular scales or ℓ of ~ 2 × 10^{4}. For clarity we do not plot the smallscale power spectrum, but only concentrate on the consistency between HFI and SPIRE CIB anisotropy at larger angular scales. We see from Fig. 18 that SPIRE measurements (green triangles) are significantly below the HFI CIB spectra (red dots). Two elements may explain this difference:

Galactic cirrus: the cirrus signal in the SPIRE field is taken fromexisting measurements in the same field with IRAS100 μm and MIPS 160 μm, and the spectrum is extrapolated from 100 μm to SPIRE wavelengths using the spectral dependence of a Galactic dust model; this procedure is less accurate than the use of Hi, and overestimates the cirrus contamination because the IRAS data contain the CIB anisotropies (Pénin et al. 2011b). Indeed, we checked with our LH2 field, that overlaps with the SPIRE SWIRELockman field to 43%, that in this region the cirrus contamination is negligible at the scales of interests for the comparison (200 < ℓ < 2000). We accordingly added back the cirrus power spectra that were removed from the Amblard et al. (2011) power spectra.

HFI/SPIRE crosscalibration: SPIRE data are calibrated for point sources, with an accuracy of 15% (Swinyard et al. 2010). The pointsource to diffuseemission calibration conversion invokes an effective beam surface that is not perfectly determined yet. Planck/HFI is directly calibrated on diffuse emission. As detailed in the Planck HFI Core Team (2011b), the accuracy is 7% at high frequencies. For now, Planck/HFI therefore has a more accurate photometric calibration for diffuse emission than SPIRE. At this stage, it appears that assuming SPIRE beam surfaces corresponding to the integral of a Gaussian beam limited by diffraction gives a more accurate SPIRE/HFI crosscalibration than the official beam surfaces given in Swinyard et al. (2010). This preliminary crosscalibration has been established by comparing the diffuse emission from several Hershel surveys (some HATLAS, SAG4, and HiGal fields) to Planck/HFI data. Compared to the beam surfaces taken in Amblard et al. (2011)^{6}, they will increase the power spectra by 10% and 20% at 857 and 545 GHz, respectively.
When corrected for the cirrus overestimate and the HFI/SPIRE crosscalibration on diffuse emission, the SPIRE and HFI data points are now compatible (blue triangles and red dots in Fig. 18, respectively). The cirrus correction is the dominant effect up to ℓ ~ 500 and ℓ ~ 1500 at 545 and 857 GHz, respectively.
5.4. A selfconsistent, cosmological, IR, galaxy evolution model
Our interpretation of the CIB anisotropy measurements with HFI relies on a model introduced in Pénin et al. (2011a). The model builds upon the halo model formalism (see Cooray & Sheth 2002, for a review) and populates dark matter halos with galaxies using a HOD, modelling the emission of dusty galaxies using the infrared evolution model of Béthermin et al. (2011, see Appendix C. Our main motivation for developing and using this parametric model is that it allows us to handle in a selfconsistent manner the observational constraints coming from galaxy clustering and the CMB with more galaxyevolutioncentered measurements such as number counts or luminosity functions at various wavelengths and redshifts. This is a key feature of our model.
Previous approaches, such as Amblard & Cooray (2007) and Viero et al. (2009), have used the Lagache et al. (2004) infraredgalaxy evolution model. Compared to Lagache et al. (2004) and Marsden et al. (2011), the parametric evolution of Béthermin et al. (2011) better reproduces the midIR to millimeter statistical observations of infrared galaxies (number counts, luminosity functions, CIB, redshift distributions). This is important because we derive from this model the mean emissivity per comoving unit volume, introduced below, which is a key quantity for interpreting CIB anisotropies.
On the scales of interest to us we can use the Limber approximation (Limber 1954) and write the angular (cross) power spectrum of infrared emission at two frequencies, ν and ν′, and at a multipole ℓ as (e.g., Knox et al. 2001) (37)where χ is the comoving angular diameter distance to redshift z, a = (1 + z)^{1} is the scale factor and is the mean emissivity per comoving unit volume at frequency ν and redshift z. The mean emissivity is derived using the empirical, parametric model of Béthermin et al. (2011)^{7}: (38)The remaining ingredient in the model is thus P_{gg}(k,z). As a foil to the HOD model for P_{gg} we begin with the simple, constant bias model in which (39)where b_{lin} is a redshift and scaleindependent bias and P_{lin}(k) is the linear theory, darkmatter power spectrum. We compute P_{lin}(k) using the fit of Eisenstein & Hu (1998). We will see that this model is not sufficient to explain the CIB anisotropies that we measure. This is not unexpected; at the mean distance of the sources we are probing Mpc scales where nonlinearities and scaledependent bias are important.
By contrast the HOD model computes P_{gg}(k,z) as the sum of the contributions of galaxies within a single dark matter halo (1 h) and galaxies belonging to two different halos (2 h): (40)The details of our assumptions for the 1 h and 2 h terms are given in Appendix C. On large scales P_{2 h} reduces to a constant bias (squared) times the linear theory power spectrum, while the 1halo term becomes a scaleindependent, shotnoise term.
Before comparing our model to Planck observations, let us identify the parameters we hope to constrain with these data. The infrared galaxy evolution model of Béthermin et al. (2011) satisfyingly reproduces current number count observations and luminosity function measurements at the price of introducing a luminosity function characterised by thirteen parameters (see Appendix B). These thirteen parameters fully define the mean emissivities, , given in Eq. (38). The standard cosmological parameters (baryon density, tilt, etc.) mostly define the shape of the linear power spectrum in Eq. (39) and the geometric functions like χ(z). The HOD formalism we introduce in Appendix C requires four more parameters. Pénin et al. (2011a) investigated this full parameter space and its degeneracies and concluded, not surprisingly, that the current generation of infrared galaxy clustering measurements will not allow us to constrain all these parameters simultaneously. Furthermore, they show that most of the constraints on the luminosity function evolution come from number counts and monochromatic luminosity function measurements. In the next section we threfore fix the luminosity function parameters to their bestfit values (from Béthermin et al. 2011) and vary only some of the HOD parameters.
5.5. Confronting the model with observations
To confront the measurements with our model, we use a LevenbergMarquardt algorithm to perform a χ^{2} minimization. The χ^{2} for our data when compared to our model is given by (41)Unless specified otherwise, we use errors including statistical and systematic photometric calibration errors (2, 2, 7 and 7% at 217, 353, 545 and 857 GHz, respectively, as defined in Sect. 2.1), added linearly, and we assume diagonal, uncorrelated error bars as justified above. To reproduce the binning performed while measuring the power spectra, is computed taking the average of at the minimum, maximum and mean ℓ of each bin. We assume a Gaussian likelihood and assume that setting the fixed parameters (e.g., the luminosity function parameters) at their bestfit value is equivalent to marginalising over them. It is important to note that the model power spectrum, includes a shotnoise term (SN) defined in Sect. 3.1. Depending on the precise configuration we study, this shotnoise level is either fitted as an extra parameter or fixed to the predicted value. We also note that the models are derived for the Planck bandpass filters and are colourcorrected to be in Jy(νI_{ν} = const.) as the data points, using the Gispert et al. (2000) CIB SED (colour corrections are (1.08)^{2}, (1.08)^{2}, (1.06)^{2} and (1.00)^{2} at 217, 353, 545 and 857 GHz, respectively).
We remark here that this approach might be limited in two ways. First, we cannot exclude the possibility that our solution is a local extremum and not the global minimum of Eq. (41). While for all the results quoted above we checked the convergence with varying starting points, this particular point would have to be validated using for example a simulated annealing technique. We also did not explore the validity of the Gaussian likelihood approximation. Some of these limitations could be resolved by implementing, for example, a Monte Carlo solver, and we defer this approach to future work.
Fig. 19 In this plot we illustrate the constant bias model at 545 GHz. The orange (solid, dashed and dotdashed) lines correspond to the bestfit linear model, its clustering component, and the shotnoise level, respectively. While this fit was performed using our fiducial emissivity defined in Eq. (38), for illustrative purpose we plot in green the analogous fit using the LDP emissivity. In both cases the required shotnoise level (dotdashed lines) is well above the 68% C.L. predicted by Béthermin et al. (2011) and given in Table 5 (yellow contour). Conversely, the solid yellow line represents the bestfit curve when the shotnoise level is fixed to the expected value (Table 3) and only the constant bias is varied. The fit is obviously unsatisfactory. These results lead us to consider the linear bias model as unphysical, despite the good fit it provides (χ^{2}/d.o.f. ≃ 0.36 (2.52/7)). For illustration purpose, we also show our bestfit powerlaw model defined in Eq. (42) (blue solid line). 
5.5.1. Linear bias model and powerlaw constraints
In this section we will illustrate the discussion looking at the 545 GHz data only, but the same conclusions are reached with the other three frequencies. The relevant results are illustrated in Fig. 19. The data points correspond to our measurements with statistical error bars only.
Fitting simultaneously for a shotnoise level and a constant bias, b_{lin}, defined in Eq. (39), we obtain a good fit (χ^{2} / d.o.f. ≃ 0.36) as visible from the solid orange curve. The best fit bias is b_{lin} = 2.45 ± 0.18 which is consistent with previous results in the literature (Lagache et al. 2007; Viero et al. 2009). The clustering contribution is plotted as the dashed orange line. However, the shotnoise level required to obtain this good fit (dotdashed orange curve) is unrealistically high: (5.6 ± 0.7) × 10^{3} Jy^{2}/sr^{8}. When compared to the expected level from our model whose 68% C.L. amplitude is displayed as the yellow shaded area, our required level is ≃ 5 times higher in power. Given that this model reproduces well all known number counts, the monochromatic luminosity function, CIB measurements (see e.g., the last line in Table 5), and Herschel/SPIRE shotnoise measurements (when no sources are removed, see Sect. 5.3), this excess level is excluded (see also the discussion in Sect. 3.1). We thus conclude that the linear bias model is not a physically realistic fit to our data. The reason for this will be made clear in the next section, as it will appear clearly that the shotnoise level we fit here absorbs a strong nonlinear component.
To further illustrate this point, we show the solid yellow curve that represents the bestfit model when we fix the shot noise to the level predicted by our model and only vary b_{lin}. The fit is obviously unsatisfactory. We also note that the specific emissivity density that comes from the underlying model is unlikely to be the culprit. The solid, dashed, and dotdashed green curves are the analogous fit when we use the emissivity coming from LDP. The conclusions remain unchanged.
Galaxy correlation functions can be reasonably well fitted, over a limited range of scales, by powerlaws. A powerlaw correlation function would project into a powerlaw angular correlation function, so we also consider a powerlaw fit to our data, (42)This simple twoparameter model (A,n) is a reasonable fit at all frequencies, giving a reduced χ^{2} / d.o.f. ≃ 0.21 (1.51/7) at 545 GHz. The bestfit values are given in Table 6 and the bestfit model at 545 GHz is displayed as the blue solid line in Fig. 19.
Powerlaw model bestfit parameters for each frequency as well as the reduced χ^{2}.
Although the powerlaw model provides a good fit to the data, it does not provide physical insight into the properties of the galaxies it describes.
5.5.2. HOD model constraints
We now consider the HOD model introduced in Sect. 5.4. In the most general configuration, our parametrisation of this model allows for four different parameters: M_{min}, M_{sat}, α_{sat} and σ_{log M} (see Appendix C). The full exploration of this parameter space turns out to be a difficult task beyond our scope in this paper. Therefore, we will restrict ourselves to only two parameters, M_{min} and α_{sat}, which describe the mass above which we expect a halo to host a CIBcontributing galaxy and the slope of the highmass end of the HOD. By varying these parameters we control the mean, galaxyweighted halo mass and the satellite fraction in the model, which in turn control the amplitudes of the 1 h and 2 h terms.
We impose M_{sat} = 3.3 M_{min} (Pénin et al. 2011a) and choose σ_{log M} = 0.65, motivated by the clustering observed in optical surveys (see Tinker & Wetzel 2010, for discussion and references). We did not find σ_{log M} to be a critical parameters for our fit, but letting it vary drives the fit to physically unrealistic region of parameter space. For each frequency, we fixed the Poisson noise level to the one expected from our model (see Sect. 3.1). To add to the robustness of our interpretation, we took one more conservative step.
As stated above, our interpretation of these data requires the knowledge of the emissivity. While our fiducial model is well tested and reproduces all relevant current observations (Béthermin et al. 2011), these do not extend beyond a redshift of about 3.5. As such, extrapolations to higher redshifts are unconstrained by previous observations except for the integral constraints provided by the CIB measurement discussed in Sect. 5.1. To let our conclusions be as modelindependent as possible and also to isolate and constrain the highz contribution, we made the extra assumption that the emissivity, j, is constant for z > 3.5 and we fitted for it simultaneously while solving for log _{10}M_{min} and α_{sat}. More precisely, we rewrite Eq. (37) as (43)where we have introduced the effective redshiftindependent emissivity, j_{eff}, that we will also solve for. We adopt the arbitrary but reasonable z = 7 cutoff of Béthermin et al. (2011).
We treated the four frequencies as independent and performed a single minimization per frequency. The results are illustrated in Fig. 20. The solid orange line represents the bestfit model per frequency. Our threeparameter model obviously fits each frequency very well. The orange dashed line represents the 2halo (linear) term, while the orange dotdashed line represents the 1halo (nonlinear) term. The green curve corresponds to the assumed Poisson noise level. Clearly, the angular scales we probe require a modelling of both the linear and the nonlinear contribution to the power spectrum for all frequencies, and the 1halo term is similar in slope to the shotnoise term, which leads to a model degeneracy.
Quantitative results are given in Table 7, where we quote a reduced χ^{2} as a goodnessoffit measure. The errors quoted in Table 7 correspond to the Gaussian errors computed from the Fisher matrix at the bestfit values. Each ℓ bin at any given frequency is considered independent from the others at all frequencies. For reference, we also give the results of a fit where we fixed j to the value given by the Béthermin et al. (2011) model and fitted for only log _{10}M_{min} and α_{sat}. While the bestfit values are consistent between the two models, it is clear from the table that allowing j_{eff} to vary degrades strongly the constraints on log _{10}M_{min} and α_{sat}. In fact, a strong degeneracy between j_{eff} and log _{10}M_{min} is observed, as might be expected.
As is the case with optical data, our data do not appear to require a departure from α_{sat} = 1. We considered different ratios of M_{sat} / M_{min}, i.e., 2, 5, 10 or 20. None of them provided a better fit, and most of them required similar values of α_{sat} (see Tinker & Wetzel 2010, for a recent summary).
Fig. 20 Angular power spectrum of CIB anisotropies at 217, 353, 545, and 857 GHz. Each panel corresponds to one frequency. For each frequency, the blue points correspond to the angular auto powerspectra, and the associated error bars include statistical and photometric calibration systematic contributions. The bestfit model per frequency (including shot noise) corresponds to the solid orange line. The dashed (dotdashed) orange lines correspond to the 2 h (1 h) contributions. The green triple dotdashed curve corresponds to the Poisson noise level, fixed to its expected value. To obtain these fits, three parameters per frequency were varied: log _{10}M_{min}, α_{sat} and j_{eff}. The fits are obviously qualitatively very good. 
Bestfit values for each frequency, as well as the reduced χ^{2}.
Taking our emissivity model at face value, the bestfit angular power spectrum at ℓ = 2000 receives the following redshift contribution at 217 (353/545/857) GHz: 5% (4/34/71) between 0 < z < 1, 7% (7/23/22) between 1 < z < 2, and 88% (89/43/7) for 2 < z.
It is obvious from our Fig. 20 that the nonlinear contribution is degenerate with the Poisson noise level. This explains the problem faced by the linear model discussed in Sect. 5.5.1. Our data by themselves are not sufficient to explore this degeneracy and we thus rely on our model. For similar reasons, we do not discuss details of the implementation of the 1h term (e.g., the truncation radius in u(k,M)) unlike Viero et al. (2009). We expect that a future joint analysis with higher angular resolution measurements from Herschel and SPT/ACT will allow us to alleviate this degeneracy.
As described in Sect. 5.1, our fiducial model predicts a mean CIB consistent with observations. While it is not possible to translate our constraints on (a weighted integral of j^{2} over redshift) into a prediction for the integral of j with the different weighting required to compute the CIB mean, rough estimates suggest that our values are consistent with the FIRAS measurements.
The relatively good consistency between the bestfit values of log _{10}M_{min} and α_{sat} observed across frequencies raises an interesting question: does a single model fit all our data? Or to put it another way, are the differences between each frequency HOD subsantial? Different frequencies, loosely speaking, correspond to different redshifts for the dominant galaxy population. As such, consistency in log _{10}M_{min} and α_{sat} could imply that the CIB fluctuations arise from a single subset of galaxies whose redshift evolution we capture well with our emissivity model, mass function, and HOD prescription, which is constant in redshift. To illustrate this hypothesis, Fig. 21 shows for each frequency’s bestfit model its prediction at 545 GHz. Even though this plot does not convey the uncertainty associated with each prediction, it clearly illustrates that each HOD leads to substantially different predictions and thus that the frequency dependence of the HOD may be significant. We postpone to future work more quantitative statements on the implications for the clustering of galaxies at high redshift.
Fig. 21 Predicted 545 GHz power spectra derived from each frequency’s bestfit model. For the bestfit model at 217, 353, 545 and 857 GHz, we plot the predicted clustering plus shotnoise power spectra at 545 GHz. Also shown are the HFI data points at 545 GHz (red diamonds). This plot suggests that the fits across frequencies are fairly different, which hints at an evolution in the population of galaxies we probed. We note, however, that the uncertainties associated with each prediction are not fully characterised by our method. 
6. Conclusion
We presented the first measurement of CIB anisotropies with Planck, detecting power from 10′ to 2°. Owing to the exceptional quality of the data, and using a complete analysis of the different steps that lead to the CIB anisotropy power spectra, we were able to measure the clustering of dusty, starforming galaxies at 217, 353, 545, and 857 GHz with unprecedented precision.
We worked on six independent fields, chosen to have high angularresolution Hi data and low foreground contamination. The CIB maps were cleaned using templates: Hi for Galactic cirrus; and the Planck 143 GHz maps for CMB. Having Hi data is necessary to cleanly separate CIB and cirrus fluctuations. Because the CIB anisotropies and Galactic dust have similar SEDs, blind component separation methods do not adequately distinguish CIB anisotropies from cirrus emission. The 143 GHz Planck channel, cleaned from sources and filtered, provides a good template for the CMB because it has low instrument noise and an angular resolution close to the higher frequency channels from which we measure the CIB. It also has the advantage of being an “internal” template, meaning its noise, data reduction process, photometric calibration, and beam are well known.
We obtained CIB anisotropy maps that reveal structures produced by the cumulative emission of highredshift, dusty, starforming galaxies. The maps are highly correlated at high frequencies. They decorrelate at lower frequency, as expected from models of the redshift distribution of sources producing the CIB anisotropies (e.g. FernandezConde et al. 2008; Pénin et al. 2011a). In these models, at 217 GHz the contribution of z ≥ 2 galaxies is becoming dominant, while at higher frequencies the dominant sources are at lower z. We computed the power spectra of the maps and their associated errors using a dedicated pipeline, based on the POKER algorithm (Ponthieu et al. 2011). After a careful examination of many systematic effects, use of the best determination of the beam window function and instrument noise determined from jackknife methods, we ended up with measurements of the angular power spectrum of the CIB anisotropy, C_{ℓ}, at 217, 353, 545, and 857 GHz, with high signaltonoise ratio over the range 200 < ℓ < 2000 (see Table 4 and Fig. 12).
The SED of CIB anisotropies is not different from the CIB mean SED, even at 217 GHz. This is expected from the model of Béthermin et al. (2011) and reflects the fact that the CIB mean and anisotropies are produced by the same population of sources. Our measurement compares very well with previous measurements at higher ℓ, at 220 GHz by SPT (Hall et al. 2010) and at 600 and 857 GHz by BLAST (Viero et al. 2009). On the contrary, Herschel/SPIRE measurements from Amblard et al. (2011), but with no sources removal to be comparable to HFI (Amblard, priv. comm.), are significantly lower than our measurements at high frequencies owing to an overestimate of the cirrus contribution. Preliminary crosscalibration between SPIRE and HFI is also increasing the Amblard et al. (2011) power spectra by 10 and 20% at 857 and 545 GHz, respectively.
From the Planck data alone we can exclude a model where galaxies trace the linear theory matter power spectrum with a scaleindependent bias: that model requires an unrealistic high level of shot noise to match the smallscale power we observe. Consequently, we developed an alternative model that couples the dusty galaxy, parametric evolution model of Béthermin et al. (2011) with a halo model approach. Characterised by only two parameters, this model provides an excellent fit to our measured anisotropy angular power spectrum for each frequency treated independently. Whereas in principle these two parameters could offer us unique insights into the clustering and the nature of dusty galaxies at high redshift, the current uncertainties in the underlying model prevent us from drawing detailed inferences. Our results suggest that a different HOD is required at each frequency, which is consistent with the fact that we expect each frequency to be dominated by contributions from different redshifts. We find that half of the contribution to the power spectrum at ℓ = 2000 comes from redshifts lower than 0.8 and 1.5 at 857 and 545 GHz, respectively. Those numbers are quite robust against exact evolution of dusty galaxies comoving emissivity at highredshift (z ≥ 3.5). This is not the case at lower frequencies and our bestfit model predicts that about 90% of the anisotropies power at ℓ = 2000 come from redshifts z > 2 at 353 GHz and 217 GHz.
Further modelling and interpretation of the CIB anisotropy will be aided by the use of crosspower spectra between bands, and by the combination of the Planck and Herschel data at 857 and 545/600 GHz and Planck and SPT/ACT data at 220 GHz. This combination will measure the CIB anisotropy power spectrum over a wide range of scales, covering the three regimes where we expect the 2halo, 1halo and shotnoise contributions to dominate. More progress could be made by measuring the CIB anisotropies over more sky and at lower frequencies (at least 143 GHz) with Planck. Going to lower frequency extends our reach in redshift, and is also important for CMB analysis and measurement of the SZ power spectrum. Additional information will be obtained by crosscorrelating the CIB maps with external tracers of the density field, like the galaxy and quasar distributions in large area catalogues (such as those from the SDSS, VIKING/VISTA and KIDS/VST surveys). This will additionally constrain 1) the populations contributing most to the CIB; and 2) the relative bias between the external tracer and the distribution of farinfrared emission. A particularly interesting crosscorrelation may be between the CMB lensing convergence and the CIB maps (Song et al. 2003). The lensing and CIB anisotropies are expected to have a high degree of overlap, and lead to a signal readily detectable by Planck. This signal will give a direct and independent measure of the bias.
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.
The convention means that the MJysr are given for a source with a spectral energy distribution . For a source with a different spectral energy distribution a colour correction has to be applied (see Planck HFI Core Team 2011b).
The Pearson correlation coefficient is (Lagache et al. 2000).
Amblard et al. (2011) use different values than those given in Swinyard et al. (2010), with beam surfaces of 1.77 × 10^{8} sr and 3.99 × 10^{8} sr at 350 and 500 μm, respectively.
Note that for illustration purpose and where specified only, we will sometimes use the older phenomenological model of Lagache et al. (2004) (LDP).
Acknowledgments
This paper has made use of modelling tools that were made available by Matthieu Béthermin and Aurélie Pénin. The Planck Collaboration acknowledges the support of: 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 DEISA (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.rssd.esa.int/Planck.
References
 Amblard, A., & Cooray, A. 2007, ApJ, 670, 903 [NASA ADS] [CrossRef] [Google Scholar]
 Amblard, A., Cooray, A., Serra, P., et al. 2011, Nature, 470, 510 [NASA ADS] [CrossRef] [Google Scholar]
 ArmitageCaplan, C., & Wandelt, B. D. 2009, ApJS, 181, 533 [NASA ADS] [CrossRef] [Google Scholar]
 Benson, A. J., Cole, S., Frenk, C. S., Baugh, C. M., & Lacey, C. G. 2000, MNRAS, 311, 793 [NASA ADS] [CrossRef] [Google Scholar]
 Berlind, A. A., & Weinberg, D. H. 2002, ApJ, 575, 587 [NASA ADS] [CrossRef] [Google Scholar]
 Bersanelli, M., Mandolesi, N., Butler, R. C., et al. 2010, A&A, 520, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berta, S., Magnelli, B., Lutz, D., et al. 2010, A&A, 518, L30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M., Dole, H., Beelen, A., & Aussel, H. 2010a, A&A, 512, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M., Dole, H., Cousin, M., & Bavouzet, N. 2010b, A&A, 516, A43 [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]
 Blagrave, K., Lockman, F. J., & Martin, P. G. 2010, in ASP Conf. Ser., 438, ed. R. Kothes, T. L. Landecker, & A. G. Willis, 156 [Google Scholar]
 Blain, A. W., Ivison, R. J., & Smail, I. 1998, MNRAS, 296, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Boulanger, F., Abergel, A., Bernard, J. P., et al. 1996, A&A, 312, 256 [NASA ADS] [Google Scholar]
 Cooray, A., & Sheth, R. K. 2002, Phys. Rept., 1 [Google Scholar]
 de Zotti, G., Ricci, R., Mesa, D., et al. 2005, A&A, 431, 893 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delabrouille, J., Cardoso, J.F., Le Jeune, M., et al. 2009, A&A, 493, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dole, H., Rieke, G. H., Lagache, G., et al. 2004, ApJS, 154, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Dole, H., Lagache, G., Puget, J. L., et al. 2006, A&A, 451, 417 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Dunkley, J., Hlozek, R., Sievers, J., et al. 2011, ApJ, 739, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605 [NASA ADS] [CrossRef] [Google Scholar]
 FernandezConde, N., Lagache, G., Puget, J. L., & Dole, H. 2008, A&A, 481, 885 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 FernandezConde, N., Lagache, G., Puget, J., & Dole, H. 2010, A&A, 515, A48 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Franceschini, A., Rodighiero, G., Vaccari, M., et al. 2010, A&A, 517, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gispert, R., Lagache, G., & Puget, J. L. 2000, A&A, 360, 1 [NASA ADS] [Google Scholar]
 Glenn, J., Conley, A., Béthermin, M., et al. 2010, MNRAS, 409, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Granato, G. L., de Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580 [Google Scholar]
 Greve, T. R., Weiss, A., Walter, F., et al. 2010, ApJ, 719, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Grossan, B., & Smoot, G. F. 2007, A&A, 474, 731 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hall, N. R., Keisler, R., Knox, L., et al. 2010, ApJ, 718, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Hauser, M. G., & Dwek, E. 2001, ARA&A, 39, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Hinshaw, G., Nolta, M. R., Bennett, C. L., et al. 2007, ApJS, 170, 288 [NASA ADS] [CrossRef] [Google Scholar]
 Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Jauzac, M., Dole, H., Le Floc’h, E., et al. 2011, A&A, 525, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Knox, L., Cooray, A., Eisenstein, D., & Haiman, Z. 2001, ApJ, 550, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2010, MNRAS, 405, 2 [NASA ADS] [Google Scholar]
 Lagache, G., Abergel, A., Boulanger, F., & Puget, J. L. 1998, A&A, 333, 709 [NASA ADS] [Google Scholar]
 Lagache, G., Abergel, A., Boulanger, F., Désert, F. X., & Puget, J. L. 1999, A&A, 344, 322 [NASA ADS] [Google Scholar]
 Lagache, G., Haffner, L. M., Reynolds, R. J., & Tufte, S. L. 2000, A&A, 354, 247 [Google Scholar]
 Lagache, G., Dole, H., & Puget, J. L. 2003, MNRAS, 338, 555 [NASA ADS] [CrossRef] [Google Scholar]
 Lagache, G., Dole, H., Puget, J. L., et al. 2004, ApJS, 154, 112 [NASA ADS] [CrossRef] [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]
 Larson, D., Dunkley, J., Hinshaw, G., et al. 2011, ApJS, 192, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Le Borgne, D., Elbaz, D., Ocvirk, P., & Pichon, C. 2009, A&A, 504, 727 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leahy, J. P., Bersanelli, M., D’Arcangelo, O., et al. 2010, A&A, 520, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Limber, D. N. 1954, ApJ, 119, 655 [NASA ADS] [CrossRef] [Google Scholar]
 Lueker, M., Reichardt, C. L., Schaffer, K. K., et al. 2010, ApJ, 719, 1045 [NASA ADS] [CrossRef] [Google Scholar]
 Mandolesi, N., Bersanelli, M., Butler, R. C., et al. 2010, A&A, 520, A3 [Google Scholar]
 Marsden, G., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 707, 1729 [NASA ADS] [CrossRef] [Google Scholar]
 Marsden, G., Chapin, E. L., Halpern, M., et al. 2011, MNRAS, 417, 1192 [NASA ADS] [CrossRef] [Google Scholar]
 Marshall, D. J., Robin, A. C., Reylé, C., Schultheis, M., & Picaud, S. 2006, A&A, 453, 635 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matsuura, S., Shirahata, M., Kawada, M., et al. 2011, ApJ, 737, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Mennella, A., Bersanelli, M., Butler, R. C., et al. 2011, A&A, 536, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mitra, S., Rocha, G., Górski, K. M., et al. 2011, ApJS, 193, 5 [NASA ADS] [CrossRef] [Google Scholar]
 MivilleDeschênes, M. A., & Lagache, G. 2005, ApJS, 157, 302 [NASA ADS] [CrossRef] [Google Scholar]
 MivilleDeschênes, M. A., Lagache, G., Boulanger, F., & Puget, J. L. 2007, A&A, 469, 595 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Negrello, M., Perrotta, F., GonzálezNuevo, J., et al. 2007, MNRAS, 377, 1557 [NASA ADS] [CrossRef] [Google Scholar]
 Nguyen, H. T., Schulz, B., Levenson, L., et al. 2010, A&A, 518, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oliver, S. J., Wang, L., Smith, A. J., et al. 2010, A&A, 518, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Patanchon, G., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 707, 1750 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Peacock, J. A., & Smith, R. E. 2000, MNRAS, 318, 1144 [NASA ADS] [CrossRef] [Google Scholar]
 Pearson, C., & Khan, S. A. 2009, MNRAS, 399, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Pénin, A., Doré, O., Lagache, G., & Béthermin, M. 2011a, A&A, in press [Google Scholar]
 Pénin, A., Lagache, G., NoriegaCrespo, A., et al. 2011b, A&A, submitted [Google Scholar]
 Planck Collaboration 2011a, A&A, 536, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011b, A&A, 536, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011c, A&A, 536, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011d, A&A, 536, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011e, A&A, 536, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011f, A&A, 536, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011g, A&A, 536, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011h, A&A, 536, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011i, A&A, 536, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011j, A&A, 536, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011k, A&A, 536, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011l, A&A, 536, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011m, A&A, 536, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011n, A&A, 536, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011o, A&A, 536, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011p, A&A, 536, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011q, A&A, 536, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011r, A&A, 536, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011s, A&A, 536, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011t, A&A, 536, A24 [Google Scholar]
 Planck Collaboration 2011u, A&A, 536, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2011v, The Explanatory Supplement to the Planck Early Release Compact Source Catalogue (ESA) [Google Scholar]
 Planck Collaboration 2011w, A&A, 536, A26 [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]
 Ponthieu, N., Grain, J., & Lagache, G. 2011, A&A, in press, DOI: 10.1051/00046361/201117098 [Google Scholar]
 Puget, J. L., Abergel, A., Bernard, J. P., et al. 1996, A&A, 308, L5 [NASA ADS] [Google Scholar]
 Richter, P., Sembach, K. R., Wakker, B. P., et al. 2001, ApJ, 559, 318 [NASA ADS] [CrossRef] [Google Scholar]
 Righi, M., HernándezMonteagudo, C., & Sunyaev, R. A. 2008, A&A, 478, 685 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rosset, C., Tristram, M., Ponthieu, N., et al. 2010, A&A, 520, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 RowanRobinson, M. 2009, MNRAS, 394, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Scott, D., & White, M. 1999, A&A, 346, 1 [NASA ADS] [Google Scholar]
 Scott, K. S., Yun, M. S., Wilson, G. W., et al. 2010, MNRAS, 405, 2260 [NASA ADS] [Google Scholar]
 Seljak, U. 2000, MNRAS, 318, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, K. M., Zahn, O., & Doré, O. 2007, Phys. Rev. D, 76, 043510 [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, I. B. 1980, ARA&A, 18, 537 [NASA ADS] [CrossRef] [Google Scholar]
 Swinyard, B. M., Ade, P., Baluteau, J.P., et al. 2010, A&A, 518, L4 [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]
 Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Tinker, J. L., Wechsler, R. H., & Zheng, Z. 2010, ApJ, 709, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Valiante, E., Lutz, D., Sturm, E., Genzel, R., & Chapin, E. L. 2009, ApJ, 701, 1814 [NASA ADS] [CrossRef] [Google Scholar]
 Vieira, J. D., Crawford, T. M., Switzer, E. R., et al. 2010, ApJ, 719, 763 [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]
 White, M., Carlstrom, J. E., Dragovan, M., & Holzapfel, W. L. 1999, ApJ, 514, 12 [NASA ADS] [CrossRef] [Google Scholar]
 White, M., Hernquist, L., & Springel, V. 2001, ApJ, 550, L129 [NASA ADS] [CrossRef] [Google Scholar]
 Wilman, R. J., Jarvis, M. J., Mauch, T., Rawlings, S., & Hickey, S. 2010, MNRAS, 405, 447 [NASA ADS] [Google Scholar]
 Zacchei, A., Maino, D., Baccigalupi, C., et al. 2011, A&A, 536, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005, ApJ, 633, 791 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: From HFI maps to CIB power spectra: flow charts of the different steps
We summarize with the two flow charts presented in Figs. A.1 and A.2 the procedures of data preparation and cleaning, and power spectra measurements and errors evaluation.
Appendix B: The parametric backward evolution model of dusty starforming galaxies
In this appendix we give some details about the dusty starforming galaxies evolution model we are using in our modelling of CIB anisotropies. The model is fully described in Béthermin et al. (2011). Two ingredients come into play: the spectral energy distribution (SED) of galaxies and the luminosity function (LF) evolution.
The SEDs are from the templates library of Lagache et al. (2004), which consists of two galaxy populations: a starforming galaxy population with SEDs that vary with IR bolometric luminosities, and a normalgalaxy population with a template SED that is fixed and is colder than the starforming galaxy templates and is scaled with IR bolometric luminosities. The normal and starforming galaxies are dominant at low and highluminosity, respectively. The fraction of each galaxy population as a function of the bolometric luminosity was given using a smooth function (B.1)where th is the hyperbolic tangent function, L_{pop} the luminosity at which the number of normal and starforming galaxies are equal, and σ_{pop} characterises the width of the transition between the two populations. At L_{IR} = L_{pop}, the starbursts fraction is 50%.
They assume that the luminosity (LF) is a classical double exponential function: (B.2)where Φ(L_{IR}) is the number of sources per logarithm of luminosity and per comoving volume unit for an infrared bolometric luminosity L_{IR}, Φ_{ ⋆ } is the normalization constant characterising the density of sources, L_{ ⋆ } is the characteristic luminosity at the break, and 1α and 1 − α − 1 / σ^{2} / ln^{2}(10) are the slope of the asymptotic powerlaw behaviour at low and high luminosity respectively.
A continuous LF redshiftevolution in luminosity and density is assumed following L^{ ⋆ } ∝ (1 + z)^{rL} and Φ^{ ⋆ } ∝ (1 + z)^{rΦ}, where r_{L} and r_{φ} are parameters driving the evolution in luminosity and density, respectively. It is impossible to reproduce the evolution of the LF with constant r_{L} and r_{φ}. They consequently authorize their value to change at two specific redshifts. The position of the first redshift break is a free parameter and converges to the same final value (z ~ 0.9) for initial values 0 < z < 2. To avoid a divergence at high redshift, the second break is fixed at z = 2. The position of the breaks are the same for both r_{L} and r_{φ}.
Fig. A.1 Cartoon illustrating the different steps from HFI frequency maps to CIB power spectra for one field. 
Fig. A.2 Cartoon illustrating the angular power spectrum measurement and errors estimate using the POKER algorithm (Ponthieu et al. 2011). 
The model has 13 free parameters that are summarized in Table B.1. The parameters were determined by fitting the model to published measurements of galaxy number counts and monochromatic LF measured at given redshifts:

Number counts: Spitzer counts at 24, 70 and160 μm, Herschel counts at 250, 350 and 500 μm, and AzTEC counts at 1.1 mm.

Monochromatic LFs: IRAS local LF at 60 μm, Spitzer LF at 24 μm at z = 0, at 15 μm at z = 0.6, at 12 μm at z = 1, and at 8 μm at z = 2.

FIRAS CIB spectrum between 200 μm and 2 mm.
Measured redshift distributions were not used because the cosmic variance and the selection effects were currently poorly quantified. The bestfit parameters as well as their uncertainties and degeneracies were obtained using a Monte Carlo Markov chain (MCMC) MetropolisHastings algorithm. The model adjusted on deep counts and monochromatic LFs at key wavelengths also reproduces recent very discriminating observations well, such as the Jauzac et al. (2011) measured redshift distribution of the CIB.
We used the socalled mean model, which is obtained using the mean value of the parameters as given in Table B.1 without the lensing contribution (dndsnudz_arr_nolensing_meanmodel_final file on the http://www.ias.upsud.fr/irgalaxies/web page).
Appendix C: The halo model
Dusty starforming galaxies evolution model parameters, fitted to selected infrared observations.
In this appendix we give the details of our halo modelling. Neglecting scaledependent halo bias, the distinction between central and satellite galaxies and halo exclusion, and assuming a Poisson distribution of galaxies, the 1 h and 2 h terms in P_{gg} have simple analytic expressions (Cooray & Sheth 2002): Here M is the halo mass, dN / dM is the halo mass function, u(k,M) is the Fourier transform of the (normalized) halo density profile, b(M) the halo bias and ⟨ N_{gal} ⟩ is the mean number of galaxies in a halo of mass M.
The mean number density of galaxies, , can be written (C.3)Note that on large scales u(k → 0,M) ≃ 1 so that we can define the “effective” bias as (C.4)and the 1halo term becomes scaleindependent (i.e., a shotnoise term).
We used the fitting function of Eisenstein & Hu (1998) to compute P_{lin}, an NFW profile (Navarro et al. 1997) truncated at the virial radius to compute u(k,M), and we relied on the mass function fit of Tinker et al. (2008) with its associated halo bias prescription (Tinker et al. 2010). All these relations were calibrated through the use of Nbody simulations. Our definition of halo mass is the mass interior to a radius within which the mean density is 200 times the mean density of the Universe.
The HOD describes the way galaxies populate the dark matter halos. While we do not distinguish between central and satellite galaxies in the above, the functional form we adopt for the mean occupation is modelled on the form frequently used in optical observations (e.g., Zheng et al. 2005) (C.5)with (C.6)and (C.7)These definitions ensure that M_{min} is the halo mass at which a halo has a probability of 50% of having a central galaxy. We introduce σ_{log M} to allow scatter in this relation between halo mass and observable, which is important on large scales. Following Zheng et al. (2005), we assume a Poisson distribution for N_{sat} and write (C.8)Within this parametrisation, halos with M ≪ M_{min} will not host any galaxies, whereas those with M ≫ M_{min} are almost certain to contain one. The satellite occupation has a similar cutoff, but the mass is chosen to be twice M_{min}, so that halos with a low probability of having a central galaxy are unlikely to contain a satellite galaxy (see Tinker & Wetzel 2010, for a further discussion of this form).
We note that this parametrisation was introduced to reproduce the observed clustering of luminositythreshold samples of optical galaxies at 0 < z < 2. We are therefore making substantial assumptions when applying this same parametric form to dusty starforming galaxies at higher redshift. As we shall see, however, our constraints on even this form of the HOD are weak enough to argue against introducing additional degrees of freedom in the model.
All Tables
CIB field description: centre (Galactic coordinates), size, mean and dispersion of Hi column density.
Flux cut from the ERCSC for our six fields, and the shotnoise power for dusty and radio galaxies appropriate to those cuts (see text).
Powerlaw model bestfit parameters for each frequency as well as the reduced χ^{2}.
Dusty starforming galaxies evolution model parameters, fitted to selected infrared observations.
All Figures
Fig. 1 From left to right and top to bottom: N1, AG, SP, LH2 and bootes fields overlaid on IRIS 100μm map (MivilleDeschênes & Lagache 2005). Fields Bootes 1 and 2 are both included in the large rectangle. All IRIS images have the same dynamic range, with a linear colour scale ranging from dark red to white from 0 to 2 MJysr^{1}. 

In the text 
Fig. 2 Wiener filter applied to the 143 GHz map for CMB subtraction. The filter essentially cuts out high multipoles where the CMBtonoise ratio of the 143 GHz map is low. Whereas this filter has to be known for estimating and subtracting the contribution of the residual CMB and 143 GHz instrument noise to the power spectrum of CMBcleaned channels, the exact value of the filter is not really critical. 

In the text 
Fig. 3 Power spectra of the different components for field SP (the figure is similar for the other fields). Power spectra of the 217, 353, 545, and 857 GHz Planck maps (continuous black line) are compared to the noise power spectra (diamonds), to the CMBcleaned power spectra (red), and to the CMB and interstellar dustcleaned power spectra (green). In this plot signal power spectra have not been corrected for the beam window function. Noise power spectra are computed using halfpointing period maps, as explained in Sect. 3.3. 

In the text 
Fig. 4 Maps of the 26 deg^{2} of the N1 field, from left to right: 217, 353, 545 and 857 GHz. From top to bottom: raw HFI maps; CMB and ERCSC sourcecleaned maps; residual maps (CMB, sources, and cirruscleaned); residual maps smoothed at 10′ to highlight the CIB anisotropies. The joint structures clearly visible (bottom row) correspond to the anisotropies of CIB. Residual point sources are also visible. They have fluxes lower than the fluxes of the ERCSC removed sources. They have no impact on our analysis. 

In the text 
Fig. 5 Hi and dust maps for two fields: SP (top)and AG (bottom). The first two maps on the left show the Hi components (Local and IVC for SP, IVC and HVC for AG), the third maps show the 857 GHz emission associated with Hi () and the maps on the right side show the HFI 857 GHz maps. Those HFI maps have been convolved by the GBT beam to allow a better comparison by eye. Hi maps are given in units of 10^{20} atomscm^{2}. Note the correlation of the dust emission with the different Hi velocity components and its variation from field to field. 

In the text 
Fig. 6 Contribution to the CIB per redshift slice, extracted from Béthermin et al. (2011). The black solid line is the CIB spectrum predicted by the model. The contribution to the CIB from 0 < z < 0.3, 0.3 < z < 1, 1 < z < 2 and z > 2 galaxies is given by the red shortdashed, green dotdashed, blue three dotdashed and purple longdashed lines, respectively. Lower limits coming from the stacking analysis at 100 μm, 160 μm (Berta et al. 2010), 250 μm, 350 μm, 500 μm (Marsden et al. 2009), 850 μm (Greve et al. 2010) and 1.1 mm (Scott et al. 2010) are shown as black arrows. The black diamonds give the Matsuura et al. (2011) absolute measurements with AKARI. The black square the Lagache et al. (2000) absolute measurements with DIRBE/WHAM and the cyan line the Lagache et al. (2000) FIRAS measurement. 

In the text 
Fig. 7 A number of recent models of dustygalaxy evolution and their associated shot noise for different flux cuts at 857 GHz. Top: comparison of the models with the Herschel and BLAST differential numbers counts. Models are from Lagache et al. (2004), Negrello et al. (2007), Le Borgne et al. (2009), Patanchon et al. (2009), Pearson & Khan (2009), Valiante et al. (2009), Béthermin et al. (2011), Franceschini et al. (2010), Lacey et al. (2010), Marsden et al. (2011), RowanRobinson (2009), Wilman et al. (2010). Data points are from Oliver et al. (2010), Béthermin et al. (2010b), Glenn et al. (2010). Bottom: shotnoise level as a function of the flux cut for the same models (same colour and line coding between the two figures). The vertical and horizontal continuous dark lines show the Planck flux cut and shotnoise level from Table 3, respectively. The Béthermin et al. (2011) model is shown by the continuous dark line. This figure shows that models predicting a very high shot noise (e.g. continuous and dashed lightblue, reddashed, continuous and dashed darkblue lines) are incompatible with the measured number counts. 

In the text 
Fig. 8 Effective beam window functions (b_{ℓ}) from FICSBell (black) and FEBeCoP (red) at 545 GHz (see Sect. 3.2 for more details). The six FEBeCoP beam window functions from each field are superimposed (red lines). Also shown for comparison is the Gaussian beam with a FWHM of 4.72′ ± 0.2′ (green lines), which is the equivalent FWHM of the beam determined on Mars. 

In the text 
Fig. 9 Three independent noisepowerspectrum measurements in the SP field at 353 GHz: red continuous line, half pointing period; green dashed, surveys I and II; black dotdashed, half focal plane array). 

In the text 
Fig. 10 Instrument noise power spectra of the six fields obtained using halfpointing period maps. From top to bottom: 217, 353, 545 and 857 GHz (continuous: N1, dotted: AG, dashed: SP, dashdotted: Bootes 1, longdash: Bootes 2, dash3 dotted: LH2). 

In the text 
Fig. 11 Contribution of residuals to the final CIB anisotropy estimate (illustrated here with the SP field at 353 GHz). The bias induced by each dust and CMB component is negligible compared to both the CIB anisotropy signal (black dots) and the statistical noise (in green, including cosmic variance on the noise estimate itself). 

In the text 
Fig. 12 CIB anisotropy power spectra of the six fields and their combined spectrum. 

In the text 
Fig. 13 Modulus of the bintobin correlation matrix derived from the simulation pipeline for the SP field at 353 GHz. Spatial frequency bins are not correlated by more than ≃ 10%. 

In the text 
Fig. 14 Fieldcombined CIB anisotropy power spectra at 217, 353, 545, and 857 GHz. The dashed line shows the expected sum of the dusty and radio galaxy shotnoise power (from Table 3). The power spectra at 217, 353, and 545 GHz were arbitrary scaled to allow for a better comparison between frequencies (they were multiplied by 2 × 10^{6}, 10^{5} and 10^{3}, respectively). 

In the text 
Fig. 15 Comparison of the observed CIB mean and anisotropy SED. The CIB measurements are from Lagache et al. (1999) (FIRAS spectrum in black) and Pénin et al. (2011b) (Spitzer and IRIS, pink diamond data points). The green and blue continuous (dashed) lines are the CIB fits from Gispert et al. (2000) and Fixsen et al. (1998) (multiplied by 0.15). The rms fluctuations of the CIB anisotropies, measured for 200 < ℓ < 2000, are shown with the red dots. Their error bars include both statistical and photometric calibration systematic errors (linearly added), as given in Table 5. This figure shows that the CIB anisotropy SED is steeper than the Fixsen et al. (1998) best fit but very close to the Gispert et al. (2000) best fit. We see no evidence for different CIB mean and anisotropy SED. 

In the text 
Fig. 16 Comparison of SPT (Hall et al. 2010, dark open diamonds) and HFI measurements (red dots) at 217 GHz. The green dashed line corresponds to the SPT shot noise and the green dotdashed line to the clustering model of Hall et al. (2010), the sum of the two is the continuous green line. The clustering model overpredicts by a factor ≃ 2.4 the HFI power at ℓ ~ 800. The blue dashdotted line shows the clustering model divided by this factor. The clustering+shot noise (blue continuous line) now underpredicts the SPT data points, which may be the signature of nonlinear contributions. 

In the text 
Fig. 17 Comparison of BLAST and HFI measurements at 545 and 857 GHz. HFI data points are the red circles; BLAST data points are the black diamonds. They were colourcorrected for the comparison (the colour was computed using the CIB SED of Gispert et al. (2000), integrated through the BLAST and HFI bandpass filters). The dashed line is the BLAST shot noise (also colourcorrected). Also shown is the BLAST bestfit clustering model (black dashdotted line) and the total contribution (shot noise plus clustering; continuous green line). It provides a good fit to the Planck data. Finally, we report in this figure a revised version of the SPIRE data points from Amblard et al. (2011) (blue triangles from Fig. 18, see text for more details). 

In the text 
Fig. 18 Comparison of SPIRE and HFI measurements at 545 and 857 GHz in the overlapping multipole range. HFI data points are the red circles; SPIRE data points from Amblard et al. (2011) are the black triangles. For these data points sources down to 50 mJy have been masked, there is consequently less power compared to HFI. The green triangles (Amblard, priv. comm.) show the SPIRE CIB measurements identical to Amblard et al. (2011), but without a flux cut applied and thus they are directly comparable to the HFI measurement. They should agree with HFI, but are a factor ~1.7 and ~1.2 below the HFI CIB data points for 400 < ℓ < 1500 at 857 and 545 GHZ, respectively. Indeed, they suffer from an overestimate of the cirrus contamination (by a factor 2). Moreover, preliminary crosscalibration between SPIRE and HFI is increasing the Amblard et al. (2011) SPIRE power spectra by 10 and 20% at 857 and 545 GHz, respectively (see Sect. 5.3 for more details). When corrected from these too factors (cirrus and crosscalibration), the SPIRE (blue triangles) and HFI measurements agree well. For this figure, all SPIRE data points were colourcorrected (colours were computed using the CIB SED of Gispert et al. (2000), integrated through the SPIRE and HFI bandpass filters). Error bars include only statistical errors (for SPIRE, error bars are only shown for the green triangles for sake of clarity). 

In the text 
Fig. 19 In this plot we illustrate the constant bias model at 545 GHz. The orange (solid, dashed and dotdashed) lines correspond to the bestfit linear model, its clustering component, and the shotnoise level, respectively. While this fit was performed using our fiducial emissivity defined in Eq. (38), for illustrative purpose we plot in green the analogous fit using the LDP emissivity. In both cases the required shotnoise level (dotdashed lines) is well above the 68% C.L. predicted by Béthermin et al. (2011) and given in Table 5 (yellow contour). Conversely, the solid yellow line represents the bestfit curve when the shotnoise level is fixed to the expected value (Table 3) and only the constant bias is varied. The fit is obviously unsatisfactory. These results lead us to consider the linear bias model as unphysical, despite the good fit it provides (χ^{2}/d.o.f. ≃ 0.36 (2.52/7)). For illustration purpose, we also show our bestfit powerlaw model defined in Eq. (42) (blue solid line). 

In the text 
Fig. 20 Angular power spectrum of CIB anisotropies at 217, 353, 545, and 857 GHz. Each panel corresponds to one frequency. For each frequency, the blue points correspond to the angular auto powerspectra, and the associated error bars include statistical and photometric calibration systematic contributions. The bestfit model per frequency (including shot noise) corresponds to the solid orange line. The dashed (dotdashed) orange lines correspond to the 2 h (1 h) contributions. The green triple dotdashed curve corresponds to the Poisson noise level, fixed to its expected value. To obtain these fits, three parameters per frequency were varied: log _{10}M_{min}, α_{sat} and j_{eff}. The fits are obviously qualitatively very good. 

In the text 
Fig. 21 Predicted 545 GHz power spectra derived from each frequency’s bestfit model. For the bestfit model at 217, 353, 545 and 857 GHz, we plot the predicted clustering plus shotnoise power spectra at 545 GHz. Also shown are the HFI data points at 545 GHz (red diamonds). This plot suggests that the fits across frequencies are fairly different, which hints at an evolution in the population of galaxies we probed. We note, however, that the uncertainties associated with each prediction are not fully characterised by our method. 

In the text 
Fig. A.1 Cartoon illustrating the different steps from HFI frequency maps to CIB power spectra for one field. 

In the text 
Fig. A.2 Cartoon illustrating the angular power spectrum measurement and errors estimate using the POKER algorithm (Ponthieu et al. 2011). 

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.