© ESO, 2011

1. Introduction

In addition to instrument noise, deep cosmological surveys in the far-infrared 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; Fernandez-Conde 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¹ (Fernandez-Conde 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 far-infrared 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 large-scale 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, star-forming 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 star-formation 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 clustering-related observable. In particular, the halo model predicts that the bias, describing the clustering of galaxies in relation to the dark matter, becomes scale-independent at large scales. This assumption of a scale-independent 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, star-forming 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 scale-independent 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 white-noise 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, star-forming 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 large-scale, linear bias, b, of dusty galaxies, which is usually assumed to be scale-independent over the relevant range. Given their limited dynamic range in scale, current measurements are equally consistent with an HOD model, a power-law correlation function or a scale-independent, linear bias. All models return a value for the large-scale bias that is 2–4 times higher than that measured for local, dusty, star-forming galaxies (where b ≃ 1).

Owing to its frequency coverage from 100 to 857 GHz, the High Frequency Instrument (HFI) on-board Planck is ideally suited to probe the dark matter – star-formation connection. Planck (Tauber et al. 2010; Planck Collaboration 2011a) is the third-generation 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 star-forming galaxies. To achieve this, we analyse small regions of sky with a total area of about 140deg², 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 power-law 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 high-redshift, 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 star-forming 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 time-ordered 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). Time-ordered 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 non-linearity 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 full-width at half-maximum (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: ring-making, ring offset estimation, and map-making. 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 signal-to-noise ratio, taking full advantage of the redundancy of observations provided by the Planck scanning strategy. The low-frequency 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 offset-corrected 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. Inter-calibration accuracy between channels is better than absolute calibration.

We made use of the so-called 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.²) or μK_CMB, the conversion between the two was exactly computed using the bandpass filters (see Planck HFI Core Team 2011b).

2.2. Extragalactic fields with high angular resolution HI data

Although Planck is an all-sky 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 21-cm Hi spectra used here were obtained with the 100-m Green Bank Telescope (GBT) over the period 2005 to 2010. Details of this high-latitude survey are presented by Martin et al. (in prep.). The total area mapped is about 825deg².

The spectra were taken with on-the-fly 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 in-scan 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 in-band 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 third-order polynomial was fitted to the emission-free 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 high-velocity 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 all-sky sight line (referenced as LH2 in the Table). The field areas are between 16 to 25deg² for a total coverage of about 140deg². Going to higher average Hi column densities (N(Hi) > 2.5 × 10²⁰ 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 shot-noise 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 gap-filling 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 high-frequency CIB anisotropy signals will not leak into lower frequency, CMB-free 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 best-fit to the effective HFI beam of the channel maps.

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(ν₂₁₇), 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 best-fit CMB model spectrum, and Y_ℓ(ν₁₄₃) 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 sub-optimal filtering rather than in biases. We checked that the CMB remaining in the CMB-cleaned 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 CMB-cleaned 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 background-cleaned 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, large-scale 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 far-infrared emission from dust and the 21-cm emission from gas³ (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.

Hi components –  The Hi data in each field show different velocity components: a local component, typical of high-latitude Hi emission, intermediate-velocity clouds (IVCs) and high-velocity 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. High-velocity clouds might be infalling clouds fueling the Galaxy with low-metallicity 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 21-cm spectrum and of the rms 21-cm 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 21-cm 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.

HI-dust correlation –  To remove the cirrus contamination from the HFI maps, we need to determine the far-IR 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 far-IR 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 χ² 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 power-law power spectrum (Miville-Deschê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 high-frequency 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 high-redshift 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 Fernandez-Conde et al. (2008, 2010) and Pénin et al. (2011a).

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 star-forming galaxies (with both shot noise, , and clustering, , components), and one from radio galaxies (with only a shot-noise 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 Sunyaev-Zeldovich (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, star-forming and radio galaxies in our case).

3.1.1. Star-forming, dusty galaxy shot noise,

We used the recent model of Béthermin et al. (2011) to compute the IR galaxy shot-noise 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/sub-mm number counts, and some mid-IR luminosity functions, with a Monte-Carlo Markov Chain (MCMC). More details on the model are given in Appendix B. The derived shot-noise 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 shot-noise 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 shot-noise 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 shot-noise power. Second, we see that the highest shot-noise level is around 13500Jy²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 shot-noise level below 8000Jy²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 mid-infrared to the millimeter, including the differential contribution of the S₂₄ ≥ 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 Lockman-SWIRE field, when none of the point sources is removed (Amblard, priv. comm.), as detailed in Sect. 5.3.

3.1.2. Radio galaxy shot noise,

The shot-noise power from radio galaxies is subdominant to that from dusty sources at the frequencies relevant to CIB anisotropy analysis. The radio galaxy shot-noise 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 scanning-related 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 (Mars-based) 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 (Armitage-Caplan & Wandelt (2009) reached a similar conclusion in their analysis of Planck-LFI beams).

• 3.

  The coefficients computed above are used to generate s-spin weighted maps for a given CMB sky realisation.

• 4.

  The spin-weighted 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.

Monte-Carlo (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 time-ordered 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 Fourier-Legendre 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 Monte-Carlo 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 FEBeCoP-generated 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 Gauss-Hermite expansion on destriped data (see Planck HFI Core Team 2011b, for more details). We see a quite good agreement between the FICSBell, all-sky and FEBeCoP, small-field 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.

3.3. Instrument noise, N_ℓ(ν)

We can use three different jack-knife difference maps to derive noise power spectra: maps made from the first and second halves of each pointing period (a half-pointing 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 half-pointing 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 half-focal 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 low-frequency noise. Removing the ERCSC sources has no impact on the noise determination.

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 signal-to-noise 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 cirrus-cleaned) 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_ℓ(ν₁₄₃) is the noise power spectrum of the 143 GHz map. It is computed as the noise in the other frequency channels, using the half-pointing 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 Hi-correlated dust emission, which is not the goal of this paper). We expect the shot-noise 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 “re-beaming” 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 shot-noise 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 CIB-correlated 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² and  <9.8 × 10^-6μK² 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 pseudo-spectrum 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 pseudo-power 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 flat-sky approximation P(ℓ). As already suggested, this final estimate makes use of the power spectrum of the masked data. This so-called pseudo-power spectrum will be denoted . In the flat-sky approximation, the standard angular frequencies labelled by their zenithal and azimuthal numbers, usually called ℓ and m respectively, are replaced by an “angular” wave-vector ℓ; 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(ℓ)|² computed from our CIB maps does not depend on the direction of ℓ.

4.1. POKER

The POKER implementation of the pseudo-spectrum 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 wave-vector ℓ, 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 square-modulus of its Fourier coefficients, i.e., P(ℓ_mn) = |D_mn|².

The direct DFT of the masked data relates the true Fourier coefficients to the pseudo-Fourier 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 pseudo-power spectrum). For a signal T plus noise N map, the ensemble averaged of the pseudo-spectrum tracing a statistically isotropic process, reads (18)where we have introduced the total transfer function F_m′n′ accounting for the beam, the “map-making” 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 “map-making” pixelisation effects are described by the power spectrum of the pixel window function for full-sky 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), time-domain filtering is included as part of the scanning beam, such that any time-domain 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 Monte-Carlo: we first generated a set of full-sky maps and computed the MC average of their pseudo-spectra. This set of maps was then re-projected onto flat maps for which MC average of their pseudo-spectra in the flat-sky approximation were computed. The ratio of the flat-sky pseudo-spectrum divided by the full-sky pseudo-spectrum gives a measurement of re-projection 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 one-dimensional 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 pseudo-power 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 re-written 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⁴.

• 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 Monte-Carlo 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 Monte-Carlo averaging. The error bar on each P_b is (27)and the bin-bin 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 extra-source 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 inter-calibration 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 pseudo-power spectrum of the instrument noise and CMB residual used to debias the simulated data pseudo-spectrum 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 “foreground-cleaned” maps to obtain a final estimate of the CIB angular power spectrum. In this step, the bias involved in Eq. (25) contains the pseudo-power spectrum of the instrument noise model and of the CMB residual model as described in the simulation pipeline;

• 4.

  the total error bars and bin-bin correlation matrix on P_b are obtained as the RMS of 1000 Monte-Carlo 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 FEBeCoP-derived beam profiles, but the estimated power spectra using either the FICSBell-derived or the FEBeCoP-derived 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 bin-bin correlation matrix, showing that two 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 Monte-Carlo 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 zero-padded. There is no general prescription as to the size of the zero-padding 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 pseudo-spectrum 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 frequency-dependent 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 MC-variance 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 inter-calibration uncertainties, which are about 2% for the CMB channels (Planck HFI Core Team 2011b). We also fitted the low-frequency channels for only CMB (32)We found that β₁ and β₂ 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 inter-calibration 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 CMB-cleaned maps changes the dust-Hi 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 cross-correlation between half-pointing period maps instead of the auto-correlation. 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 bin-bin covariance of the “combined-field” estimator, , is a function of the bin-bin covariance of each “single-field” estimator, , and the covariance between two “single-field” 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 field-to-field and bin-to-bin 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 bin-bin covariance matrix could be computed⁵. However, the forthcoming cosmological interpretation of the derived power spectra assumes a zero bin-to-bin 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 signal-to-noise compared to the single-field measurements.

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 shot-noise 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 shot-noise 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 best-fit 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.

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 over-plot 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, non-evolving, 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 smaller-scale 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 non-linear corrections.

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 GOODS-South field. They detect CIB anisotropy and shot-noise power in the range 940 ≤ ℓ ≤ 10800. The measured correlations are well fitted by a power-law 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 star-forming 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 largest-scale 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 shot-noise powers measured by BLAST. Once the colour corrections are applied, they are 1843 and 7326Jy²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 best-fit 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!

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 sub-arcminute angular scales or ℓ of  ~ 2 × 10⁴. For clarity we do not plot the small-scale 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 SWIRE-Lockman 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 cross-calibration: SPIRE data are calibrated for point sources, with an accuracy of 15% (Swinyard et al. 2010). The point-source to diffuse-emission 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 cross-calibration than the official beam surfaces given in Swinyard et al. (2010). This preliminary cross-calibration has been established by comparing the diffuse emission from several Hershel surveys (some HATLAS, SAG4, and Hi-Gal fields) to Planck/HFI data. Compared to the beam surfaces taken in Amblard et al. (2011)⁶, 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 cross-calibration 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 self-consistent, 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 self-consistent manner the observational constraints coming from galaxy clustering and the CMB with more galaxy-evolution-centered 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) infrared-galaxy 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 mid-IR 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)⁷: (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 scale-independent bias and P_lin(k) is the linear theory, dark-matter 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 non-linearities and scale-dependent 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 1-halo term becomes a scale-independent, shot-noise 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 best-fit 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 Levenberg-Marquardt algorithm to perform a χ² minimization. The χ² 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 best-fit value is equivalent to marginalising over them. It is important to note that the model power spectrum, includes a shot-noise term (SN) defined in Sect. 3.1. Depending on the precise configuration we study, this shot-noise 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 colour-corrected to be in Jy(νI_ν = const.) as the data points, using the Gispert et al. (2000) CIB SED (colour corrections are (1.08)², (1.08)², (1.06)² and (1.00)² 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.

5.5.1. Linear bias model and power-law 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 shot-noise level and a constant bias, b_lin, defined in Eq. (39), we obtain a good fit (χ² / 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 shot-noise level required to obtain this good fit (dot-dashed orange curve) is unrealistically high: (5.6 ± 0.7) × 10³ Jy²/sr⁸. 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 shot-noise 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 shot-noise level we fit here absorbs a strong non-linear component.

To further illustrate this point, we show the solid yellow curve that represents the best-fit 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 dot-dashed 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 power-laws. A power-law correlation function would project into a power-law angular correlation function, so we also consider a power-law fit to our data, (42)This simple two-parameter model (A,n) is a reasonable fit at all frequencies, giving a reduced χ² / d.o.f. ≃ 0.21 (1.51/7) at 545 GHz. The best-fit values are given in Table 6 and the best-fit model at 545 GHz is displayed as the blue solid line in Fig. 19.

Although the power-law 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 CIB-contributing galaxy and the slope of the high-mass end of the HOD. By varying these parameters we control the mean, galaxy-weighted 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 model-independent as possible and also to isolate and constrain the high-z 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 ₁₀M_min and α_sat. More precisely, we rewrite Eq. (37) as (43)where we have introduced the effective redshift-independent emissivity, j_eff, that we will also solve for. We adopt the arbitrary but reasonable z = 7 cut-off 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 best-fit model per frequency. Our three-parameter model obviously fits each frequency very well. The orange dashed line represents the 2-halo (linear) term, while the orange dot-dashed line represents the 1-halo (non-linear) 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 non-linear contribution to the power spectrum for all frequencies, and the 1-halo term is similar in slope to the shot-noise term, which leads to a model degeneracy.

Quantitative results are given in Table 7, where we quote a reduced χ² as a goodness-of-fit measure. The errors quoted in Table 7 correspond to the Gaussian errors computed from the Fisher matrix at the best-fit 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 ₁₀M_min and α_sat. While the best-fit 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 ₁₀M_min and α_sat. In fact, a strong degeneracy between j_eff and log ₁₀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).

Taking our emissivity model at face value, the best-fit 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 non-linear 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² 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 best-fit values of log ₁₀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 ₁₀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 best-fit 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.

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, star-forming galaxies at 217, 353, 545, and 857 GHz with unprecedented precision.

We worked on six independent fields, chosen to have high angular-resolution 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 high-redshift, dusty, star-forming 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. Fernandez-Conde 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 jack-knife 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 signal-to-noise 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 cross-calibration 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 scale-independent bias: that model requires an unrealistic high level of shot noise to match the small-scale 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 high-redshift (z ≥ 3.5). This is not the case at lower frequencies and our best-fit 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 cross-power 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 2-halo, 1-halo and shot-noise 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 cross-correlating 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 far-infrared emission. A particularly interesting cross-correlation 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.

---

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/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN and JA (Spain); Tekes, AoF and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and 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

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 star-forming galaxies

In this appendix we give some details about the dusty star-forming 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 star-forming galaxy population with SEDs that vary with IR bolometric luminosities, and a normal-galaxy population with a template SED that is fixed and is colder than the star-forming galaxy templates and is scaled with IR bolometric luminosities. The normal and star-forming galaxies are dominant at low- and high-luminosity, 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 star-forming 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 / σ² / ln²(10) are the slope of the asymptotic power-law behaviour at low and high luminosity respectively.

A continuous LF redshift-evolution in luminosity and density is assumed following L ^⋆  ∝ (1 + z)^r_L 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_φ.

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 best-fit parameters as well as their uncertainties and degeneracies were obtained using a Monte Carlo Markov chain (MCMC) Metropolis-Hastings 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 so-called 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.u-psud.fr/irgalaxies/web page).

Appendix C: The halo model

In this appendix we give the details of our halo modelling. Neglecting scale-dependent 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 1-halo term becomes scale-independent (i.e., a shot-noise 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 N-body 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 cut-off, 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 luminosity-threshold samples of optical galaxies at 0 < z < 2. We are therefore making substantial assumptions when applying this same parametric form to dusty star-forming 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

All Figures